src/Plugins/MEDWriter/Test/TestMEDWriter0.py

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