Interpolation

About interpolation

PLaSK provides a mechanism to calculate (interpolate) a field of some physical quantity in arbitrary requested points, if values of this field are known in different points. Both sets of points are described by meshes and associated with the vectors of corresponding values.

Let us denote:

• src_mesh - set of points in which the field values are known
• src_vec - vector of known field values, in points described by the sec_mesh
• dst_mesh - requested set of points, in which field values should be calculated (interpolated)
• dst_vec - vector of field values in points described by the dst_mesh, to calculate

plask::interpolate method calculates and returns dst_vec for a given src_mesh, src_vec, dst_mesh, interpolation method, and interpolation flags that specify behavior outside of the source mesh.

plask::interpolate can return a newly created vector or the src_vec if src_mesh and dst_mesh are the same. Furthermore, the lifespan of both source and destination data cannot be determined in advance. For this reason src_vec is passed and dst_vec is returned through a DataVector class, which is responsible for deleting the data in the proper time (i.e. when all the existing solvers delete their copies of the pointer, indicating they are not going to use this data any more). However, for this mechanism to work efficiently, all the solvers should allocate the data using DataVector, as described in Solvers.

Typically, plask::interpolate is called inside providers of the fields of scalars or vectors (see Providers and receivers). Note that the exact src_mesh type must be known at the compile time by the solver providing the data (usually this is the native mesh of the solver). Interpolation algorithms depend on this type. On the other hand the dst_mesh can be very generic and needs only to provide the iterator over its points. This mechanism allows to connect the solver providing some particular data with any solver requesting it, regardless of its mesh.

How to write a new interpolation algorithm?

To implement an interpolation method (which you must usually do in any case where your mesh is the source mesh) you have to write a specialization or a partial specialization of the plask::InterpolationAlgorithm class template for specific: source mesh type, data type, and/or interpolation method. Your specialization must contain an implementation of the static method plask::InterpolationAlgorithm::interpolate.

For example to implement linear interpolation for MyMeshType source mesh type using the same code for all possible data types, you write:

template <typename SrcT, typename DstT>    // for any data type
struct plask::InterpolationAlgorithm<MyMeshType, SrcT, DstT, plask::INTERPOLATION_LINEAR> {
    static void interpolate(MyMeshType& src_mesh, const DataVector<const SrcT>& src_vec, const plask::MeshD<MyMeshType::DIM>& dst_mesh, DataVector<DstT>& dst_vec) {
 
        // here comes your interpolation code
    }
};

Note that above code is template and must be placed in a header file.

The next example, shows how to implement algorithm for a particular data type. To implement the interpolation version for the 'double' type, you should write:

template <>
struct plask::InterpolationAlgorithm<MyMeshType, double, plask::INTERPOLATION_LINEAR> {
    static void interpolate(MyMeshType& src_mesh, const DataVector<const double>& src_vec, const plask::MeshD<MyMeshType::DIM>& dst_mesh, DataVector<double>& dst_vec) {
 
        // interpolation code for vectors of doubles
    }
};

You can simultaneously have codes from both examples. In such case, when linear interpolation from the source mesh MyMeshType is requested, the compiler will use the second implementation to interpolate vectors of doubles and the first one in all other cases.

Typically, the code of the function should iterate over all the points of the dst_mesh and fill the dst_vec with the interpolated values in the respective points.