Version: 6.3.1

Interpolation tools

Introduction

The main purpose of this module is to propose a set of algorithms for mesh interpolation fully independant of the mesh datastructure to support several type of format. This component is parameterized as much as possible using C++ templates. For the moment only interpolators for unstructured meshes are present in the interpolation kernel.

Main architecture of interpolation kernel.

In the interpolation kernel, algorithms that computes the intersection $ T_i\cap S_j$ given the locations and geometries of source cell $ S_j $ and target cell $ T_i $ are called Intersectors.

As can be seen in the theory of interpolation, all the proposed interpolators aim at filling the interpolation matrix W (which is generally sparse). For each pair (i,j), $ W_{ij} $ is obtained by calling the desired intersector. The problem is that each call to this algorithm is CPU-expensive. To reduce the computational time, a first filtering is done to detect pairs (i,j) $ W_{ij} $ is obviously equal to 0. It is typically the case when a cell in the source mesh is too far from an another cell in the target mesh each.

So for a given type of interpolation, the computation of W is performed in two steps :

  1. A filtering process reduces the number of pairs of elements for which the calculation must be carried out by eliminating the pairs whose bounding boxes do not intersect.
  2. For all remaining pairs calling for each intersector (see Special features of 2D intersectors, Special features of 3D surface intersectors or Special features of 3D volumes intersectors).

Whatever its dimension and type, each interpolator inherits from INTERP_KERNEL::Interpolation which is a template (CRTP) class than enable an easy access to the main API without useless CPU cost.

class MeshType

Each Interpolators and Intersectors are parameterized (templated in C++ langage) with class MeshType . This type of generalization has been chosen to reduce at maximum overhead.
Thanks to this principle intersectors and interpolators are usable with mesh formats such as MED or VTK, without preformance loss. MeshType is a concept that should strictly fulfilled the following rules :

  • Const values / Types
    • MyConnType : represents type of connectivity index. This is typically int or long int .
    • MY_SPACEDIM : space dimension. Dimension relative to coordinates.
    • MY_MESHDIM : the dimension of all cells in meshes.
    • My_numPol : policy of numbering. C Format ( $ [0,n-1] $ ) or FORTRAN ( $ [1,n] $ ).
  • Methods
    1.  void getBoundingBox(double *boundingBox) const 

    2.  INTERP_KERNEL::NormalizedCellType getTypeOfElement(MyConnType eltId) const 

    3.  unsigned char getNumberOfNodesOfElement(MyConnType eltId) const 

    4.  unsigned long getNumberOfNodes() const 

    5.  unsigned long getNumberOfElements() const 

    6.  const MyConnType *getConnectivityPtr() const 

    7.  const double *getCoordinatesPtr() const 

    8.  const MyConnType *getConnectivityIndexPtr() const 

    9.  void releaseTempArrays()
  • Formats of arrays
    • the array returned by getCoordinatesPtr must be a full interlace array.
    • the arrays returned by getConnectivityPtr and getConnectivityIndexPtr must be with the same principle as it is in medmem. Of course the numbering format may change according to My_numPol policy.

Note that the array format for connectivity is kept close to MED. It is close to VTK format too but slightly different. So it may require for the VTK side a copy on wrap. To avoid this copy of a part of the connectivity structure, an iterator should be used.

class MatrixType

As already said, the matrix returned by interpolator is typically a sparse matrix. Instances of class MatrixType are used to store the resulting interpolation matrix. To be able to be filled by the interpolator the MatrixType class has to match the following concept :

  • Methods
    1.  void resize(uint nbrows) 

    2.  Row &operator [] (uint irow)

class Row has to match at least the following concept :

  • Methods
    •  void insert(const std::pair<int,T>& myPair)

Note that std::vector< std::map<int,double> > is a candidate for MatrixType.

Usage of interpolation tools: the REMAPPER classes.

high-level usage

The simplest way of using the interpolation tools is in sequential mode to use the REMAPPER classes. These classes fulfill HXX2SALOME rules and may be used in coupling graphs. Two sequential REMAPPERS exist, ParaMEDMEM::MEDCouplingRemapper and MEDMEM::MEDMEM_REMAPPER . These classes are strongly linked to their corresponding data structure, respectively MEDCoupling and MEDMEM.

