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