doc/tutorial/medcouplingremapper_en.rst

   1
   2 Interpolate with MEDCouplingRemapper
   3 ------------------------------------
   4
   5 The purpose of this exercise is to interpolate between two meshes "srcMesh" and "trgMesh".
   6 To make the reader aware of some subtleties about interpolation, a special case
   7 is considered: "srcMesh" is a refinement of "trgMesh".
   8
   9
  10 Implementation start
  11 ~~~~~~~~~~~~~~~~~~~~
  12
  13 To start the exercise import the whole Python module  MEDCouplingRemapper. ::
  14
  15         from MEDCouplingRemapper import *
  16
  17
  18 Create a MEDCouplingUMesh 2D instance, built from a Cartesian mesh
  19 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  20
  21 Build the unstructured mesh "trgMesh" from a 2D Cartesian mesh with 10x10 cells,
  22 starting at point [0.,0.] and having a step of 1.0 along the X and Y directions.
  23 ::
  24
  25         arr=DataArrayDouble(11) ; arr.iota(0)
  26         trgMesh=MEDCouplingCMesh() ; trgMesh.setCoords(arr,arr) ; trgMesh=trgMesh.buildUnstructured()
  27
  28 Create the source mesh "srcMesh"
  29 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  30
  31 Create a mesh "srcMesh" from a 2D Cartesian mesh having 20x20 cells and starting
  32 at point [0.,0.] with a step of 0.5 along the X and Y directions.
  33 ::
  34
  35         arr=DataArrayDouble(21) ; arr.iota(0) ; arr*=0.5
  36         srcMesh=MEDCouplingCMesh() ; srcMesh.setCoords(arr,arr) ; srcMesh=srcMesh.buildUnstructured()
  37
  38 Triangulate the 20 first cells of source mesh (using MEDCouplingUMesh.simplexize()).
  39 Store the result in "srcMesh".
  40 ::
  41
  42         tmp=srcMesh[:20] ; tmp.simplexize(0)
  43         srcMesh=MEDCouplingUMesh.MergeUMeshes([tmp,srcMesh[20:]])
  44
  45 Interpolate using MEDCouplingRemapper
  46 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  47
  48 Compute the first part of the interpolation matrix with the following considerations:
  49 "srcMesh" is regarded as a discretization at cell points, and so is "trgMesh".
  50 To this end invoke prepare() on an instance of the MEDCouplingRemapper class ("remap").
  51 ::
  52
  53         remap=MEDCouplingRemapper()
  54         remap.prepare(srcMesh,trgMesh,"P0P0")
  55
  56 Check that the computed matrix is correct in this trivial case: get the internal
  57 interpolation matrix by calling MEDCouplingRemapper.getCrudeMatrix() and save it in
  58 "myMatrix". This matrix gives for each cell in "trgMesh" the cell IDs of "srcMesh"
  59 which are intersecting the target and the area of intersection. Check that for each
  60 cell in "trgMesh" the sum of the areas is always equal to 1.0.
  61 ::
  62
  63         myMatrix=remap.getCrudeMatrix()
  64         print myMatrix # to see what it looks like
  65         sumByRows=DataArrayDouble(len(myMatrix))
  66         for i,wIt in enumerate(sumByRows):
  67           su=0.
  68           for it in myMatrix[i]:
  69             su+=myMatrix[i][it]
  70           wIt[0]=su
  71         print "Does interpolation look OK? %s"%(str(sumByRows.isUniform(1.,1e-12)))
  72
  73 .. note:: Some triangles were added into "srcMesh" to make "myMatrix" less boring. "myMatrix".
  74
  75 Create a field at cell points "srcField" built from the following analytical formula:
  76 "7-sqrt((x-5.)*(x-5.)+(y-5.)*(y-5.))" where x and y represent the usual space coordinates.
  77 ::
  78
  79         srcField=MEDCouplingFieldDouble(ON_CELLS,ONE_TIME) ; srcField.setMesh(srcMesh)
  80         srcField.fillFromAnalytic(1,"7-sqrt((x-5.)*(x-5.)+(y-5.)*(y-5.))") ; srcField.getArray().setInfoOnComponent(0,"powercell [W]")
  81
  82 .. image:: images/Remapper1.png
  83
  84 Apply the interpolation using MEDCouplingRemapper.transferField(): ::
  85
  86         remap.transferField(srcField,1e300)
  87
  88 .. note:: 1e300 is a default value. It will be systematically assigned to all cells
  89         in "trgField" which do not intersect any cell in "srcMesh". The common usage is to assign
  90         a huge value to identify what is often considered as a bug. But some other use cases
  91         indicate using 0 here, e.g. when considering a parallel interpolation merged ultimately
  92         with an addition.
  93
  94 .. note:: An exception is raised since "srcField" hasn't got an explicitly defined nature.
  95         In what follows the impact of this attribute on the final result will be explained.
  96
  97 Set the nature of "srcField" to IntensiveMaximum (intensive field, e.g. a temperature). ::
  98
  99         srcField.setNature(IntensiveMaximum)
 100         trgFieldCV=remap.transferField(srcField,1e300)
 101
 102 Check that with this nature the field integral is conserved. On the other side
 103 the sum on cells (accumulation) is NOT conserved. ::
 104
 105         print "IntensiveMaximum %lf == %lf"%(srcField.integral(True)[0],trgFieldCV.integral(True)[0])
 106         print "IntensiveMaximum %lf != %lf"%(srcField.getArray().accumulate()[0],trgFieldCV.getArray().accumulate()[0])
 107
 108 Set the nature of "srcField" to ExtensiveConservation (extensive field, e.g. a power). ::
 109
 110         srcField.setNature(ExtensiveConservation)
 111         trgFieldI=remap.transferField(srcField,1e300)
 112
 113 Check that given this nature the field integral is NOT conserved. On the other side the
 114 cumulative sum on cells is conserved. ::
 115 ::
 116
 117         print "ExtensiveConservation %lf != %lf"%(srcField.integral(True)[0],trgFieldI.integral(True)[0])
 118         print "ExtensiveConservation %lf == %lf"%(srcField.getArray().accumulate()[0],trgFieldI.getArray().accumulate()[0])
 119
 120 Visualize the fields using ParaViS.
 121
 122 Solution
 123 ~~~~~~~~
 124
 125 :ref:`python_testMEDCouplingremapper1_solution`