Med Memory Users' Guide 5.1.3
|
The main purpose of this module is to propose a set of algorithms for mesh interpolation fully independant of 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 interpolation kernel.
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 coincidant and some general formulae in next sub section are simpler. As far as possible in the next sub sections the formulae for coincidant and non-coincidant meshes will be given.
For fields with polynomial representation on each cell, the components of the discretized field on the source side can be expressed as linear combinations of the components of the discretized field
on the target side, in terms of a matrix-vector product :
All the aim of interpolators is to compute W depending on a physical quantities and the type of interpolation wanted (P0, P1, P1d etc...).
For the intensive field from the source mesh
to the target mesh
the interpolation should preserve the integral of
on any domain. At the discrete level, for any target cell
, the following general interpolation equation has to be always verified :
To compute 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 remap fields with P0 representations.
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.
In the general interpolation equation the left hand side becomes :
In the general interpolation equation the right hand side becomes :
In the case where the intensive field values are constant on each cell, the coefficients of the linear remapping matrix are given by the formula :
and in case of perfect overlapping :
Where Vol represents the volume with mesh dimension of interpolation equals to 3, the area when mesh dimension equals to 2, and length when mesh dimension equals to 1.
In code coupling from neutronics to hydraulic code extensive field of power is exchanged and the all power as to be kept the same. The principle is to 'intensify' the field to move on from extensive field P to an intensive one .
In the general interpolation equation the left hand side becomes :
In the general interpolation equation the right hand side becomes :
In the case where the extensive field values are constant on each cell, the coefficients of the linear remapping matrix are given by the formula :
and in case of perfect overlapping :
In interpolation kernel, the algorithm that computes given cell i in source mesh and cell j in target mesh is called intersector.
As seen in the theory of interpolation, for all interpolation the aim is to fill the W matrix (which is generally a sparse matrix). Fundamatally for each pair (i,j) is obtained by calling each time the wanted intersector. The problem is that each call to algorithm is CPU-expensive. To reduce the computation time a first filtering is done to found a maximim of (i,j) pairs where
is obviously equal to 0. Typically it is the case when a cell in source mesh is too far from an another cell in target mesh each other.
So for a given type of interpolation, the computation of W is performed in two steps :
Each interpolator inherits from INTERP_KERNEL::Interpolation ( whatever its dimension and its type ) that is a CRTP class in order to clearly see the main API without useless CPU cost.
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 several formats without preformance loss. For example MED, VTK...
MeshType
is a concept that should strictly fulfilled the following rules :
int
or long
int
.void getBoundingBox(double *boundingBox) const
INTERP_KERNEL::NormalizedCellType getTypeOfElement(MyConnType eltId) const
unsigned char getNumberOfNodesOfElement(MyConnType eltId) const
unsigned long getNumberOfNodes() const
unsigned long getNumberOfElements() const
const MyConnType *getConnectivityPtr() const
const double *getCoordinatesPtr() const
const MyConnType *getConnectivityIndexPtr() const
void ReleaseTempArrays()
getCoordinatesPtr
must be a full interlace array.getConnectivityPtr
and
As already said, the matrix returned by interpolator is typically a sparse matrix. Instances of class
MatrixType
are used to stores these results of interpolation. To be able to be filled by the interpolator the MatrixType
class has to match the following concept :
void resize(uint nbrows)
Row &operator [] (uint irow)
class
Row
has to match at least following concept :
void insert(const std::pair<int,T>& myPair)
std::vector<
std::map<int,double>
>
is a candidate for MatrixType
.This the simplest mode of usage of interpolator. This way is strongly linked with MED data-structure. All interpolators is completely hidden to you. Even sparse interpolation matrix is hidden. An exemple of usage :
... 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); mapper.transfer(sourceField,targetField); //use targetField ...
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.
... 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(...) ...
... 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(); ...