interpolation.hpp

Go to the documentation of this file.

/*
 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
 * Copyright (c) 2022 Lodz University of Technology
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 3.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 */
#ifndef PLASK__INTERPOLATION_H
#define PLASK__INTERPOLATION_H
 
#include <typeinfo>  // for 'typeid'
 
#include "mesh.hpp"
#include "../exceptions.hpp"
#include "../memory.hpp"
#include "../lazydata.hpp"
#include "../log/log.hpp"
#include "plask/geometry/space.hpp"
 
namespace plask {
 
enum InterpolationMethod: unsigned {
    INTERPOLATION_DEFAULT = 0,  
    INTERPOLATION_NEAREST,      
    INTERPOLATION_LINEAR,       
    INTERPOLATION_SPLINE,       
    INTERPOLATION_SMOOTH_SPLINE,
    INTERPOLATION_PERIODIC_SPLINE,
    INTERPOLATION_FOURIER,      
    // ...add new interpolation algorithms here...
#   ifndef DOXYGEN
    __ILLEGAL_INTERPOLATION_METHOD__  // necessary for metaprogram loop and automatic Python enums
#   endif // DOXYGEN
};
 
PLASK_API extern const char* interpolationMethodNames[__ILLEGAL_INTERPOLATION_METHOD__+1];
 
template <InterpolationMethod default_method>
inline InterpolationMethod getInterpolationMethod(InterpolationMethod method) {
    return (method == INTERPOLATION_DEFAULT)? default_method : method;
}
 
class InterpolationFlags {
    unsigned char sym[3];   
    unsigned char per;      
    double lo[3], hi[3];    
 
public:
 
    enum class Symmetry: unsigned char {
        NO = 0,
        POSITIVE = 1,
        PP = 1,
        NP = 3,
        PN = 5,
        NN = 7,
        PPP = 1,
        NPP = 3,
        PNP = 5,
        NNP = 7,
        PPN = 9,
        NPN = 11,
        PNN = 13,
        NNN = 15,
        NEGATIVE = 15,
    };
 
    InterpolationFlags(): sym{0,0,0}, per(0), lo{0.,0.,0.} , hi{0., 0., 0.}  {}
 
    InterpolationFlags(shared_ptr<const GeometryD<2>> geometry, Symmetry sym0, Symmetry sym1):
        sym{geometry->isSymmetric(Geometry::DIRECTION_TRAN)? (unsigned char)sym0 : (unsigned char)0, geometry->isSymmetric(Geometry::DIRECTION_VERT)? (unsigned char)sym1 : (unsigned char)0},
        per((unsigned char)((geometry->isPeriodic(Geometry::DIRECTION_TRAN)? 1 : 0) + (geometry->isPeriodic(Geometry::DIRECTION_VERT)? 2 : 0))),
        lo{geometry->getChildBoundingBox().left(), geometry->getChildBoundingBox().bottom(), 0.},
        hi{geometry->getChildBoundingBox().right(), geometry->getChildBoundingBox().top(), 0.}
    {
        if (geometry->isSymmetric(Geometry::DIRECTION_TRAN)) {
            if (lo[0] < 0. && hi[0] > 0.)
                throw Exception("interpolation: Symmetric geometry spans at both sides of transverse axis");
            if (!sym[0]) {
                hi[0] = max(-lo[0], hi[0]); lo[0] = -hi[0];
            }
        }
        if (geometry->isSymmetric(Geometry::DIRECTION_VERT)) {
            if (lo[1] < 0. && hi[1] > 0.)
                throw Exception("interpolation: Symmetric geometry spans at both sides of vertical axis");
            if (!sym[1]) {
                hi[1] = max(-lo[1], hi[1]); lo[1] = -hi[1];
            }
        }
    }
 
    InterpolationFlags(shared_ptr<const Geometry2DCartesian> geometry): InterpolationFlags(geometry, Symmetry::POSITIVE, Symmetry::POSITIVE) {}
    InterpolationFlags(shared_ptr<Geometry2DCartesian> geometry): InterpolationFlags(geometry, Symmetry::POSITIVE, Symmetry::POSITIVE) {}
    InterpolationFlags(shared_ptr<const Geometry2DCylindrical> geometry): InterpolationFlags(geometry, Symmetry::POSITIVE, Symmetry::POSITIVE) {}
    InterpolationFlags(shared_ptr<Geometry2DCylindrical> geometry): InterpolationFlags(geometry, Symmetry::POSITIVE, Symmetry::POSITIVE) {}
 
