BUG: mantis #1373: eigen value/vector calculation now handles repeated eigen values

applications/test/tensor/Test-tensor.C

@@ -69,5 +69,121 @@ int main()
     Info<< (symm(t7) && t7) - (0.5*(t7 + t7.T()) && t7) << endl;
     Info<< (t7 && symm(t7)) - (t7 && 0.5*(t7 + t7.T())) << endl;
 
+
+    /*
+    // Lots of awkward eigenvector tests ...
+
+    tensor T_rand_real
+    (
+        0.9999996423721313, 0.3330855667591095, 0.6646450161933899,
+        0.9745196104049683, 0.0369445420801640, 0.0846728682518005,
+        0.6474838852882385, 0.1617118716239929, 0.2041363865137100
+    );
+    Debug(T_rand_real);
+    vector L_rand_real(eigenValues(T_rand_real));
+    Debug(L_rand_real);
+    tensor U_rand_real(eigenVectors(T_rand_real));
+    Debug(U_rand_real);
+
+    Info << endl << endl;
+
+    tensor T_rand_imag
+    (
+        0.8668024539947510, 0.1664607226848602, 0.8925783634185791,
+        0.9126510620117188, 0.7408077120780945, 0.1499115079641342,
+        0.0936608463525772, 0.7615650296211243, 0.8953040242195129
+    );
+    Debug(T_rand_imag);
+    vector L_rand_imag(eigenValues(T_rand_imag));
+    Debug(L_rand_imag);
+    tensor U_rand_imag(eigenVectors(T_rand_imag));
+    Debug(U_rand_imag);
+
+    Info << endl << endl;
+
+    tensor T_rand_symm
+    (
+        1.9999992847442627, 1.3076051771640778, 1.3121289014816284,
+        1.3076051771640778, 0.0738890841603279, 0.2463847398757935,
+        1.3121289014816284, 0.2463847398757935, 0.4082727730274200
+    );
+    Debug(T_rand_symm);
+    vector L_rand_symm(eigenValues(T_rand_symm));
+    Debug(L_rand_symm);
+    tensor U_rand_symm(eigenVectors(T_rand_symm));
+    Debug(U_rand_symm);
+
+    Info << endl << endl;
+
+    symmTensor T_rand_Symm
+    (
+        1.9999992847442627, 1.3076051771640778, 1.3121289014816284,
+                            0.0738890841603279, 0.2463847398757935,
+                                                0.4082727730274200
+    );
+    Debug(T_rand_Symm);
+    vector L_rand_Symm(eigenValues(T_rand_Symm));
+    Debug(L_rand_Symm);
+    tensor U_rand_Symm(eigenVectors(T_rand_Symm));
+    Debug(U_rand_Symm);
+
+    Info << endl << endl;
+
+    tensor T_rand_diag
+    (
+        0.8668024539947510, 0, 0,
+        0, 0.7408077120780945, 0,
+        0, 0, 0.8953040242195129
+    );
+    Debug(T_rand_diag);
+    vector L_rand_diag(eigenValues(T_rand_diag));
+    Debug(L_rand_diag);
+    tensor U_rand_diag(eigenVectors(T_rand_diag));
+    Debug(U_rand_diag);
+
+    Info << endl << endl;
+
+    tensor T_repeated
+    (
+        0, 1, 1,
+        1, 0, 1,
+        1, 1, 0
+    );
+    Debug(T_repeated);
+    vector L_repeated(eigenValues(T_repeated));
+    Debug(L_repeated);
+    tensor U_repeated(eigenVectors(T_repeated));
+    Debug(U_repeated);
+
+    Info << endl << endl;
+
+    tensor T_repeated_zero
+    (
+        1, 1, 1,
+        1, 1, 1,
+        1, 1, 1
+    );
+    Debug(T_repeated_zero);
+    vector L_repeated_zero(eigenValues(T_repeated_zero));
+    Debug(L_repeated_zero);
+    tensor U_repeated_zero(eigenVectors(T_repeated_zero));
+    Debug(U_repeated_zero);
+
+    Info << endl << endl;
+
+    tensor T_triple
+    (
+        2, 0, 0,
+        0, 2, 0,
+        0, 0, 2
+    );
+    Debug(T_triple);
+    vector L_triple(eigenValues(T_triple));
+    Debug(L_triple);
+    tensor U_triple(eigenVectors(T_triple));
+    Debug(U_triple);
+    */
+
+
     return 0;
 }

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

