59 tensor2D::uniform(VGREAT)
65 tensor2D::uniform(-VGREAT)
71 tensor2D::uniform(ROOTVGREAT)
77 tensor2D::uniform(-ROOTVGREAT)
92 const scalar a =
T.xx();
93 const scalar
b =
T.xy();
94 const scalar
c =
T.yx();
95 const scalar d =
T.yy();
98 if ((
sqr(
b) +
sqr(c)) < ROOTSMALL)
103 const scalar trace = a + d;
107 scalar
e = std::fma(-
b, c, w);
108 scalar
f = std::fma(a, d, -w);
109 const scalar determinant =
f +
e;
112 const scalar gapSqr = std::fma(-4.0, determinant,
sqr(trace));
120 if (
mag(firstRoot) < SMALL)
123 <<
"Zero-valued root is found. Adding SMALL to the root "
124 <<
"to avoid floating-point exception." <<
nl;
128 Vector2D<complex> eVals
131 complex(determinant/firstRoot, 0)
135 if (eVals.x().real() > eVals.y().real())
137 Swap(eVals.x(), eVals.y());
147 return Vector2D<complex>
176 if (
mag(eVec) < SMALL)
179 <<
"Eigenvector magnitude should be non-zero:"
180 <<
"mag(eigenvector) = " <<
mag(eVec)
185 return eVec/
mag(eVec);
187 else if (
mag(
A.xx()) > SMALL)
192 if (
mag(eVec) < SMALL)
195 <<
"Eigenvector magnitude should be non-zero:"
196 <<
"mag(eigenvector) = " <<
mag(eVec)
201 return eVec/
mag(eVec);
204 else if (
mag(
T.yx()) >
mag(
T.xy()) &&
mag(
T.yx()) > SMALL)
209 if (
mag(eVec) < SMALL)
212 <<
"Eigenvector magnitude should be non-zero:"
213 <<
"mag(eigenvector) = " <<
mag(eVec)
218 return eVec/
mag(eVec);
220 else if (
mag(
T.xy()) > SMALL)
225 if (
mag(eVec) < SMALL)
228 <<
"Eigenvector magnitude should be non-zero:"
229 <<
"mag(eigenvector) = " <<
mag(eVec)
234 return eVec/
mag(eVec);
260 const Vector2D<complex> eVals(eigenValues(
T));
262 return eigenVectors(
T, eVals);
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
A templated (2 x 2) tensor of objects of <T> derived from VectorSpace.
Templated 2D Vector derived from VectorSpace adding construction from 2 components,...
const Cmpt & y() const
Access to the vector y component.
const Cmpt & x() const
Access to the vector x component.
A complex number, similar to the C++ complex type.
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.
const dimensionedScalar c
Speed of light in a vacuum.
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
dimensionedScalar sign(const dimensionedScalar &ds)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
vector eigenVector(const symmTensor &T, const scalar eVal, const vector &standardBasis1, const vector &standardBasis2)
void Swap(DynamicList< T, SizeMinA > &a, DynamicList< T, SizeMinB > &b)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
errorManip< error > abort(error &err)
constexpr char nl
The newline '\n' character (0x0a)
static const char *const typeName
The type name used in ensight case files.
A non-counting (dummy) refCount.