src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
   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
  28 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  29
  30 namespace Foam
  31 {
  32
  33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  34
  35 template<>
  36 const char* const tensor2D::typeName = "tensor2D";
  37
  38 template<>
  39 const char* tensor2D::componentNames[] =
  40 {
  41     "xx", "xy",
  42     "yx", "yy"
  43 };
  44
  45 template<>
  46 const tensor2D tensor2D::zero
  47 (
  48     0, 0,
  49     0, 0
  50 );
  51
  52 template<>
  53 const tensor2D tensor2D::one
  54 (
  55     1, 1,
  56     1, 1
  57 );
  58
  59 template<>
  60 const tensor2D tensor2D::max
  61 (
  62     VGREAT, VGREAT,
  63     VGREAT, VGREAT
  64 );
  65
  66 template<>
  67 const tensor2D tensor2D::min
  68 (
  69     -VGREAT, -VGREAT,
  70     -VGREAT, -VGREAT
  71 );
  72
  73
  74 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  75
  76 // Return eigenvalues in ascending order of absolute values
  77 vector2D eigenValues(const tensor2D& t)
  78 {
  79     scalar i = 0;
  80     scalar ii = 0;
  81
  82     if (mag(t.xy()) < SMALL && mag(t.yx()) < SMALL)
  83     {
  84         i = t.xx();
  85         ii = t.yy();
  86     }
  87     else
  88     {
  89         scalar mb = t.xx() + t.yy();
  90         scalar c = t.xx()*t.yy() - t.xy()*t.yx();
  91
  92         // If there is a zero root
  93         if (mag(c) < SMALL)
  94         {
  95             i = 0;
  96             ii = mb;
  97         }
  98         else
  99         {
 100             scalar disc = sqr(mb) - 4*c;
 101
 102             if (disc > 0)
 103             {
 104                 scalar q = sqrt(disc);
 105
 106                 i = 0.5*(mb - q);
 107                 ii = 0.5*(mb + q);
 108             }
 109             else
 110             {
 111                 FatalErrorIn("eigenValues(const tensor2D&)")
 112                     << "zero and complex eigenvalues in tensor2D: " << t
 113                     << abort(FatalError);
 114             }
 115         }
 116     }
 117
 118     // Sort the eigenvalues into ascending order
 119     if (mag(i) > mag(ii))
 120     {
 121         Swap(i, ii);
 122     }
 123
 124     return vector2D(i, ii);
 125 }
 126
 127
 128 vector2D eigenVector(const tensor2D& t, const scalar lambda)
 129 {
 130     if (lambda < SMALL)
 131     {
 132         return vector2D::zero;
 133     }
 134
 135     if (mag(t.xy()) < SMALL && mag(t.yx()) < SMALL)
 136     {
 137         if (lambda > min(t.xx(), t.yy()))
 138         {
 139             return vector2D(1, 0);
 140         }
 141         else
 142         {
 143             return vector2D(0, 1);
 144         }
 145     }
 146     else if (mag(t.xy()) < SMALL)
 147     {
 148         return vector2D(lambda - t.yy(), t.yx());
 149     }
 150     else
 151     {
 152         return vector2D(t.xy(), lambda - t.yy());
 153     }
 154 }
 155
 156
 157 tensor2D eigenVectors(const tensor2D& t)
 158 {
 159     vector2D evals(eigenValues(t));
 160
 161     tensor2D evs
 162     (
 163         eigenVector(t, evals.x()),
 164         eigenVector(t, evals.y())
 165     );
 166
 167     return evs;
 168 }
 169
 170
 171 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 172
 173 } // End namespace Foam
 174
 175 // ************************************************************************* //