...
const char sourceFileName[]="source.med";
MEDCouplingFieldDouble *sourceField=MEDLoader::ReadFieldCell(sourceFileName,"Source_Mesh",0,"Density",/*iteration*/0,/*order*/0);
const char targetFileName[]="target.med";
MEDCouplingUMesh *med_target_mesh=MEDLoader::ReadUMeshFromFile(targetFileName,"Target_Mesh",0);
//
sourceField->setNature(ConservativeVolumic);//Specify nature is needed to allow remapper object to apply correct policy for denominator computation !
MEDCouplingRemapper remapper;
remapper.setPrecision(1e-12);
remapper.setIntersectionType(INTERP_KERNEL::Triangulation);
remapper.prepare(sourceField->getMesh(),med_target_mesh,"P0P0");
MEDCouplingFieldDouble *targetField=remapper.transferField(sourceField,/*default_value*/4.57);//Any target cell not intercepted by any source cell will have value set to 4.57.
...
// clean-up
targetField->decrRef();
sourceField->decrRef();
med_target_mesh->decrRef();
  • If you intend to use MEDMEM data structure, medmemremapper class should be used :
...
std::string sourceFileName("source.med");
MEDMEM::MESH med_source_mesh(MED_DRIVER,sourceFileName,"Source_Mesh");
std::string targetFileName("target.med");
MEDMEM::MESH med_target_mesh(MED_DRIVER,targetFileName,"Target_Mesh");
FIELD<double> sourceField(MED_DRIVER,sourceFileName,"Density",0,0);
FIELD<double> targetField;
Remapper mapper;
mapper.prepare(med_source_mesh,med_target_mesh,"P0P1");
mapper.transfer(sourceField,targetField);
//use targetField
...

middle-level usage

This mode is the mode that needs the minimum of prerequisites (algorithms and the datastructure you intend to use). On the other hand it is needed to specify precisely nature of interpolator.

As consequence of genericity of interpolators, they are usable only by instanciating an underneath mesh data structure. The two following examples show how to use interpolator at this level.

  • The simplest way to use interpolator with MEDCoupling datastruture is put in the following example. Note that this code is close to those used by ParaMEDMEM to perform synchronization of meshes between processes of a MPI communicator :
...
MEDCouplingUMesh *med_source_mesh=MEDLoader::ReadUMeshFromFile("source.med","Source_mesh",0);
MEDCouplingUMesh *med_target_mesh=MEDLoader::ReadUMeshFromFile("target.med","Target_mesh",0);
MEDCouplingNormalizedUnstructuredMesh<2,2> wrap_source_mesh(med_source_mesh);
MEDCouplingNormalizedUnstructuredMesh<2,2> wrap_target_mesh(med_target_mesh);
// Go for interpolation...
INTERP_KERNEL::Interpolation2D myInterpolator;
myInterpolator.setPrecision(1e-7);
myInterpolator.setIntersectionType(INTERP_KERNEL::Geometric2D);
std::vector<std::map<int,double> > resultMatrix;
INTERP_KERNEL::Matrix<double,ALL_C_MODE> resultMatrix2;
// here the interpolation is performed twice for this code to show the capability of storing data of out matrix in 2 different data structures.
myInterpolator.interpolateMeshes(wrap_source_mesh,wrap_target_mesh,resultMatrix,"P0P0");
myInterpolator.interpolateMeshes(wrap_source_mesh,wrap_target_mesh,resultMatrix2,"P0P0");
//Ok resultMatrix and resultMatrix2 contain matrix now
...
  • An another way to use the interpolator with MEDMEM datastructure is :
...
std::string sourceFileName("source.med");
MEDMEM::MESH med_source_mesh(MED_DRIVER,sourceFileName,"Source_Mesh");
std::string targetFileName("target.med");
MEDMEM::MESH med_target_mesh(MED_DRIVER,targetFileName,"Target_Mesh");
// Ok at this point we have our mesh in MED-Memory format.
// Go to wrap med_source_mesh and med_target_mesh.
MEDNormalizedUnstructuredMesh<2,2> wrap_source_mesh(&med_source_mesh);
MEDNormalizedUnstructuredMesh<2,2> wrap_target_mesh(&med_target_mesh);
// Go for interpolation...
INTERP_KERNEL::Interpolation2D myInterpolator;
//optionnal call to parametrize your interpolation. First precision, tracelevel, intersector wanted.
myInterpolator.setOptions(1e-7,0,Geometric2D);
INTERP_KERNEL::Matrix<double,ALL_FORTRAN_MODE> resultMatrix;
myInterpolator.interpolateMeshes(wrap_source_mesh,wrap_target_mesh,resultMatrix,"P0P0");
//Ok let's multiply resultMatrix by source field to interpolate to target field.
resultMatrix.multiply(...)
...
  • Same with VTK datastructure :