@@ -89,99 +89,80 @@ namespace Foam
 
 Foam::vector Foam::eigenValues(const tensor& t)
 {
-    scalar i = 0;
+    // The eigenvalues
-    scalar ii = 0;
+    scalar i, ii, iii;
-    scalar iii = 0;
 
-    if
+    // Coefficients of the characteristic polynmial
-    (
+    // x^3 + a*x^2 + b*x + c = 0
-        (
+    scalar a =
-            mag(t.xy()) + mag(t.xz()) + mag(t.yx())
+       - t.xx() - t.yy() - t.zz();
-          + mag(t.yz()) + mag(t.zx()) + mag(t.zy())
-        )
-      < SMALL
-    )
-    {
-        // diagonal matrix
-        i = t.xx();
-        ii = t.yy();
-        iii = t.zz();
-    }
-    else
-    {
-        scalar a = -t.xx() - t.yy() - t.zz();
 
-        scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
+    scalar b =
-            - t.xy()*t.yx() - t.xz()*t.zx() - t.yz()*t.zy();
+        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()
+    scalar c =
-            - t.xz()*t.yx()*t.zy() + t.xz()*t.yy()*t.zx()
+      - t.xx()*t.yy()*t.zz()
-            + t.xy()*t.yx()*t.zz() + t.xx()*t.yz()*t.zy();
+      - 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();
 
-        // If there is a zero root
+    // Auxillary variables
-        if (mag(c) < 1e-100)
-        {
-            scalar disc = sqr(a) - 4*b;
-
-            if (disc >= -SMALL)
-            {
-                scalar q = -0.5*sqrt(max(0.0, disc));
-
-                i = 0;
-                ii = -0.5*a + q;
-                iii = -0.5*a - q;
-            }
-            else
-            {
-                FatalErrorIn("eigenValues(const tensor&)")
-                    << "zero and complex eigenvalues in tensor: " << t
-                    << abort(FatalError);
-            }
-        }
-        else
-        {
-            scalar Q = (a*a - 3*b)/9;
-            scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
-
-            scalar R2 = sqr(R);
-            scalar Q3 = pow3(Q);
-
-            // Three different real roots
-            if (R2 < Q3)
-            {
-                scalar sqrtQ = sqrt(Q);
-                scalar theta = acos(min(1.0, max(-1.0, R/(Q*sqrtQ))));
-
-                scalar m2SqrtQ = -2*sqrtQ;
     scalar aBy3 = a/3;
 
-                i = m2SqrtQ*cos(theta/3) - aBy3;
+    scalar P = (a*a - 3*b)/9; // == -p_wikipedia/3
-                ii = m2SqrtQ*cos((theta + twoPi)/3) - aBy3;
+    scalar PPP = P*P*P;
-                iii = m2SqrtQ*cos((theta - twoPi)/3) - aBy3;
-            }
-            else
-            {
-                scalar A = cbrt(R + sqrt(R2 - Q3));
 
-                // Three equal real roots
+    scalar Q = (2*a*a*a - 9*a*b + 27*c)/54; // == q_wikipedia/2
-                if (A < SMALL)
+    scalar QQ = Q*Q;
+
+    // Three identical roots
+    if (mag(P) < SMALL && mag(Q) < SMALL)
     {
-                    scalar root = -a/3;
+        return vector(- aBy3, - aBy3, - aBy3);
-                    return vector(root, root, root);
     }
+
+    // Two identical roots and one distinct root
+    else if (mag(PPP/QQ - 1) < SMALL)
+    {
+        scalar sqrtP = sqrt(P);
+        scalar signQ = sign(Q);
+
+        i = ii = signQ*sqrtP - aBy3;
+        iii = - 2*signQ*sqrtP - aBy3;
+    }
+
+    // Three distinct roots
+    else if (PPP > QQ)
+    {
+        scalar sqrtP = sqrt(P);
+        scalar value = cos(acos(Q/sqrt(PPP))/3);
+        scalar delta = sqrt(3 - 3*value*value);
+
+        i = - 2*sqrtP*value - aBy3;
+        ii = sqrtP*(value + delta) - aBy3;
+        iii = sqrtP*(value - delta) - aBy3;
+    }
+
+    // One real root, two imaginary roots
+    // based on the above logic, PPP must be less than QQ
     else
     {
-                    // Complex roots
         WarningIn("eigenValues(const tensor&)")
             << "complex eigenvalues detected for tensor: " << t
             << endl;
 
-                    return vector::zero;
+        if (mag(P) < SMALL)
-                }
+        {
-            }
+            i = cbrt(QQ/2);
         }
+        else
+        {
+            scalar w = cbrt(- Q - sqrt(QQ - PPP));
+            i = w + P/w - aBy3;
         }
 
+        return vector(-VGREAT, i, VGREAT);
+    }
 
     // Sort the eigenvalues into ascending order
     if (i > ii)
@@ -203,24 +184,35 @@ Foam::vector Foam::eigenValues(const tensor& t)
 }
 
 
-Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
+Foam::vector Foam::eigenVector
+(
+    const tensor& t,
+    const scalar lambda
+)
 {
-    if (mag(lambda) < SMALL)
+    // Constantly rotating direction ensures different eigenvectors are
-    {
+    // generated when called sequentially with a multiple eigenvalue
-        return vector::zero;
+    static vector direction(0,0,1);
-    }
+    vector oldDirection(direction);
+    scalar temp = direction[2];
+    direction[2] = direction[1];
+    direction[1] = direction[0];
+    direction[0] = temp;
 
-    // Construct the matrix for the eigenvector problem
+    // Construct the linear system for this eigenvalue
     tensor A(t - lambda*I);
 
-    // Calculate the sub-determinants of the 3 components
+    // Determinants of the 2x2 sub-matrices used to find the eigenvectors
-    scalar sd0 = A.yy()*A.zz() - A.yz()*A.zy();
+    scalar sd0, sd1, sd2;
-    scalar sd1 = A.xx()*A.zz() - A.xz()*A.zx();
+    scalar magSd0, magSd1, magSd2;
-    scalar sd2 = A.xx()*A.yy() - A.xy()*A.yx();
 
-    scalar magSd0 = mag(sd0);
+    // Sub-determinants for a unique eivenvalue
-    scalar magSd1 = mag(sd1);
+    sd0 = A.yy()*A.zz() - A.yz()*A.zy();
-    scalar magSd2 = mag(sd2);
+    sd1 = A.zz()*A.xx() - A.zx()*A.xz();
+    sd2 = A.xx()*A.yy() - A.xy()*A.yx();
+    magSd0 = mag(sd0);
+    magSd1 = mag(sd1);
+    magSd2 = mag(sd2);
 
     // Evaluate the eigenvector using the largest sub-determinant
     if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
@@ -231,9 +223,8 @@ Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
             (A.yz()*A.zx() - A.zz()*A.yx())/sd0,
             (A.zy()*A.yx() - A.yy()*A.zx())/sd0
         );
-        ev /= mag(ev);
 
-        return ev;
+        return ev/mag(ev);
     }
     else if (magSd1 >= magSd2 && magSd1 > SMALL)
     {
@@ -243,9 +234,8 @@ Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
             1,
             (A.zx()*A.xy() - A.xx()*A.zy())/sd1
         );
