doc/tutorial/atestMEDCouplingFieldDouble1.rst

   1
   2 .. _python_testMEDCouplingfielddouble1_solution:
   3
   4 Playing with fields
   5 ~~~~~~~~~~~~~~~~~~~
   6
   7 ::
   8
   9         import MEDCoupling as mc
  10
  11         # Create an unstructured mesh from a Cartesian one
  12         xarr = mc.DataArrayDouble.New(11,1)
  13         xarr.iota(0.)
  14         cmesh = mc.MEDCouplingCMesh.New()
  15         cmesh.setCoords(xarr,xarr,xarr)
  16         mesh = cmesh.buildUnstructured()
  17         mesh.convertToPolyTypes(mc.DataArrayInt.Range(0,mesh.getNumberOfCells(),2))
  18         # Create a field
  19         f = mesh.fillFromAnalytic(mc.ON_CELLS,1,"(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)")  # 1 means that the field should have one component
  20         f.setName("MyField")
  21         # A variant:
  22         f2 = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
  23         f2.setMesh(mesh)
  24         f2.setName("MyField2")
  25         f2.fillFromAnalytic(1,"(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)")    # 1 means that the field should have one component
  26         print "Are f and f2 equal?", f.isEqualWithoutConsideringStr(f2,1e-12,1e-12)
  27         #
  28         da1 = f.getArray()              # a DataArrayDouble, which is a direct reference (not a copy) of the field's values
  29         ids1 = da1.getIdsInRange(0.,5.)
  30         fPart1 = f.buildSubPart(ids1)
  31         fPart1.writeVTK("ExoField_fPart1.vtu")
  32         ids2 = f.getArray().getIdsInRange(50.,1.e300)
  33         fPart2 = f.buildSubPart(ids2)
  34         # Renumbering cells to follow MED file rules
  35         fPart1Cpy = fPart1.deepCpy()
  36         o2n = fPart1Cpy.getMesh().sortCellsInMEDFileFrmt()
  37         fPart1Cpy.getArray().renumberInPlace(o2n)
  38         # Check that fPart1Cpy and fPart1 are the same
  39         fPart1Cpy.substractInPlaceDM(fPart1,12,1e-12)
  40         fPart1Cpy.getArray().abs()
  41         print "Are the fields equal?", (fPart1Cpy.getArray().accumulate()[0]<1e-12)
  42         # Aggregate fields
  43         fPart12 = mc.MEDCouplingFieldDouble.MergeFields([fPart1,fPart2])
  44         fPart12.writeVTK("ExoField_fPart12.vtu")
  45         # Evaluation on points
  46         bary = fPart12.getMesh().getBarycenterAndOwner()
  47         arr1 = fPart12.getValueOnMulti(bary)
  48         arr2 = f.getValueOnMulti(bary)
  49         delta = arr1-arr2
  50         delta.abs()
  51         print "Is field evaluation matching?", (delta.accumulate()[0]<1e-12)
  52         # Integral computations
  53         integ1 = fPart12.integral(0,True)
  54         integ1_bis = fPart12.getArray().accumulate()[0]
  55         print "First integral matching ?", ( abs(integ1 - integ1_bis) < 1e-8 )
  56         fPart12.getMesh().scale([0.,0.,0.], 1.2)
  57         integ2 = fPart12.integral(0,True)
  58         print "Second integral matching ?", ( abs(integ2-integ1_bis*1.2*1.2*1.2) < 1e-8 )
  59         # Explosion of field
  60         fVec = mesh.fillFromAnalytic(mc.ON_CELLS,3,"(x-5.)*IVec+(y-5.)*JVec+(z-5.)*KVec")
  61         fVecPart1 = fVec.buildSubPart(ids1)
  62         fVecPart1.setName("fVecPart1")
  63         cells = fPart1.getMesh().getNumberOfCells() * [None]
  64         for icell,vec in enumerate(fVecPart1.getArray()):
  65           m = fPart1.getMesh()[[icell]]
  66           m.zipCoords()      # Not mandatory but saves memory
  67           m.translate(vec)
  68           cells[icell] = m
  69           pass
  70         meshFVecPart1Exploded = mc.MEDCouplingUMesh.MergeUMeshes(cells)
  71         fPart1.setMesh(meshFVecPart1Exploded)
  72         fPart1.writeVTK("ExoField_fPart1_explo.vtu")
  73
  74