tensor.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) 2011-2018 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 "cubicEqn.H"
28 #include "mathematicalConstants.H"
29 
30 using namespace Foam::constant::mathematical;
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 template<>
35 const char* const Foam::tensor::vsType::typeName = "tensor";
36 
37 template<>
38 const char* const Foam::tensor::vsType::componentNames[] =
39 {
40  "xx", "xy", "xz",
41  "yx", "yy", "yz",
42  "zx", "zy", "zz"
43 };
44 
45 template<>
47 
48 template<>
50 
51 template<>
53 
54 template<>
56 
57 template<>
59 
60 template<>
62 
63 template<>
65 (
66  1, 0, 0,
67  0, 1, 0,
68  0, 0, 1
69 );
70 
71 
72 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
73 
75 {
76  // Coefficients of the characteristic cubic polynomial (a = 1)
77  const scalar b =
78  - t.xx() - t.yy() - t.zz();
79  const scalar c =
80  t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
81  - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz();
82  const scalar d =
83  - t.xx()*t.yy()*t.zz()
84  - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx()
85  + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx();
86 
87  // Solve
88  Roots<3> roots = cubicEqn(1, b, c, d).roots();
89 
90  // Check the root types
91  vector lambda = vector::zero;
92  forAll(roots, i)
93  {
94  switch (roots.type(i))
95  {
96  case rootType::real:
97  lambda[i] = roots[i];
98  break;
99  case rootType::complex:
101  << "Complex eigenvalues detected for tensor: " << t
102  << endl;
103  lambda[i] = 0;
104  break;
105  case rootType::posInf:
106  lambda[i] = vGreat;
107  break;
108  case rootType::negInf:
109  lambda[i] = - vGreat;
110  break;
111  case rootType::nan:
113  << "Eigenvalue calculation failed for tensor: " << t
114  << exit(FatalError);
115  }
116  }
117 
118  // Sort the eigenvalues into ascending order
119  if (lambda.x() > lambda.y())
120  {
121  Swap(lambda.x(), lambda.y());
122  }
123  if (lambda.y() > lambda.z())
124  {
125  Swap(lambda.y(), lambda.z());
126  }
127  if (lambda.x() > lambda.y())
128  {
129  Swap(lambda.x(), lambda.y());
130  }
131 
132  return lambda;
133 }
134 
135 
137 (
138  const tensor& T,
139  const scalar lambda,
140  const vector& direction1,
141  const vector& direction2
142 )
143 {
144  // Construct the linear system for this eigenvalue
145  tensor A(T - lambda*I);
146 
147  // Determinants of the 2x2 sub-matrices used to find the eigenvectors
148  scalar sd0, sd1, sd2;
149  scalar magSd0, magSd1, magSd2;
150 
151  // Sub-determinants for a unique eivenvalue
152  sd0 = A.yy()*A.zz() - A.yz()*A.zy();
153  sd1 = A.zz()*A.xx() - A.zx()*A.xz();
154  sd2 = A.xx()*A.yy() - A.xy()*A.yx();
155  magSd0 = mag(sd0);
156  magSd1 = mag(sd1);
157  magSd2 = mag(sd2);
158 
159  // Evaluate the eigenvector using the largest sub-determinant
160  if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > small)
161  {
162  vector ev
163  (
164  1,
165  (A.yz()*A.zx() - A.zz()*A.yx())/sd0,
166  (A.zy()*A.yx() - A.yy()*A.zx())/sd0
167  );
168 
169  return ev/mag(ev);
170  }
171  else if (magSd1 >= magSd2 && magSd1 > small)
172  {
173  vector ev
174  (
175  (A.xz()*A.zy() - A.zz()*A.xy())/sd1,
176  1,
177  (A.zx()*A.xy() - A.xx()*A.zy())/sd1
178  );
179 
180  return ev/mag(ev);
181  }
182  else if (magSd2 > small)
183  {
184  vector ev
185  (
186  (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
187  (A.yx()*A.xz() - A.xx()*A.yz())/sd2,
188  1
189  );
190 
191  return ev/mag(ev);
192  }
193 
194  // Sub-determinants for a repeated eigenvalue
195  sd0 = A.yy()*direction1.z() - A.yz()*direction1.y();
196  sd1 = A.zz()*direction1.x() - A.zx()*direction1.z();
197  sd2 = A.xx()*direction1.y() - A.xy()*direction1.x();
198  magSd0 = mag(sd0);
199  magSd1 = mag(sd1);
200  magSd2 = mag(sd2);
201 
202  // Evaluate the eigenvector using the largest sub-determinant
203  if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > small)
204  {
205  vector ev
206  (
207  1,
208  (A.yz()*direction1.x() - direction1.z()*A.yx())/sd0,
209  (direction1.y()*A.yx() - A.yy()*direction1.x())/sd0
210  );
211 
212  return ev/mag(ev);
213  }
214  else if (magSd1 >= magSd2 && magSd1 > small)
215  {
216  vector ev
217  (
218  (direction1.z()*A.zy() - A.zz()*direction1.y())/sd1,
219  1,
220  (A.zx()*direction1.y() - direction1.x()*A.zy())/sd1
221  );
222 
223  return ev/mag(ev);
224  }
225  else if (magSd2 > small)
226  {
227  vector ev
228  (
229  (A.xy()*direction1.z() - direction1.y()*A.xz())/sd2,
230  (direction1.x()*A.xz() - A.xx()*direction1.z())/sd2,
231  1
232  );
233 
234  return ev/mag(ev);
235  }
236 
237  // Triple eigenvalue
238  return direction1^direction2;
239 }
240 
241 
242 Foam::tensor Foam::eigenVectors(const tensor& T, const vector& lambdas)
243 {
244  vector Ux(1, 0, 0), Uy(0, 1, 0), Uz(0, 0, 1);
245 
246  Ux = eigenVector(T, lambdas.x(), Uy, Uz);
247  Uy = eigenVector(T, lambdas.y(), Uz, Ux);
248  Uz = eigenVector(T, lambdas.z(), Ux, Uy);
249 
250  return tensor(Ux, Uy, Uz);
251 }
252 
253 
255 {
256  const vector lambdas(eigenValues(T));
257 
258  return eigenVectors(T, lambdas);
259 }
260 
261 
263 {
264  return eigenValues(tensor(T));
265 }
266 
267 
269 (
270  const symmTensor& T,
271  const scalar lambda,
272  const vector& direction1,
273  const vector& direction2
274 )
275 {
276  return eigenVector(tensor(T), lambda, direction1, direction2);
277 }
278 
279 
281 {
282  return eigenVectors(tensor(T), lambdas);
283 }
284 
285 
287 {
288  return eigenVectors(tensor(T));
289 }
290 
291 
292 // ************************************************************************* //
vector eigenVector(const tensor &T, const scalar lambda, const vector &direction1, const vector &direction2)
Definition: tensor.C:137
const Cmpt & yx() const
Definition: TensorI.H:174
static const char *const typeName
Definition: VectorSpace.H:111
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const Cmpt & zy() const
Definition: TensorI.H:202
const Cmpt & xz() const
Definition: TensorI.H:167
static const Form max
Definition: VectorSpace.H:115
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
static const char *const componentNames[]
Definition: VectorSpace.H:112
const Cmpt & yy() const
Definition: TensorI.H:181
const Cmpt & xx() const
Definition: TensorI.H:153
dimensionedVector eigenValues(const dimensionedTensor &dt)
const Cmpt & yz() const
Definition: TensorI.H:188
static const Form rootMin
Definition: VectorSpace.H:118
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
const Cmpt & z() const
Definition: VectorI.H:87
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
Definition: Roots.H:64
const Cmpt & y() const
Definition: VectorI.H:81
static const Form min
Definition: VectorSpace.H:116
Cubic equation of the form a*x^3 + b*x^2 + c*x + d = 0.
Definition: cubicEqn.H:49
void Swap(T &a, T &b)
Definition: Swap.H:43
static const Identity< scalar > I
Definition: Identity.H:93
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
mathematical constants.
const Cmpt & x() const
Definition: VectorI.H:75
static const Tensor I
Definition: Tensor.H:80
dimensionedTensor eigenVectors(const dimensionedTensor &dt)
const Cmpt & zx() const
Definition: TensorI.H:195
Roots< 3 > roots() const
Get the roots.
Definition: cubicEqn.C:32
static const Form rootMax
Definition: VectorSpace.H:117
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar c
Speed of light in a vacuum.
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
static const Form one
Definition: VectorSpace.H:114
static Tensor< Cmpt > uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
Definition: VectorSpaceI.H:168
dimensioned< scalar > mag(const dimensioned< Type > &)
const Cmpt & xy() const
Definition: TensorI.H:160
const Cmpt & zz() const
Definition: TensorI.H:209
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
static const Form zero
Definition: VectorSpace.H:113
void type(const direction i, const rootType t)
Set the type of the i-th root.
Definition: RootsI.H:120