    InterpolationFlags(shared_ptr<const GeometryD<3>> geometry, Symmetry sym0, Symmetry sym1, Symmetry sym2):
        sym{geometry->isSymmetric(Geometry::DIRECTION_LONG)? (unsigned char)sym0 : (unsigned char)0, geometry->isSymmetric(Geometry::DIRECTION_TRAN)? (unsigned char)sym1 : (unsigned char)0, geometry->isSymmetric(Geometry::DIRECTION_VERT)? (unsigned char)sym2 : (unsigned char)0},
        per((unsigned char)((geometry->isPeriodic(Geometry::DIRECTION_LONG)? 1 : 0) + (geometry->isPeriodic(Geometry::DIRECTION_TRAN)? 2 : 0) + (geometry->isPeriodic(Geometry::DIRECTION_VERT)? 4 : 0))),
        lo{geometry->getChildBoundingBox().back(), geometry->getChildBoundingBox().left(), geometry->getChildBoundingBox().bottom()},
        hi{geometry->getChildBoundingBox().front(), geometry->getChildBoundingBox().right(), geometry->getChildBoundingBox().top()}
    {
        if (geometry->isSymmetric(Geometry::DIRECTION_LONG)) {
            if (lo[0] < 0. && hi[0] > 0.)
                throw Exception("interpolation: Symmetric geometry spans at both sides of longitudinal axis");
            if (!sym[0]) {
                hi[0] = max(-lo[0], hi[0]); lo[0] = -hi[0];
            }
        }
        if (geometry->isSymmetric(Geometry::DIRECTION_TRAN)) {
            if (lo[1] < 0. && hi[1] > 0.)
                throw Exception("interpolation: Symmetric geometry spans at both sides of transverse axis");
            if (!sym[1]) {
                hi[1] = max(-lo[1], hi[1]); lo[1] = -hi[1];
            }
        }
        if (geometry->isSymmetric(Geometry::DIRECTION_VERT)) {
            if (lo[2] < 0. && hi[2] > 0.)
                throw Exception("interpolation: Symmetric geometry spans at both sides of vertical axis");
            if (!sym[2]) {
                hi[2] = max(-lo[2], hi[2]); lo[2] = -hi[2];
            }
        }
    }
 
    InterpolationFlags(shared_ptr<const Geometry3D> geometry): InterpolationFlags(geometry, Symmetry::POSITIVE, Symmetry::POSITIVE, Symmetry::POSITIVE) {}
    InterpolationFlags(shared_ptr<Geometry3D> geometry): InterpolationFlags(geometry, Symmetry::POSITIVE, Symmetry::POSITIVE, Symmetry::POSITIVE) {}
 
    unsigned char symmetric(int axis) const { return sym[axis]; }
 
    bool periodic(int axis) const { return (per & (1 << axis)) != 0; }
 
    double low(int axis) const {
        if (sym[axis]) return min(lo[axis], -hi[axis]);
        else return lo[axis];
    }
 
    double high(int axis) const {
        if (sym[axis]) return max(lo[axis], hi[axis]);
        else return hi[axis];
    }
 
    double wrap(int ax, double x) const {
        double d = hi[ax] - lo[ax];
        if (periodic(ax)) {
            if (sym[ax]) {
                x = std::fmod(abs(x), 2.*d);
                if (x == 0.) x = 1e-12;
                if (x > d) x = - (x - 2.*d);
                if (hi[ax] < 0) x = - x;
            } else {
                x = std::fmod(x-lo[ax], d);
                x += (x >= 0)? lo[ax] : hi[ax];
            }
        } else if (sym[ax]) {
            if (x == 0.) x = 1e-12;
            if (lo[ax] >= 0) x = abs(x);
            else x = - abs(x);
        }
        return x;
    }
 
    template <int dim>
    Vec<dim> wrap(Vec<dim> pos) const {
        for (int i = 0; i != dim; ++i) {
            pos[i] = wrap(i, pos[i]);
        }
        return pos;
    }
 
    template <int dim, typename DataT>
    Vec<dim,DataT> reflect(int ax, Vec<dim,DataT> vec) const {
        for (int i = 0; i != dim; ++i) if (sym[ax] & (2 << i)) vec[i] = -vec[i];
        return vec;
    }
 
    template <typename DataT>
    DataT reflect(int ax, DataT val) const {
        if (sym[ax] & 14) return -val;
        else return val;
    }
 
