circular_to_cartesian.hpp

// SPDX-License-Identifier: MIT
#pragma once
 
#include <cassert>
#include <cmath>
 
#include <ddc/ddc.hpp>
 
#include <sll/view.hpp>
 
#include "coordinate_converter.hpp"
#include "curvilinear2d_to_cartesian.hpp"
#include "mapping_tools.hpp"
#include "pseudo_cartesian_compatible_mapping.hpp"
 
template <class X, class Y, class R, class Theta>
class CircularToCartesian
    : public CoordinateConverter<ddc::Coordinate<X, Y>, ddc::Coordinate<R, Theta>>
    , public PseudoCartesianCompatibleMapping
    , public Curvilinear2DToCartesian<X, Y, R, Theta>
{
public:
    using cartesian_tag_x = typename Curvilinear2DToCartesian<X, Y, R, Theta>::cartesian_tag_x;
    using cartesian_tag_y = typename Curvilinear2DToCartesian<X, Y, R, Theta>::cartesian_tag_y;
    using curvilinear_tag_r = typename Curvilinear2DToCartesian<X, Y, R, Theta>::curvilinear_tag_r;
    using curvilinear_tag_theta =
            typename Curvilinear2DToCartesian<X, Y, R, Theta>::curvilinear_tag_theta;
 
    using CoordArg = ddc::Coordinate<R, Theta>;
    using CoordResult = ddc::Coordinate<X, Y>;
 
public:
    CircularToCartesian() = default;
 
    KOKKOS_FUNCTION CircularToCartesian(CircularToCartesian const& other) {}
 
    CircularToCartesian(CircularToCartesian&& x) = default;
 
    ~CircularToCartesian() = default;
 
    CircularToCartesian& operator=(CircularToCartesian const& x) = default;
 
    CircularToCartesian& operator=(CircularToCartesian&& x) = default;
 
    KOKKOS_FUNCTION ddc::Coordinate<X, Y> operator()(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
        const double x = r * Kokkos::cos(theta);
        const double y = r * Kokkos::sin(theta);
        return ddc::Coordinate<X, Y>(x, y);
    }
 
    KOKKOS_FUNCTION ddc::Coordinate<R, Theta> operator()(
            ddc::Coordinate<X, Y> const& coord) const final
    {
        const double x = ddc::get<X>(coord);
        const double y = ddc::get<Y>(coord);
        const double r = Kokkos::sqrt(x * x + y * y);
        const double theta = Kokkos::atan2(y, x);
        return ddc::Coordinate<R, Theta>(r, theta);
    }
 
    KOKKOS_FUNCTION double jacobian(ddc::Coordinate<R, Theta> const& coord) const
    {
        double r = ddc::get<R>(coord);
        return r;
    }
 
 
    KOKKOS_FUNCTION void jacobian_matrix(ddc::Coordinate<R, Theta> const& coord, Matrix_2x2& matrix)
            const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
        matrix[0][0] = Kokkos::cos(theta);
        matrix[0][1] = -r * Kokkos::sin(theta);
        matrix[1][0] = Kokkos::sin(theta);
        matrix[1][1] = r * Kokkos::cos(theta);
    }
 
    KOKKOS_FUNCTION double jacobian_11(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double theta = ddc::get<Theta>(coord);
        return Kokkos::cos(theta);
    }
 
    KOKKOS_FUNCTION double jacobian_12(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
        return -r * Kokkos::sin(theta);
    }
 
    KOKKOS_FUNCTION double jacobian_21(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double theta = ddc::get<Theta>(coord);
        return Kokkos::sin(theta);
    }
 
    KOKKOS_FUNCTION double jacobian_22(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
        return r * Kokkos::cos(theta);
    }
 
 
    KOKKOS_FUNCTION void inv_jacobian_matrix(
            ddc::Coordinate<R, Theta> const& coord,
            Matrix_2x2& matrix) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
        assert(fabs(r) >= 1e-15);
        matrix[0][0] = Kokkos::cos(theta);
        matrix[0][1] = Kokkos::sin(theta);
        matrix[1][0] = -1 / r * Kokkos::sin(theta);
        matrix[1][1] = 1 / r * Kokkos::cos(theta);
    }
 
    KOKKOS_FUNCTION double inv_jacobian_11(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double theta = ddc::get<Theta>(coord);
        return Kokkos::cos(theta);
    }
 
    KOKKOS_FUNCTION double inv_jacobian_12(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double theta = ddc::get<Theta>(coord);
        return Kokkos::sin(theta);
    }
 
    KOKKOS_FUNCTION double inv_jacobian_21(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
        assert(fabs(r) >= 1e-15);
        return -1 / r * Kokkos::sin(theta);
    }
 
    KOKKOS_FUNCTION double inv_jacobian_22(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
        assert(fabs(r) >= 1e-15);
        return 1 / r * Kokkos::cos(theta);
    }
 
 
 
    KOKKOS_FUNCTION void to_pseudo_cartesian_jacobian_center_matrix(Matrix_2x2& matrix) const final
    {
        matrix[0][0] = 1.;
        matrix[0][1] = 0.;
        matrix[1][0] = 0.;
        matrix[1][1] = 1.;
    }
 
    KOKKOS_FUNCTION double to_pseudo_cartesian_jacobian_11_center() const final
    {
        return 1.;
    }
 
    KOKKOS_FUNCTION double to_pseudo_cartesian_jacobian_12_center() const final
    {
        return 0.;
    }
 
    KOKKOS_FUNCTION double to_pseudo_cartesian_jacobian_21_center() const final
    {
        return 0.;
    }
 
    KOKKOS_FUNCTION double to_pseudo_cartesian_jacobian_22_center() const final
    {
        return 1.;
    }
};
 
 
namespace mapping_detail {
template <class X, class Y, class R, class Theta, class ExecSpace>
struct MappingAccessibility<ExecSpace, CircularToCartesian<X, Y, R, Theta>> : std::true_type
{
};
} // namespace mapping_detail