spline_interpolator_2d_rp.hpp

// SPDX-License-Identifier: MIT
#pragma once
#include "ddc_alias_inline_functions.hpp"
#include "ddc_aliases.hpp"
#include "geometry.hpp"
#include "i_interpolator_2d_rp.hpp"
 
 
template <class RadialExtrapolationRule>
class SplineInterpolatorRTheta : public IInterpolatorRTheta
{
public:
    using evaluator_type = ddc::SplineEvaluator2D<
            Kokkos::DefaultHostExecutionSpace,
            Kokkos::HostSpace,
            BSplinesR,
            BSplinesTheta,
            GridR,
            GridTheta,
            RadialExtrapolationRule,
            RadialExtrapolationRule,
            ddc::PeriodicExtrapolationRule<Theta>,
            ddc::PeriodicExtrapolationRule<Theta>,
            GridR,
            GridTheta>;
 
private:
    SplineRThetaBuilder const& m_builder;
 
    evaluator_type const& m_evaluator;
 
    mutable host_t<DFieldMem<IdxRangeBSRTheta>> m_coefs;
 
    using r_deriv_type = host_t<DConstField<SplineRThetaBuilder::batched_derivs_domain_type1>>;
    using p_deriv_type = host_t<DConstField<SplineRThetaBuilder::batched_derivs_domain_type2>>;
    using mixed_deriv_type = host_t<DConstField<SplineRThetaBuilder::batched_derivs_domain_type>>;
 
public:
    SplineInterpolatorRTheta(SplineRThetaBuilder const& builder, evaluator_type const& evaluator)
        : m_builder(builder)
        , m_evaluator(evaluator)
        , m_coefs(get_spline_idx_range(builder))
    {
    }
 
    ~SplineInterpolatorRTheta() override = default;
 
    host_t<DFieldRTheta> operator()(
            host_t<DFieldRTheta> const inout_data,
            host_t<Field<CoordRTheta const, IdxRangeRTheta>> const coordinates) const override
    {
#ifndef NDEBUG
        // To ensure that the interpolator is C0, we ensure that
        // the value at (r=0,theta) is the same for all theta.
        IdxRangeR r_idx_range = get_idx_range<GridR>(inout_data);
        IdxRangeTheta theta_idx_range = get_idx_range<GridTheta>(inout_data);
        if (ddc::coordinate(r_idx_range.front()) == 0) {
            ddc::for_each(theta_idx_range, [&](IdxTheta const ip) {
                bool const unicity_center_point
                        = inout_data(r_idx_range.front(), ip)
                          == inout_data(r_idx_range.front(), theta_idx_range.front());
                if (!unicity_center_point) {
                    std::printf("Unicity of the value at the center point is not verified.");
                    assert(unicity_center_point);
                }
            });
        }
#endif
 
        m_builder(get_field(m_coefs), get_const_field(inout_data));
        m_evaluator(get_field(inout_data), coordinates, get_const_field(m_coefs));
        return inout_data;
    }
};
 
 
 
template <class RadialExtrapolationRule>
class PreallocatableSplineInterpolatorRTheta : public IPreallocatableInterpolatorRTheta
{
public:
    using evaluator_type = ddc::SplineEvaluator2D<
            Kokkos::DefaultHostExecutionSpace,
            Kokkos::HostSpace,
            BSplinesR,
            BSplinesTheta,
            GridR,
            GridTheta,
            RadialExtrapolationRule,
            RadialExtrapolationRule,
            ddc::PeriodicExtrapolationRule<Theta>,
            ddc::PeriodicExtrapolationRule<Theta>,
            GridR,
            GridTheta>;
 
private:
    SplineRThetaBuilder const& m_builder;
 
    evaluator_type const& m_evaluator;
 
public:
    PreallocatableSplineInterpolatorRTheta(
            SplineRThetaBuilder const& builder,
            evaluator_type const& evaluator)
        : m_builder(builder)
        , m_evaluator(evaluator)
    {
    }
 
    ~PreallocatableSplineInterpolatorRTheta() override = default;
 
    std::unique_ptr<IInterpolatorRTheta> preallocate() const override
    {
        return std::make_unique<
                SplineInterpolatorRTheta<RadialExtrapolationRule>>(m_builder, m_evaluator);
    }
};