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) 2011-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 "tensor.H"
27 #include "mathematicalConstants.H"
28 
29 using namespace Foam::constant::mathematical;
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35  template<>
36  const char* const tensor::typeName = "tensor";
37 
38  template<>
39  const char* tensor::componentNames[] =
40  {
41  "xx", "xy", "xz",
42  "yx", "yy", "yz",
43  "zx", "zy", "zz"
44  };
45 
46  template<>
47  const tensor tensor::zero
48  (
49  0, 0, 0,
50  0, 0, 0,
51  0, 0, 0
52  );
53 
54  template<>
55  const tensor tensor::one
56  (
57  1, 1, 1,
58  1, 1, 1,
59  1, 1, 1
60  );
61 
62  template<>
63  const tensor tensor::max
64  (
65  VGREAT, VGREAT, VGREAT,
66  VGREAT, VGREAT, VGREAT,
67  VGREAT, VGREAT, VGREAT
68  );
69 
70  template<>
71  const tensor tensor::min
72  (
73  -VGREAT, -VGREAT, -VGREAT,
74  -VGREAT, -VGREAT, -VGREAT,
75  -VGREAT, -VGREAT, -VGREAT
76  );
77 
78  template<>
79  const tensor tensor::I
80  (
81  1, 0, 0,
82  0, 1, 0,
83  0, 0, 1
84  );
85 }
86 
87 
88 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
89 
91 {
92  // The eigenvalues
93  scalar i, ii, iii;
94 
95  // diagonal matrix
96  if
97  (
98  (
99  mag(t.xy()) + mag(t.xz()) + mag(t.yx())
100  + mag(t.yz()) + mag(t.zx()) + mag(t.zy())
101  )
102  < SMALL
103  )
104  {
105  i = t.xx();
106  ii = t.yy();
107  iii = t.zz();
108  }
109 
110  // non-diagonal matrix
111  else
112  {
113  // Coefficients of the characteristic polynmial
114  // x^3 + a*x^2 + b*x + c = 0
115  scalar a =
116  - t.xx() - t.yy() - t.zz();
117 
118  scalar b =
119  t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
120  - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz();
121 
122  scalar c =
123  - t.xx()*t.yy()*t.zz()
124  - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx()
125  + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx();
126 
127  // Auxillary variables
128  scalar aBy3 = a/3;
129 
130  scalar P = (a*a - 3*b)/9; // == -p_wikipedia/3
131  scalar PPP = P*P*P;
132 
133  scalar Q = (2*a*a*a - 9*a*b + 27*c)/54; // == q_wikipedia/2
134  scalar QQ = Q*Q;
135 
136  // Three identical roots
137  if (mag(P) < SMALL && mag(Q) < SMALL)
138  {
139  return vector(- aBy3, - aBy3, - aBy3);
140  }
141 
142  // Two identical roots and one distinct root
143  else if (mag(PPP/QQ - 1) < SMALL)
144  {
145  scalar sqrtP = sqrt(P);
146  scalar signQ = sign(Q);
147 
148  i = ii = signQ*sqrtP - aBy3;
149  iii = - 2*signQ*sqrtP - aBy3;
150  }
151 
152  // Three distinct roots
153  else if (PPP > QQ)
154  {
155  scalar sqrtP = sqrt(P);
156  scalar value = cos(acos(Q/sqrt(PPP))/3);
157  scalar delta = sqrt(3 - 3*value*value);
158 
159  i = - 2*sqrtP*value - aBy3;
160  ii = sqrtP*(value + delta) - aBy3;
161  iii = sqrtP*(value - delta) - aBy3;
162  }
163 
164  // One real root, two imaginary roots
165  // based on the above logic, PPP must be less than QQ
166  else
167  {
168  WarningIn("eigenValues(const tensor&)")
169  << "complex eigenvalues detected for tensor: " << t
170  << endl;
171 
172  if (mag(P) < SMALL)
173  {
174  i = cbrt(QQ/2);
175  }
176  else
177  {
178  scalar w = cbrt(- Q - sqrt(QQ - PPP));
179  i = w + P/w - aBy3;
180  }
181 
182  return vector(-VGREAT, i, VGREAT);
183  }
184  }
185 
186  // Sort the eigenvalues into ascending order
187  if (i > ii)
188  {
189  Swap(i, ii);
190  }
191 
192  if (ii > iii)
193  {
194  Swap(ii, iii);
195  }
196 
197  if (i > ii)
198  {
199  Swap(i, ii);
200  }
201 
202  return vector(i, ii, iii);
203 }
204 
205 
207 (
208  const tensor& t,
209  const scalar lambda
210 )
211 {
212  // Constantly rotating direction ensures different eigenvectors are
213  // generated when called sequentially with a multiple eigenvalue
214  static vector direction(1,0,0);
215  vector oldDirection(direction);
216  scalar temp = direction[2];
217  direction[2] = direction[1];
218  direction[1] = direction[0];
219  direction[0] = temp;
220 
221  // Construct the linear system for this eigenvalue
222  tensor A(t - lambda*I);
223 
224  // Determinants of the 2x2 sub-matrices used to find the eigenvectors
225  scalar sd0, sd1, sd2;
226  scalar magSd0, magSd1, magSd2;
227 
228  // Sub-determinants for a unique eivenvalue
229  sd0 = A.yy()*A.zz() - A.yz()*A.zy();
230  sd1 = A.zz()*A.xx() - A.zx()*A.xz();
231  sd2 = A.xx()*A.yy() - A.xy()*A.yx();
232  magSd0 = mag(sd0);
233  magSd1 = mag(sd1);
234  magSd2 = mag(sd2);
235 
236  // Evaluate the eigenvector using the largest sub-determinant
237  if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
238  {
239  vector ev
240  (
241  1,
242  (A.yz()*A.zx() - A.zz()*A.yx())/sd0,
243  (A.zy()*A.yx() - A.yy()*A.zx())/sd0
244  );
245 
246  return ev/mag(ev);
247  }
248  else if (magSd1 >= magSd2 && magSd1 > SMALL)
249  {
250  vector ev
251  (
252  (A.xz()*A.zy() - A.zz()*A.xy())/sd1,
253  1,
254  (A.zx()*A.xy() - A.xx()*A.zy())/sd1
255  );
256 
257  return ev/mag(ev);
258  }
259  else if (magSd2 > SMALL)
260  {
261  vector ev
262  (
263  (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
264  (A.yx()*A.xz() - A.xx()*A.yz())/sd2,
265  1
266  );
267 
268  return ev/mag(ev);
269  }
270 
271  // Sub-determinants for a repeated eigenvalue
272  sd0 = A.yy()*direction.z() - A.yz()*direction.y();
273  sd1 = A.zz()*direction.x() - A.zx()*direction.z();
274  sd2 = A.xx()*direction.y() - A.xy()*direction.x();
275  magSd0 = mag(sd0);
276  magSd1 = mag(sd1);
277  magSd2 = mag(sd2);
278 
279  // Evaluate the eigenvector using the largest sub-determinant
280  if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
281  {
282  vector ev
283  (
284  1,
285  (A.yz()*direction.x() - direction.z()*A.yx())/sd0,
286  (direction.y()*A.yx() - A.yy()*direction.x())/sd0
287  );
288 
289  return ev/mag(ev);
290  }
291  else if (magSd1 >= magSd2 && magSd1 > SMALL)
292  {
293  vector ev
294  (
295  (direction.z()*A.zy() - A.zz()*direction.y())/sd1,
296  1,
297  (A.zx()*direction.y() - direction.x()*A.zy())/sd1
298  );
299 
300  return ev/mag(ev);
301  }
302  else if (magSd2 > SMALL)
303  {
304  vector ev
305  (
306  (A.xy()*direction.z() - direction.y()*A.xz())/sd2,
307  (direction.x()*A.xz() - A.xx()*direction.z())/sd2,
308  1
309  );
310 
311  return ev/mag(ev);
312  }
313 
314  // Triple eigenvalue
315  return oldDirection;
316 }
317 
318 
320 {
321  vector evals(eigenValues(t));
322 
323  tensor evs
324  (
325  eigenVector(t, evals.x()),
326  eigenVector(t, evals.y()),
327  eigenVector(t, evals.z())
328  );
329 
330  return evs;
331 }
332 
333 
335 {
336  return eigenValues(tensor(t));
337 }
338 
339 
340 Foam::vector Foam::eigenVector(const symmTensor& t, const scalar lambda)
341 {
342  return eigenVector(tensor(t), lambda);
343 }
344 
345 
347 {
348  return eigenVectors(tensor(t));
349 }
350 
351 
352 // ************************************************************************* //
static const char *const typeName
Definition: Tensor.H:77
dimensionedScalar sqrt(const dimensionedScalar &ds)
const Cmpt & xy() const
Definition: TensorI.H:167
unsigned char direction
Definition: direction.H:43
dimensionedVector eigenValues(const dimensionedTensor &dt)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
const Cmpt & zz() const
Definition: TensorI.H:216
dimensioned< scalar > mag(const dimensioned< Type > &)
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
Templated 3D symmetric tensor derived from VectorSpace adding construction from 6 components...
Definition: SymmTensor.H:53
static const Tensor zero
Definition: Tensor.H:80
static const sphericalTensor I(1)
dimensionedScalar cbrt(const dimensionedScalar &ds)
const Cmpt & yy() const
Definition: TensorI.H:188
dimensionedScalar sign(const dimensionedScalar &ds)
static const Tensor min
Definition: Tensor.H:83
dimensionedTensor eigenVectors(const dimensionedTensor &dt)
Namespace for OpenFOAM.
static const char * componentNames[]
Definition: Tensor.H:78
const Cmpt & y() const
Definition: VectorI.H:71
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define WarningIn(functionName)
Report a warning using Foam::Warning.
dimensionedScalar cos(const dimensionedScalar &ds)
const Cmpt & x() const
Definition: VectorI.H:65
const Cmpt & yx() const
Definition: TensorI.H:181
scalar delta
dimensionedScalar acos(const dimensionedScalar &ds)
const Cmpt & z() const
Definition: VectorI.H:77
void Swap(T &a, T &b)
Definition: Swap.H:43
const Cmpt & zy() const
Definition: TensorI.H:209
const dimensionedScalar c
Speed of light in a vacuum.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
vector eigenVector(const tensor &, const scalar lambda)
Definition: tensor.C:207
static const Tensor max
Definition: Tensor.H:82
static const Tensor I
Definition: Tensor.H:84
const Cmpt & xx() const
Definition: TensorI.H:160
const Cmpt & yz() const
Definition: TensorI.H:195
mathematical constants.
const Cmpt & xz() const
Definition: TensorI.H:174
This function object calculates and outputs the second invariant of the velocity gradient tensor [1/s...
Definition: Q.H:66
const Cmpt & zx() const
Definition: TensorI.H:202
static const Tensor one
Definition: Tensor.H:81