triad.C
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) 2012-2014 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 \*---------------------------------------------------------------------------*/
25 
26 #include "triad.H"
27 #include "transform.H"
28 #include "quaternion.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 template<>
36 const char* const triad::Vector<vector>::typeName = "triad";
37 
38 template<>
39 const char* triad::Vector<vector>::componentNames[] = {"x", "y", "z"};
40 
41 const triad triad::zero
42 (
43  vector(0, 0, 0),
44  vector(0, 0, 0),
45  vector(0, 0, 0)
46 );
47 
48 const triad triad::one
49 (
50  vector(1, 1, 1),
51  vector(1, 1, 1),
52  vector(1, 1, 1)
53 );
54 
55 const triad triad::max
56 (
57  vector(VGREAT, VGREAT, VGREAT),
58  vector(VGREAT, VGREAT, VGREAT),
59  vector(VGREAT, VGREAT, VGREAT)
60 );
61 
62 const triad triad::min
63 (
64  vector(-VGREAT, -VGREAT, -VGREAT),
65  vector(-VGREAT, -VGREAT, -VGREAT),
66  vector(-VGREAT, -VGREAT, -VGREAT)
67 );
68 
69 const triad triad::rootMax
70 (
71  vector(ROOTVGREAT, ROOTVGREAT, ROOTVGREAT),
72  vector(ROOTVGREAT, ROOTVGREAT, ROOTVGREAT),
73  vector(ROOTVGREAT, ROOTVGREAT, ROOTVGREAT)
74 );
75 
76 const triad triad::rootMin
77 (
78  vector(-ROOTVGREAT, -ROOTVGREAT, -ROOTVGREAT),
79  vector(-ROOTVGREAT, -ROOTVGREAT, -ROOTVGREAT),
80  vector(-ROOTVGREAT, -ROOTVGREAT, -ROOTVGREAT)
81 );
82 
83 const triad triad::I
84 (
85  vector(1, 0, 0),
86  vector(0, 1, 0),
87  vector(0, 0, 1)
88 );
89 
90 const triad triad::unset
91 (
92  vector(VGREAT, VGREAT, VGREAT),
93  vector(VGREAT, VGREAT, VGREAT),
94  vector(VGREAT, VGREAT, VGREAT)
95 );
96 
97 }
98 
99 
100 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
101 
103 {
104  tensor Rt(q.R().T());
105  x() = Rt.x();
106  y() = Rt.y();
107  z() = Rt.z();
108 }
109 
110 
112 {
113  x() = t.x();
114  y() = t.y();
115  z() = t.z();
116 }
117 
118 
119 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
120 
122 {
123  // Hack for 2D z-slab cases
124  // if (!set(2))
125  // {
126  // operator[](2) = vector(0, 0, 1);
127  // }
128 
129  // If only two of the axes are set, set the third
130  if (set(0) && set(1) && !set(2))
131  {
132  operator[](2) = orthogonal(operator[](0), operator[](1));
133  }
134  else if (set(0) && set(2) && !set(1))
135  {
136  operator[](1) = orthogonal(operator[](0), operator[](2));
137  }
138  else if (set(1) && set(2) && !set(0))
139  {
140  operator[](0) = orthogonal(operator[](1), operator[](2));
141  }
142 
143  // If all the axes are set
144  if (set())
145  {
146  for (int i=0; i<2; i++)
147  {
148  scalar o01 = mag(operator[](0) & operator[](1));
149  scalar o02 = mag(operator[](0) & operator[](2));
150  scalar o12 = mag(operator[](1) & operator[](2));
151 
152  if (o01 < o02 && o01 < o12)
153  {
154  operator[](2) = orthogonal(operator[](0), operator[](1));
155 
156  // if (o02 < o12)
157  // {
158  // operator[](1) = orthogonal(operator[](0), operator[](2));
159  // }
160  // else
161  // {
162  // operator[](0) = orthogonal(operator[](1), operator[](2));
163  // }
164  }
165  else if (o02 < o12)
166  {
167  operator[](1) = orthogonal(operator[](0), operator[](2));
168 
169  // if (o01 < o12)
170  // {
171  // operator[](2) = orthogonal(operator[](0), operator[](1));
172  // }
173  // else
174  // {
175  // operator[](0) = orthogonal(operator[](1), operator[](2));
176  // }
177  }
178  else
179  {
180  operator[](0) = orthogonal(operator[](1), operator[](2));
181 
182  // if (o02 < o01)
183  // {
184  // operator[](1) = orthogonal(operator[](0), operator[](2));
185  // }
186  // else
187  // {
188  // operator[](2) = orthogonal(operator[](0), operator[](1));
189  // }
190  }
191  }
192  }
193 }
194 
195 
197 {
198  bool preset[3];
199 
200  for (direction i=0; i<3; i++)
201  {
202  if (t2.set(i) && !set(i))
203  {
204  operator[](i) = t2.operator[](i);
205  preset[i] = true;
206  }
207  else
208  {
209  preset[i] = false;
210  }
211  }
212 
213  if (set() && t2.set())
214  {
215  direction correspondance[3];
216  short signd[3];
217 
218  for (direction i=0; i<3; i++)
219  {
220  if (preset[i])
221  {
222  signd[i] = 0;
223  continue;
224  }
225 
226  scalar mostAligned = -1;
227  for (direction j=0; j<3; j++)
228  {
229  bool set = false;
230  for (direction k=0; k<i; k++)
231  {
232  if (correspondance[k] == j)
233  {
234  set = true;
235  break;
236  }
237  }
238 
239  if (!set)
240  {
241  scalar a = operator[](i) & t2.operator[](j);
242  scalar maga = mag(a);
243 
244  if (maga > mostAligned)
245  {
246  correspondance[i] = j;
247  mostAligned = maga;
248  signd[i] = sign(a);
249  }
250  }
251  }
252 
253  operator[](i) += signd[i]*t2.operator[](correspondance[i]);
254  }
255  }
256 }
257 
258 
260 {
261  if (set())
262  {
263  vector mostAligned
264  (
265  mag(v & operator[](0)),
266  mag(v & operator[](1)),
267  mag(v & operator[](2))
268  );
269 
270  scalar mav;
271 
272  if
273  (
274  mostAligned.x() > mostAligned.y()
275  && mostAligned.x() > mostAligned.z()
276  )
277  {
278  mav = mostAligned.x();
279  mostAligned = operator[](0);
280  }
281  else if (mostAligned.y() > mostAligned.z())
282  {
283  mav = mostAligned.y();
284  mostAligned = operator[](1);
285  }
286  else
287  {
288  mav = mostAligned.z();
289  mostAligned = operator[](2);
290  }
291 
292  if (mav < 0.99)
293  {
294  tensor R(rotationTensor(mostAligned, v));
295 
296  operator[](0) = transform(R, operator[](0));
297  operator[](1) = transform(R, operator[](1));
298  operator[](2) = transform(R, operator[](2));
299  }
300  }
301 }
302 
303 
305 {
306  if (!this->set())
307  {
308  return *this;
309  }
310 
311  triad t;
312 
313  if
314  (
315  mag(operator[](0).x()) > mag(operator[](1).x())
316  && mag(operator[](0).x()) > mag(operator[](2).x())
317  )
318  {
319  t[0] = operator[](0);
320 
321  if (mag(operator[](1).y()) > mag(operator[](2).y()))
322  {
323  t[1] = operator[](1);
324  t[2] = operator[](2);
325  }
326  else
327  {
328  t[1] = operator[](2);
329  t[2] = operator[](1);
330  }
331  }
332  else if
333  (
334  mag(operator[](1).x()) > mag(operator[](2).x())
335  )
336  {
337  t[0] = operator[](1);
338 
339  if (mag(operator[](0).y()) > mag(operator[](2).y()))
340  {
341  t[1] = operator[](0);
342  t[2] = operator[](2);
343  }
344  else
345  {
346  t[1] = operator[](2);
347  t[2] = operator[](0);
348  }
349  }
350  else
351  {
352  t[0] = operator[](2);
353 
354  if (mag(operator[](0).y()) > mag(operator[](1).y()))
355  {
356  t[1] = operator[](0);
357  t[2] = operator[](1);
358  }
359  else
360  {
361  t[1] = operator[](1);
362  t[2] = operator[](0);
363  }
364  }
365 
366  if (t[0].x() < 0) t[0] *= -1;
367  if (t[1].y() < 0) t[1] *= -1;
368  if (t[2].z() < 0) t[2] *= -1;
369 
370  return t;
371 }
372 
373 
374 
376 {
377  tensor R;
378 
379  R.xx() = x().x();
380  R.xy() = y().x();
381  R.xz() = z().x();
382 
383  R.yx() = x().y();
384  R.yy() = y().y();
385  R.yz() = z().y();
386 
387  R.zx() = x().z();
388  R.zy() = y().z();
389  R.zz() = z().z();
390 
391  return quaternion(R);
392 }
393 
394 
395 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
396 
398 {
399  x() = t.x();
400  y() = t.y();
401  z() = t.z();
402 }
403 
404 
405 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
406 
407 Foam::scalar Foam::diff(const triad& A, const triad& B)
408 {
409  triad tmpA = A.sortxyz();
410  triad tmpB = B.sortxyz();
411 
412  scalar sumDifference = 0;
413 
414  for (direction dir = 0; dir < 3; dir++)
415  {
416  if (!tmpA.set(dir) || !tmpB.set(dir))
417  {
418  continue;
419  }
420 
421  scalar cosPhi =
422  (tmpA[dir] & tmpB[dir])
423  /(mag(tmpA[dir])*mag(tmpA[dir]) + SMALL);
424 
425  cosPhi = min(max(cosPhi, -1), 1);
426 
427  sumDifference += mag(cosPhi - 1);
428  }
429 
430  return (sumDifference/3);
431 }
432 
433 
434 // ************************************************************************* //
void align(const vector &v)
Align this triad with the given vector v.
Definition: triad.C:259
const Cmpt & xy() const
Definition: TensorI.H:167
unsigned char direction
Definition: direction.H:43
bool set(const direction d) const
Is the vector in the direction d set.
Definition: triadI.H:61
static const triad I
Definition: triad.H:104
const Cmpt & zz() const
Definition: TensorI.H:216
dimensioned< scalar > mag(const dimensioned< Type > &)
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
Definition: triad.C:407
#define R(A, B, C, D, E, F, K, M)
Quaternion class used to perform rotations in 3D space.
Definition: quaternion.H:60
static const triad rootMin
Definition: triad.H:103
triad sortxyz() const
Sort the axes such that they are closest to the x, y and z axes.
Definition: triad.C:304
const Cmpt & yy() const
Definition: TensorI.H:188
dimensionedScalar sign(const dimensionedScalar &ds)
triad()
Construct null.
Definition: triadI.H:28
Vector< Cmpt > x() const
Definition: TensorI.H:121
Namespace for OpenFOAM.
static const triad one
Definition: triad.H:99
Vector< Cmpt > y() const
Definition: TensorI.H:128
tensor R() const
The rotation tensor corresponding the quaternion.
Definition: quaternionI.H:235
void operator+=(const triad &t2)
Add the triad t2 to this triad.
Definition: triad.C:196
static const triad zero
Definition: triad.H:98
const Cmpt & y() const
Definition: VectorI.H:71
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:465
T cosPhi(const Vector< T > &a, const Vector< T > &b, const T &tolerance=SMALL)
Calculate angle between a and b in radians.
Definition: vectorTools.H:105
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
label k
Boltzmann constant.
Vector< Cmpt > z() const
Definition: TensorI.H:135
static const triad min
Definition: triad.H:101
const Cmpt & x() const
Definition: VectorI.H:65
const Cmpt & yx() const
Definition: TensorI.H:181
void operator=(const Vector< vector > &)
Definition: triadI.H:121
Representation of a 3D Cartesian coordinate system as a Vector of vectors.
Definition: triad.H:65
const Cmpt & z() const
Definition: VectorI.H:77
static const triad unset
Definition: triad.H:105
static const triad max
Definition: triad.H:100
friend Ostream & operator(Ostream &, const VectorSpace< Vector< vector >, vector, nCmpt > &)
scalar y
const Cmpt & zy() const
Definition: TensorI.H:209
3D tensor transformation operations.
Tensor< Cmpt > T() const
Transpose.
Definition: TensorI.H:286
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
void orthogonalize()
Orthogonalize this triad so that it is ortho-normal.
Definition: triad.C:121
tensor rotationTensor(const vector &n1, const vector &n2)
Definition: transform.H:46
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const Cmpt & xx() const
Definition: TensorI.H:160
const Cmpt & yz() const
Definition: TensorI.H:195
const Cmpt & xz() const
Definition: TensorI.H:174
const Cmpt & zx() const
Definition: TensorI.H:202
static const triad rootMax
Definition: triad.H:102