    template <int dim, typename DataT>
    DataT postprocess(Vec<dim> pos, DataT data) const {
        for (int i = 0; i != dim; ++i) {
            if (sym[i]) {
                if (periodic(i)) {
                    double d = hi[i] - lo[i];
                    pos[i] = std::fmod(pos[i], 2.*d);
                    if (pos[i] > d || (pos[i] < 0. && pos[i] > -d)) data = reflect(i, data);
                } else {
                    if (lo[i] >= 0.) { if (pos[i] < 0.) data = reflect(i, data); }
                    else { if (pos[i] > 0.) data = reflect(i, data); }
                }
            }
        }
        return data;
    }
};
 
template <typename SrcMeshT, typename SrcT, typename DstT, InterpolationMethod method>
struct InterpolationAlgorithm
{
    static LazyData<DstT> interpolate(const shared_ptr<const SrcMeshT>& src_mesh, const DataVector<const SrcT>&,
                                      const shared_ptr<const MeshD<SrcMeshT::DIM>>&, const InterpolationFlags&) {
        std::string msg = "interpolate (source mesh type: ";
        msg += typeid(*src_mesh).name();
        msg += ", interpolation method: ";
        msg += interpolationMethodNames[method];
        msg += ")";
        throw NotImplemented(msg);
    }
};
 
template <typename SrcMeshT, typename SrcT, typename DstT>
struct InterpolationAlgorithm<SrcMeshT, SrcT, DstT, INTERPOLATION_DEFAULT>
{
    static LazyData<DstT> interpolate(const shared_ptr<const SrcMeshT>&, const DataVector<const SrcT>&,
                                      const shared_ptr<const MeshD<SrcMeshT::DIM>>&, const InterpolationFlags& PLASK_UNUSED(flags)) {
        throw CriticalException("interpolate(...) called for INTERPOLATION_DEFAULT method. Contact solver author to fix this issue."
#ifndef NDEBUG
                                "\n\nINFO FOR SOLVER AUTHOR: To avoid this error use 'getInterpolationMethod<YOUR_DEFAULT_METHOD>(interpolation_method) in C++ code of the provider in your solver.\n"
#endif
                               );
    }
};
 
#ifndef DOXYGEN
// The following structures are solely used for metaprogramming
template <typename SrcMeshT, typename SrcT, typename DstT, int iter>
struct __InterpolateMeta__
{
    inline static LazyData<typename std::remove_const<DstT>::type> interpolate(
            const shared_ptr<const SrcMeshT>& src_mesh, const DataVector<const SrcT>& src_vec,
            const shared_ptr<const MeshD<SrcMeshT::DIM>>& dst_mesh, InterpolationMethod method, const InterpolationFlags& flags) {
        if (int(method) == iter)
            return InterpolationAlgorithm<SrcMeshT, SrcT, typename std::remove_const<DstT>::type, (InterpolationMethod)iter>
                    ::interpolate(src_mesh, DataVector<const SrcT>(src_vec), dst_mesh, flags);
        else
            return __InterpolateMeta__<SrcMeshT, SrcT, DstT, iter+1>::interpolate(src_mesh, src_vec, dst_mesh, method, flags);
    }
};
template <typename SrcMeshT, typename SrcT, typename DstT>
struct __InterpolateMeta__<SrcMeshT, SrcT, DstT, __ILLEGAL_INTERPOLATION_METHOD__>
{
    inline static LazyData<typename std::remove_const<DstT>::type> interpolate(const shared_ptr<const SrcMeshT>&, const DataVector<const SrcT>&,
                                   const shared_ptr<const MeshD<SrcMeshT::DIM>>&, InterpolationMethod, const InterpolationFlags&) {
        throw CriticalException("no such interpolation method");
    }
};
#endif // DOXYGEN
 
 
template <typename SrcMeshT, typename SrcT, typename DstT=SrcT>
LazyData<typename std::remove_const<DstT>::type> interpolate(
                                        shared_ptr<const SrcMeshT> src_mesh,
                                        DataVector<const SrcT> src_vec,
                                        shared_ptr<const MeshD<SrcMeshT::DIM>> dst_mesh,
                                        InterpolationMethod method=INTERPOLATION_DEFAULT,
                                        const InterpolationFlags& flags=InterpolationFlags(),
                                        bool verbose=true)
{
    if (src_mesh->size() != src_vec.size())
        throw BadMesh("interpolate", "mesh size ({1}) and values size ({0}) do not match", src_vec.size(), src_mesh->size());
    if (src_mesh == dst_mesh) return new LazyDataFromVectorImpl<typename std::remove_const<DstT>::type>(src_vec); // meshes are identical, so just return src_vec
    if (verbose && method < __ILLEGAL_INTERPOLATION_METHOD__) writelog(LOG_DEBUG, "interpolate: Running {0} interpolation", interpolationMethodNames[method]);
    return __InterpolateMeta__<SrcMeshT, SrcT, DstT, 0>::interpolate(src_mesh, src_vec, dst_mesh, method, flags);
}
 
template <typename SrcMeshT, typename SrcT, typename DstT=SrcT, typename DstMeshT>
LazyData<typename std::remove_const<DstT>::type> interpolate(
                            shared_ptr<SrcMeshT> src_mesh,
                            DataVector<SrcT> src_vec,
                            shared_ptr<DstMeshT> dst_mesh,
                            InterpolationMethod method=INTERPOLATION_DEFAULT,
                            const InterpolationFlags& flags=InterpolationFlags(),
                            bool verbose=true)
{
    return interpolate(shared_ptr<const SrcMeshT>(src_mesh), DataVector<const SrcT>(src_vec),
                       shared_ptr<const MeshD<SrcMeshT::DIM>>(dst_mesh), method, flags, verbose);
}
 
 
/*template <typename SrcMeshT, typename SrcT, typename DstT=SrcT>
LazyData<typename std::remove_const<DstT>::type> interpolate(shared_ptr<SrcMeshT> src_mesh, DataVector<const SrcT> src_vec,
                             shared_ptr<MeshD<SrcMeshT::DIM>> dst_mesh,
                             InterpolationMethod method=INTERPOLATION_DEFAULT, bool verbose=true)
{
    return interpolate(shared_ptr<const SrcMeshT>(src_mesh), src_vec, shared_ptr<const MeshD<SrcMeshT::DIM>>(dst_mesh), method, verbose);
}
 
template <typename SrcMeshT, typename SrcT, typename DstT=SrcT>
LazyData<typename std::remove_const<DstT>::type> interpolate(shared_ptr<const SrcMeshT> src_mesh, DataVector<const SrcT> src_vec,
                             shared_ptr<MeshD<SrcMeshT::DIM>> dst_mesh,
                             InterpolationMethod method=INTERPOLATION_DEFAULT, bool verbose=true)
{
    return interpolate(src_mesh, src_vec, shared_ptr<const MeshD<SrcMeshT::DIM>>(dst_mesh), method, verbose);
}*/
 
template <typename DstT, typename SrcMeshType, typename SrcT = DstT>
struct InterpolatedLazyDataImpl: public LazyDataImpl<DstT> {
 
