FreeFOAM The Cross-Platform CFD Toolkit
tensor.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) 1991-2010 OpenCFD Ltd.
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 "tensor.H"
28 
29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 template<>
37 const char* const tensor::typeName = "tensor";
38 
39 template<>
40 const char* tensor::componentNames[] =
41 {
42  "xx", "xy", "xz",
43  "yx", "yy", "yz",
44  "zx", "zy", "zz"
45 };
46 
47 template<>
48 const tensor tensor::zero
49 (
50  0, 0, 0,
51  0, 0, 0,
52  0, 0, 0
53 );
54 
55 template<>
56 const tensor tensor::one
57 (
58  1, 1, 1,
59  1, 1, 1,
60  1, 1, 1
61 );
62 
63 template<>
64 const tensor tensor::max
65 (
66  VGREAT, VGREAT, VGREAT,
67  VGREAT, VGREAT, VGREAT,
68  VGREAT, VGREAT, VGREAT
69 );
70 
71 template<>
72 const tensor tensor::min
73 (
74  -VGREAT, -VGREAT, -VGREAT,
75  -VGREAT, -VGREAT, -VGREAT,
76  -VGREAT, -VGREAT, -VGREAT
77 );
78 
79 
80 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
81 
82 // Return eigenvalues in ascending order of absolute values
84 {
85  scalar i = 0;
86  scalar ii = 0;
87  scalar iii = 0;
88 
89  if
90  (
91  (
92  mag(t.xy()) + mag(t.xz()) + mag(t.yx())
93  + mag(t.yz()) + mag(t.zx()) + mag(t.zy())
94  )
95  < SMALL
96  )
97  {
98  // diagonal matrix
99  i = t.xx();
100  ii = t.yy();
101  iii = t.zz();
102  }
103  else
104  {
105  scalar a = -t.xx() - t.yy() - t.zz();
106 
107  scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
108  - t.xy()*t.yx() - t.xz()*t.zx() - t.yz()*t.zy();
109 
110  scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.zx()
111  - t.xz()*t.yx()*t.zy() + t.xz()*t.yy()*t.zx()
112  + t.xy()*t.yx()*t.zz() + t.xx()*t.yz()*t.zy();
113 
114  // If there is a zero root
115  if (mag(c) < 1.0e-100)
116  {
117  scalar disc = sqr(a) - 4*b;
118 
119  if (disc >= -SMALL)
120  {
121  scalar q = -0.5*sqrt(max(0.0, disc));
122 
123  i = 0;
124  ii = -0.5*a + q;
125  iii = -0.5*a - q;
126  }
127  else
128  {
129  FatalErrorIn("eigenValues(const tensor&)")
130  << "zero and complex eigenvalues in tensor: " << t
131  << abort(FatalError);
132  }
133  }
134  else
135  {
136  scalar Q = (a*a - 3*b)/9;
137  scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
138 
139  scalar R2 = sqr(R);
140  scalar Q3 = pow3(Q);
141 
142  // Three different real roots
143  if (R2 < Q3)
144  {
145  scalar sqrtQ = sqrt(Q);
146  scalar theta = acos(R/(Q*sqrtQ));
147 
148  scalar m2SqrtQ = -2*sqrtQ;
149  scalar aBy3 = a/3;
150 
151  i = m2SqrtQ*cos(theta/3) - aBy3;
152  ii = m2SqrtQ*cos((theta + mathematicalConstant::twoPi)/3)
153  - aBy3;
154  iii = m2SqrtQ*cos((theta - mathematicalConstant::twoPi)/3)
155  - aBy3;
156  }
157  else
158  {
159  scalar A = cbrt(R + sqrt(R2 - Q3));
160 
161  // Three equal real roots
162  if (A < SMALL)
163  {
164  scalar root = -a/3;
165  return vector(root, root, root);
166  }
167  else
168  {
169  // Complex roots
170  WarningIn("eigenValues(const tensor&)")
171  << "complex eigenvalues detected for tensor: " << t
172  << endl;
173 
174  return vector::zero;
175  }
176  }
177  }
178  }
179 
180 
181  // Sort the eigenvalues into ascending order
182  if (mag(i) > mag(ii))
183  {
184  Swap(i, ii);
185  }
186 
187  if (mag(ii) > mag(iii))
188  {
189  Swap(ii, iii);
190  }
191 
192  if (mag(i) > mag(ii))
193  {
194  Swap(i, ii);
195  }
196 
197  return vector(i, ii, iii);
198 }
199 
200 
201 vector eigenVector(const tensor& t, const scalar lambda)
202 {
203  if (mag(lambda) < SMALL)
204  {
205  return vector::zero;
206  }
207 
208  // Construct the matrix for the eigenvector problem
209  tensor A(t - lambda*I);
210 
211  // Calculate the sub-determinants of the 3 components
212  scalar sd0 = A.yy()*A.zz() - A.yz()*A.zy();
213  scalar sd1 = A.xx()*A.zz() - A.xz()*A.zx();
214  scalar sd2 = A.xx()*A.yy() - A.xy()*A.yx();
215 
216  scalar magSd0 = mag(sd0);
217  scalar magSd1 = mag(sd1);
218  scalar magSd2 = mag(sd2);
219 
220  // Evaluate the eigenvector using the largest sub-determinant
221  if (magSd0 > magSd1 && magSd0 > magSd2 && magSd0 > SMALL)
222  {
223  vector ev
224  (
225  1,
226  (A.yz()*A.zx() - A.zz()*A.yx())/sd0,
227  (A.zy()*A.yx() - A.yy()*A.zx())/sd0
228  );
229  ev /= mag(ev);
230 
231  return ev;
232  }
233  else if (magSd1 > magSd2 && magSd1 > SMALL)
234  {
235  vector ev
236  (
237  (A.xz()*A.zy() - A.zz()*A.xy())/sd1,
238  1,
239  (A.zx()*A.xy() - A.xx()*A.zy())/sd1
240  );
241  ev /= mag(ev);
242 
243  return ev;
244  }
245  else if (magSd2 > SMALL)
246  {
247  vector ev
248  (
249  (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
250  (A.yx()*A.xz() - A.xx()*A.yz())/sd2,
251  1
252  );
253  ev /= mag(ev);
254 
255  return ev;
256  }
257  else
258  {
259  return vector::zero;
260  }
261 }
262 
263 
265 {
266  vector evals(eigenValues(t));
267 
268  tensor evs
269  (
270  eigenVector(t, evals.x()),
271  eigenVector(t, evals.y()),
272  eigenVector(t, evals.z())
273  );
274 
275  return evs;
276 }
277 
278 
279 // Return eigenvalues in ascending order of absolute values
281 {
282  scalar i = 0;
283  scalar ii = 0;
284  scalar iii = 0;
285 
286  if
287  (
288  (
289  mag(t.xy()) + mag(t.xz()) + mag(t.xy())
290  + mag(t.yz()) + mag(t.xz()) + mag(t.yz())
291  )
292  < SMALL
293  )
294  {
295  // diagonal matrix
296  i = t.xx();
297  ii = t.yy();
298  iii = t.zz();
299  }
300  else
301  {
302  scalar a = -t.xx() - t.yy() - t.zz();
303 
304  scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
305  - t.xy()*t.xy() - t.xz()*t.xz() - t.yz()*t.yz();
306 
307  scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.xz()
308  - t.xz()*t.xy()*t.yz() + t.xz()*t.yy()*t.xz()
309  + t.xy()*t.xy()*t.zz() + t.xx()*t.yz()*t.yz();
310 
311  // If there is a zero root
312  if (mag(c) < 1.0e-100)
313  {
314  scalar disc = sqr(a) - 4*b;
315 
316  if (disc >= -SMALL)
317  {
318  scalar q = -0.5*sqrt(max(0.0, disc));
319 
320  i = 0;
321  ii = -0.5*a + q;
322  iii = -0.5*a - q;
323  }
324  else
325  {
326  FatalErrorIn("eigenValues(const tensor&)")
327  << "zero and complex eigenvalues in tensor: " << t
328  << abort(FatalError);
329  }
330  }
331  else
332  {
333  scalar Q = (a*a - 3*b)/9;
334  scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
335 
336  scalar R2 = sqr(R);
337  scalar Q3 = pow3(Q);
338 
339  // Three different real roots
340  if (R2 < Q3)
341  {
342  scalar sqrtQ = sqrt(Q);
343  scalar theta = acos(R/(Q*sqrtQ));
344 
345  scalar m2SqrtQ = -2*sqrtQ;
346  scalar aBy3 = a/3;
347 
348  i = m2SqrtQ*cos(theta/3) - aBy3;
349  ii = m2SqrtQ*cos((theta + mathematicalConstant::twoPi)/3)
350  - aBy3;
351  iii = m2SqrtQ*cos((theta - mathematicalConstant::twoPi)/3)
352  - aBy3;
353  }
354  else
355  {
356  scalar A = cbrt(R + sqrt(R2 - Q3));
357 
358  // Three equal real roots
359  if (A < SMALL)
360  {
361  scalar root = -a/3;
362  return vector(root, root, root);
363  }
364  else
365  {
366  // Complex roots
367  WarningIn("eigenValues(const symmTensor&)")
368  << "complex eigenvalues detected for symmTensor: " << t
369  << endl;
370 
371  return vector::zero;
372  }
373  }
374  }
375  }
376 
377 
378  // Sort the eigenvalues into ascending order
379  if (mag(i) > mag(ii))
380  {
381  Swap(i, ii);
382  }
383 
384  if (mag(ii) > mag(iii))
385  {
386  Swap(ii, iii);
387  }
388 
389  if (mag(i) > mag(ii))
390  {
391  Swap(i, ii);
392  }
393 
394  return vector(i, ii, iii);
395 }
396 
397 
398 vector eigenVector(const symmTensor& t, const scalar lambda)
399 {
400  if (mag(lambda) < SMALL)
401  {
402  return vector::zero;
403  }
404 
405  // Construct the matrix for the eigenvector problem
406  symmTensor A(t - lambda*I);
407 
408  // Calculate the sub-determinants of the 3 components
409  scalar sd0 = A.yy()*A.zz() - A.yz()*A.yz();
410  scalar sd1 = A.xx()*A.zz() - A.xz()*A.xz();
411  scalar sd2 = A.xx()*A.yy() - A.xy()*A.xy();
412 
413  scalar magSd0 = mag(sd0);
414  scalar magSd1 = mag(sd1);
415  scalar magSd2 = mag(sd2);
416 
417  // Evaluate the eigenvector using the largest sub-determinant
418  if (magSd0 > magSd1 && magSd0 > magSd2 && magSd0 > SMALL)
419  {
420  vector ev
421  (
422  1,
423  (A.yz()*A.xz() - A.zz()*A.xy())/sd0,
424  (A.yz()*A.xy() - A.yy()*A.xz())/sd0
425  );
426  ev /= mag(ev);
427 
428  return ev;
429  }
430  else if (magSd1 > magSd2 && magSd1 > SMALL)
431  {
432  vector ev
433  (
434  (A.xz()*A.yz() - A.zz()*A.xy())/sd1,
435  1,
436  (A.xz()*A.xy() - A.xx()*A.yz())/sd1
437  );
438  ev /= mag(ev);
439 
440  return ev;
441  }
442  else if (magSd2 > SMALL)
443  {
444  vector ev
445  (
446  (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
447  (A.xy()*A.xz() - A.xx()*A.yz())/sd2,
448  1
449  );
450  ev /= mag(ev);
451 
452  return ev;
453  }
454  else
455  {
456  return vector::zero;
457  }
458 }
459 
460 
462 {
463  vector evals(eigenValues(t));
464 
465  tensor evs
466  (
467  eigenVector(t, evals.x()),
468  eigenVector(t, evals.y()),
469  eigenVector(t, evals.z())
470  );
471 
472  return evs;
473 }
474 
475 
476 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
477 
478 } // End namespace Foam
479 
480 // ************************ vim: set sw=4 sts=4 et: ************************ //