Skip to content

Commit

Permalink
Merge branch 'OrthonormalCheckForLocalCoordinateSystem' into 'master'
Browse files Browse the repository at this point in the history 
Additional check for coordinate system and small refactorings.

See merge request ogs/ogs!5033
  • Loading branch information
@endJunction
endJunction committed Jul 1, 2024
2 parents eb7959b + f9cf152 commit 86c4348
Show file tree
Hide file tree
Showing 3 changed files with 67 additions and 81 deletions.
139 changes: 63 additions & 76 deletions ParameterLib/CoordinateSystem.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,26 +12,57 @@

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

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

namespace ParameterLib
{
static double const tolerance = std::numeric_limits<double>::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 <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)
Expand Down Expand Up @@ -100,28 +131,16 @@ CoordinateSystem::CoordinateSystem(Parameter<double> const& e0,
Eigen::Matrix<double, 2, 2> getTransformationFromSingleBase2D(
Parameter<double> const& unit_direction, SpatialPosition const& pos)
{
Eigen::Matrix<double, 2, 2> t;
#ifndef NDEBUG
if (std::abs(Eigen::Map<Eigen::Vector2d>(
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<Eigen::Vector2d const>(normal.data()),
unit_direction.name);

auto normal = unit_direction(0 /* time independent */, pos);
Eigen::Matrix<double, 2, 2> 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;
}

Expand All @@ -147,39 +166,31 @@ Eigen::Matrix<double, 2, 2> CoordinateSystem::transformation<2>(
t.col(1) = Eigen::Map<Eigen::Vector2d>(
(*_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<double, 3, 3> getTransformationFromSingleBase3D(
Parameter<double> const& unit_direction, SpatialPosition const& pos)
{
auto e2 = unit_direction(0 /* time independent */, pos);
#ifndef NDEBUG
if (std::abs(Eigen::Map<Eigen::Vector3d>(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<double, 3, 3> t;
auto e2 = t.col(2);
e2 = Eigen::Map<Eigen::Vector3d const>(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)
{
Expand All @@ -197,26 +208,12 @@ Eigen::Matrix<double, 3, 3> getTransformationFromSingleBase3D(

e1.normalize();

auto eigen_mapped_e2 = Eigen::Map<Vector3>(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<double, 3, 3> 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;
}
Expand Down Expand Up @@ -245,12 +242,7 @@ Eigen::Matrix<double, 3, 3> CoordinateSystem::transformation<3>(
t.col(2) = Eigen::Map<Eigen::Vector3d>(
(*_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;
}

Expand Down Expand Up @@ -288,12 +280,7 @@ Eigen::Matrix<double, 3, 3> CoordinateSystem::transformation_3d(
Eigen::Matrix<double, 3, 3> t = Eigen::Matrix<double, 3, 3>::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;
}

Expand Down
5 changes: 2 additions & 3 deletions Tests/MaterialLib/SolidModels/TestLinearElasticTransverseIsotropic.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -241,9 +241,8 @@ TEST_F(LinearElasticTransverseIsotropic, test_agaist_ElasticOrthotropic)
}
{ // 3D
ParameterLib::ConstantParameter<double> const e1{
"e1", {-0.8191520442889918, 0.0, -0.573576436351046}};
ParameterLib::ConstantParameter<double> const e2{"e2",
{0.0, -1.0, 0.0}};
"e1", {0.8191520442889918, 0.0, 0.573576436351046}};
ParameterLib::ConstantParameter<double> const e2{"e2", {0.0, 1.0, 0.0}};
ParameterLib::ConstantParameter<double> const e3{
"e3", {-0.573576436351046, 0.0, 0.8191520442889918}};
ParameterLib::CoordinateSystem const coordinate_system{e1, e2, e3};
Expand Down
4 changes: 2 additions & 2 deletions Tests/ParameterLib/TestCoordinateSystem.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -214,10 +214,10 @@ TEST(CoordinateSystem, test3D)
{
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> parameters;

std::vector<double> e0{-0.8191520442889918, 0.0, -0.573576436351046};
std::vector<double> e0{0.8191520442889918, 0.0, 0.573576436351046};
parameters.emplace_back(
std::make_unique<ParameterLib::ConstantParameter<double>>("e0", e0));
std::vector<double> e1{0.0, -1.0, 0.0};
std::vector<double> e1{0.0, 1.0, 0.0};
parameters.emplace_back(
std::make_unique<ParameterLib::ConstantParameter<double>>("e1", e1));
std::vector<double> e2{-0.573576436351046, 0.0, 0.8191520442889918};
Expand Down

0 comments on commit 86c4348

Please sign in to comment.