src/MEDLoader/Swig/MEDLoaderCouplingTrainingSession.py

   1 #  -*- coding: iso-8859-1 -*-
   2 # Copyright (C) 2007-2020  CEA/DEN, EDF R&D
   3 #
   4 # This library is free software; you can redistribute it and/or
   5 # modify it under the terms of the GNU Lesser General Public
   6 # License as published by the Free Software Foundation; either
   7 # version 2.1 of the License, or (at your option) any later version.
   8 #
   9 # This library is distributed in the hope that it will be useful,
  10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 # Lesser General Public License for more details.
  13 #
  14 # You should have received a copy of the GNU Lesser General Public
  15 # License along with this library; if not, write to the Free Software
  16 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 #
  18 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 #
  20 # Author : Anthony GEAY (CEA/DEN/DM2S/STMF/LGLS)
  21
  22 from MEDLoader import *
  23 from MEDCouplingRemapper import *
  24 import math, os
  25 import tempfile,os,shutil
  26
  27 zeDir = tempfile.mkdtemp()
  28 os.chdir(zeDir)
  29
  30 d=DataArrayDouble.New(6,2)
  31 d[:,0]=3.
  32 d[:, 1] = list(range(6))
  33 d[:,1]*=math.pi/3.
  34 d=d.fromPolarToCart()
  35 d.setInfoOnComponents(["X [m]","Y [m]"])
  36 print(d.getValues())
  37 print(d)
  38 print(d.magnitude().isUniform(3.,1e-12))
  39 #
  40 radius=3.
  41 translationToPerform=[[0.,0.],[3./2.*radius,-radius*math.sqrt(3.)/2],[3./2.*radius,radius*math.sqrt(3.)/2],[0.,radius*math.sqrt(3.)],[-3./2.*radius,radius*math.sqrt(3.)/2],[-3./2.*radius,-radius*math.sqrt(3.)/2],[0.,-radius*math.sqrt(3.)]]
  42 ds=len(translationToPerform)*[None]
  43 for pos,t in enumerate(translationToPerform):
  44   ds[pos]=d[:]
  45   ds[pos]+=t
  46   pass
  47 #
  48 d2=DataArrayDouble.Aggregate(ds)
  49 oldNbOfTuples=d2.getNumberOfTuples()
  50 c,cI=d2.findCommonTuples(1e-12)
  51 tmp=c[cI[0]:cI[0+1]]
  52 print(tmp)
  53 a=cI.deltaShiftIndex()
  54 b=a-1
  55 myNewNbOfTuples=oldNbOfTuples-sum(b.getValues())
  56 o2n,newNbOfTuples=DataArrayInt.ConvertIndexArrayToO2N(oldNbOfTuples,c,cI)
  57 print("Ai je trouve le bon resultat ? %s"%(str(myNewNbOfTuples==newNbOfTuples))) ; assert myNewNbOfTuples==newNbOfTuples
  58 #
  59 d3=d2.renumberAndReduce(o2n,newNbOfTuples)
  60 n2o=o2n.invertArrayO2N2N2O(newNbOfTuples)
  61 d3_bis=d2[n2o]
  62 print("Ai je trouve le bon resultat (2) ? %s"%(str(d3.isEqual(d3_bis,1e-12)))) ; assert d3.isEqual(d3_bis,1e-12)
  63 #
  64 d3+=[3.3,4.4]
  65 # d3 contains coordinates
  66 m=MEDCouplingUMesh.New("My7hexagons",2)
  67 m.setCoords(d3)
  68 m.allocateCells(7)
  69 for i in range(7):
  70   m.insertNextCell(NORM_POLYGON,o2n[6*i:6*(i+1)].getValues())
  71   pass
  72 m.finishInsertingCells()
  73 m.checkConsistencyLight()
  74 #
  75 m.writeVTK("My7hexagons.vtu")
  76
  77 ########
  78
  79 coords=[0.,0.,0., 1.,1.,0., 1.,1.25,0., 1.,0.,0., 1.,1.5,0., 2.,0.,0., 2.,1.,0., 1.,2.,0., 0.,2.,0., 3.,1.,0.,
  80         3.,2.,0., 0.,1.,0., 1.,3.,0., 2.,2.,0., 2.,3.,0.,
  81         0.,0.,1., 1.,1.,1., 1.,1.25,1., 1.,0.,1., 1.,1.5,1., 2.,0.,1., 2.,1.,1., 1.,2.,1., 0.,2.,1., 3.,1.,1.,
  82         3.,2.,1., 0.,1.,1., 1.,3.,1., 2.,2.,1., 2.,3.,1.,
  83         0.,0.,2., 1.,1.,2., 1.,1.25,2., 1.,0.,2., 1.,1.5,2., 2.,0.,2., 2.,1.,2., 1.,2.,2., 0.,2.,2., 3.,1.,2.,
  84         3.,2.,2., 0.,1.,2., 1.,3.,2., 2.,2.,2., 2.,3.,2.,
  85         0.,0.,3., 1.,1.,3., 1.,1.25,3., 1.,0.,3., 1.,1.5,3., 2.,0.,3., 2.,1.,3., 1.,2.,3., 0.,2.,3., 3.,1.,3.,
  86         3.,2.,3., 0.,1.,3., 1.,3.,3., 2.,2.,3., 2.,3.,3.]
  87 conn=[0,11,1,3,15,26,16,18,   1,2,4,7,13,6,-1,1,16,21,6,-1,6,21,28,13,-1,13,7,22,28,-1,7,4,19,22,-1,4,2,17,19,-1,2,1,16,17,-1,16,21,28,22,19,17,
  88       1,6,5,3,16,21,20,18,   13,10,9,6,28,25,24,21, 11,8,7,4,2,1,-1,11,26,16,1,-1,1,16,17,2,-1,2,17,19,4,-1,4,19,22,7,-1,7,8,23,22,-1,8,11,26,23,-1,26,16,17,19,22,23,
  89       7,12,14,13,22,27,29,28,  15,26,16,18,30,41,31,33, 16,17,19,22,28,21,-1,16,31,36,21,-1,21,36,43,28,-1,28,22,37,43,-1,22,19,34,37,-1,19,17,32,34,-1,17,16,31,32,-1,31,36,43,37,34,32,
  90       16,21,20,18,31,36,35,33,   28,25,24,21,43,40,39,36, 26,23,22,19,17,16,-1,26,41,31,16,-1,16,31,32,17,-1,17,32,34,19,-1,19,34,37,22,-1,22,23,38,37,-1,23,26,41,38,-1,41,31,32,34,37,38,
  91       22,27,29,28,37,42,44,43, 30,41,31,33,45,56,46,48,  31,32,34,37,43,36,-1,31,46,51,36,-1,36,51,58,43,-1,43,37,52,58,-1,37,34,49,52,-1,34,32,47,49,-1,32,31,46,47,-1,46,51,58,52,49,47,
  92       31,36,35,33,46,51,50,48,  43,40,39,36,58,55,54,51, 41,38,37,34,32,31,-1,41,56,46,31,-1,31,46,47,32,-1,32,47,49,34,-1,34,49,52,37,-1,37,38,53,52,-1,38,41,56,53,-1,56,46,47,49,52,53,
  93       37,42,44,43,52,57,59,58]
  94 mesh3D=MEDCouplingUMesh.New("mesh3D",3);
  95 mesh3D.allocateCells(18);
  96 mesh3D.insertNextCell(NORM_HEXA8,conn[0:8]); mesh3D.insertNextCell(NORM_POLYHED,conn[8:51]); mesh3D.insertNextCell(NORM_HEXA8,conn[51:59]); mesh3D.insertNextCell(NORM_HEXA8,conn[59:67]); mesh3D.insertNextCell(NORM_POLYHED,conn[67:110]); mesh3D.insertNextCell(NORM_HEXA8,conn[110:118]);
  97 mesh3D.insertNextCell(NORM_HEXA8,conn[118:126]); mesh3D.insertNextCell(NORM_POLYHED,conn[126:169]); mesh3D.insertNextCell(NORM_HEXA8,conn[169:177]); mesh3D.insertNextCell(NORM_HEXA8,conn[177:185]); mesh3D.insertNextCell(NORM_POLYHED,conn[185:228]); mesh3D.insertNextCell(NORM_HEXA8,conn[228:236]);
  98 mesh3D.insertNextCell(NORM_HEXA8,conn[236:244]); mesh3D.insertNextCell(NORM_POLYHED,conn[244:287]); mesh3D.insertNextCell(NORM_HEXA8,conn[287:295]); mesh3D.insertNextCell(NORM_HEXA8,conn[295:303]); mesh3D.insertNextCell(NORM_POLYHED,conn[303:346]); mesh3D.insertNextCell(NORM_HEXA8,conn[346:354]);
  99 mesh3D.finishInsertingCells();
 100 myCoords=DataArrayDouble.New(coords,60,3);
 101 myCoords.setInfoOnComponents(["X [m]","Y [m]","Z [m]"])
 102 mesh3D.setCoords(myCoords);
 103 mesh3D.orientCorrectlyPolyhedrons()
 104 mesh3D.sortCellsInMEDFileFrmt()
 105 mesh3D.checkConsistencyLight()
 106 renum = DataArrayInt.New(60) ; renum[:15] = list(range(15, 30)) ; renum[15:30] = list(range(15)) ; renum[30:45] = list(range(45, 60)) ; renum[45:] = list(range(30, 45))
 107 mesh3D.renumberNodes(renum,60)
 108 #
 109 mesh3D.getCoords()[:]*=100.
 110 mesh3D.getCoords().setInfoOnComponents(["X [cm]","Y [cm]","Z [cm]"])
 111 #
 112 zLev=mesh3D.getCoords()[:,2]
 113 zLev = zLev.getDifferentValues(1e-12)
 114 zLev.sort()
 115 #
 116 tmp,cellIdsSol1=mesh3D.buildSlice3D([0.,0.,(zLev[1]+zLev[2])/2],[0.,0.,1.],1e-12)
 117 bary=mesh3D.computeCellCenterOfMass()
 118 baryZ=bary[:,2]
 119 cellIdsSol2=baryZ.findIdsInRange(zLev[1],zLev[2])
 120 nodeIds=mesh3D.findNodesOnPlane([0.,0.,zLev[0]],[0.,0.,1.],1e-10)
 121 mesh2D=mesh3D.buildFacePartOfMySelfNode(nodeIds,True)
 122 extMesh=MEDCouplingMappedExtrudedMesh.New(mesh3D,mesh2D,0)
 123 cellIdsSol3=extMesh.getMesh3DIds()[mesh2D.getNumberOfCells():2*mesh2D.getNumberOfCells()]
 124 for i in range(3):
 125   exec("print( cellIdsSol%s.getValues())"%(i+1))
 126 #
 127 mesh3DPart=mesh3D[cellIdsSol2] # equivalent to mesh3DPart=mesh3D.buildPartOfMySelf(cellIdsSol2,True)
 128 mesh3DPart.zipCoords()
 129 print(mesh3DPart.checkConsecutiveCellTypesAndOrder([NORM_HEXA8,NORM_POLYHED])) ; assert mesh3DPart.checkConsecutiveCellTypesAndOrder([NORM_HEXA8,NORM_POLYHED])
 130 print(mesh3DPart.checkConsecutiveCellTypes()) ; assert mesh3DPart.checkConsecutiveCellTypes()
 131 #print mesh3DPart.advancedRepr()
 132 #
 133 baryXY=bary[:,[0,1]]
 134 baryXY-=[250.,150.]
 135 magn=baryXY.magnitude()
 136 cellIds2Sol1=magn.findIdsInRange(0.,1e-12)
 137 #
 138 bary2=mesh2D.computeCellCenterOfMass()[:,[0,1]]
 139 bary2-=[250.,150.]
 140 magn=bary2.magnitude()
 141 ids=magn.findIdsInRange(0.,1e-12)
 142 idStart=int(ids) # ids is assumed to contain only one value, if not an exception is thrown
 143 cellIds2Sol2 = extMesh.getMesh3DIds()[list(range(idStart, mesh3D.getNumberOfCells(), mesh2D.getNumberOfCells()))]
 144 #
 145 mesh3DSlice2=mesh3D[cellIds2Sol1]
 146 mesh3DSlice2.zipCoords()
 147 #
 148 mesh3DSlice2bis=mesh3DSlice2.deepCopy()
 149 mesh3DSlice2bis.translate([0.,1000.,0.])
 150 mesh3DSlice2All=MEDCouplingUMesh.MergeUMeshes([mesh3DSlice2,mesh3DSlice2bis])
 151 mesh3DSlice2All.writeVTK("mesh3DSlice2All.vtu")
 152 #
 153 mesh3DSurf,desc,descIndx,revDesc,revDescIndx=mesh3D.buildDescendingConnectivity()
 154 numberOf3DCellSharing=revDescIndx.deltaShiftIndex()
 155 cellIds=numberOf3DCellSharing.findIdsNotEqual(1)
 156 mesh3DSurfInside=mesh3DSurf[cellIds]
 157 mesh3DSurfInside.writeVTK("mesh3DSurfInside.vtu")
 158
 159 ######
 160
 161 xarr=DataArrayDouble.New(11,1)
 162 xarr.iota(0.)
 163 cmesh=MEDCouplingCMesh.New()
 164 cmesh.setCoords(xarr,xarr,xarr)
 165 mesh=cmesh.buildUnstructured()
 166 mesh.convertToPolyTypes(DataArrayInt.Range(0,mesh.getNumberOfCells(),2))
 167 #
 168
 169 f=mesh.fillFromAnalytic(ON_CELLS,1,"(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)")
 170 f.setName("MyField")
 171 #
 172 f2=MEDCouplingFieldDouble.New(ON_CELLS,ONE_TIME)
 173 f2.setMesh(mesh)
 174 f2.setName("MyField2")
 175 f2.fillFromAnalytic(1,"(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)")
 176 print("f and f2 are equal : %s"%(f.isEqualWithoutConsideringStr(f2,1e-13,1e-12))) ; assert f.isEqualWithoutConsideringStr(f2,1e-13,1e-12)
 177 #
 178 ids1=f.getArray().findIdsInRange(0.,5.)
 179 fPart1=f.buildSubPart(ids1)
 180 ids2=f.getArray().findIdsInRange(50.,1.e300)
 181 fPart2=f.buildSubPart(ids2)
 182 #Renumbering cells to follow MED file
 183 fPart1Cpy=fPart1.deepCopy()
 184 o2n=fPart1Cpy.getMesh().sortCellsInMEDFileFrmt()
 185 fPart1Cpy.getArray().renumberInPlace(o2n)
 186 #Check that fPart1Cpy and fPart1 are the same
 187 fPart1Cpy.substractInPlaceDM(fPart1,12,1e-12)
 188 fPart1Cpy.getArray().abs()
 189 print("Fields are the same ? %s"%(fPart1Cpy.getArray().accumulate()[0]<1e-12)) ; assert fPart1Cpy.getArray().accumulate()[0]<1e-12
 190 #
 191 fPart12=MEDCouplingFieldDouble.MergeFields([fPart1,fPart2])
 192 # evaluation on points
 193 bary=fPart12.getMesh().computeCellCenterOfMass()
 194 arr1=fPart12.getValueOnMulti(bary)
 195 arr2=f.getValueOnMulti(bary)
 196 delta=arr1-arr2
 197 delta.abs()
 198 print("Check OK : %s"%(delta.accumulate()[0]<1e-12)) ; assert delta.accumulate()[0]<1e-12
 199 #
 200 print(abs(fPart12.integral(0,True)-fPart12.getArray().accumulate()[0])<1e-10) ; assert abs(fPart12.integral(0,True)-fPart12.getArray().accumulate()[0])<1e-10
 201 fPart12.getMesh().scale([0.,0.,0.],1.2)
 202 print(abs(fPart12.integral(0,True)-fPart12.getArray().accumulate()[0]*1.2*1.2*1.2)<1e-8) ; assert abs(fPart12.integral(0,True)-fPart12.getArray().accumulate()[0]*1.2*1.2*1.2)<1e-8
 203 # Explosion of field
 204 fVec=mesh.fillFromAnalytic(ON_CELLS,3,"(x-5.)*IVec+(y-5.)*JVec+(z-5.)*KVec")
 205 fVecPart1=fVec.buildSubPart(ids1)
 206 fVecPart1.setName("fVecPart1")
 207 cells=fPart1.getMesh().getNumberOfCells()*[None]
 208 for icell,vec in enumerate(fVecPart1.getArray()):
 209   m=fPart1.getMesh()[[icell]]
 210   m.zipCoords()
 211   m.translate(vec)
 212   cells[icell]=m
 213   pass
 214 meshFVecPart1Exploded=MEDCouplingUMesh.MergeUMeshes(cells)
 215 fPart1.setMesh(meshFVecPart1Exploded)
 216
 217 ####
 218
 219 targetCoords=[-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 ];
 220 targetConn=[0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4];
 221 targetMesh=MEDCouplingUMesh.New("MyMesh",2);
 222 targetMesh.allocateCells(5);
 223 targetMesh.insertNextCell(NORM_TRI3,3,targetConn[4:7]);
 224 targetMesh.insertNextCell(NORM_TRI3,3,targetConn[7:10]);
 225 targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[0:4]);
 226 targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[10:14]);
 227 targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[14:18]);
 228 targetMesh.finishInsertingCells();
 229 myCoords=DataArrayDouble.New(targetCoords,9,2);
 230 myCoords.setInfoOnComponents(["X [km]","YY [mm]"])
 231 targetMesh.setCoords(myCoords);
 232 #
 233 WriteUMesh("TargetMesh.med",targetMesh,True)
 234 #
 235 meshRead=ReadUMeshFromFile("TargetMesh.med",targetMesh.getName(),0)
 236 print("Is the mesh read in file equals targetMesh ? %s"%(meshRead.isEqual(targetMesh,1e-12))) ; assert meshRead.isEqual(targetMesh,1e-12)
 237 #
 238 f=MEDCouplingFieldDouble.New(ON_CELLS,ONE_TIME)
 239 f.setTime(5.6,7,8)
 240 f.setArray(targetMesh.computeCellCenterOfMass())
 241 f.setMesh(targetMesh)
 242 f.setName("AFieldName")
 243 WriteField("MyFirstField.med",f,True)
 244 #
 245 f2=ReadFieldCell("MyFirstField.med",f.getMesh().getName(),0,f.getName(),7,8)
 246 print("Is the field read in file equals f ? %s"%(f2.isEqual(f,1e-12,1e-12))) ; assert f2.isEqual(f,1e-12,1e-12)
 247 #
 248 WriteUMesh("MySecondField.med",f.getMesh(),True)
 249 WriteFieldUsingAlreadyWrittenMesh("MySecondField.med",f)
 250 #
 251 f2=f.clone(True)
 252 f2.getArray()[:]*=2.0
 253 f2.setTime(7.8,9,10)
 254 WriteFieldUsingAlreadyWrittenMesh("MySecondField.med",f2)
 255 #
 256 f3=ReadFieldCell("MySecondField.med",f.getMesh().getName(),0,f.getName(),7,8)
 257 print("Is the field read in file equals f ? %s"%(f.isEqual(f3,1e-12,1e-12))) ; assert f.isEqual(f3,1e-12,1e-12)
 258 f4=ReadFieldCell("MySecondField.med",f.getMesh().getName(),0,f.getName(),9,10)
 259 print("Is the field read in file equals f ? %s"%(f2.isEqual(f4,1e-12,1e-12))) ; assert f2.isEqual(f4,1e-12,1e-12)
 260
 261 #####
 262
 263 targetCoords=[-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 ];
 264 targetConn=[0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4];
 265 targetMesh=MEDCouplingUMesh.New("MyMesh",2);
 266 targetMesh.allocateCells(5);
 267 targetMesh.insertNextCell(NORM_TRI3,3,targetConn[4:7]);
 268 targetMesh.insertNextCell(NORM_TRI3,3,targetConn[7:10]);
 269 targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[0:4]);
 270 targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[10:14]);
 271 targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[14:18]);
 272 targetMesh.finishInsertingCells();
 273 myCoords=DataArrayDouble.New(targetCoords,9,2);
 274 myCoords.setInfoOnComponents(["X [km]","YY [mm]"])
 275 targetMesh.setCoords(myCoords);
 276 #
 277 targetMeshConsti=targetMesh.buildDescendingConnectivity()[0]
 278 targetMesh1=targetMeshConsti[[3,4,7,8]]
 279 targetMesh1.setName(targetMesh.getName())
 280 #
 281 meshMEDFile=MEDFileUMesh.New()
 282 meshMEDFile.setMeshAtLevel(0,targetMesh)
 283 meshMEDFile.setMeshAtLevel(-1,targetMesh1)
 284 # Some groups on cells Level 0
 285 grp0_0=DataArrayInt.New([0,1,3]) ; grp0_0.setName("grp0_Lev0")
 286 grp1_0=DataArrayInt.New([1,2,3,4]) ; grp1_0.setName("grp1_Lev0")
 287 meshMEDFile.setGroupsAtLevel(0,[grp0_0,grp1_0])
 288 # Some groups on cells Level -1
 289 grp0_M1=DataArrayInt.New([0,1]) ; grp0_M1.setName("grp0_LevM1")
 290 grp1_M1=DataArrayInt.New([0,1,2]) ; grp1_M1.setName("grp1_LevM1")
 291 grp2_M1=DataArrayInt.New([1,2,3]) ; grp2_M1.setName("grp2_LevM1")
 292 meshMEDFile.setGroupsAtLevel(-1,[grp0_M1,grp1_M1,grp2_M1])
 293 #
 294 meshMEDFile.write("TargetMesh2.med",2) # 2 stands for write from scratch
 295 #
 296 meshMEDFileRead=MEDFileMesh.New("TargetMesh2.med")
 297 meshRead0=meshMEDFileRead.getMeshAtLevel(0)
 298 meshRead1=meshMEDFileRead.getMeshAtLevel(-1)
 299 print("Is the mesh at level 0 read in file equals targetMesh ? %s"%(meshRead0.isEqual(targetMesh,1e-12))) ; assert meshRead0.isEqual(targetMesh,1e-12)
 300 print("Is the mesh at level -1 read in file equals targetMesh ? %s"%(meshRead1.isEqual(targetMesh1,1e-12))) ; assert meshRead1.isEqual(targetMesh1,1e-12)
 301 #
 302 print(meshMEDFileRead.getGrpNonEmptyLevels("grp0_Lev0"))
 303 grp0_0_read=meshMEDFileRead.getGroupArr(0,"grp0_Lev0")
 304 print("Is group \"grp0_Lev0\" are the same ? %s"%(grp0_0_read.isEqual(grp0_0))) ; assert grp0_0_read.isEqual(grp0_0)
 305 #
 306 # Fields
 307 #
 308 f=MEDCouplingFieldDouble.New(ON_CELLS,ONE_TIME)
 309 f.setTime(5.6,7,8)
 310 f.setArray(targetMesh.computeCellCenterOfMass())
 311 f.setMesh(targetMesh)
 312 f.setName("AFieldName")
 313 #
 314 fMEDFile=MEDFileField1TS.New()
 315 fMEDFile.setFieldNoProfileSBT(f)
 316 #
 317 fMEDFile.write("TargetMesh2.med",0) # 0 is very important here because we want to append to TargetMesh2.med and not to scratch it
 318 #
 319 fMEDFileRead=MEDFileField1TS.New("TargetMesh2.med",f.getName(),7,8)
 320 fRead1=fMEDFileRead.getFieldOnMeshAtLevel(ON_CELLS,0,meshMEDFileRead) # fastest method. No read in file.
 321 fRead2=fMEDFileRead.getFieldAtLevel(ON_CELLS,0) # basic method like, mesh is reread in file...
 322 print("Does the field f remains the same using fast method ? %s"%(fRead1.isEqual(f,1e-12,1e-12))) ; assert fRead1.isEqual(f,1e-12,1e-12)
 323 print("Does the field f remains the same using slow method ? %s"%(fRead2.isEqual(f,1e-12,1e-12))) ; assert fRead2.isEqual(f,1e-12,1e-12)
 324 #
 325 # Writing and Reading fields on profile using MEDLoader advanced API
 326 #
 327 pfl=DataArrayInt.New([1,2,3]) ; pfl.setName("My1stPfl")
 328 fPart=f.buildSubPart(pfl)
 329 fPart.setName("fPart")
 330 #
 331 fMEDFile2=MEDFileField1TS.New()
 332 fMEDFile2.setFieldProfile(fPart,meshMEDFileRead,0,pfl)
 333 fMEDFile2.write("TargetMesh2.med",0) # 0 is very important here because we want to append to TargetMesh2.med and not to scratch it
 334 #
 335 fMEDFileRead2=MEDFileField1TS.New("TargetMesh2.med",fPart.getName(),7,8)
 336 fPartRead,pflRead=fMEDFileRead2.getFieldWithProfile(ON_CELLS,0,meshMEDFileRead)
 337 print(fPartRead.isEqualWithoutConsideringStr(fPart.getArray(),1e-12)) ; assert fPartRead.isEqualWithoutConsideringStr(fPart.getArray(),1e-12)
 338 print(pflRead.isEqualWithoutConsideringStr(pfl)) ; assert pflRead.isEqualWithoutConsideringStr(pfl)
 339
 340 #####
 341
 342 m0=MEDCouplingCMesh()
 343 arr=DataArrayDouble(31,1) ; arr.iota(0.)
 344 m0.setCoords(arr,arr)
 345 m0=m0.buildUnstructured()
 346 m00=m0[::2] ; m00.simplexize(0) ; m01=m0[1::2]
 347 m0=MEDCouplingUMesh.MergeUMeshes([m00,m01])
 348 m0.getCoords()[:]*=1/15.
 349 m0.setName("mesh")
 350 #
 351 CellField=MEDCouplingFieldDouble(ON_CELLS,ONE_TIME) ; CellField.setTime(5.6,5,6) ; CellField.setMesh(m0)
 352 CellField.setName("CellField")
 353 CellField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))") ; CellField.getArray().setInfoOnComponent(0,"powercell [W]")
 354 NodeField=MEDCouplingFieldDouble(ON_NODES,ONE_TIME) ; NodeField.setTime(5.6,5,6) ; NodeField.setMesh(m0)
 355 NodeField.setName("NodeField")
 356 NodeField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))") ; NodeField.getArray().setInfoOnComponent(0,"powernode [W]")
 357 #
 358 proc0=m0.getCellsInBoundingBox([(0.,0.4),(0.,0.4)],1e-10)
 359 proc1=proc0.buildComplement(m0.getNumberOfCells())
 360 #
 361 NodeField0=NodeField[proc0] ; NodeField0.getMesh().setName(m0.getName()) ; CellField0=CellField[proc0] ; CellField0.setMesh(NodeField0.getMesh())
 362 NodeField1=NodeField[proc1] ; NodeField1.getMesh().setName(m0.getName()) ; CellField1=CellField[proc1] ; CellField1.setMesh(NodeField1.getMesh())
 363 #
 364 proc0_fname="proc0.med"
 365 WriteField(proc0_fname,NodeField0,True)
 366 WriteFieldUsingAlreadyWrittenMesh(proc0_fname,CellField0)
 367 proc1_fname="proc1.med"
 368 WriteField(proc1_fname,NodeField1,True)
 369 WriteFieldUsingAlreadyWrittenMesh(proc1_fname,CellField1)
 370 #
 371 CellField0_read=ReadFieldCell("proc0.med","mesh",0,"CellField",5,6)
 372 CellField1_read=ReadFieldCell("proc1.med","mesh",0,"CellField",5,6)
 373 CellField_read=MEDCouplingFieldDouble.MergeFields([CellField0_read,CellField1_read])
 374 CellFieldCpy=CellField.deepCopy()
 375 CellFieldCpy.substractInPlaceDM(CellField_read,10,1e-12)
 376 CellFieldCpy.getArray().abs()
 377 print(CellFieldCpy.getArray().isUniform(0.,1e-12))
 378 #
 379 NodeField0_read=ReadFieldNode("proc0.med","mesh",0,"NodeField",5,6)
 380 NodeField1_read=ReadFieldNode("proc1.med","mesh",0,"NodeField",5,6)
 381 NodeField_read=MEDCouplingFieldDouble.MergeFields([NodeField0_read,NodeField1_read])
 382 NodeField_read.mergeNodes(1e-10)
 383 NodeFieldCpy=NodeField.deepCopy()
 384 NodeFieldCpy.mergeNodes(1e-10)
 385 NodeFieldCpy.substractInPlaceDM(NodeField_read,10,1e-12)
 386 print(NodeFieldCpy.getArray().isUniform(0.,1e-12)) ; assert NodeFieldCpy.getArray().isUniform(0.,1e-12)
 387 #
 388 fileNames=["proc0.med","proc1.med"]
 389 msML=[MEDFileMesh.New(fname) for fname in fileNames]
 390 fsML=[MEDFileFields.New(fname) for fname in fileNames]
 391 mergeMLMesh=MEDFileUMesh()
 392 mergeMLFields=MEDFileFields()
 393 for lev in msML[0].getNonEmptyLevels():
 394     o2nML=len(msML[0].getNonEmptyLevels())*[None]
 395     cs=[mML.getCoords() for mML in msML]
 396     mergeMLMesh.setCoords(DataArrayDouble.Aggregate(cs))
 397     ms=[mML.getMeshAtLevel(lev) for mML in msML]
 398     m=MEDCouplingUMesh.MergeUMeshes(ms) ; m.setCoords(mergeMLMesh.getCoords())
 399     o2nML[lev]=m.sortCellsInMEDFileFrmt()
 400     mergeMLMesh.setMeshAtLevel(lev,m)
 401     pass
 402 #
 403 for fieldName in fsML[0].getFieldsNames():
 404     fmts=[fML[fieldName] for fML in fsML]
 405     mergeField=MEDFileFieldMultiTS()
 406     for dt,it,tim in fmts[0].getTimeSteps():
 407         fts=[fmt[dt,it] for fmt in fmts]
 408         arrs=len(fts)*[None]
 409         for typp in fts[0].getTypesOfFieldAvailable():
 410             arr1s=[]
 411             if typp==ON_CELLS:
 412                for ft in fts:
 413                    for geoTyp,smth in ft.getFieldSplitedByType():
 414                        if geoTyp!=NORM_ERROR:
 415                            smth1=[x for x in smth if x[0]==ON_CELLS]
 416                            arr2s=[ft.getUndergroundDataArray()[elt[1][0]:elt[1][1]] for elt in smth1]
 417                            arr1s.append(DataArrayDouble.Aggregate(arr2s))
 418                            pass
 419                        pass
 420                    pass
 421                pass
 422             else:
 423                 for ft in fts:
 424                     smth=[x for x in ft.getFieldSplitedByType() if x[0]==NORM_ERROR]
 425                     arr2=DataArrayDouble.Aggregate([ft.getUndergroundDataArray()[elt[1][0][1][0]:elt[1][0][1][1]] for elt in smth])
 426                     arr1s.append(arr2)
 427                     pass
 428                 pass
 429             arr=DataArrayDouble.Aggregate(arr1s)
 430             if typp==ON_CELLS:
 431                arr.renumberInPlace(o2nML[lev])
 432             mcf=MEDCouplingFieldDouble(typp,ONE_TIME) ; mcf.setName(fieldName) ; mcf.setTime(tim,dt,it) ; mcf.setArray(arr)
 433             mcf.setMesh(mergeMLMesh.getMeshAtLevel(lev)) ; mcf.checkConsistencyLight()
 434             mergeField.appendFieldNoProfileSBT(mcf)
 435             pass
 436         pass
 437     mergeMLFields.pushField(mergeField)
 438     pass
 439 mergeMLMesh.write("merge.med",2)
 440 mergeMLFields.write("merge.med",0)
 441
 442 #####
 443
 444 arr=DataArrayDouble(11) ; arr.iota(0)
 445 trgMesh=MEDCouplingCMesh() ; trgMesh.setCoords(arr,arr) ; trgMesh=trgMesh.buildUnstructured()
 446 #
 447 arr=DataArrayDouble(21) ; arr.iota(0) ; arr*=0.5
 448 srcMesh=MEDCouplingCMesh() ; srcMesh.setCoords(arr,arr) ; srcMesh=srcMesh.buildUnstructured()
 449 #
 450 tmp=srcMesh[:20] ; tmp.simplexize(0)
 451 srcMesh=MEDCouplingUMesh.MergeUMeshes([tmp,srcMesh[20:]])
 452 #
 453 remap=MEDCouplingRemapper()
 454 remap.prepare(srcMesh,trgMesh,"P0P0")
 455 #
 456 myMatrix=remap.getCrudeMatrix()
 457 print(myMatrix) # pour voir a quoi elle ressemble
 458 sumByRows=DataArrayDouble(len(myMatrix))
 459 for i,wIt in enumerate(sumByRows):
 460   su=0.
 461   for it in myMatrix[i]:
 462     su+=myMatrix[i][it]
 463   wIt[0]=su
 464 print("Does interpolation look OK ? %s"%(str(sumByRows.isUniform(1.,1e-12)))) ; assert sumByRows.isUniform(1.,1e-12)
 465 #
 466 srcField=MEDCouplingFieldDouble(ON_CELLS,ONE_TIME) ; srcField.setMesh(srcMesh)
 467 srcField.fillFromAnalytic(1,"7-sqrt((x-5.)*(x-5.)+(y-5.)*(y-5.))") ; CellField.getArray().setInfoOnComponent(0,"powercell [W]")
 468 #
 469 #remap.transferField(srcField,1e300)
 470 srcField.setNature(IntensiveMaximum)
 471 trgFieldCV=remap.transferField(srcField,1e300)
 472 #
 473 print("IntensiveMaximum %lf == %lf"%(srcField.integral(True)[0],trgFieldCV.integral(True)[0])) ; assert abs(srcField.integral(True)[0]-trgFieldCV.integral(True)[0])<1e-6
 474 print("IntensiveMaximum %lf != %lf"%(srcField.getArray().accumulate()[0],trgFieldCV.getArray().accumulate()[0])) ; assert abs(srcField.getArray().accumulate()[0]-trgFieldCV.getArray().accumulate()[0])>1e-6
 475 #
 476 srcField.setNature(ExtensiveMaximum)
 477 trgFieldI=remap.transferField(srcField,1e300)
 478 #
 479 print("ExtensiveConservation %lf != %lf"%(srcField.integral(True)[0],trgFieldI.integral(True)[0])) ; assert abs(srcField.integral(True)[0]-trgFieldI.integral(True)[0])>1e-6
 480 print("ExtensiveConservation %lf == %lf"%(srcField.getArray().accumulate()[0],trgFieldI.getArray().accumulate()[0])) ; assert abs(srcField.getArray().accumulate()[0]-trgFieldI.getArray().accumulate()[0])<1e-6
 481
 482 ######
 483
 484 from numpy import *
 485 from math import acos
 486
 487 med_root_dir=os.getenv("MEDCOUPLING_ROOT_DIR")
 488 agitateur_file = ""
 489 if med_root_dir:
 490   agitateur_file = os.path.join(os.getenv("MEDCOUPLING_ROOT_DIR"),"share","resources","med","agitateur.med")
 491 if not os.path.exists(agitateur_file):
 492   current_dir = os.path.dirname(os.path.realpath(__file__))
 493   agitateur_file=os.path.join(current_dir, "..", "..", "..", "resources","agitateur.med")
 494 pass
 495 data=MEDFileData(agitateur_file)
 496 ts=data.getFields()[0].getTimeSteps()
 497 print(ts)
 498 #
 499 fMts=data.getFields()["DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"]
 500 f1ts=fMts[(2,-1)]
 501 fMc=f1ts.getFieldAtLevel(ON_CELLS,0)
 502 arr=fMc.getArray()
 503 arr.getMinMaxPerComponent() # juste pour voir la plage de variation du champ par compo
 504 ids=arr.findIdsInRange(0.,1.)
 505 f2Mc=fMc[ids]
 506 #
 507 pressMts=data.getFields()["PRESSION_ELEM_DOM"]
 508 press1ts=pressMts[(2,-1)]
 509 pressMc=press1ts.getFieldAtLevel(ON_CELLS,0)
 510 pressOnAgitateurMc=pressMc[ids]
 511 #
 512 pressOnAgitateurMc.getMesh().zipCoords()
 513 #
 514 agitateurMesh3DMc=pressOnAgitateurMc.getMesh()
 515 m3DSurf,desc,descI,revDesc,revDescI=agitateurMesh3DMc.buildDescendingConnectivity()
 516 nbOf3DCellSharing=revDescI.deltaShiftIndex()
 517 ids2=nbOf3DCellSharing.findIdsEqual(1)
 518 agitateurSkinMc=m3DSurf[ids2]
 519 OffsetsOfTupleIdsInField=revDescI[ids2]
 520 tupleIdsInField=revDesc[OffsetsOfTupleIdsInField]
 521 pressOnSkinAgitateurMc=pressOnAgitateurMc[tupleIdsInField]
 522 pressOnSkinAgitateurMc.setMesh(agitateurSkinMc)
 523 #
 524 pressSkin=pressOnSkinAgitateurMc.getArray()
 525 pressSkin*=1e5
 526 areaSkin=agitateurSkinMc.getMeasureField(True).getArray()
 527 forceSkin=pressSkin*areaSkin
 528 normalSkin=agitateurSkinMc.buildOrthogonalField().getArray()
 529 forceVectSkin=forceSkin*normalSkin
 530 #
 531 singlePolyhedron=agitateurMesh3DMc.buildSpreadZonesWithPoly()
 532 singlePolyhedron.orientCorrectlyPolyhedrons()
 533 centerOfMass=singlePolyhedron.computeCellCenterOfMass()
 534
 535 barySkin=agitateurSkinMc.computeCellCenterOfMass()
 536 posSkin=barySkin-centerOfMass
 537
 538 torquePerCellOnSkin=DataArrayDouble.CrossProduct(posSkin,forceVectSkin)
 539
 540 zeTorque=torquePerCellOnSkin.accumulate()
 541 print("couple = %r N.m"%(zeTorque[2])) ; assert abs(zeTorque[2]-0.37)<1e-2
 542
 543 speedMts=data.getFields()["VITESSE_ELEM_DOM"]
 544 speed1ts=speedMts[(2,-1)]
 545 speedMc=speed1ts.getFieldAtLevel(ON_CELLS,0)
 546 speedOnSkin=speedMc.getArray()[tupleIdsInField]
 547 powerSkin=DataArrayDouble.Dot(forceVectSkin,speedOnSkin)
 548 power=powerSkin.accumulate()[0]
 549 print("power = %r W"%(power)) ; assert abs(power-4.22)<1e-2
 550
 551 x2=posSkin[:,0]*posSkin[:,0] ; x2=x2.accumulate()[0]
 552 y2=posSkin[:,1]*posSkin[:,1] ; y2=y2.accumulate()[0]
 553 xy=posSkin[:,0]*posSkin[:,1] ; xy=xy.accumulate()[0]
 554 inertiaSkin=matrix([[x2,xy],[xy,y2]])
 555 inertiaSkinValues,inertiaSkinVects=linalg.eig(inertiaSkin)
 556 pos=max(enumerate(inertiaSkinValues),key=lambda x: x[1])[0]
 557 vect0=inertiaSkinVects[pos].tolist()[0]
 558 print(vect0)
 559
 560 def computeAngle(locAgitateur1ts):
 561     fMc=locAgitateur1ts.getFieldAtLevel(ON_CELLS,0)
 562     arr=fMc.getArray()
 563     ids=arr.findIdsInRange(0.,1.)
 564     f2Mc=fMc[ids]
 565     m3DSurf,desc,descI,revDesc,revDescI=f2Mc.getMesh().buildDescendingConnectivity()
 566     nbOf3DCellSharing=revDescI.deltaShiftIndex()
 567     ids2=nbOf3DCellSharing.findIdsEqual(1)
 568     agitateurSkinMc=m3DSurf[ids2]
 569     #
 570     singlePolyhedron=agitateurMesh3DMc.buildSpreadZonesWithPoly()
 571     singlePolyhedron.orientCorrectlyPolyhedrons()
 572     centerOfMass=singlePolyhedron.computeCellCenterOfMass()
 573     bary=agitateurSkinMc.computeCellCenterOfMass()
 574     posSkin=bary-centerOfMass
 575     x2=posSkin[:,0]*posSkin[:,0] ; x2=x2.accumulate()[0]
 576     y2=posSkin[:,1]*posSkin[:,1] ; y2=y2.accumulate()[0]
 577     xy=posSkin[:,0]*posSkin[:,1] ; xy=xy.accumulate()[0]
 578     inertiaSkin=matrix([[x2,xy],[xy,y2]])
 579     inertiaSkinValues,inertiaSkinVects=linalg.eig(inertiaSkin)
 580     pos=max(enumerate(inertiaSkinValues),key=lambda x: x[1])[0]
 581     vect0=inertiaSkinVects[pos].tolist()[0]
 582     return vect0
 583
 584 vects=len(ts)*[None]
 585 for itts,locAgitateur1ts in zip(ts,data.getFields()["DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"]):
 586     angle=computeAngle(locAgitateur1ts)
 587     vects[itts[0]]=angle
 588     pass
 589
 590 angle2=len(ts)*[0.]
 591 for pos in range(2, len(vects)):
 592     norm1=sqrt(vects[pos-1][0]*vects[pos-1][0]+vects[pos-1][1]*vects[pos-1][1])
 593     norm2=sqrt(vects[pos][0]*vects[pos][0]+vects[pos][1]*vects[pos][1])
 594     crs=vects[pos-1][0]*vects[pos][0]+vects[pos-1][1]*vects[pos][1]
 595     crs/=norm1 ; crs/=norm2 ; crs=min(crs,1.)
 596     angle2[pos]=acos(crs)#/(ts[pos][2]-ts[pos-1][2])
 597     pass
 598
 599 omega=sum(angle2)/(ts[-1][2]-ts[0][2])
 600 print(sum(angle2)) ; assert abs(sum(angle2)-1.12)<1e-2
 601 print("Au pdt (%d,%d) a %r s le couple est de : %r N.m, power/omega=%r N.m"%(ts[2][0],ts[2][1],ts[2][2],zeTorque[2],power/omega))
 602 assert abs(power/omega-0.37)<1e-2
 603 shutil.rmtree(zeDir)