2 \page interpkernel Interpolation kernel toolkit
4 \section InterpKerIntro Introduction
6 The main purpose of this module is to propose a set of algorithms for
7 mesh interpolation \b fully \b independant \b of \b mesh \b datastructure to
8 support several type of format. This component is parameterized as
9 much as possible using C++ templates.
10 For the moment only interpolators for unstructured meshes are present in
11 %interpolation kernel.
13 \section InterpKerTheory Theory of interpolation
15 \subsection InterpKerPerfOverl Mesh overlapping
17 When interpolation is performed between a source mesh S and a target
18 mesh T the aspect of overlapping is important. In fact if any cell of
19 of S is fully overlapped by cells of T and inversely any cell of T is
20 fully overlapped by cells of S the meshes S and T are said \b
21 coincidant and some general formulae in next sub section are simpler.
22 As far as possible in the next sub sections the formulae for
23 coincidant and non-coincidant meshes will be given.
25 \subsection InterpKerRemapGlobal Global conservative remapping
27 For fields with polynomial representation on each cell, the components of the discretized field \f$ \phi_s \f$ on the source side can be expressed as linear combinations of the components of the discretized field \f$ \phi_t \f$ on the target side, in terms of a matrix-vector product :
33 All the aim of interpolators is to compute W depending on a physical
34 quantities and the type of interpolation wanted (P0, P1, P1d etc...).
36 For the \b intensive \b field \f$ \phi \f$ from the
37 source mesh \f$ S \f$ to the target mesh \f$ T \f$ the interpolation should preserve the
38 integral of \f$ \phi \f$ on any domain. At the discrete level, for any
39 target cell \f$ T_i \f$, the following \b general \b interpolation \b
40 equation \anchor InterpKerGenralEq has to
44 \int_{T_i} \phi = \sum_{S_j\cap T_i \neq \emptyset} \int_{T_i\cap S_j} \phi.
47 To compute \f$ W_{ij} \f$
48 this equation is used. The evaluation of integrals depends on the source and target meshes and on the nature of interpolation chosen : P0, P1, P1d etc... For the moment it is only possible to
49 remap fields with P0 representations.
51 \subsection InterpKerRemapInt Conservative remapping of intensive physical quantities
53 At the basis of many CFD numerical schemes is the fact that physical
54 quantities such as density, momentum per unit volume or energy per
55 unit volume obey some balance laws that should be preserved at the
56 discrete level on every cell. This property is critical for example to
57 accurately capture shockwaves.
59 \subsection InterpKerP0P0Int cell-cell (P0->P0) conservative remapping of intensive physical quantities
61 In the \ref InterpKerGenralEq "general interpolation equation" the
62 left hand side becomes :
65 \int_{T_i} \phi = (\sum_{S_j} Vol(T_i\cap S_j)).\phi_{T_i}.
68 \note \f$ \sum_{S_j} Vol(T_i\cap S_j) = Vol(T_i) \f$ \ref InterpKerPerfOverl "in case of perfect overlapping".
70 In the \ref InterpKerGenralEq "general interpolation equation" the
71 right hand side becomes :
74 \sum_{S_j\cap T_i \neq \emptyset} \int_{T_i\cap S_j} \phi = \sum_{S_j\cap T_i \neq \emptyset} {Vol(T_i\cap S_j)}.\phi_{S_j}.
77 In the case where the \b intensive field values are constant on each
78 cell, the coefficients of the linear remapping matrix \f$ W \f$ are
79 given by the formula :
82 W_{ij}=\frac{Vol(T_i\cap S_j)}{ \sum_{S_j} Vol(T_i\cap S_j) }.
85 and \ref InterpKerPerfOverl "in case of perfect overlapping" :
88 W_{ij}=\frac{Vol(T_i\cap S_j)}{ Vol(T_i) }.
91 Where Vol represents the volume with mesh dimension of interpolation equals to 3, the
92 area when mesh dimension equals to 2, and length when mesh dimension equals to 1.
94 \subsection InterpKerRemapExt Conservative remapping of extensive physical quantities
96 In code coupling from neutronics to hydraulic code \b extensive field
97 of power is exchanged and the all power as to be kept the same. The
98 principle is to 'intensify' the field to move on from extensive field
99 \e P to an intensive one \f$ \phi \f$.
101 \subsection InterpKerP0P0Ext cell-cell (P0->P0) conservative remapping of extensive physical quantities
103 In the \ref InterpKerGenralEq "general interpolation equation" the
104 left hand side becomes :
107 \int_{T_i} \phi = P_{T_i}.
110 In the \ref InterpKerGenralEq "general interpolation equation" the
111 right hand side becomes :
114 \sum_{S_j\cap T_i \neq \emptyset} \int_{T_i\cap S_j} \phi =
115 \sum_{S_j\cap T_i \neq \emptyset} \frac{Vol(T_i\cap S_j)}{\sum_{T_i} Vol(T_i \cap S_j)}.P_{S_j}.
118 \note \f$ \sum_{T_i} Vol(T_i \cap S_j) = Vol(S_j) \f$ \ref InterpKerPerfOverl "in case of perfect overlapping".
120 In the case where the \b extensive field values are constant on each
121 cell, the coefficients of the linear remapping matrix \f$ W \f$ are
122 given by the formula :
125 W_{ij}=\frac{Vol(T_i\cap S_j)}{ \sum_{T_i} Vol(T_i \cap S_j) }.
128 and \ref InterpKerPerfOverl "in case of perfect overlapping" :
131 W_{ij}=\frac{Vol(T_i\cap S_j)}{ Vol(S_j) }.
134 \section InterpKerMainArchitecture Main architecture of interpolation kernel.
136 In %interpolation kernel, the algorithm that computes \f$ W_{ij} \f$ given cell i
137 in source mesh and cell j in target mesh is called intersector.
139 As seen in \ref InterpKerTheory "the theory of interpolation", for all interpolation the aim is to
140 fill the W matrix (which is generally a sparse matrix). Fundamatally for each pair (i,j) \f$ W_{ij} \f$ is obtained
141 by calling each time the wanted intersector. The problem is that each call to algorithm
143 To reduce the computation time a first filtering is done to found a
144 maximim of (i,j) pairs where \f$ W_{ij} \f$ is obviously equal to 0. Typically it
145 is the case when a cell in source mesh is too far from an another cell
146 in target mesh each other.
148 So for a given type of interpolation, the computation of W is
149 performed in two steps :
151 -# A filtering process reduces the number of pairs of
152 elements for which the calculation must be carried out by
153 eliminating pairs that do not intersect through a comparison of
154 their bounding boxes. It reduces as much as possible call to intersector.
155 -# For all remaining pairs calling for each intersector (see \ref interpolation2D, \ref interpolation3Dsurf or \ref interpolation3D).
157 Each interpolator inherits from INTERP_KERNEL::Interpolation ( whatever
158 its dimension and its type ) that is a
159 CRTP class in order to clearly see the main API without useless CPU cost.
161 \subsection InterpKerMeshType class MeshType
163 Each Interpolators and Intersectors are parameterized (templated in
164 C++ langage) with \c class \c MeshType . This type of generalization
165 has been chosen to reduce at maximum overhead. \n
166 Thanks to this principle \b intersectors \b and \b interpolators \b are \b usable
167 \b with \b several \b formats \b without \b preformance \b loss. For example MED, VTK...\n
168 \c MeshType is a concept that should strictly fulfilled the following
171 - Const values / Types
172 - MyConnType : represents type of connectivity index. This is typically \c int or \c long \c int .
173 - MY_SPACEDIM : space dimension. Dimension relative to coordinates.
174 - MY_MESHDIM : the dimension of all cells in meshes.
175 - My_numPol : policy of numbering. C Format ( \f$ [0,n-1] \f$ ) or FORTRAN ( \f$ [1,n] \f$ ).
177 -# \code void getBoundingBox(double *boundingBox) const \endcode
178 -# \code INTERP_KERNEL::NormalizedCellType getTypeOfElement(MyConnType eltId) const \endcode
179 -# \code unsigned char getNumberOfNodesOfElement(MyConnType eltId) const \endcode
180 -# \code unsigned long getNumberOfNodes() const \endcode
181 -# \code unsigned long getNumberOfElements() const \endcode
182 -# \code const MyConnType *getConnectivityPtr() const \endcode
183 -# \code const double *getCoordinatesPtr() const \endcode
184 -# \code const MyConnType *getConnectivityIndexPtr() const \endcode
185 -# \code void ReleaseTempArrays() \endcode
187 - the array returned by \c getCoordinatesPtr must be a \b full \b interlace array.
188 - the arrays returned by \c getConnectivityPtr and \c
189 - getConnectivityIndexPtr must be with the same principle as it is in \ref medmemConnArrays "medmem". Of course the numbering format may change according to \a My_numPol policy.
191 \note The arrays formats of connectivity is kept close to MED. It is
192 close to VTK too but slightly different. So it needs VTK side a copy
193 on wrap. To avoid this copy of a part of connectivity structure, iterator should be used.
195 \subsection InterpKerMatrixType class MatrixType
197 As already said, the matrix returned by interpolator is typically a sparse matrix. Instances of
198 \c class \c MatrixType are used to stores these results of
199 interpolation. To be able to be filled by the interpolator the \c MatrixType class has to match the following concept :
202 -# \code void resize(uint nbrows) \endcode
203 -# \code Row &operator [] (uint irow) \endcode
205 \c class \c Row has to match at least following concept :
208 - \code void insert(const std::pair<int,T>& myPair) \endcode
210 \note \c std::vector\c < \c std::map<int,double> > is a candidate for
213 \section InterpKerGenUsage Usage of interpolation kernel.
215 \subsection InterpKerHighLevUsage high-level usage
217 This the simplest mode of usage of interpolator. This way is strongly
218 linked with MED data-structure. All interpolators is completely hidden
219 to you. Even sparse interpolation matrix is hidden. An exemple of
224 std::string sourceFileName("source.med");
225 MEDMEM::MESH med_source_mesh(MED_DRIVER,sourceFileName,"Source_Mesh");
226 std::string targetFileName("target.med");
227 MEDMEM::MESH med_target_mesh(MED_DRIVER,targetFileName,"Target_Mesh");
228 FIELD<double> sourceField(MED_DRIVER,sourceFileName,"Density",0,0);
229 FIELD<double> targetField;
231 mapper.prepare(med_source_mesh,med_target_mesh);
232 mapper.transfer(sourceField,targetField);
237 \subsection InterpKerMidLevUsage middle-level usage
239 This mode is the mode that needs the minimum of prerequisites
240 (algorithms and the datastructure you intend to use). On the other
241 hand it is needed to specify precisely nature of interpolator.
243 As consequence of genericity of interpolators, they are usable only by
244 instanciating an underneath mesh data structure. The two following
245 examples show how to use interpolator at this level.
247 - The simplest way to use the interpolator with MED datastructure is :
251 std::string sourceFileName("source.med");
252 MEDMEM::MESH med_source_mesh(MED_DRIVER,sourceFileName,"Source_Mesh");
253 std::string targetFileName("target.med");
254 MEDMEM::MESH med_target_mesh(MED_DRIVER,targetFileName,"Target_Mesh");
255 // Ok at this point we have our mesh in MED-Memory format.
256 // Go to wrap med_source_mesh and med_target_mesh.
257 MEDNormalizedUnstructuredMesh<2,2> wrap_source_mesh(&med_source_mesh);
258 MEDNormalizedUnstructuredMesh<2,2> wrap_target_mesh(&med_target_mesh);
259 // Go for interpolation...
260 INTERP_KERNEL::Interpolation2D myInterpolator;
261 //optionnal call to parametrize your interpolation. First precision, tracelevel, intersector wanted.
262 myInterpolator.setOptions(1e-7,0,Geometric2D);
263 INTERP_KERNEL::Matrix<double,ALL_FORTRAN_MODE> resultMatrix;
264 myInterpolator.interpolateMeshes(wrap_source_mesh,wrap_target_mesh,resultMatrix,"P0P0");
265 //Ok let's multiply resultMatrix by source field to interpolate to target field.
266 resultMatrix.multiply(...)
270 - Same with VTK datastructure :
274 vtkXMLUnstructuredGridReader *readerSource=vtkXMLUnstructuredGridReader::New();
275 readerSource->SetFileName("source.vtu");
276 vtkUnstructuredGrid *vtk_source_mesh=readerSource->GetOutput();
277 readerSource->Update();
278 vtkXMLUnstructuredGridReader *readerTarget=vtkXMLUnstructuredGridReader::New();
279 readerTarget->SetFileName("target.vtu");
280 vtkUnstructuredGrid *vtk_target_mesh=readerTarget->GetOutput();
281 readerTarget->Update();
282 // Ok at this point we have our mesh in VTK format.
283 // Go to wrap vtk_source_mesh and vtk_target_mesh.
284 VTKNormalizedUnstructuredMesh<2> wrap_source_mesh(vtk_source_mesh);
285 VTKNormalizedUnstructuredMesh<2> wrap_target_mesh(vtk_target_mesh);
286 // Go for interpolation...
287 INTERP_KERNEL::Interpolation2D myInterpolator;
288 //optionnal call to parametrize your interpolation. First precision, tracelevel, intersector wanted.
289 myInterpolator.setOptions(1e-7,0,Geometric2D);
290 INTERP_KERNEL::Matrix<double,ALL_C_MODE> resultMatrix;
291 myInterpolator.interpolateMeshes(wrap_source_mesh,wrap_target_mesh,resultMatrix,"P0P0");
292 //Ok let's multiply resultMatrix by source field to interpolate to target field.
293 resultMatrix.multiply(...)
295 readerSource->Delete();
296 readerTarget->Delete();