src/OpenFOAM/primitives/tensor2D/tensor2D.C

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