doc/tutorial/medcoupling_fielddouble1_en.rst

   1
   2 Playing with fields
   3 -------------------
   4
   5 In MEDCoupling fields are directly supported by a mesh. That is a major difference
   6 with the field notion in MED file.
   7
   8 Fields instances are well suited:
   9   * for high level services where an interaction with the mesh is requested as in getValueOn(), getValueOnMulti(), integral(), normL1(), normL2(), fillFromAnalytic(), changeUnderlyingMesh(), ...
  10   * to precisely transmit information between components when coupling codes.
  11
  12 For information, the MEDCouplingFieldDouble implementation is tiny since
  13 it delegates most of its work to MEDCouplingMesh, DataArrayDouble, MEDCouplingSpatialDiscretization classes.
  14 MEDCouplingFieldDouble insure the coherency in this model.
  15 It is often possible and even advised to manipulate the array and mesh of a MEDCouplingFieldDouble directly.
  16
  17 Objective
  18 ~~~~~~~~~
  19
  20 This exercise focuses on the relationship between a mesh and the values of a field.
  21
  22 * Create a field
  23 * Aggregate fields
  24 * Build parts of a field
  25 * Renumbering a field
  26 * Compare two fields coming from 2 sources
  27 * Evaluation of a field on a point set
  28 * Explode a field
  29
  30 Implementation start
  31 ~~~~~~~~~~~~~~~~~~~~
  32
  33 Import the MEDCoupling Python module. ::
  34
  35         from MEDCoupling import *
  36
  37 We are going to create a MEDCouplingUMesh from a 3D cartesian mesh. Each direction will contain 10 cells and 11 nodes. The generated MEDCouplingUMesh
  38 will contain 1000 cells. ::
  39
  40         xarr=DataArrayDouble.New(11,1)
  41         xarr.iota(0.)
  42         cmesh=MEDCouplingCMesh.New()
  43         cmesh.setCoords(xarr,xarr,xarr)
  44         mesh=cmesh.buildUnstructured()
  45
  46 In order to put the focus on mixed types, cells with an even id will be converted to polyhedrons ::
  47
  48         mesh.convertToPolyTypes(DataArrayInt.Range(0,mesh.getNumberOfCells(),2))
  49
  50 Creation of a field
  51 ~~~~~~~~~~~~~~~~~~~
  52
  53 Create a field called "MyField" by application of an analytic function "(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)".
  54 Two possibilities :
  55
  56 * Directly by calling fillFromAnalytic on a mesh ::
  57
  58         f=mesh.fillFromAnalytic(ON_CELLS,1,"(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)")
  59         f.setName("MyField")
  60
  61 * Or by building a new instance ::
  62
  63         f2=MEDCouplingFieldDouble.New(ON_CELLS,ONE_TIME)
  64         f2.setMesh(mesh)
  65         f2.setName("MyField2")
  66         f2.fillFromAnalytic(1,"(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)")
  67
  68 Compare the two fields:
  69 Compare f and f2 with a precision of 1e-12 on coordinates and 1e-12 on values. ::
  70
  71         print "f and f2 are equal: %s"%(f.isEqualWithoutConsideringStr(f2,1e-12,1e-12))
  72
  73 Builing of a subpart of a field
  74 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  75
  76 Store in ids1 the list of tuple ids whose value is within [0.0,5.0] (DataArrayDouble.getIdsInRange)     . From ids1 build the sub-part fPart1 of the field "f". ::
  77
  78         ids1=f.getArray().getIdsInRange(0.,5.)
  79         fPart1=f.buildSubPart(ids1)
  80
  81 .. image:: images/FieldDouble1_1.png
  82
  83 Select the part "fPart2" of the field "f" whose values are in [50.,infinity). ::
  84
  85         ids2=f.getArray().getIdsInRange(50.,1.e300)
  86         fPart2=f.buildSubPart(ids2)
  87
  88 Renumbering cells
  89 ~~~~~~~~~~~~~~~~~
  90
  91 The generated file "fPart1" is valid for MEDCoupling, but its cells are not sorted by geometric type: it is not valid from a MED file point of view. By using MEDCouplingUMesh.sortCellsInMEDFileFrmt and DataArrayDouble.renumberInPlace
  92 renumber manually fPart1 starting from a deep copy of fPart1. ::
  93
  94         fPart1Cpy=fPart1.deepCpy()
  95         o2n=fPart1Cpy.getMesh().sortCellsInMEDFileFrmt()
  96         fPart1Cpy.getArray().renumberInPlace(o2n)
  97
  98 "fPart1Cpy" is now normalized to be stored in a MED file (we will tackle this later).
  99 Check that fPart1Cpy and fPart1 are the same (discarding any permutation): ::
 100
 101         fPart1Cpy.substractInPlaceDM(fPart1,12,1e-12)
 102         fPart1Cpy.getArray().abs()
 103         print "Fields are the same? %s"%(fPart1Cpy.getArray().accumulate()[0]<1e-12)
 104
 105 .. note:: This is in fact a very special case of interpolation. Except that here
 106         we assume that the supports of "fPart1" and "fPart1Cpy" are equal, discarding any
 107         cell and/or node permutation.
 108
 109 Aggregate Fields
 110 ~~~~~~~~~~~~~~~~
 111
 112 Aggregate fields "fPart1" and "fPart2". The result is stored in "fPart12". ::
 113
 114         fPart12=MEDCouplingFieldDouble.MergeFields([fPart1,fPart2])
 115
 116 .. image:: images/FieldDouble1_2.png
 117
 118 .. note:: Apologies for the name MEDCouplingFieldDouble.MergeFields instead of
 119                 AggregateFields.
 120
 121 Evaluation of a MEDCouplingFieldDouble on given space points
 122 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 123
 124 Evaluate the values of the computed field "fPart12" on the barycenters of its mesh.
 125 Evaluate the field "f" on the same barycenters. The method used is MEDCouplingFieldDouble.getValueOnMulti(). ::
 126
 127         bary=fPart12.getMesh().getBarycenterAndOwner()
 128         arr1=fPart12.getValueOnMulti(bary)
 129         arr2=f.getValueOnMulti(bary)
 130         delta=arr1-arr2
 131         delta.abs()
 132         print "Check OK: %s"%(delta.accumulate()[0]<1e-12)
 133
 134 .. note:: In this context and for example for a field on cells, "evaluate" at a point means returning the value of the cell containing the point.
 135 .. note:: This technique can be used to quickly assess the quality of an interpolation.
 136
 137 Operations on a field
 138 ~~~~~~~~~~~~~~~~~~~~~
 139
 140 Compute the integral of the field "fPart12" and compute it a second time by using
 141 DataArrayDouble.accumulate on the underlying DataArrayDouble of this "fPart12" (remember that the cell volumes are all 1.0).
 142 To show the link with the underlying mesh, scale the underlying mesh (fPart12.getMesh()) by 1.2 and centered at [0.,0.,0.].
 143 Recompute the integral.
 144 ::
 145
 146         fPart12.integral(0,True)
 147         fPart12.getArray().accumulate()
 148         fPart12.getMesh().scale([0.,0.,0.],1.2)
 149         abs(fPart12.integral(0,True)-fPart12.getArray().accumulate()[0]*1.2*1.2*1.2)<1e-8
 150
 151 Exploding a field
 152 ~~~~~~~~~~~~~~~~~
 153
 154 Starting from "mesh", create a vector field on cells "fVec" with 3 components representing the displacement between each cell's barycenter and the point [5.,5.,5.]. Use MEDCouplingMesh.fillFromAnalytic(). ::
 155
 156         fVec=mesh.fillFromAnalytic(ON_CELLS,3,"(x-5.)*IVec+(y-5.)*JVec+(z-5.)*KVec")
 157
 158 Create the reduction of "fVec" ("fVecPart1") on cell IDs "ids1" (previously obtained). ::
 159
 160         fVecPart1=fVec.buildSubPart(ids1)
 161         fVecPart1.setName("fVecPart1")
 162
 163 Build the scalar field fPart1Exploded having the same values as "fPart1" but supported by an exploded mesh (in comparison to fPart1.getMesh()).
 164 To explode the underlying mesh fPart1.getMesh(), use the vectorial displacement field "fVecPart1" in order to apply to each cell the proper translation. ::
 165
 166         cells=fPart1.getMesh().getNumberOfCells()*[None]
 167         for icell,vec in enumerate(fVecPart1.getArray()):
 168           m=fPart1.getMesh()[[icell]]
 169           m.zipCoords() # pas absolument nécessaire mais permet d'être économe en mémoire
 170           m.translate(vec)
 171           cells[icell]=m
 172           pass
 173         meshFVecPart1Exploded=MEDCouplingUMesh.MergeUMeshes(cells)
 174         fPart1.setMesh(meshFVecPart1Exploded)
 175
 176 .. image:: images/FieldDouble1_1_exploded.png
 177
 178 Solution
 179 ~~~~~~~~
 180
 181 :ref:`python_testMEDCouplingfielddouble1_solution`