56 const scalar
b = -
T.xx() -
T.yy() -
T.zz();
58 T.xx()*
T.yy() +
T.xx()*
T.zz() +
T.yy()*
T.zz()
59 -
T.xy()*
T.yx() -
T.yz()*
T.zy() -
T.zx()*
T.xz();
61 -
T.xx()*
T.yy()*
T.zz()
62 -
T.xy()*
T.yz()*
T.zx() -
T.xz()*
T.zy()*
T.yx()
63 +
T.xx()*
T.yz()*
T.zy() +
T.yy()*
T.zx()*
T.xz() +
T.zz()*
T.xy()*
T.yx();
68 bool isComplex =
false;
71 switch (roots.
type(i))
80 <<
"Eigenvalue computation fails for tensor: " <<
T 81 <<
"due to the not-a-number root = " << roots[i]
113 const Vector<complex>& standardBasis1,
114 const Vector<complex>& standardBasis2
118 Tensor<complex>
A(
Zero);
132 scalar magSd0 =
mag(sd0);
133 scalar magSd1 =
mag(sd1);
134 scalar magSd2 =
mag(sd2);
137 if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
139 const Vector<complex> eVec
142 (
A.yz()*
A.zx() -
A.zz()*
A.yx())/sd0,
143 (
A.zy()*
A.yx() -
A.yy()*
A.zx())/sd0
147 if (
mag(eVec) < SMALL)
150 <<
"Eigenvector magnitude should be non-zero:" 151 <<
"mag(eigenvector) = " <<
mag(eVec)
156 return eVec/
mag(eVec);
158 else if (magSd1 >= magSd2 && magSd1 > SMALL)
160 const Vector<complex> eVec
162 (
A.xz()*
A.zy() -
A.zz()*
A.xy())/sd1,
164 (
A.zx()*
A.xy() -
A.xx()*
A.zy())/sd1
168 if (
mag(eVec) < SMALL)
171 <<
"Eigenvector magnitude should be non-zero:" 172 <<
"mag(eigenvector) = " <<
mag(eVec)
177 return eVec/
mag(eVec);
179 else if (magSd2 > SMALL)
181 const Vector<complex> eVec
183 (
A.xy()*
A.yz() -
A.yy()*
A.xz())/sd2,
184 (
A.yx()*
A.xz() -
A.xx()*
A.yz())/sd2,
189 if (
mag(eVec) < SMALL)
192 <<
"Eigenvector magnitude should be non-zero:" 193 <<
"mag(eigenvector) = " <<
mag(eVec)
198 return eVec/
mag(eVec);
202 sd0 =
A.yy()*standardBasis1.z() -
A.yz()*standardBasis1.y();
203 sd1 =
A.zz()*standardBasis1.x() -
A.zx()*standardBasis1.z();
204 sd2 =
A.xx()*standardBasis1.y() -
A.xy()*standardBasis1.x();
210 if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
212 const Vector<complex> eVec
215 (
A.yz()*standardBasis1.x() - standardBasis1.z()*
A.yx())/sd0,
216 (standardBasis1.y()*
A.yx() -
A.yy()*standardBasis1.x())/sd0
220 if (
mag(eVec) < SMALL)
223 <<
"Eigenvector magnitude should be non-zero:" 224 <<
"mag(eigenvector) = " <<
mag(eVec)
229 return eVec/
mag(eVec);
231 else if (magSd1 >= magSd2 && magSd1 > SMALL)
233 const Vector<complex> eVec
235 (standardBasis1.z()*
A.zy() -
A.zz()*standardBasis1.y())/sd1,
237 (
A.zx()*standardBasis1.y() - standardBasis1.x()*
A.zy())/sd1
241 if (
mag(eVec) < SMALL)
244 <<
"Eigenvector magnitude should be non-zero:" 245 <<
"mag(eigenvector) = " <<
mag(eVec)
250 return eVec/
mag(eVec);
252 else if (magSd2 > SMALL)
254 const Vector<complex> eVec
256 (
A.xy()*standardBasis1.z() - standardBasis1.y()*
A.xz())/sd2,
257 (standardBasis1.x()*
A.xz() -
A.xx()*standardBasis1.z())/sd2,
262 if (
mag(eVec) < SMALL)
265 <<
"Eigenvector magnitude should be non-zero:" 266 <<
"mag(eigenvector) = " <<
mag(eVec)
271 return eVec/
mag(eVec);
275 return standardBasis1^standardBasis2;
282 const Vector<complex>& eVals
307 const scalar detval = t.
det();
309 if (detval < ROOTVSMALL)
314 mat(0,0) = t.xx(); mat(0,1) = t.xy(); mat(0,2) = t.xz();
315 mat(1,0) = t.yx(); mat(1,1) = t.yy(); mat(1,2) = t.yz();
316 mat(2,0) = t.zx(); mat(2,1) = t.zy(); mat(2,2) = t.zz();
322 mat(0,0), mat(0,1), mat(0,2),
323 mat(1,0), mat(1,1), mat(1,2),
324 mat(2,0), mat(2,1), mat(2,2)
vector eigenVector(const symmTensor &T, const scalar eVal, const vector &standardBasis1, const vector &standardBasis2)
Return a real eigenvector corresponding to a given real eigenvalue of a given symmTensor.
scalarRectangularMatrix SVDinv(const scalarRectangularMatrix &A, scalar minCondition=0)
Return the inverse of matrix A using SVD.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
Cmpt det() const
The determinate.
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
#define forAll(list, i)
Loop across all elements in list.
Container to encapsulate various operations for cubic equation of the forms with real coefficients: ...
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
symmTensor pinv(const symmTensor &st)
Return inverse of a given symmTensor, and fall back to pseudo-inverse if the symmTensor is singular...
errorManip< error > abort(error &err)
Templated 3D Vector derived from VectorSpace adding construction from 3 components, element access using x(), y() and z() member functions and the inner-product (dot-product) and cross-product operators.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
#define WarningInFunction
Report a warning using Foam::Warning.
void type(const direction i, const roots::type t)
Set the type of the i-th root.
const dimensionedScalar c
Speed of light in a vacuum.
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
A templated (3 x 3) tensor of objects of <T> derived from MatrixSpace.
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Tensor of scalars, i.e. Tensor<scalar>.
RectangularMatrix< scalar > scalarRectangularMatrix
static constexpr const zero Zero
Global zero (0)