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

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