src/MEDMEM_SWIG/dumpMEDMEM.py

   1 # Copyright (C) 2005  OPEN CASCADE, CEA, EDF R&D, LEG
   2 #           PRINCIPIA R&D, EADS CCR, Lip6, BV, CEDRAT
   3 # This library is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU Lesser General Public
   5 # License as published by the Free Software Foundation; either
   6 # version 2.1 of the License.
   7 #
   8 # This library is distributed in the hope that it will be useful
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 # Lesser General Public License for more details.
  12 #
  13 # You should have received a copy of the GNU Lesser General Public
  14 # License along with this library; if not, write to the Free Software
  15 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 #
  17 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 #
  19 ############################################################################
  20 #
  21 # This file provides utilities to dump content of MEDMEM objects: mesh,
  22 # field, group, family, nodal connectivity, node coordinates.
  23 #
  24 ############################################################################
  25
  26 from libMEDMEM_Swig import *
  27
  28 import os
  29
  30 theEntityName = { MED_CELL        :"CELL",
  31                   MED_FACE        :"FACE",
  32                   MED_EDGE        :"EDGE",
  33                   MED_NODE        :"NODE",
  34                   MED_ALL_ENTITIES:"ALL_ENTITIES" }
  35
  36 theTypeName = {MED_NONE        :"NONE",
  37                MED_POINT1      :"POINT1",
  38                MED_SEG2        :"SEG2",
  39                MED_SEG3        :"SEG3",
  40                MED_TRIA3       :"TRIA3",
  41                MED_QUAD4       :"QUAD4",
  42                MED_TRIA6       :"TRIA6",
  43                MED_QUAD8       :"QUAD8",
  44                MED_TETRA4      :"TETRA4",
  45                MED_PYRA5       :"PYRA5",
  46                MED_PENTA6      :"PENTA6",
  47                MED_HEXA8       :"HEXA8",
  48                MED_TETRA10     :"TETRA10",
  49                MED_PYRA13      :"PYRA13",
  50                MED_PENTA15     :"PENTA15",
  51                MED_HEXA20      :"HEXA20",
  52                MED_POLYGON     :"POLYGON",
  53                MED_POLYHEDRA   :"POLYHEDRA",
  54                MED_ALL_ELEMENTS:"ALL_ELEMENTS"}
  55
  56 medModeSwitch = { 0:"FULL_INTERLACE",
  57                   1:"NO_INTERLACE",
  58                   3:"UNDEFINED_INTERLACE" };
  59
  60 med_type_champ = { 6 : "REEL64",
  61                    24: "INT32",
  62                    26: "INT64",
  63                    0 : "UNDEFINED_TYPE" } ;
  64
  65 tab="   "
  66
  67 debugShowConn=True
  68 debugShowConn=False
  69
  70 SHOW_ALL = -1
  71
  72 # private
  73 def _showNodalConnectivity(mesh,entity,type,elems,tablevel,showOnly=SHOW_ALL):
  74     if debugShowConn: print "ELEMENS:",elems
  75     tab1 = tab*tablevel
  76     tab2 = tab*(tablevel+1)
  77     typeName = theTypeName[type]
  78     nbElem = len( elems )
  79     if showOnly > 0:
  80         elems = elems[:showOnly]
  81     nbShow = len( elems )
  82     if type == MED_POLYGON:
  83         connectivity = mesh.getPolygonsConnectivity(MED_NODAL,entity)
  84         index = mesh.getPolygonsConnectivityIndex(MED_NODAL,entity)
  85         if debugShowConn: print "CONN:",connectivity,"\nIND:",index
  86         d = mesh.getNumberOfElements(entity,MED_ALL_ELEMENTS)
  87         for i in elems:
  88             j = i - d
  89             print tab1,typeName,i,":",connectivity[ index[j-1]-1 : index[j]-1 ]
  90             pass
  91         pass
  92     elif type == MED_POLYHEDRA:
  93         connectivity = mesh.getPolyhedronConnectivity(MED_NODAL)
  94         fIndex = mesh.getPolyhedronFacesIndex()
  95         index = mesh.getPolyhedronIndex(MED_NODAL)
  96         d = mesh.getNumberOfElements(entity,MED_ALL_ELEMENTS)
  97         if debugShowConn: print "CONN:",connectivity,"\nIND:",index,"\nFIND:",fIndex
  98         for i in elems:
  99             j = i - d
 100             print tab1,typeName, i
 101             iF1, iF2 = index[ j-1 ]-1, index[ j ]-1
 102             for f in range( iF2 - iF1 ):
 103                 iN1, iN2 = fIndex[ iF1+f ]-1, fIndex[ iF1+f+1 ]-1
 104                 print tab2,"Face",(f+1),":",connectivity[ iN1 : iN2 ]
 105                 pass
 106             pass
 107         pass
 108     else:
 109         elemShift = 0
 110         types = mesh.getTypes( entity )
 111         for t in types:
 112             if t != type:
 113                 elemShift += mesh.getNumberOfElements(entity,t)
 114             else:
 115                 break
 116             pass
 117         connectivity = mesh.getConnectivity(MED_FULL_INTERLACE, MED_NODAL, entity,MED_ALL_ELEMENTS)
 118         index = mesh.getConnectivityIndex(MED_FULL_INTERLACE, entity)
 119         if debugShowConn: print "CONN: %s\n INX: %s" % (connectivity,index)
 120         nbNodesPerCell = type%100
 121         for i in elems:
 122             elem = i + elemShift
 123             n = index[elem-1]-1
 124             print tab1,typeName,i,":",connectivity[n:n+nbNodesPerCell]
 125             pass
 126         pass
 127     nbSkip = nbElem - nbShow
 128     if nbSkip > 0:
 129         print tab1,"...",nbSkip,"elements not shown"
 130         pass
 131     pass
 132
 133 #private
 134 def _showSupport(support, tablevel,showElems=0):
 135     tab1 = tab*(tablevel+0)
 136     tab2 = tab*(tablevel+1)
 137     tab3 = tab*(tablevel+3)
 138     entity    = support.getEntity()
 139     types     = support.getTypes()
 140     nbOfTypes = support.getNumberOfTypes()
 141     onAll     = support.isOnAllElements()
 142     print tab1,"-Entity:",theEntityName[entity]
 143     print tab1,"-Types :",types
 144     print tab1,"-Elements"
 145     if onAll:
 146         print tab2,"<< Is on all elements >>"
 147     else:
 148         for type in types:
 149             nbOfElmtsOfType = support.getNumberOfElements(type)
 150             number = support.getNumber(type)
 151             number.sort()
 152             print tab2,"* Type:",theTypeName[type]
 153             print     tab3,". Nb elements:",nbOfElmtsOfType
 154             print     tab3,". Numbers:",number[:min(100,nbOfElmtsOfType)],
 155             if nbOfElmtsOfType > 100:
 156                 print "...skip", nbOfElmtsOfType-100
 157             else:
 158                 print
 159                 pass
 160             if entity != MED_NODE and showElems:
 161                 print tab3,". Nodal connectivity"
 162                 _showNodalConnectivity(support.getMesh(),entity,type,number,tablevel+4,showElems)
 163                 pass
 164             pass
 165         pass
 166     print
 167     return
 168
 169 ## Dump i-th family of given entity in mesh.
 170 ## Optionally dump nodal connectivity of <showElems> first elements.
 171 ## Use showElems=SHOW_ALL to dump connectivity of all elements.
 172
 173 def ShowFamily(mesh, entity, i, showElems=0):
 174     family = mesh.getFamily(entity,i)
 175     familyName = family.getName()
 176     familyDescription = family.getDescription()
 177     entity = family.getEntity()
 178     familyBool = family.isOnAllElements()
 179     print "\nFAMILY", i, "on", theEntityName[entity]
 180     print tab*1,"-Name       :",familyName
 181     print tab*1,"-Description:",familyDescription
 182     familyIdentifier = family.getIdentifier()
 183     nbOfAtt = family.getNumberOfAttributes()
 184     print tab*1,"-Identifier :",familyIdentifier
 185     print tab*1,"-Attributes"
 186     attributesids = family.getAttributesIdentifiers()
 187     attributesvals = family.getAttributesValues()
 188     for k in range(nbOfAtt):
 189         print tab*2,"* Id         :",attributesids[k]
 190         print tab*2,"* Value      :",attributesvals[k]
 191         print tab*2,"* Description:",family.getAttributeDescription(k+1)
 192         pass
 193     nbOfGrp = family.getNumberOfGroups()
 194     print tab*1,"-Nb Of Groups:",nbOfGrp
 195     print tab*1,"-Groups:"
 196     for k in range(nbOfGrp):
 197         print tab*2,k+1,":",family.getGroupName(k+1)
 198         pass
 199     _showSupport(family,1,showElems)
 200     return
 201
 202 ## Dump all families in mesh.
 203 ## Optionally dump nodal connectivity of <showElems> first elements of each family.
 204 ## Use showElems=SHOW_ALL to dump connectivity of all elements.
 205
 206 def ShowFamilies(mesh, showElems=0):
 207     line = "families in mesh <" + mesh.getName() + ">"
 208     print "\n",line,"\n","-"*len(line)
 209     for entity in [MED_NODE,MED_CELL,MED_FACE,MED_EDGE]:
 210         nbFam = mesh.getNumberOfFamilies(entity)
 211         for i in range(nbFam):
 212             ShowFamily(mesh,entity,i+1,showElems)
 213         pass
 214     print
 215
 216 ## Dump a GROUP.
 217 ## Optionally dump nodal connectivity of <showElems> first elements of the group.
 218 ## Use showElems=SHOW_ALL to dump connectivity of all elements.
 219
 220 def ShowGroup(group, showElems):
 221     groupName = group.getName()
 222     groupDescription = group.getDescription()
 223     nbOfFam = group.getNumberOfFamilies()
 224     print "\nGROUP on",theEntityName[group.getEntity()]
 225     print tab*1,"-Name          :",groupName
 226     print tab*1,"-Description   :",groupDescription
 227     print tab*1,"-Nb Of Families:",nbOfFam
 228     print tab*1,"-Families"
 229     for k in range(nbOfFam):
 230         print tab*2,k+1,":",group.getFamily(k+1).getName()
 231         pass
 232     _showSupport(group,1,showElems)
 233     return
 234
 235 ## Dump all GROUP's in mesh.
 236 ## Optionally dump nodal connectivity of <showElems> first elements of each group.
 237 ## Use showElems=SHOW_ALL to dump connectivity of all elements.
 238
 239 def ShowGroups(mesh, showElems=0):
 240     line = "groups in mesh <" + mesh.getName() + ">"
 241     print "\n",line,"\n","-"*len(line)
 242     for entity in [MED_NODE,MED_CELL,MED_FACE,MED_EDGE]:
 243         nbGrp = mesh.getNumberOfGroups(entity)
 244         if nbGrp > 0:
 245             for j in range(nbGrp):
 246                 group = mesh.getGroup(entity,j+1)
 247                 ShowGroup(group,showElems)
 248             pass
 249         pass
 250     pass
 251
 252 ## Dump mesh general information.
 253 ## Optionally dump node coordinates of first <nodes2Show> nodes.
 254 ## <entity2Show> gives number of elements to dump nodal connectivity by
 255 ## entities: [<nb CELLs to show>, <nb FACEs>, <nb EDGEs> ].
 256 ## Use SHOW_ALL to dump all elements or node coordinates.
 257
 258 def ShowMesh(mesh, nodes2Show=0, entity2Show=[0,0,0]):
 259     print "---------------------- MESH -------------------------"
 260     meshName   = mesh.getName()
 261     spaceDim   = mesh.getSpaceDimension()
 262     meshDim    = mesh.getMeshDimension()
 263     print "The mesh <%s> is a %dD mesh on a %dD geometry" % (meshName,meshDim,spaceDim)
 264     nbNodes    = mesh.getNumberOfNodes()
 265     print "There are",nbNodes,"MED_NODE's"
 266     coordSyst  = mesh.getCoordinatesSystem()
 267     print "The Coordinates :"
 268     coordNames = []
 269     coordUnits = []
 270     for isd in range(spaceDim):
 271         coordNames.append(mesh.getCoordinateName(isd))
 272         coordUnits.append(mesh.getCoordinateUnit(isd))
 273         pass
 274     print tab,"system:",coordSyst
 275     print tab,"names:", coordNames
 276     print tab,"units:", coordUnits
 277     ## coordinates
 278     if nodes2Show:
 279         print tab,"values:"
 280         coordinates = mesh.getCoordinates(MED_FULL_INTERLACE)
 281         nbCoord = nodes2Show
 282         maxNbCoord = len( coordinates ) / spaceDim
 283         if nbCoord < 0: nbCoord = maxNbCoord
 284         else:           nbCoord = min( nbCoord, maxNbCoord )
 285         for k in range( nbCoord ):
 286             n = k*spaceDim
 287             print tab*2,k+1,coordinates[n:n+spaceDim]
 288             pass
 289         if nbCoord < maxNbCoord: print tab*2,"... %d nodes skipped" % (maxNbCoord-nbCoord)
 290         pass
 291     # elem types
 292     print "The Elements :"
 293     i = -1
 294     for entity in [MED_CELL,MED_FACE,MED_EDGE]:
 295         i += 1
 296         entityName = theEntityName[ entity ]
 297         nbTypes = mesh.getNumberOfTypesWithPoly( entity )
 298         if nbTypes == 0 : continue
 299         try:
 300             types = mesh.getTypesWithPoly( entity )
 301         except:
 302             continue
 303         print tab,"%s types:" % entityName
 304         for type in types:
 305             nbElemType = mesh.getNumberOfElementsWithPoly(entity,type)
 306             print tab*2,"%s: \t %d elements" % ( theTypeName[ type ], nbElemType )
 307             pass
 308         # nodal connectivity
 309         if i >= len( entity2Show ): break
 310         if not entity2Show[ i ]: continue
 311         print tab,"%s nodal connectivity:" % entityName
 312         for type in types:
 313             typeName = theTypeName[ type ]
 314             nbElemType = mesh.getNumberOfElementsWithPoly(entity,type)
 315             if nbElemType == 0:
 316                 continue
 317             d = 1
 318             if type==MED_POLYGON or type==MED_POLYHEDRA:
 319                 d += mesh.getNumberOfElements(entity,MED_ALL_ELEMENTS)
 320             number = range (d, nbElemType+d)
 321             _showNodalConnectivity(mesh,entity,type,number,2,entity2Show[ i ])
 322             pass
 323         pass
 324
 325     print "----------------------Groups, Families-------------------------"
 326     nbF = 0
 327     nbG = 0
 328     for entity in [MED_NODE,MED_CELL,MED_FACE,MED_EDGE]:
 329         nbFam = mesh.getNumberOfFamilies(entity)
 330         nbGrp = mesh.getNumberOfGroups(entity)
 331         nbElem= mesh.getNumberOfElements(entity, MED_ALL_ELEMENTS);
 332         nbF += nbFam
 333         nbG += nbGrp
 334         if (entity == MED_NODE) :
 335             if (nbFam > 0) : print "This mesh has",nbFam,"Node Family(ies)"
 336             if (nbGrp > 0) : print "This mesh has",nbGrp,"Node Group(s)"
 337             if (nbElem> 0) : print "This mesh has",nbElem,"Node Element(s)"
 338             pass
 339         elif (entity == MED_CELL) :
 340             if (nbFam > 0) : print "This mesh has",nbFam,"Cell Family(ies)"
 341             if (nbGrp > 0) : print "This mesh has",nbGrp,"Cell Group(s)"
 342             if (nbElem> 0) : print "This mesh has",nbElem,"Cell Element(s)"
 343             pass
 344         elif (entity == MED_FACE) :
 345             if (nbFam > 0) : print "This mesh has",nbFam,"Face Family(ies)"
 346             if (nbGrp > 0) : print "This mesh has",nbGrp,"Face Group(s)"
 347             if (nbElem> 0) : print "This mesh has",nbElem,"Face Element(s)"
 348             pass
 349         elif (entity == MED_EDGE) :
 350             if (nbFam > 0) : print "This mesh has",nbFam,"Edge Family(ies)"
 351             if (nbGrp > 0) : print "This mesh has",nbGrp,"Edge Group(s)"
 352             if (nbElem> 0) : print "This mesh has",nbElem,"Edge Element(s)"
 353             pass
 354         pass
 355     print "Total nbF", nbF,"nbG",nbG
 356
 357 ## Dump all FIELD's in MED.
 358 ## Optionally dump <showValues> first values.
 359 ## Use showValues=SHOW_ALL to dump all values.
 360
 361 def ShowFields( med, showValues=0 ):
 362     nbFields = med.getNumberOfFields()
 363     print "---------------------- Fields-------------------------"
 364     print "Nb fields", nbFields
 365     for iField in range( med.getNumberOfFields() ):
 366         fName = med.getFieldName( iField )
 367         for it in range( med.getFieldNumberOfIteration( fName ) ):
 368             dt_it_    = med.getFieldIteration( fName, it )
 369             f         = med.getField( fName, dt_it_.dt, dt_it_.it )
 370             sup       = f.getSupport()
 371             name      = f.getName()
 372             desc      = f.getDescription()
 373             itnb      = f.getIterationNumber()
 374             time      = f.getTime()
 375             order     = f.getOrderNumber()
 376             ftype     = f.getValueType()
 377             mode      = f.getInterlacingType()
 378             nbcomp    = f.getNumberOfComponents()
 379             nbval     = f.getNumberOfValues()
 380             nbelem    = sup.getNumberOfElements(MED_ALL_ELEMENTS)
 381             nbOfTypes = sup.getNumberOfTypes()
 382             types     = sup.getTypes()
 383             isOnAll   = sup.isOnAllElements()
 384             print '\nFIELD',iField
 385             print tab*1,'-Name             : "%s"' % name
 386             print tab*1,'-Description      : "%s"' % desc
 387             print tab*1,'-IterationNumber  :  %s' % itnb
 388             print tab*1,'-Time             :  %s' % time
 389             print tab*1,'-OrderNumber      :  %s' % order
 390             print tab*1,'-Nb Values        :  %s' % nbval
 391             print tab*1,'-Nb Supp. Elements:  %s' % nbelem
 392             print tab*1,'-Nb Componenets   :  %s' % nbcomp
 393             print tab*1,'-ValueType        :  %s' % med_type_champ[ftype]
 394             print tab*1,'-Interlace        :  %s' % medModeSwitch[mode]
 395             print tab*1,'-Conponents'
 396             for k in range(nbcomp):
 397                 kp1 = k+1
 398                 compName = f.getComponentName(kp1)
 399                 compDesc = f.getComponentDescription(kp1)
 400                 compUnit = f.getMEDComponentUnit(kp1)
 401                 print     tab*2,kp1,'*Name       : "%s"' % compName
 402                 try:
 403                     print tab*2,'  *Description: "%s"' % compDesc
 404                 except:
 405                     print 'INVALID'
 406                 try:
 407                     print tab*2,'  *Unit       : "%s"' % compUnit
 408                 except:
 409                     print 'INVALID'
 410                 pass
 411             print tab*1,'-SUPPORT          : "%s"' % sup.getName()
 412             print tab*1,'-On all elements  :  %s' % bool(isOnAll)
 413             print tab*1,'-Types            :  %s'  % types
 414
 415             if ftype == MED_REEL64:
 416                 if mode == MED_FULL_INTERLACE:
 417                     f = createFieldDoubleFromField(f)
 418                 else:
 419                     f = createFieldDoubleNoInterlaceFromField( f )
 420             elif ftype == MED_INT32:
 421                 if mode == MED_FULL_INTERLACE:
 422                     f = createFieldIntFromField(f)
 423                 else:
 424                     f = createFieldIntNoInterlaceFromField( f )
 425             else:
 426                 print tab*1,'<< Unknown field type >>:',ftype
 427                 continue
 428             nbGauss = 1
 429             hasGauss = False
 430             if nbcomp == 0:
 431                 nbGauss = 0
 432             else:
 433                 hasGauss = f.getGaussPresence()
 434             if hasGauss:
 435                 nbGaussByType = f.getNumberOfGaussPoints()
 436             for k in range(nbOfTypes):
 437                 type = types[k]
 438                 nbOfElmtsOfType = sup.getNumberOfElements(type)
 439                 if hasGauss: nbGauss = nbGaussByType[ k ]
 440                 if type == 0: type = MED_POINT1
 441                 print tab*2,k+1,theTypeName[type],':',nbOfElmtsOfType, 'elements,',\
 442                       nbGauss,'gauss point(s)'
 443                 pass
 444             nbOf = sup.getNumberOfElements(MED_ALL_ELEMENTS)
 445             elements = []
 446             if not isOnAll:
 447                 elements = sup.getNumber(MED_ALL_ELEMENTS)
 448             if nbcomp == 0:
 449                 nbOf = 0
 450             print tab*1,'-Nb Values        :',nbOf
 451             #value = f.getValue(MED_FULL_INTERLACE)
 452             #print value[0: min( 100, len(value)-1 )]
 453
 454             toShow = min( nbOf, showValues )
 455             if toShow < 0: toShow = nbOf
 456             for I in range( toShow ):
 457                 if elements:
 458                     i = elements[ I ]
 459                 else:
 460                     i = I+1
 461                 if mode == MED_FULL_INTERLACE:
 462                     valueI = f.getRow(i)
 463                 else:
 464                     valueI = []
 465                     for j in range( nbcomp ):
 466                         for k in range( f.getNbGaussI( i ) ):
 467                             valueI.append( f.getValueIJK(i,j+1,k+1) )
 468                 print '         ',i,' - ',valueI #[:nbcomp]
 469                 pass
 470             if nbOf > toShow:
 471                 print '            ...skip',nbOf - toShow,'values'
 472                 pass
 473             pass
 474         pass
 475     pass
 476
 477 ## Read all fields in MED
 478
 479 def ReadFields(med):
 480     nbFields = med.getNumberOfFields()
 481     if (nbFields>0):
 482         print 'READ FIELDs ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
 483         med.updateSupport()
 484         for i in range(nbFields):
 485             field_name = med.getFieldName(i)
 486             nbOfIt = med.getFieldNumberOfIteration(field_name)
 487             print '  The field is',field_name,'with',nbOfIt,'iteration(s)'
 488             for j in range(nbOfIt):
 489                 dtitfield = med.getFieldIteration(field_name,j)
 490                 dt = dtitfield.getdt()
 491                 it = dtitfield.getit()
 492                 field = med.getField(field_name,dt,it)
 493                 type = field.getValueType()
 494                 print '     * Iteration:',dt,'Order number:',it,'Type:',type
 495                 mode = field.getInterlacingType()
 496                 if type == MED_INT32:
 497                     if mode == MED_FULL_INTERLACE:
 498                         fieldint = createFieldIntFromField(field)
 499                     else:
 500                         fieldint = createFieldIntNoInterlaceFromField( field )
 501                     print '     Reading',fieldint.getName(),'...'
 502                     fieldint.read()
 503                 elif type == MED_REEL64:
 504                     if mode == MED_FULL_INTERLACE:
 505                         f = createFieldDoubleFromField(field)
 506                     else:
 507                         f = createFieldDoubleNoInterlaceFromField( field )
 508                     print '     Reading',f.getName(),'...'
 509                     f.read()
 510                 else:
 511                     print '  !!!! Bad type of Field !!!!'
 512
 513 # Remove a file if it exists
 514
 515 def supprimer(f):
 516     if os.access(f, os.F_OK):
 517         os.remove(f)
 518
 519 # Remove a file if it exists
 520
 521 def deleteFile( f ):
 522     supprimer( f )