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) 1991-2009 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 template<>
  62 const tensor2D tensor2D::max
  63 (
  64     VGREAT, VGREAT,
  65     VGREAT, VGREAT
  66 );
  67
  68 template<>
  69 const tensor2D tensor2D::min
  70 (
  71     -VGREAT, -VGREAT,
  72     -VGREAT, -VGREAT
  73 );
  74
  75
  76 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  77
  78 // Return eigenvalues in ascending order of absolute values
  79 vector2D eigenValues(const tensor2D& t)
  80 {
  81     scalar i = 0;
  82     scalar ii = 0;
  83
  84     if (mag(t.xy()) < SMALL && mag(t.yx()) < SMALL)
  85     {
  86         i = t.xx();
  87         ii = t.yy();
  88     }
  89     else
  90     {
  91         scalar mb = t.xx() + t.yy();
  92         scalar c = t.xx()*t.yy() - t.xy()*t.yx();
  93
  94         // If there is a zero root
  95         if (mag(c) < SMALL)
  96         {
  97             i = 0;
  98             ii = mb;
  99         }
 100         else
 101         {
 102             scalar disc = sqr(mb) - 4*c;
 103
 104             if (disc > 0)
 105             {
 106                 scalar q = sqrt(disc);
 107
 108                 i = 0.5*(mb - q);
 109                 ii = 0.5*(mb + q);
 110             }
 111             else
 112             {
 113                 FatalErrorIn("eigenValues(const tensor2D&)")
 114                     << "zero and complex eigenvalues in tensor2D: " << t
 115                     << abort(FatalError);
 116             }
 117         }
 118     }
 119
 120     // Sort the eigenvalues into ascending order
 121     if (mag(i) > mag(ii))
 122     {
 123         Swap(i, ii);
 124     }
 125
 126     return vector2D(i, ii);
 127 }
 128
 129
 130 vector2D eigenVector(const tensor2D& t, const scalar lambda)
 131 {
 132     if (lambda < SMALL)
 133     {
 134         return vector2D::zero;
 135     }
 136
 137     if (mag(t.xy()) < SMALL && mag(t.yx()) < SMALL)
 138     {
 139         if (lambda > min(t.xx(), t.yy()))
 140         {
 141             return vector2D(1, 0);
 142         }
 143         else
 144         {
 145             return vector2D(0, 1);
 146         }
 147     }
 148     else if (mag(t.xy()) < SMALL)
 149     {
 150         return vector2D(lambda - t.yy(), t.yx());
 151     }
 152     else
 153     {
 154         return vector2D(t.xy(), lambda - t.yy());
 155     }
 156 }
 157
 158
 159 tensor2D eigenVectors(const tensor2D& t)
 160 {
 161     vector2D evals(eigenValues(t));
 162
 163     tensor2D evs;
 164     evs.x() = eigenVector(t, evals.x());
 165     evs.y() = eigenVector(t, evals.y());
 166
 167     return evs;
 168 }
 169
 170
 171 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 172
 173 } // End namespace Foam
 174
 175 // ************************************************************************* //