diff --git a/src/geometryRTheta/advection_field/advection_field_rp.hpp b/src/geometryRTheta/advection_field/advection_field_rp.hpp
index 578dce584..80d5b6c79 100644
--- a/src/geometryRTheta/advection_field/advection_field_rp.hpp
+++ b/src/geometryRTheta/advection_field/advection_field_rp.hpp
@@ -117,7 +117,7 @@ class AdvectionFieldFinder
      *      The parameter @f$ \varepsilon @f$ for the linearization of the
      *      electric field.
      */
-    AdvectionFieldFinder(Mapping const& mapping, double const epsilon = 1e-12)
+    explicit AdvectionFieldFinder(Mapping const& mapping, double const epsilon = 1e-12)
         : m_mapping(mapping)
         , m_polar_spline_evaluator(ddc::NullExtrapolationRule())
         , m_spline_evaluator {ddc::NullExtrapolationRule(), ddc::NullExtrapolationRule(), ddc::PeriodicExtrapolationRule<Theta>(), ddc::PeriodicExtrapolationRule<Theta>()}
@@ -245,6 +245,8 @@ class AdvectionFieldFinder
                 get_const_field(coords),
                 get_const_field(electrostatic_potential_coef));
 
+        InverseJacobianMatrix<Mapping, CoordRTheta> inv_jacobian_matrix(m_mapping);
+
         // > computation of the electric field
         ddc::for_each(grid, [&](IdxRTheta const irp) {
             double const r = ddc::coordinate(ddc::select<GridR>(irp));
@@ -253,8 +255,7 @@ class AdvectionFieldFinder
             if (r > m_epsilon) {
                 CoordRTheta const coord_rp(r, th);
 
-                Matrix_2x2 inv_J; // Inverse Jacobian matrix
-                m_mapping.inv_jacobian_matrix(coord_rp, inv_J);
+                Matrix_2x2 inv_J = inv_jacobian_matrix(coord_rp);
 
                 // Gradiant of phi in the physical index range (Cartesian index range)
                 double const deriv_x_phi
@@ -304,8 +305,7 @@ class AdvectionFieldFinder
                 // --- Value at r = m_epsilon:
                 CoordRTheta const coord_rp_epsilon(m_epsilon, th);
 
-                Matrix_2x2 inv_J_eps; // Jacobian matrix
-                m_mapping.inv_jacobian_matrix(coord_rp_epsilon, inv_J_eps);
+                Matrix_2x2 inv_J_eps = inv_jacobian_matrix(coord_rp_epsilon);
 
                 double const deriv_r_phi_epsilon = evaluator.deriv_dim_1(
                         coord_rp_epsilon,
diff --git a/tests/geometryRTheta/polar_poisson/test_cases.hpp b/tests/geometryRTheta/polar_poisson/test_cases.hpp
index 9a423a065..d74d5d098 100644
--- a/tests/geometryRTheta/polar_poisson/test_cases.hpp
+++ b/tests/geometryRTheta/polar_poisson/test_cases.hpp
@@ -116,6 +116,9 @@ class CurvilinearSolution : public PoissonSolution<CurvilinearToCartesian>
 template <class CurvilinearToCartesian>
 class CartesianSolution : public PoissonSolution<CurvilinearToCartesian>
 {
+private:
+    InverseJacobianMatrix<CurvilinearToCartesian, Coord<R, Theta>> m_inverse_jacobian;
+
 public:
     /**
      * @brief Instantiate a CartesianSolution.
@@ -126,6 +129,7 @@ class CartesianSolution : public PoissonSolution<CurvilinearToCartesian>
      */
     explicit CartesianSolution(CurvilinearToCartesian const& coordinate_converter)
         : PoissonSolution<CurvilinearToCartesian>(coordinate_converter)
+        , m_inverse_jacobian(coordinate_converter)
     {
     }
 
@@ -156,9 +160,7 @@ class CartesianSolution : public PoissonSolution<CurvilinearToCartesian>
         const double y = ddc::get<Y>(cart_coord);
 
         const double r = ddc::get<R>(coord);
-        const double dx_r
-                = PoissonSolution<CurvilinearToCartesian>::m_coordinate_converter.inv_jacobian_11(
-                        coord);
+        const double dx_r = m_inverse_jacobian.inv_jacobian_11(coord);
 
         const double C = 1e-4 * ipow(2., 12);
         return C * std::sin(2 * M_PI * y)
@@ -182,9 +184,7 @@ class CartesianSolution : public PoissonSolution<CurvilinearToCartesian>
         const double y = ddc::get<Y>(cart_coord);
 
         const double r = ddc::get<R>(coord);
-        const double dy_r
-                = PoissonSolution<CurvilinearToCartesian>::m_coordinate_converter.inv_jacobian_12(
-                        coord);
+        const double dy_r = m_inverse_jacobian.inv_jacobian_12(coord);
 
         const double C = 1e-4 * ipow(2., 12);
         return C * std::cos(2 * M_PI * x)
diff --git a/vendor/sll/include/sll/mapping/cartesian_to_pseudo_cartesian.hpp b/vendor/sll/include/sll/mapping/cartesian_to_pseudo_cartesian.hpp
index 924a2954d..b496cae6c 100644
--- a/vendor/sll/include/sll/mapping/cartesian_to_pseudo_cartesian.hpp
+++ b/vendor/sll/include/sll/mapping/cartesian_to_pseudo_cartesian.hpp
@@ -3,6 +3,9 @@
 
 #include <sll/view.hpp>
 
+#include "inverse_jacobian_matrix.hpp"
+#include "mapping_tools.hpp"
+
 /**
  * @brief A class describing a mapping from a Cartesian domain to a pseudo-Cartesian.
  * This operation is carried out using 2 curvilinear to Cartesian mappings.
@@ -48,6 +51,11 @@ class CartesianToPseudoCartesian
     using Theta = typename MappingToCartesian::curvilinear_tag_theta;
     using CoordRTheta = ddc::Coordinate<R, Theta>;
 
+    static_assert(is_2d_mapping_v<MappingToCartesian>);
+    static_assert(is_2d_mapping_v<MappingToPseudoCartesian>);
+    static_assert(has_2d_jacobian_v<MappingToCartesian, CoordRTheta>);
+    static_assert(has_2d_jacobian_v<MappingToPseudoCartesian, CoordRTheta>);
+
 private:
     MappingToCartesian m_mapping_to_cartesian;
     MappingToPseudoCartesian m_mapping_to_pseudo_cartesian;
@@ -91,13 +99,14 @@ class CartesianToPseudoCartesian
 
         Matrix_2x2 inv_jacobian_from_cartesian;
         Matrix_2x2 jacobian_from_pseudo_cartesian;
+        InverseJacobianMatrix inv_jacobian_cartesian(m_mapping_to_cartesian);
         if (r < m_epsilon) {
             Matrix_2x2 J_0;
             m_mapping_to_cartesian.to_pseudo_cartesian_jacobian_center_matrix(J_0);
 
             CoordRTheta coord_eps(m_epsilon, ddc::get<Theta>(coord_rtheta));
 
-            m_mapping_to_cartesian.inv_jacobian_matrix(coord_eps, inv_jacobian_from_cartesian);
+            inv_jacobian_from_cartesian = inv_jacobian_cartesian(coord_eps);
             m_mapping_to_pseudo_cartesian
                     .jacobian_matrix(coord_eps, jacobian_from_pseudo_cartesian);
             Matrix_2x2 J_eps = mat_mul(jacobian_from_pseudo_cartesian, inv_jacobian_from_cartesian);
@@ -107,7 +116,7 @@ class CartesianToPseudoCartesian
             J[1][0] = (1 - r / m_epsilon) * J_0[1][0] + r / m_epsilon * J_eps[1][0];
             J[1][1] = (1 - r / m_epsilon) * J_0[1][1] + r / m_epsilon * J_eps[1][1];
         } else {
-            m_mapping_to_cartesian.inv_jacobian_matrix(coord_rtheta, inv_jacobian_from_cartesian);
+            inv_jacobian_from_cartesian = inv_jacobian_cartesian(coord_rtheta);
             m_mapping_to_pseudo_cartesian
                     .jacobian_matrix(coord_rtheta, jacobian_from_pseudo_cartesian);
             J = mat_mul(jacobian_from_pseudo_cartesian, inv_jacobian_from_cartesian);
@@ -198,6 +207,7 @@ class CartesianToPseudoCartesian
 
         Matrix_2x2 jacobian_from_cartesian;
         Matrix_2x2 inv_jacobian_from_pseudo_cartesian;
+        InverseJacobianMatrix inv_jacobian_pseudo_cartesian(m_mapping_to_pseudo_cartesian);
         if (r < m_epsilon) {
             Matrix_2x2 J_0;
             m_mapping_to_cartesian.to_pseudo_cartesian_jacobian_center_matrix(J_0);
@@ -206,8 +216,7 @@ class CartesianToPseudoCartesian
             CoordRTheta coord_eps(m_epsilon, ddc::get<Theta>(coord_rtheta));
 
             m_mapping_to_cartesian.jacobian_matrix(coord_eps, jacobian_from_cartesian);
-            m_mapping_to_pseudo_cartesian
-                    .inv_jacobian_matrix(coord_eps, inv_jacobian_from_pseudo_cartesian);
+            inv_jacobian_from_pseudo_cartesian = inv_jacobian_pseudo_cartesian(coord_eps);
             Matrix_2x2 J_eps = mat_mul(inv_jacobian_from_pseudo_cartesian, jacobian_from_cartesian);
 
             J[0][0] = (1 - r / m_epsilon) * J_0_inv[0][0] + r / m_epsilon * J_eps[0][0];
@@ -216,8 +225,7 @@ class CartesianToPseudoCartesian
             J[1][1] = (1 - r / m_epsilon) * J_0_inv[1][1] + r / m_epsilon * J_eps[1][1];
         } else {
             m_mapping_to_cartesian.jacobian_matrix(coord_rtheta, jacobian_from_cartesian);
-            m_mapping_to_pseudo_cartesian
-                    .inv_jacobian_matrix(coord_rtheta, inv_jacobian_from_pseudo_cartesian);
+            inv_jacobian_from_pseudo_cartesian = inv_jacobian_pseudo_cartesian(coord_rtheta);
             J = mat_mul(jacobian_from_cartesian, inv_jacobian_from_pseudo_cartesian);
         }
     }
diff --git a/vendor/sll/include/sll/mapping/discrete_to_cartesian.hpp b/vendor/sll/include/sll/mapping/discrete_to_cartesian.hpp
index 02ed36ce8..6517f865d 100644
--- a/vendor/sll/include/sll/mapping/discrete_to_cartesian.hpp
+++ b/vendor/sll/include/sll/mapping/discrete_to_cartesian.hpp
@@ -4,11 +4,11 @@
 
 #include <ddc/ddc.hpp>
 
+#include <sll/math_tools.hpp>
 #include <sll/view.hpp>
 
 #include "coordinate_converter.hpp"
 #include "curvilinear2d_to_cartesian.hpp"
-#include "jacobian.hpp"
 #include "mapping_tools.hpp"
 #include "pseudo_cartesian_compatible_mapping.hpp"
 
@@ -35,7 +35,6 @@ template <
         class MemorySpace = typename SplineEvaluator::memory_space>
 class DiscreteToCartesian
     : public CoordinateConverter<ddc::Coordinate<R, Theta>, ddc::Coordinate<X, Y>>
-    , public NonAnalyticalJacobian<ddc::Coordinate<R, Theta>>
     , public PseudoCartesianCompatibleMapping
     , public Curvilinear2DToCartesian<X, Y, R, Theta>
 {
@@ -166,7 +165,7 @@ class DiscreteToCartesian
      */
     KOKKOS_FUNCTION void jacobian_matrix(
             ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord,
-            Matrix_2x2& matrix) const final
+            Matrix_2x2& matrix) const
     {
         matrix[0][0]
                 = m_spline_evaluator.deriv_dim_1(coord, m_x_spline_representation.span_cview());
@@ -195,7 +194,7 @@ class DiscreteToCartesian
      * @see SplineEvaluator2D
      */
     KOKKOS_FUNCTION double jacobian_11(
-            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const final
+            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());
     }
@@ -217,7 +216,7 @@ class DiscreteToCartesian
      * @see SplineEvaluator2D
      */
     KOKKOS_FUNCTION double jacobian_12(
-            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const final
+            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());
     }
@@ -239,7 +238,7 @@ class DiscreteToCartesian
      * @see SplineEvaluator2D
      */
     KOKKOS_FUNCTION double jacobian_21(
-            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const final
+            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());
     }
@@ -261,11 +260,26 @@ class DiscreteToCartesian
      * @see SplineEvaluator2D
      */
     KOKKOS_FUNCTION double jacobian_22(
-            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const final
+            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());
     }
 
+    /**
+     * @brief Compute the Jacobian, the determinant of the Jacobian matrix of the mapping.
+     *
+     * @param[in] coord
+     *          The coordinate where we evaluate the Jacobian.
+     *
+     * @return A double with the value of the determinant of the Jacobian matrix.
+     */
+    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);
+    }
 
 
     /**
diff --git a/vendor/sll/include/sll/mapping/inverse_jacobian_matrix.hpp b/vendor/sll/include/sll/mapping/inverse_jacobian_matrix.hpp
new file mode 100644
index 000000000..d6f877c59
--- /dev/null
+++ b/vendor/sll/include/sll/mapping/inverse_jacobian_matrix.hpp
@@ -0,0 +1,148 @@
+// SPDX-License-Identifier: MIT
+#pragma once
+
+#include <sll/view.hpp>
+
+#include "mapping_tools.hpp"
+
+/**
+ * A class to calculate the inverse of the Jacobian matrix.
+ * If specialised methods are available then these are used by the class.
+ * Otherwise the inverse is calculated from the Jacobian matrix.
+ *
+ * @tparam Mapping The mapping whose inverse we are interested in.
+ * @tparam PositionCoordinate The coordinate system in which the inverse should be calculated.
+ */
+template <class Mapping, class PositionCoordinate = typename Mapping::CoordArg>
+class InverseJacobianMatrix
+{
+private:
+    Mapping m_mapping;
+
+public:
+    /**
+     * @brief A constructor for the InverseJacobianMatrix.
+     * @param[in] mapping The mapping whose inverse we are interested in.
+     */
+    KOKKOS_FUNCTION explicit InverseJacobianMatrix(Mapping const& mapping) : m_mapping(mapping) {}
+
+    /**
+     * @brief Compute full inverse Jacobian matrix.
+     *
+     * For some computations, we need the complete inverse Jacobian matrix or just the
+     * coefficients.
+     * The coefficients can be given indendently with the functions
+     * inv_jacobian_11, inv_jacobian_12, inv_jacobian_21 and inv_jacobian_22.
+     *
+     * @param[in] coord The coordinate where we evaluate the Jacobian matrix.
+     * @returns The inverse Jacobian matrix returned.
+     *
+     * @see inv_jacobian_11
+     * @see inv_jacobian_12
+     * @see inv_jacobian_21
+     * @see inv_jacobian_22
+     */
+    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;
+    }
+
+    /**
+     * @brief Compute the (1,1) coefficient of the inverse Jacobian matrix.
+     *
+     * Be careful because not all mappings are invertible, especially at the center point.
+     *
+     * @param[in] coord
+     * 				The coordinate where we evaluate the inverse Jacobian matrix.
+     *
+     * @return A double with the value of the (1,1) coefficient of the inverse Jacobian 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;
+        }
+    }
+
+    /**
+     * @brief Compute the (1,2) coefficient of the inverse Jacobian matrix.
+     *
+     * Be careful because not all mappings are invertible, especially at the center point.
+     *
+     * @param[in] coord
+     * 				The coordinate where we evaluate the inverse Jacobian matrix.
+     *
+     * @return A double with the value of the (1,2) coefficient of the inverse Jacobian matrix.
+     */
+    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;
+        }
+    }
+    /**
+     * @brief Compute the (2,1) coefficient of the inverse Jacobian matrix.
+     *
+     * Be careful because not all mappings are invertible, especially at the center point.
+     *
+     * @param[in] coord
+     * 				The coordinate where we evaluate the inverse Jacobian matrix.
+     *
+     * @return A double with the value of the (2,1) coefficient of the inverse Jacobian matrix.
+     */
+    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;
+        }
+    }
+
+    /**
+     * @brief Compute the (2,2) coefficient of the inverse Jacobian matrix.
+     *
+     * Be careful because not all mappings are invertible, especially at the center point.
+     *
+     * @param[in] coord
+     * 				The coordinate where we evaluate the inverse Jacobian matrix.
+     *
+     * @return A double with the value of the (2,2) coefficient of the inverse Jacobian matrix.
+     */
+    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;
+        }
+    }
+};
diff --git a/vendor/sll/include/sll/mapping/jacobian.hpp b/vendor/sll/include/sll/mapping/jacobian.hpp
deleted file mode 100644
index 07031dc14..000000000
--- a/vendor/sll/include/sll/mapping/jacobian.hpp
+++ /dev/null
@@ -1,300 +0,0 @@
-// SPDX-License-Identifier: MIT
-#pragma once
-#include <array>
-#include <cassert>
-
-#include <sll/view.hpp>
-
-/**
- * An operator to calculate the Jacobian matrix and its inverse.
- * All operators which can calculate terms of the Jacobian matrix should inherit from this class.
- *
- * @tparam PositionCoordinate The type of the coordinate at which the Jacobian matrix can be calculated.
- */
-template <class PositionCoordinate>
-class Jacobian
-{
-public:
-    /**
-     * @brief Compute the Jacobian, the determinant of the Jacobian matrix of the mapping.
-     *
-     * @param[in] coord
-     * 			The coordinate where we evaluate the Jacobian.
-     *
-     * @return A double with the value of the determinant of the Jacobian matrix.
-     */
-    virtual KOKKOS_FUNCTION double jacobian(PositionCoordinate const& coord) const = 0;
-
-    /**
-     * @brief Compute full Jacobian matrix.
-     *
-     * For some computations, we need the complete Jacobian matrix or just the
-     * coefficients.
-     * The coefficients can be given indendently with the functions
-     * jacobian_11, jacobian_12,  jacobian_21 and jacobian_22.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the Jacobian matrix.
-     * @param[out] matrix
-     * 				The Jacobian matrix returned.
-     *
-     *
-     * @see Jacobian::jacobian_11
-     * @see Jacobian::jacobian_12
-     * @see Jacobian::jacobian_21
-     * @see Jacobian::jacobian_22
-     */
-    virtual KOKKOS_FUNCTION void jacobian_matrix(
-            PositionCoordinate const& coord,
-            Matrix_2x2& matrix) const = 0;
-
-    /**
-     * @brief Compute the (1,1) coefficient of the Jacobian matrix.
-     *
-     * For a mapping given by @f$ \mathcal{F} : (r,\theta)\mapsto (x,y) @f$, the
-     * (1,1) coefficient of the Jacobian matrix is given by @f$ \frac{\partial x}{\partial r} @f$.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the Jacobian matrix.
-     *
-     * @return A double with the value of the (1,1) coefficient of the Jacobian matrix.
-     */
-    virtual KOKKOS_FUNCTION double jacobian_11(PositionCoordinate const& coord) const = 0;
-
-    /**
-     * @brief Compute the (1,2) coefficient of the Jacobian matrix.
-     *
-     * For a mapping given by @f$ \mathcal{F} : (r,\theta)\mapsto (x,y) @f$, the
-     * (1,2) coefficient of the Jacobian matrix is given by @f$ \frac{\partial x}{\partial \theta} @f$.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the Jacobian matrix.
-     *
-     * @return A double with the value of the (1,2) coefficient of the Jacobian matrix.
-     */
-    virtual KOKKOS_FUNCTION double jacobian_12(PositionCoordinate const& coord) const = 0;
-
-    /**
-     * @brief Compute the (2,1) coefficient of the Jacobian matrix.
-     *
-     * For a mapping given by @f$ \mathcal{F} : (r,\theta)\mapsto (x,y) @f$, the
-     * (2,1) coefficient of the Jacobian matrix is given by @f$ \frac{\partial y}{\partial r} @f$.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the Jacobian matrix. .
-     *
-     * @return A double with the value of the (2,1) coefficient of the Jacobian matrix.
-     */
-    virtual KOKKOS_FUNCTION double jacobian_21(PositionCoordinate const& coord) const = 0;
-
-    /**
-     * @brief Compute the (2,2) coefficient of the Jacobian matrix.
-     *
-     * For a mapping given by @f$ \mathcal{F} : (r,\theta)\mapsto (x,y) @f$, the
-     * (2,2) coefficient of the Jacobian matrix is given by @f$ \frac{\partial y}{\partial \theta} @f$.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the Jacobian matrix.
-     *
-     * @return A double with the value of the (2,2) coefficient of the Jacobian matrix.
-     */
-    virtual KOKKOS_FUNCTION double jacobian_22(PositionCoordinate const& coord) const = 0;
-
-    /**
-     * @brief Compute full inverse Jacobian matrix.
-     *
-     * For some computations, we need the complete inverse Jacobian matrix or just the
-     * coefficients.
-     * The coefficients can be given indendently with the functions
-     * inv_jacobian_11, inv_jacobian_12, inv_jacobian_21 and inv_jacobian_22.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the Jacobian matrix.
-     * @param[out] matrix
-     * 				The inverse Jacobian matrix returned.
-     *
-     *
-     * @see Jacobian::inv_jacobian_11
-     * @see Jacobian::inv_jacobian_12
-     * @see Jacobian::inv_jacobian_21
-     * @see Jacobian::inv_jacobian_22
-     */
-    virtual KOKKOS_FUNCTION void inv_jacobian_matrix(
-            PositionCoordinate const& coord,
-            Matrix_2x2& matrix) const = 0;
-
-    /**
-     * @brief Compute the (1,1) coefficient of the inverse Jacobian matrix.
-     *
-     * Be careful because not all mappings are invertible, especially at the center point.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the inverse Jacobian matrix.
-     *
-     * @return A double with the value of the (1,1) coefficient of the inverse Jacobian matrix.
-     */
-    virtual KOKKOS_FUNCTION double inv_jacobian_11(PositionCoordinate const& coord) const = 0;
-
-    /**
-     * @brief Compute the (1,2) coefficient of the inverse Jacobian matrix.
-     *
-     * Be careful because not all mappings are invertible, especially at the center point.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the inverse Jacobian matrix.
-     *
-     * @return A double with the value of the (1,2) coefficient of the inverse Jacobian matrix.
-     */
-    virtual KOKKOS_FUNCTION double inv_jacobian_12(PositionCoordinate const& coord) const = 0;
-
-    /**
-     * @brief Compute the (2,1) coefficient of the inverse Jacobian matrix.
-     *
-     * Be careful because not all mappings are invertible, especially at the center point.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the inverse Jacobian matrix.
-     *
-     * @return A double with the value of the (2,1) coefficient of the inverse Jacobian matrix.
-     */
-    virtual KOKKOS_FUNCTION double inv_jacobian_21(PositionCoordinate const& coord) const = 0;
-
-    /**
-     * @brief Compute the (2,2) coefficient of the inverse Jacobian matrix.
-     *
-     * Be careful because not all mappings are invertible, especially at the center point.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the inverse Jacobian matrix.
-     *
-     * @return A double with the value of the (2,2) coefficient of the inverse Jacobian matrix.
-     */
-    virtual KOKKOS_FUNCTION double inv_jacobian_22(PositionCoordinate const& coord) const = 0;
-};
-
-/**
- * A specialisation of Jacobian to handle non-analytical terms. In this case the inverse and the determinant
- * are calculated from the Jacobian matrix in the same way regardless of the implementation of the calculation
- * of the Jacobian itself.
- *
- * @tparam PositionCoordinate The type of the coordinate at which the Jacobian matrix can be calculated.
- */
-template <class PositionCoordinate>
-class NonAnalyticalJacobian : public Jacobian<PositionCoordinate>
-{
-public:
-    /**
-     * @brief Compute the Jacobian, the determinant of the Jacobian matrix of the mapping.
-     *
-     * @param[in] coord
-     * 			The coordinate where we evaluate the Jacobian.
-     *
-     * @return A double with the value of the determinant of the Jacobian matrix.
-     */
-    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;
-    }
-
-    /**
-     * @brief Compute full inverse Jacobian matrix.
-     *
-     * For some computations, we need the complete inverse Jacobian matrix or just the
-     * coefficients.
-     * The coefficients can be given indendently with the functions
-     * inv_jacobian_11, inv_jacobian_12, inv_jacobian_21 and inv_jacobian_22.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the Jacobian matrix.
-     * @param[out] matrix
-     * 				The inverse Jacobian matrix returned.
-     *
-     *
-     * @see Jacobian::inv_jacobian_11
-     * @see Jacobian::inv_jacobian_12
-     * @see Jacobian::inv_jacobian_21
-     * @see Jacobian::inv_jacobian_22
-     */
-    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;
-    }
-
-    /**
-     * @brief Compute the (1,1) coefficient of the inverse Jacobian matrix.
-     *
-     * Be careful because not all mappings are invertible, especially at the center point.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the inverse Jacobian matrix.
-     *
-     * @return A double with the value of the (1,1) coefficient of the inverse Jacobian matrix.
-     */
-    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;
-    }
-
-    /**
-     * @brief Compute the (1,2) coefficient of the inverse Jacobian matrix.
-     *
-     * Be careful because not all mappings are invertible, especially at the center point.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the inverse Jacobian matrix.
-     *
-     * @return A double with the value of the (1,2) coefficient of the inverse Jacobian matrix.
-     */
-    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;
-    }
-
-    /**
-     * @brief Compute the (2,1) coefficient of the inverse Jacobian matrix.
-     *
-     * Be careful because not all mappings are invertible, especially at the center point.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the inverse Jacobian matrix.
-     *
-     * @return A double with the value of the (2,1) coefficient of the inverse Jacobian matrix.
-     */
-    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;
-    }
-
-    /**
-     * @brief Compute the (2,2) coefficient of the inverse Jacobian matrix.
-     *
-     * Be careful because not all mappings are invertible, especially at the center point.
-     *
-     * @param[in] coord
-     * 				The coordinate where we evaluate the inverse Jacobian matrix.
-     *
-     * @return A double with the value of the (2,2) coefficient of the inverse Jacobian matrix.
-     */
-    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;
-    }
-};
diff --git a/vendor/sll/include/sll/mapping/metric_tensor.hpp b/vendor/sll/include/sll/mapping/metric_tensor.hpp
index 12ab08101..3ab52033e 100644
--- a/vendor/sll/include/sll/mapping/metric_tensor.hpp
+++ b/vendor/sll/include/sll/mapping/metric_tensor.hpp
@@ -3,7 +3,6 @@
 
 #include <sll/view.hpp>
 
-#include "jacobian.hpp"
 #include "mapping_tools.hpp"
 
 /**
diff --git a/vendor/sll/tests/mapping_jacobian.cpp b/vendor/sll/tests/mapping_jacobian.cpp
index 12ca8b593..1d45f703e 100644
--- a/vendor/sll/tests/mapping_jacobian.cpp
+++ b/vendor/sll/tests/mapping_jacobian.cpp
@@ -8,6 +8,7 @@
 #include "sll/mapping/czarny_to_cartesian.hpp"
 #include "sll/mapping/discrete_mapping_builder.hpp"
 #include "sll/mapping/discrete_to_cartesian.hpp"
+#include "sll/mapping/inverse_jacobian_matrix.hpp"
 
 #include "test_utils.hpp"
 
@@ -142,13 +143,13 @@ void check_inverse(Matrix_2x2 matrix, Matrix_2x2 inv_matrix)
 /**
  * @brief A class for the Google tests.
  */
-class InverseJacobianMatrix : public testing::TestWithParam<std::tuple<std::size_t, std::size_t>>
+class InvJacobianMatrix : public testing::TestWithParam<std::tuple<std::size_t, std::size_t>>
 {
 };
 
 
 
-TEST_P(InverseJacobianMatrix, InverseMatrixCircMap)
+TEST_P(InvJacobianMatrix, InverseMatrixCircMap)
 {
     auto const [Nr, Nt] = GetParam();
     const CircularToCartesian<X, Y, R, Theta> mapping;
@@ -179,22 +180,21 @@ TEST_P(InverseJacobianMatrix, InverseMatrixCircMap)
     });
 
     static_assert(has_2d_jacobian_v<CircularToCartesian<X, Y, R, Theta>, CoordRTheta>);
