jacobian.hpp

// SPDX-License-Identifier: MIT
#pragma once
#include <array>
#include <cassert>
 
#include <sll/view.hpp>
 
template <class PositionCoordinate>
class Jacobian
{
public:
    virtual KOKKOS_FUNCTION double jacobian(PositionCoordinate const& coord) const = 0;
 
    virtual KOKKOS_FUNCTION void jacobian_matrix(
            PositionCoordinate const& coord,
            Matrix_2x2& matrix) const = 0;
 
    virtual KOKKOS_FUNCTION double jacobian_11(PositionCoordinate const& coord) const = 0;
 
    virtual KOKKOS_FUNCTION double jacobian_12(PositionCoordinate const& coord) const = 0;
 
    virtual KOKKOS_FUNCTION double jacobian_21(PositionCoordinate const& coord) const = 0;
 
    virtual KOKKOS_FUNCTION double jacobian_22(PositionCoordinate const& coord) const = 0;
 
    virtual KOKKOS_FUNCTION void inv_jacobian_matrix(
            PositionCoordinate const& coord,
            Matrix_2x2& matrix) const = 0;
 
    virtual KOKKOS_FUNCTION double inv_jacobian_11(PositionCoordinate const& coord) const = 0;
 
    virtual KOKKOS_FUNCTION double inv_jacobian_12(PositionCoordinate const& coord) const = 0;
 
    virtual KOKKOS_FUNCTION double inv_jacobian_21(PositionCoordinate const& coord) const = 0;
 
    virtual KOKKOS_FUNCTION double inv_jacobian_22(PositionCoordinate const& coord) const = 0;
};
 
template <class PositionCoordinate>
class NonAnalyticalJacobian : public Jacobian<PositionCoordinate>
{
public:
    KOKKOS_FUNCTION double jacobian(PositionCoordinate const& coord) const final
    {
        const double j_rr = this->jacobian_11(coord);
        const double j_rp = this->jacobian_12(coord);
        const double j_pr = this->jacobian_21(coord);
        const double j_pp = this->jacobian_22(coord);
        return j_rr * j_pp - j_rp * j_pr;
    }
 
    KOKKOS_FUNCTION void inv_jacobian_matrix(PositionCoordinate const& coord, Matrix_2x2& matrix)
            const final
    {
        double jacob = jacobian(coord);
        assert(fabs(jacob) > 1e-15);
        matrix[0][0] = this->jacobian_22(coord) / jacob;
        matrix[0][1] = -this->jacobian_12(coord) / jacob;
        matrix[1][0] = -this->jacobian_21(coord) / jacob;
        matrix[1][1] = this->jacobian_11(coord) / jacob;
    }
 
    KOKKOS_FUNCTION double inv_jacobian_11(PositionCoordinate const& coord) const final
    {
        double jacob = jacobian(coord);
        assert(fabs(jacob) > 1e-15);
        return this->jacobian_22(coord) / jacob;
    }
 
    KOKKOS_FUNCTION double inv_jacobian_12(PositionCoordinate const& coord) const final
    {
        double jacob = jacobian(coord);
        assert(fabs(jacob) > 1e-15);
        return -this->jacobian_12(coord) / jacob;
    }
 
    KOKKOS_FUNCTION double inv_jacobian_21(PositionCoordinate const& coord) const final
    {
        double jacob = jacobian(coord);
        assert(fabs(jacob) > 1e-15);
        return -this->jacobian_21(coord) / jacob;
    }
 
    KOKKOS_FUNCTION double inv_jacobian_22(PositionCoordinate const& coord) const final
    {
        double jacob = jacobian(coord);
        assert(fabs(jacob) > 1e-15);
        return this->jacobian_11(coord) / jacob;
    }
};