14#include <Eigen/Geometry>
24static double const tolerance = std::numeric_limits<double>::epsilon();
27 "The determinant of the coordinate system transformation matrix is '{:g}', "
28 "which is not sufficiently close to unity with the tolerance of '{:g}'. "
29 "Please adjust the accuracy of the local system bases";
32 "The direction vector given by parameter {:s} for local_coordinate_system "
33 "is not normalized to unit length";
37 : _base{nullptr, nullptr, &unit_direction}, _has_implicit_base(true)
39 if (
_base[2]->isTimeDependent())
42 "The unit_normal parameter named {} must not be time dependent.",
49 : _base{&e0, &e1, nullptr}, _has_implicit_base(false)
54 "The parameter types for the basis must be equal but they are "
56 typeid(
_base[0]).name(),
57 typeid(
_base[1]).name());
59 if (
_base[0]->isTimeDependent() ||
_base[1]->isTimeDependent())
61 OGS_FATAL(
"The parameters for the basis must not be time dependent.");
63 if (
_base[0]->getNumberOfGlobalComponents() != 2 ||
64 _base[1]->getNumberOfGlobalComponents() != 2)
66 OGS_FATAL(
"The parameters for the 2D basis must have two components.");
73 : _base{&e0, &e1, &e2}, _has_implicit_base(false)
80 "The parameter types for the basis must be equal but they are "
81 "'{:s}', '{:s}', and '{:s}'.",
82 typeid(
_base[0]).name(),
83 typeid(
_base[1]).name(),
84 typeid(
_base[2]).name());
86 if (
_base[0]->isTimeDependent() ||
_base[1]->isTimeDependent(),
87 _base[2]->isTimeDependent())
89 OGS_FATAL(
"The parameters for the basis must not be time dependent.");
91 if (
_base[0]->getNumberOfGlobalComponents() != 3 ||
92 _base[1]->getNumberOfGlobalComponents() != 3 ||
93 _base[2]->getNumberOfGlobalComponents() != 3)
96 "The parameters for the 3D basis must have three components.");
103 Eigen::Matrix<double, 2, 2> t;
105 if (std::abs(Eigen::Map<Eigen::Vector2d>(
106 unit_direction(0 , pos).data())
114 auto normal = unit_direction(0 , pos);
117 t << normal[1], normal[0], -normal[0], normal[1];
120 if (std::abs(t.determinant() - 1) >
tolerance)
129Eigen::Matrix<double, 2, 2> CoordinateSystem::transformation<2>(
137 if (
_base[2] !=
nullptr)
140 "The coordinate system is 3D but a transformation for 2D case is "
144 Eigen::Matrix<double, 2, 2> t;
145 t.col(0) = Eigen::Map<Eigen::Vector2d>(
146 (*
_base[0])(0 , pos).data());
147 t.col(1) = Eigen::Map<Eigen::Vector2d>(
148 (*
_base[1])(0 , pos).data());
151 if (std::abs(t.determinant() - 1) >
tolerance)
162 auto e2 = unit_direction(0 , pos);
164 if (std::abs(Eigen::Map<Eigen::Vector3d>(e2.data()).norm() - 1.0) >
172 auto const it = std::max_element(e2.begin(), e2.end(),
173 [](
const double a,
const double b)
174 { return std::abs(a) < std::abs(b); });
175 const auto id = std::distance(e2.begin(), it);
177 const auto id_a = (
id + 1) % 3;
178 const auto id_b = (
id + 2) % 3;
181 using Vector3 = Eigen::Vector3d;
182 Vector3 e1(0.0, 0.0, 0.0);
195 e1[id_b] = -e2[id_a] / e2[id_b];
200 auto eigen_mapped_e2 = Eigen::Map<Vector3>(e2.data());
202 auto e0 = e1.cross(eigen_mapped_e2);
207 Eigen::Matrix<double, 3, 3> t;
210 t.col(2) = eigen_mapped_e2;
214 if (std::abs(t.determinant() - 1) >
tolerance)
225Eigen::Matrix<double, 3, 3> CoordinateSystem::transformation<3>(
233 if (
_base[2] ==
nullptr)
236 "The coordinate system is 2D but a transformation for 3D case is "
240 Eigen::Matrix<double, 3, 3> t;
241 t.col(0) = Eigen::Map<Eigen::Vector3d>(
242 (*
_base[0])(0 , pos).data());
243 t.col(1) = Eigen::Map<Eigen::Vector3d>(
244 (*
_base[1])(0 , pos).data());
245 t.col(2) = Eigen::Map<Eigen::Vector3d>(
246 (*
_base[2])(0 , pos).data());
249 if (std::abs(t.determinant() - 1) >
tolerance)
261 if ((*
_base[2])(0 , pos).size() == 3)
267 Eigen::Matrix<double, 3, 3> t = Eigen::Matrix<double, 3, 3>::Identity();
268 t.template topLeftCorner<2, 2>() = transformation<2>(pos);
281 if (
_base[2] !=
nullptr)
283 return transformation<3>(pos);
286 auto e0 = (*
_base[0])(0 , pos);
287 auto e1 = (*
_base[1])(0 , pos);
288 Eigen::Matrix<double, 3, 3> t = Eigen::Matrix<double, 3, 3>::Identity();
289 t.template topLeftCorner<2, 2>() << e0[0], e1[0], e0[1], e1[1];
292 if (std::abs(t.determinant() - 1) >
tolerance)
300template <
int Dimension>
304 assert(values.size() == Dimension * Dimension ||
305 "Input vector has wrong dimension; expected 4 or 9 entries.");
307 Eigen::Map<Eigen::Matrix<double, Dimension, Dimension>
const>(
308 values.data(), Dimension, Dimension);
309 auto const R = transformation<Dimension>(pos);
310 return R * tensor * R.transpose();
313template <
int Dimension>
314Eigen::Matrix<double, Dimension, Dimension>
318 assert(values.size() == Dimension ||
319 "Input vector has wrong dimension; expected 2 or 3 entries.");
320 auto const tensor = Eigen::Map<Eigen::Matrix<double, Dimension, 1>
const>(
321 values.data(), Dimension, 1);
322 auto const R = transformation<Dimension>(pos);
323 return R * tensor.asDiagonal() * R.transpose();
326template Eigen::Matrix<double, 2, 2> CoordinateSystem::rotateTensor<2>(
328template Eigen::Matrix<double, 3, 3> CoordinateSystem::rotateTensor<3>(
330template Eigen::Matrix<double, 2, 2> CoordinateSystem::rotateDiagonalTensor<2>(
332template Eigen::Matrix<double, 3, 3> CoordinateSystem::rotateDiagonalTensor<3>(
Eigen::Matrix< double, 2, 2 > getTransformationFromSingleBase2D(Parameter< double > const &unit_direction, SpatialPosition const &pos)
static double const tolerance
static constexpr char error_info[]
static constexpr char normalization_error_info[]
Eigen::Matrix< double, 3, 3 > getTransformationFromSingleBase3D(Parameter< double > const &unit_direction, SpatialPosition const &pos)
std::array< Parameter< double > const *, 3 > _base
Eigen::Matrix< double, 3, 3 > transformationFromSingleBase_3d(SpatialPosition const &pos) const
Eigen::Matrix< double, 3, 3 > transformation_3d(SpatialPosition const &pos) const
Eigen::Matrix< double, Dimension, Dimension > rotateTensor(std::vector< double > const &values, SpatialPosition const &pos) const
Eigen::Matrix< double, Dimension, Dimension > rotateDiagonalTensor(std::vector< double > const &values, SpatialPosition const &pos) const
CoordinateSystem(Parameter< double > const &unit_direction)