src/OpenFOAM/primitives/tensor/tensor.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 "tensor.H"
  28 #include "mathematicalConstants.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 template<>
  38 const char* const tensor::typeName = "tensor";
  39
  40 template<>
  41 const char* tensor::componentNames[] =
  42 {
  43     "xx", "xy", "xz",
  44     "yx", "yy", "yz",
  45     "zx", "zy", "zz"
  46 };
  47
  48 template<>
  49 const tensor tensor::zero
  50 (
  51     0, 0, 0,
  52     0, 0, 0,
  53     0, 0, 0
  54 );
  55
  56 template<>
  57 const tensor tensor::one
  58 (
  59     1, 1, 1,
  60     1, 1, 1,
  61     1, 1, 1
  62 );
  63
  64
  65 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  66
  67 // Return eigenvalues in ascending order of absolute values
  68 vector eigenValues(const tensor& t)
  69 {
  70     scalar i = 0;
  71     scalar ii = 0;
  72     scalar iii = 0;
  73
  74     if
  75     (
  76         (
  77             mag(t.xy()) + mag(t.xz()) + mag(t.yx())
  78           + mag(t.yz()) + mag(t.zx()) + mag(t.zy())
  79         )
  80       < SMALL
  81     )
  82     {
  83         // diagonal matrix
  84         i = t.xx();
  85         ii = t.yy();
  86         iii = t.zz();
  87     }
  88     else
  89     {
  90         scalar a = -t.xx() - t.yy() - t.zz();
  91
  92         scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
  93             - t.xy()*t.yx() - t.xz()*t.zx() - t.yz()*t.zy();
  94
  95         scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.zx()
  96             - t.xz()*t.yx()*t.zy() + t.xz()*t.yy()*t.zx()
  97             + t.xy()*t.yx()*t.zz() + t.xx()*t.yz()*t.zy();
  98
  99         // If there is a zero root
 100         if (mag(c) < 1.0e-100)
 101         {
 102             scalar disc = sqr(a) - 4*b;
 103
 104             if (disc >= -SMALL)
 105             {
 106                 scalar q = -0.5*sqrt(max(0.0, disc));
 107
 108                 i = 0;
 109                 ii = -0.5*a + q;
 110                 iii = -0.5*a - q;
 111             }
 112             else
 113             {
 114                 FatalErrorIn("eigenValues(const tensor&)")
 115                     << "zero and complex eigenvalues in tensor: " << t
 116                     << abort(FatalError);
 117             }
 118         }
 119         else
 120         {
 121             scalar Q = (a*a - 3*b)/9;
 122             scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
 123
 124             scalar R2 = sqr(R);
 125             scalar Q3 = pow3(Q);
 126
 127             // Three different real roots
 128             if (R2 < Q3)
 129             {
 130                 scalar sqrtQ = sqrt(Q);
 131                 scalar theta = acos(R/(Q*sqrtQ));
 132
 133                 scalar m2SqrtQ = -2*sqrtQ;
 134                 scalar aBy3 = a/3;
 135
 136                 i = m2SqrtQ*cos(theta/3) - aBy3;
 137                 ii = m2SqrtQ*cos((theta + mathematicalConstant::twoPi)/3)
 138                     - aBy3;
 139                 iii = m2SqrtQ*cos((theta - mathematicalConstant::twoPi)/3)
 140                     - aBy3;
 141             }
 142             else
 143             {
 144                 scalar A = cbrt(R + sqrt(R2 - Q3));
 145
 146                 // Three equal real roots
 147                 if (A < SMALL)
 148                 {
 149                     scalar root = -a/3;
 150                     return vector(root, root, root);
 151                 }
 152                 else
 153                 {
 154                     // Complex roots
 155                     WarningIn("eigenValues(const tensor&)")
 156                         << "complex eigenvalues detected for tensor: " << t
 157                         << endl;
 158
 159                     return vector::zero;
 160                 }
 161             }
 162         }
 163     }
 164
 165
 166     // Sort the eigenvalues into ascending order
 167     if (mag(i) > mag(ii))
 168     {
 169         Swap(i, ii);
 170     }
 171
 172     if (mag(ii) > mag(iii))
 173     {
 174         Swap(ii, iii);
 175     }
 176
 177     if (mag(i) > mag(ii))
 178     {
 179         Swap(i, ii);
 180     }
 181
 182     return vector(i, ii, iii);
 183 }
 184
 185
 186 vector eigenVector(const tensor& t, const scalar lambda)
 187 {
 188     if (mag(lambda) < SMALL)
 189     {
 190         return vector::zero;
 191     }
 192
 193     // Construct the matrix for the eigenvector problem
 194     tensor A(t - lambda*I);
 195
 196     // Calculate the sub-determinants of the 3 components
 197     scalar sd0 = A.yy()*A.zz() - A.yz()*A.zy();
 198     scalar sd1 = A.xx()*A.zz() - A.xz()*A.zx();
 199     scalar sd2 = A.xx()*A.yy() - A.xy()*A.yx();
 200
 201     scalar magSd0 = mag(sd0);
 202     scalar magSd1 = mag(sd1);
 203     scalar magSd2 = mag(sd2);
 204
 205     // Evaluate the eigenvector using the largest sub-determinant
 206     if (magSd0 > magSd1 && magSd0 > magSd2 && magSd0 > SMALL)
 207     {
 208         vector ev
 209         (
 210             1,
 211             (A.yz()*A.zx() - A.zz()*A.yx())/sd0,
 212             (A.zy()*A.yx() - A.yy()*A.zx())/sd0
 213         );
 214         ev /= mag(ev);
 215
 216         return ev;
 217     }
 218     else if (magSd1 > magSd2 && magSd1 > SMALL)
 219     {
 220         vector ev
 221         (
 222             (A.xz()*A.zy() - A.zz()*A.xy())/sd1,
 223             1,
 224             (A.zx()*A.xy() - A.xx()*A.zy())/sd1
 225         );
 226         ev /= mag(ev);
 227
 228         return ev;
 229     }
 230     else if (magSd2 > SMALL)
 231     {
 232         vector ev
 233         (
 234             (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
 235             (A.yx()*A.xz() - A.xx()*A.yz())/sd2,
 236             1
 237         );
 238         ev /= mag(ev);
 239
 240         return ev;
 241     }
 242     else
 243     {
 244         return vector::zero;
 245     }
 246 }
 247
 248
 249 tensor eigenVectors(const tensor& t)
 250 {
 251     vector evals(eigenValues(t));
 252
 253     tensor evs;
 254     evs.x() = eigenVector(t, evals.x());
 255     evs.y() = eigenVector(t, evals.y());
 256     evs.z() = eigenVector(t, evals.z());
 257
 258     return evs;
 259 }
 260
 261
 262 // Return eigenvalues in ascending order of absolute values
 263 vector eigenValues(const symmTensor& t)
 264 {
 265     scalar i = 0;
 266     scalar ii = 0;
 267     scalar iii = 0;
 268
 269     if
 270     (
 271         (
 272             mag(t.xy()) + mag(t.xz()) + mag(t.xy())
 273           + mag(t.yz()) + mag(t.xz()) + mag(t.yz())
 274         )
 275       < SMALL
 276     )
 277     {
 278         // diagonal matrix
 279         i = t.xx();
 280         ii = t.yy();
 281         iii = t.zz();
 282     }
 283     else
 284     {
 285         scalar a = -t.xx() - t.yy() - t.zz();
 286
 287         scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
 288             - t.xy()*t.xy() - t.xz()*t.xz() - t.yz()*t.yz();
 289
 290         scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.xz()
 291             - t.xz()*t.xy()*t.yz() + t.xz()*t.yy()*t.xz()
 292             + t.xy()*t.xy()*t.zz() + t.xx()*t.yz()*t.yz();
 293
 294         // If there is a zero root
 295         if (mag(c) < 1.0e-100)
 296         {
 297             scalar disc = sqr(a) - 4*b;
 298
 299             if (disc >= -SMALL)
 300             {
 301                 scalar q = -0.5*sqrt(max(0.0, disc));
 302
 303                 i = 0;
 304                 ii = -0.5*a + q;
 305                 iii = -0.5*a - q;
 306             }
 307             else
 308             {
 309                 FatalErrorIn("eigenValues(const tensor&)")
 310                     << "zero and complex eigenvalues in tensor: " << t
 311                     << abort(FatalError);
 312             }
 313         }
 314         else
 315         {
 316             scalar Q = (a*a - 3*b)/9;
 317             scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
 318
 319             scalar R2 = sqr(R);
 320             scalar Q3 = pow3(Q);
 321
 322             // Three different real roots
 323             if (R2 < Q3)
 324             {
 325                 scalar sqrtQ = sqrt(Q);
 326                 scalar theta = acos(R/(Q*sqrtQ));
 327
 328                 scalar m2SqrtQ = -2*sqrtQ;
 329                 scalar aBy3 = a/3;
 330
 331                 i = m2SqrtQ*cos(theta/3) - aBy3;
 332                 ii = m2SqrtQ*cos((theta + mathematicalConstant::twoPi)/3)
 333                     - aBy3;
 334                 iii = m2SqrtQ*cos((theta - mathematicalConstant::twoPi)/3)
 335                     - aBy3;
 336             }
 337             else
 338             {
 339                 scalar A = cbrt(R + sqrt(R2 - Q3));
 340
 341                 // Three equal real roots
 342                 if (A < SMALL)
 343                 {
 344                     scalar root = -a/3;
 345                     return vector(root, root, root);
 346                 }
 347                 else
 348                 {
 349                     // Complex roots
 350                     WarningIn("eigenValues(const symmTensor&)")
 351                         << "complex eigenvalues detected for symmTensor: " << t
 352                         << endl;
 353
 354                     return vector::zero;
 355                 }
 356             }
 357         }
 358     }
 359
 360
 361     // Sort the eigenvalues into ascending order
 362     if (mag(i) > mag(ii))
 363     {
 364         Swap(i, ii);
 365     }
 366
 367     if (mag(ii) > mag(iii))
 368     {
 369         Swap(ii, iii);
 370     }
 371
 372     if (mag(i) > mag(ii))
 373     {
 374         Swap(i, ii);
 375     }
 376
 377     return vector(i, ii, iii);
 378 }
 379
 380
 381 vector eigenVector(const symmTensor& t, const scalar lambda)
 382 {
 383     if (mag(lambda) < SMALL)
 384     {
 385         return vector::zero;
 386     }
 387
 388     // Construct the matrix for the eigenvector problem
 389     symmTensor A(t - lambda*I);
 390
 391     // Calculate the sub-determinants of the 3 components
 392     scalar sd0 = A.yy()*A.zz() - A.yz()*A.yz();
 393     scalar sd1 = A.xx()*A.zz() - A.xz()*A.xz();
 394     scalar sd2 = A.xx()*A.yy() - A.xy()*A.xy();
 395
 396     scalar magSd0 = mag(sd0);
 397     scalar magSd1 = mag(sd1);
 398     scalar magSd2 = mag(sd2);
 399
 400     // Evaluate the eigenvector using the largest sub-determinant
 401     if (magSd0 > magSd1 && magSd0 > magSd2 && magSd0 > SMALL)
 402     {
 403         vector ev
 404         (
 405             1,
 406             (A.yz()*A.xz() - A.zz()*A.xy())/sd0,
 407             (A.yz()*A.xy() - A.yy()*A.xz())/sd0
 408         );
 409         ev /= mag(ev);
 410
 411         return ev;
 412     }
 413     else if (magSd1 > magSd2 && magSd1 > SMALL)
 414     {
 415         vector ev
 416         (
 417             (A.xz()*A.yz() - A.zz()*A.xy())/sd1,
 418             1,
 419             (A.xz()*A.xy() - A.xx()*A.yz())/sd1
 420         );
 421         ev /= mag(ev);
 422
 423         return ev;
 424     }
 425     else if (magSd2 > SMALL)
 426     {
 427         vector ev
 428         (
 429             (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
 430             (A.xy()*A.xz() - A.xx()*A.yz())/sd2,
 431             1
 432         );
 433         ev /= mag(ev);
 434
 435         return ev;
 436     }
 437     else
 438     {
 439         return vector::zero;
 440     }
 441 }
 442
 443
 444 tensor eigenVectors(const symmTensor& t)
 445 {
 446     vector evals(eigenValues(t));
 447
 448     tensor evs;
 449     evs.x() = eigenVector(t, evals.x());
 450     evs.y() = eigenVector(t, evals.y());
 451     evs.z() = eigenVector(t, evals.z());
 452
 453     return evs;
 454 }
 455
 456
 457 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 458
 459 } // End namespace Foam
 460
 461 // ************************************************************************* //