From c5c66068800a7508989c205d66aaf491dbbedf76 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 27 Jun 2024 11:23:07 +0200 Subject: [PATCH 1/8] [T] "Flip" the coordinate system to more usual one --- .../SolidModels/TestLinearElasticTransverseIsotropic.cpp | 5 ++--- Tests/ParameterLib/TestCoordinateSystem.cpp | 4 ++-- 2 files changed, 4 insertions(+), 5 deletions(-) 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}; From dfcd800ed77e1dd3c2092320985d8c6c485d73de Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 27 Jun 2024 12:05:24 +0200 Subject: [PATCH 2/8] [ParL] Extract transformation's det(T)==1 check --- ParameterLib/CoordinateSystem.cpp | 47 +++++++++++++------------------ 1 file changed, 20 insertions(+), 27 deletions(-) diff --git a/ParameterLib/CoordinateSystem.cpp b/ParameterLib/CoordinateSystem.cpp index c66fdc7c2b3..b293cac8db6 100644 --- a/ParameterLib/CoordinateSystem.cpp +++ b/ParameterLib/CoordinateSystem.cpp @@ -23,16 +23,26 @@ 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); + } +} + CoordinateSystem::CoordinateSystem(Parameter const& unit_direction) : _base{nullptr, nullptr, &unit_direction}, _has_implicit_base(true) { @@ -117,10 +127,7 @@ Eigen::Matrix getTransformationFromSingleBase2D( 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); - } + checkTransformationIsSON(t); #endif // NDEBUG return t; } @@ -148,10 +155,7 @@ Eigen::Matrix CoordinateSystem::transformation<2>( (*_base[1])(0 /* time independent */, pos).data()); #ifndef NDEBUG - if (std::abs(t.determinant() - 1) > tolerance) - { - OGS_FATAL(error_info, t.determinant(), tolerance); - } + checkTransformationIsSON(t); #endif // NDEBUG return t; } @@ -210,12 +214,7 @@ Eigen::Matrix getTransformationFromSingleBase3D( t.col(2) = eigen_mapped_e2; #ifndef NDEBUG - - if (std::abs(t.determinant() - 1) > tolerance) - { - OGS_FATAL(error_info, t.determinant(), tolerance); - } - + checkTransformationIsSON(t); #endif // NDEBUG return t; @@ -246,10 +245,7 @@ Eigen::Matrix CoordinateSystem::transformation<3>( (*_base[2])(0 /* time independent */, pos).data()); #ifndef NDEBUG - if (std::abs(t.determinant() - 1) > tolerance) - { - OGS_FATAL(error_info, t.determinant(), tolerance); - } + checkTransformationIsSON(t); #endif // NDEBUG return t; } @@ -289,10 +285,7 @@ Eigen::Matrix CoordinateSystem::transformation_3d( 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); - } + checkTransformationIsSON(t); #endif // NDEBUG return t; } From b0cf16b9679538143e080316a32e59ae92bf8971 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 27 Jun 2024 12:05:24 +0200 Subject: [PATCH 3/8] [ParL] Add orthogonality check of transformations --- ParameterLib/CoordinateSystem.cpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/ParameterLib/CoordinateSystem.cpp b/ParameterLib/CoordinateSystem.cpp index b293cac8db6..2d227da80ac 100644 --- a/ParameterLib/CoordinateSystem.cpp +++ b/ParameterLib/CoordinateSystem.cpp @@ -17,6 +17,7 @@ #include #include +#include "MathLib/FormattingUtils.h" #include "Parameter.h" namespace ParameterLib @@ -41,6 +42,18 @@ static void checkTransformationIsSON( "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); + } } CoordinateSystem::CoordinateSystem(Parameter const& unit_direction) From 61b105235b5978229fe92625d42b4bead444f4b7 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 28 Jun 2024 13:37:27 +0200 Subject: [PATCH 4/8] [ParL] Extract unity check for vectors --- ParameterLib/CoordinateSystem.cpp | 47 ++++++++++++++++--------------- 1 file changed, 24 insertions(+), 23 deletions(-) diff --git a/ParameterLib/CoordinateSystem.cpp b/ParameterLib/CoordinateSystem.cpp index 2d227da80ac..8cfb3b39475 100644 --- a/ParameterLib/CoordinateSystem.cpp +++ b/ParameterLib/CoordinateSystem.cpp @@ -23,11 +23,6 @@ namespace ParameterLib { static double const tolerance = std::numeric_limits::epsilon(); -#ifndef NDEBUG -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( @@ -56,6 +51,20 @@ static void checkTransformationIsSON( } } +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) { @@ -123,18 +132,14 @@ CoordinateSystem::CoordinateSystem(Parameter const& e0, Eigen::Matrix getTransformationFromSingleBase2D( Parameter const& unit_direction, SpatialPosition const& pos) { - Eigen::Matrix t; + auto const& normal = unit_direction(0 /* time independent */, pos); + #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); - } + checkNormalization(Eigen::Map(normal.data()), + unit_direction.name); #endif // NDEBUG - 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]; @@ -176,14 +181,12 @@ Eigen::Matrix CoordinateSystem::transformation<2>( Eigen::Matrix getTransformationFromSingleBase3D( Parameter const& unit_direction, SpatialPosition const& pos) { - auto e2 = unit_direction(0 /* time independent */, pos); + auto const& e2 = unit_direction(0 /* time independent */, pos); + + auto const eigen_mapped_e2 = Eigen::Map(e2.data()); #ifndef NDEBUG - if (std::abs(Eigen::Map(e2.data()).norm() - 1.0) > - tolerance) - { - OGS_FATAL(normalization_error_info, unit_direction.name); - } -#endif + checkNormalization(eigen_mapped_e2, unit_direction.name); +#endif // NDEBUG // Find the id of the first non-zero component of e2: auto const it = std::max_element(e2.begin(), e2.end(), @@ -214,8 +217,6 @@ 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 From 93f28439c6913119389ec416028e2a106e2f63d2 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 27 Jun 2024 14:24:48 +0200 Subject: [PATCH 5/8] [ParL] Enable transformation tests for non-debug No measurable run-time difference. --- ParameterLib/CoordinateSystem.cpp | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/ParameterLib/CoordinateSystem.cpp b/ParameterLib/CoordinateSystem.cpp index 8cfb3b39475..80d0abb61ff 100644 --- a/ParameterLib/CoordinateSystem.cpp +++ b/ParameterLib/CoordinateSystem.cpp @@ -133,20 +133,15 @@ Eigen::Matrix getTransformationFromSingleBase2D( Parameter const& unit_direction, SpatialPosition const& pos) { auto const& normal = unit_direction(0 /* time independent */, pos); - -#ifndef NDEBUG checkNormalization(Eigen::Map(normal.data()), unit_direction.name); -#endif // NDEBUG 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 checkTransformationIsSON(t); -#endif // NDEBUG return t; } @@ -172,9 +167,7 @@ Eigen::Matrix CoordinateSystem::transformation<2>( t.col(1) = Eigen::Map( (*_base[1])(0 /* time independent */, pos).data()); -#ifndef NDEBUG checkTransformationIsSON(t); -#endif // NDEBUG return t; } @@ -184,9 +177,7 @@ Eigen::Matrix getTransformationFromSingleBase3D( auto const& e2 = unit_direction(0 /* time independent */, pos); auto const eigen_mapped_e2 = Eigen::Map(e2.data()); -#ifndef NDEBUG checkNormalization(eigen_mapped_e2, unit_direction.name); -#endif // NDEBUG // Find the id of the first non-zero component of e2: auto const it = std::max_element(e2.begin(), e2.end(), @@ -227,9 +218,7 @@ Eigen::Matrix getTransformationFromSingleBase3D( t.col(1) = e1; t.col(2) = eigen_mapped_e2; -#ifndef NDEBUG checkTransformationIsSON(t); -#endif // NDEBUG return t; } @@ -258,9 +247,7 @@ Eigen::Matrix CoordinateSystem::transformation<3>( t.col(2) = Eigen::Map( (*_base[2])(0 /* time independent */, pos).data()); -#ifndef NDEBUG checkTransformationIsSON(t); -#endif // NDEBUG return t; } @@ -298,9 +285,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 checkTransformationIsSON(t); -#endif // NDEBUG return t; } From 2dd751c209c6967ee9adc8bc9daa85425d804206 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 28 Jun 2024 13:39:33 +0200 Subject: [PATCH 6/8] [ParL] Use Eigen's functions for Eigen vectors avoiding STL iterators. --- ParameterLib/CoordinateSystem.cpp | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/ParameterLib/CoordinateSystem.cpp b/ParameterLib/CoordinateSystem.cpp index 80d0abb61ff..4fd9f0850b7 100644 --- a/ParameterLib/CoordinateSystem.cpp +++ b/ParameterLib/CoordinateSystem.cpp @@ -12,8 +12,7 @@ #include #include -#include -#include +#include #include #include @@ -180,17 +179,15 @@ Eigen::Matrix getTransformationFromSingleBase3D( checkNormalization(eigen_mapped_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; + eigen_mapped_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); + Eigen::Vector3d e1 = Eigen::Vector3d::Zero(); if (std::abs(e2[id_a]) < tolerance) { From 2c95179dc26667345108ca230e45d3ac0729abb5 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 28 Jun 2024 13:53:39 +0200 Subject: [PATCH 7/8] [ParL] Use result transformation memory space --- ParameterLib/CoordinateSystem.cpp | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/ParameterLib/CoordinateSystem.cpp b/ParameterLib/CoordinateSystem.cpp index 4fd9f0850b7..76ec7d9bbc8 100644 --- a/ParameterLib/CoordinateSystem.cpp +++ b/ParameterLib/CoordinateSystem.cpp @@ -173,21 +173,24 @@ Eigen::Matrix CoordinateSystem::transformation<2>( Eigen::Matrix getTransformationFromSingleBase3D( Parameter const& unit_direction, SpatialPosition const& pos) { - auto const& e2 = unit_direction(0 /* time independent */, pos); + auto const& normal = unit_direction(0 /* time independent */, pos); - auto const eigen_mapped_e2 = Eigen::Map(e2.data()); - checkNormalization(eigen_mapped_e2, unit_direction.name); + 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: int id; - eigen_mapped_e2.cwiseAbs().maxCoeff(&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 - Eigen::Vector3d e1 = Eigen::Vector3d::Zero(); + auto e1 = t.col(1); + e1 = Eigen::Vector3d::Zero(); if (std::abs(e2[id_a]) < tolerance) { @@ -205,15 +208,10 @@ Eigen::Matrix getTransformationFromSingleBase3D( e1.normalize(); - 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. - - Eigen::Matrix t; - t.col(0) = e0; - t.col(1) = e1; - t.col(2) = eigen_mapped_e2; + t.col(0) = e1.cross(e2); checkTransformationIsSON(t); From f9cf152a174dd682fc30bcbf8a6a730e53118c6c Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 28 Jun 2024 15:51:01 +0200 Subject: [PATCH 8/8] [ParL] More template parameters for MSVC Current error: /ParameterLib/CoordinateSystem.cpp(280): error C2672: 'checkTransformationIsOrthonormal': no matching overloaded function found /ParameterLib/CoordinateSystem.cpp(280): error C2784: 'void ParameterLib::checkTransformationIsOrthonormal(const Eigen::Matrix &)': could not deduce template argument for 'const Eigen::Matrix &' from 'Eigen::Matrix' /ParameterLib/CoordinateSystem.cpp(26): note: see declaration of 'ParameterLib::checkTransformationIsOrthonormal' [609/644] Building CXX object Tests\CMakeFiles\testrunner.dir\Unity\unity_14_cxx.cxx.obj --- ParameterLib/CoordinateSystem.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ParameterLib/CoordinateSystem.cpp b/ParameterLib/CoordinateSystem.cpp index 76ec7d9bbc8..678249c6c94 100644 --- a/ParameterLib/CoordinateSystem.cpp +++ b/ParameterLib/CoordinateSystem.cpp @@ -25,7 +25,7 @@ static double const tolerance = std::numeric_limits::epsilon(); template static void checkTransformationIsSON( - Eigen::Matrix const& t) + Eigen::Matrix const& t) { if (std::abs(t.determinant() - 1) > tolerance) {