triad.C
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) 2012-2021 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  const scalar o01
136  (
137  (set(0) && set(1)) ? mag(operator[](0) & operator[](1)) : vGreat
138  );
139  const scalar o02
140  (
141  (set(0) && set(2)) ? mag(operator[](0) & operator[](2)) : vGreat
142  );
143  const scalar o12
144  (
145  (set(1) && set(2)) ? mag(operator[](1) & operator[](2)) : vGreat
146  );
147 
148  if (o01 < o02 && o01 < o12)
149  {
150  operator[](2) = orthogonal(operator[](0), operator[](1));
151 
152  // if (o02 < o12)
153  // {
154  // operator[](1) = orthogonal(operator[](0), operator[](2));
155  // }
156  // else
157  // {
158  // operator[](0) = orthogonal(operator[](1), operator[](2));
159  // }
160  }
161  else if (o02 < o12)
162  {
163  operator[](1) = orthogonal(operator[](0), operator[](2));
164 
165  // if (o01 < o12)
166  // {
167  // operator[](2) = orthogonal(operator[](0), operator[](1));
168  // }
169  // else
170  // {
171  // operator[](0) = orthogonal(operator[](1), operator[](2));
172  // }
173  }
174  else
175  {
176  operator[](0) = orthogonal(operator[](1), operator[](2));
177 
178  // if (o02 < o01)
179  // {
180  // operator[](1) = orthogonal(operator[](0), operator[](2));
181  // }
182  // else
183  // {
184  // operator[](2) = orthogonal(operator[](0), operator[](1));
185  // }
186  }
187  }
188  }
189 }
190 
191 
193 {
194  bool preset[3];
195 
196  for (direction i=0; i<3; i++)
197  {
198  if (t2.set(i) && !set(i))
199  {
200  operator[](i) = t2.operator[](i);
201  preset[i] = true;
202  }
203  else
204  {
205  preset[i] = false;
206  }
207  }
208 
209  if (set() && t2.set())
210  {
211  direction correspondence[3]{0, 0, 0};
212  short signd[3];
213 
214  for (direction i=0; i<3; i++)
215  {
216  if (preset[i])
217  {
218  signd[i] = 0;
219  continue;
220  }
221 
222  scalar mostAligned = -1;
223  for (direction j=0; j<3; j++)
224  {
225  bool set = false;
226  for (direction k=0; k<i; k++)
227  {
228  if (correspondence[k] == j)
229  {
230  set = true;
231  break;
232  }
233  }
234 
235  if (!set)
236  {
237  scalar a = operator[](i) & t2.operator[](j);
238  scalar maga = mag(a);
239 
240  if (maga > mostAligned)
241  {
242  correspondence[i] = j;
243  mostAligned = maga;
244  signd[i] = sign(a);
245  }
246  }
247  }
248 
249  operator[](i) += signd[i]*t2.operator[](correspondence[i]);
250  }
251  }
252 }
253 
254 
256 {
257  if (set())
258  {
259  vector mostAligned
260  (
261  mag(v & operator[](0)),
262  mag(v & operator[](1)),
263  mag(v & operator[](2))
264  );
265 
266  scalar mav;
267 
268  if
269  (
270  mostAligned.x() > mostAligned.y()
271  && mostAligned.x() > mostAligned.z()
272  )
273  {
274  mav = mostAligned.x();
275  mostAligned = operator[](0);
276  }
277  else if (mostAligned.y() > mostAligned.z())
278  {
279  mav = mostAligned.y();
280  mostAligned = operator[](1);
281  }
282  else
283  {
284  mav = mostAligned.z();
285  mostAligned = operator[](2);
286  }
287 
288  if (mav < 0.99)
289  {
290  tensor R(rotationTensor(mostAligned, v));
291 
292  operator[](0) = transform(R, operator[](0));
293  operator[](1) = transform(R, operator[](1));
294  operator[](2) = transform(R, operator[](2));
295  }
296  }
297 }
298 
299 
301 {
302  if (!this->set())
303  {
304  return *this;
305  }
306 
307  triad t;
308 
309  if
310  (
311  mag(operator[](0).x()) > mag(operator[](1).x())
312  && mag(operator[](0).x()) > mag(operator[](2).x())
313  )
314  {
315  t[0] = operator[](0);
316 
317  if (mag(operator[](1).y()) > mag(operator[](2).y()))
318  {
319  t[1] = operator[](1);
320  t[2] = operator[](2);
321  }
322  else
323  {
324  t[1] = operator[](2);
325  t[2] = operator[](1);
326  }
327  }
328  else if
329  (
330  mag(operator[](1).x()) > mag(operator[](2).x())
331  )
332  {
333  t[0] = operator[](1);
334 
335  if (mag(operator[](0).y()) > mag(operator[](2).y()))
336  {
337  t[1] = operator[](0);
338  t[2] = operator[](2);
339  }
340  else
341  {
342  t[1] = operator[](2);
343  t[2] = operator[](0);
344  }
345  }
346  else
347  {
348  t[0] = operator[](2);
349 
350  if (mag(operator[](0).y()) > mag(operator[](1).y()))
351  {
352  t[1] = operator[](0);
353  t[2] = operator[](1);
354  }
355  else
356  {
357  t[1] = operator[](1);
358  t[2] = operator[](0);
359  }
360  }
361 
362  if (t[0].x() < 0) t[0] *= -1;
363  if (t[1].y() < 0) t[1] *= -1;
364  if (t[2].z() < 0) t[2] *= -1;
365 
366  return t;
367 }
368 
369 
370 
372 {
373  tensor R;
374 
375  R.xx() = x().x();
376  R.xy() = y().x();
377  R.xz() = z().x();
378 
379  R.yx() = x().y();
380  R.yy() = y().y();
381  R.yz() = z().y();
382 
383  R.zx() = x().z();
384  R.zy() = y().z();
385  R.zz() = z().z();
386 
387  return quaternion(R);
388 }
389 
390 
391 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
392 
394 {
395  x() = t.x();
396  y() = t.y();
397  z() = t.z();
398 }
399 
400 
401 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
402 
403 Foam::scalar Foam::diff(const triad& A, const triad& B)
404 {
405  triad tmpA = A.sortxyz();
406  triad tmpB = B.sortxyz();
407 
408  scalar sumDifference = 0;
409 
410  for (direction dir = 0; dir < 3; dir++)
411  {
412  if (!tmpA.set(dir) || !tmpB.set(dir))
413  {
414  continue;
415  }
416 
417  scalar cosPhi =
418  (tmpA[dir] & tmpB[dir])
419  /(mag(tmpA[dir])*mag(tmpA[dir]) + small);
420 
421  cosPhi = min(max(cosPhi, -1), 1);
422 
423  sumDifference += mag(cosPhi - 1);
424  }
425 
426  return (sumDifference/3);
427 }
428 
429 
430 // ************************************************************************* //
const Cmpt & yx() const
Definition: TensorI.H:184
dimensionedScalar sign(const dimensionedScalar &ds)
static const char *const typeName
Definition: VectorSpace.H:111
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
Definition: triad.C:403
tensor R() const
The rotation tensor corresponding the quaternion.
Definition: quaternionI.H:323
const Cmpt & zy() const
Definition: TensorI.H:212
const Cmpt & xz() const
Definition: TensorI.H:177
static const Form max
Definition: VectorSpace.H:115
void align(const vector &v)
Align this triad with the given vector v.
Definition: triad.C:255
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:331
Vector< Cmpt > x() const
Definition: TensorI.H:289
static const char *const componentNames[]
Definition: VectorSpace.H:112
const Cmpt & yy() const
Definition: TensorI.H:191
const Cmpt & xx() const
Definition: TensorI.H:163
const Cmpt & yz() const
Definition: TensorI.H:198
static const Form rootMin
Definition: VectorSpace.H:118
Vector< Cmpt > y() const
Definition: TensorI.H:296
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.
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from unit vector n1 to n2.
Definition: transform.H:47
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
const vector & y() const
Definition: VectorI.H:81
static const Form min
Definition: VectorSpace.H:116
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
void orthogonalise()
Orthogonalise this triad so that it is ortho-normal.
Definition: triad.C:108
Vector< Cmpt > z() const
Definition: TensorI.H:303
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:300
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:205
static const Form rootMax
Definition: VectorSpace.H:117
void operator+=(const triad &t2)
Add the triad t2 to this triad.
Definition: triad.C:192
#define R(A, B, C, D, E, F, K, M)
static const Form one
Definition: VectorSpace.H:114
dimensioned< scalar > mag(const dimensioned< Type > &)
const vector & operator[](const direction) const
Definition: VectorSpaceI.H:190
const Cmpt & xy() const
Definition: TensorI.H:170
const Cmpt & zz() const
Definition: TensorI.H:219
static const Form zero
Definition: VectorSpace.H:113
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:477