FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

tensor.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "tensor.H"
00027 #include <OpenFOAM/mathematicalConstants.H>
00028 
00029 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00030 
00031 namespace Foam
00032 {
00033 
00034 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00035 
00036 template<>
00037 const char* const tensor::typeName = "tensor";
00038 
00039 template<>
00040 const char* tensor::componentNames[] =
00041 {
00042     "xx", "xy", "xz",
00043     "yx", "yy", "yz",
00044     "zx", "zy", "zz"
00045 };
00046 
00047 template<>
00048 const tensor tensor::zero
00049 (
00050     0, 0, 0,
00051     0, 0, 0,
00052     0, 0, 0
00053 );
00054 
00055 template<>
00056 const tensor tensor::one
00057 (
00058     1, 1, 1,
00059     1, 1, 1,
00060     1, 1, 1
00061 );
00062 
00063 template<>
00064 const tensor tensor::max
00065 (
00066     VGREAT, VGREAT, VGREAT,
00067     VGREAT, VGREAT, VGREAT,
00068     VGREAT, VGREAT, VGREAT
00069 );
00070 
00071 template<>
00072 const tensor tensor::min
00073 (
00074     -VGREAT, -VGREAT, -VGREAT,
00075     -VGREAT, -VGREAT, -VGREAT,
00076     -VGREAT, -VGREAT, -VGREAT
00077 );
00078 
00079 
00080 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00081 
00082 // Return eigenvalues in ascending order of absolute values
00083 vector eigenValues(const tensor& t)
00084 {
00085     scalar i = 0;
00086     scalar ii = 0;
00087     scalar iii = 0;
00088 
00089     if
00090     (
00091         (
00092             mag(t.xy()) + mag(t.xz()) + mag(t.yx())
00093           + mag(t.yz()) + mag(t.zx()) + mag(t.zy())
00094         )
00095       < SMALL
00096     )
00097     {
00098         // diagonal matrix
00099         i = t.xx();
00100         ii = t.yy();
00101         iii = t.zz();
00102     }
00103     else
00104     {
00105         scalar a = -t.xx() - t.yy() - t.zz();
00106 
00107         scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
00108             - t.xy()*t.yx() - t.xz()*t.zx() - t.yz()*t.zy();
00109 
00110         scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.zx()
00111             - t.xz()*t.yx()*t.zy() + t.xz()*t.yy()*t.zx()
00112             + t.xy()*t.yx()*t.zz() + t.xx()*t.yz()*t.zy();
00113 
00114         // If there is a zero root
00115         if (mag(c) < 1.0e-100)
00116         {
00117             scalar disc = sqr(a) - 4*b;
00118 
00119             if (disc >= -SMALL)
00120             {
00121                 scalar q = -0.5*sqrt(max(0.0, disc));
00122 
00123                 i = 0;
00124                 ii = -0.5*a + q;
00125                 iii = -0.5*a - q;
00126             }
00127             else
00128             {
00129                 FatalErrorIn("eigenValues(const tensor&)")
00130                     << "zero and complex eigenvalues in tensor: " << t
00131                     << abort(FatalError);
00132             }
00133         }
00134         else
00135         {
00136             scalar Q = (a*a - 3*b)/9;
00137             scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
00138 
00139             scalar R2 = sqr(R);
00140             scalar Q3 = pow3(Q);
00141 
00142             // Three different real roots
00143             if (R2 < Q3)
00144             {
00145                 scalar sqrtQ = sqrt(Q);
00146                 scalar theta = acos(R/(Q*sqrtQ));
00147 
00148                 scalar m2SqrtQ = -2*sqrtQ;
00149                 scalar aBy3 = a/3;
00150 
00151                 i = m2SqrtQ*cos(theta/3) - aBy3;
00152                 ii = m2SqrtQ*cos((theta + mathematicalConstant::twoPi)/3)
00153                     - aBy3;
00154                 iii = m2SqrtQ*cos((theta - mathematicalConstant::twoPi)/3)
00155                     - aBy3;
00156             }
00157             else
00158             {
00159                 scalar A = cbrt(R + sqrt(R2 - Q3));
00160 
00161                 // Three equal real roots
00162                 if (A < SMALL)
00163                 {
00164                     scalar root = -a/3;
00165                     return vector(root, root, root);
00166                 }
00167                 else
00168                 {
00169                     // Complex roots
00170                     WarningIn("eigenValues(const tensor&)")
00171                         << "complex eigenvalues detected for tensor: " << t
00172                         << endl;
00173 
00174                     return vector::zero;
00175                 }
00176             }
00177         }
00178     }
00179 
00180 
00181     // Sort the eigenvalues into ascending order
00182     if (mag(i) > mag(ii))
00183     {
00184         Swap(i, ii);
00185     }
00186 
00187     if (mag(ii) > mag(iii))
00188     {
00189         Swap(ii, iii);
00190     }
00191 
00192     if (mag(i) > mag(ii))
00193     {
00194         Swap(i, ii);
00195     }
00196 
00197     return vector(i, ii, iii);
00198 }
00199 
00200 
00201 vector eigenVector(const tensor& t, const scalar lambda)
00202 {
00203     if (mag(lambda) < SMALL)
00204     {
00205         return vector::zero;
00206     }
00207 
00208     // Construct the matrix for the eigenvector problem
00209     tensor A(t - lambda*I);
00210 
00211     // Calculate the sub-determinants of the 3 components
00212     scalar sd0 = A.yy()*A.zz() - A.yz()*A.zy();
00213     scalar sd1 = A.xx()*A.zz() - A.xz()*A.zx();
00214     scalar sd2 = A.xx()*A.yy() - A.xy()*A.yx();
00215 
00216     scalar magSd0 = mag(sd0);
00217     scalar magSd1 = mag(sd1);
00218     scalar magSd2 = mag(sd2);
00219 
00220     // Evaluate the eigenvector using the largest sub-determinant
00221     if (magSd0 > magSd1 && magSd0 > magSd2 && magSd0 > SMALL)
00222     {
00223         vector ev
00224         (
00225             1,
00226             (A.yz()*A.zx() - A.zz()*A.yx())/sd0,
00227             (A.zy()*A.yx() - A.yy()*A.zx())/sd0
00228         );
00229         ev /= mag(ev);
00230 
00231         return ev;
00232     }
00233     else if (magSd1 > magSd2 && magSd1 > SMALL)
00234     {
00235         vector ev
00236         (
00237             (A.xz()*A.zy() - A.zz()*A.xy())/sd1,
00238             1,
00239             (A.zx()*A.xy() - A.xx()*A.zy())/sd1
00240         );
00241         ev /= mag(ev);
00242 
00243         return ev;
00244     }
00245     else if (magSd2 > SMALL)
00246     {
00247         vector ev
00248         (
00249             (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
00250             (A.yx()*A.xz() - A.xx()*A.yz())/sd2,
00251             1
00252         );
00253         ev /= mag(ev);
00254 
00255         return ev;
00256     }
00257     else
00258     {
00259         return vector::zero;
00260     }
00261 }
00262 
00263 
00264 tensor eigenVectors(const tensor& t)
00265 {
00266     vector evals(eigenValues(t));
00267 
00268     tensor evs
00269     (
00270         eigenVector(t, evals.x()),
00271         eigenVector(t, evals.y()),
00272         eigenVector(t, evals.z())
00273     );
00274 
00275     return evs;
00276 }
00277 
00278 
00279 // Return eigenvalues in ascending order of absolute values
00280 vector eigenValues(const symmTensor& t)
00281 {
00282     scalar i = 0;
00283     scalar ii = 0;
00284     scalar iii = 0;
00285 
00286     if
00287     (
00288         (
00289             mag(t.xy()) + mag(t.xz()) + mag(t.xy())
00290           + mag(t.yz()) + mag(t.xz()) + mag(t.yz())
00291         )
00292       < SMALL
00293     )
00294     {
00295         // diagonal matrix
00296         i = t.xx();
00297         ii = t.yy();
00298         iii = t.zz();
00299     }
00300     else
00301     {
00302         scalar a = -t.xx() - t.yy() - t.zz();
00303 
00304         scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
00305             - t.xy()*t.xy() - t.xz()*t.xz() - t.yz()*t.yz();
00306 
00307         scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.xz()
00308             - t.xz()*t.xy()*t.yz() + t.xz()*t.yy()*t.xz()
00309             + t.xy()*t.xy()*t.zz() + t.xx()*t.yz()*t.yz();
00310 
00311         // If there is a zero root
00312         if (mag(c) < 1.0e-100)
00313         {
00314             scalar disc = sqr(a) - 4*b;
00315 
00316             if (disc >= -SMALL)
00317             {
00318                 scalar q = -0.5*sqrt(max(0.0, disc));
00319 
00320                 i = 0;
00321                 ii = -0.5*a + q;
00322                 iii = -0.5*a - q;
00323             }
00324             else
00325             {
00326                 FatalErrorIn("eigenValues(const tensor&)")
00327                     << "zero and complex eigenvalues in tensor: " << t
00328                     << abort(FatalError);
00329             }
00330         }
00331         else
00332         {
00333             scalar Q = (a*a - 3*b)/9;
00334             scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
00335 
00336             scalar R2 = sqr(R);
00337             scalar Q3 = pow3(Q);
00338 
00339             // Three different real roots
00340             if (R2 < Q3)
00341             {
00342                 scalar sqrtQ = sqrt(Q);
00343                 scalar theta = acos(R/(Q*sqrtQ));
00344 
00345                 scalar m2SqrtQ = -2*sqrtQ;
00346                 scalar aBy3 = a/3;
00347 
00348                 i = m2SqrtQ*cos(theta/3) - aBy3;
00349                 ii = m2SqrtQ*cos((theta + mathematicalConstant::twoPi)/3)
00350                     - aBy3;
00351                 iii = m2SqrtQ*cos((theta - mathematicalConstant::twoPi)/3)
00352                     - aBy3;
00353             }
00354             else
00355             {
00356                 scalar A = cbrt(R + sqrt(R2 - Q3));
00357 
00358                 // Three equal real roots
00359                 if (A < SMALL)
00360                 {
00361                     scalar root = -a/3;
00362                     return vector(root, root, root);
00363                 }
00364                 else
00365                 {
00366                     // Complex roots
00367                     WarningIn("eigenValues(const symmTensor&)")
00368                         << "complex eigenvalues detected for symmTensor: " << t
00369                         << endl;
00370 
00371                     return vector::zero;
00372                 }
00373             }
00374         }
00375     }
00376 
00377 
00378     // Sort the eigenvalues into ascending order
00379     if (mag(i) > mag(ii))
00380     {
00381         Swap(i, ii);
00382     }
00383 
00384     if (mag(ii) > mag(iii))
00385     {
00386         Swap(ii, iii);
00387     }
00388 
00389     if (mag(i) > mag(ii))
00390     {
00391         Swap(i, ii);
00392     }
00393 
00394     return vector(i, ii, iii);
00395 }
00396 
00397 
00398 vector eigenVector(const symmTensor& t, const scalar lambda)
00399 {
00400     if (mag(lambda) < SMALL)
00401     {
00402         return vector::zero;
00403     }
00404 
00405     // Construct the matrix for the eigenvector problem
00406     symmTensor A(t - lambda*I);
00407 
00408     // Calculate the sub-determinants of the 3 components
00409     scalar sd0 = A.yy()*A.zz() - A.yz()*A.yz();
00410     scalar sd1 = A.xx()*A.zz() - A.xz()*A.xz();
00411     scalar sd2 = A.xx()*A.yy() - A.xy()*A.xy();
00412 
00413     scalar magSd0 = mag(sd0);
00414     scalar magSd1 = mag(sd1);
00415     scalar magSd2 = mag(sd2);
00416 
00417     // Evaluate the eigenvector using the largest sub-determinant
00418     if (magSd0 > magSd1 && magSd0 > magSd2 && magSd0 > SMALL)
00419     {
00420         vector ev
00421         (
00422             1,
00423             (A.yz()*A.xz() - A.zz()*A.xy())/sd0,
00424             (A.yz()*A.xy() - A.yy()*A.xz())/sd0
00425         );
00426         ev /= mag(ev);
00427 
00428         return ev;
00429     }
00430     else if (magSd1 > magSd2 && magSd1 > SMALL)
00431     {
00432         vector ev
00433         (
00434             (A.xz()*A.yz() - A.zz()*A.xy())/sd1,
00435             1,
00436             (A.xz()*A.xy() - A.xx()*A.yz())/sd1
00437         );
00438         ev /= mag(ev);
00439 
00440         return ev;
00441     }
00442     else if (magSd2 > SMALL)
00443     {
00444         vector ev
00445         (
00446             (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
00447             (A.xy()*A.xz() - A.xx()*A.yz())/sd2,
00448             1
00449         );
00450         ev /= mag(ev);
00451 
00452         return ev;
00453     }
00454     else
00455     {
00456         return vector::zero;
00457     }
00458 }
00459 
00460 
00461 tensor eigenVectors(const symmTensor& t)
00462 {
00463     vector evals(eigenValues(t));
00464 
00465     tensor evs
00466     (
00467         eigenVector(t, evals.x()),
00468         eigenVector(t, evals.y()),
00469         eigenVector(t, evals.z())
00470     );
00471 
00472     return evs;
00473 }
00474 
00475 
00476 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00477 
00478 } // End namespace Foam
00479 
00480 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines