metric_tensor.hpp

#pragma once
#include <type_traits>
 
#include <sll/view.hpp>
 
#include "mapping_tools.hpp"
 
template <class Mapping, class PositionCoordinate>
class MetricTensor
{
    static_assert(is_2d_mapping_v<Mapping>);
    static_assert(has_2d_jacobian_v<Mapping, PositionCoordinate>);
 
public:
    using Matrix_2x2 = std::array<std::array<double, 2>, 2>;
 
private:
    Mapping m_mapping;
 
public:
    KOKKOS_FUNCTION MetricTensor(Mapping mapping) : m_mapping(mapping) {}
 
    KOKKOS_FUNCTION void operator()(Matrix_2x2& matrix, PositionCoordinate const& coord) const
    {
        const double J_11 = m_mapping.jacobian_11(coord);
        const double J_12 = m_mapping.jacobian_12(coord);
        const double J_21 = m_mapping.jacobian_21(coord);
        const double J_22 = m_mapping.jacobian_22(coord);
        matrix[0][0] = (J_11 * J_11 + J_21 * J_21);
        matrix[0][1] = (J_11 * J_12 + J_21 * J_22);
        matrix[1][0] = (J_11 * J_12 + J_21 * J_22);
        matrix[1][1] = (J_12 * J_12 + J_22 * J_22);
    }
 
    KOKKOS_FUNCTION void inverse(Matrix_2x2& matrix, PositionCoordinate const& coord) const
    {
        const double J_11 = m_mapping.jacobian_11(coord);
        const double J_12 = m_mapping.jacobian_12(coord);
        const double J_21 = m_mapping.jacobian_21(coord);
        const double J_22 = m_mapping.jacobian_22(coord);
        const double jacob_2 = m_mapping.jacobian(coord) * m_mapping.jacobian(coord);
        matrix[0][0] = (J_12 * J_12 + J_22 * J_22) / jacob_2;
        matrix[0][1] = (-J_11 * J_12 - J_21 * J_22) / jacob_2;
        matrix[1][0] = (-J_11 * J_12 - J_21 * J_22) / jacob_2;
        matrix[1][1] = (J_11 * J_11 + J_21 * J_21) / jacob_2;
    }
 
    KOKKOS_FUNCTION std::array<double, 2> to_covariant(
            std::array<double, 2> const& contravariant_vector,
            PositionCoordinate const& coord) const
    {
        Matrix_2x2 inv_metric_tensor;
        inverse(inv_metric_tensor, coord);
        std::array<double, 2> covariant_vector;
        covariant_vector[0] = inv_metric_tensor[0][0] * contravariant_vector[0]
                              + inv_metric_tensor[0][1] * contravariant_vector[1];
        covariant_vector[1] = inv_metric_tensor[1][0] * contravariant_vector[0]
                              + inv_metric_tensor[1][1] * contravariant_vector[1];
        return covariant_vector;
    }
};