-    static_assert(has_2d_inv_jacobian_v<CircularToCartesian<X, Y, R, Theta>, CoordRTheta>);
+    InverseJacobianMatrix inv_jacobian(mapping);
 
     // Test for each coordinates if the inv_Jacobian_matrix is the inverse of the Jacobian_matrix
     ddc::for_each(grid, [&](IdxRTheta const irp) {
         Matrix_2x2 Jacobian_matrix;
-        Matrix_2x2 inv_Jacobian_matrix;
+        Matrix_2x2 inv_Jacobian_matrix = inv_jacobian(coords(irp));
 
         mapping.jacobian_matrix(coords(irp), Jacobian_matrix);
-        mapping.inv_jacobian_matrix(coords(irp), inv_Jacobian_matrix);
 
         check_inverse(Jacobian_matrix, inv_Jacobian_matrix);
     });
 }
 
 
-TEST_P(InverseJacobianMatrix, InverseMatrixCzarMap)
+TEST_P(InvJacobianMatrix, InverseMatrixCzarMap)
 {
     auto const [Nr, Nt] = GetParam();
     const CzarnyToCartesian<X, Y, R, Theta> mapping(0.3, 1.4);
@@ -240,7 +240,7 @@ TEST_P(InverseJacobianMatrix, InverseMatrixCzarMap)
 }
 
 
