doc/user/input/data_conversion.rst

   1 Data conversion
   2 ===============
   3
   4 medcoupling includes MEDCouplingRemapper object that allows given a
   5 source field and a target mesh to produce a field lying on target
   6 mesh.
   7
   8 MEDCouplingRemapper goal is to compute a projection by minimizing
   9 diffusive effects and conserve at most the magnitude of key
  10 variables.
  11
  12 The MEDCouplingRemapper object computes a remap function that allows
  13 to transform source field into target field.
  14
  15 The function is represented using sparse matrix given:
  16
  17   - a source mesh
  18   - a target mesh
  19   - a source spatial discretization ("P0" means ON_CELLS and "P1" means ON_NODES)
  20   - a target spatial discretization
  21
  22 Apply a projection is then equivalent to apply matrix vector
  23 multiply.
  24
  25 As projection requires mathematically to compute for every cell in
  26 target mesh, the contribution relative to each of cell in source
  27 mesh. So projection is computation-intensive task.
  28
  29 To optimize time spent by MEDCouplingRemapper, it implements a fast
  30 algorithm that select coarsely the pairs of cells that have a chance
  31 to be in interaction each other.
  32 It allows to drop from O(n^2) to O(n) the complexity of algorithm.
  33 The algorithm is based on rectilinear bounding boxes of cells along
  34 axis. 2 cells in interaction have their respective bounding boxes in
  35 interaction.
  36
  37 Unfortunately it works well mathematically but in the real life
  38 projections especially when spacedim and meshdim are not the same it
  39 can be harder to set.
  40
  41 .. _projection_P0P0:
  42
  43 Projection cells/cells with meshes having same mesh dimension/space dimension
  44 -----------------------------------------------------------------------------
  45
  46 Let us create meshes for sample projection: a source mesh consisting of two quadrangles and a target mesh consisting of a triangle:
  47
  48 .. figure:: ../images/projection_P0P1_meshes.png
  49    :align: center
  50
  51    Two source quadrangles in blue and a target triangle in red
  52
  53 .. literalinclude:: ../../../src/MEDCoupling_Swig/UsersGuideExamplesTest.py
  54    :start-after: UG_Projection_0
  55    :end-before:  UG_Projection_0
  56
  57 Now perform projection:
  58
  59 .. literalinclude:: ../../../src/MEDCoupling_Swig/UsersGuideExamplesTest.py
  60    :start-after: UG_Projection_1
  61    :end-before:  UG_Projection_1
  62
  63 Here "P0P0" specifies spatial discretization of both source and target field. "P0" is for cell discretization (ON_CELLS). "P1" is for node discretization (ON_NODES). For example, "P0P1" means from a source cell field (lying on a source mesh) to a target node field (lying on a target mesh).
  64
  65 Remapper *rem* has just performed intersection of *src* and *tgt* meshes. Result of intersection is stored in a matrix which is returned by getCrudeMatrix method. Each element of the matrix is size of intersection of two cells i and j of the source and target meshes correspondingly.
  66
  67 Generally the matrix is sparse. To reduce memory usage, the space matrix is represented by an array whose i-th element is a map storing indices and intersection sizes of source cells intersecting i-th target cell.
  68
  69 Here is a sample sparse matrix:
  70
  71 .. math::
  72
  73    \begin{bmatrix}
  74        11 & 0  & 0  & 0  \\
  75        0  & 0  & 32 & 0  \\
  76        0  & 0  & 33 & 43 \\
  77        14 & 24 & 0  & 0  \\
  78        0  & 25 & 0  & 45
  79    \end{bmatrix}
  80
  81 and a corresponding array, which stores only non-zero matrix elements:
  82
  83 ::
  84
  85   [{0: 11, 3: 14}, {3: 24, 4: 25}, {1: 32, 2: 33}, {2: 43, 4: 45}]
  86
  87 Now we create a field on all cells of *src* mesh and call rem.transferField() to transfer the field to the cells of *tgt* mesh
  88
  89 .. literalinclude:: ../../../src/MEDCoupling_Swig/UsersGuideExamplesTest.py
  90    :start-after: UG_Projection_2
  91    :end-before:  UG_Projection_2
  92
  93 The last argument -1 is a default value that will be assigned in *trgF* to each entity of *tgt* mesh that is not intercepted by any entity of *src* mesh.
  94
  95 It is crucial to setNature() of the field by which an adequate interpolation formula is chosen. The field nature is expressed by one of the following values:  ExtensiveMaximum, IntensiveMaximum, ExtensiveConservation, IntensiveConservation.
  96
  97 The field is represented by a vector with a discrete value on each cell. This value can represent either
  98
  99 * an average value of the field in the cell (average density, velocity or temperature in the cell) in which case the representation is said to be **intensive**,
 100 * an integrated value over the cell (total mass, power of the cell) in which case the representation is said to be **extensive**.
 101
 102 In the case where the source and target meshes are not fully overlapping, it is important to make a correct choice between Maximum and Conservation.
 103
 104 * **Conservation** means that the interpolation algorithm preserves the integral of the field over a domain.
 105 * **Maximum** means that the field values resulting from the interpolation remain between the upper and lower bounds of the original field.
 106
 107 Projection cells/nodes nodes/cells nodes/nodes with meshes having same mesh dimension/space dimension
 108 -----------------------------------------------------------------------------------------------------
 109
 110 Suppose we have the same source and target meshes as in the :ref:`previous chapter <projection_P0P0>`. I.e. the source mesh consists of two quadrangles and the target mesh consists of a triangle.
 111
 112 .. figure:: ../images/projection_P0P1_meshes.png
 113    :align: center
 114
 115    Two source quadrangles in blue and a target triangle in red
 116
 117 The source field is also on cells but the target field is on **nodes** of target mesh. In such a case we do the following:
 118
 119 .. literalinclude:: ../../../src/MEDCoupling_Swig/UsersGuideExamplesTest.py
 120    :start-after: UG_Projection_3
 121    :end-before:  UG_Projection_3
 122
 123 Here "P0P1" specifies support of the source (P0=ON_CELLS) and the target (P1=ON_NODES) fields. The computed martix is::
 124
 125   [{0: 0.125, 1: 0.04166666666666663}, {1:0.16666666666666669}, {1: 0.041666666666666685}]
 126
 127 The result matrix contains values for each of the three nodes of the target mesh. In order to be able to compute these values by means of intersection of cells, an algorithm constructs an auxiliary target mesh. This auxiliary mesh is a dual mesh of the actual target mesh, hence it includes cells corresponding to each actual target node.
 128
 129 .. figure:: ../images/projection_P0P1_dual.png
 130    :align: center
 131
 132    The target mesh used for projection is a dual one of the actual target mesh
 133
 134 If the case is reverse, i.e. the source field is on nodes and the target one should be on cells, we change "P0P1" to "P1P0":
 135
 136 .. literalinclude:: ../../../src/MEDCoupling_Swig/UsersGuideExamplesTest.py
 137    :start-after: UG_Projection_4
 138    :end-before:  UG_Projection_4
 139
 140 The computed martix is::
 141
 142   [{1: 0.07291666666666663, 4: 0.375}]
 143
 144 In the same way as in "P0P1" case, an auxiliary dual mesh is internally constructed for a mesh supporting a field on nodes (the source field in our case). But since the source cells are not simplexes, simplexization of quadrangles is perfrormed on the fly to be able to construct the dual mesh. Cell #0 (0,3,4,1) is split into (0,3,4) and (0,4,1)
 145
 146 .. figure:: ../images/projection_P1P0_dual.png
 147    :align: center
 148
 149    The source mesh used for projection is a dual one of simplexized actual source mesh
 150
 151 .. |br| raw:: html
 152
 153   <br />
 154
 155 .. note:: in present implementation of Remapper, the projection of fields on nodes lying on quadratic meshes is not available. If you are in this case you have 2 choices :|br|
 156    1) you accept to lose precision by converting your mesh (MEDCouplingUMesh.convertLinearCellsToQuadratic) |br|
 157    2) you subdivide your mesh.
 158
 159 Projection cell/cell when meshdim != spacedim
 160 ---------------------------------------------
 161
 162 When meshdim is not equal to spacedim some additional important parameters should be considered.
 163
 164 Consider a simple case with spacedim=3 and meshdim=2 where :
 165
 166  - source mesh contains one triangle cell [0.,0.,0.], [1,0,0], [0,1,0]
 167  - target mesh also contains one triange cell [0.,0.,0.], [1,0,0], [0,1,0], which coincides with the source triangle.
 168
 169 Only "P0P0" projection (ON_CELLS->ON_CELLS) will be considered here but all elements presented here are applicable for other spatial discretizations projections.
 170
 171 .. literalinclude:: ../../../src/MEDCoupling_Swig/UsersGuideExamplesTest.py
 172    :start-after: UG_Projection_5
 173    :end-before:  UG_Projection_5
 174
 175 A remapper returns the following correct matrix::
 176
 177   [{0: 0.5}]
 178
 179 where 0.5 is equal to area of the source (or target) triangle.
 180
 181 Let us see what the remapper will return if the meshes are not overlapping. To this end we translate  the source mesh along its normal by a small distance (1e-3):
 182
 183 .. literalinclude:: ../../../src/MEDCoupling_Swig/UsersGuideExamplesTest.py
 184    :start-after: UG_Projection_6
 185    :end-before:  UG_Projection_6
 186
 187 The remapper returns the following matrix, which means that it have not found any intersection of the meshes::
 188
 189   [{}]
 190
 191 This happens because bounding boxes of cells used to detect possibly intersecting cells do not intersect. In order to get a correct result of intersection we can increase those bounding boxes. The remapper provides two ways to do that:
 192
 193 - rem.setBoundingBoxAdjustmentAbs( *thickness* ) expands a bounding box by *thickness*.
 194 - rem.setBoundingBoxAdjustment( *factor* ) expands a bounding box by thickness computed by multiplying a maximal bounding box size by *factor*.
 195
 196 .. figure:: ../images/projection_bnd_box.png
 197    :align: center
 198
 199    If source and target meshes do not overlap we can increase their bounding boxes to help remapper to detect intersecting cells
 200
 201 .. literalinclude:: ../../../src/MEDCoupling_Swig/UsersGuideExamplesTest.py
 202    :start-after: UG_Projection_7
 203    :end-before:  UG_Projection_7
 204
 205 With bounding box tuning we obtain the correct result matrix::
 206
 207   [{0: 0.5}]
 208
 209 The bounding box tuning has a disadvantage in the case where a cell has large angle with either of planes OXY, OYZ or OXZ. In such a case cells that are rather far one from another can be detected as possibly intersecting since their bounding boxes are too large. setMaxDistance3DSurfIntersect method helps avoiding excess intersections. If a distance between fast barycenter of target cell and medium source plane is more than the value specified via setMaxDistance3DSurfIntersect then cell intersection is not performed.
 210
 211 .. figure:: ../images/projection_setMaxDistance3DSurfIntersect.png
 212    :align: center
 213
 214    If the tuned bounding box is too large we can use setMaxDistance3DSurfIntersect to exclude unnecessary intersections
 215
 216 Let us rotate the source mesh by pi/4 angle around Y axis. The result matrix is
 217
 218 .. literalinclude:: ../../../src/MEDCoupling_Swig/UsersGuideExamplesTest.py
 219    :start-after: UG_Projection_8
 220    :end-before:  UG_Projection_8
 221
 222 ::
 223
 224   [{0: 0.46155716207961217}]
 225
 226 If we use setMaxDistance3DSurfIntersect, the result matrix shows that no intersection found:
 227
 228 .. literalinclude:: ../../../src/MEDCoupling_Swig/UsersGuideExamplesTest.py
 229    :start-after: UG_Projection_9
 230    :end-before:  UG_Projection_9
 231
 232 ::
 233
 234   [{}]
 235
 236 Another function useful to avoid unnecessary intersection is setMinDotBtwPlane3DSurfIntersect. It tells the remapper not to intersect cells whose normals are too much different. setMinDotBtwPlane3DSurfIntersect specifies minimal dot product of unitary normal vectors of a source and target cells.
 237
 238 Let us switch off setMaxDistance3DSurfIntersect by passing a negative value. And use setMinDotBtwPlane3DSurfIntersect to avoid intersecting cells that has 45 degrees angle between their normals in our case.
 239
 240 .. literalinclude:: ../../../src/MEDCoupling_Swig/UsersGuideExamplesTest.py
 241    :start-after: UG_Projection_10
 242    :end-before:  UG_Projection_10
 243
 244 Result matrix is::
 245
 246   [{}]