Go to the documentation of this file.
99 if ((
sqr(
T.xy()) +
sqr(
T.xz()) +
sqr(
T.yz())) < ROOTSMALL)
109 const scalar
b = -
T.xx() -
T.yy() -
T.zz();
111 T.xx()*
T.yy() +
T.xx()*
T.zz() +
T.yy()*
T.zz()
112 -
T.xy()*
T.yx() -
T.yz()*
T.zy() -
T.zx()*
T.xz();
114 -
T.xx()*
T.yy()*
T.zz()
115 -
T.xy()*
T.yz()*
T.zx() -
T.xz()*
T.zy()*
T.yx()
116 +
T.xx()*
T.yz()*
T.zy() +
T.yy()*
T.zx()*
T.xz() +
T.zz()*
T.xy()*
T.yx();
126 switch (roots.
type(i))
136 <<
"Eigenvalue computation fails for symmTensor: " <<
T
137 <<
"due to the non-real root = " << roots[i]
155 const vector& standardBasis1,
156 const vector& standardBasis2
165 const scalar sd0 =
A.yy()*
A.zz() -
A.yz()*
A.zy();
166 const scalar sd1 =
A.zz()*
A.xx() -
A.zx()*
A.xz();
167 const scalar sd2 =
A.xx()*
A.yy() -
A.xy()*
A.yx();
168 const scalar magSd0 =
mag(sd0);
169 const scalar magSd1 =
mag(sd1);
170 const scalar magSd2 =
mag(sd2);
173 if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
178 (
A.yz()*
A.zx() -
A.zz()*
A.yx())/sd0,
179 (
A.zy()*
A.yx() -
A.yy()*
A.zx())/sd0
183 if (
mag(eVec) < SMALL)
186 <<
"Eigenvector magnitude should be non-zero:"
187 <<
"mag(eigenvector) = " <<
mag(eVec)
192 return eVec/
mag(eVec);
194 else if (magSd1 >= magSd2 && magSd1 > SMALL)
198 (
A.xz()*
A.zy() -
A.zz()*
A.xy())/sd1,
200 (
A.zx()*
A.xy() -
A.xx()*
A.zy())/sd1
204 if (
mag(eVec) < SMALL)
207 <<
"Eigenvector magnitude should be non-zero:"
208 <<
"mag(eigenvector) = " <<
mag(eVec)
213 return eVec/
mag(eVec);
215 else if (magSd2 > SMALL)
219 (
A.xy()*
A.yz() -
A.yy()*
A.xz())/sd2,
220 (
A.yx()*
A.xz() -
A.xx()*
A.yz())/sd2,
225 if (
mag(eVec) < SMALL)
228 <<
"Eigenvector magnitude should be non-zero:"
229 <<
"mag(eigenvector) = " <<
mag(eVec)
234 return eVec/
mag(eVec);
239 const scalar sd0 =
A.yy()*standardBasis1.
z() -
A.yz()*standardBasis1.
y();
240 const scalar sd1 =
A.zz()*standardBasis1.
x() -
A.zx()*standardBasis1.
z();
241 const scalar sd2 =
A.xx()*standardBasis1.
y() -
A.xy()*standardBasis1.
x();
242 const scalar magSd0 =
mag(sd0);
243 const scalar magSd1 =
mag(sd1);
244 const scalar magSd2 =
mag(sd2);
247 if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
252 (
A.yz()*standardBasis1.
x() - standardBasis1.
z()*
A.yx())/sd0,
253 (standardBasis1.
y()*
A.yx() -
A.yy()*standardBasis1.
x())/sd0
257 if (
mag(eVec) < SMALL)
260 <<
"Eigenvector magnitude should be non-zero:"
261 <<
"mag(eigenvector) = " <<
mag(eVec)
266 return eVec/
mag(eVec);
268 else if (magSd1 >= magSd2 && magSd1 > SMALL)
272 (standardBasis1.
z()*
A.zy() -
A.zz()*standardBasis1.
y())/sd1,
274 (
A.zx()*standardBasis1.
y() - standardBasis1.
x()*
A.zy())/sd1
278 if (
mag(eVec) < SMALL)
281 <<
"Eigenvector magnitude should be non-zero:"
282 <<
"mag(eigenvector) = " <<
mag(eVec)
287 return eVec/
mag(eVec);
289 else if (magSd2 > SMALL)
293 (
A.xy()*standardBasis1.
z() - standardBasis1.
y()*
A.xz())/sd2,
294 (standardBasis1.
x()*
A.xz() -
A.xx()*standardBasis1.
z())/sd2,
299 if (
mag(eVec) < SMALL)
302 <<
"Eigenvector magnitude should be non-zero:"
303 <<
"mag(eigenvector) = " <<
mag(eVec)
308 return eVec/
mag(eVec);
312 return standardBasis1^standardBasis2;
322 vector Ux(1, 0, 0), Uy(0, 1, 0), Uz(0, 0, 1);
328 return tensor(Ux, Uy, Uz);
static const Form rootMax
void type(const direction i, const roots::type t)
Set the type of the i-th root.
const Cmpt & x() const
Access to the vector x component.
A templated (3 x 3) symmetric tensor of objects of <T>, effectively containing 6 elements,...
static SymmTensor< scalar > uniform(const scalar &s)
Return a VectorSpace with all elements = s.
static const char *const componentNames[]
static constexpr const zero Zero
Global zero (0)
vector eigenVector(const symmTensor &T, const scalar eVal, const vector &standardBasis1, const vector &standardBasis2)
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
iterator end() noexcept
Return an iterator to end of VectorSpace.
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
iterator begin() noexcept
Return an iterator to begin of VectorSpace.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const Cmpt & z() const
Access to the vector z component.
#define forAll(list, i)
Loop across all elements in list.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
static const SymmTensor I
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
errorManip< error > abort(error &err)
static const Form rootMin
const Cmpt & y() const
Access to the vector y component.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
const dimensionedScalar c
Speed of light in a vacuum.
Container to encapsulate various operations for cubic equation of the forms with real coefficients:
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
#define WarningInFunction
Report a warning using Foam::Warning.
Tensor< scalar > tensor
Tensor of scalars, i.e. Tensor<scalar>.
static const Identity< scalar > I
static const char *const typeName