55 const scalar
b = -
T.xx() -
T.yy() -
T.zz();
57 T.xx()*
T.yy() +
T.xx()*
T.zz() +
T.yy()*
T.zz()
58 -
T.xy()*
T.yx() -
T.yz()*
T.zy() -
T.zx()*
T.xz();
60 -
T.xx()*
T.yy()*
T.zz()
61 -
T.xy()*
T.yz()*
T.zx() -
T.xz()*
T.zy()*
T.yx()
62 +
T.xx()*
T.yz()*
T.zy() +
T.yy()*
T.zx()*
T.xz() +
T.zz()*
T.xy()*
T.yx();
67 bool isComplex =
false;
70 switch (roots.
type(i))
79 <<
"Eigenvalue computation fails for tensor: " <<
T
80 <<
"due to the not-a-number root = " << roots[i]
131 scalar magSd0 =
mag(sd0);
132 scalar magSd1 =
mag(sd1);
133 scalar magSd2 =
mag(sd2);
136 if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
141 (
A.yz()*
A.zx() -
A.zz()*
A.yx())/sd0,
142 (
A.zy()*
A.yx() -
A.yy()*
A.zx())/sd0
146 if (
mag(eVec) < SMALL)
149 <<
"Eigenvector magnitude should be non-zero:"
150 <<
"mag(eigenvector) = " <<
mag(eVec)
155 return eVec/
mag(eVec);
157 else if (magSd1 >= magSd2 && magSd1 > SMALL)
161 (
A.xz()*
A.zy() -
A.zz()*
A.xy())/sd1,
163 (
A.zx()*
A.xy() -
A.xx()*
A.zy())/sd1
167 if (
mag(eVec) < SMALL)
170 <<
"Eigenvector magnitude should be non-zero:"
171 <<
"mag(eigenvector) = " <<
mag(eVec)
176 return eVec/
mag(eVec);
178 else if (magSd2 > SMALL)
182 (
A.xy()*
A.yz() -
A.yy()*
A.xz())/sd2,
183 (
A.yx()*
A.xz() -
A.xx()*
A.yz())/sd2,
188 if (
mag(eVec) < SMALL)
191 <<
"Eigenvector magnitude should be non-zero:"
192 <<
"mag(eigenvector) = " <<
mag(eVec)
197 return eVec/
mag(eVec);
201 sd0 =
A.yy()*standardBasis1.
z() -
A.yz()*standardBasis1.
y();
202 sd1 =
A.zz()*standardBasis1.
x() -
A.zx()*standardBasis1.
z();
203 sd2 =
A.xx()*standardBasis1.
y() -
A.xy()*standardBasis1.
x();
209 if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
214 (
A.yz()*standardBasis1.
x() - standardBasis1.
z()*
A.yx())/sd0,
215 (standardBasis1.
y()*
A.yx() -
A.yy()*standardBasis1.
x())/sd0
219 if (
mag(eVec) < SMALL)
222 <<
"Eigenvector magnitude should be non-zero:"
223 <<
"mag(eigenvector) = " <<
mag(eVec)
228 return eVec/
mag(eVec);
230 else if (magSd1 >= magSd2 && magSd1 > SMALL)
234 (standardBasis1.
z()*
A.zy() -
A.zz()*standardBasis1.
y())/sd1,
236 (
A.zx()*standardBasis1.
y() - standardBasis1.
x()*
A.zy())/sd1
240 if (
mag(eVec) < SMALL)
243 <<
"Eigenvector magnitude should be non-zero:"
244 <<
"mag(eigenvector) = " <<
mag(eVec)
249 return eVec/
mag(eVec);
251 else if (magSd2 > SMALL)
255 (
A.xy()*standardBasis1.
z() - standardBasis1.
y()*
A.xz())/sd2,
256 (standardBasis1.
x()*
A.xz() -
A.xx()*standardBasis1.
z())/sd2,
261 if (
mag(eVec) < SMALL)
264 <<
"Eigenvector magnitude should be non-zero:"
265 <<
"mag(eigenvector) = " <<
mag(eVec)
270 return eVec/
mag(eVec);
274 return standardBasis1^standardBasis2;
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
void type(const direction i, const roots::type t)
Set the type of the i-th root.
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
const Cmpt & z() const
Access to the vector z component.
const Cmpt & y() const
Access to the vector y component.
const Cmpt & x() const
Access to the vector x component.
A complex number, similar to the C++ complex type.
Container to encapsulate various operations for cubic equation of the forms with real coefficients:
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
#define WarningInFunction
Report a warning using Foam::Warning.
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
vector eigenVector(const symmTensor &T, const scalar eVal, const vector &standardBasis1, const vector &standardBasis2)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
errorManip< error > abort(error &err)
#define forAll(list, i)
Loop across all elements in list.