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-2013 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 
111 Foam::triad::triad(const tensor& t)
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 
259 void Foam::triad::align(const vector& v)
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 
304 Foam::triad Foam::triad::sortxyz() const
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 
375 Foam::triad::operator quaternion() const
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 
397 void Foam::triad::operator=(const tensor& t)
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 // ************************************************************************* //