-TEST_P(InverseJacobianMatrix, InverseMatrixDiscCzarMap)
+TEST_P(InvJacobianMatrix, InverseMatrixDiscCzarMap)
 {
     auto const [Nr, Nt] = GetParam();
     const CzarnyToCartesian<X, Y, R, Theta> analytical_mapping(0.3, 1.4);
@@ -293,7 +293,7 @@ TEST_P(InverseJacobianMatrix, InverseMatrixDiscCzarMap)
     DiscreteToCartesian mapping = mapping_builder();
 
     static_assert(has_2d_jacobian_v<decltype(mapping), CoordRTheta>);
-    static_assert(has_2d_inv_jacobian_v<decltype(mapping), CoordRTheta>);
+    InverseJacobianMatrix inv_jacobian(mapping);
 
     // Test for each coordinates if the inv_Jacobian_matrix is the inverse of the Jacobian_matrix
     ddc::for_each(grid, [&](IdxRTheta const irp) {
@@ -301,10 +301,9 @@ TEST_P(InverseJacobianMatrix, InverseMatrixDiscCzarMap)
         const double r = ddc::get<R>(coord_rp);
         if (fabs(r) > 1e-15) {
             Matrix_2x2 Jacobian_matrix;
-            Matrix_2x2 inv_Jacobian_matrix;
+            Matrix_2x2 inv_Jacobian_matrix = inv_jacobian(coord_rp);
 
             mapping.jacobian_matrix(coord_rp, Jacobian_matrix);
-            mapping.inv_jacobian_matrix(coord_rp, inv_Jacobian_matrix);
 
             check_inverse(Jacobian_matrix, inv_Jacobian_matrix);
         }
@@ -315,5 +314,5 @@ TEST_P(InverseJacobianMatrix, InverseMatrixDiscCzarMap)
 
 INSTANTIATE_TEST_SUITE_P(
         MyGroup,
-        InverseJacobianMatrix,
+        InvJacobianMatrix,
         testing::Combine(testing::Values<std::size_t>(64), testing::Values<std::size_t>(128)));
diff --git a/vendor/sll/tests/mapping_jacobian_matrix_coef.cpp b/vendor/sll/tests/mapping_jacobian_matrix_coef.cpp
index fd9c1d142..e9ad28378 100644
--- a/vendor/sll/tests/mapping_jacobian_matrix_coef.cpp
+++ b/vendor/sll/tests/mapping_jacobian_matrix_coef.cpp
@@ -8,6 +8,7 @@
 #include "sll/mapping/czarny_to_cartesian.hpp"
 #include "sll/mapping/discrete_mapping_builder.hpp"
 #include "sll/mapping/discrete_to_cartesian.hpp"
+#include "sll/mapping/inverse_jacobian_matrix.hpp"
 
 #include "test_utils.hpp"
 
@@ -244,7 +245,6 @@ TEST_P(JacobianMatrixAndJacobianCoefficients, MatrixCzarMap)
     });
 
     static_assert(has_2d_jacobian_v<decltype(mapping), CoordRTheta>);
-    static_assert(has_2d_inv_jacobian_v<decltype(mapping), CoordRTheta>);
 
     // Test for each coordinates if the coefficients defined by the coefficients functions
     //are the same as the coefficients in the matrix function.
@@ -262,16 +262,16 @@ TEST_P(JacobianMatrixAndJacobianCoefficients, MatrixCzarMap)
         check_matrix(Jacobian_matrix, Jacobian_matrix_coeff);
     });
 