-        ev /= mag(ev);
 
-        return ev;
+        return ev/mag(ev);
     }
     else if (magSd2 > SMALL)
     {
@@ -255,14 +245,55 @@ Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
             (A.yx()*A.xz() - A.xx()*A.yz())/sd2,
             1
         );
-        ev /= mag(ev);
 
-        return ev;
+        return ev/mag(ev);
     }
-    else
+
+    // 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();
+    magSd0 = mag(sd0);
+    magSd1 = mag(sd1);
+    magSd2 = mag(sd2);
+
+    // Evaluate the eigenvector using the largest sub-determinant
+    if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
     {
-        return vector::zero;
+        vector ev
+        (
+            1,
+            (A.yz()*direction.x() - direction.z()*A.yx())/sd0,
+            (direction.y()*A.yx() - A.yy()*direction.x())/sd0
+        );
+
+        return ev/mag(ev);
     }
+    else if (magSd1 >= magSd2 && magSd1 > SMALL)
+    {
+        vector ev
+        (
+            (direction.z()*A.zy() - A.zz()*direction.y())/sd1,
+            1,
+            (A.zx()*direction.y() - direction.x()*A.zy())/sd1
+        );
+
+        return ev/mag(ev);
+    }
+    else if (magSd2 > SMALL)
+    {
+        vector ev
+        (
+            (A.xy()*direction.z() - direction.y()*A.xz())/sd2,
+            (direction.x()*A.xz() - A.xx()*direction.z())/sd2,
+            1
+        );
+
+        return ev/mag(ev);
+    }
+
+    // Triple eigenvalue
+    return oldDirection;
 }
 
 
@@ -281,198 +312,21 @@ Foam::tensor Foam::eigenVectors(const tensor& t)
 }
 
 
