OGS
CoordinateSystem.cpp
Go to the documentation of this file.
1
11#include "CoordinateSystem.h"
12
13#include <Eigen/Dense>
14#include <Eigen/Geometry>
15#include <cmath>
16#include <limits>
17#include <typeinfo>
18
20#include "Parameter.h"
21
22namespace ParameterLib
23{
24static double const tolerance = std::numeric_limits<double>::epsilon();
25
26template <int Dim>
28 Eigen::Matrix<double, Dim, Dim, Eigen::ColMajor, Dim, Dim> const& t)
29{
30 if (std::abs(t.determinant() - 1) > tolerance)
31 {
33 "The determinant of the coordinate system transformation matrix is "
34 "'{:g}', which is not sufficiently close to unity with the "
35 "tolerance of '{:g}'. Please adjust the accuracy of the local "
36 "system bases",
37 t.determinant(), tolerance);
38 }
39 if (((t * t.transpose() - Eigen::Matrix<double, Dim, Dim>::Identity())
40 .array() > 2 * tolerance)
41 .any())
42 {
44 "The transformation is not orthogonal because the difference "
45 "A*A^T - I is:\n{}\nand at least one component deviates from zero "
46 "by more then '{:g}'.",
47 (t * t.transpose() - Eigen::Matrix<double, Dim, Dim>::Identity())
48 .eval(),
49 2 * tolerance);
50 }
51}
52
53template <typename Derived>
54static void checkNormalization(Eigen::MatrixBase<Derived> const& vec,
55 std::string_view const parmeter_name)
56{
57 if (std::abs(vec.squaredNorm() - 1.0) > tolerance)
58 {
60 "The direction vector given by parameter {:s} for the "
61 "local_coordinate_system is not normalized to unit length. Vector "
62 "norm is {:g} and |v|^2-1 = {:g}.",
63 parmeter_name, vec.norm(), vec.squaredNorm() - 1.0);
64 }
65}
66
68 : _base{nullptr, nullptr, &unit_direction}, _has_implicit_base(true)
69{
70 if (_base[2]->isTimeDependent())
71 {
73 "The unit_normal parameter named {} must not be time dependent.",
74 unit_direction.name);
75 }
76}
77
79 Parameter<double> const& e1)
80 : _base{&e0, &e1, nullptr}, _has_implicit_base(false)
81{
82 if (typeid(_base[0]) != typeid(_base[1]))
83 {
85 "The parameter types for the basis must be equal but they are "
86 "'{:s}' and '{:s}'.",
87 typeid(_base[0]).name(),
88 typeid(_base[1]).name());
89 }
90 if (_base[0]->isTimeDependent() || _base[1]->isTimeDependent())
91 {
92 OGS_FATAL("The parameters for the basis must not be time dependent.");
93 }
94 if (_base[0]->getNumberOfGlobalComponents() != 2 ||
95 _base[1]->getNumberOfGlobalComponents() != 2)
96 {
97 OGS_FATAL("The parameters for the 2D basis must have two components.");
98 }
99}
100
102 Parameter<double> const& e1,
103 Parameter<double> const& e2)
104 : _base{&e0, &e1, &e2}, _has_implicit_base(false)
105{
106 if ((typeid(_base[0]) != typeid(_base[1])) ||
107 (typeid(_base[1]) != typeid(_base[2])) ||
108 (typeid(_base[2]) != typeid(_base[0])))
109 {
110 OGS_FATAL(
111 "The parameter types for the basis must be equal but they are "
112 "'{:s}', '{:s}', and '{:s}'.",
113 typeid(_base[0]).name(),
114 typeid(_base[1]).name(),
115 typeid(_base[2]).name());
116 }
117 if (_base[0]->isTimeDependent() || _base[1]->isTimeDependent(),
118 _base[2]->isTimeDependent())
119 {
120 OGS_FATAL("The parameters for the basis must not be time dependent.");
121 }
122 if (_base[0]->getNumberOfGlobalComponents() != 3 ||
123 _base[1]->getNumberOfGlobalComponents() != 3 ||
124 _base[2]->getNumberOfGlobalComponents() != 3)
125 {
126 OGS_FATAL(
127 "The parameters for the 3D basis must have three components.");
128 }
129}
130
131Eigen::Matrix<double, 2, 2> getTransformationFromSingleBase2D(
132 Parameter<double> const& unit_direction, SpatialPosition const& pos)
133{
134 auto const& normal = unit_direction(0 /* time independent */, pos);
135 checkNormalization(Eigen::Map<Eigen::Vector2d const>(normal.data()),
136 unit_direction.name);
137
138 Eigen::Matrix<double, 2, 2> t;
139 // base 0: ( normal[1], -normal[0])^T
140 // base 1: ( normal[0], normal[1])^T
141 t << normal[1], normal[0], -normal[0], normal[1];
142
144 return t;
145}
146
147template <>
148Eigen::Matrix<double, 2, 2> CoordinateSystem::transformation<2>(
149 SpatialPosition const& pos) const
150{
152 {
154 }
155
156 if (_base[2] != nullptr)
157 {
158 OGS_FATAL(
159 "The coordinate system is 3D but a transformation for 2D case is "
160 "requested.");
161 }
162
163 Eigen::Matrix<double, 2, 2> t;
164 t.col(0) = Eigen::Map<Eigen::Vector2d>(
165 (*_base[0])(0 /* time independent */, pos).data());
166 t.col(1) = Eigen::Map<Eigen::Vector2d>(
167 (*_base[1])(0 /* time independent */, pos).data());
168
170 return t;
171}
172
173Eigen::Matrix<double, 3, 3> getTransformationFromSingleBase3D(
174 Parameter<double> const& unit_direction, SpatialPosition const& pos)
175{
176 auto const& normal = unit_direction(0 /* time independent */, pos);
177
178 Eigen::Matrix<double, 3, 3> t;
179 auto e2 = t.col(2);
180 e2 = Eigen::Map<Eigen::Vector3d const>(normal.data());
181 checkNormalization(e2, unit_direction.name);
182
183 // Find the id of the first non-zero component of e2:
184 int id;
185 e2.cwiseAbs().maxCoeff(&id);
186
187 // Get other two component ids:
188 const auto id_a = (id + 1) % 3;
189 const auto id_b = (id + 2) % 3;
190
191 // Compute basis vector e1 orthogonal to e2
192 auto e1 = t.col(1);
193 e1 = Eigen::Vector3d::Zero();
194
195 if (std::abs(e2[id_a]) < tolerance)
196 {
197 e1[id_a] = 1.0;
198 }
199 else if (std::abs(e2[id_b]) < tolerance)
200 {
201 e1[id_b] = 1.0;
202 }
203 else
204 {
205 e1[id_a] = 1.0;
206 e1[id_b] = -e2[id_a] / e2[id_b];
207 }
208
209 e1.normalize();
210
211 // |e0| = |e1 x e2| = |e1||e2|sin(theta) with theta the angle between e1 and
212 // e2. Since |e1| = |e2| = 1.0, and theta = pi/2, we have |e0|=1. Therefore
213 // e0 is normalized by nature.
214 t.col(0) = e1.cross(e2);
215
217
218 return t;
219}
220
221template <>
222Eigen::Matrix<double, 3, 3> CoordinateSystem::transformation<3>(
223 SpatialPosition const& pos) const
224{
226 {
228 }
229
230 if (_base[2] == nullptr)
231 {
232 OGS_FATAL(
233 "The coordinate system is 2D but a transformation for 3D case is "
234 "requested.");
235 }
236
237 Eigen::Matrix<double, 3, 3> t;
238 t.col(0) = Eigen::Map<Eigen::Vector3d>(
239 (*_base[0])(0 /* time independent */, pos).data());
240 t.col(1) = Eigen::Map<Eigen::Vector3d>(
241 (*_base[1])(0 /* time independent */, pos).data());
242 t.col(2) = Eigen::Map<Eigen::Vector3d>(
243 (*_base[2])(0 /* time independent */, pos).data());
244
246 return t;
247}
248
250 SpatialPosition const& pos) const
251{
252 // If base 2, which stores the unit direction vector, has three components:
253 if ((*_base[2])(0 /* time independent */, pos).size() == 3)
254 {
256 }
257
258 // If base 2, which stores the unit direction vector, has two components
259 Eigen::Matrix<double, 3, 3> t = Eigen::Matrix<double, 3, 3>::Identity();
260 t.template topLeftCorner<2, 2>() = transformation<2>(pos);
261
262 return t;
263}
264
265Eigen::Matrix<double, 3, 3> CoordinateSystem::transformation_3d(
266 SpatialPosition const& pos) const
267{
269 {
271 }
272
273 if (_base[2] != nullptr)
274 {
275 return transformation<3>(pos);
276 }
277
278 auto e0 = (*_base[0])(0 /* time independent */, pos);
279 auto e1 = (*_base[1])(0 /* time independent */, pos);
280 Eigen::Matrix<double, 3, 3> t = Eigen::Matrix<double, 3, 3>::Identity();
281 t.template topLeftCorner<2, 2>() << e0[0], e1[0], e0[1], e1[1];
282
284 return t;
285}
286
287template <int Dimension>
288Eigen::Matrix<double, Dimension, Dimension> CoordinateSystem::rotateTensor(
289 std::vector<double> const& values, SpatialPosition const& pos) const
290{
291 assert(values.size() == Dimension * Dimension ||
292 "Input vector has wrong dimension; expected 4 or 9 entries.");
293 auto const tensor =
294 Eigen::Map<Eigen::Matrix<double, Dimension, Dimension> const>(
295 values.data(), Dimension, Dimension);
296 return rotateTensor<Dimension>(tensor, pos);
297}
298
299template <int Dimension>
300Eigen::Matrix<double, Dimension, Dimension>
301CoordinateSystem::rotateDiagonalTensor(std::vector<double> const& values,
302 SpatialPosition const& pos) const
303{
304 assert(values.size() == Dimension ||
305 "Input vector has wrong dimension; expected 2 or 3 entries.");
306 auto const tensor = Eigen::Map<Eigen::Matrix<double, Dimension, 1> const>(
307 values.data(), Dimension, 1);
308 auto const R = transformation<Dimension>(pos);
309 return R * tensor.asDiagonal() * R.transpose();
310}
311
312template Eigen::Matrix<double, 2, 2> CoordinateSystem::rotateTensor<2>(
313 std::vector<double> const& values, SpatialPosition const& pos) const;
314template Eigen::Matrix<double, 3, 3> CoordinateSystem::rotateTensor<3>(
315 std::vector<double> const& values, SpatialPosition const& pos) const;
316template Eigen::Matrix<double, 2, 2> CoordinateSystem::rotateDiagonalTensor<2>(
317 std::vector<double> const& values, SpatialPosition const& pos) const;
318template Eigen::Matrix<double, 3, 3> CoordinateSystem::rotateDiagonalTensor<3>(
319 std::vector<double> const& values, SpatialPosition const& pos) const;
320} // namespace ParameterLib
#define OGS_FATAL(...)
Definition Error.h:26
Eigen::Matrix< double, 2, 2 > getTransformationFromSingleBase2D(Parameter< double > const &unit_direction, SpatialPosition const &pos)
static double const tolerance
static void checkNormalization(Eigen::MatrixBase< Derived > const &vec, std::string_view const parmeter_name)
static void checkTransformationIsSON(Eigen::Matrix< double, Dim, Dim, Eigen::ColMajor, Dim, Dim > const &t)
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)
std::string const name
Definition Parameter.h:73