src/Plugins/MEDWriter/plugin/Test/TestMEDWriter0.py

   1 # Copyright (C) 2016-2022  CEA/DEN, EDF R&D
   2 #
   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, or (at your option) any later version.
   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 # Author : Anthony Geay (EDF R&D)
  20
  21 #### import the simple module from the paraview
  22 from paraview.simple import *
  23 import MEDLoader as ml
  24 import os
  25 from math import pi,sqrt
  26 import tempfile
  27
  28 #### disable automatic camera reset on 'Show'
  29 paraview.simple._DisableFirstRenderCameraReset()
  30
  31 with tempfile.TemporaryDirectory(prefix="MEDWriter_") as tmpdir:
  32     pat=os.path.join(tmpdir, 'testMEDWriter_%i.med')
  33     fname0=pat%0
  34     fname1=pat%1
  35     fname2=pat%2
  36     fname3=pat%3
  37     fname4=pat%4
  38     fname4_vtp=os.path.splitext(pat%4)[0]+".vtp"
  39     fname5=pat%5
  40     fname6_vtu=os.path.splitext(pat%6)[0]+".vtu"
  41     fname6=pat%6
  42     fname7_vtu=os.path.splitext(pat%7)[0]+".vtu"
  43     fname7=pat%7
  44     fname8_vtr=os.path.splitext(pat%8)[0]+".vtr"
  45     fname8=pat%8
  46
  47     ##### First test with a simple sphere
  48
  49     plane1 = Sphere()
  50     SaveData(fname0,proxy=plane1,WriteAllTimeSteps=1)
  51     #
  52     totomed=MEDReader(FileName=fname0)
  53     totomed.AllArrays=['TS1/Mesh/ComSup0/Mesh@@][@@P0']
  54     totomed.AllTimeSteps=['0000']
  55     SaveData(fname1,proxy=totomed,WriteAllTimeSteps=1)
  56     # Sphere has been written. Try to check to write it in MED file !
  57     mfd=ml.MEDFileData(fname0)
  58     mm=mfd.getMeshes()[0] ; m0=mm[0]
  59     area=m0.getMeasureField(True).accumulate()[0]
  60     assert(abs(sqrt(area/(4*pi))-0.975/2.)<0.01) # 4*pi*radius**2
  61     f=mfd.getFields()[0][0].getFieldOnMeshAtLevel(ml.ON_NODES,0,mm)
  62     f_array_double = f.getArray().convertToDblArr() # f is MEDCouplingFieldFloat
  63     assert(abs(ml.DataArrayDouble(f_array_double.accumulate(),1,3).magnitude()[0])<1e-12) # sum of all normal vector should be 0
  64
  65
  66     ##### Build a MED file from scratch
  67
  68     fieldName0="F0"
  69     fieldName1="F1"
  70     c=ml.DataArrayDouble([0.,0.,1.,1.,2.,2.,1.,0.,2.,0.,0.,2.,1.,2.,0.,1.,2.,1.],9,2)
  71     c.setInfoOnComponents(["X abc","Y defg"])
  72     #
  73     mName="mesh"
  74     m0=ml.MEDCouplingUMesh(mName,2) ; m0.allocateCells() ; m0.setCoords(c)
  75     m0.insertNextCell(ml.NORM_TRI3,[3,1,4])
  76     m0.insertNextCell(ml.NORM_TRI3,[1,8,4])
  77     m0.insertNextCell(ml.NORM_QUAD4,[0,7,1,3])
  78     m0.insertNextCell(ml.NORM_QUAD4,[1,6,2,4])
  79     m0.insertNextCell(ml.NORM_QUAD4,[7,5,6,1])
  80     m1=ml.MEDCouplingUMesh(mName,1) ; m1.allocateCells() ; m1.setCoords(c)
  81     m1.insertNextCell(ml.NORM_SEG2,[0,7])
  82     m1.insertNextCell(ml.NORM_SEG2,[7,5])
  83     m1.insertNextCell(ml.NORM_SEG2,[5,6])
  84     mm=ml.MEDFileUMesh() ; mm[0]=m0 ; mm[-1]=m1
  85     mm.setFamilyFieldArr(1,ml.DataArrayInt([200,201,202,203,204,205,206,207,208]))
  86     mm.setFamilyFieldArr(0,ml.DataArrayInt([100,101,102,103,104]))
  87     mm.setFamilyFieldArr(-1,ml.DataArrayInt([105,106,107]))
  88     mm.setFamilyId("Fam3",3) ; mm.setFamilyId("Fam5",7)
  89     mm.setFamiliesOnGroup("gr0",["Fam3"])
  90     mm.setFamiliesOnGroup("gr1",["Fam5"])
  91     mm.setFamiliesOnGroup("gr2",["Fam3","Fam5"])
  92     mm.write(fname2,2)
  93     #
  94     f1ts0=ml.MEDFileField1TS()
  95     f0=ml.MEDCouplingFieldDouble(ml.ON_CELLS) ; f0.setName(fieldName0)
  96     f0.setMesh(m0) ; f0.setArray(ml.DataArrayDouble([8,7,6,5,4]))
  97     f1ts0.setFieldNoProfileSBT(f0)
  98     f0=ml.MEDCouplingFieldDouble(ml.ON_CELLS) ; f0.setName(fieldName0)
  99     f0.setMesh(m1) ; f0.setArray(ml.DataArrayDouble([3,2,1]))
 100     f0.setTime(0,0,0)
 101     f1ts0.setFieldNoProfileSBT(f0)
 102     f1ts0.write(fname2,0)
 103     #
 104     f1ts1=ml.MEDFileField1TS()
 105     f0=ml.MEDCouplingFieldDouble(ml.ON_NODES) ; f0.setName(fieldName1)
 106     arr=ml.DataArrayDouble([9,109,8,108,7,107,6,106,5,105,4,104,3,103,2,102,1,101],9,2)
 107     arr.setInfoOnComponents(["aa","bbb"])
 108     f0.setMesh(m0) ; f0.setArray(arr)
 109     f0.setTime(0,0,0)
 110     f1ts1.setFieldNoProfileSBT(f0)
 111     f1ts1.write(fname2,0)
 112     #
 113     test3=MEDReader(FileName=fname2)
 114     test3.AllArrays=['TS0/%s/ComSup0/%s@@][@@P0'%(mName,fieldName0),'TS0/%s/ComSup0/%s@@][@@P1'%(mName,fieldName1)]
 115     test3.AllTimeSteps = ['0000']
 116     SaveData(fname3,proxy=test3,WriteAllTimeSteps=1)
 117     ### test content of fname3
 118     mfd2=ml.MEDFileData(fname3)
 119     mm2=mfd2.getMeshes()[0]
 120     c1=mm2.getCoords()
 121     assert(c.isEqualWithoutConsideringStr(c1[:,:2],1e-12))
 122     fs2=ml.MEDFileFields(fname3)
 123     assert(len(fs2)==2)
 124     assert(mm2.getSpaceDimension()==3) ; assert(mm2.getCoords()[:,2].isUniform(0.,0.))
 125     m2_0=mm2[0].deepCopy() ; m2_0.changeSpaceDimension(2,0.) ; m2_0.getCoords().setInfoOnComponents(mm[0].getCoords().getInfoOnComponents())
 126     assert(m2_0.isEqual(mm[0],1e-12))
 127     m2_1=mm2[-1].deepCopy() ; m2_1.changeSpaceDimension(2,0.) ; m2_1.getCoords().setInfoOnComponents(mm[0].getCoords().getInfoOnComponents())
 128     assert(m2_1.isEqual(mm[-1],1e-12))
 129     f2_0=mfd2.getFields()[fieldName0][0].getFieldOnMeshAtLevel(ml.ON_CELLS,0,mm2) ; f2_0.setMesh(m2_0)
 130     assert(f1ts0.getFieldOnMeshAtLevel(ml.ON_CELLS,0,mm).isEqual(f2_0,1e-12,1e-12))
 131     f2_1=mfd2.getFields()[fieldName1][0].getFieldOnMeshAtLevel(ml.ON_NODES,0,mm2) ; f2_1.setMesh(m2_0)
 132     assert(f1ts1.getFieldOnMeshAtLevel(ml.ON_NODES,0,mm).isEqual(f2_1,1e-12,1e-12))
 133     assert(mm2.getGroupsNames()==('gr0','gr1','gr2'))
 134     assert(mm2.getFamiliesOnGroup("gr0")==("Fam3",))
 135     assert(mm2.getFamiliesOnGroup("gr1")==("Fam5",))
 136     assert(mm2.getFamiliesOnGroup("gr2")==("Fam3","Fam5"))
 137     assert(mm2.getFamiliesNames()==('FAMILLE_ZERO','Fam3','Fam5'))
 138     assert([mm2.getFamilyId(elt) for elt in ['FAMILLE_ZERO','Fam3','Fam5']]==[0,3,7])
 139     assert(mm2.getFamilyFieldAtLevel(0).isEqual(mm.getFamilyFieldAtLevel(0)))
 140     assert(mm2.getFamilyFieldAtLevel(-1).isEqual(mm.getFamilyFieldAtLevel(-1)))
 141     assert(mm2.getFamilyFieldAtLevel(1).isEqual(mm.getFamilyFieldAtLevel(1)))
 142     # Write a polydata mesh
 143     mergeBlocks1 = MergeBlocks(Input=test3)
 144     extractSurface1 = ExtractSurface(Input=mergeBlocks1)
 145     SaveData(fname4_vtp,proxy=extractSurface1)
 146     test4vtp = XMLPolyDataReader(FileName=[fname4_vtp])
 147     test4vtp.CellArrayStatus = ['F0', 'FamilyIdCell']
 148     SaveData(fname5,proxy=test4vtp,WriteAllTimeSteps=1)
 149     ### test content of fname5
 150     mfd5=ml.MEDFileData(fname5)
 151     m5=mfd5.getMeshes()[0][0].deepCopy()
 152     assert(m5.getSpaceDimension()==3) #
 153     m5.setName(mm.getName()) ; m5.changeSpaceDimension(2,0.) ; m5.getCoords().setInfoOnComponents(mm[0].getCoords().getInfoOnComponents())
 154     bary5=m5.computeCellCenterOfMass()
 155     bary=mm[0].computeCellCenterOfMass()
 156     a,b=bary5.areIncludedInMe(bary,1e-12) ; assert(a)
 157     a,c=mm[0].getCoords().areIncludedInMe(m5.getCoords(),1e-12) ; assert(a)
 158     m5.renumberNodes(c,len(c))#c.invertArrayO2N2N2O(len(c)))
 159     assert(m5.unPolyze())
 160     assert(m5.getCoords().isEqual(mm[0].getCoords(),1e-12))
 161     assert(m5.isEqual(mm[0],1e-12))
 162     f5_0=mfd5.getFields()[fieldName0][0].getFieldOnMeshAtLevel(ml.ON_CELLS,0,mfd5.getMeshes()[0]) ; f5_0.setMesh(m5)
 163     assert(f1ts0.getFieldOnMeshAtLevel(ml.ON_CELLS,0,mm).isEqual(f5_0,1e-12,1e-12))
 164     f5_1=mfd5.getFields()[fieldName1][0].getFieldOnMeshAtLevel(ml.ON_NODES,0,mfd5.getMeshes()[0]) ; f5_1.setMesh(m5)
 165     f5_1.setArray(f5_1.getArray()[c.invertArrayO2N2N2O(len(c))])
 166     assert(f1ts1.getFieldOnMeshAtLevel(ml.ON_NODES,0,mm).isEqual(f5_1,1e-12,1e-12))
 167
 168     ### test with a non geo types non sorted
 169
 170     c=ml.DataArrayDouble([0.,0.,1.,1.,2.,2.,1.,0.,2.,0.,0.,2.,1.,2.,0.,1.,2.,1.],9,2)
 171     c.setInfoOnComponents(["X abc","Y defg"])
 172     m6=ml.MEDCouplingUMesh(mName,2) ; m6.allocateCells() ; m6.setCoords(c)
 173     m6.insertNextCell(ml.NORM_TRI3,[3,1,4])
 174     m6.insertNextCell(ml.NORM_QUAD4,[0,7,1,3])
 175     m6.insertNextCell(ml.NORM_TRI3,[1,8,4])
 176     m6.insertNextCell(ml.NORM_QUAD4,[1,6,2,4])
 177     m6.insertNextCell(ml.NORM_QUAD4,[7,5,6,1])
 178     fieldName6="F6"
 179     f6=ml.MEDCouplingFieldDouble(ml.ON_CELLS) ; f6.setMesh(m6) ; f6.setName(fieldName6)
 180     f6.setArray(ml.DataArrayDouble([20,21,22,23,24]))
 181     f6.writeVTK(fname6_vtu)
 182     test6vtu=XMLUnstructuredGridReader(FileName=[fname6_vtu])
 183     SaveData(fname6,proxy=test6vtu,WriteAllTimeSteps=1)
 184     mfd7=ml.MEDFileData(fname6)
 185     assert(len(mfd7.getMeshes())==1)
 186     m7=mfd7.getMeshes()[0][0]
 187     assert(len(mfd7.getFields())==1)
 188     f7=mfd7.getFields()[0][0].getFieldOnMeshAtLevel(ml.ON_CELLS,0,mfd7.getMeshes()[0])
 189     assert(f7.getMesh().isEqual(m7,1e-12))
 190     assert(m7.getCoords()[:,:2].isEqualWithoutConsideringStr(m6.getCoords(),1e-12))
 191     assert(m7.getNodalConnectivity().isEqual(ml.DataArrayInt([3,3,1,4,3,1,8,4,4,0,7,1,3,4,1,6,2,4,4,7,5,6,1]))) # there is a permutation of cells
 192     assert(m7.getNodalConnectivityIndex().isEqual(ml.DataArrayInt([0,4,8,13,18,23]))) # there is a permutation of cells
 193     assert(f7.getArray().isEqual(ml.DataArrayFloat([20,22,21,23,24]),1e-12)) # there is a permutation of cells
 194
 195     ### test with polyhedron
 196
 197     m8=ml.MEDCouplingCMesh() ; m8.setCoords(ml.DataArrayDouble([0,1,2,3]),ml.DataArrayDouble([0,1]),ml.DataArrayDouble([0,1]))
 198     m8=m8.buildUnstructured()
 199     m8.getCoords().setInfoOnComponents(['X', 'Y', 'Z'])
 200     m8_0=m8[0] ; m8_0.simplexize(ml.PLANAR_FACE_5)
 201     m8_1=m8[1:]
 202     m8_1.convertAllToPoly()
 203     m8=ml.MEDCouplingUMesh.MergeUMeshesOnSameCoords([m8_0,m8_1])
 204     m8=m8[[0,5,1,2,6,3,4]]
 205     fieldName8="F8"
 206     f8=ml.MEDCouplingFieldDouble(ml.ON_CELLS) ; f8.setMesh(m8) ; f8.setName(fieldName8)
 207     f8.setArray(ml.DataArrayDouble([20,21,22,23,24,25,26]))
 208     f8.writeVTK(fname7_vtu)
 209     test8vtu=XMLUnstructuredGridReader(FileName=[fname7_vtu])
 210     SaveData(fname7,proxy=test8vtu,WriteAllTimeSteps=1)
 211     mfd9=ml.MEDFileData(fname7)
 212     assert(len(mfd9.getMeshes())==1)
 213     m9=mfd9.getMeshes()[0][0]
 214     assert(len(mfd9.getFields())==1)
 215     assert(m9.getCoords().isEqual(m8.getCoords(),1e-12))
 216     c9=ml.DataArrayInt([0,2,3,5,6,1,4])
 217     assert(m8[c9].isEqualWithoutConsideringStr(m9,1e-12))
 218     f9=mfd9.getFields()[0][0].getFieldOnMeshAtLevel(ml.ON_CELLS,0,mfd9.getMeshes()[0])
 219     f9_array_double = f9.getArray().convertToDblArr() # f9 is MEDCouplingFieldFloat
 220     assert(f9_array_double.isEqual(f8.getArray()[c9],1e-12))
 221
 222     ### test with cartesian
 223
 224     FieldName10="F10"
 225     FieldName10_n="F10_n"
 226     m10=ml.MEDCouplingCMesh()
 227     m10.setCoordsAt(0,ml.DataArrayDouble([0,1,2]))
 228     m10.setCoordsAt(1,ml.DataArrayDouble([1,2,3,4]))
 229     m10.setCoordsAt(2,ml.DataArrayDouble([3,5,6,7,8]))
 230     f10=ml.MEDCouplingFieldDouble(ml.ON_CELLS) ; f10.setMesh(m10)
 231     f10.setName(FieldName10)
 232     f10.setArray(ml.DataArrayInt.Range(0,m10.getNumberOfCells(),1).convertToDblArr()) ; f10.checkConsistencyLight()
 233     f10_n=ml.MEDCouplingFieldDouble(ml.ON_NODES) ; f10_n.setMesh(m10)
 234     f10_n.setName(FieldName10_n)
 235     f10_n.setArray(ml.DataArrayInt.Range(0,m10.getNumberOfNodes(),1).convertToDblArr()) ; f10_n.checkConsistencyLight()
 236     ml.MEDCouplingFieldDouble.WriteVTK(fname8_vtr,[f10,f10_n])
 237     test10vtr=XMLRectilinearGridReader(FileName=[fname8_vtr])
 238     SaveData(fname8,proxy=test10vtr,WriteAllTimeSteps=1)
 239     mfd11=ml.MEDFileData(fname8)
 240     assert(len(mfd11.getMeshes())==1)
 241     assert(len(mfd11.getFields())==2)
 242     mfd11=ml.MEDFileData(fname8)
 243     m11=mfd11.getMeshes()[0]
 244     assert(isinstance(m11,ml.MEDFileCMesh))
 245     f11=mfd11.getFields()[FieldName10][0].getFieldOnMeshAtLevel(ml.ON_CELLS,0,m11)
 246     f11_n=mfd11.getFields()[FieldName10_n][0].getFieldOnMeshAtLevel(ml.ON_NODES,0,m11)
 247     f11_array_double = f11.getArray().convertToDblArr() # f11 is MEDCouplingFieldFloat
 248     assert(f11_array_double.isEqualWithoutConsideringStr(f10.getArray(),1e-12))
 249     f11_n_array_double = f11_n.getArray().convertToDblArr() # f11_n is MEDCouplingFieldFloat
 250     assert(f11_n_array_double.isEqualWithoutConsideringStr(f10_n.getArray(),1e-12))
 251