transform.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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 vector n1 to n2
47 (
48  const vector& n1,
49  const vector& n2
50 )
51 {
52  const scalar s = 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
60  return
61  s*I
62  + (1 - s)*sqr(n3)/magSqrN3
63  + (n2*n1 - n1*n2);
64  }
65  // n1 and n2 are contradirectional
66  else if (s < 0)
67  {
68  // Return mirror transformation tensor
69  return I + 2*n1*n2;
70  }
71  // n1 and n2 are codirectional
72  else
73  {
74  // Return null transformation tensor
75  return I;
76  }
77 }
78 
79 
80 //- Rotational transformation tensor about the x-axis by omega radians
81 inline tensor Rx(const scalar& omega)
82 {
83  const scalar s = sin(omega);
84  const scalar c = cos(omega);
85  return tensor
86  (
87  1, 0, 0,
88  0, c, s,
89  0, -s, c
90  );
91 }
92 
93 
94 //- Rotational transformation tensor about the y-axis by omega radians
95 inline tensor Ry(const scalar& omega)
96 {
97  const scalar s = sin(omega);
98  const scalar c = cos(omega);
99  return tensor
100  (
101  c, 0, -s,
102  0, 1, 0,
103  s, 0, c
104  );
105 }
106 
107 
108 //- Rotational transformation tensor about the z-axis by omega radians
109 inline tensor Rz(const scalar& omega)
110 {
111  const scalar s = sin(omega);
112  const scalar c = cos(omega);
113  return tensor
114  (
115  c, s, 0,
116  -s, c, 0,
117  0, 0, 1
118  );
119 }
120 
121 
122 //- Rotational transformation tensor about axis a by omega radians
123 inline tensor Ra(const vector& a, const scalar omega)
124 {
125  const scalar s = sin(omega);
126  const scalar c = cos(omega);
127 
128  return tensor
129  (
130  sqr(a.x())*(1 - c) + c,
131  a.y()*a.x()*(1 - c) + a.z()*s,
132  a.x()*a.z()*(1 - c) - a.y()*s,
133 
134  a.x()*a.y()*(1 - c) - a.z()*s,
135  sqr(a.y())*(1 - c) + c,
136  a.y()*a.z()*(1 - c) + a.x()*s,
137 
138  a.x()*a.z()*(1 - c) + a.y()*s,
139  a.y()*a.z()*(1 - c) - a.x()*s,
140  sqr(a.z())*(1 - c) + c
141  );
142 }
143 
144 
145 inline label transform(const tensor&, const bool i)
146 {
147  return i;
148 }
149 
150 
151 inline label transform(const tensor&, const label i)
152 {
153  return i;
154 }
155 
156 
157 inline scalar transform(const tensor&, const scalar s)
158 {
159  return s;
160 }
161 
162 
163 template<class Cmpt>
164 inline Vector<Cmpt> transform(const tensor& tt, const Vector<Cmpt>& v)
165 {
166  return tt & v;
167 }
168 
169 
170 template<class Cmpt>
171 inline Tensor<Cmpt> transform(const tensor& tt, const Tensor<Cmpt>& t)
172 {
173  return Tensor<Cmpt>
174  (
175  (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.xx()
176  + (tt.xx()*t.xy() + tt.xy()*t.yy() + tt.xz()*t.zy())*tt.xy()
177  + (tt.xx()*t.xz() + tt.xy()*t.yz() + tt.xz()*t.zz())*tt.xz(),
178 
179  (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.yx()
180  + (tt.xx()*t.xy() + tt.xy()*t.yy() + tt.xz()*t.zy())*tt.yy()
181  + (tt.xx()*t.xz() + tt.xy()*t.yz() + tt.xz()*t.zz())*tt.yz(),
182 
183  (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.zx()
184  + (tt.xx()*t.xy() + tt.xy()*t.yy() + tt.xz()*t.zy())*tt.zy()
185  + (tt.xx()*t.xz() + tt.xy()*t.yz() + tt.xz()*t.zz())*tt.zz(),
186 
187  (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.xx()
188  + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.xy()
189  + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.xz(),
190 
191  (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.yx()
192  + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.yy()
193  + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.yz(),
194 
195  (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.zx()
196  + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.zy()
197  + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.zz(),
198 
199  (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.xx()
200  + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.xy()
201  + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.xz(),
202 
203  (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.yx()
204  + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.yy()
205  + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.yz(),
206 
207  (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.zx()
208  + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.zy()
209  + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.zz()
210  );
211 }
212 
213 
214 template<class Cmpt>
216 (
217  const tensor& tt,
218  const SphericalTensor<Cmpt>& st
219 )
220 {
221  return st;
222 }
223 
224 
225 template<class Cmpt>
226 inline SymmTensor<Cmpt> transform(const tensor& tt, const SymmTensor<Cmpt>& st)
227 {
228  return SymmTensor<Cmpt>
229  (
230  (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.xx()
231  + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.xy()
232  + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.xz(),
233 
234  (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.yx()
235  + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.yy()
236  + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.yz(),
237 
238  (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.zx()
239  + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.zy()
240  + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.zz(),
241 
242  (tt.yx()*st.xx() + tt.yy()*st.xy() + tt.yz()*st.xz())*tt.yx()
243  + (tt.yx()*st.xy() + tt.yy()*st.yy() + tt.yz()*st.yz())*tt.yy()
244  + (tt.yx()*st.xz() + tt.yy()*st.yz() + tt.yz()*st.zz())*tt.yz(),
245 
246  (tt.yx()*st.xx() + tt.yy()*st.xy() + tt.yz()*st.xz())*tt.zx()
247  + (tt.yx()*st.xy() + tt.yy()*st.yy() + tt.yz()*st.yz())*tt.zy()
248  + (tt.yx()*st.xz() + tt.yy()*st.yz() + tt.yz()*st.zz())*tt.zz(),
249 
250  (tt.zx()*st.xx() + tt.zy()*st.xy() + tt.zz()*st.xz())*tt.zx()
251  + (tt.zx()*st.xy() + tt.zy()*st.yy() + tt.zz()*st.yz())*tt.zy()
252  + (tt.zx()*st.xz() + tt.zy()*st.yz() + tt.zz()*st.zz())*tt.zz()
253  );
254 }
255 
256 
257 template<class Type1, class Type2>
258 inline Type1 transformMask(const Type2& t)
259 {
260  return t;
261 }
262 
263 
264 template<>
266 {
267  return sph(t);
268 }
269 
270 
271 template<>
273 {
274  return symm(t);
275 }
276 
277 
278 //- Estimate angle of vec in coordinate system (e0, e1, e0^e1).
279 // Is guaranteed to return increasing number but is not correct
280 // angle. Used for sorting angles. All input vectors need to be normalized.
281 //
282 // Calculates scalar which increases with angle going from e0 to vec in
283 // the coordinate system e0, e1, e0^e1
284 //
285 // Jumps from 2*pi -> 0 at -SMALL so hopefully parallel vectors with small
286 // rounding errors should still get the same quadrant.
287 //
288 inline scalar pseudoAngle
289 (
290  const vector& e0,
291  const vector& e1,
292  const vector& vec
293 )
294 {
295  scalar cos = vec & e0;
296  scalar sin = vec & e1;
297 
298  if (sin < -SMALL)
299  {
300  return (3.0 + cos)*constant::mathematical::piByTwo;
301  }
302  else
303  {
304  return (1.0 - cos)*constant::mathematical::piByTwo;
305  }
306 }
307 
308 
309 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
310 
311 } // End namespace Foam
312 
313 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
314 
315 #endif
316 
317 // ************************************************************************* //
const Cmpt & yz() const
Definition: SymmTensorI.H:111
const Cmpt & yx() const
Definition: TensorI.H:174
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:202
Templated 3D symmetric tensor derived from VectorSpace adding construction from 6 components...
Definition: SymmTensor.H:53
const Cmpt & xz() const
Definition: TensorI.H:167
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:289
const Cmpt & xx() const
Definition: SymmTensorI.H:87
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const Cmpt & yy() const
Definition: TensorI.H:181
const Cmpt & xx() const
Definition: TensorI.H:153
const Cmpt & yz() const
Definition: TensorI.H:188
const Cmpt & z() const
Definition: VectorI.H:87
sphericalTensor transformMask< sphericalTensor >(const symmTensor &st)
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from 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:95
const Cmpt & zx() const
Definition: TensorI.H:195
Type1 transformMask(const Type2 &t)
Definition: transform.H:258
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
tensor Rz(const scalar &omega)
Rotational transformation tensor about the z-axis by omega radians.
Definition: transform.H:109
const dimensionedScalar c
Speed of light in a vacuum.
tensor Ra(const vector &a, const scalar omega)
Rotational transformation tensor about axis a by omega radians.
Definition: transform.H:123
const Cmpt & xy() const
Definition: SymmTensorI.H:93
const scalar piByTwo(0.5 *pi)
const Cmpt & xy() const
Definition: TensorI.H:160
const Cmpt & zz() const
Definition: TensorI.H:209
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:81
Namespace for OpenFOAM.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:477