ENH: Migrated eigen vector/value updates from old development line and updated dependencies

src/OpenFOAM/primitives/Tensor/tensor/tensor.C

@@ -70,24 +70,30 @@ const Foam::tensor Foam::tensor::I
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-Foam::vector Foam::eigenValues(const tensor& t)
+Foam::vector Foam::eigenValues(const tensor& T)
 {
     // The eigenvalues
     scalar i, ii, iii;
 
     // diagonal matrix
-    if
-    (
+    const scalar onDiagMagSum =
         (
-            mag(t.xy()) + mag(t.xz()) + mag(t.yx())
-            + mag(t.yz()) + mag(t.zx()) + mag(t.zy())
-        )
-        < SMALL
-    )
+            mag(T.xx()) + mag(T.yy()) + mag(T.zz())
+        );
+
+    const scalar offDiagMagSum =
+        (
+            mag(T.xy()) + mag(T.xz()) + mag(T.yx())
+          + mag(T.yz()) + mag(T.zx()) + mag(T.zy())
+        );
+
+    const scalar magSum = onDiagMagSum + offDiagMagSum;
+
+    if (offDiagMagSum < max(VSMALL, SMALL*magSum))
     {
-        i = t.xx();
-        ii = t.yy();
-        iii = t.zz();
+        i = T.xx();
+        ii = T.yy();
+        iii = T.zz();
     }
 
     // non-diagonal matrix
@@ -96,16 +102,16 @@ Foam::vector Foam::eigenValues(const tensor& t)
         // Coefficients of the characteristic polynmial
         // x^3 + a*x^2 + b*x + c = 0
         scalar a =
-           - t.xx() - t.yy() - t.zz();
+           - T.xx() - T.yy() - T.zz();
 
         scalar b =
-            t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
-          - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz();
+            T.xx()*T.yy() + T.xx()*T.zz() + T.yy()*T.zz()
+          - T.xy()*T.yx() - T.yz()*T.zy() - T.zx()*T.xz();
 
         scalar c =
-          - t.xx()*t.yy()*t.zz()
-          - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx()
-          + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx();
+          - T.xx()*T.yy()*T.zz()
+          - T.xy()*T.yz()*T.zx() - T.xz()*T.zy()*T.yx()
+          + T.xx()*T.yz()*T.zy() + T.yy()*T.zx()*T.xz() + T.zz()*T.xy()*T.yx();
 
         // Auxillary variables
         scalar aBy3 = a/3;
@@ -117,13 +123,13 @@ Foam::vector Foam::eigenValues(const tensor& t)
         scalar QQ = Q*Q;
 
         // Three identical roots
-        if (mag(P) < SMALL && mag(Q) < SMALL)
+        if (mag(P) < SMALL*sqr(magSum) && mag(Q) < SMALL*pow3(magSum))
         {
             return vector(- aBy3, - aBy3, - aBy3);
         }
 
         // Two identical roots and one distinct root