+    InverseJacobianMatrix inv_jacobian(mapping);
     // --- for the inverseJacobian matrix:
     ddc::for_each(grid, [&](IdxRTheta const irp) {
-        Matrix_2x2 inv_Jacobian_matrix;
+        Matrix_2x2 inv_Jacobian_matrix = inv_jacobian(coords(irp));
         Matrix_2x2 inv_Jacobian_matrix_coeff;
 
-        mapping.inv_jacobian_matrix(coords(irp), inv_Jacobian_matrix);
-        inv_Jacobian_matrix_coeff[0][0] = mapping.inv_jacobian_11(coords(irp));
-        inv_Jacobian_matrix_coeff[0][1] = mapping.inv_jacobian_12(coords(irp));
-        inv_Jacobian_matrix_coeff[1][0] = mapping.inv_jacobian_21(coords(irp));
-        inv_Jacobian_matrix_coeff[1][1] = mapping.inv_jacobian_22(coords(irp));
+        inv_Jacobian_matrix_coeff[0][0] = inv_jacobian.inv_jacobian_11(coords(irp));
+        inv_Jacobian_matrix_coeff[0][1] = inv_jacobian.inv_jacobian_12(coords(irp));
+        inv_Jacobian_matrix_coeff[1][0] = inv_jacobian.inv_jacobian_21(coords(irp));
+        inv_Jacobian_matrix_coeff[1][1] = inv_jacobian.inv_jacobian_22(coords(irp));
 
         check_matrix(inv_Jacobian_matrix, inv_Jacobian_matrix_coeff);
     });
