CoordinateSystem.cpp

/**
 * \file
 * \copyright
 * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 *
 */

#include "CoordinateSystem.h"

#include <Eigen/Dense>
#include <Eigen/Geometry>
#include <cmath>
#include <limits>
#include <typeinfo>

#include "MathLib/FormattingUtils.h"
#include "Parameter.h"

namespace ParameterLib
{
static double const tolerance = std::numeric_limits<double>::epsilon();

template <int Dim>
static void checkTransformationIsSON(
    Eigen::Matrix<double, Dim, Dim, Eigen::ColMajor, Dim, Dim> const& t)
{
    if (std::abs(t.determinant() - 1) > tolerance)
    {
        OGS_FATAL(
            "The determinant of the coordinate system transformation matrix is "
            "'{:g}', which is not sufficiently close to unity with the "
            "tolerance of '{:g}'. Please adjust the accuracy of the local "
            "system bases",
            t.determinant(), tolerance);
    }
    if (((t * t.transpose() - Eigen::Matrix<double, Dim, Dim>::Identity())
             .array() > 2 * tolerance)
            .any())
    {
        OGS_FATAL(
            "The transformation is not orthogonal because the difference "
            "A*A^T - I is:\n{}\nand at least one component deviates from zero "
            "by more then '{:g}'.",
            (t * t.transpose() - Eigen::Matrix<double, Dim, Dim>::Identity())
                .eval(),
            2 * tolerance);
    }
}

template <typename Derived>
static void checkNormalization(Eigen::MatrixBase<Derived> const& vec,
                               std::string_view const parmeter_name)
{
    if (std::abs(vec.squaredNorm() - 1.0) > tolerance)
    {
        OGS_FATAL(
            "The direction vector given by parameter {:s} for the "
            "local_coordinate_system is not normalized to unit length. Vector "
            "norm is {:g} and |v|^2-1 = {:g}.",
            parmeter_name, vec.norm(), vec.squaredNorm() - 1.0);
    }
}

CoordinateSystem::CoordinateSystem(Parameter<double> const& unit_direction)
    : _base{nullptr, nullptr, &unit_direction}, _has_implicit_base(true)
{
    if (_base[2]->isTimeDependent())