src/MEDLoader/Swig/MEDLoaderCouplingTrainingSession.py

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