-        else if (mag(PPP/QQ - 1) < SMALL)
+        else if (mag(PPP - QQ) < SMALL*pow6(magSum))
         {
             scalar sqrtP = sqrt(P);
             scalar signQ = sign(Q);
@@ -148,11 +154,11 @@ Foam::vector Foam::eigenValues(const tensor& t)
         // based on the above logic, PPP must be less than QQ
         else
         {
-            WarningInFunction
-                << "complex eigenvalues detected for tensor: " << t
+            WarningIn("eigenValues(const tensor&)")
+                << "complex eigenvalues detected for tensor: " << T
                 << endl;
 
-            if (mag(P) < SMALL)
+            if (mag(P) < SMALL*sqr(magSum))
             {
                 i = cbrt(QQ/2);
             }
@@ -188,21 +194,14 @@ Foam::vector Foam::eigenValues(const tensor& t)
 
 Foam::vector Foam::eigenVector
 (
-    const tensor& t,
-    const scalar lambda
+    const tensor& T,
+    const scalar lambda,
+    const vector& direction1,
+    const vector& direction2
 )
 {
-    // Constantly rotating direction ensures different eigenvectors are
-    // generated when called sequentially with a multiple eigenvalue
-    static vector direction(1,0,0);
-    vector oldDirection(direction);
-    scalar temp = direction[2];
-    direction[2] = direction[1];
-    direction[1] = direction[0];
-    direction[0] = temp;
-
     // Construct the linear system for this eigenvalue
-    tensor A(t - lambda*I);
+    tensor A(T - lambda*I);
 
     // Determinants of the 2x2 sub-matrices used to find the eigenvectors
     scalar sd0, sd1, sd2;
@@ -252,9 +251,9 @@ Foam::vector Foam::eigenVector
     }
 
     // Sub-determinants for a repeated eigenvalue
-    sd0 = A.yy()*direction.z() - A.yz()*direction.y();
-    sd1 = A.zz()*direction.x() - A.zx()*direction.z();
-    sd2 = A.xx()*direction.y() - A.xy()*direction.x();
+    sd0 = A.yy()*direction1.z() - A.yz()*direction1.y();
+    sd1 = A.zz()*direction1.x() - A.zx()*direction1.z();
+    sd2 = A.xx()*direction1.y() - A.xy()*direction1.x();
     magSd0 = mag(sd0);
     magSd1 = mag(sd1);
     magSd2 = mag(sd2);
@@ -265,8 +264,8 @@ Foam::vector Foam::eigenVector
         vector ev
         (
             1,
-            (A.yz()*direction.x() - direction.z()*A.yx())/sd0,
-            (direction.y()*A.yx() - A.yy()*direction.x())/sd0
+            (A.yz()*direction1.x() - direction1.z()*A.yx())/sd0,
+            (direction1.y()*A.yx() - A.yy()*direction1.x())/sd0
         );
 
         return ev/mag(ev);
@@ -275,9 +274,9 @@ Foam::vector Foam::eigenVector
     {
         vector ev
         (
-            (direction.z()*A.zy() - A.zz()*direction.y())/sd1,
+            (direction1.z()*A.zy() - A.zz()*direction1.y())/sd1,
             1,
-            (A.zx()*direction.y() - direction.x()*A.zy())/sd1
+            (A.zx()*direction1.y() - direction1.x()*A.zy())/sd1
         );
 
         return ev/mag(ev);
@@ -286,8 +285,8 @@ Foam::vector Foam::eigenVector
     {
         vector ev
         (
-            (A.xy()*direction.z() - direction.y()*A.xz())/sd2,
-            (direction.x()*A.xz() - A.xx()*direction.z())/sd2,
+            (A.xy()*direction1.z() - direction1.y()*A.xz())/sd2,
+            (direction1.x()*A.xz() - A.xx()*direction1.z())/sd2,
             1
         );
 
@@ -295,40 +294,57 @@ Foam::vector Foam::eigenVector
     }
 
     // Triple eigenvalue
-    return oldDirection;
+    return direction1^direction2;
 }
 
 
-Foam::tensor Foam::eigenVectors(const tensor& t)
+Foam::tensor Foam::eigenVectors(const tensor& T, const vector& lambdas)
 {
-    vector evals(eigenValues(t));
+    vector Ux(1, 0, 0), Uy(0, 1, 0), Uz(0, 0, 1);
 
-    tensor evs
-    (
-        eigenVector(t, evals.x()),
-        eigenVector(t, evals.y()),
-        eigenVector(t, evals.z())
-    );
+    Ux = eigenVector(T, lambdas.x(), Uy, Uz);
+    Uy = eigenVector(T, lambdas.y(), Uz, Ux);
+    Uz = eigenVector(T, lambdas.z(), Ux, Uy);
 
-    return evs;
+    return tensor(Ux, Uy, Uz);
 }
 
 
-Foam::vector Foam::eigenValues(const symmTensor& t)
+Foam::tensor Foam::eigenVectors(const tensor& T)
 {
-    return eigenValues(tensor(t));
+    const vector lambdas(eigenValues(T));
+
+    return eigenVectors(T, lambdas);
 }
 
 
-Foam::vector Foam::eigenVector(const symmTensor& t, const scalar lambda)
+Foam::vector Foam::eigenValues(const symmTensor& T)
 {
-    return eigenVector(tensor(t), lambda);
+    return eigenValues(tensor(T));
 }
 
 
-Foam::tensor Foam::eigenVectors(const symmTensor& t)
+Foam::vector Foam::eigenVector
+(
+    const symmTensor& T,
+    const scalar lambda,
+    const vector& direction1,
+    const vector& direction2
+)
 {
-    return eigenVectors(tensor(t));
+    return eigenVector(tensor(T), lambda, direction1, direction2);
+}
+
+
+Foam::tensor Foam::eigenVectors(const symmTensor& T, const vector& lambdas)
+{
+    return eigenVectors(tensor(T), lambdas);
+}
+
+
+Foam::tensor Foam::eigenVectors(const symmTensor& T)
+{
+    return eigenVectors(tensor(T));
 }
 
 

src/OpenFOAM/primitives/Tensor/tensor/tensor.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,13 +50,27 @@ namespace Foam
 
 typedef Tensor<scalar> tensor;
 
-vector eigenValues(const tensor&);
-vector eigenVector(const tensor&, const scalar lambda);
-tensor eigenVectors(const tensor&);
+vector eigenValues(const tensor& T);
+vector eigenVector
+(
+    const tensor& T,
+    const scalar lambda,
+    const vector& direction1,
+    const vector& direction2
+);
+tensor eigenVectors(const tensor& T, const vector& lambdas);
+tensor eigenVectors(const tensor& T);
 
-vector eigenValues(const symmTensor&);
-vector eigenVector(const symmTensor&, const scalar lambda);
-tensor eigenVectors(const symmTensor&);
+vector eigenValues(const symmTensor& T);
+vector eigenVector
+(
+    const symmTensor& T,
+    const scalar lambda,
+    const vector& direction1,
+    const vector& direction2
+);
+tensor eigenVectors(const symmTensor& T, const vector& lambdas);
+tensor eigenVectors(const symmTensor& T);
 
 //- Data associated with tensor type are contiguous
 template<>

src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C

@@ -458,7 +458,7 @@ void Foam::edgeCollapser::faceCollapseAxisAndAspectRatio
             // normal, as it has the greatest value.  The minimum eigenvalue
             // is the dominant collapse axis for high aspect ratio faces.
 
-            collapseAxis = eigenVector(J, eVals.x());
+            collapseAxis = eigenVectors(J, eVals).x();
 
             // The inertia calculation describes the mass distribution as a
             // function of distance squared to the axis, so the square root of