src/MEDMEM_SWIG/med_test2.py

   1 #  -*- coding: iso-8859-1 -*-
   2 # Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   3 #
   4 # Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   5 # CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   6 #
   7 # This library is free software; you can redistribute it and/or
   8 # modify it under the terms of the GNU Lesser General Public
   9 # License as published by the Free Software Foundation; either
  10 # version 2.1 of the License.
  11 #
  12 # This library is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15 # Lesser General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU Lesser General Public
  18 # License along with this library; if not, write to the Free Software
  19 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  20 #
  21 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  22 #
  23
  24 ###################################################################################
  25 # This Python script is parsing a MED file using MED Memory from SALOME platform:
  26 # It analyses all meshes in the MED file (coordinates, connectivity of d-cells as
  27 # well as (d-1)-cells, families), it tests all fields generated in the MESH class
  28 # and write them in a 2 new med files (V2.1 and V2.2), it gives only the number of
  29 # fields stored in the MED file (d is the mesh/space dimension).
  30 ###################################################################################
  31 #
  32 from libMEDMEM_Swig import *
  33 from random import *
  34
  35 #==============================================================================
  36
  37 def AnalyzeField(field):
  38     name = field.getName()
  39     desc = field.getDescription()
  40     nbComp = field.getNumberOfComponents()
  41     itNum = field.getIterationNumber()
  42     ordNum = field.getOrderNumber()
  43     type = field.getValueType()
  44
  45     print "Analysis of the field ",name," with the description ",desc
  46     print "iteration number ",itNum," order Number ",ordNum
  47     print "It has ",nbComp," component(s) with the type ",type
  48
  49     fieldValue = field.getValue()
  50     fieldSupport = field.getSupport()
  51     fieldMesh = fieldSupport.getMesh()
  52     fieldEntity = fieldSupport.getEntity()
  53     bool = fieldSupport.isOnAllElements()
  54
  55     if bool:
  56         print "The support of this field is on all entities ",fieldEntity," of the mesh ",fieldMesh.getName()
  57         if fieldEntity == MED_NODE:
  58             nbValByComp = fieldMesh.getNumberOfNodes()
  59         else:
  60             nbValByComp = fieldMesh.getNumberOfElements(fieldEntity,MED_ALL_ELEMENTS)
  61         print "and its dimension (number of values by component of the field) is ",nbValByComp
  62     else:
  63         print "The support of this field is partially on entities ",fieldEntity," of the mesh ",fieldMesh.getName()
  64         nbValByComp = fieldSupport.getNumberOfElements(MED_ALL_ELEMENTS)
  65         print "and its dimension (number of values by component of the field) is ",nbValByComp
  66
  67     for i in range(nbComp):
  68         ip1 = i + 1
  69         compName = field.getComponentName(ip1)
  70         compDesc = field.getComponentDescription(ip1)
  71         compUnit = field.getMEDComponentUnit(ip1)
  72         print "The ",(i+1),"-th  component ",compName," with the dexription ",compDesc," and the unit ",compUnit
  73
  74     for i in range(nbValByComp):
  75         print "  * ",fieldValue[i*nbComp:(i+1)*nbComp]
  76
  77 #==============================================================================
  78 import os
  79 #
  80 #befor running this script, please be sure about the path the file fileName
  81 #
  82 filePath=os.environ["MED_ROOT_DIR"]
  83 filePath=os.path.join(filePath, "share", "salome", "resources", "med")
  84
  85 medFile = os.path.join(filePath, "carre_en_quad4_seg2.med")
  86 #medFile = os.path.join(filePath, "cube_hexa8_quad4.med")
  87
  88 def print_ord(i):
  89     if i == 0:
  90         return 'first'
  91     elif i == 1:
  92         return 'second'
  93     elif i == 2:
  94         return 'third'
  95     else:
  96         return `i`+'th'
  97
  98 md = MEDFILEBROWSER(medFile)
  99
 100 nbMeshes = md.getNumberOfMeshes()
 101
 102 nbFields = md.getNumberOfFields()
 103
 104 print "The med file", medFile, "contains", nbMeshes, "mesh(es) and", nbFields, "field(s)"
 105
 106 if (nbMeshes>0):
 107     print "Mesh(es) Name(s) is(are) "
 108
 109     for i in range(nbMeshes):
 110         mesh_name = md.getMeshName(i)
 111         print "   - ",mesh_name
 112
 113 if (nbFields>0):
 114     print "Field(s) Name(s) is(are) "
 115
 116     for i in range(nbFields):
 117         field_name = md.getFieldName(i)
 118         print "   - ",field_name
 119
 120 print ""
 121
 122 if (nbMeshes>0):
 123     print "Mesh(es) Analysis "
 124     for i in range(nbMeshes):
 125         mesh_name = md.getMeshName(i)
 126         mesh = MESH(MED_DRIVER,md.getFileName(),mesh_name)
 127         spaceDim = mesh.getSpaceDimension()
 128         meshDim = mesh.getMeshDimension()
 129         print "The",print_ord(i), "mesh, '",mesh_name,"', is a",spaceDim,"D mesh on a",meshDim,"D geometry"
 130         nbNodes = mesh.getNumberOfNodes()
 131         print "The mesh has",nbNodes,"Nodes"
 132         coordSyst = mesh.getCoordinatesSystem()
 133         print "The coordinates system is",coordSyst
 134         print "The Coordinates :"
 135         coordNames = []
 136         coordUnits = []
 137         for isd in range(spaceDim):
 138             coordNames.append(mesh.getCoordinateName(isd))
 139             coordUnits.append(mesh.getCoordinateUnit(isd))
 140
 141         print "names:", coordNames
 142         print "units", coordUnits
 143         print "values:"
 144         coordinates = mesh.getCoordinates(MED_FULL_INTERLACE)
 145         for k in range(nbNodes):
 146             kp1 = k+1
 147             coords = []
 148             for isd in range(spaceDim):
 149                 isdp1 = isd+1
 150                 coords.append(mesh.getCoordinate(kp1,isdp1))
 151
 152             print coords," ---- ", coordinates[k*spaceDim:((k+1)*spaceDim)]
 153
 154         print ""
 155         print "Show the Nodal Connectivity:"
 156         nbTypesCell = mesh.getNumberOfTypes(MED_CELL)
 157         print ""
 158         if (nbTypesCell>0):
 159             print "The Mesh has",nbTypesCell,"Type(s) of Cell"
 160             types = mesh.getTypes(MED_CELL)
 161             for k in range(nbTypesCell):
 162                 type = types[k]
 163                 nbElemType = mesh.getNumberOfElements(MED_CELL,type)
 164                 print "For the type:",type,"there is(are)",nbElemType,"elemnt(s)"
 165                 connectivity = mesh.getConnectivity(MED_NODAL,MED_CELL,type)
 166                 nbNodesPerCell = type%100
 167                 for j in range(nbElemType):
 168                     print "Element",(j+1)," ",connectivity[j*nbNodesPerCell:(j+1)*nbNodesPerCell]
 169
 170         print ""
 171         print "Show the Reverse Nodal Connectivity:"
 172         ReverseConnectivity = mesh.getReverseConnectivity(MED_NODAL)
 173         ReverseConnectivityIndex = mesh.getReverseConnectivityIndex(MED_NODAL)
 174         print ""
 175         for j in range(nbNodes):
 176             begin = ReverseConnectivityIndex[j]-1
 177             end = ReverseConnectivityIndex[j+1]-1
 178             print "Node",(j+1),"-->",ReverseConnectivity[begin:end]
 179
 180         print ""
 181         print "Show the Descending Connectivity:"
 182         mesh.calculateConnectivity(MED_DESCENDING,MED_CELL)
 183         nbElemts = mesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)
 184         Connectivity = mesh.getConnectivity(MED_DESCENDING,MED_CELL,MED_ALL_ELEMENTS)
 185         ConnectivityIndex = mesh.getConnectivityIndex(MED_DESCENDING,MED_CELL)
 186         print ""
 187         for j in range(nbElemts):
 188             begin = ConnectivityIndex[j]-1
 189             end = ConnectivityIndex[j+1]-1
 190             print "Element",(j+1),"-->",Connectivity[begin:end]
 191
 192         print ""
 193         for entity in [MED_NODE,MED_CELL,MED_FACE,MED_EDGE]:
 194             nbFam = mesh.getNumberOfFamilies(entity)
 195             if (entity == MED_NODE) & (nbFam > 0):
 196                 print "This mesh has",nbFam,"Node Family(ies)"
 197             elif (entity == MED_CELL) & (nbFam > 0):
 198                 print "This mesh has",nbFam,"Cell Family(ies)"
 199             elif (entity == MED_FACE) & (nbFam > 0):
 200                 print "This mesh has",nbFam,"Face Family(ies)"
 201             elif (entity == MED_EDGE) & (nbFam > 0):
 202                 print "This mesh has",nbFam,"Edge Family(ies)"
 203
 204             if nbFam > 0:
 205                 for j in range(nbFam):
 206                     print ""
 207                     family = mesh.getFamily(entity,j+1)
 208                     familyName = family.getName()
 209                     familyDescription = family.getDescription()
 210                     familyEntity = family.getEntity()
 211                     familyBool = family.isOnAllElements()
 212                     print "  -Name:",familyName
 213                     print "  -Description:",familyDescription
 214                     print "  -Entity:",familyEntity
 215                     familyIdentifier = family.getIdentifier()
 216                     nbOfAtt = family.getNumberOfAttributes()
 217                     print "  -Identifier:",familyIdentifier
 218                     print "  -Number Of Attributes:",nbOfAtt
 219                     attributesids = family.getAttributesIdentifiers()
 220                     attributesvals = family.getAttributesValues()
 221                     for k in range(nbOfAtt):
 222                         print "    * Attributes:",attributesids[k],":",attributesvals[k],",",family.getAttributeDescription(k+1)
 223                     nbOfGrp = family.getNumberOfGroups()
 224                     print "  -Number Of Groups:",nbOfGrp
 225                     for k in range(nbOfGrp):
 226                         print "    * Group:",family.getGroupName(k+1)
 227                     print "  -Entities list:"
 228                     if (familyBool):
 229                         print "  -Is on all entities"
 230                     else:
 231                         nbOfTypes = family.getNumberOfTypes()
 232                         types = family.getTypes()
 233                         print "  -Number Of Types:",nbOfTypes
 234                         for k in range(nbOfTypes):
 235                             type = types[k]
 236                             nbOfElmtsOfType = family.getNumberOfElements(type)
 237                             number = family.getNumber(type)
 238                             print "    * Type",type
 239                             print "    * Number",number[0:nbOfElmtsOfType]
 240                         print ""
 241                         numberFamily = family.getNumber(MED_ALL_ELEMENTS)
 242                         print "    * Getting an Integer Field on the family ",familyName
 243                         fieldFamilyIntg = FIELDINT(family,spaceDim)
 244                         fieldFamilyIntg.setIterationNumber(0)
 245                         fieldFamilyIntg.setOrderNumber(0)
 246                         fieldFamilyIntg.setTime(0.0)
 247                         for kcomp in range(spaceDim):
 248                             kcomp1 = kcomp+1
 249                             if kcomp == 0:
 250                                 fieldCompName = "comp1"
 251                                 fieldCompDesc = "desc1"
 252                                 fieldCompUnit = "unit1"
 253                             if kcomp == 1:
 254                                 fieldCompName = "comp2"
 255                                 fieldCompDesc = "desc2"
 256                                 fieldCompUnit = "unit2"
 257                             if kcomp == 2:
 258                                 fieldCompName = "comp2"
 259                                 fieldCompDesc = "desc2"
 260                                 fieldCompUnit = "unit2"
 261
 262                             fieldFamilyIntg.setComponentName(kcomp1,fieldCompName)
 263                             fieldFamilyIntg.setComponentDescription(kcomp1,fieldCompDesc)
 264                             fieldFamilyIntg.setMEDComponentUnit(kcomp1,fieldCompUnit)
 265                         fieldFamilyName = "Integer Field on "+familyName
 266                         fieldFamilyIntg.setName(fieldFamilyName)
 267                         field_name = fieldFamilyIntg.getName()
 268                         type_field = fieldFamilyIntg.getValueType()
 269                         nbOfComp = fieldFamilyIntg.getNumberOfComponents()
 270                         print "      The field",field_name,"is with the type",type_field
 271                         print "      It has",nbOfComp,"Component(s)"
 272                         for kcomp in range(nbOfComp):
 273                             kcomp1 = kcomp+1
 274                             compName = fieldFamilyIntg.getComponentName(kcomp1)
 275                             compDesc = fieldFamilyIntg.getComponentDescription(kcomp1)
 276                             compUnit = fieldFamilyIntg.getMEDComponentUnit(kcomp1)
 277                             print "      * Component:",kcomp1
 278                             print "          Name:",compName
 279                             print "          Description:",compDesc
 280                             print "          Unit:",compUnit
 281
 282                         nbOf = fieldFamilyIntg.getSupport().getNumberOfElements(MED_ALL_ELEMENTS)
 283                         print "      Values:",nbOf
 284                         print "      Randomly set and get to check ..!"
 285                         for k in range(nbOf):
 286                             valueI = []
 287                             for kcomp in range(nbOfComp):
 288                                 valueI.append(randint(0,100))
 289
 290 #                            print "     Set Entry *",(k+1)," ",valueI[:nbOfComp]
 291                             valInd = numberFamily[k]
 292                             fieldFamilyIntg.setRow(valInd,valueI)
 293                             valueIverif = fieldFamilyIntg.getRow(valInd)
 294                             print "     Set/Get Entry *",(k+1)," ",valueI[:nbOfComp],"  /  ",valueIverif[:nbOfComp]
 295                         print "    * Getting a Real Field"
 296                         fieldFamilyDble = FIELDDOUBLE(family,spaceDim)
 297                         fieldFamilyDble.setIterationNumber(0)
 298                         fieldFamilyDble.setOrderNumber(0)
 299                         fieldFamilyDble.setTime(0.0)
 300                         for kcomp in range(spaceDim):
 301                             kcomp1 = kcomp+1
 302                             if kcomp == 0:
 303                                 fieldCompName = "comp1"
 304                                 fieldCompDesc = "desc1"
 305                                 fieldCompUnit = "unit1"
 306                             if kcomp == 1:
 307                                 fieldCompName = "comp2"
 308                                 fieldCompDesc = "desc2"
 309                                 fieldCompUnit = "unit2"
 310                             if kcomp == 2:
 311                                 fieldCompName = "comp2"
 312                                 fieldCompDesc = "desc2"
 313                                 fieldCompUnit = "unit2"
 314
 315                             fieldFamilyDble.setComponentName(kcomp1,fieldCompName)
 316                             fieldFamilyDble.setComponentDescription(kcomp1,fieldCompDesc)
 317                             fieldFamilyDble.setMEDComponentUnit(kcomp1,fieldCompUnit)
 318
 319                         fieldFamilyName = "Real Field on "+familyName
 320                         fieldFamilyDble.setName(fieldFamilyName)
 321                         field_name = fieldFamilyDble.getName()
 322                         type_field = fieldFamilyDble.getValueType()
 323                         nbOfComp = fieldFamilyDble.getNumberOfComponents()
 324                         print "      The field",field_name,"is with the type",type_field
 325                         print "      It has",nbOfComp,"Component(s)"
 326                         for kcomp in range(nbOfComp):
 327                             kcomp1 = kcomp+1
 328                             compName = fieldFamilyDble.getComponentName(kcomp1)
 329                             compDesc = fieldFamilyDble.getComponentDescription(kcomp1)
 330                             compUnit = fieldFamilyDble.getMEDComponentUnit(kcomp1)
 331                             print "      * Component:",kcomp1
 332                             print "          Name:",compName
 333                             print "          Description:",compDesc
 334                             print "          Unit:",compUnit
 335
 336                         nbOf = fieldFamilyDble.getSupport().getNumberOfElements(MED_ALL_ELEMENTS)
 337                         print "      Values:",nbOf
 338                         print "      Randomly set and get to check ..!"
 339                         for k in range(nbOf):
 340                             valueI = []
 341                             for kcomp in range(nbOfComp):
 342                                 valueI.append(random())
 343
 344 #                            print "     Set Entry *",(k+1)," ",valueI[:nbOfComp]
 345                             valInd = numberFamily[k]
 346                             fieldFamilyDble.setRow(valInd,valueI)
 347                             valueIverif = fieldFamilyDble.getRow(valInd)
 348                             print "     Set/Get Entry *",(k+1)," ",valueI[:nbOfComp],"  /  ",valueIverif[:nbOfComp]
 349                     if (entity != MED_NODE):
 350                         print ""
 351                         print "Getting barycenter on this family"
 352                         barycenterfamily = mesh.getBarycenter(family)
 353                         if (not familyBool): numberFamily = family.getNumber(MED_ALL_ELEMENTS)
 354                         nbVal = barycenterfamily.getSupport().getNumberOfElements(MED_ALL_ELEMENTS)
 355                         nbComp = barycenterfamily.getNumberOfComponents()
 356                         for j in range(nbVal):
 357                             valInd = j+1
 358                             if (not familyBool): valInd = numberFamily[j]
 359                             barycenterfamilyentity = barycenterfamily.getRow(valInd)
 360                             print "    * ",barycenterfamilyentity[:nbComp]
 361                 print ""
 362
 363         print "Building of the support on all Cells of the mesh."
 364         supportCell = mesh.getSupportOnAll( MED_CELL )
 365         print ""
 366         print "Getting barycenter of all Cells of the mesh"
 367         barycenter = mesh.getBarycenter(supportCell)
 368         for j in range(nbElemts):
 369             barycenterCell = barycenter.getRow(j+1)
 370             print "    * ",barycenterCell[:spaceDim]
 371
 372         print "Writing on file the mesh"
 373         writeMed21File = medFile[0:(len(medFile)-4)]+"_V21_fields.med"
 374         writeMed22File = medFile[0:(len(medFile)-4)]+"_V22_fields.med"
 375         fieldsMesh = barycenter.getSupport().getMesh()
 376         fieldsMeshName = "Fields Mesh"
 377         fieldsMesh.setName(fieldsMeshName)
 378
 379         index22Mesh = fieldsMesh.addDriver(MED_DRIVER,writeMed22File,fieldsMeshName)
 380         fieldsMesh.write(index22Mesh)
 381
 382         AnalyzeField(barycenter)
 383
 384         print "Writing on file the cells barycenter field"
 385
 386         barycenterName = barycenter.getName()
 387
 388         index22FieldBarycenter = barycenter.addDriver(MED_DRIVER,writeMed22File,barycenterName)
 389         barycenter.write(index22FieldBarycenter)
 390
 391         print ""
 392         if spaceDim == 3 :
 393             print "Getting volume of all Cells of the mesh:"
 394             volume = mesh.getVolume(supportCell)
 395             voltot = 0.
 396             for j in range(nbElemts):
 397                 volumeCell = volume.getValueIJ(j+1,1)
 398                 print "    * ",volumeCell
 399                 voltot = voltot + volumeCell
 400             print "Volume of the mesh:",voltot
 401             print ""
 402
 403             AnalyzeField(volume)
 404
 405             print "Writing on file the cells volume field"
 406
 407             volumeName = volume.getName()
 408
 409             index22FieldVolume = volume.addDriver(MED_DRIVER,writeMed22File,volumeName)
 410             volume.write(index22FieldVolume)
 411
 412             print ""
 413             print "Building of the support on all Faces of the mesh."
 414             supportFace = SUPPORT(mesh,"Support on all faces of the mesh",MED_FACE)
 415             supportFace.update()
 416             nbFace = mesh.getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS)
 417             print ""
 418             print "Getting normal of each face of this support",nbFace
 419             nbTypeFace = mesh.getNumberOfTypes(MED_FACE)
 420             TypeFace = mesh.getTypes(MED_FACE)
 421             print "nbTypeFace:",nbTypeFace,"----",TypeFace[:nbTypeFace]
 422             normal = mesh.getNormal(supportFace)
 423             for j in range(nbFace):
 424                 normalFace = normal.getRow(j+1)
 425                 value1 = normalFace[0]
 426                 value2 = normalFace[1]
 427                 value3 = normalFace[2]
 428                 norm = (value1*value1 + value2*value2 + value3*value3)**(0.5)
 429                 print "    * ",normalFace[:spaceDim],"norm:",norm
 430             print ""
 431
 432             AnalyzeField(normal)
 433
 434             print "Writing on file the face normal field"
 435
 436             normalName = normal.getName()
 437
 438             index22FieldNormal = normal.addDriver(MED_DRIVER,writeMed22File,normalName)
 439             normal.write(index22FieldNormal)
 440
 441         elif spaceDim == 2:
 442             print "Getting area on all Cells of the mesh:"
 443             area = mesh.getArea(supportCell)
 444             areatot = 0.
 445             for j in range(nbElemts):
 446                 areaCell = area.getValueIJ(j+1,1)
 447                 print "    * ",areaCell
 448                 areatot = areatot + areaCell
 449             print "Area of the mesh:",areatot
 450             print ""
 451
 452             AnalyzeField(area)
 453
 454             print "Writing on file the cells area field"
 455
 456             areaName = area.getName()
 457
 458             index22FieldArea = area.addDriver(MED_DRIVER,writeMed22File,areaName)
 459             area.write(index22FieldArea)
 460
 461             print ""
 462             print "Getting the support on all Edges of the mesh."
 463             supportEdge = mesh.getSupportOnAll(MED_EDGE)
 464             nbEdge = mesh.getNumberOfElements(MED_EDGE,MED_ALL_ELEMENTS)
 465             print ""
 466             print "Getting normal of each edge of this support",nbEdge
 467             nbTypeEdge = mesh.getNumberOfTypes(MED_EDGE)
 468             TypeEdge = mesh.getTypes(MED_EDGE)
 469             print "nbTypeEdge:",nbTypeEdge,"----",TypeEdge[:nbTypeEdge]
 470             normal = mesh.getNormal(supportEdge)
 471             for j in range(nbEdge):
 472                 normalEdge = normal.getRow(j+1)
 473                 value1 = normalEdge[0]
 474                 value2 = normalEdge[1]
 475                 norm = (value1*value1 + value2*value2)**(0.5)
 476                 print "    * ",normalEdge[:spaceDim],"norm:",norm
 477             print ""
 478
 479             AnalyzeField(normal)
 480
 481             print "Writing on file the edge normal field"
 482
 483             normalName = normal.getName()
 484
 485             index22FieldNormal = normal.addDriver(MED_DRIVER,writeMed22File,normalName)
 486             normal.write(index22FieldNormal)
 487         print ""
 488
 489 print "END of the Pyhton script ..... Ctrl D to exit"