    shared_ptr<const SrcMeshType> src_mesh;
    shared_ptr<const MeshD<SrcMeshType::DIM>> dst_mesh;
    DataVector<const SrcT> src_vec;
    InterpolationFlags flags;
 
    InterpolatedLazyDataImpl(const shared_ptr<const SrcMeshType>& src_mesh, const DataVector<const SrcT>& src_vec,
                             const shared_ptr<const MeshD<SrcMeshType::DIM>>& dst_mesh, const InterpolationFlags& flags)
        : src_mesh(src_mesh), dst_mesh(dst_mesh), src_vec(src_vec), flags(flags) {}
 
    std::size_t size() const override { return dst_mesh->size(); }
 
};
 
template <typename DstT, typename SrcMeshType, typename SrcT = DstT>
struct LinearInterpolatedLazyDataImpl: public InterpolatedLazyDataImpl<DstT, SrcMeshType, SrcT> {
 
    LinearInterpolatedLazyDataImpl(shared_ptr<const SrcMeshType> src_mesh, const DataVector<const SrcT>& src_vec,
                                   shared_ptr<const MeshD<SrcMeshType::DIM>> dst_mesh, const InterpolationFlags& flags):
        InterpolatedLazyDataImpl<DstT, SrcMeshType>(src_mesh, src_vec, dst_mesh, flags) {}
 
    DstT at(std::size_t index) const override {
        return this->src_mesh->interpolateLinear(this->src_vec, this->dst_mesh->at(index), this->flags);
    }
 
};
 
template <typename DstT, typename SrcMeshType, typename SrcT=DstT>
struct NearestNeighborInterpolatedLazyDataImpl: public InterpolatedLazyDataImpl<DstT, SrcMeshType, SrcT> {
 
    NearestNeighborInterpolatedLazyDataImpl(shared_ptr<const SrcMeshType> src_mesh, const DataVector<const SrcT>& src_vec,
                                            shared_ptr<const MeshD<SrcMeshType::DIM>> dst_mesh, const InterpolationFlags& flags):
        InterpolatedLazyDataImpl<DstT, SrcMeshType, DstT>(src_mesh, src_vec, dst_mesh, flags) {}
 
    SrcT at(std::size_t index) const override {
        return this->src_mesh->interpolateNearestNeighbor(this->src_vec, this->dst_mesh->at(index), this->flags);
    }
 
};
 
 
 
} // namespace plask
 
 
namespace boost {
 
    template <>
    inline plask::InterpolationMethod lexical_cast<plask::InterpolationMethod, std::string>(const std::string& arg) {
        std::string val = arg; to_upper(val);
        for (unsigned i = plask::INTERPOLATION_DEFAULT+1; i != plask::__ILLEGAL_INTERPOLATION_METHOD__; ++i) {
            if (val == plask::interpolationMethodNames[i]) return (plask::InterpolationMethod)i;
        }
        throw bad_lexical_cast(typeid(std::string), typeid(plask::InterpolationMethod));
    }
 
}
 
 
 
#endif  //PLASK__INTERPOLATION_H