OpenLB 1.7
Loading...
Searching...
No Matches
geometricOperations.h
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2021 Julius Jessberger
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22 */
23
28#ifndef GEOMETRIC_OPERATIONS_H
29#define GEOMETRIC_OPERATIONS_H
30
31#include "core/vector.h"
32#include "dimensionConverter.h"
33#include "matrix.h"
34
35namespace olb {
36
37namespace util {
38
39// angle conversions
40template <typename T, unsigned D>
41decltype(Vector<decltype(util::sqrt(T())), D>())
43{
44 constexpr BaseType<decltype(util::sqrt(T()))> conversionFactor = M_PI / 180.;
45 return conversionFactor * angle;
46}
47
48template <typename T>
49decltype(util::sqrt(T())) degreeToRadian(T angle)
50{
51 return degreeToRadian(Vector<T, 1>(angle))[0];
52}
53
54template <typename T, unsigned D>
57{
58 constexpr BaseType<decltype(util::sqrt(T()))> conversionFactor = 180. / M_PI;
59 return conversionFactor * angle;
60}
61
62template <typename T>
63decltype(util::sqrt(T())) radianToDegree(T angle)
64{
65 return radianToDegree(Vector<T, 1>(angle))[0];
66}
67
68template <typename T, unsigned D>
71{
73
74 if constexpr (D == 2) {
75 T const cos = util::cos(angle[0]);
76 T const sin = util::sin(angle[0]);
77
78 // row 1
79 rotationMatrix[0] = cos;
80 rotationMatrix[1] = -sin;
81 // row 2
82 rotationMatrix[2] = sin;
83 rotationMatrix[3] = cos;
84 }
85 else {
86 T const cos[3] = {util::cos(angle[0]), util::cos(angle[1]),
87 util::cos(angle[2])};
88 T const sin[3] = {util::sin(angle[0]), util::sin(angle[1]),
89 util::sin(angle[2])};
90
91 // |x0| / 0 1 2 \ |x1|
92 // |y0| | 3 4 5 | = |x1|
93 // |z0| \ 6 7 8 / |x1|
94
95 // row 1
96 rotationMatrix[0] = cos[1] * cos[2];
97 rotationMatrix[1] = sin[0] * sin[1] * cos[2] - cos[0] * sin[2];
98 rotationMatrix[2] = cos[0] * sin[1] * cos[2] + sin[0] * sin[2];
99 // row 2
100 rotationMatrix[3] = cos[1] * sin[2];
101 rotationMatrix[4] = sin[0] * sin[1] * sin[2] + cos[0] * cos[2];
102 rotationMatrix[5] = cos[0] * sin[1] * sin[2] - sin[0] * cos[2];
103 // row 3
104 rotationMatrix[6] = -sin[1];
105 rotationMatrix[7] = sin[0] * cos[1];
106 rotationMatrix[8] = cos[0] * cos[1];
107 }
108
109 return rotationMatrix;
110}
111
112template <typename T, unsigned D>
114 const Vector<T, utilities::dimensions::convert<D>::matrix>& rotationMatrix)
115{
116 if constexpr (D == 2) {
117 return Vector<T, 4>(rotationMatrix[0], rotationMatrix[2], rotationMatrix[1],
118 rotationMatrix[3]);
119 }
120 else {
121
122 //Individual Entries ( source: https://doi.org/10.1016/j.compfluid.2018.02.027 )
123 // |x0| / 0 3 6 \ |x1|
124 // |y0| | 1 4 7 | = |x1|
125 // |z0| \ 2 5 8 / |x1|
126
127 return Vector<T, 9>(rotationMatrix[0], rotationMatrix[3], rotationMatrix[6],
128 rotationMatrix[1], rotationMatrix[4], rotationMatrix[7],
129 rotationMatrix[2], rotationMatrix[5],
130 rotationMatrix[8]);
131 }
132 __builtin_unreachable();
133}
134
135template <typename T, unsigned D>
139{
140 return invertRotationMatrix<T, D>(calculateRotationMatrix<T, D>(angle));
141}
142
144template <typename T, unsigned D,
145 bool OUTPUT_USES_ROTATION_CENTER_AS_ORIGIN = false>
147 const Vector<T, D>& input,
148 const Vector<T, utilities::dimensions::convert<D>::matrix>& rotationMatrix,
149 const Vector<T, D>& rotationCenter = Vector<T, D>(0.))
150{
151 const Vector<T, D> dist = input - rotationCenter;
152 Vector<T, D> rotated;
153
154 if constexpr (D == 2) {
155 rotated =
156 Vector<T, 2>(dist[0] * rotationMatrix[0] + dist[1] * rotationMatrix[1],
157 dist[0] * rotationMatrix[2] + dist[1] * rotationMatrix[3]);
158 }
159 else {
160 rotated =
161 Vector<T, 3>(rotationMatrix[0] * dist[0] + rotationMatrix[1] * dist[1] +
162 rotationMatrix[2] * dist[2],
163 rotationMatrix[3] * dist[0] + rotationMatrix[4] * dist[1] +
164 rotationMatrix[5] * dist[2],
165 rotationMatrix[6] * dist[0] + rotationMatrix[7] * dist[1] +
166 rotationMatrix[8] * dist[2]);
167 }
168
169 if constexpr (!OUTPUT_USES_ROTATION_CENTER_AS_ORIGIN) {
170 return rotationCenter + rotated;
171 }
172 else {
173 return rotated;
174 }
175
176 __builtin_unreachable();
177}
178
180template <typename T, unsigned D,
181 bool OUTPUT_USES_ROTATION_CENTER_AS_ORIGIN = false>
183 const Vector<T, D>& input,
184 const Vector<T, utilities::dimensions::convert<D>::matrix>& rotationMatrix,
185 const Vector<T, D>& rotationCenter = Vector<T, D>(0.))
186{
188 invertRotationMatrix<T, D>(rotationMatrix);
189 return executeRotation(input, invRotationMatrix, rotationCenter);
190}
191
193template <typename T, unsigned D>
195 const Vector<T, D>& rotationCenter, const Vector<T, D>& velocity,
197 angularVelocity,
198 const Vector<T, D>& position)
199{
200 if constexpr (D == 2) {
201 // two dimensions: u = U + w x r = (Ux, Uy, 0) + (0,0,w) x (X,Y,0) = (Ux, Uy, 0) + (-w*Y, w*X, 0)
202 return Vector<T, 2>(
203 velocity[0] - angularVelocity[0] * (position[1] - rotationCenter[1]),
204 velocity[1] + angularVelocity[0] * (position[0] - rotationCenter[0]));
205 }
206 else {
207 // three dimensions: u = U + w x r = (Ux, Uy, Uz) + (wx,wy,wz) x (X,Y,Z) = (Ux, Uy, Uz) + (wy*Z-wz*Y, wz*X-wx*Z, wx*Y-wy*X)
208 return velocity +
209 crossProduct3D(angularVelocity, position - rotationCenter);
210 }
211 __builtin_unreachable();
212}
213
215template <typename T>
217 const Vector<T, 9>& rotationMatrix)
218{
219 // I' = R(angle) * I * R(angle)^T
220 // TODO: The inertia tensor is symmetric - make use of it
221 T data[3][3] = {
222 {mofi[0], T {}, T {}},
223 { T {}, mofi[1], T {}},
224 { T {}, T {}, mofi[2]}
225 };
226 Matrix<T, 3, 3> inertiaTensor(data);
227 const Matrix<T, 3, 3> _rotationMatrix(rotationMatrix);
228 inertiaTensor =
229 _rotationMatrix * inertiaTensor * (_rotationMatrix.transpose());
230
231 return inertiaTensor;
232}
233
234} // namespace util
235
236} // namespace olb
237
238#endif
#define M_PI
Matrix with a defined number of ROWS and columns (COLS)
Definition matrix.h:31
constexpr Matrix< T, COLS, ROWS > transpose() const
Definition matrix.h:125
Plain old scalar vector.
Definition vector.h:47
decltype(Vector< decltype(util::sqrt(T())), D >()) radianToDegree(const Vector< T, D > &angle)
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
Vector< T, D > invertRotation(const Vector< T, D > &input, const Vector< T, utilities::dimensions::convert< D >::matrix > &rotationMatrix, const Vector< T, D > &rotationCenter=Vector< T, D >(0.))
Rotates the input around the rotationCenter with a given rotationMatrix in the opposite direction.
constexpr Matrix< T, 3, 3 > rotateMofi(const Vector< T, 3 > &mofi, const Vector< T, 9 > &rotationMatrix)
Rotate moment of inertia (mofi)
ADf< T, DIM > sin(const ADf< T, DIM > &a)
Definition aDiff.h:569
Vector< T, utilities::dimensions::convert< D >::matrix > calculateRotationMatrix(const Vector< T, utilities::dimensions::convert< D >::rotation > &angle)
decltype(Vector< decltype(util::sqrt(T())), D >()) degreeToRadian(const Vector< T, D > &angle)
Vector< T, D > executeRotation(const Vector< T, D > &input, const Vector< T, utilities::dimensions::convert< D >::matrix > &rotationMatrix, const Vector< T, D > &rotationCenter=Vector< T, D >(0.))
Rotates the input around the rotationCenter with a given rotationMatrix.
constexpr Vector< T, D > calculateLocalVelocity(const Vector< T, D > &rotationCenter, const Vector< T, D > &velocity, const Vector< T, utilities::dimensions::convert< D >::rotation > &angularVelocity, const Vector< T, D > &position)
Calculate local velocity.
Vector< T, utilities::dimensions::convert< D >::matrix > invertRotationMatrix(const Vector< T, utilities::dimensions::convert< D >::matrix > &rotationMatrix)
Vector< T, utilities::dimensions::convert< D >::matrix > calculateInverseRotationMatrix(const Vector< T, utilities::dimensions::convert< D >::rotation > &angle)
ADf< T, DIM > cos(const ADf< T, DIM > &a)
Definition aDiff.h:578
Top level namespace for all of OpenLB.
typename util::BaseTypeHelper< T >::type BaseType
Definition baseType.h:59
constexpr Vector< T, 3 > crossProduct3D(const ScalarVector< T, 3, IMPL > &a, const ScalarVector< T, 3, IMPL_ > &b) any_platform
Definition vector.h:224
Converts dimensions by deriving from given cartesian dimension D.
efficient implementation of a vector class