Go to the documentation of this file.
38 const char*
const Foam::tensor::vsType::typeName =
"tensor";
41 const char*
const Foam::tensor::vsType::componentNames[] =
49 const Foam::tensor Foam::tensor::vsType::zero(tensor::uniform(0));
52 const Foam::tensor Foam::tensor::vsType::one(tensor::uniform(1));
61 const Foam::tensor Foam::tensor::vsType::rootMax(tensor::uniform(ROOTVGREAT));
64 const Foam::tensor Foam::tensor::vsType::rootMin(tensor::uniform(-ROOTVGREAT));
95 const scalar
b = -
T.xx() -
T.yy() -
T.zz();
97 T.xx()*
T.yy() +
T.xx()*
T.zz() +
T.yy()*
T.zz()
98 -
T.xy()*
T.yx() -
T.yz()*
T.zy() -
T.zx()*
T.xz();
100 -
T.xx()*
T.yy()*
T.zz()
101 -
T.xy()*
T.yz()*
T.zx() -
T.xz()*
T.zy()*
T.yx()
102 +
T.xx()*
T.yz()*
T.zy() +
T.yy()*
T.zx()*
T.xz() +
T.zz()*
T.xy()*
T.yx();
107 bool isComplex =
false;
110 switch (roots.
type(i))
119 <<
"Eigenvalue computation fails for tensor: " <<
T
120 <<
"due to the not-a-number root = " << roots[i]
171 scalar magSd0 =
mag(sd0);
172 scalar magSd1 =
mag(sd1);
173 scalar magSd2 =
mag(sd2);
176 if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
181 (
A.yz()*
A.zx() -
A.zz()*
A.yx())/sd0,
182 (
A.zy()*
A.yx() -
A.yy()*
A.zx())/sd0
186 if (
mag(eVec) < SMALL)
189 <<
"Eigenvector magnitude should be non-zero:"
190 <<
"mag(eigenvector) = " <<
mag(eVec)
195 return eVec/
mag(eVec);
197 else if (magSd1 >= magSd2 && magSd1 > SMALL)
201 (
A.xz()*
A.zy() -
A.zz()*
A.xy())/sd1,
203 (
A.zx()*
A.xy() -
A.xx()*
A.zy())/sd1
207 if (
mag(eVec) < SMALL)
210 <<
"Eigenvector magnitude should be non-zero:"
211 <<
"mag(eigenvector) = " <<
mag(eVec)
216 return eVec/
mag(eVec);
218 else if (magSd2 > SMALL)
222 (
A.xy()*
A.yz() -
A.yy()*
A.xz())/sd2,
223 (
A.yx()*
A.xz() -
A.xx()*
A.yz())/sd2,
228 if (
mag(eVec) < SMALL)
231 <<
"Eigenvector magnitude should be non-zero:"
232 <<
"mag(eigenvector) = " <<
mag(eVec)
237 return eVec/
mag(eVec);
241 sd0 =
A.yy()*standardBasis1.
z() -
A.yz()*standardBasis1.
y();
242 sd1 =
A.zz()*standardBasis1.
x() -
A.zx()*standardBasis1.
z();
243 sd2 =
A.xx()*standardBasis1.
y() -
A.xy()*standardBasis1.
x();
249 if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
254 (
A.yz()*standardBasis1.
x() - standardBasis1.
z()*
A.yx())/sd0,
255 (standardBasis1.
y()*
A.yx() -
A.yy()*standardBasis1.
x())/sd0
259 if (
mag(eVec) < SMALL)
262 <<
"Eigenvector magnitude should be non-zero:"
263 <<
"mag(eigenvector) = " <<
mag(eVec)
268 return eVec/
mag(eVec);
270 else if (magSd1 >= magSd2 && magSd1 > SMALL)
274 (standardBasis1.
z()*
A.zy() -
A.zz()*standardBasis1.
y())/sd1,
276 (
A.zx()*standardBasis1.
y() - standardBasis1.
x()*
A.zy())/sd1
280 if (
mag(eVec) < SMALL)
283 <<
"Eigenvector magnitude should be non-zero:"
284 <<
"mag(eigenvector) = " <<
mag(eVec)
289 return eVec/
mag(eVec);
291 else if (magSd2 > SMALL)
295 (
A.xy()*standardBasis1.
z() - standardBasis1.
y()*
A.xz())/sd2,
296 (standardBasis1.
x()*
A.xz() -
A.xx()*standardBasis1.
z())/sd2,
301 if (
mag(eVec) < SMALL)
304 <<
"Eigenvector magnitude should be non-zero:"
305 <<
"mag(eigenvector) = " <<
mag(eVec)
310 return eVec/
mag(eVec);
314 return standardBasis1^standardBasis2;
void type(const direction i, const roots::type t)
Set the type of the i-th root.
A templated (3 x 3) tensor of objects of <T> derived from MatrixSpace.
const Cmpt & x() const
Access to the vector x component.
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)
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Ostream & endl(Ostream &os)
Add newline and flush stream.
const Cmpt & z() const
Access to the vector z component.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
#define forAll(list, i)
Loop across all elements in list.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
A complex number, similar to the C++ complex type.
errorManip< error > abort(error &err)
const Cmpt & y() const
Access to the vector y component.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
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.