geometry_helper.h

#pragma once

#include "../InputFunctions/domainGeometry.h"

#include <cmath>

inline void compute_jacobian_elements(const DomainGeometry& domain_geometry, double r, double theta, double coeff_alpha,
                                      double& arr, double& att, double& art, double& detDF)
{
    /* Calculate the elements of the Jacobian matrix for the transformation mapping */
    /* The Jacobian matrix is: */
    /* [Jrr, Jrt] */
    /* [Jtr, Jtt] */
    const double Jrr = domain_geometry.dFx_dr(r, theta);
    const double Jtr = domain_geometry.dFy_dr(r, theta);
    const double Jrt = domain_geometry.dFx_dt(r, theta);
    const double Jtt = domain_geometry.dFy_dt(r, theta);
    /* Compute the determinant of the Jacobian matrix */
    detDF = Jrr * Jtt - Jrt * Jtr;
    /* Compute the elements of the symmetric matrix: */
    /* 0.5 * alpha * DF^{-1} * DF^{-T} * |det(DF)| */
    /* which is represented by: */
    /*  [arr, 0.5*art]  */
    /*  [0.5*atr, att]  */
    arr = 0.5 * (Jtt * Jtt + Jrt * Jrt) * coeff_alpha / std::fabs(detDF);
    att = 0.5 * (Jtr * Jtr + Jrr * Jrr) * coeff_alpha / std::fabs(detDF);
    art = (-Jtt * Jtr - Jrt * Jrr) * coeff_alpha / std::fabs(detDF);
    /* Note that the inverse Jacobian matrix DF^{-1} is: */
    /*  1.0 / det(DF) *  */
    /*  [Jtt, -Jrt]      */
    /*  [-Jtr, Jrr]      */
}