tensor2D.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "tensor2D.H"
27 #include "quadraticEqn.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 template<>
32 const char* const Foam::tensor2D::vsType::typeName = "tensor2D";
33 
34 template<>
35 const char* const Foam::tensor2D::vsType::componentNames[] =
36 {
37  "xx", "xy",
38  "yx", "yy"
39 };
40 
41 template<>
43 (
44  tensor2D::uniform(0)
45 );
46 
47 template<>
49 (
50  tensor2D::uniform(1)
51 );
52 
53 template<>
55 (
56  tensor2D::uniform(vGreat)
57 );
58 
59 template<>
61 (
62  tensor2D::uniform(-vGreat)
63 );
64 
65 template<>
67 (
68  tensor2D::uniform(rootVGreat)
69 );
70 
71 template<>
73 (
74  tensor2D::uniform(-rootVGreat)
75 );
76 
77 template<>
79 (
80  1, 0,
81  0, 1
82 );
83 
84 
85 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
86 
88 {
89  // Coefficients of the characteristic quadratic polynomial (a = 1)
90  const scalar b = - t.xx() - t.yy();
91  const scalar c = t.xx()*t.yy() - t.xy()*t.yx();
92 
93  // Solve
94  Roots<2> roots = quadraticEqn(1, b, c).roots();
95 
96  // Check the root types
97  vector2D lambda = vector2D::zero;
98  forAll(roots, i)
99  {
100  switch (roots.type(i))
101  {
102  case rootType::real:
103  lambda[i] = roots[i];
104  break;
105  case rootType::complex:
107  << "Complex eigenvalues detected for tensor: " << t
108  << endl;
109  lambda[i] = 0;
110  break;
111  case rootType::posInf:
112  lambda[i] = vGreat;
113  break;
114  case rootType::negInf:
115  lambda[i] = - vGreat;
116  break;
117  case rootType::nan:
119  << "Eigenvalue calculation failed for tensor: " << t
120  << exit(FatalError);
121  }
122  }
123 
124  // Sort the eigenvalues into ascending order
125  if (lambda.x() > lambda.y())
126  {
127  Swap(lambda.x(), lambda.y());
128  }
129 
130  return lambda;
131 }
132 
133 
135 (
136  const tensor2D& T,
137  const scalar lambda,
138  const vector2D& direction1
139 )
140 {
141  // Construct the linear system for this eigenvalue
142  tensor2D A(T - lambda*tensor2D::I);
143 
144  // Evaluate the eigenvector using the largest divisor
145  if (mag(A.yy()) > mag(A.xx()) && mag(A.yy()) > small)
146  {
147  vector2D ev(1, - A.yx()/A.yy());
148 
149  return ev/mag(ev);
150  }
151  else if (mag(A.xx()) > small)
152  {
153  vector2D ev(- A.xy()/A.xx(), 1);
154 
155  return ev/mag(ev);
156  }
157 
158  // Repeated eigenvalue
159  return vector2D(- direction1.y(), direction1.x());
160 }
161 
162 
163 Foam::tensor2D Foam::eigenVectors(const tensor2D& T, const vector2D& lambdas)
164 {
165  vector2D Ux(1, 0), Uy(0, 1);
166 
167  Ux = eigenVector(T, lambdas.x(), Uy);
168  Uy = eigenVector(T, lambdas.y(), Ux);
169 
170  return tensor2D(Ux, Uy);
171 }
172 
173 
175 {
176  const vector2D lambdas(eigenValues(T));
177 
178  return eigenVectors(T, lambdas);
179 }
180 
181 
182 // ************************************************************************* //
vector eigenVector(const tensor &T, const scalar lambda, const vector &direction1, const vector &direction2)
Definition: tensor.C:137
static const char *const typeName
Definition: VectorSpace.H:111
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
static const Form max
Definition: VectorSpace.H:115
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
Vector2D< scalar > vector2D
vector2D obtained from generic Vector2D
Definition: vector2D.H:49
static const char *const componentNames[]
Definition: VectorSpace.H:112
dimensionedVector eigenValues(const dimensionedTensor &dt)
static const Form rootMin
Definition: VectorSpace.H:118
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
dimensionedScalar lambda(viscosity->lookup("lambda"))
Templated 2D tensor derived from VectorSpace adding construction from 4 components, element access using xx(), xy(), yx() and yy() member functions and the iner-product (dot-product) and outer-product of two Vector2Ds (tensor-product) operators.
Definition: Tensor2D.H:56
static const Form min
Definition: VectorSpace.H:116
void Swap(T &a, T &b)
Definition: Swap.H:43
static const Identity< scalar > I
Definition: Identity.H:93
static const Tensor2D I
Definition: Tensor2D.H:75
dimensionedTensor eigenVectors(const dimensionedTensor &dt)
static const Form rootMax
Definition: VectorSpace.H:117
#define WarningInFunction
Report a warning using Foam::Warning.
static const Form one
Definition: VectorSpace.H:114
dimensioned< scalar > mag(const dimensioned< Type > &)
Templated 2D Vector derived from VectorSpace adding construction from 2 components, element access using x() and y() member functions and the inner-product (dot-product).
Definition: Vector2D.H:51
Tensor2D< scalar > tensor2D
Tensor2D or scalars.
Definition: tensor2D.H:49
static const Form zero
Definition: VectorSpace.H:113