-// Return eigenvalues in ascending order of absolute values
 Foam::vector Foam::eigenValues(const symmTensor& t)
 {
-    scalar i = 0;
+    return eigenValues(tensor(t));
-    scalar ii = 0;
-    scalar iii = 0;
-
-    if
-    (
-        (
-            mag(t.xy()) + mag(t.xz()) + mag(t.xy())
-          + mag(t.yz()) + mag(t.xz()) + mag(t.yz())
-        )
-      < SMALL
-    )
-    {
-        // diagonal matrix
-        i = t.xx();
-        ii = t.yy();
-        iii = t.zz();
-    }
-    else
-    {
-        scalar a = -t.xx() - t.yy() - t.zz();
-
-        scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
-            - t.xy()*t.xy() - t.xz()*t.xz() - t.yz()*t.yz();
-
-        scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.xz()
-            - t.xz()*t.xy()*t.yz() + t.xz()*t.yy()*t.xz()
-            + t.xy()*t.xy()*t.zz() + t.xx()*t.yz()*t.yz();
-
-        // If there is a zero root
-        if (mag(c) < 1e-100)
-        {
-            scalar disc = sqr(a) - 4*b;
-
-            if (disc >= -SMALL)
-            {
-                scalar q = -0.5*sqrt(max(0.0, disc));
-
-                i = 0;
-                ii = -0.5*a + q;
-                iii = -0.5*a - q;
-            }
-            else
-            {
-                FatalErrorIn("eigenValues(const tensor&)")
-                    << "zero and complex eigenvalues in tensor: " << t
-                    << abort(FatalError);
-            }
-        }
-        else
-        {
-            scalar Q = (a*a - 3*b)/9;
-            scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
-
-            scalar R2 = sqr(R);
-            scalar Q3 = pow3(Q);
-
-            // Three different real roots
-            if (R2 < Q3)
-            {
-                scalar sqrtQ = sqrt(Q);
-                scalar theta = acos(min(1.0, max(-1.0, R/(Q*sqrtQ))));
-
-                scalar m2SqrtQ = -2*sqrtQ;
-                scalar aBy3 = a/3;
-
-                i = m2SqrtQ*cos(theta/3) - aBy3;
-                ii = m2SqrtQ*cos((theta + twoPi)/3) - aBy3;
-                iii = m2SqrtQ*cos((theta - twoPi)/3) - aBy3;
-            }
-            else
-            {
-                scalar A = cbrt(R + sqrt(R2 - Q3));
-
-                // Three equal real roots
-                if (A < SMALL)
-                {
-                    scalar root = -a/3;
-                    return vector(root, root, root);
-                }
-                else
-                {
-                    // Complex roots
-                    WarningIn("eigenValues(const symmTensor&)")
-                        << "complex eigenvalues detected for symmTensor: " << t
-                        << endl;
-
-                    return vector::zero;
-                }
-            }
-        }
-    }
-
-
-    // Sort the eigenvalues into ascending order
-    if (i > ii)
-    {
-        Swap(i, ii);
-    }
-
-    if (ii > iii)
-    {
-        Swap(ii, iii);
-    }
-
-    if (i > ii)
-    {
-        Swap(i, ii);
-    }
-
-    return vector(i, ii, iii);
 }
 
 
 Foam::vector Foam::eigenVector(const symmTensor& t, const scalar lambda)
 {
-    if (mag(lambda) < SMALL)
+    return eigenVector(tensor(t), lambda);
-    {
-        return vector::zero;
-    }
-
-    // Construct the matrix for the eigenvector problem
-    symmTensor A(t - lambda*I);
-
-    // Calculate the sub-determinants of the 3 components
-    scalar sd0 = A.yy()*A.zz() - A.yz()*A.yz();
-    scalar sd1 = A.xx()*A.zz() - A.xz()*A.xz();
-    scalar sd2 = A.xx()*A.yy() - A.xy()*A.xy();
-
-    scalar magSd0 = mag(sd0);
-    scalar magSd1 = mag(sd1);
-    scalar magSd2 = mag(sd2);
-
-    // Evaluate the eigenvector using the largest sub-determinant
-    if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
-    {
-        vector ev
-        (
-            1,
-            (A.yz()*A.xz() - A.zz()*A.xy())/sd0,
-            (A.yz()*A.xy() - A.yy()*A.xz())/sd0
-        );
-        ev /= mag(ev);
-
-        return ev;
-    }
-    else if (magSd1 >= magSd2 && magSd1 > SMALL)
-    {
-        vector ev
-        (
-            (A.xz()*A.yz() - A.zz()*A.xy())/sd1,
-            1,
-            (A.xz()*A.xy() - A.xx()*A.yz())/sd1
-        );
-        ev /= mag(ev);
-
-        return ev;
-    }
-    else if (magSd2 > SMALL)
-    {
-        vector ev
-        (
-            (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
-            (A.xy()*A.xz() - A.xx()*A.yz())/sd2,
-            1
-        );
-        ev /= mag(ev);
-
-        return ev;
-    }
-    else
-    {
-        return vector::zero;
-    }
 }
 
 
 Foam::tensor Foam::eigenVectors(const symmTensor& t)
 {
-    vector evals(eigenValues(t));
+    return eigenVectors(tensor(t));
-
-    tensor evs
-    (
-        eigenVector(t, evals.x()),
-        eigenVector(t, evals.y()),
-        eigenVector(t, evals.z())
-    );
-
-    return evs;
 }