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-2020 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 inline tensor Rx(const scalar& omega)
87 {
88  const scalar s = sin(omega);
89  const scalar c = cos(omega);
90  return tensor
91  (
92  1, 0, 0,
93  0, c, s,
94  0, -s, c
95  );
96 }
97 
98 
99 //- Rotational transformation tensor about the y-axis by omega radians
100 inline tensor Ry(const scalar& omega)
101 {
102  const scalar s = sin(omega);
103  const scalar c = cos(omega);
104  return tensor
105  (
106  c, 0, -s,
107  0, 1, 0,
108  s, 0, c
109  );
110 }
111 
112 
113 //- Rotational transformation tensor about the z-axis by omega radians
114 inline tensor Rz(const scalar& omega)
115 {
116  const scalar s = sin(omega);
117  const scalar c = cos(omega);
118  return tensor
119  (
120  c, s, 0,
121  -s, c, 0,
122  0, 0, 1
123  );
124 }
125 
126 
127 //- Rotational transformation tensor about axis a by omega radians
128 inline tensor Ra(const vector& a, const scalar omega)
129 {
130  const scalar s = sin(omega);
131  const scalar c = cos(omega);
132 
133  return tensor
134  (
135  sqr(a.x())*(1 - c) + c,
136  a.y()*a.x()*(1 - c) + a.z()*s,
137  a.x()*a.z()*(1 - c) - a.y()*s,
138 
139  a.x()*a.y()*(1 - c) - a.z()*s,
140  sqr(a.y())*(1 - c) + c,
141  a.y()*a.z()*(1 - c) + a.x()*s,
142 
143  a.x()*a.z()*(1 - c) + a.y()*s,
144  a.y()*a.z()*(1 - c) - a.x()*s,
145  sqr(a.z())*(1 - c) + c
146  );
147 }
148 
149 
150 inline label transform(const tensor&, const bool i)
151 {
152  return i;
153 }
154 
155 
156 inline label transform(const tensor&, const label i)
157 {
158  return i;
159 }
160 
161 
162 inline scalar transform(const tensor&, const scalar s)
163 {
164  return s;
165 }
166 
167 
168 template<class Cmpt>
169 inline Vector<Cmpt> transform(const tensor& tt, const Vector<Cmpt>& v)
170 {
171  return tt & v;
172 }
173 
174 
175 template<class Cmpt>
176 inline Tensor<Cmpt> transform(const tensor& tt, const Tensor<Cmpt>& t)
177 {
178  return Tensor<Cmpt>
179  (
180  (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.xx()
181  + (tt.xx()*t.xy() + tt.xy()*t.yy() + tt.xz()*t.zy())*tt.xy()
182  + (tt.xx()*t.xz() + tt.xy()*t.yz() + tt.xz()*t.zz())*tt.xz(),
183 
184  (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.yx()
185  + (tt.xx()*t.xy() + tt.xy()*t.yy() + tt.xz()*t.zy())*tt.yy()
186  + (tt.xx()*t.xz() + tt.xy()*t.yz() + tt.xz()*t.zz())*tt.yz(),
187 
188  (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.zx()
189  + (tt.xx()*t.xy() + tt.xy()*t.yy() + tt.xz()*t.zy())*tt.zy()
190  + (tt.xx()*t.xz() + tt.xy()*t.yz() + tt.xz()*t.zz())*tt.zz(),
191 
192  (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.xx()
193  + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.xy()
194  + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.xz(),
195 
196  (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.yx()
197  + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.yy()
198  + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.yz(),
199 
200  (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.zx()
201  + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.zy()
202  + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.zz(),
203 
204  (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.xx()
205  + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.xy()
206  + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.xz(),
207 
208  (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.yx()
209  + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.yy()
210  + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.yz(),
211 
212  (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.zx()
213  + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.zy()
214  + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.zz()
215  );
216 }
217 
218 
219 template<class Cmpt>
221 (
222  const tensor& tt,
223  const SphericalTensor<Cmpt>& st
224 )
225 {
226  return st;
227 }
228 
229 
230 template<class Cmpt>
231 inline SymmTensor<Cmpt> transform(const tensor& tt, const SymmTensor<Cmpt>& st)
232 {
233  return SymmTensor<Cmpt>
234  (
235  (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.xx()
236  + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.xy()
237  + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.xz(),
238 
239  (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.yx()
240  + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.yy()
241  + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.yz(),
242 
243  (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.zx()
244  + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.zy()
245  + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.zz(),
246 
247  (tt.yx()*st.xx() + tt.yy()*st.xy() + tt.yz()*st.xz())*tt.yx()
248  + (tt.yx()*st.xy() + tt.yy()*st.yy() + tt.yz()*st.yz())*tt.yy()
249  + (tt.yx()*st.xz() + tt.yy()*st.yz() + tt.yz()*st.zz())*tt.yz(),
250 
251  (tt.yx()*st.xx() + tt.yy()*st.xy() + tt.yz()*st.xz())*tt.zx()
252  + (tt.yx()*st.xy() + tt.yy()*st.yy() + tt.yz()*st.yz())*tt.zy()
253  + (tt.yx()*st.xz() + tt.yy()*st.yz() + tt.yz()*st.zz())*tt.zz(),
254 
255  (tt.zx()*st.xx() + tt.zy()*st.xy() + tt.zz()*st.xz())*tt.zx()
256  + (tt.zx()*st.xy() + tt.zy()*st.yy() + tt.zz()*st.yz())*tt.zy()
257  + (tt.zx()*st.xz() + tt.zy()*st.yz() + tt.zz()*st.zz())*tt.zz()
258  );
259 }
260 
261 
262 template<class Type1, class Type2>
263 inline Type1 transformMask(const Type2& t)
264 {
265  return t;
266 }
267 
268 
269 template<>
271 {
272  return sph(t);
273 }
274 
275 
276 template<>
278 {
279  return symm(t);
280 }
281 
282 
283 //- Estimate angle of vec in coordinate system (e0, e1, e0^e1).
284 // Is guaranteed to return increasing number but is not correct
285 // angle. Used for sorting angles. All input vectors need to be normalized.
286 //
287 // Calculates scalar which increases with angle going from e0 to vec in
288 // the coordinate system e0, e1, e0^e1
289 //
290 // Jumps from 2*pi -> 0 at -small so hopefully parallel vectors with small
291 // rounding errors should still get the same quadrant.
292 //
293 inline scalar pseudoAngle
294 (
295  const vector& e0,
296  const vector& e1,
297  const vector& vec
298 )
299 {
300  scalar cos = vec & e0;
301  scalar sin = vec & e1;
302 
303  if (sin < -small)
304  {
305  return (3.0 + cos)*constant::mathematical::piByTwo;
306  }
307  else
308  {
309  return (1.0 - cos)*constant::mathematical::piByTwo;
310  }
311 }
312 
313 
314 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
315 
316 } // End namespace Foam
317 
318 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
319 
320 #endif
321 
322 // ************************************************************************* //
const Cmpt & yz() const
Definition: SymmTensorI.H:111
const Cmpt & yx() const
Definition: TensorI.H:184
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
const Cmpt & zy() const
Definition: TensorI.H:212
Templated 3D symmetric tensor derived from VectorSpace adding construction from 6 components...
Definition: SymmTensor.H:53
const Cmpt & xz() const
Definition: TensorI.H:177
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:294
const Cmpt & xx() const
Definition: SymmTensorI.H:87
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const Cmpt & yy() const
Definition: TensorI.H:191
const Cmpt & xx() const
Definition: TensorI.H:163
const Cmpt & yz() const
Definition: TensorI.H:198
const Cmpt & z() const
Definition: VectorI.H:87
sphericalTensor transformMask< sphericalTensor >(const symmTensor &st)
const dimensionedScalar & c
Speed of light in a vacuum.
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from unit vector n1 to n2.
Definition: transform.H:47
SphericalTensor< Cmpt > sph(const DiagTensor< Cmpt > &dt)
Return the spherical part of a diagonal tensor.
Definition: DiagTensorI.H:288
const Cmpt & yy() const
Definition: SymmTensorI.H:105
const Cmpt & y() const
Definition: VectorI.H:81
Templated 3D SphericalTensor derived from VectorSpace adding construction from 1 component, element access using th ii() member function and the inner-product (dot-product) and outer-product operators.
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.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const Cmpt & xz() const
Definition: SymmTensorI.H:99
dimensionedScalar cos(const dimensionedScalar &ds)
symmTensor transformMask< symmTensor >(const symmTensor &st)
const Cmpt & zz() const
Definition: SymmTensorI.H:117
static const Identity< scalar > I
Definition: Identity.H:93
const Cmpt & x() const
Definition: VectorI.H:75
dimensioned< scalar > magSqr(const dimensioned< Type > &)
dimensionedScalar sin(const dimensionedScalar &ds)
tensor Ry(const scalar &omega)
Rotational transformation tensor about the y-axis by omega radians.
Definition: transform.H:100
Vector< Cmpt > perpendicular(const Vector< Cmpt > &v)
Definition: VectorI.H:166
const Cmpt & zx() const
Definition: TensorI.H:205
Type1 transformMask(const Type2 &t)
Definition: transform.H:263
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
tensor Rz(const scalar &omega)
Rotational transformation tensor about the z-axis by omega radians.
Definition: transform.H:114
tensor Ra(const vector &a, const scalar omega)
Rotational transformation tensor about axis a by omega radians.
Definition: transform.H:128
const Cmpt & xy() const
Definition: SymmTensorI.H:93
const scalar piByTwo(0.5 *pi)
const Cmpt & xy() const
Definition: TensorI.H:170
const Cmpt & zz() const
Definition: TensorI.H:219
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
tensor Rx(const scalar &omega)
Rotational transformation tensor about the x-axis by omega radians.
Definition: transform.H:86
Namespace for OpenFOAM.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:477