Back to index

salome-med  6.5.0
med_test2.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 # Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00005 # CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00006 #
00007 # This library is free software; you can redistribute it and/or
00008 # modify it under the terms of the GNU Lesser General Public
00009 # License as published by the Free Software Foundation; either
00010 # version 2.1 of the License.
00011 #
00012 # This library is distributed in the hope that it will be useful,
00013 # but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00015 # Lesser General Public License for more details.
00016 #
00017 # You should have received a copy of the GNU Lesser General Public
00018 # License along with this library; if not, write to the Free Software
00019 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00020 #
00021 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00022 #
00023 
00024 ###################################################################################
00025 # This Python script is parsing a MED file using MED Memory from SALOME platform:
00026 # It analyses all meshes in the MED file (coordinates, connectivity of d-cells as
00027 # well as (d-1)-cells, families), it tests all fields generated in the MESH class
00028 # and write them in a 2 new med files (V2.1 and V2.2), it gives only the number of
00029 # fields stored in the MED file (d is the mesh/space dimension).
00030 ###################################################################################
00031 #
00032 from libMEDMEM_Swig import *
00033 from random import *
00034 
00035 #==============================================================================
00036 
00037 def AnalyzeField(field):
00038     name = field.getName()
00039     desc = field.getDescription()
00040     nbComp = field.getNumberOfComponents()
00041     itNum = field.getIterationNumber()
00042     ordNum = field.getOrderNumber()
00043     type = field.getValueType()
00044 
00045     print "Analysis of the field ",name," with the description ",desc
00046     print "iteration number ",itNum," order Number ",ordNum
00047     print "It has ",nbComp," component(s) with the type ",type
00048 
00049     fieldValue = field.getValue()
00050     fieldSupport = field.getSupport()
00051     fieldMesh = fieldSupport.getMesh()
00052     fieldEntity = fieldSupport.getEntity()
00053     bool = fieldSupport.isOnAllElements()
00054 
00055     if bool:
00056         print "The support of this field is on all entities ",fieldEntity," of the mesh ",fieldMesh.getName()
00057         if fieldEntity == MED_NODE:
00058             nbValByComp = fieldMesh.getNumberOfNodes()
00059         else:
00060             nbValByComp = fieldMesh.getNumberOfElements(fieldEntity,MED_ALL_ELEMENTS)
00061         print "and its dimension (number of values by component of the field) is ",nbValByComp
00062     else:
00063         print "The support of this field is partially on entities ",fieldEntity," of the mesh ",fieldMesh.getName()
00064         nbValByComp = fieldSupport.getNumberOfElements(MED_ALL_ELEMENTS)
00065         print "and its dimension (number of values by component of the field) is ",nbValByComp
00066 
00067     for i in range(nbComp):
00068         ip1 = i + 1
00069         compName = field.getComponentName(ip1)
00070         compDesc = field.getComponentDescription(ip1)
00071         compUnit = field.getMEDComponentUnit(ip1)
00072         print "The ",(i+1),"-th  component ",compName," with the dexription ",compDesc," and the unit ",compUnit
00073 
00074     for i in range(nbValByComp):
00075         print "  * ",fieldValue[i*nbComp:(i+1)*nbComp]
00076 
00077 #==============================================================================
00078 import os
00079 #
00080 #befor running this script, please be sure about the path the file fileName
00081 #
00082 filePath=os.environ["MED_ROOT_DIR"]
00083 filePath=os.path.join(filePath, "share", "salome", "resources", "med")
00084 
00085 medFile = os.path.join(filePath, "carre_en_quad4_seg2.med")
00086 #medFile = os.path.join(filePath, "cube_hexa8_quad4.med")
00087 
00088 def print_ord(i):
00089     if i == 0:
00090         return 'first'
00091     elif i == 1:
00092         return 'second'
00093     elif i == 2:
00094         return 'third'
00095     else:
00096         return `i`+'th'
00097 
00098 md = MEDFILEBROWSER(medFile)
00099 
00100 nbMeshes = md.getNumberOfMeshes()
00101 
00102 nbFields = md.getNumberOfFields()
00103 
00104 print "The med file", medFile, "contains", nbMeshes, "mesh(es) and", nbFields, "field(s)"
00105 
00106 if (nbMeshes>0):
00107     print "Mesh(es) Name(s) is(are) "
00108 
00109     for i in range(nbMeshes):
00110         mesh_name = md.getMeshName(i)
00111         print "   - ",mesh_name
00112 
00113 if (nbFields>0):
00114     print "Field(s) Name(s) is(are) "
00115 
00116     for i in range(nbFields):
00117         field_name = md.getFieldName(i)
00118         print "   - ",field_name
00119 
00120 print ""
00121 
00122 if (nbMeshes>0):
00123     print "Mesh(es) Analysis "
00124     for i in range(nbMeshes):
00125         mesh_name = md.getMeshName(i)
00126         mesh = MESH(MED_DRIVER,md.getFileName(),mesh_name)
00127         spaceDim = mesh.getSpaceDimension()
00128         meshDim = mesh.getMeshDimension()
00129         print "The",print_ord(i), "mesh, '",mesh_name,"', is a",spaceDim,"D mesh on a",meshDim,"D geometry"
00130         nbNodes = mesh.getNumberOfNodes()
00131         print "The mesh has",nbNodes,"Nodes"
00132         coordSyst = mesh.getCoordinatesSystem()
00133         print "The coordinates system is",coordSyst
00134         print "The Coordinates :"
00135         coordNames = []
00136         coordUnits = []
00137         for isd in range(spaceDim):
00138             coordNames.append(mesh.getCoordinateName(isd))
00139             coordUnits.append(mesh.getCoordinateUnit(isd))
00140 
00141         print "names:", coordNames
00142         print "units", coordUnits
00143         print "values:"
00144         coordinates = mesh.getCoordinates(MED_FULL_INTERLACE)
00145         for k in range(nbNodes):
00146             kp1 = k+1
00147             coords = []
00148             for isd in range(spaceDim):
00149                 isdp1 = isd+1
00150                 coords.append(mesh.getCoordinate(kp1,isdp1))
00151 
00152             print coords," ---- ", coordinates[k*spaceDim:((k+1)*spaceDim)]
00153 
00154         print ""
00155         print "Show the Nodal Connectivity:"
00156         nbTypesCell = mesh.getNumberOfTypes(MED_CELL)
00157         print ""
00158         if (nbTypesCell>0):
00159             print "The Mesh has",nbTypesCell,"Type(s) of Cell"
00160             types = mesh.getTypes(MED_CELL)
00161             for k in range(nbTypesCell):
00162                 type = types[k]
00163                 nbElemType = mesh.getNumberOfElements(MED_CELL,type)
00164                 print "For the type:",type,"there is(are)",nbElemType,"elemnt(s)"
00165                 connectivity = mesh.getConnectivity(MED_NODAL,MED_CELL,type)
00166                 nbNodesPerCell = type%100
00167                 for j in range(nbElemType):
00168                     print "Element",(j+1)," ",connectivity[j*nbNodesPerCell:(j+1)*nbNodesPerCell]
00169 
00170         print ""
00171         print "Show the Reverse Nodal Connectivity:"
00172         ReverseConnectivity = mesh.getReverseConnectivity(MED_NODAL)
00173         ReverseConnectivityIndex = mesh.getReverseConnectivityIndex(MED_NODAL)
00174         print ""
00175         for j in range(nbNodes):
00176             begin = ReverseConnectivityIndex[j]-1
00177             end = ReverseConnectivityIndex[j+1]-1
00178             print "Node",(j+1),"-->",ReverseConnectivity[begin:end]
00179 
00180         print ""
00181         print "Show the Descending Connectivity:"
00182         mesh.calculateConnectivity(MED_DESCENDING,MED_CELL)
00183         nbElemts = mesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)
00184         Connectivity = mesh.getConnectivity(MED_DESCENDING,MED_CELL,MED_ALL_ELEMENTS)
00185         ConnectivityIndex = mesh.getConnectivityIndex(MED_DESCENDING,MED_CELL)
00186         print ""
00187         for j in range(nbElemts):
00188             begin = ConnectivityIndex[j]-1
00189             end = ConnectivityIndex[j+1]-1
00190             print "Element",(j+1),"-->",Connectivity[begin:end]
00191 
00192         print ""
00193         for entity in [MED_NODE,MED_CELL,MED_FACE,MED_EDGE]:
00194             nbFam = mesh.getNumberOfFamilies(entity)
00195             if (entity == MED_NODE) & (nbFam > 0):
00196                 print "This mesh has",nbFam,"Node Family(ies)"
00197             elif (entity == MED_CELL) & (nbFam > 0):
00198                 print "This mesh has",nbFam,"Cell Family(ies)"
00199             elif (entity == MED_FACE) & (nbFam > 0):
00200                 print "This mesh has",nbFam,"Face Family(ies)"
00201             elif (entity == MED_EDGE) & (nbFam > 0):
00202                 print "This mesh has",nbFam,"Edge Family(ies)"
00203 
00204             if nbFam > 0:
00205                 for j in range(nbFam):
00206                     print ""
00207                     family = mesh.getFamily(entity,j+1)
00208                     familyName = family.getName()
00209                     familyDescription = family.getDescription()
00210                     familyEntity = family.getEntity()
00211                     familyBool = family.isOnAllElements()
00212                     print "  -Name:",familyName
00213                     print "  -Description:",familyDescription
00214                     print "  -Entity:",familyEntity
00215                     familyIdentifier = family.getIdentifier()
00216                     nbOfAtt = family.getNumberOfAttributes()
00217                     print "  -Identifier:",familyIdentifier
00218                     print "  -Number Of Attributes:",nbOfAtt
00219                     attributesids = family.getAttributesIdentifiers()
00220                     attributesvals = family.getAttributesValues()
00221                     for k in range(nbOfAtt):
00222                         print "    * Attributes:",attributesids[k],":",attributesvals[k],",",family.getAttributeDescription(k+1)
00223                     nbOfGrp = family.getNumberOfGroups()
00224                     print "  -Number Of Groups:",nbOfGrp
00225                     for k in range(nbOfGrp):
00226                         print "    * Group:",family.getGroupName(k+1)
00227                     print "  -Entities list:"
00228                     if (familyBool):
00229                         print "  -Is on all entities"
00230                     else:
00231                         nbOfTypes = family.getNumberOfTypes()
00232                         types = family.getTypes()
00233                         print "  -Number Of Types:",nbOfTypes
00234                         for k in range(nbOfTypes):
00235                             type = types[k]
00236                             nbOfElmtsOfType = family.getNumberOfElements(type)
00237                             number = family.getNumber(type)
00238                             print "    * Type",type
00239                             print "    * Number",number[0:nbOfElmtsOfType]
00240                         print ""
00241                         numberFamily = family.getNumber(MED_ALL_ELEMENTS)
00242                         print "    * Getting an Integer Field on the family ",familyName
00243                         fieldFamilyIntg = FIELDINT(family,spaceDim)
00244                         fieldFamilyIntg.setIterationNumber(0)
00245                         fieldFamilyIntg.setOrderNumber(0)
00246                         fieldFamilyIntg.setTime(0.0)
00247                         for kcomp in range(spaceDim):
00248                             kcomp1 = kcomp+1
00249                             if kcomp == 0:
00250                                 fieldCompName = "comp1"
00251                                 fieldCompDesc = "desc1"
00252                                 fieldCompUnit = "unit1"
00253                             if kcomp == 1:
00254                                 fieldCompName = "comp2"
00255                                 fieldCompDesc = "desc2"
00256                                 fieldCompUnit = "unit2"
00257                             if kcomp == 2:
00258                                 fieldCompName = "comp2"
00259                                 fieldCompDesc = "desc2"
00260                                 fieldCompUnit = "unit2"
00261 
00262                             fieldFamilyIntg.setComponentName(kcomp1,fieldCompName)
00263                             fieldFamilyIntg.setComponentDescription(kcomp1,fieldCompDesc)
00264                             fieldFamilyIntg.setMEDComponentUnit(kcomp1,fieldCompUnit)
00265                         fieldFamilyName = "Integer Field on "+familyName
00266                         fieldFamilyIntg.setName(fieldFamilyName)
00267                         field_name = fieldFamilyIntg.getName()
00268                         type_field = fieldFamilyIntg.getValueType()
00269                         nbOfComp = fieldFamilyIntg.getNumberOfComponents()
00270                         print "      The field",field_name,"is with the type",type_field
00271                         print "      It has",nbOfComp,"Component(s)"
00272                         for kcomp in range(nbOfComp):
00273                             kcomp1 = kcomp+1
00274                             compName = fieldFamilyIntg.getComponentName(kcomp1)
00275                             compDesc = fieldFamilyIntg.getComponentDescription(kcomp1)
00276                             compUnit = fieldFamilyIntg.getMEDComponentUnit(kcomp1)
00277                             print "      * Component:",kcomp1
00278                             print "          Name:",compName
00279                             print "          Description:",compDesc
00280                             print "          Unit:",compUnit
00281 
00282                         nbOf = fieldFamilyIntg.getSupport().getNumberOfElements(MED_ALL_ELEMENTS)
00283                         print "      Values:",nbOf
00284                         print "      Randomly set and get to check ..!"
00285                         for k in range(nbOf):
00286                             valueI = []
00287                             for kcomp in range(nbOfComp):
00288                                 valueI.append(randint(0,100))
00289 
00290 #                            print "     Set Entry *",(k+1)," ",valueI[:nbOfComp]
00291                             valInd = numberFamily[k]                           
00292                             fieldFamilyIntg.setRow(valInd,valueI)
00293                             valueIverif = fieldFamilyIntg.getRow(valInd)
00294                             print "     Set/Get Entry *",(k+1)," ",valueI[:nbOfComp],"  /  ",valueIverif[:nbOfComp]
00295                         print "    * Getting a Real Field"
00296                         fieldFamilyDble = FIELDDOUBLE(family,spaceDim)
00297                         fieldFamilyDble.setIterationNumber(0)
00298                         fieldFamilyDble.setOrderNumber(0)
00299                         fieldFamilyDble.setTime(0.0)
00300                         for kcomp in range(spaceDim):
00301                             kcomp1 = kcomp+1
00302                             if kcomp == 0:
00303                                 fieldCompName = "comp1"
00304                                 fieldCompDesc = "desc1"
00305                                 fieldCompUnit = "unit1"
00306                             if kcomp == 1:
00307                                 fieldCompName = "comp2"
00308                                 fieldCompDesc = "desc2"
00309                                 fieldCompUnit = "unit2"
00310                             if kcomp == 2:
00311                                 fieldCompName = "comp2"
00312                                 fieldCompDesc = "desc2"
00313                                 fieldCompUnit = "unit2"
00314 
00315                             fieldFamilyDble.setComponentName(kcomp1,fieldCompName)
00316                             fieldFamilyDble.setComponentDescription(kcomp1,fieldCompDesc)
00317                             fieldFamilyDble.setMEDComponentUnit(kcomp1,fieldCompUnit)
00318 
00319                         fieldFamilyName = "Real Field on "+familyName
00320                         fieldFamilyDble.setName(fieldFamilyName)
00321                         field_name = fieldFamilyDble.getName()
00322                         type_field = fieldFamilyDble.getValueType()
00323                         nbOfComp = fieldFamilyDble.getNumberOfComponents()
00324                         print "      The field",field_name,"is with the type",type_field
00325                         print "      It has",nbOfComp,"Component(s)"
00326                         for kcomp in range(nbOfComp):
00327                             kcomp1 = kcomp+1
00328                             compName = fieldFamilyDble.getComponentName(kcomp1)
00329                             compDesc = fieldFamilyDble.getComponentDescription(kcomp1)
00330                             compUnit = fieldFamilyDble.getMEDComponentUnit(kcomp1)
00331                             print "      * Component:",kcomp1
00332                             print "          Name:",compName
00333                             print "          Description:",compDesc
00334                             print "          Unit:",compUnit
00335 
00336                         nbOf = fieldFamilyDble.getSupport().getNumberOfElements(MED_ALL_ELEMENTS)
00337                         print "      Values:",nbOf
00338                         print "      Randomly set and get to check ..!"
00339                         for k in range(nbOf):
00340                             valueI = []
00341                             for kcomp in range(nbOfComp):
00342                                 valueI.append(random())
00343 
00344 #                            print "     Set Entry *",(k+1)," ",valueI[:nbOfComp]
00345                             valInd = numberFamily[k]
00346                             fieldFamilyDble.setRow(valInd,valueI)
00347                             valueIverif = fieldFamilyDble.getRow(valInd)
00348                             print "     Set/Get Entry *",(k+1)," ",valueI[:nbOfComp],"  /  ",valueIverif[:nbOfComp]
00349                     if (entity != MED_NODE):
00350                         print ""
00351                         print "Getting barycenter on this family"
00352                         barycenterfamily = mesh.getBarycenter(family)
00353                         if (not familyBool): numberFamily = family.getNumber(MED_ALL_ELEMENTS)
00354                         nbVal = barycenterfamily.getSupport().getNumberOfElements(MED_ALL_ELEMENTS)
00355                         nbComp = barycenterfamily.getNumberOfComponents()
00356                         for j in range(nbVal):
00357                             valInd = j+1
00358                             if (not familyBool): valInd = numberFamily[j]
00359                             barycenterfamilyentity = barycenterfamily.getRow(valInd)
00360                             print "    * ",barycenterfamilyentity[:nbComp]
00361                 print ""
00362 
00363         print "Building of the support on all Cells of the mesh."
00364         supportCell = mesh.getSupportOnAll( MED_CELL )
00365         print ""
00366         print "Getting barycenter of all Cells of the mesh"
00367         barycenter = mesh.getBarycenter(supportCell)
00368         for j in range(nbElemts):
00369             barycenterCell = barycenter.getRow(j+1)
00370             print "    * ",barycenterCell[:spaceDim]
00371 
00372         print "Writing on file the mesh"
00373         writeMed21File = medFile[0:(len(medFile)-4)]+"_V21_fields.med"
00374         writeMed22File = medFile[0:(len(medFile)-4)]+"_V22_fields.med"
00375         fieldsMesh = barycenter.getSupport().getMesh()
00376         fieldsMeshName = "Fields Mesh"
00377         fieldsMesh.setName(fieldsMeshName)
00378 
00379         index22Mesh = fieldsMesh.addDriver(MED_DRIVER,writeMed22File,fieldsMeshName)
00380         fieldsMesh.write(index22Mesh)
00381 
00382         AnalyzeField(barycenter)
00383 
00384         print "Writing on file the cells barycenter field"
00385 
00386         barycenterName = barycenter.getName()
00387 
00388         index22FieldBarycenter = barycenter.addDriver(MED_DRIVER,writeMed22File,barycenterName)
00389         barycenter.write(index22FieldBarycenter)
00390 
00391         print ""
00392         if spaceDim == 3 :
00393             print "Getting volume of all Cells of the mesh:"
00394             volume = mesh.getVolume(supportCell)
00395             voltot = 0.
00396             for j in range(nbElemts):
00397                 volumeCell = volume.getValueIJ(j+1,1)
00398                 print "    * ",volumeCell
00399                 voltot = voltot + volumeCell
00400             print "Volume of the mesh:",voltot
00401             print ""
00402 
00403             AnalyzeField(volume)
00404 
00405             print "Writing on file the cells volume field"
00406 
00407             volumeName = volume.getName()
00408 
00409             index22FieldVolume = volume.addDriver(MED_DRIVER,writeMed22File,volumeName)
00410             volume.write(index22FieldVolume)
00411 
00412             print ""
00413             print "Building of the support on all Faces of the mesh."
00414             supportFace = SUPPORT(mesh,"Support on all faces of the mesh",MED_FACE)
00415             supportFace.update()
00416             nbFace = mesh.getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS)
00417             print ""
00418             print "Getting normal of each face of this support",nbFace
00419             nbTypeFace = mesh.getNumberOfTypes(MED_FACE)
00420             TypeFace = mesh.getTypes(MED_FACE)
00421             print "nbTypeFace:",nbTypeFace,"----",TypeFace[:nbTypeFace]
00422             normal = mesh.getNormal(supportFace)
00423             for j in range(nbFace):
00424                 normalFace = normal.getRow(j+1)
00425                 value1 = normalFace[0]
00426                 value2 = normalFace[1]
00427                 value3 = normalFace[2]
00428                 norm = (value1*value1 + value2*value2 + value3*value3)**(0.5)
00429                 print "    * ",normalFace[:spaceDim],"norm:",norm
00430             print ""
00431 
00432             AnalyzeField(normal)
00433 
00434             print "Writing on file the face normal field"
00435 
00436             normalName = normal.getName()
00437 
00438             index22FieldNormal = normal.addDriver(MED_DRIVER,writeMed22File,normalName)
00439             normal.write(index22FieldNormal)
00440 
00441         elif spaceDim == 2:
00442             print "Getting area on all Cells of the mesh:"
00443             area = mesh.getArea(supportCell)
00444             areatot = 0.
00445             for j in range(nbElemts):
00446                 areaCell = area.getValueIJ(j+1,1)
00447                 print "    * ",areaCell
00448                 areatot = areatot + areaCell
00449             print "Area of the mesh:",areatot
00450             print ""
00451 
00452             AnalyzeField(area)
00453 
00454             print "Writing on file the cells area field"
00455 
00456             areaName = area.getName()
00457 
00458             index22FieldArea = area.addDriver(MED_DRIVER,writeMed22File,areaName)
00459             area.write(index22FieldArea)
00460 
00461             print ""
00462             print "Getting the support on all Edges of the mesh."
00463             supportEdge = mesh.getSupportOnAll(MED_EDGE)
00464             nbEdge = mesh.getNumberOfElements(MED_EDGE,MED_ALL_ELEMENTS)
00465             print ""
00466             print "Getting normal of each edge of this support",nbEdge
00467             nbTypeEdge = mesh.getNumberOfTypes(MED_EDGE)
00468             TypeEdge = mesh.getTypes(MED_EDGE)
00469             print "nbTypeEdge:",nbTypeEdge,"----",TypeEdge[:nbTypeEdge]
00470             normal = mesh.getNormal(supportEdge)
00471             for j in range(nbEdge):
00472                 normalEdge = normal.getRow(j+1)
00473                 value1 = normalEdge[0]
00474                 value2 = normalEdge[1]
00475                 norm = (value1*value1 + value2*value2)**(0.5)
00476                 print "    * ",normalEdge[:spaceDim],"norm:",norm
00477             print ""
00478 
00479             AnalyzeField(normal)
00480 
00481             print "Writing on file the edge normal field"
00482 
00483             normalName = normal.getName()
00484 
00485             index22FieldNormal = normal.addDriver(MED_DRIVER,writeMed22File,normalName)
00486             normal.write(index22FieldNormal)
00487         print ""
00488 
00489 print "END of the Pyhton script ..... Ctrl D to exit"