OGS
CoordinateSystem.cpp
Go to the documentation of this file.
1
11#include "CoordinateSystem.h"
12
13#include <typeinfo>
14
15#include "Parameter.h"
16
17namespace ParameterLib
18{
19static double const tolerance = 1.e-15;
20static const char* const error_info =
21 "The determinant of the coordinate system transformation matrix is '{:g}', "
22 "which is not sufficiently close to unity with the tolerance of '{:g}'. "
23 "Please adjust the accuracy of the local system bases";
24
26 Parameter<double> const& e1)
27 : _base{&e0, &e1, nullptr}
28{
29 if (typeid(_base[0]) != typeid(_base[1]))
30 {
32 "The parameter types for the basis must be equal but they are "
33 "'{:s}' and '{:s}'.",
34 typeid(_base[0]).name(),
35 typeid(_base[1]).name());
36 }
37 if (_base[0]->isTimeDependent() || _base[1]->isTimeDependent())
38 {
39 OGS_FATAL("The parameters for the basis must not be time dependent.");
40 }
41 if (_base[0]->getNumberOfGlobalComponents() != 2 ||
42 _base[1]->getNumberOfGlobalComponents() != 2)
43 {
44 OGS_FATAL("The parameters for the 2D basis must have two components.");
45 }
46}
47
49 Parameter<double> const& e1,
50 Parameter<double> const& e2)
51 : _base{&e0, &e1, &e2}
52{
53 if ((typeid(_base[0]) != typeid(_base[1])) ||
54 (typeid(_base[1]) != typeid(_base[2])) ||
55 (typeid(_base[2]) != typeid(_base[0])))
56 {
58 "The parameter types for the basis must be equal but they are "
59 "'{:s}', '{:s}', and '{:s}'.",
60 typeid(_base[0]).name(),
61 typeid(_base[1]).name(),
62 typeid(_base[2]).name());
63 }
64 if (_base[0]->isTimeDependent() || _base[1]->isTimeDependent(),
65 _base[2]->isTimeDependent())
66 {
67 OGS_FATAL("The parameters for the basis must not be time dependent.");
68 }
69 if (_base[0]->getNumberOfGlobalComponents() != 3 ||
70 _base[1]->getNumberOfGlobalComponents() != 3 ||
71 _base[2]->getNumberOfGlobalComponents() != 3)
72 {
74 "The parameters for the 3D basis must have three components.");
75 }
76}
77
78template <>
79Eigen::Matrix<double, 2, 2> CoordinateSystem::transformation<2>(
80 SpatialPosition const& pos) const
81{
82 if (_base[2] != nullptr)
83 {
85 "The coordinate system is 3D but a transformation for 2D case is "
86 "requested.");
87 }
88
89 auto e0 = (*_base[0])(0 /* time independent */, pos);
90 auto e1 = (*_base[1])(0 /* time independent */, pos);
91 Eigen::Matrix<double, 2, 2> t;
92 t << e0[0], e1[0], e0[1], e1[1];
93
94#ifndef NDEBUG
95 if (std::abs(t.determinant() - 1) > tolerance)
96 {
97 OGS_FATAL(error_info, t.determinant(), tolerance);
98 }
99#endif // NDEBUG
100 return t;
101}
102
103template <>
104Eigen::Matrix<double, 3, 3> CoordinateSystem::transformation<3>(
105 SpatialPosition const& pos) const
106{
107 if (_base[2] == nullptr)
108 {
109 OGS_FATAL(
110 "The coordinate system is 2D but a transformation for 3D case is "
111 "requested.");
112 }
113
114 auto e0 = (*_base[0])(0 /* time independent */, pos);
115 auto e1 = (*_base[1])(0 /* time independent */, pos);
116 auto e2 = (*_base[2])(0 /* time independent */, pos);
117 Eigen::Matrix<double, 3, 3> t;
118 t << e0[0], e1[0], e2[0], e0[1], e1[1], e2[1], e0[2], e1[2], e2[2];
119
120#ifndef NDEBUG
121 if (std::abs(t.determinant() - 1) > tolerance)
122 {
123 OGS_FATAL(error_info, t.determinant(), tolerance);
124 }
125#endif // NDEBUG
126 return t;
127}
128
129Eigen::Matrix<double, 3, 3> CoordinateSystem::transformation_3d(
130 SpatialPosition const& pos) const
131{
132 if (_base[2] != nullptr)
133 {
134 return transformation<3>(pos);
135 }
136
137 auto e0 = (*_base[0])(0 /* time independent */, pos);
138 auto e1 = (*_base[1])(0 /* time independent */, pos);
139 Eigen::Matrix<double, 3, 3> t = Eigen::Matrix<double, 3, 3>::Identity();
140 t.template topLeftCorner<2, 2>() << e0[0], e1[0], e0[1], e1[1];
141
142#ifndef NDEBUG
143 if (std::abs(t.determinant() - 1) > tolerance)
144 {
145 OGS_FATAL(error_info, t.determinant(), tolerance);
146 }
147#endif // NDEBUG
148 return t;
149}
150
151template <int Dimension>
152Eigen::Matrix<double, Dimension, Dimension> CoordinateSystem::rotateTensor(
153 std::vector<double> const& values, SpatialPosition const& pos) const
154{
155 assert(values.size() == Dimension * Dimension ||
156 "Input vector has wrong dimension; expected 4 or 9 entries.");
157 auto const tensor =
158 Eigen::Map<Eigen::Matrix<double, Dimension, Dimension> const>(
159 values.data(), Dimension, Dimension);
160 auto const R = transformation<Dimension>(pos);
161 return R * tensor * R.transpose();
162}
163
164template <int Dimension>
165Eigen::Matrix<double, Dimension, Dimension>
166CoordinateSystem::rotateDiagonalTensor(std::vector<double> const& values,
167 SpatialPosition const& pos) const
168{
169 assert(values.size() == Dimension ||
170 "Input vector has wrong dimension; expected 2 or 3 entries.");
171 auto const tensor = Eigen::Map<Eigen::Matrix<double, Dimension, 1> const>(
172 values.data(), Dimension, 1);
173 auto const R = transformation<Dimension>(pos);
174 return R * tensor.asDiagonal() * R.transpose();
175}
176
177template Eigen::Matrix<double, 2, 2> CoordinateSystem::rotateTensor<2>(
178 std::vector<double> const& values, SpatialPosition const& pos) const;
179template Eigen::Matrix<double, 3, 3> CoordinateSystem::rotateTensor<3>(
180 std::vector<double> const& values, SpatialPosition const& pos) const;
181template Eigen::Matrix<double, 2, 2> CoordinateSystem::rotateDiagonalTensor<2>(
182 std::vector<double> const& values, SpatialPosition const& pos) const;
183template Eigen::Matrix<double, 3, 3> CoordinateSystem::rotateDiagonalTensor<3>(
184 std::vector<double> const& values, SpatialPosition const& pos) const;
185} // namespace ParameterLib
#define OGS_FATAL(...)
Definition: Error.h:26
static const double t
static double const tolerance
static const char *const error_info
std::array< Parameter< double > const *, 3 > _base
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
CoordinateSystem(Parameter< double > const &e0, Parameter< double > const &e1)
Eigen::Matrix< double, Dimension, Dimension > rotateDiagonalTensor(std::vector< double > const &values, SpatialPosition const &pos) const