src/MEDMEM_SWIG/dumpMEDMEM.py

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