@@ -333,7 +333,7 @@ TEST_P(JacobianMatrixAndJacobianCoefficients, MatrixDiscCzarMap)
     DiscreteToCartesian mapping = mapping_builder();
 
     static_assert(has_2d_jacobian_v<decltype(mapping), CoordRTheta>);
-    static_assert(has_2d_inv_jacobian_v<decltype(mapping), CoordRTheta>);
+    InverseJacobianMatrix inv_jacobian(mapping);
 
     // Test for each coordinates if the coefficients defined by the coefficients functions
     //are the same as the coefficients in the matrix function.
@@ -355,14 +355,15 @@ TEST_P(JacobianMatrixAndJacobianCoefficients, MatrixDiscCzarMap)
 
 
             // --- for the inverse Jacobian matrix:
-            Matrix_2x2 inv_Jacobian_matrix;
+            Matrix_2x2 inv_Jacobian_matrix = inv_jacobian(coord_rp);
             Matrix_2x2 inv_Jacobian_matrix_coeff;
 
-            mapping.inv_jacobian_matrix(coord_rp, inv_Jacobian_matrix);
-            inv_Jacobian_matrix_coeff[0][0] = mapping.inv_jacobian_11(coord_rp);
-            inv_Jacobian_matrix_coeff[0][1] = mapping.inv_jacobian_12(coord_rp);
-            inv_Jacobian_matrix_coeff[1][0] = mapping.inv_jacobian_21(coord_rp);
-            inv_Jacobian_matrix_coeff[1][1] = mapping.inv_jacobian_22(coord_rp);
+            inv_Jacobian_matrix_coeff[0][0] = inv_jacobian.inv_jacobian_11(coord_rp);
+            inv_Jacobian_matrix_coeff[0][1] = inv_jacobian.inv_jacobian_12(coord_rp);
+            inv_Jacobian_matrix_coeff[1][0] = inv_jacobian.inv_jacobian_21(coord_rp);
+            inv_Jacobian_matrix_coeff[1][1] = inv_jacobian.inv_jacobian_22(coord_rp);
+
+            check_matrix(inv_Jacobian_matrix, inv_Jacobian_matrix_coeff);
         }
     });
 }