Back to index

salome-med  6.5.0
dumpMEDMEM.py
Go to the documentation of this file.
00001 #  -*- coding: iso-8859-1 -*-
00002 # Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00003 #
00004 # This library is free software; you can redistribute it and/or
00005 # modify it under the terms of the GNU Lesser General Public
00006 # License as published by the Free Software Foundation; either
00007 # version 2.1 of the License.
00008 #
00009 # This library is distributed in the hope that it will be useful,
00010 # but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012 # Lesser General Public License for more details.
00013 #
00014 # You should have received a copy of the GNU Lesser General Public
00015 # License along with this library; if not, write to the Free Software
00016 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00017 #
00018 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00019 #
00020 
00021 ############################################################################
00022 # This file provides utilities to dump content of MEDMEM objects: mesh,
00023 # field, group, family, nodal connectivity, node coordinates.
00024 ############################################################################
00025 #
00026 from libMEDMEM_Swig import *
00027 
00028 import os
00029 
00030 theEntityName = { MED_CELL        :"CELL",
00031                   MED_FACE        :"FACE",
00032                   MED_EDGE        :"EDGE",
00033                   MED_NODE        :"NODE",
00034                   MED_ALL_ENTITIES:"ALL_ENTITIES" }
00035 
00036 theTypeName = {MED_NONE        :"NONE",
00037                MED_POINT1      :"POINT1",
00038                MED_SEG2        :"SEG2",
00039                MED_SEG3        :"SEG3",
00040                MED_TRIA3       :"TRIA3",
00041                MED_QUAD4       :"QUAD4",
00042                MED_TRIA6       :"TRIA6",
00043                MED_QUAD8       :"QUAD8",
00044                MED_TETRA4      :"TETRA4",
00045                MED_PYRA5       :"PYRA5",
00046                MED_PENTA6      :"PENTA6",
00047                MED_HEXA8       :"HEXA8",
00048                MED_TETRA10     :"TETRA10",
00049                MED_PYRA13      :"PYRA13",
00050                MED_PENTA15     :"PENTA15",
00051                MED_HEXA20      :"HEXA20",
00052                MED_POLYGON     :"POLYGON",
00053                MED_POLYHEDRA   :"POLYHEDRA",
00054                MED_ALL_ELEMENTS:"ALL_ELEMENTS"}
00055 
00056 medModeSwitch = { 0:"FULL_INTERLACE",
00057                   1:"NO_INTERLACE",
00058                   3:"UNDEFINED_INTERLACE" };
00059 
00060 med_type_champ = { 6 : "REEL64",
00061                    24: "INT32",
00062                    26: "INT64",
00063                    0 : "UNDEFINED_TYPE" } ;
00064 
00065 tab="   "
00066 
00067 debugShowConn=True
00068 debugShowConn=False
00069 
00070 SHOW_ALL = -1
00071 
00072 # private
00073 def _showNodalConnectivity(mesh,entity,type,elems,tablevel,showOnly=SHOW_ALL):
00074     if debugShowConn: print "ELEMENTS:",elems
00075     tab1 = tab*tablevel
00076     tab2 = tab*(tablevel+1)
00077     typeName = theTypeName[type]
00078     nbElem = len( elems )
00079     if showOnly > 0:
00080         elems = elems[:showOnly]
00081     nbShow = len( elems )
00082     connectivity = mesh.getConnectivity(MED_NODAL,entity,MED_ALL_ELEMENTS)
00083     index = mesh.getConnectivityIndex(MED_NODAL,entity)
00084     if debugShowConn: print "CONN:",connectivity,"\nIND:",index
00085     elemShift = 0
00086     types = mesh.getTypes( entity )
00087     for t in types:
00088         if t != type:
00089             elemShift += mesh.getNumberOfElements(entity,t)
00090         else:
00091             break
00092         pass
00093     for i in elems:
00094         elem = i + elemShift
00095         print tab1,typeName,i,":",connectivity[index[elem-1]-1 : index[elem]-1]
00096         pass
00097     nbSkip = nbElem - nbShow
00098     if nbSkip > 0:
00099         print tab1,"...",nbSkip,"elements not shown"
00100         pass
00101     pass
00102 
00103 #private
00104 def _showSupport(support, tablevel,showElems=0):
00105     tab1 = tab*(tablevel+0)
00106     tab2 = tab*(tablevel+1)
00107     tab3 = tab*(tablevel+3)
00108     entity    = support.getEntity()
00109     types     = support.getTypes()
00110     nbOfTypes = support.getNumberOfTypes()
00111     onAll     = support.isOnAllElements()
00112     print tab1,"-Entity:",theEntityName[entity]
00113     print tab1,"-Types :",types
00114     print tab1,"-Elements"
00115     if onAll:
00116         print tab2,"<< Is on all elements >>"
00117     else:
00118         for type in types:
00119             nbOfElmtsOfType = support.getNumberOfElements(type)
00120             number = support.getNumber(type)
00121             number.sort()
00122             print tab2,"* Type:",theTypeName[type]
00123             print     tab3,". Nb elements:",nbOfElmtsOfType
00124             print     tab3,". Numbers:",number[:min(100,nbOfElmtsOfType)],
00125             if nbOfElmtsOfType > 100:
00126                 print "...skip", nbOfElmtsOfType-100
00127             else:
00128                 print
00129                 pass
00130             if entity != MED_NODE and showElems:
00131                 print tab3,". Nodal connectivity"
00132                 _showNodalConnectivity(support.getMesh(),entity,type,number,tablevel+4,showElems)
00133                 pass
00134             pass
00135         pass
00136     print
00137     return
00138 
00139 ## Dump i-th family of given entity in mesh.
00140 ## Optionally dump nodal connectivity of <showElems> first elements.
00141 ## Use showElems=SHOW_ALL to dump connectivity of all elements.
00142 
00143 def ShowFamily(mesh, entity, i, showElems=0):
00144     family = mesh.getFamily(entity,i)
00145     familyName = family.getName()
00146     familyDescription = family.getDescription()
00147     entity = family.getEntity()
00148     familyBool = family.isOnAllElements()
00149     print "\nFAMILY", i, "on", theEntityName[entity]
00150     print tab*1,"-Name       :",familyName
00151     print tab*1,"-Description:",familyDescription
00152     familyIdentifier = family.getIdentifier()
00153     nbOfAtt = family.getNumberOfAttributes()
00154     print tab*1,"-Identifier :",familyIdentifier
00155     print tab*1,"-Attributes"
00156     attributesids = family.getAttributesIdentifiers()
00157     attributesvals = family.getAttributesValues()
00158     for k in range(nbOfAtt):
00159         print tab*2,"* Id         :",attributesids[k]
00160         print tab*2,"* Value      :",attributesvals[k]
00161         print tab*2,"* Description:",family.getAttributeDescription(k+1)
00162         pass
00163     nbOfGrp = family.getNumberOfGroups()
00164     print tab*1,"-Nb Of Groups:",nbOfGrp
00165     print tab*1,"-Groups:"
00166     for k in range(nbOfGrp):
00167         print tab*2,k+1,":",family.getGroupName(k+1)
00168         pass
00169     _showSupport(family,1,showElems)
00170     return
00171 
00172 ## Dump all families in mesh.
00173 ## Optionally dump nodal connectivity of <showElems> first elements of each family.
00174 ## Use showElems=SHOW_ALL to dump connectivity of all elements.
00175 
00176 def ShowFamilies(mesh, showElems=0):
00177     line = "families in mesh <" + mesh.getName() + ">"
00178     print "\n",line,"\n","-"*len(line)
00179     for entity in [MED_NODE,MED_CELL,MED_FACE,MED_EDGE]:
00180         nbFam = mesh.getNumberOfFamilies(entity)
00181         for i in range(nbFam):
00182             ShowFamily(mesh,entity,i+1,showElems)
00183         pass
00184     print
00185 
00186 ## Dump a GROUP.
00187 ## Optionally dump nodal connectivity of <showElems> first elements of the group.
00188 ## Use showElems=SHOW_ALL to dump connectivity of all elements.
00189 
00190 def ShowGroup(group, showElems):
00191     groupName = group.getName()
00192     groupDescription = group.getDescription()
00193     nbOfFam = group.getNumberOfFamilies()
00194     print "\nGROUP on",theEntityName[group.getEntity()]
00195     print tab*1,"-Name          :",groupName
00196     print tab*1,"-Description   :",groupDescription
00197     print tab*1,"-Nb Of Families:",nbOfFam
00198     print tab*1,"-Families"
00199     for k in range(nbOfFam):
00200         print tab*2,k+1,":",group.getFamily(k+1).getName()
00201         pass
00202     _showSupport(group,1,showElems)
00203     return
00204 
00205 ## Dump all GROUP's in mesh.
00206 ## Optionally dump nodal connectivity of <showElems> first elements of each group.
00207 ## Use showElems=SHOW_ALL to dump connectivity of all elements.
00208 
00209 def ShowGroups(mesh, showElems=0):
00210     line = "groups in mesh <" + mesh.getName() + ">"
00211     print "\n",line,"\n","-"*len(line)
00212     for entity in [MED_NODE,MED_CELL,MED_FACE,MED_EDGE]:
00213         nbGrp = mesh.getNumberOfGroups(entity)
00214         if nbGrp > 0:
00215             for j in range(nbGrp):
00216                 group = mesh.getGroup(entity,j+1)
00217                 ShowGroup(group,showElems)
00218             pass
00219         pass
00220     pass
00221 
00222 ## Dump mesh general information.
00223 ## Optionally dump node coordinates of first <nodes2Show> nodes.
00224 ## <entity2Show> gives number of elements to dump nodal connectivity by
00225 ## entities: [<nb CELLs to show>, <nb FACEs>, <nb EDGEs> ].
00226 ## Use SHOW_ALL to dump all elements or node coordinates.
00227 
00228 def ShowMesh(mesh, nodes2Show=0, entity2Show=[0,0,0]):
00229     print "---------------------- MESH -------------------------"
00230     meshName   = mesh.getName()
00231     spaceDim   = mesh.getSpaceDimension()
00232     meshDim    = mesh.getMeshDimension()
00233     print "The mesh <%s> is a %dD mesh on a %dD geometry" % (meshName,meshDim,spaceDim)
00234     nbNodes    = mesh.getNumberOfNodes()
00235     print "There are",nbNodes,"MED_NODE's"
00236     coordSyst  = mesh.getCoordinatesSystem()
00237     print "The Coordinates :"
00238     coordNames = []
00239     coordUnits = []
00240     for isd in range(spaceDim):
00241         coordNames.append(mesh.getCoordinateName(isd))
00242         coordUnits.append(mesh.getCoordinateUnit(isd))
00243         pass
00244     print tab,"system:",coordSyst
00245     print tab,"names:", coordNames
00246     print tab,"units:", coordUnits
00247     ## coordinates
00248     if nodes2Show:
00249         print tab,"values:"
00250         coordinates = mesh.getCoordinates(MED_FULL_INTERLACE)
00251         nbCoord = nodes2Show
00252         maxNbCoord = len( coordinates ) / spaceDim
00253         if nbCoord < 0: nbCoord = maxNbCoord
00254         else:           nbCoord = min( nbCoord, maxNbCoord )
00255         for k in range( nbCoord ):
00256             n = k*spaceDim
00257             print tab*2,k+1,coordinates[n:n+spaceDim]
00258             pass
00259         if nbCoord < maxNbCoord: print tab*2,"... %d nodes skipped" % (maxNbCoord-nbCoord)
00260         pass
00261     # elem types
00262     print "The Elements :"
00263     i = -1
00264     for entity in [MED_CELL,MED_FACE,MED_EDGE]:
00265         i += 1
00266         entityName = theEntityName[ entity ]
00267         if mesh.getNumberOfElements(entity,MED_ALL_ELEMENTS) < 1: continue
00268         nbTypes = mesh.getNumberOfTypes( entity )
00269         try:
00270             types = mesh.getTypes( entity )
00271         except:
00272             continue
00273         print tab,"%s types:" % entityName
00274         for type in types:
00275             nbElemType = mesh.getNumberOfElements(entity,type)
00276             print tab*2,"%s: \t %d elements" % ( theTypeName[ type ], nbElemType )
00277             pass
00278         # nodal connectivity
00279         if i >= len( entity2Show ): break
00280         if not entity2Show[ i ]: continue
00281         print tab,"%s nodal connectivity:" % entityName
00282         for type in types:
00283             typeName = theTypeName[ type ]
00284             nbElemType = mesh.getNumberOfElements(entity,type)
00285             if nbElemType == 0:
00286                 continue
00287             d = 1
00288             number = range (d, nbElemType+d)
00289             _showNodalConnectivity(mesh,entity,type,number,2,entity2Show[ i ])
00290             pass
00291         pass
00292 
00293     print "----------------------Groups, Families-------------------------"
00294     nbF = 0
00295     nbG = 0
00296     for entity in [MED_NODE,MED_CELL,MED_FACE,MED_EDGE]:
00297         nbFam = mesh.getNumberOfFamilies(entity)
00298         nbGrp = mesh.getNumberOfGroups(entity)
00299         nbElem= mesh.getNumberOfElements(entity, MED_ALL_ELEMENTS);
00300         nbF += nbFam
00301         nbG += nbGrp
00302         if (entity == MED_NODE) :
00303             if (nbFam > 0) : print "This mesh has",nbFam,"Node Family(ies)"
00304             if (nbGrp > 0) : print "This mesh has",nbGrp,"Node Group(s)"
00305             if (nbElem> 0) : print "This mesh has",nbElem,"Node Element(s)"
00306             pass
00307         elif (entity == MED_CELL) :
00308             if (nbFam > 0) : print "This mesh has",nbFam,"Cell Family(ies)"
00309             if (nbGrp > 0) : print "This mesh has",nbGrp,"Cell Group(s)"
00310             if (nbElem> 0) : print "This mesh has",nbElem,"Cell Element(s)"
00311             pass
00312         elif (entity == MED_FACE) :
00313             if (nbFam > 0) : print "This mesh has",nbFam,"Face Family(ies)"
00314             if (nbGrp > 0) : print "This mesh has",nbGrp,"Face Group(s)"
00315             if (nbElem> 0) : print "This mesh has",nbElem,"Face Element(s)"
00316             pass
00317         elif (entity == MED_EDGE) :
00318             if (nbFam > 0) : print "This mesh has",nbFam,"Edge Family(ies)"
00319             if (nbGrp > 0) : print "This mesh has",nbGrp,"Edge Group(s)"
00320             if (nbElem> 0) : print "This mesh has",nbElem,"Edge Element(s)"
00321             pass
00322         pass
00323     print "Total nbF", nbF,"nbG",nbG
00324     return
00325 
00326 ## Dump all FIELD's in MED.
00327 ## Optionally dump <showValues> first values.
00328 ## Use showValues=SHOW_ALL to dump all values.
00329 
00330 def ShowFields( fields, showValues=0 ):
00331     nbFields = len(fields)
00332     print "---------------------- Fields-------------------------"
00333     print "Nb fields", nbFields
00334     for (iField, f ) in enumerate( fields ):
00335         sup       = f.getSupport()
00336         name      = f.getName()
00337         desc      = f.getDescription()
00338         itnb      = f.getIterationNumber()
00339         time      = f.getTime()
00340         order     = f.getOrderNumber()
00341         ftype     = f.getValueType()
00342         mode      = f.getInterlacingType()
00343         nbcomp    = f.getNumberOfComponents()
00344         nbval     = f.getNumberOfValues()
00345         nbelem    = sup.getNumberOfElements(MED_ALL_ELEMENTS)
00346         nbOfTypes = sup.getNumberOfTypes()
00347         types     = sup.getTypes()
00348         isOnAll   = sup.isOnAllElements()
00349         print '\nFIELD',iField
00350         print tab*1,'-Name             : "%s"' % name
00351         print tab*1,'-Description      : "%s"' % desc
00352         print tab*1,'-IterationNumber  :  %s' % itnb
00353         print tab*1,'-Time             :  %s' % time
00354         print tab*1,'-OrderNumber      :  %s' % order
00355         print tab*1,'-Nb Values        :  %s' % nbval
00356         print tab*1,'-Nb Supp. Elements:  %s' % nbelem
00357         print tab*1,'-Nb Componenets   :  %s' % nbcomp
00358         print tab*1,'-ValueType        :  %s' % med_type_champ[ftype]
00359         print tab*1,'-Interlace        :  %s' % medModeSwitch[mode]
00360         print tab*1,'-Conponents'
00361         for k in range(nbcomp):
00362             kp1 = k+1
00363             compName = f.getComponentName(kp1)
00364             compDesc = f.getComponentDescription(kp1)
00365             compUnit = f.getMEDComponentUnit(kp1)
00366             print     tab*2,kp1,'*Name       : "%s"' % compName
00367             try:
00368                 print tab*2,'  *Description: "%s"' % compDesc
00369             except:
00370                 print 'INVALID'
00371                 pass
00372             try:
00373                 print tab*2,'  *Unit       : "%s"' % compUnit
00374             except:
00375                 print 'INVALID'
00376                 pass
00377             pass
00378         print tab*1,'-MESH             : "%s"' % sup.getMeshName()
00379         print tab*1,'-SUPPORT          : "%s"' % sup.getName()
00380         print tab*1,'-On all elements  :  %s' % bool(isOnAll)
00381         print tab*1,'-Types            :  %s'  % types
00382 
00383         if ftype == MED_REEL64:
00384             if mode == MED_FULL_INTERLACE:
00385                 f = createFieldDoubleFromField(f)
00386             else:
00387                 f = createFieldDoubleNoInterlaceFromField( f )
00388                 pass
00389             pass
00390         elif ftype == MED_INT32:
00391             if mode == MED_FULL_INTERLACE:
00392                 f = createFieldIntFromField(f)
00393             else:
00394                 f = createFieldIntNoInterlaceFromField( f )
00395                 pass
00396             pass
00397         else:
00398             print tab*1,'<< Unknown field type >>:',ftype
00399             continue
00400         nbGauss = 1
00401         hasGauss = False
00402         if nbcomp == 0:
00403             nbGauss = 0
00404         else:
00405             hasGauss = f.getGaussPresence()
00406             pass
00407         if hasGauss:
00408             nbGaussByType = f.getNumberOfGaussPoints()
00409             pass
00410         for k in range(nbOfTypes):
00411             type = types[k]
00412             nbOfElmtsOfType = sup.getNumberOfElements(type)
00413             if hasGauss: nbGauss = nbGaussByType[ k ]
00414             if type == 0: type = MED_POINT1
00415             print tab*2,k+1,theTypeName[type],':',nbOfElmtsOfType, 'elements,',\
00416                   nbGauss,'gauss point(s)'
00417             pass
00418         nbOf = sup.getNumberOfElements(MED_ALL_ELEMENTS)
00419         elements = []
00420         if not isOnAll:
00421             elements = sup.getNumber(MED_ALL_ELEMENTS)
00422             pass
00423         if nbcomp == 0:
00424             nbOf = 0
00425             print tab*1,'-Nb Values        :',nbOf
00426             #value = f.getValue(MED_FULL_INTERLACE)
00427             #print value[0: min( 100, len(value)-1 )]
00428 
00429         toShow = min( nbOf, showValues )
00430         if toShow < 0: toShow = nbOf
00431         for I in range( toShow ):
00432             if elements:
00433                 i = elements[ I ]
00434             else:
00435                 i = I+1
00436             if mode == MED_FULL_INTERLACE:
00437                 valueI = f.getRow(i)
00438             else:
00439                 valueI = []
00440                 for j in range( nbcomp ):
00441                     for k in range( f.getNbGaussI( i ) ):
00442                         valueI.append( f.getValueIJK(i,j+1,k+1) )
00443             print '         ',i,' - ',valueI #[:nbcomp]
00444             pass
00445         if nbOf > toShow:
00446             print '            ...skip',nbOf - toShow,'values'
00447             pass
00448         pass
00449     pass
00450 
00451 ## Read all fields in MED
00452 
00453 def ReadFields(med):
00454     nbFields = med.getNumberOfFields()
00455     if (nbFields>0):
00456         print 'READ FIELDs ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
00457         med.updateSupport()
00458         for i in range(nbFields):
00459             field_name = med.getFieldName(i)
00460             nbOfIt = med.getFieldNumberOfIteration(field_name)
00461             print '  The field is',field_name,'with',nbOfIt,'iteration(s)'
00462             for j in range(nbOfIt):
00463                 dtitfield = med.getFieldIteration(field_name,j)
00464                 dt = dtitfield.getdt()
00465                 it = dtitfield.getit()
00466                 field = med.getField(field_name,dt,it)
00467                 type = field.getValueType()
00468                 print '     * Iteration:',dt,'Order number:',it,'Type:',type
00469                 mode = field.getInterlacingType()
00470                 if type == MED_INT32:
00471                     if mode == MED_FULL_INTERLACE:
00472                         fieldint = createFieldIntFromField(field)
00473                     else:
00474                         fieldint = createFieldIntNoInterlaceFromField( field )
00475                     print '     Reading',fieldint.getName(),'...'
00476                     fieldint.read()
00477                 elif type == MED_REEL64:
00478                     if mode == MED_FULL_INTERLACE:
00479                         f = createFieldDoubleFromField(field)
00480                     else:
00481                         f = createFieldDoubleNoInterlaceFromField( field )
00482                     print '     Reading',f.getName(),'...'
00483                     f.read()
00484                 else:
00485                     print '  !!!! Bad type of Field !!!!'
00486 
00487 # Remove a file if it exists
00488 
00489 def supprimer(f):
00490     if os.access(f, os.F_OK):
00491         os.remove(f)
00492 
00493 # Remove a file if it exists
00494 
00495 def deleteFile( f ):
00496     supprimer( f )