src/foam/primitives/Tensor/tensor/tensor.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  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(scalar(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 = m2SqrtQ*cos((theta + mathematicalConstant::twoPi)/3)
 153                     - aBy3;
 154                 iii = m2SqrtQ*cos((theta - mathematicalConstant::twoPi)/3)
 155                     - aBy3;
 156             }
 157             else
 158             {
 159                 scalar A = cbrt(R + sqrt(R2 - Q3));
 160
 161                 // Three equal real roots
 162                 if (A < SMALL)
 163                 {
 164                     scalar root = -a/3;
 165                     return vector(root, root, root);
 166                 }
 167                 else
 168                 {
 169                     // Complex roots
 170                     WarningIn("eigenValues(const tensor&)")
 171                         << "complex eigenvalues detected for tensor: " << t
 172                         << endl;
 173
 174                     return vector::zero;
 175                 }
 176             }
 177         }
 178     }
 179
 180
 181     // Sort the eigenvalues into ascending order
 182     if (mag(i) > mag(ii))
 183     {
 184         Swap(i, ii);
 185     }
 186
 187     if (mag(ii) > mag(iii))
 188     {
 189         Swap(ii, iii);
 190     }
 191
 192     if (mag(i) > mag(ii))
 193     {
 194         Swap(i, ii);
 195     }
 196
 197     return vector(i, ii, iii);
 198 }
 199
 200
 201 vector eigenVector(const tensor& t, const scalar lambda)
 202 {
 203     // Construct the matrix for the eigenvector problem
 204     tensor A(t - lambda*I);
 205
 206     // Calculate the sub-determinants of the 3 components
 207     scalar sd0 = A.yy()*A.zz() - A.yz()*A.zy();
 208     scalar sd1 = A.xx()*A.zz() - A.xz()*A.zx();
 209     scalar sd2 = A.xx()*A.yy() - A.xy()*A.yx();
 210
 211     scalar magSd0 = mag(sd0);
 212     scalar magSd1 = mag(sd1);
 213     scalar magSd2 = mag(sd2);
 214
 215     // Evaluate the eigenvector using the largest sub-determinant
 216     if (magSd0 > magSd1 && magSd0 > magSd2 && magSd0 > SMALL)
 217     {
 218         vector ev
 219         (
 220             1,
 221             (A.yz()*A.zx() - A.zz()*A.yx())/sd0,
 222             (A.zy()*A.yx() - A.yy()*A.zx())/sd0
 223         );
 224         ev /= mag(ev);
 225
 226         return ev;
 227     }
 228     else if (magSd1 > magSd2 && magSd1 > SMALL)
 229     {
 230         vector ev
 231         (
 232             (A.xz()*A.zy() - A.zz()*A.xy())/sd1,
 233             1,
 234             (A.zx()*A.xy() - A.xx()*A.zy())/sd1
 235         );
 236         ev /= mag(ev);
 237
 238         return ev;
 239     }
 240     else if (magSd2 > SMALL)
 241     {
 242         vector ev
 243         (
 244             (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
 245             (A.yx()*A.xz() - A.xx()*A.yz())/sd2,
 246             1
 247         );
 248         ev /= mag(ev);
 249
 250         return ev;
 251     }
 252     else
 253     {
 254         // Double identical eigen-value
 255         if (mag(A.xx()) > SMALL)
 256         {
 257             return vector(0, 1, 0);
 258         }
 259         else if (mag(A.yy()) > SMALL)
 260         {
 261             return vector(0, 0, 1);
 262         }
 263         else if (mag(A.zz()) > SMALL)
 264         {
 265             return vector(1, 0, 0);
 266         }
 267         else
 268         {
 269             // Everything is zero.  Return arbitrary vector
 270             return vector(1, 0, 0);
 271         }
 272     }
 273 }
 274
 275
 276 tensor eigenVectors(const tensor& t)
 277 {
 278     vector evals(eigenValues(t));
 279
 280     // Modification for strict-aliasing compliance.
 281     // Sandeep Menon, 01/Dec/2010
 282
 283     // Test for null eigen values to return a not null eigen vector.
 284     // Jovani Favero, 18/Nov/2009. Adjusted by HJ, to correct for multiple zero
 285     // entries in eigenvalues
 286
 287     // Start with largest eigen-value: if this is zero, all are zero
 288     // and xyz tensor is returned
 289     if (mag(evals.z()) < SMALL)
 290     {
 291         return I;
 292     }
 293
 294     // One valid eigen-vector.  Manufacture second and third direction
 295     // as orthogonal vectors onto the first one with arbitrary orientation
 296     if (mag(evals.y()) < SMALL)
 297     {
 298         // Take the z eigenvector and find a non-zero component
 299         vector zVec = eigenVector(t, evals.z());
 300
 301         vector yVec;
 302
 303         if (mag(zVec.z()) > 0)
 304         {
 305             // Rotate z into y
 306             yVec = vector(zVec.x(), -zVec.z(), zVec.y());
 307         }
 308         else if (mag(zVec.y()) > 0)
 309         {
 310             // Rotate y into x
 311             yVec = vector(-zVec.y(), zVec.x(), zVec.z());
 312         }
 313         else
 314         {
 315             // Rotate x into z
 316             yVec = vector(zVec.z(), zVec.y(), -zVec.x());
 317         }
 318
 319         vector xVec = yVec ^ zVec;
 320
 321         return tensor(xVec, yVec, zVec);
 322     }
 323
 324     if (mag(evals.x()) < SMALL)
 325     {
 326         vector xVec = eigenVector(t, evals.x());
 327         vector yVec = eigenVector(t, evals.y());
 328         vector zVec = eigenVector(t, evals.z());
 329
 330         // If second and third eigen value is the same, the two vectors
 331         // will be identical.  Fix this
 332         if (mag(evals.z() - evals.y()) < SMALL)
 333         {
 334             zVec = xVec ^ yVec;
 335         }
 336
 337         return tensor(xVec, yVec, zVec);
 338     }
 339
 340     return tensor
 341     (
 342         eigenVector(t, evals.x()),
 343         eigenVector(t, evals.y()),
 344         eigenVector(t, evals.z())
 345     );
 346 }
 347
 348
 349 // Return eigenvalues in ascending order of absolute values
 350 vector eigenValues(const symmTensor& t)
 351 {
 352     scalar i = 0;
 353     scalar ii = 0;
 354     scalar iii = 0;
 355
 356     if
 357     (
 358         (
 359             mag(t.xy()) + mag(t.xz()) + mag(t.xy())
 360           + mag(t.yz()) + mag(t.xz()) + mag(t.yz())
 361         )
 362       < SMALL
 363     )
 364     {
 365         // diagonal matrix
 366         i = t.xx();
 367         ii = t.yy();
 368         iii = t.zz();
 369     }
 370     else
 371     {
 372         scalar a = -t.xx() - t.yy() - t.zz();
 373
 374         scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
 375             - t.xy()*t.xy() - t.xz()*t.xz() - t.yz()*t.yz();
 376
 377         scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.xz()
 378             - t.xz()*t.xy()*t.yz() + t.xz()*t.yy()*t.xz()
 379             + t.xy()*t.xy()*t.zz() + t.xx()*t.yz()*t.yz();
 380
 381         // If there is a zero root
 382         if (mag(c) < SMALL)
 383         {
 384             const scalar disc = Foam::max(sqr(a) - 4*b, scalar(0));
 385
 386             scalar q = -0.5*sqrt(max(scalar(0), disc));
 387
 388             i = 0;
 389             ii = -0.5*a + q;
 390             iii = -0.5*a - q;
 391         }
 392         else
 393         {
 394             scalar Q = (a*a - 3*b)/9;
 395             scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
 396
 397             scalar R2 = sqr(R);
 398             scalar Q3 = pow3(Q);
 399
 400             // Three different real roots
 401             if (R2 < Q3)
 402             {
 403                 scalar sqrtQ = sqrt(Q);
 404                 scalar theta = acos(R/(Q*sqrtQ));
 405
 406                 scalar m2SqrtQ = -2*sqrtQ;
 407                 scalar aBy3 = a/3;
 408
 409                 i = m2SqrtQ*cos(theta/3) - aBy3;
 410                 ii = m2SqrtQ*cos((theta + mathematicalConstant::twoPi)/3)
 411                     - aBy3;
 412                 iii = m2SqrtQ*cos((theta - mathematicalConstant::twoPi)/3)
 413                     - aBy3;
 414             }
 415             else
 416             {
 417                 scalar A = cbrt(R + sqrt(R2 - Q3));
 418
 419                 // Three equal real roots
 420                 if (A < SMALL)
 421                 {
 422                     scalar root = -a/3;
 423                     return vector(root, root, root);
 424                 }
 425                 else
 426                 {
 427                     // Complex roots
 428                     WarningIn("eigenValues(const symmTensor&)")
 429                         << "complex eigenvalues detected for symmTensor: " << t
 430                         << endl;
 431
 432                     return vector::zero;
 433                 }
 434             }
 435         }
 436     }
 437
 438
 439     // Sort the eigenvalues into ascending order
 440     if (mag(i) > mag(ii))
 441     {
 442         Swap(i, ii);
 443     }
 444
 445     if (mag(ii) > mag(iii))
 446     {
 447         Swap(ii, iii);
 448     }
 449
 450     if (mag(i) > mag(ii))
 451     {
 452         Swap(i, ii);
 453     }
 454
 455     return vector(i, ii, iii);
 456 }
 457
 458
 459 vector eigenVector(const symmTensor& t, const scalar lambda)
 460 {
 461     // Construct the matrix for the eigenvector problem
 462     symmTensor A(t - lambda*I);
 463
 464     // Calculate the sub-determinants of the 3 components
 465     scalar sd0 = A.yy()*A.zz() - A.yz()*A.yz();
 466     scalar sd1 = A.xx()*A.zz() - A.xz()*A.xz();
 467     scalar sd2 = A.xx()*A.yy() - A.xy()*A.xy();
 468
 469     scalar magSd0 = mag(sd0);
 470     scalar magSd1 = mag(sd1);
 471     scalar magSd2 = mag(sd2);
 472
 473     // Evaluate the eigenvector using the largest sub-determinant
 474     if (magSd0 > magSd1 && magSd0 > magSd2 && magSd0 > SMALL)
 475     {
 476         vector ev
 477         (
 478             1,
 479             (A.yz()*A.xz() - A.zz()*A.xy())/sd0,
 480             (A.yz()*A.xy() - A.yy()*A.xz())/sd0
 481         );
 482         ev /= mag(ev);
 483
 484         return ev;
 485     }
 486     else if (magSd1 > magSd2 && magSd1 > SMALL)
 487     {
 488         vector ev
 489         (
 490             (A.xz()*A.yz() - A.zz()*A.xy())/sd1,
 491             1,
 492             (A.xz()*A.xy() - A.xx()*A.yz())/sd1
 493         );
 494         ev /= mag(ev);
 495
 496         return ev;
 497     }
 498     else if (magSd2 > SMALL)
 499     {
 500         vector ev
 501         (
 502             (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
 503             (A.xy()*A.xz() - A.xx()*A.yz())/sd2,
 504             1
 505         );
 506         ev /= mag(ev);
 507
 508         return ev;
 509     }
 510     else
 511     {
 512         // Double identical eigen-value
 513         if (mag(A.xx()) > SMALL)
 514         {
 515             return vector(0, 1, 0);
 516         }
 517         else if (mag(A.yy()) > SMALL)
 518         {
 519             return vector(0, 0, 1);
 520         }
 521         else if (mag(A.zz()) > SMALL)
 522         {
 523             return vector(1, 0, 0);
 524         }
 525         else
 526         {
 527             // Everything is zero.  Return arbitrary vector
 528             return vector(1, 0, 0);
 529         }
 530     }
 531 }
 532
 533
 534 tensor eigenVectors(const symmTensor& t)
 535 {
 536     vector evals(eigenValues(t));
 537
 538     // Modification for strict-aliasing compliance.
 539     // Sandeep Menon, 01/Dec/2010
 540
 541     // Test for null eigen values to return a not null eigen vector
 542     // Jovani Favero, 18/Nov/2009
 543
 544     // Start with largest eigen-value: if this is zero, all are zero
 545     // and original tensor is returned
 546     if (mag(evals.z()) < SMALL)
 547     {
 548         return I;
 549     }
 550
 551     // One valid eigen-vector.  Manufacture second and third direction
 552     // as orthogonal vectors onto the first one with arbitrary orientation
 553     if (mag(evals.y()) < SMALL)
 554     {
 555         // Take the z eigenvector and find a non-zero component
 556         vector zVec = eigenVector(t, evals.z());
 557
 558         vector yVec;
 559
 560         if (mag(zVec.z()) > 0)
 561         {
 562             // Rotate z into y
 563             yVec = vector(zVec.x(), -zVec.z(), zVec.y());
 564         }
 565         else if (mag(zVec.y()) > 0)
 566         {
 567             // Rotate y into x
 568             yVec = vector(-zVec.y(), zVec.x(), zVec.z());
 569         }
 570         else
 571         {
 572             // Rotate x into z
 573             yVec = vector(zVec.z(), zVec.y(), -zVec.x());
 574         }
 575
 576         vector xVec = yVec ^ zVec;
 577
 578         // Note different return
 579         return tensor(xVec, yVec, zVec);
 580     }
 581
 582     if (mag(evals.x()) < SMALL)
 583     {
 584         vector xVec = eigenVector(t, evals.x());
 585         vector yVec = eigenVector(t, evals.y());
 586         vector zVec = eigenVector(t, evals.z());
 587
 588         // If second and third eigen value is the same, the two vectors
 589         // will be identical.  Fix this
 590         if (mag(evals.z() - evals.y()) < SMALL)
 591         {
 592             zVec = xVec ^ yVec;
 593         }
 594
 595         return tensor(xVec, yVec, zVec);
 596     }
 597
 598     return tensor
 599     (
 600         eigenVector(t, evals.x()),
 601         eigenVector(t, evals.y()),
 602         eigenVector(t, evals.z())
 603     );
 604 }
 605
 606
 607 // Matrix inversion with singular value decomposition
 608 tensor hinv(const tensor& t)
 609 {
 610     static const scalar large = 1e10;
 611     static const scalar small = 1e-10;
 612
 613     if (det(t) > small)
 614     {
 615         return inv(t);
 616     }
 617     else
 618     {
 619         vector eig = eigenValues(t);
 620         tensor eigVecs = eigenVectors(t);
 621
 622         tensor zeroInv = tensor::zero;
 623
 624         // Test if all eigen values are zero.
 625         // If this happens then eig.z() = SMALL, and hinv(t)
 626         // returns a zero tensor.
 627         // Jovani Favero, 18/Nov/2009
 628         // Further bug fix: replace > with == and add SMALL to zeroInv
 629         // Dominik Christ, 7/Aug/2012
 630         if (mag(eig.z()) == large*mag(eig.z()))
 631         {
 632             // Three zero eigen values (z is largest in magnitude).
 633             // Return zero inverse
 634             return zeroInv;
 635         }
 636
 637         // Compare smaller eigen values and if out of range of large
 638         // consider them singular
 639
 640         if (mag(eig.z()) > large*mag(eig.x()))
 641         {
 642             // Make a tensor out of symmTensor sqr.  HJ, 24/Oct/2009
 643             zeroInv += tensor(sqr(eigVecs.x()));
 644         }
 645
 646         if (mag(eig.z()) > large*mag(eig.y()))
 647         {
 648             // Make a tensor out of symmTensor sqr.  HJ, 24/Oct/2009
 649             zeroInv += tensor(sqr(eigVecs.y()));
 650         }
 651
 652         return inv(t + zeroInv) - zeroInv;
 653     }
 654 }
 655
 656
 657 symmTensor hinv(const symmTensor& t)
 658 {
 659     static const scalar large = 1e10;
 660     static const scalar small = 1e-10;
 661
 662     if (det(t) > small)
 663     {
 664         return inv(t);
 665     }
 666     else
 667     {
 668         vector eig = eigenValues(t);
 669         tensor eigVecs = eigenVectors(t);
 670
 671         symmTensor zeroInv = symmTensor::zero;
 672
 673         // Test if all eigen values are zero,
 674         // If this happens then eig.z() = SMALL
 675         // and hinv(t) return a zero tensor.
 676         // Jovani Favero, 18/Nov/2009
 677         // Further bug fix: replace > with == and add SMALL to zeroInv
 678         // Dominik Christ, 7/Aug/2012
 679         if (mag(eig.z()) == large*mag(eig.z()))
 680         {
 681             // Three zero eigen values (z is largest in magnitude).
 682             // Return zero inverse
 683             return zeroInv;
 684         }
 685
 686         // Compare smaller eigen values and if out of range of large
 687         // consider them singular
 688
 689         if (mag(eig.z()) > large*mag(eig.x()))
 690         {
 691             zeroInv += sqr(eigVecs.x());
 692         }
 693
 694         if (mag(eig.z()) > large*mag(eig.y()))
 695         {
 696             zeroInv += sqr(eigVecs.y());
 697         }
 698
 699         return inv(t + zeroInv) - zeroInv;
 700     }
 701 }
 702
 703
 704 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 705
 706 } // End namespace Foam
 707
 708 // ************************************************************************* //