99 if ((
sqr(
T.xy()) +
sqr(
T.xz()) +
sqr(
T.yz())) < ROOTSMALL)
103 std::sort(eVals.
begin(), eVals.
end());
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]
145 std::sort(eVals.
begin(), eVals.
end());
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 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.
static const SymmTensor I
iterator end() noexcept
Return an iterator to end of VectorSpace.
iterator begin() noexcept
Return an iterator to begin of VectorSpace.
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.
Container to encapsulate various operations for cubic equation of the forms with real coefficients:
A class representing the concept of 1 (one) that can be used to avoid manipulating objects known to b...
static const char *const componentNames[]
static const complex rootMax
complex (ROOTVGREAT, ROOTVGREAT)
static const complex min
complex (-VGREAT,-VGREAT)
static const complex max
complex (VGREAT,VGREAT)
static const complex rootMin
complex (-ROOTVGREAT, -ROOTVGREAT)
#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.
static const char *const typeName
The type name used in ensight case files.
A non-counting (dummy) refCount.