transform.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2024 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 InNamespace
25  Foam
26 
27 Description
28  3D tensor transformation operations.
29 
30 \*---------------------------------------------------------------------------*/
31 
32 #ifndef transform_H
33 #define transform_H
34 
35 #include "tensor.H"
36 #include "mathematicalConstants.H"
37 
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42 
43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
44 
45 //- Rotational transformation tensor from unit vector n1 to n2
47 (
48  const vector& n1,
49  const vector& n2
50 )
51 {
52  const scalar c = n1 & n2;
53  const vector n3 = n1 ^ n2;
54  const scalar magSqrN3 = magSqr(n3);
55 
56  // n1 and n2 define a plane n3
57  if (magSqrN3 > small)
58  {
59  // Return rotational transformation tensor in the n3-plane using
60  // Rodrigues' rotation formula
61  return
62  c*I
63  + (n2*n1 - n1*n2)
64  + (1 - c)*sqr(n3)/magSqrN3;
65  }
66  // n1 and n2 are contradirectional
67  else if (c < 0)
68  {
69  // Return rotational transformation tensor in a plane with arbitrary
70  // normal perpendicular to n1 and n2 using a subset of Rodrigues'
71  // rotation formula
72  const vector n4 = perpendicular(n1);
73  const scalar magSqrN4 = magSqr(n4);
74  return - I + 2*sqr(n4)/magSqrN4;
75  }
76  // n1 and n2 are codirectional
77  else
78  {
79  // Return null transformation tensor
80  return I;
81  }
82 }
83 
84 
85 //- Rotational transformation tensor about the x-axis by omega radians
86 // The rotation is defined in a right-handed coordinate system
87 // i.e. clockwise with respect to the axis from -ve to +ve
88 // (looking along the axis).
89 inline tensor Rx(const scalar& omega)
90 {
91  const scalar s = sin(omega);
92  const scalar c = cos(omega);
93  return tensor
94  (
95  1, 0, 0,
96  0, c, -s,
97  0, s, c
98  );
99 }
100 
101 
102 //- Rotational transformation tensor about the y-axis by omega radians
103 // The rotation is defined in a right-handed coordinate system
104 // i.e. clockwise with respect to the axis from -ve to +ve
105 // (looking along the axis).
106 inline tensor Ry(const scalar& omega)
107 {
108  const scalar s = sin(omega);
109  const scalar c = cos(omega);
110  return tensor
111  (
112  c, 0, s,
113  0, 1, 0,
114  -s, 0, c
115  );
116 }
117 
118 
119 //- Rotational transformation tensor about the z-axis by omega radians
120 // The rotation is defined in a right-handed coordinate system
121 // i.e. clockwise with respect to the axis from -ve to +ve
122 // (looking along the axis).
123 inline tensor Rz(const scalar& omega)
124 {
125  const scalar s = sin(omega);
126  const scalar c = cos(omega);
127  return tensor
128  (
129  c, -s, 0,
130  s, c, 0,
131  0, 0, 1
132  );
133 }
134 
135 
136 //- Rotational transformation tensor about axis a by omega radians
137 // The rotation is defined in a right-handed coordinate system
138 // i.e. clockwise with respect to the axis from -ve to +ve
139 // (looking along the axis).
140 inline tensor Ra(const vector& a, const scalar omega)
141 {
142  const scalar s = sin(omega);
143  const scalar c = cos(omega);
144 
145  return tensor
146  (
147  sqr(a.x())*(1 - c) + c,
148  a.y()*a.x()*(1 - c) - a.z()*s,
149  a.x()*a.z()*(1 - c) + a.y()*s,
150 
151  a.x()*a.y()*(1 - c) + a.z()*s,
152  sqr(a.y())*(1 - c) + c,
153  a.y()*a.z()*(1 - c) - a.x()*s,
154 
155  a.x()*a.z()*(1 - c) - a.y()*s,
156  a.y()*a.z()*(1 - c) + a.x()*s,
157  sqr(a.z())*(1 - c) + c
158  );
159 }
160 
161 
162 inline bool transform(const tensor&, const bool i)
163 {
164  return i;
165 }
166 
167 
168 inline label transform(const tensor&, const label i)
169 {
170  return i;
171 }
172 
173 
174 inline scalar transform(const tensor&, const scalar s)
175 {
176  return s;
177 }
178 
179 
180 template<class Cmpt>
181 inline Vector<Cmpt> transform(const tensor& tt, const Vector<Cmpt>& v)
182 {
183  return tt & v;
184 }
185 
186 
187 template<class Cmpt>
188 inline Tensor<Cmpt> transform(const tensor& tt, const Tensor<Cmpt>& t)
189 {
190  return Tensor<Cmpt>
191  (
192  (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.xx()
193  + (tt.xx()*t.xy() + tt.xy()*t.yy() + tt.xz()*t.zy())*tt.xy()
194  + (tt.xx()*t.xz() + tt.xy()*t.yz() + tt.xz()*t.zz())*tt.xz(),
195 
196  (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.yx()
197  + (tt.xx()*t.xy() + tt.xy()*t.yy() + tt.xz()*t.zy())*tt.yy()
198  + (tt.xx()*t.xz() + tt.xy()*t.yz() + tt.xz()*t.zz())*tt.yz(),
199 
200  (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.zx()
201  + (tt.xx()*t.xy() + tt.xy()*t.yy() + tt.xz()*t.zy())*tt.zy()
202  + (tt.xx()*t.xz() + tt.xy()*t.yz() + tt.xz()*t.zz())*tt.zz(),
203 
204  (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.xx()
205  + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.xy()
206  + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.xz(),
207 
208  (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.yx()
209  + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.yy()
210  + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.yz(),
211 
212  (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.zx()
213  + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.zy()
214  + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.zz(),
215 
216  (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.xx()
217  + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.xy()
218  + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.xz(),
219 
220  (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.yx()
221  + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.yy()
222  + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.yz(),
223 
224  (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.zx()
225  + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.zy()
226  + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.zz()
227  );
228 }
229 
230 
231 template<class Cmpt>
233 (
234  const tensor& tt,
235  const SphericalTensor<Cmpt>& st
236 )
237 {
238  return st;
239 }
240 
241 
242 template<class Cmpt>
243 inline SymmTensor<Cmpt> transform(const tensor& tt, const SymmTensor<Cmpt>& st)
244 {
245  return SymmTensor<Cmpt>
246  (
247  (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.xx()
248  + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.xy()
249  + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.xz(),
250 
251  (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.yx()
252  + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.yy()
253  + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.yz(),
254 
255  (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.zx()
256  + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.zy()
257  + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.zz(),
258 
259  (tt.yx()*st.xx() + tt.yy()*st.xy() + tt.yz()*st.xz())*tt.yx()
260  + (tt.yx()*st.xy() + tt.yy()*st.yy() + tt.yz()*st.yz())*tt.yy()
261  + (tt.yx()*st.xz() + tt.yy()*st.yz() + tt.yz()*st.zz())*tt.yz(),
262 
263  (tt.yx()*st.xx() + tt.yy()*st.xy() + tt.yz()*st.xz())*tt.zx()
264  + (tt.yx()*st.xy() + tt.yy()*st.yy() + tt.yz()*st.yz())*tt.zy()
265  + (tt.yx()*st.xz() + tt.yy()*st.yz() + tt.yz()*st.zz())*tt.zz(),
266 
267  (tt.zx()*st.xx() + tt.zy()*st.xy() + tt.zz()*st.xz())*tt.zx()
268  + (tt.zx()*st.xy() + tt.zy()*st.yy() + tt.zz()*st.yz())*tt.zy()
269  + (tt.zx()*st.xz() + tt.zy()*st.yz() + tt.zz()*st.zz())*tt.zz()
270  );
271 }
272 
273 
274 template<class Type1, class Type2>
275 inline Type1 transformMask(const Type2& t)
276 {
277  return t;
278 }
279 
280 
281 template<>
283 {
284  return sph(t);
285 }
286 
287 
288 template<>
290 {
291  return symm(t);
292 }
293 
294 
295 //- Estimate angle of vec in coordinate system (e0, e1, e0^e1).
296 // Is guaranteed to return increasing number but is not correct
297 // angle. Used for sorting angles. All input vectors need to be normalised.
298 //
299 // Calculates scalar which increases with angle going from e0 to vec in
300 // the coordinate system e0, e1, e0^e1
301 //
302 // Jumps from 2*pi -> 0 at -small so hopefully parallel vectors with small
303 // rounding errors should still get the same quadrant.
304 //
305 inline scalar pseudoAngle
306 (
307  const vector& e0,
308  const vector& e1,
309  const vector& vec
310 )
311 {
312  scalar cos = vec & e0;
313  scalar sin = vec & e1;
314 
315  if (sin < -small)
316  {
317  return (3.0 + cos)*constant::mathematical::piByTwo;
318  }
319  else
320  {
321  return (1.0 - cos)*constant::mathematical::piByTwo;
322  }
323 }
324 
325 
326 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
327 
328 } // End namespace Foam
329 
330 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
331 
332 #endif
333 
334 // ************************************************************************* //
Templated 3D SphericalTensor derived from VectorSpace adding construction from 1 component,...
Templated 3D symmetric tensor derived from VectorSpace adding construction from 6 components,...
Definition: SymmTensor.H:56
const Cmpt & xx() const
Definition: SymmTensorI.H:87
const Cmpt & yz() const
Definition: SymmTensorI.H:111
const Cmpt & xz() const
Definition: SymmTensorI.H:99
const Cmpt & zz() const
Definition: SymmTensorI.H:117
const Cmpt & xy() const
Definition: SymmTensorI.H:93
const Cmpt & yy() const
Definition: SymmTensorI.H:105
const Cmpt & xx() const
Definition: TensorI.H:163
const Cmpt & yx() const
Definition: TensorI.H:184
const Cmpt & yz() const
Definition: TensorI.H:198
const Cmpt & xz() const
Definition: TensorI.H:177
const Cmpt & zz() const
Definition: TensorI.H:219
const Cmpt & xy() const
Definition: TensorI.H:170
const Cmpt & zx() const
Definition: TensorI.H:205
const Cmpt & zy() const
Definition: TensorI.H:212
const Cmpt & yy() const
Definition: TensorI.H:191
const Cmpt & z() const
Definition: VectorI.H:87
const Cmpt & y() const
Definition: VectorI.H:81
const Cmpt & x() const
Definition: VectorI.H:75
Convenience class to handle the input of constant rotational speed. Reads an omega entry with default...
Definition: omega.H:54
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const scalar piByTwo(0.5 *pi)
const dimensionedScalar c
Speed of light in a vacuum.
Namespace for OpenFOAM.
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from unit vector n1 to n2.
Definition: transform.H:47
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
tensor Ra(const vector &a, const scalar omega)
Rotational transformation tensor about axis a by omega radians.
Definition: transform.H:140
tensor Rx(const scalar &omega)
Rotational transformation tensor about the x-axis by omega radians.
Definition: transform.H:89
dimensionedScalar sin(const dimensionedScalar &ds)
SphericalTensor< Cmpt > sph(const DiagTensor< Cmpt > &dt)
Return the spherical part of a diagonal tensor.
Definition: DiagTensorI.H:353
static const Identity< scalar > I
Definition: Identity.H:93
tensor Ry(const scalar &omega)
Rotational transformation tensor about the y-axis by omega radians.
Definition: transform.H:106
Type1 transformMask(const Type2 &t)
Definition: transform.H:275
sphericalTensor transformMask< sphericalTensor >(const symmTensor &st)
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:504
scalar pseudoAngle(const vector &e0, const vector &e1, const vector &vec)
Estimate angle of vec in coordinate system (e0, e1, e0^e1).
Definition: transform.H:306
tensor Rz(const scalar &omega)
Rotational transformation tensor about the z-axis by omega radians.
Definition: transform.H:123
dimensionSet perpendicular(const dimensionSet &)
Definition: dimensionSet.C:516
symmTensor transformMask< symmTensor >(const symmTensor &st)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
dimensionedScalar cos(const dimensionedScalar &ds)