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 <algorithm>
16#include <iterator>
17#include <limits>
18#include <typeinfo>
19
20#include "Parameter.h"
21
22namespace ParameterLib
23{
24static double const tolerance = std::numeric_limits<double>::epsilon();
25#ifndef NDEBUG
26static constexpr char error_info[] =
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";
30
31static constexpr char normalization_error_info[] =
32 "The direction vector given by parameter {:s} for local_coordinate_system "
33 "is not normalized to unit length";
34#endif // NDEBUG
35
37 : _base{nullptr, nullptr, &unit_direction}, _has_implicit_base(true)
38{
39 if (_base[2]->isTimeDependent())
40 {
42 "The unit_normal parameter named {} must not be time dependent.",
43 unit_direction.name);
44 }
45}
46
48 Parameter<double> const& e1)
49 : _base{&e0, &e1, nullptr}, _has_implicit_base(false)
50{
51 if (typeid(_base[0]) != typeid(_base[1]))
52 {
54 "The parameter types for the basis must be equal but they are "
55 "'{:s}' and '{:s}'.",
56 typeid(_base[0]).name(),
57 typeid(_base[1]).name());
58 }
59 if (_base[0]->isTimeDependent() || _base[1]->isTimeDependent())
60 {
61 OGS_FATAL("The parameters for the basis must not be time dependent.");
62 }
63 if (_base[0]->getNumberOfGlobalComponents() != 2 ||
64 _base[1]->getNumberOfGlobalComponents() != 2)
65 {
66 OGS_FATAL("The parameters for the 2D basis must have two components.");
67 }
68}
69
71 Parameter<double> const& e1,
72 Parameter<double> const& e2)
73 : _base{&e0, &e1, &e2}, _has_implicit_base(false)
74{
75 if ((typeid(_base[0]) != typeid(_base[1])) ||
76 (typeid(_base[1]) != typeid(_base[2])) ||
77 (typeid(_base[2]) != typeid(_base[0])))
78 {
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());
85 }
86 if (_base[0]->isTimeDependent() || _base[1]->isTimeDependent(),
87 _base[2]->isTimeDependent())
88 {
89 OGS_FATAL("The parameters for the basis must not be time dependent.");
90 }
91 if (_base[0]->getNumberOfGlobalComponents() != 3 ||
92 _base[1]->getNumberOfGlobalComponents() != 3 ||
93 _base[2]->getNumberOfGlobalComponents() != 3)
94 {
96 "The parameters for the 3D basis must have three components.");
97 }
98}
99
100Eigen::Matrix<double, 2, 2> getTransformationFromSingleBase2D(
101 Parameter<double> const& unit_direction, SpatialPosition const& pos)
102{
103 Eigen::Matrix<double, 2, 2> t;
104#ifndef NDEBUG
105 if (std::abs(Eigen::Map<Eigen::Vector2d>(
106 unit_direction(0 /* time independent */, pos).data())
107 .norm() -
108 1) > tolerance)
109 {
110 OGS_FATAL(normalization_error_info, unit_direction.name);
111 }
112#endif // NDEBUG
113
114 auto normal = unit_direction(0 /* time independent */, pos);
115 // base 0: ( normal[1], -normal[0])^T
116 // base 1: ( normal[0], normal[1])^T
117 t << normal[1], normal[0], -normal[0], normal[1];
118
119#ifndef NDEBUG
120 if (std::abs(t.determinant() - 1) > tolerance)
121 {
122 OGS_FATAL(error_info, t.determinant(), tolerance);
123 }
124#endif // NDEBUG
125 return t;
126}
127
128template <>
129Eigen::Matrix<double, 2, 2> CoordinateSystem::transformation<2>(
130 SpatialPosition const& pos) const
131{
133 {
135 }
136
137 if (_base[2] != nullptr)
138 {
139 OGS_FATAL(
140 "The coordinate system is 3D but a transformation for 2D case is "
141 "requested.");
142 }
143
144 Eigen::Matrix<double, 2, 2> t;
145 t.col(0) = Eigen::Map<Eigen::Vector2d>(
146 (*_base[0])(0 /* time independent */, pos).data());
147 t.col(1) = Eigen::Map<Eigen::Vector2d>(
148 (*_base[1])(0 /* time independent */, pos).data());
149
150#ifndef NDEBUG
151 if (std::abs(t.determinant() - 1) > tolerance)
152 {
153 OGS_FATAL(error_info, t.determinant(), tolerance);
154 }
155#endif // NDEBUG
156 return t;
157}
158
159Eigen::Matrix<double, 3, 3> getTransformationFromSingleBase3D(
160 Parameter<double> const& unit_direction, SpatialPosition const& pos)
161{
162 auto e2 = unit_direction(0 /* time independent */, pos);
163#ifndef NDEBUG
164 if (std::abs(Eigen::Map<Eigen::Vector3d>(e2.data()).norm() - 1.0) >
165 tolerance)
166 {
167 OGS_FATAL(normalization_error_info, unit_direction.name);
168 }
169#endif
170
171 // Find the id of the first non-zero component of e2:
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);
176 // Get other two component ids:
177 const auto id_a = (id + 1) % 3;
178 const auto id_b = (id + 2) % 3;
179
180 // Compute basis vector e1 orthogonal to e2
181 using Vector3 = Eigen::Vector3d;
182 Vector3 e1(0.0, 0.0, 0.0);
183
184 if (std::abs(e2[id_a]) < tolerance)
185 {
186 e1[id_a] = 1.0;
187 }
188 else if (std::abs(e2[id_b]) < tolerance)
189 {
190 e1[id_b] = 1.0;
191 }
192 else
193 {
194 e1[id_a] = 1.0;
195 e1[id_b] = -e2[id_a] / e2[id_b];
196 }
197
198 e1.normalize();
199
200 auto eigen_mapped_e2 = Eigen::Map<Vector3>(e2.data());
201
202 auto e0 = e1.cross(eigen_mapped_e2);
203 // |e0| = |e1 x e2| = |e1||e2|sin(theta) with theta the angle between e1 and
204 // e2. Since |e1| = |e2| = 1.0, and theta = pi/2, we have |e0|=1. Therefore
205 // e0 is normalized by nature.
206
207 Eigen::Matrix<double, 3, 3> t;
208 t.col(0) = e0;
209 t.col(1) = e1;
210 t.col(2) = eigen_mapped_e2;
211
212#ifndef NDEBUG
213
214 if (std::abs(t.determinant() - 1) > tolerance)
215 {
216 OGS_FATAL(error_info, t.determinant(), tolerance);
217 }
218
219#endif // NDEBUG
220
221 return t;
222}
223
224template <>
225Eigen::Matrix<double, 3, 3> CoordinateSystem::transformation<3>(
226 SpatialPosition const& pos) const
227{
229 {
231 }
232
233 if (_base[2] == nullptr)
234 {
235 OGS_FATAL(
236 "The coordinate system is 2D but a transformation for 3D case is "
237 "requested.");
238 }
239
240 Eigen::Matrix<double, 3, 3> t;
241 t.col(0) = Eigen::Map<Eigen::Vector3d>(
242 (*_base[0])(0 /* time independent */, pos).data());
243 t.col(1) = Eigen::Map<Eigen::Vector3d>(
244 (*_base[1])(0 /* time independent */, pos).data());
245 t.col(2) = Eigen::Map<Eigen::Vector3d>(
246 (*_base[2])(0 /* time independent */, pos).data());
247
248#ifndef NDEBUG
249 if (std::abs(t.determinant() - 1) > tolerance)
250 {
251 OGS_FATAL(error_info, t.determinant(), tolerance);
252 }
253#endif // NDEBUG
254 return t;
255}
256
258 SpatialPosition const& pos) const
259{
260 // If base 2, which stores the unit direction vector, has three components:
261 if ((*_base[2])(0 /* time independent */, pos).size() == 3)
262 {
264 }
265
266 // If base 2, which stores the unit direction vector, has two components
267 Eigen::Matrix<double, 3, 3> t = Eigen::Matrix<double, 3, 3>::Identity();
268 t.template topLeftCorner<2, 2>() = transformation<2>(pos);
269
270 return t;
271}
272
273Eigen::Matrix<double, 3, 3> CoordinateSystem::transformation_3d(
274 SpatialPosition const& pos) const
275{
277 {
279 }
280
281 if (_base[2] != nullptr)
282 {
283 return transformation<3>(pos);
284 }
285
286 auto e0 = (*_base[0])(0 /* time independent */, pos);
287 auto e1 = (*_base[1])(0 /* time independent */, 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];
290
291#ifndef NDEBUG
292 if (std::abs(t.determinant() - 1) > tolerance)
293 {
294 OGS_FATAL(error_info, t.determinant(), tolerance);
295 }
296#endif // NDEBUG
297 return t;
298}
299
300template <int Dimension>
301Eigen::Matrix<double, Dimension, Dimension> CoordinateSystem::rotateTensor(
302 std::vector<double> const& values, SpatialPosition const& pos) const
303{
304 assert(values.size() == Dimension * Dimension ||
305 "Input vector has wrong dimension; expected 4 or 9 entries.");
306 auto const tensor =
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();
311}
312
313template <int Dimension>
314Eigen::Matrix<double, Dimension, Dimension>
315CoordinateSystem::rotateDiagonalTensor(std::vector<double> const& values,
316 SpatialPosition const& pos) const
317{
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();
324}
325
326template Eigen::Matrix<double, 2, 2> CoordinateSystem::rotateTensor<2>(
327 std::vector<double> const& values, SpatialPosition const& pos) const;
328template Eigen::Matrix<double, 3, 3> CoordinateSystem::rotateTensor<3>(
329 std::vector<double> const& values, SpatialPosition const& pos) const;
330template Eigen::Matrix<double, 2, 2> CoordinateSystem::rotateDiagonalTensor<2>(
331 std::vector<double> const& values, SpatialPosition const& pos) const;
332template Eigen::Matrix<double, 3, 3> CoordinateSystem::rotateDiagonalTensor<3>(
333 std::vector<double> const& values, SpatialPosition const& pos) const;
334} // 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 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)
std::string const name
Definition: Parameter.h:73