grid_interpolation.hh

#ifndef VTESTBED_CORE_GRID_INTERPOLATION_HH
#define VTESTBED_CORE_GRID_INTERPOLATION_HH

#include <vtestbed/core/grid.hh>

namespace vtb {

    namespace core {

        template <class T, int N, class Function, class Return=T>
        class Grid_interpolation {

        public:
            using vertex_type = blitz::TinyVector<T,N>;
            using grid_type = Grid<T,N>;
            using index_type = blitz::TinyVector<int,N>;
            using function_type = Function;
            using return_type = Return;

        private:
            grid_type _grid;
            function_type _function;

        public:

            inline explicit
            Grid_interpolation(const grid_type& grid, function_type function):
            _grid(grid), _function(function) {}

            inline return_type
            operator()(const vertex_type& x) {
                #if defined(VTB_DEBUG_INTERPOLATION)
                this->detect_extrapolation(x);
                #endif
                const auto& grid = this->_grid;
                index_type idx_max{grid.shape()-1};
                index_type i00{floor(grid.index(x))};
                i00 = blitz::clamp(i00, grid.begin(), index_type(idx_max-1));
                return this->_function(x, i00);
            }

            inline const grid_type& grid() const { return this->_grid; }
            inline const function_type& function() const { return this->_function; }

        private:

            #if defined(VTB_DEBUG_INTERPOLATION)
            inline void
            detect_extrapolation(const vertex_type& x) {
                const auto& grid = this->_grid;
                if (!grid.contains(x)) {
                    bool success = true;
                    for (int i=0; i<2; ++i) {
                        using std::abs;
                        T x_min = grid.lbound(i);
                        T x_max = grid.ubound(i);
                        if (x(i) < x_min && std::abs(x(i)-x_min) > T{1e-3}) {
                            success = false;
                        }
                        if (x(i) > x_max && std::abs(x(i)-x_max) > T{1e-3}) {
                            success = false;
                        }
                    }
                    if (!success) {
                        std::clog << "_grid=" << _grid << std::endl;
                        std::clog << "x=" << x << std::endl;
                        throw std::runtime_error("extrapolation");
                    }
                }
            }
            #endif

        };

        template <class T>
        class Barycentric_function {

        public:
            using vec3 = blitz::TinyVector<T,3>;
            using vec2 = blitz::TinyVector<T,2>;
            using vertex_matrix = blitz::Array<vec3,2>;
            using index_type = blitz::TinyVector<int,2>;

        private:
            vertex_matrix _points;

        public:

            inline explicit
            Barycentric_function(const vertex_matrix& points): _points(points) {}

            T operator()(const vec2& x, const index_type& i00);

        };

        template <class T,int N>
        class Linear_function {

        public:
            using array_type = blitz::Array<T,N>;
            using vertex_type = blitz::TinyVector<T,N>;
            using index_type = blitz::TinyVector<int,N>;
            using grid_type = Grid<T,N>;

        private:
            array_type _func;
            grid_type _grid;

        public:

            inline explicit
            Linear_function(array_type func, const grid_type& grid):
            _func(func), _grid(grid) {}

            T operator()(const vertex_type& x, const index_type& i00);

        };

        template <class T>
        class Radial_basis_function {

        public:
            using vec3 = blitz::TinyVector<T,3>;
            using vec2 = blitz::TinyVector<T,2>;
            using vertex_matrix = blitz::Array<vec3,2>;
            using index_type = blitz::TinyVector<int,2>;

        private:
            vertex_matrix _points;
            T _order{1};
            T _sigma{1};

        public:
            inline explicit
            Radial_basis_function(
                const vertex_matrix& points,
                T order = T{1},
                T sigma = T{1}
            ): _points(points), _order(order), _sigma(sigma) {}

            T operator()(const vec2& x, const index_type& i00);

        };

        template <class T, int N>
        class Linear_interpolation:
            public Grid_interpolation<T,N,Linear_function<T,N>> {

        public:
            using function_type = Linear_function<T,N>;
            using array_type = blitz::Array<T,N>;
            using grid_type = Grid<T,N>;

        private:
            using base_type = Grid_interpolation<T,N,function_type>;

        public:
            using base_type::base_type;
            using base_type::operator();

            inline explicit
            Linear_interpolation(array_type points, const grid_type& grid):
            base_type(grid, function_type(points,grid)) {}

        };

        template <class T, class Return=T>
        class Linear_function_irregular {

        public:
            using return_type = Return;
            using vec3 = blitz::TinyVector<T,3>;
            using vertex_field = blitz::Array<return_type,3>;
            using index_type = blitz::TinyVector<int,3>;
            using grid_type = Grid<T,3>;

        private:
            vertex_field _points;
            grid_type _grid;

        public:
            inline explicit
            Linear_function_irregular(const vertex_field& points, const grid_type& grid):
            _points(points), _grid(grid) {}

            return_type operator()(const vec3& x, const index_type& i00);

        };

    }

}

#endif // vim:filetype=cpp