Further tinkering and additional test of the eigenvaues of symmTensors.

applications/test/tensor/tensorTest.C

@@ -38,6 +38,12 @@ int main()
         << (eigenVector(t6, e[2]) & t6) << e[2]*eigenVector(t6, e[2])
         << endl;
 
+    Info<< "Check eigenvalues for symmTensor "
+        << eigenValues(symm(t6)) - eigenValues(tensor(symm(t6))) << endl;
+
+    Info<< "Check eigenvectors for symmTensor "
+        << eigenVectors(symm(t6)) - eigenVectors(tensor(symm(t6))) << endl;
+
     tensor t7(1, 2, 3, 2, 4, 5, 3, 5, 6);
 
     Info<< "Check transformation "
@@ -49,7 +55,7 @@ int main()
         << transform(t1, st1) << endl;
 
     vector v1(1, 2, 3);
-    
+
     Info<< sqr(v1) << endl;
     Info<< symm(t7) << endl;
     Info<< twoSymm(t7) << endl;

src/OpenFOAM/primitives/tensor/tensor.C

@@ -185,7 +185,7 @@ vector eigenValues(const tensor& t)
 
 vector eigenVector(const tensor& t, const scalar lambda)
 {
-    if (lambda < SMALL)
+    if (mag(lambda) < SMALL)
     {
         return vector::zero;
     }
@@ -292,22 +292,22 @@ vector eigenValues(const symmTensor& t)
             + t.xy()*t.xy()*t.zz() + t.xx()*t.yz()*t.yz();
 
         // If there is a zero root
-        if (mag(c) < SMALL)
+        if (mag(c) < 1.0e-100)
         {
             scalar disc = sqr(a) - 4*b;
 
-            if (disc > 0)
+            if (disc >= -SMALL)
             {
-                scalar q = -0.5*(a + sign(a)*sqrt(disc));
+                scalar q = -0.5*sqrt(max(0.0, disc));
 
                 i = 0;
-                ii = q;
-                iii = b/q;
+                ii = -0.5*a + q;
+                iii = -0.5*a - q;
             }
             else
             {
-                FatalErrorIn("eigenValues(const symmTensor&)")
-                    << "zero and complex eigenvalues in symmTensor: " << t
+                FatalErrorIn("eigenValues(const tensor&)")
+                    << "zero and complex eigenvalues in tensor: " << t
                     << abort(FatalError);
             }
         }
@@ -380,7 +380,7 @@ vector eigenValues(const symmTensor& t)
 
 vector eigenVector(const symmTensor& t, const scalar lambda)
 {
-    if (lambda < SMALL)
+    if (mag(lambda) < SMALL)
     {
         return vector::zero;
     }