...
vtkXMLUnstructuredGridReader *readerSource=vtkXMLUnstructuredGridReader::New();
readerSource->SetFileName("source.vtu");
vtkUnstructuredGrid *vtk_source_mesh=readerSource->GetOutput();
readerSource->Update();
vtkXMLUnstructuredGridReader *readerTarget=vtkXMLUnstructuredGridReader::New();
readerTarget->SetFileName("target.vtu");
vtkUnstructuredGrid *vtk_target_mesh=readerTarget->GetOutput();
readerTarget->Update();
// Ok at this point we have our mesh in VTK format.
// Go to wrap vtk_source_mesh and vtk_target_mesh.
VTKNormalizedUnstructuredMesh<2> wrap_source_mesh(vtk_source_mesh);
VTKNormalizedUnstructuredMesh<2> wrap_target_mesh(vtk_target_mesh);
// Go for interpolation...
INTERP_KERNEL::Interpolation2D myInterpolator;
//optionnal call to parametrize your interpolation. First precision, tracelevel, intersector wanted.
myInterpolator.setOptions(1e-7,0,Geometric2D);
INTERP_KERNEL::Matrix<double,ALL_C_MODE> resultMatrix;
myInterpolator.interpolateMeshes(wrap_source_mesh,wrap_target_mesh,resultMatrix,"P0P0");
//Ok let's multiply resultMatrix by source field to interpolate to target field.
resultMatrix.multiply(...)
//clean-up
readerSource->Delete();
readerTarget->Delete();
...

Theory of conservative interpolation

At the basis of many CFD numerical schemes is the fact that physical quantities such as density, momentum per unit volume or energy per unit volume obey some balance laws that should be preserved at the discrete level on every cell. This property is critical for example to accurately capture shockwaves.

Mesh overlapping

When interpolation is performed between a source mesh S and a target mesh T the aspect of overlapping is important. In fact if any cell of of S is fully overlapped by cells of T and inversely any cell of T is fully overlapped by cells of S the meshes S and T are said to be coincident and some general formulae in next sub section are simpler. As far as possible in the next sub sections the formulae for coincident and non-coincident meshes will be given.

Linear conservative remapping

For fields with polynomial representation on each cell, the components of the discretized field $ \phi_s $ on the source side can be expressed as linear combinations of the components of the discretized field $ \phi_t $ on the target side, in terms of a matrix-vector product :

\[ \phi_t=W.\phi_s. \]

The objectives of interpolators is to compute the matrix W depending on their physical properties (intensive or extensive field) and their mesh discretisation (P0, P1,...).

It is often desired that the process interpolation preserve the integral of $ \phi $ on any domain. At the discrete level, for any target cell $ T_i $, the following general interpolation equation has to be always verified :

\[ \int_{T_i} \phi = \sum_{S_j\cap T_i \neq \emptyset} \int_{T_i\cap S_j} \phi. \]

This equation is used to compute $ W_{ij} $, based on the fields representation ( P0, P1, P1d etc..) and the geometry of source and target mesh cells. :

Conservative remapping of P0 (cell based) fields

We assume that the field is represented by a vector with a discrete value on each cell. This value can represent either

  • 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,
  • an integrated value over the cell (total mass, power of the cell) in which case the representation is said to be extensive

cell-cell (P0->P0) conservative remapping of intensive fields

In the general interpolation equation the left hand side becomes :

\[ \int_{T_i} \phi = (\sum_{S_j} Vol(T_i\cap S_j)).\phi_{T_i}. \]

Here Vol represents the volume when the mesh dimension is equal to 3, the area when mesh dimension is equal to 2, and length when mesh dimension is equal to 1.

Note that $ \sum_{S_j} Vol(T_i\cap S_j) = Vol(T_i) $ in case of perfect overlapping.

In the general interpolation equation the right hand side becomes :

\[ \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}. \]

As the field values are constant on each cell, the coefficients of the linear remapping matrix $ W $ are given by the formula :

\[ W_{ij}=\frac{Vol(T_i\cap S_j)}{ \sum_{S_j} Vol(T_i\cap S_j) }. \]

and in case of perfect overlapping :

\[ W_{ij}=\frac{Vol(T_i\cap S_j)}{ Vol(T_i) }. \]

cell-cell (P0->P0) conservative remapping of extensive physical quantities

In code coupling from neutronics to hydraulics, extensive field of power is exchanged and the total power should remain the same. The discrete values of the field represent the total power contained in the cell. Hence in the general interpolation equation the left hand side becomes :

\[ \int_{T_i} \phi = P_{T_i}, \]

while the right hand side is now :

\[ \sum_{S_j\cap T_i \neq \emptyset} \int_{T_i\cap S_j} \phi = \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}. \]

Note $ \sum_{T_i} Vol(T_i \cap S_j) = Vol(S_j) $ in case of perfect overlapping.

The coefficients of the linear remapping matrix $ W $ are then given by the formula :

\[ W_{ij}=\frac{Vol(T_i\cap S_j)}{ \sum_{T_i} Vol(T_i \cap S_j) }. \]

and in case of perfect overlapping :

\[ W_{ij}=\frac{Vol(T_i\cap S_j)}{ Vol(S_j) }. \]

Copyright © 2007-2011 CEA/DEN, EDF R&D, OPEN CASCADE
Copyright © 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS