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