100 if ((
sqr(
T.xy()) +
sqr(
T.xz()) +
sqr(
T.yz())) < ROOTSMALL)
110 const scalar
b = -
T.xx() -
T.yy() -
T.zz();
112 T.xx()*
T.yy() +
T.xx()*
T.zz() +
T.yy()*
T.zz()
113 -
T.xy()*
T.yx() -
T.yz()*
T.zy() -
T.zx()*
T.xz();
115 -
T.xx()*
T.yy()*
T.zz()
116 -
T.xy()*
T.yz()*
T.zx() -
T.xz()*
T.zy()*
T.yx()
117 +
T.xx()*
T.yz()*
T.zy() +
T.yy()*
T.zx()*
T.xz() +
T.zz()*
T.xy()*
T.yx();
120 Roots<3> roots(cubicEqn(1,
b,
c, d).roots());
127 switch (roots.type(i))
137 <<
"Eigenvalue computation fails for symmTensor: " <<
T 138 <<
"due to the non-real root = " << roots[i]
156 const vector& standardBasis1,
157 const vector& standardBasis2
166 const scalar sd0 =
A.yy()*
A.zz() -
A.yz()*
A.zy();
167 const scalar sd1 =
A.zz()*
A.xx() -
A.zx()*
A.xz();
168 const scalar sd2 =
A.xx()*
A.yy() -
A.xy()*
A.yx();
169 const scalar magSd0 =
mag(sd0);
170 const scalar magSd1 =
mag(sd1);
171 const scalar magSd2 =
mag(sd2);
174 if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
179 (
A.yz()*
A.zx() -
A.zz()*
A.yx())/sd0,
180 (
A.zy()*
A.yx() -
A.yy()*
A.zx())/sd0
184 if (
mag(eVec) < SMALL)
187 <<
"Eigenvector magnitude should be non-zero:" 188 <<
"mag(eigenvector) = " <<
mag(eVec)
193 return eVec/
mag(eVec);
195 else if (magSd1 >= magSd2 && magSd1 > SMALL)
199 (
A.xz()*
A.zy() -
A.zz()*
A.xy())/sd1,
201 (
A.zx()*
A.xy() -
A.xx()*
A.zy())/sd1
205 if (
mag(eVec) < SMALL)
208 <<
"Eigenvector magnitude should be non-zero:" 209 <<
"mag(eigenvector) = " <<
mag(eVec)
214 return eVec/
mag(eVec);
216 else if (magSd2 > SMALL)
220 (
A.xy()*
A.yz() -
A.yy()*
A.xz())/sd2,
221 (
A.yx()*
A.xz() -
A.xx()*
A.yz())/sd2,
226 if (
mag(eVec) < SMALL)
229 <<
"Eigenvector magnitude should be non-zero:" 230 <<
"mag(eigenvector) = " <<
mag(eVec)
235 return eVec/
mag(eVec);
240 const scalar sd0 =
A.yy()*standardBasis1.z() -
A.yz()*standardBasis1.y();
241 const scalar sd1 =
A.zz()*standardBasis1.x() -
A.zx()*standardBasis1.z();
242 const scalar sd2 =
A.xx()*standardBasis1.y() -
A.xy()*standardBasis1.x();
243 const scalar magSd0 =
mag(sd0);
244 const scalar magSd1 =
mag(sd1);
245 const scalar magSd2 =
mag(sd2);
248 if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
253 (
A.yz()*standardBasis1.x() - standardBasis1.z()*
A.yx())/sd0,
254 (standardBasis1.y()*
A.yx() -
A.yy()*standardBasis1.x())/sd0
258 if (
mag(eVec) < SMALL)
261 <<
"Eigenvector magnitude should be non-zero:" 262 <<
"mag(eigenvector) = " <<
mag(eVec)
267 return eVec/
mag(eVec);
269 else if (magSd1 >= magSd2 && magSd1 > SMALL)
273 (standardBasis1.z()*
A.zy() -
A.zz()*standardBasis1.y())/sd1,
275 (
A.zx()*standardBasis1.y() - standardBasis1.x()*
A.zy())/sd1
279 if (
mag(eVec) < SMALL)
282 <<
"Eigenvector magnitude should be non-zero:" 283 <<
"mag(eigenvector) = " <<
mag(eVec)
288 return eVec/
mag(eVec);
290 else if (magSd2 > SMALL)
294 (
A.xy()*standardBasis1.z() - standardBasis1.y()*
A.xz())/sd2,
295 (standardBasis1.x()*
A.xz() -
A.xx()*standardBasis1.z())/sd2,
300 if (
mag(eVec) < SMALL)
303 <<
"Eigenvector magnitude should be non-zero:" 304 <<
"mag(eigenvector) = " <<
mag(eVec)
309 return eVec/
mag(eVec);
313 return standardBasis1^standardBasis2;
323 vector Ux(1, 0, 0), Uy(0, 1, 0), Uz(0, 0, 1);
329 return tensor(Ux, Uy, Uz);
343 const scalar dt =
det(st);
350 mat(0,0) = st.xx(); mat(0,1) = st.xy(); mat(0,2) = st.xz();
351 mat(1,0) = st.yx(); mat(1,1) = st.yy(); mat(1,2) = st.yz();
352 mat(2,0) = st.zx(); mat(2,1) = st.zy(); mat(2,2) = st.zz();
358 mat(0,0), mat(0,1), mat(0,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.
static const char *const typeName
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)
static const char *const componentNames[]
static const Form rootMin
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
#define forAll(list, i)
Loop across all elements in list.
static const SymmTensor I
static const Identity< scalar > I
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
SymmTensor< scalar > symmTensor
SymmTensor of scalars, i.e. SymmTensor<scalar>.
void sort(UList< T > &list)
Sort the list.
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)
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
static const Form rootMax
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar c
Speed of light in a vacuum.
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
static Form uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
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)