inverse_jacobian_matrix.hpp

// SPDX-License-Identifier: MIT
#pragma once
 
#include <sll/view.hpp>
 
#include "mapping_tools.hpp"
 
template <class Mapping, class PositionCoordinate = typename Mapping::CoordArg>
class InverseJacobianMatrix
{
private:
    Mapping m_mapping;
 
public:
    KOKKOS_FUNCTION explicit InverseJacobianMatrix(Mapping const& mapping) : m_mapping(mapping) {}
 
    KOKKOS_INLINE_FUNCTION Matrix_2x2 operator()(PositionCoordinate const& coord) const
    {
        Matrix_2x2 matrix;
        if constexpr (has_2d_inv_jacobian_v<Mapping, PositionCoordinate>) {
            m_mapping.inv_jacobian_matrix(coord, matrix);
        } else {
            static_assert(has_2d_jacobian_v<Mapping, PositionCoordinate>);
            double jacob = m_mapping.jacobian(coord);
            assert(fabs(jacob) > 1e-15);
            matrix[0][0] = m_mapping.jacobian_22(coord) / jacob;
            matrix[0][1] = -m_mapping.jacobian_12(coord) / jacob;
            matrix[1][0] = -m_mapping.jacobian_21(coord) / jacob;
            matrix[1][1] = m_mapping.jacobian_11(coord) / jacob;
        }
        return matrix;
    }
 
    KOKKOS_INLINE_FUNCTION double inv_jacobian_11(PositionCoordinate const& coord) const
    {
        if constexpr (has_2d_inv_jacobian_v<Mapping, PositionCoordinate>) {
            return m_mapping.inv_jacobian_11(coord);
        } else {
            static_assert(has_2d_jacobian_v<Mapping, PositionCoordinate>);
            double jacob = m_mapping.jacobian(coord);
            assert(fabs(jacob) > 1e-15);
            return m_mapping.jacobian_22(coord) / jacob;
        }
    }
 
    KOKKOS_INLINE_FUNCTION double inv_jacobian_12(PositionCoordinate const& coord) const
    {
        if constexpr (has_2d_inv_jacobian_v<Mapping, PositionCoordinate>) {
            return m_mapping.inv_jacobian_12(coord);
        } else {
            static_assert(has_2d_jacobian_v<Mapping, PositionCoordinate>);
            double jacob = m_mapping.jacobian(coord);
            assert(fabs(jacob) > 1e-15);
            return -m_mapping.jacobian_12(coord) / jacob;
        }
    }
    KOKKOS_INLINE_FUNCTION double inv_jacobian_21(PositionCoordinate const& coord) const
    {
        if constexpr (has_2d_inv_jacobian_v<Mapping, PositionCoordinate>) {
            return m_mapping.inv_jacobian_21(coord);
        } else {
            static_assert(has_2d_jacobian_v<Mapping, PositionCoordinate>);
            double jacob = m_mapping.jacobian(coord);
            assert(fabs(jacob) > 1e-15);
            return -m_mapping.jacobian_21(coord) / jacob;
        }
    }
 
    KOKKOS_INLINE_FUNCTION double inv_jacobian_22(PositionCoordinate const& coord) const
    {
        if constexpr (has_2d_inv_jacobian_v<Mapping, PositionCoordinate>) {
            return m_mapping.inv_jacobian_22(coord);
        } else {
            static_assert(has_2d_jacobian_v<Mapping, PositionCoordinate>);
            double jacob = m_mapping.jacobian(coord);
            assert(fabs(jacob) > 1e-15);
            return m_mapping.jacobian_11(coord) / jacob;
        }
    }
};