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