src/foam/primitives/Tensor2D/tensor2D/tensor2D.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "tensor2D.H"
  27 #include "mathematicalConstants.H"
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 template<>
  37 const char* const tensor2D::typeName = "tensor2D";
  38
  39 template<>
  40 const char* tensor2D::componentNames[] =
  41 {
  42     "xx", "xy",
  43     "yx", "yy"
  44 };
  45
  46 template<>
  47 const tensor2D tensor2D::zero
  48 (
  49     0, 0,
  50     0, 0
  51 );
  52
  53 template<>
  54 const tensor2D tensor2D::one
  55 (
  56     1, 1,
  57     1, 1
  58 );
  59
  60 template<>
  61 const tensor2D tensor2D::max
  62 (
  63     VGREAT, VGREAT,
  64     VGREAT, VGREAT
  65 );
  66
  67 template<>
  68 const tensor2D tensor2D::min
  69 (
  70     -VGREAT, -VGREAT,
  71     -VGREAT, -VGREAT
  72 );
  73
  74
  75 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  76
  77 // Return eigenvalues in ascending order of absolute values
  78 vector2D eigenValues(const tensor2D& t)
  79 {
  80     scalar i = 0;
  81     scalar ii = 0;
  82
  83     if (mag(t.xy()) < SMALL && mag(t.yx()) < SMALL)
  84     {
  85         i = t.xx();
  86         ii = t.yy();
  87     }
  88     else
  89     {
  90         scalar mb = t.xx() + t.yy();
  91         scalar c = t.xx()*t.yy() - t.xy()*t.yx();
  92
  93         // If there is a zero root
  94         if (mag(c) < SMALL)
  95         {
  96             i = 0;
  97             ii = mb;
  98         }
  99         else
 100         {
 101             scalar disc = sqr(mb) - 4*c;
 102
 103             if (disc > 0)
 104             {
 105                 scalar q = sqrt(disc);
 106
 107                 i = 0.5*(mb - q);
 108                 ii = 0.5*(mb + q);
 109             }
 110             else
 111             {
 112                 FatalErrorIn("eigenValues(const tensor2D&)")
 113                     << "zero and complex eigenvalues in tensor2D: " << t
 114                     << abort(FatalError);
 115             }
 116         }
 117     }
 118
 119     // Sort the eigenvalues into ascending order
 120     if (mag(i) > mag(ii))
 121     {
 122         Swap(i, ii);
 123     }
 124
 125     return vector2D(i, ii);
 126 }
 127
 128
 129 vector2D eigenVector(const tensor2D& t, const scalar lambda)
 130 {
 131     if (lambda < SMALL)
 132     {
 133         return vector2D::zero;
 134     }
 135
 136     if (mag(t.xy()) < SMALL && mag(t.yx()) < SMALL)
 137     {
 138         if (lambda > min(t.xx(), t.yy()))
 139         {
 140             return vector2D(1, 0);
 141         }
 142         else
 143         {
 144             return vector2D(0, 1);
 145         }
 146     }
 147     else if (mag(t.xy()) < SMALL)
 148     {
 149         return vector2D(lambda - t.yy(), t.yx());
 150     }
 151     else
 152     {
 153         return vector2D(t.xy(), lambda - t.yy());
 154     }
 155 }
 156
 157
 158 tensor2D eigenVectors(const tensor2D& t)
 159 {
 160     vector2D evals(eigenValues(t));
 161
 162     tensor2D evs
 163     (
 164         eigenVector(t, evals.x()),
 165         eigenVector(t, evals.y())
 166     );
 167
 168     return evs;
 169 }
 170
 171
 172 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 173
 174 } // End namespace Foam
 175
 176 // ************************************************************************* //