OGS
KelvinVector.cpp
Go to the documentation of this file.
1
11#include "KelvinVector.h"
12
13#include "BaseLib/Error.h"
14
15namespace MathLib
16{
17namespace KelvinVector
18{
19template <>
20double Invariants<6>::determinant(Eigen::Matrix<double, 6, 1> const& v)
21{
22 return v(0) * v(1) * v(2) + v(3) * v(4) * v(5) / std::sqrt(2.) -
23 v(3) * v(3) * v(2) / 2. - v(4) * v(4) * v(0) / 2. -
24 v(5) * v(5) * v(1) / 2.;
25}
26
27template <>
28double Invariants<4>::determinant(Eigen::Matrix<double, 4, 1> const& v)
29{
30 return v(2) * (v(0) * v(1) - v(3) * v(3) / 2.);
31}
32
33template <>
34Eigen::Matrix<double, 4, 1, Eigen::ColMajor, 4, 1> inverse(
35 Eigen::Matrix<double, 4, 1, Eigen::ColMajor, 4, 1> const& v)
36{
37 assert(Invariants<4>::determinant(v) != 0);
38
39 Eigen::Matrix<double, 4, 1, Eigen::ColMajor, 4, 1> inv;
40 inv(0) = v(1) * v(2);
41 inv(1) = v(0) * v(2);
42 inv(2) = v(0) * v(1) - v(3) * v(3) / 2.;
43 inv(3) = -v(3) * v(2);
44 return inv / Invariants<4>::determinant(v);
45}
46
47template <>
48Eigen::Matrix<double, 6, 1, Eigen::ColMajor, 6, 1> inverse(
49 Eigen::Matrix<double, 6, 1, Eigen::ColMajor, 6, 1> const& v)
50{
51 assert(Invariants<6>::determinant(v) != 0);
52
53 Eigen::Matrix<double, 6, 1, Eigen::ColMajor, 6, 1> inv;
54 inv(0) = v(1) * v(2) - v(4) * v(4) / 2.;
55 inv(1) = v(0) * v(2) - v(5) * v(5) / 2.;
56 inv(2) = v(0) * v(1) - v(3) * v(3) / 2.;
57 inv(3) = v(4) * v(5) / std::sqrt(2.) - v(3) * v(2);
58 inv(4) = v(3) * v(5) / std::sqrt(2.) - v(4) * v(0);
59 inv(5) = v(4) * v(3) / std::sqrt(2.) - v(1) * v(5);
60 return inv / Invariants<6>::determinant(v);
61}
62
63template <>
64Eigen::Matrix<double, 3, 3> kelvinVectorToTensor(
65 Eigen::Matrix<double, 4, 1, Eigen::ColMajor, 4, 1> const& v)
66{
67 Eigen::Matrix<double, 3, 3> m;
68 m << v[0], v[3] / std::sqrt(2.), 0, v[3] / std::sqrt(2.), v[1], 0, 0, 0,
69 v[2];
70 return m;
71}
72
73template <>
74Eigen::Matrix<double, 3, 3> kelvinVectorToTensor(
75 Eigen::Matrix<double, 6, 1, Eigen::ColMajor, 6, 1> const& v)
76{
77 Eigen::Matrix<double, 3, 3> m;
78 m << v[0], v[3] / std::sqrt(2.), v[5] / std::sqrt(2.), v[3] / std::sqrt(2.),
79 v[1], v[4] / std::sqrt(2.), v[5] / std::sqrt(2.), v[4] / std::sqrt(2.),
80 v[2];
81 return m;
82}
83
84template <>
85Eigen::Matrix<double, 3, 3> kelvinVectorToTensor(Eigen::Matrix<double,
86 Eigen::Dynamic,
87 1,
88 Eigen::ColMajor,
89 Eigen::Dynamic,
90 1> const& v)
91{
92 if (v.size() == 4)
93 {
94 Eigen::Matrix<double, 4, 1, Eigen::ColMajor, 4, 1> v4;
95 v4 << v[0], v[1], v[2], v[3];
96 return kelvinVectorToTensor(v4);
97 }
98 if (v.size() == 6)
99 {
100 Eigen::Matrix<double, 6, 1, Eigen::ColMajor, 6, 1> v6;
101 v6 << v[0], v[1], v[2], v[3], v[4], v[5];
102 return kelvinVectorToTensor(v6);
103 }
104 OGS_FATAL(
105 "Conversion of dynamic Kelvin vector of size {:d} to a tensor is not "
106 "possible. Kelvin vector must be of size 4 or 6.",
107 v.size());
108}
109
110template <>
111KelvinVectorType<2> tensorToKelvin<2>(Eigen::Matrix<double, 3, 3> const& m)
112{
113 assert(std::abs(m(0, 1) - m(1, 0)) <
114 std::numeric_limits<double>::epsilon());
115 assert(m(0, 2) == 0);
116 assert(m(1, 2) == 0);
117 assert(m(2, 0) == 0);
118 assert(m(2, 1) == 0);
119
120 KelvinVectorType<2> v;
121 v << m(0, 0), m(1, 1), m(2, 2), m(0, 1) * std::sqrt(2.);
122 return v;
123}
124
125template <>
126KelvinVectorType<3> tensorToKelvin<3>(Eigen::Matrix<double, 3, 3> const& m)
127{
128 assert(std::abs(m(0, 1) - m(1, 0)) <
129 std::numeric_limits<double>::epsilon());
130 assert(std::abs(m(1, 2) - m(2, 1)) <
131 std::numeric_limits<double>::epsilon());
132 assert(std::abs(m(0, 2) - m(2, 0)) <
133 std::numeric_limits<double>::epsilon());
134
135 KelvinVectorType<3> v;
136 v << m(0, 0), m(1, 1), m(2, 2), m(0, 1) * std::sqrt(2.),
137 m(1, 2) * std::sqrt(2.), m(0, 2) * std::sqrt(2.);
138 return v;
139}
140
141template <>
142Eigen::Matrix<double, 4, 1> kelvinVectorToSymmetricTensor(
143 Eigen::Matrix<double, 4, 1, Eigen::ColMajor, 4, 1> const& v)
144{
145 Eigen::Matrix<double, 4, 1> m;
146 m << v[0], v[1], v[2], v[3] / std::sqrt(2.);
147 return m;
148}
149
150template <>
151Eigen::Matrix<double, 6, 1> kelvinVectorToSymmetricTensor(
152 Eigen::Matrix<double, 6, 1, Eigen::ColMajor, 6, 1> const& v)
153{
154 Eigen::Matrix<double, 6, 1> m;
155 m << v[0], v[1], v[2], v[3] / std::sqrt(2.), v[4] / std::sqrt(2.),
156 v[5] / std::sqrt(2.);
157 return m;
158}
159
160template <>
161Eigen::Matrix<double, Eigen::Dynamic, 1, Eigen::ColMajor, Eigen::Dynamic, 1>
162kelvinVectorToSymmetricTensor(Eigen::Matrix<double,
163 Eigen::Dynamic,
164 1,
165 Eigen::ColMajor,
166 Eigen::Dynamic,
167 1> const& v)
168{
169 if (v.size() == 4)
170 {
171 return kelvinVectorToSymmetricTensor<4>(v);
172 }
173 if (v.size() == 6)
174 {
175 return kelvinVectorToSymmetricTensor<6>(v);
176 }
177 OGS_FATAL(
178 "Kelvin vector to tensor conversion expected an input vector of size 4 "
179 "or 6, but a vector of size {:d} was given.",
180 v.size());
181}
182
183template <>
184Eigen::Matrix<double, 2, 4> liftVectorToKelvin<2>(
185 Eigen::Matrix<double, 2, 1> const& v)
186{
187 Eigen::Matrix<double, 2, 4> m;
188 m << v[0], 0, 0, v[1] / std::sqrt(2.), 0, v[1], 0, v[0] / std::sqrt(2.);
189 return m;
190}
191
192template <>
193Eigen::Matrix<double, 3, 6> liftVectorToKelvin<3>(
194 Eigen::Matrix<double, 3, 1> const& v)
195{
196 Eigen::Matrix<double, 3, 6> m;
197 m << v[0], 0, 0, v[1] / std::sqrt(2.), 0, v[2] / std::sqrt(2.), 0, v[1], 0,
198 v[0] / std::sqrt(2.), v[2] / std::sqrt(2.), 0, 0, 0, v[2], 0,
199 v[1] / std::sqrt(2.), v[0] / std::sqrt(2.);
200 return m;
201}
202
203template <>
204Eigen::Matrix<double, 2, 1> reduceKelvinToVector<2>(
205 Eigen::Matrix<double, 2, 4> const& m)
206{
207 assert(m(0, 1) == 0);
208 assert(m(0, 2) == 0);
209 assert(m(1, 0) == 0);
210 assert(m(1, 2) == 0);
211 Eigen::Matrix<double, 2, 1> v;
212 v << m(0, 0), m(1, 1);
213 return v;
214}
215
216template <>
217Eigen::Matrix<double, 3, 1> reduceKelvinToVector<3>(
218 Eigen::Matrix<double, 3, 6> const& m)
219{
220 assert(m(0, 1) == 0);
221 assert(m(0, 2) == 0);
222 assert(m(0, 4) == 0);
223 assert(m(1, 0) == 0);
224 assert(m(1, 2) == 0);
225 assert(m(1, 5) == 0);
226 assert(m(2, 0) == 0);
227 assert(m(2, 1) == 0);
228 assert(m(2, 3) == 0);
229 assert(std::abs(m(0, 3) - m(2, 4)) <
230 std::numeric_limits<double>::epsilon());
231 assert(std::abs(m(0, 5) - m(1, 4)) <
232 std::numeric_limits<double>::epsilon());
233 assert(std::abs(m(1, 3) - m(2, 5)) <
234 std::numeric_limits<double>::epsilon());
235 Eigen::Matrix<double, 3, 1> v;
236 v << m(0, 0), m(1, 1), m(2, 2);
237 return v;
238}
239
240template <>
241KelvinMatrixType<2> fourthOrderRotationMatrix<2>(
242 Eigen::Matrix<double, 2, 2, Eigen::ColMajor, 2, 2> const& transformation)
243{
244 // 1-based index access for convenience.
245 auto Q = [&](int const i, int const j)
246 { return transformation(i - 1, j - 1); };
247
249 R << Q(1, 1) * Q(1, 1), Q(1, 2) * Q(1, 2), 0,
250 std::sqrt(2) * Q(1, 1) * Q(1, 2), Q(2, 1) * Q(2, 1), Q(2, 2) * Q(2, 2),
251 0, std::sqrt(2) * Q(2, 1) * Q(2, 2), 0, 0, 1, 0,
252 std::sqrt(2) * Q(1, 1) * Q(2, 1), std::sqrt(2) * Q(1, 2) * Q(2, 2), 0,
253 Q(1, 1) * Q(2, 2) + Q(1, 2) * Q(2, 1);
254 return R;
255}
256
257template <>
258KelvinMatrixType<3> fourthOrderRotationMatrix<3>(
259 Eigen::Matrix<double, 3, 3, Eigen::ColMajor, 3, 3> const& transformation)
260{
261 // 1-based index access for convenience.
262 auto Q = [&](int const i, int const j)
263 { return transformation(i - 1, j - 1); };
264
266 R << Q(1, 1) * Q(1, 1), Q(1, 2) * Q(1, 2), Q(1, 3) * Q(1, 3),
267 std::sqrt(2) * Q(1, 1) * Q(1, 2), std::sqrt(2) * Q(1, 2) * Q(1, 3),
268 std::sqrt(2) * Q(1, 1) * Q(1, 3), Q(2, 1) * Q(2, 1), Q(2, 2) * Q(2, 2),
269 Q(2, 3) * Q(2, 3), std::sqrt(2) * Q(2, 1) * Q(2, 2),
270 std::sqrt(2) * Q(2, 2) * Q(2, 3), std::sqrt(2) * Q(2, 1) * Q(2, 3),
271 Q(3, 1) * Q(3, 1), Q(3, 2) * Q(3, 2), Q(3, 3) * Q(3, 3),
272 std::sqrt(2) * Q(3, 1) * Q(3, 2), std::sqrt(2) * Q(3, 2) * Q(3, 3),
273 std::sqrt(2) * Q(3, 1) * Q(3, 3), std::sqrt(2) * Q(1, 1) * Q(2, 1),
274 std::sqrt(2) * Q(1, 2) * Q(2, 2), std::sqrt(2) * Q(1, 3) * Q(2, 3),
275 Q(1, 1) * Q(2, 2) + Q(1, 2) * Q(2, 1),
276 Q(1, 2) * Q(2, 3) + Q(1, 3) * Q(2, 2),
277 Q(1, 1) * Q(2, 3) + Q(1, 3) * Q(2, 1), std::sqrt(2) * Q(2, 1) * Q(3, 1),
278 std::sqrt(2) * Q(2, 2) * Q(3, 2), std::sqrt(2) * Q(2, 3) * Q(3, 3),
279 Q(2, 1) * Q(3, 2) + Q(2, 2) * Q(3, 1),
280 Q(2, 2) * Q(3, 3) + Q(2, 3) * Q(3, 2),
281 Q(2, 1) * Q(3, 3) + Q(2, 3) * Q(3, 1), std::sqrt(2) * Q(1, 1) * Q(3, 1),
282 std::sqrt(2) * Q(1, 2) * Q(3, 2), std::sqrt(2) * Q(1, 3) * Q(3, 3),
283 Q(1, 1) * Q(3, 2) + Q(1, 2) * Q(3, 1),
284 Q(1, 2) * Q(3, 3) + Q(1, 3) * Q(3, 2),
285 Q(1, 1) * Q(3, 3) + Q(1, 3) * Q(3, 1);
286 return R;
287}
288} // namespace KelvinVector
289} // namespace MathLib
#define OGS_FATAL(...)
Definition Error.h:26
KelvinMatrixType< 3 > fourthOrderRotationMatrix< 3 >(Eigen::Matrix< double, 3, 3, Eigen::ColMajor, 3, 3 > const &transformation)
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
KelvinVectorType< 2 > tensorToKelvin< 2 >(Eigen::Matrix< double, 3, 3 > const &m)
Eigen::Matrix< double, 2, 4 > liftVectorToKelvin< 2 >(Eigen::Matrix< double, 2, 1 > const &v)
Eigen::Matrix< double, 2, 1 > reduceKelvinToVector< 2 >(Eigen::Matrix< double, 2, 4 > const &m)
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Eigen::Matrix< double, 3, 6 > liftVectorToKelvin< 3 >(Eigen::Matrix< double, 3, 1 > const &v)
KelvinMatrixType< 2 > fourthOrderRotationMatrix< 2 >(Eigen::Matrix< double, 2, 2, Eigen::ColMajor, 2, 2 > const &transformation)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
KelvinVectorType< 3 > tensorToKelvin< 3 >(Eigen::Matrix< double, 3, 3 > const &m)
Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > inverse(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Eigen::Matrix< double, 3, 1 > reduceKelvinToVector< 3 >(Eigen::Matrix< double, 3, 6 > const &m)
static const double v
KV::KelvinVectorType< DisplacementDim > KelvinVector
Definition Base.h:27
static double determinant(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Determinant of a matrix in Kelvin vector representation.