evaluator_2d.hpp

#pragma once
#include <ddc/ddc.hpp>
 
#include "ddc_aliases.hpp"
 
struct Evaluator2D
{
    template <class Eval1, class Eval2>
    class Evaluator
    {
    private:
        Eval1 eval_func1;
        Eval2 eval_func2;
 
    public:
        template <class IdxRange>
        Evaluator(IdxRange idx_range)
            : eval_func1(ddc::select<typename Eval1::Dim>(idx_range))
            , eval_func2(ddc::select<typename Eval2::Dim>(idx_range))
        {
        }
 
        double operator()(double const x, double const y) const noexcept
        {
            return eval_func1(x) * eval_func2(y);
        }
 
        template <class Grid1, class Grid2>
        double operator()(ddc::Coordinate<Grid1, Grid2> const x) const noexcept
        {
            return eval_func1(ddc::get<Grid1>(x)) * eval_func2(ddc::get<Grid2>(x));
        }
 
        template <class Grid1, class Grid2>
        void operator()(ddc::ChunkSpan<double, ddc::DiscreteDomain<Grid1, Grid2>> chunk) const
        {
            auto const& idx_range = chunk.idx_range();
 
            ddc::for_each(idx_range, [=](ddc::DiscreteElement<Grid1, Grid2> const i) {
                chunk(i) = eval_func1(ddc::coordinate(ddc::select<Grid1>(i)))
                           * eval_func2(ddc::coordinate(ddc::select<Grid2>(i)));
            });
        }
 
        double deriv(double const x, double const y, int const derivative_x, int const derivative_y)
                const noexcept
        {
            return eval_func1.deriv(x, derivative_x) * eval_func2.deriv(y, derivative_y);
        }
 
        template <class Grid1, class Grid2>
        double deriv(
                ddc::Coordinate<Grid1, Grid2> const x,
                int const derivative_x,
                int const derivative_y) const noexcept
        {
            return eval_func1.deriv(ddc::get<Grid1>(x), derivative_x)
                   * eval_func2.deriv(ddc::get<Grid2>(x), derivative_y);
        }
 
        template <class Grid1, class Grid2>
        void deriv(
                ddc::ChunkSpan<double, ddc::DiscreteDomain<Grid1, Grid2>> chunk,
                int const derivative_x,
                int const derivative_y) const
        {
            auto const& idx_range = chunk.idx_range();
 
            ddc::for_each(idx_range, [=](ddc::DiscreteElement<Grid1, Grid2> const i) {
                chunk(i) = eval_func1.deriv(ddc::coordinate(ddc::select<Grid1>(i)), derivative_x)
                           * eval_func2.deriv(ddc::coordinate(ddc::select<Grid2>(i)), derivative_y);
            });
        }
 
        double max_norm(int diff1 = 0, int diff2 = 0) const
        {
            return eval_func1.max_norm(diff1) * eval_func2.max_norm(diff2);
        }
    };
};