diff --git a/ParameterLib/CoordinateSystem.cpp b/ParameterLib/CoordinateSystem.cpp index c66fdc7c2b3..678249c6c94 100644 --- a/ParameterLib/CoordinateSystem.cpp +++ b/ParameterLib/CoordinateSystem.cpp @@ -12,26 +12,57 @@ #include #include -#include -#include +#include #include #include +#include "MathLib/FormattingUtils.h" #include "Parameter.h" namespace ParameterLib { static double const tolerance = std::numeric_limits::epsilon(); -#ifndef NDEBUG -static constexpr char error_info[] = - "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"; -static constexpr char normalization_error_info[] = - "The direction vector given by parameter {:s} for local_coordinate_system " - "is not normalized to unit length"; -#endif // NDEBUG +template +static void checkTransformationIsSON( + Eigen::Matrix 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::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::Identity()) + .eval(), + 2 * tolerance); + } +} + +template +static void checkNormalization(Eigen::MatrixBase 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 const& unit_direction) : _base{nullptr, nullptr, &unit_direction}, _has_implicit_base(true) @@ -100,28 +131,16 @@ CoordinateSystem::CoordinateSystem(Parameter const& e0, Eigen::Matrix getTransformationFromSingleBase2D( Parameter const& unit_direction, SpatialPosition const& pos) { - Eigen::Matrix t; -#ifndef NDEBUG - if (std::abs(Eigen::Map( - unit_direction(0 /* time independent */, pos).data()) - .norm() - - 1) > tolerance) - { - OGS_FATAL(normalization_error_info, unit_direction.name); - } -#endif // NDEBUG + auto const& normal = unit_direction(0 /* time independent */, pos); + checkNormalization(Eigen::Map(normal.data()), + unit_direction.name); - auto normal = unit_direction(0 /* time independent */, pos); + Eigen::Matrix t; // base 0: ( normal[1], -normal[0])^T // base 1: ( normal[0], normal[1])^T t << normal[1], normal[0], -normal[0], normal[1]; -#ifndef NDEBUG - if (std::abs(t.determinant() - 1) > tolerance) - { - OGS_FATAL(error_info, t.determinant(), tolerance); - } -#endif // NDEBUG + checkTransformationIsSON(t); return t; } @@ -147,39 +166,31 @@ Eigen::Matrix CoordinateSystem::transformation<2>( t.col(1) = Eigen::Map( (*_base[1])(0 /* time independent */, pos).data()); -#ifndef NDEBUG - if (std::abs(t.determinant() - 1) > tolerance) - { - OGS_FATAL(error_info, t.determinant(), tolerance); - } -#endif // NDEBUG + checkTransformationIsSON(t); return t; } Eigen::Matrix getTransformationFromSingleBase3D( Parameter const& unit_direction, SpatialPosition const& pos) { - auto e2 = unit_direction(0 /* time independent */, pos); -#ifndef NDEBUG - if (std::abs(Eigen::Map(e2.data()).norm() - 1.0) > - tolerance) - { - OGS_FATAL(normalization_error_info, unit_direction.name); - } -#endif + auto const& normal = unit_direction(0 /* time independent */, pos); + + Eigen::Matrix t; + auto e2 = t.col(2); + e2 = Eigen::Map(normal.data()); + checkNormalization(e2, unit_direction.name); // Find the id of the first non-zero component of e2: - auto const it = std::max_element(e2.begin(), e2.end(), - [](const double a, const double b) - { return std::abs(a) < std::abs(b); }); - const auto id = std::distance(e2.begin(), it); + int id; + e2.cwiseAbs().maxCoeff(&id); + // Get other two component ids: const auto id_a = (id + 1) % 3; const auto id_b = (id + 2) % 3; // Compute basis vector e1 orthogonal to e2 - using Vector3 = Eigen::Vector3d; - Vector3 e1(0.0, 0.0, 0.0); + auto e1 = t.col(1); + e1 = Eigen::Vector3d::Zero(); if (std::abs(e2[id_a]) < tolerance) { @@ -197,26 +208,12 @@ Eigen::Matrix getTransformationFromSingleBase3D( e1.normalize(); - auto eigen_mapped_e2 = Eigen::Map(e2.data()); - - auto e0 = e1.cross(eigen_mapped_e2); // |e0| = |e1 x e2| = |e1||e2|sin(theta) with theta the angle between e1 and // e2. Since |e1| = |e2| = 1.0, and theta = pi/2, we have |e0|=1. Therefore // e0 is normalized by nature. + t.col(0) = e1.cross(e2); - Eigen::Matrix t; - t.col(0) = e0; - t.col(1) = e1; - t.col(2) = eigen_mapped_e2; - -#ifndef NDEBUG - - if (std::abs(t.determinant() - 1) > tolerance) - { - OGS_FATAL(error_info, t.determinant(), tolerance); - } - -#endif // NDEBUG + checkTransformationIsSON(t); return t; } @@ -245,12 +242,7 @@ Eigen::Matrix CoordinateSystem::transformation<3>( t.col(2) = Eigen::Map( (*_base[2])(0 /* time independent */, pos).data()); -#ifndef NDEBUG - if (std::abs(t.determinant() - 1) > tolerance) - { - OGS_FATAL(error_info, t.determinant(), tolerance); - } -#endif // NDEBUG + checkTransformationIsSON(t); return t; } @@ -288,12 +280,7 @@ Eigen::Matrix CoordinateSystem::transformation_3d( Eigen::Matrix t = Eigen::Matrix::Identity(); t.template topLeftCorner<2, 2>() << e0[0], e1[0], e0[1], e1[1]; -#ifndef NDEBUG - if (std::abs(t.determinant() - 1) > tolerance) - { - OGS_FATAL(error_info, t.determinant(), tolerance); - } -#endif // NDEBUG + checkTransformationIsSON(t); return t; } diff --git a/Tests/MaterialLib/SolidModels/TestLinearElasticTransverseIsotropic.cpp b/Tests/MaterialLib/SolidModels/TestLinearElasticTransverseIsotropic.cpp index 287847f838c..ba568a3a382 100644 --- a/Tests/MaterialLib/SolidModels/TestLinearElasticTransverseIsotropic.cpp +++ b/Tests/MaterialLib/SolidModels/TestLinearElasticTransverseIsotropic.cpp @@ -241,9 +241,8 @@ TEST_F(LinearElasticTransverseIsotropic, test_agaist_ElasticOrthotropic) } { // 3D ParameterLib::ConstantParameter const e1{ - "e1", {-0.8191520442889918, 0.0, -0.573576436351046}}; - ParameterLib::ConstantParameter const e2{"e2", - {0.0, -1.0, 0.0}}; + "e1", {0.8191520442889918, 0.0, 0.573576436351046}}; + ParameterLib::ConstantParameter const e2{"e2", {0.0, 1.0, 0.0}}; ParameterLib::ConstantParameter const e3{ "e3", {-0.573576436351046, 0.0, 0.8191520442889918}}; ParameterLib::CoordinateSystem const coordinate_system{e1, e2, e3}; diff --git a/Tests/ParameterLib/TestCoordinateSystem.cpp b/Tests/ParameterLib/TestCoordinateSystem.cpp index 0d69d520a38..4cdc728e4ce 100644 --- a/Tests/ParameterLib/TestCoordinateSystem.cpp +++ b/Tests/ParameterLib/TestCoordinateSystem.cpp @@ -214,10 +214,10 @@ TEST(CoordinateSystem, test3D) { std::vector> parameters; - std::vector e0{-0.8191520442889918, 0.0, -0.573576436351046}; + std::vector e0{0.8191520442889918, 0.0, 0.573576436351046}; parameters.emplace_back( std::make_unique>("e0", e0)); - std::vector e1{0.0, -1.0, 0.0}; + std::vector e1{0.0, 1.0, 0.0}; parameters.emplace_back( std::make_unique>("e1", e1)); std::vector e2{-0.573576436351046, 0.0, 0.8191520442889918};