tensor2D.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "tensor2D.H"
30 #include "quadraticEqn.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 template<>
35 const char* const Foam::tensor2D::vsType::typeName = "tensor2D";
36 
37 template<>
38 const char* const Foam::tensor2D::vsType::componentNames[] =
39 {
40  "xx", "xy",
41  "yx", "yy"
42 };
43 
44 template<>
46 (
47  tensor2D::uniform(0)
48 );
49 
50 template<>
52 (
53  tensor2D::uniform(1)
54 );
55 
56 template<>
58 (
59  tensor2D::uniform(VGREAT)
60 );
61 
62 template<>
64 (
65  tensor2D::uniform(-VGREAT)
66 );
67 
68 template<>
70 (
71  tensor2D::uniform(ROOTVGREAT)
72 );
73 
74 template<>
76 (
77  tensor2D::uniform(-ROOTVGREAT)
78 );
79 
80 template<>
82 (
83  1, 0,
84  0, 1
85 );
86 
87 
88 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
89 
91 {
92  const scalar a = T.xx();
93  const scalar b = T.xy();
94  const scalar c = T.yx();
95  const scalar d = T.yy();
96 
97  // Return diagonal if T is effectively diagonal tensor
98  if ((sqr(b) + sqr(c)) < ROOTSMALL)
99  {
100  return Vector2D<complex>(complex(a), complex(d));
101  }
102 
103  const scalar trace = a + d;
104 
105  // (JLM:p. 2246)
106  scalar w = b*c;
107  scalar e = std::fma(-b, c, w);
108  scalar f = std::fma(a, d, -w);
109  const scalar determinant = f + e;
110 
111  // Square-distance between two eigenvalues
112  const scalar gapSqr = std::fma(-4.0, determinant, sqr(trace));
113 
114  // (F:Sec. 8.4.2.)
115  // Eigenvalues are effectively real
116  if (0 <= gapSqr)
117  {
118  scalar firstRoot = 0.5*(trace - sign(-trace)*Foam::sqrt(gapSqr));
119 
120  if (mag(firstRoot) < SMALL)
121  {
123  << "Zero-valued root is found. Adding SMALL to the root "
124  << "to avoid floating-point exception." << nl;
125  firstRoot = SMALL;
126  }
127 
128  Vector2D<complex> eVals
129  (
130  complex(firstRoot, 0),
131  complex(determinant/firstRoot, 0)
132  );
133 
134  // Sort the eigenvalues into ascending order
135  if (eVals.x().real() > eVals.y().real())
136  {
137  Swap(eVals.x(), eVals.y());
138  }
139 
140  return eVals;
141  }
142  // Eigenvalues are complex
143  else
144  {
145  const complex eVal(0.5*trace, 0.5*Foam::sqrt(mag(gapSqr)));
146 
147  return Vector2D<complex>
148  (
149  eVal,
150  eVal.conjugate()
151  );
152  }
153 }
154 
155 
157 (
158  const tensor2D& T,
159  const complex& eVal,
160  const Vector2D<complex>& standardBasis
161 )
162 {
163  // Construct the linear system for this eigenvalue
164  const Tensor2D<complex> A
165  (
166  complex(T.xx()) - eVal, complex(T.xy()),
167  complex(T.yx()), complex(T.yy()) - eVal
168  );
169 
170  // Evaluate the eigenvector using the largest divisor
171  if (mag(A.yy()) > mag(A.xx()) && mag(A.yy()) > SMALL)
172  {
173  Vector2D<complex> eVec(complex(1), -A.yx()/A.yy());
174 
175  #ifdef FULLDEBUG
176  if (mag(eVec) < SMALL)
177  {
179  << "Eigenvector magnitude should be non-zero:"
180  << "mag(eigenvector) = " << mag(eVec)
181  << abort(FatalError);
182  }
183  #endif
184 
185  return eVec/mag(eVec);
186  }
187  else if (mag(A.xx()) > SMALL)
188  {
189  Vector2D<complex> eVec(-A.xy()/A.xx(), complex(1));
190 
191  #ifdef FULLDEBUG
192  if (mag(eVec) < SMALL)
193  {
195  << "Eigenvector magnitude should be non-zero:"
196  << "mag(eigenvector) = " << mag(eVec)
197  << abort(FatalError);
198  }
199  #endif
200 
201  return eVec/mag(eVec);
202  }
203  // (K:p. 47-48)
204  else if (mag(T.yx()) > mag(T.xy()) && mag(T.yx()) > SMALL)
205  {
206  const Vector2D<complex> eVec(eVal - T.yy(), complex(T.yx()));
207 
208  #ifdef FULLDEBUG
209  if (mag(eVec) < SMALL)
210  {
212  << "Eigenvector magnitude should be non-zero:"
213  << "mag(eigenvector) = " << mag(eVec)
214  << abort(FatalError);
215  }
216  #endif
217 
218  return eVec/mag(eVec);
219  }
220  else if (mag(T.xy()) > SMALL)
221  {
222  const Vector2D<complex> eVec(complex(T.xy()), eVal - T.xx());
223 
224  #ifdef FULLDEBUG
225  if (mag(eVec) < SMALL)
226  {
228  << "Eigenvector magnitude should be non-zero:"
229  << "mag(eigenvector) = " << mag(eVec)
230  << abort(FatalError);
231  }
232  #endif
233 
234  return eVec/mag(eVec);
235  }
236 
237  // Repeated eigenvalue
238  return Vector2D<complex>(-standardBasis.y(), standardBasis.x());
239 }
240 
241 
243 (
244  const tensor2D& T,
245  const Vector2D<complex>& eVals
246 )
247 {
248  Vector2D<complex> Ux(pTraits<complex>::one, Zero);
249  Vector2D<complex> Uy(Zero, pTraits<complex>::one);
250 
251  Ux = eigenVector(T, eVals.x(), Uy);
252  Uy = eigenVector(T, eVals.y(), Ux);
253 
254  return Tensor2D<complex>(Ux, Uy);
255 }
256 
257 
259 {
260  const Vector2D<complex> eVals(eigenValues(T));
261 
262  return eigenVectors(T, eVals);
263 }
264 
265 
266 // ************************************************************************* //
Foam::VectorSpace::rootMax
static const Form rootMax
Definition: VectorSpace.H:119
Foam::roots::complex
Definition: Roots.H:57
Foam::VectorSpace::one
static const Form one
Definition: VectorSpace.H:116
Foam::VectorSpace::componentNames
static const char *const componentNames[]
Definition: VectorSpace.H:114
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::eigenVector
vector eigenVector(const symmTensor &T, const scalar eVal, const vector &standardBasis1, const vector &standardBasis2)
Definition: symmTensor.C:152
Foam::eigenValues
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
Definition: dimensionedTensor.C:149
Foam::eigenVectors
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
Definition: dimensionedTensor.C:160
tensor2D.H
A
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Foam::tensor2D
Tensor2D< scalar > tensor2D
Tensor2D of scalars, i.e. Tensor2D<scalar>.
Definition: symmTensor2D.H:68
Foam::Vector2D
Templated 2D Vector derived from VectorSpace adding construction from 2 components,...
Definition: Vector2D.H:55
Foam::sign
dimensionedScalar sign(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:166
Foam::Tensor2D::I
static const Tensor2D I
Definition: Tensor2D.H:79
Foam::Swap
void Swap(DynamicList< T, SizeMin1 > &a, DynamicList< T, SizeMin2 > &b)
Definition: DynamicListI.H:909
Foam::Tensor2D
A templated (2 x 2) tensor of objects of <T> derived from VectorSpace.
Definition: Tensor2D.H:59
Foam::VectorSpace::min
static const Form min
Definition: VectorSpace.H:118
Foam::constant::physicoChemical::b
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Foam::FatalError
error FatalError
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:137
T
const volScalarField & T
Definition: createFieldRefs.H:2
Foam::VectorSpace::rootMin
static const Form rootMin
Definition: VectorSpace.H:120
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:372
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::nl
constexpr char nl
Definition: Ostream.H:385
Foam::VectorSpace::max
static const Form max
Definition: VectorSpace.H:117
f
labelList f(nPoints)
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::VectorSpace::zero
static const Form zero
Definition: VectorSpace.H:115
quadraticEqn.H
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:298
Foam::VectorSpace::typeName
static const char *const typeName
Definition: VectorSpace.H:113