discrete_to_cartesian.hpp

// SPDX-License-Identifier: MIT
#pragma once
#include <iostream>
 
#include <ddc/ddc.hpp>
 
#include <sll/math_tools.hpp>
#include <sll/view.hpp>
 
#include "mapping_tools.hpp"
 
template <
        class X,
        class Y,
        class SplineEvaluator,
        class R = typename SplineEvaluator::continuous_dimension_type1,
        class Theta = typename SplineEvaluator::continuous_dimension_type2,
        class MemorySpace = typename SplineEvaluator::memory_space>
class DiscreteToCartesian
{
    static_assert(std::is_same_v<MemorySpace, typename SplineEvaluator::memory_space>);
 
public:
    using BSplineR = typename SplineEvaluator::bsplines_type1;
    using BSplineTheta = typename SplineEvaluator::bsplines_type2;
 
    using cartesian_tag_x = X;
    using cartesian_tag_y = Y;
    using curvilinear_tag_r = R;
    using curvilinear_tag_theta = Theta;
 
    using CoordArg = ddc::Coordinate<R, Theta>;
    using CoordResult = ddc::Coordinate<X, Y>;
 
private:
    using spline_idx_range = ddc::DiscreteDomain<BSplineR, BSplineTheta>;
 
    using SplineType = ddc::ChunkView<double, spline_idx_range, Kokkos::layout_right, MemorySpace>;
 
    using IdxRangeRTheta = typename SplineEvaluator::evaluation_domain_type;
    using IdxRangeTheta = typename SplineEvaluator::evaluation_domain_type2;
    using IdxTheta = typename IdxRangeTheta::discrete_element_type;
 
private:
    SplineType m_x_spline_representation;
    SplineType m_y_spline_representation;
    SplineEvaluator m_spline_evaluator;
    IdxRangeRTheta m_idx_range_singular_point;
 
public:
    KOKKOS_FUNCTION DiscreteToCartesian(
            SplineType curvilinear_to_x,
            SplineType curvilinear_to_y,
            SplineEvaluator const& evaluator,
            IdxRangeRTheta idx_range_singular_point)
        : m_x_spline_representation(curvilinear_to_x)
        , m_y_spline_representation(curvilinear_to_y)
        , m_spline_evaluator(evaluator)
        , m_idx_range_singular_point(idx_range_singular_point)
    {
    }
 
    KOKKOS_FUNCTION ddc::Coordinate<X, Y> operator()(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const
    {
        const double x = m_spline_evaluator(coord, m_x_spline_representation.span_cview());
        const double y = m_spline_evaluator(coord, m_y_spline_representation.span_cview());
        return ddc::Coordinate<X, Y>(x, y);
    }
 
    KOKKOS_FUNCTION void jacobian_matrix(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord,
            Matrix_2x2& matrix) const
    {
        matrix[0][0]
                = m_spline_evaluator.deriv_dim_1(coord, m_x_spline_representation.span_cview());
        matrix[0][1]
                = m_spline_evaluator.deriv_dim_2(coord, m_x_spline_representation.span_cview());
        matrix[1][0]
                = m_spline_evaluator.deriv_dim_1(coord, m_y_spline_representation.span_cview());
        matrix[1][1]
                = m_spline_evaluator.deriv_dim_2(coord, m_y_spline_representation.span_cview());
    }
 
    KOKKOS_FUNCTION double jacobian_11(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const
    {
        return m_spline_evaluator.deriv_dim_1(coord, m_x_spline_representation.span_cview());
    }
 
    KOKKOS_FUNCTION double jacobian_12(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const
    {
        return m_spline_evaluator.deriv_dim_2(coord, m_x_spline_representation.span_cview());
    }
 
    KOKKOS_FUNCTION double jacobian_21(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const
    {
        return m_spline_evaluator.deriv_dim_1(coord, m_y_spline_representation.span_cview());
    }
 
    KOKKOS_FUNCTION double jacobian_22(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const
    {
        return m_spline_evaluator.deriv_dim_2(coord, m_y_spline_representation.span_cview());
    }
 
    KOKKOS_FUNCTION double jacobian(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const
    {
        Matrix_2x2 J;
        jacobian_matrix(coord, J);
        return determinant(J);
    }
 
    KOKKOS_INLINE_FUNCTION Matrix_2x2 first_order_jacobian_matrix_r_rtheta(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const
    {
        Matrix_2x2 matrix;
        matrix[0][0]
                = m_spline_evaluator.deriv_dim_1(coord, m_x_spline_representation.span_cview());
        matrix[0][1]
                = m_spline_evaluator.deriv_1_and_2(coord, m_x_spline_representation.span_cview());
        matrix[1][0]
                = m_spline_evaluator.deriv_dim_1(coord, m_y_spline_representation.span_cview());
        matrix[1][1]
                = m_spline_evaluator.deriv_1_and_2(coord, m_y_spline_representation.span_cview());
        return matrix;
    }
 
    KOKKOS_INLINE_FUNCTION IdxRangeRTheta idx_range_singular_point() const
    {
        return m_idx_range_singular_point;
    }
 
    KOKKOS_INLINE_FUNCTION const ddc::Coordinate<X, Y> control_point(
            ddc::DiscreteElement<BSplineR, BSplineTheta> const& el) const
    {
        return ddc::Coordinate<X, Y>(m_x_spline_representation(el), m_y_spline_representation(el));
    }
};
 
 
namespace mapping_detail {
template <
        class X,
        class Y,
        class SplineEvaluator,
        class R,
        class Theta,
        class MemorySpace,
        class ExecSpace>
struct MappingAccessibility<
        ExecSpace,
        DiscreteToCartesian<X, Y, SplineEvaluator, R, Theta, MemorySpace>>
{
    static constexpr bool value = Kokkos::SpaceAccessibility<ExecSpace, MemorySpace>::accessible;
};
 
template <class X, class Y, class SplineEvaluator, class R, class Theta, class MemorySpace>
struct IsCurvilinear2DMapping<DiscreteToCartesian<X, Y, SplineEvaluator, R, Theta, MemorySpace>>
    : std::true_type
{
};
 
template <class X, class Y, class SplineEvaluator, class R, class Theta, class MemorySpace>
struct SingularOPointInvJacobian<DiscreteToCartesian<X, Y, SplineEvaluator, R, Theta, MemorySpace>>
    : std::true_type
{
};
 
} // namespace mapping_detail