ENH: triangle and tetrahedron robustness improvements.

Currently very verbose when a degenerate shape is encountered.

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H

@@ -156,29 +156,12 @@ inline Point tetrahedron<Point, PointRef>::circumCentre() const
 
     if (Foam::mag(denom) < ROOTVSMALL)
     {
-        // Degenerate tet. Use test of the individual triangles.
+        WarningIn("Point tetrahedron<Point, PointRef>::circumCentre() const")
-        {
+            << "Degenerate tetrahedron:" << nl << *this << nl
-            point triCentre = triPointRef(a_, b_, c_).circumCentre();
+            <<"Returning centre instead of circumCentre."
-            if (magSqr(d_ - triCentre) < magSqr(a_ - triCentre))
+            << endl;
-            {
+
-                return triCentre;
+        return centre();
-            }
-        }
-        {
-            point triCentre = triPointRef(a_, b_, d_).circumCentre();
-            if (magSqr(c_ - triCentre) < magSqr(a_ - triCentre))
-            {
-                return triCentre;
-            }
-        }
-        {
-            point triCentre = triPointRef(a_, c_, d_).circumCentre();
-            if (magSqr(b_ - triCentre) < magSqr(a_ - triCentre))
-            {
-                return triCentre;
-            }
-        }
-        return triPointRef(b_, c_, d_).circumCentre();
     }
 
     return a_ + 0.5*(a + num/denom);
@@ -188,16 +171,43 @@ inline Point tetrahedron<Point, PointRef>::circumCentre() const
 template<class Point, class PointRef>
 inline scalar tetrahedron<Point, PointRef>::circumRadius() const
 {
-    return Foam::mag(a_ - circumCentre());
+    vector a = b_ - a_;
+    vector b = c_ - a_;
+    vector c = d_ - a_;
+
+    scalar lambda = magSqr(c) - (a & c);
+    scalar mu = magSqr(b) - (a & b);
+
+    vector ba = b ^ a;
+    vector ca = c ^ a;
+
+    vector num = lambda*ba - mu*ca;
+    scalar denom = (c & ba);
+
+    if (Foam::mag(denom) < ROOTVSMALL)
+    {
+        WarningIn("Point tetrahedron<Point, PointRef>::circumCentre() const")
+            << "Degenerate tetrahedron:" << nl << *this << nl
+            << "Returning GREAT for circumRadius."
+            << endl;
+
+        return GREAT;
+    }
+
+    return Foam::mag(0.5*(a + num/denom));
 }
 
 
 template<class Point, class PointRef>
 inline scalar tetrahedron<Point, PointRef>::quality() const
 {
-    // Note: 8/(9*sqrt(3)) = 0.5132002393
+    return
-
+        mag()
-    return mag()/(0.5132002393*pow3(min(circumRadius(), GREAT)) + ROOTVSMALL);
+       /(
+            8.0/(9.0*sqrt(3.0))
+           *pow3(min(circumRadius(), GREAT))
+          + ROOTVSMALL
+        );
 }
 
 
@@ -266,7 +276,8 @@ scalar tetrahedron<Point, PointRef>::barycentric
                 "const point& pt"
             ") const"
         )
-            << "Degenerate tetrahedron - returning 1/4 barycentric coordinates."
+            << "Degenerate tetrahedron:" << nl << *this << nl
+            << "Returning 1/4 barycentric coordinates."
             << endl;
 
         bary = List<scalar>(4, 0.25);

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -109,9 +109,9 @@ inline vector triangle<Point, PointRef>::normal() const
 template<class Point, class PointRef>
 inline Point triangle<Point, PointRef>::circumCentre() const
 {
-    scalar d1 = (c_ - a_)&(b_ - a_);
+    scalar d1 =  (c_ - a_) & (b_ - a_);
-    scalar d2 = -(c_ - b_)&(b_ - a_);
+    scalar d2 = -(c_ - b_) & (b_ - a_);
-    scalar d3 = (c_ - a_)&(c_ - b_);
+    scalar d3 =  (c_ - a_) & (c_ - b_);
 
     scalar c1 = d2*d3;
     scalar c2 = d3*d1;
@@ -119,6 +119,16 @@ inline Point triangle<Point, PointRef>::circumCentre() const
 
     scalar c = c1 + c2 + c3;
 
+    if (Foam::mag(c) < ROOTVSMALL)
+    {
+        WarningIn("Point triangle<Point, PointRef>::circumCentre() const")
+            << "Degenerate triangle:" << nl << *this << nl
+            << "Returning centre instead of circumCentre."
+            << endl;
+
+        return centre();
+    }
+
     return
     (
         ((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*c)
@@ -129,14 +139,19 @@ inline Point triangle<Point, PointRef>::circumCentre() const
 template<class Point, class PointRef>
 inline scalar triangle<Point, PointRef>::circumRadius() const
 {
-    scalar d1 = (c_ - a_) & (b_ - a_);
+    scalar d1 =  (c_ - a_) & (b_ - a_);
-    scalar d2 = - (c_ - b_) & (b_ - a_);
+    scalar d2 = -(c_ - b_) & (b_ - a_);
-    scalar d3 = (c_ - a_) & (c_ - b_);
+    scalar d3 =  (c_ - a_) & (c_ - b_);
 
     scalar denom = d2*d3 + d3*d1 + d1*d2;
 
     if (Foam::mag(denom) < VSMALL)
     {
+        WarningIn("scalar triangle<Point, PointRef>::circumRadius() const")
+            << "Degenerate triangle:" << nl << *this << nl
+            << "Returning GREAT for circumRadius."
+            << endl;
+
         return GREAT;
     }
     else
@@ -151,16 +166,7 @@ inline scalar triangle<Point, PointRef>::circumRadius() const
 template<class Point, class PointRef>
 inline scalar triangle<Point, PointRef>::quality() const
 {
-    // Note: 3*sqr(3)/(4*pi) = 0.4134966716
+    return mag()/(Foam::sqr(circumRadius())*3.0*sqrt(3.0)/4.0 + VSMALL);
-
-    return
-        mag()
-      / (
-            constant::mathematical::pi
-           *Foam::sqr(circumRadius())
-           *0.4134966716
-          + VSMALL
-        );
 }
 
 
@@ -264,7 +270,8 @@ scalar triangle<Point, PointRef>::barycentric
                 "const point& pt"
             ") const"
         )
-            << "Degenerate triangle - returning 1/3 barycentric coordinates."
+            << "Degenerate triangle:" << nl << *this << nl
+            << "Returning 1/3 barycentric coordinates."
             << endl;
 
         bary = List<scalar>(3, 1.0/3.0);
@@ -490,7 +497,7 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
 ) const
 {
     // Adapted from:
-    // Real-time collision detection, Christer Ericson, 2005, 136-142
+    // Real-time collision detection, Christer Ericson, 2005, p136-142
 
     // Check if P in vertex region outside A
     vector ab = b_ - a_;
@@ -528,6 +535,27 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
 
     if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0)
     {
+        if ((d1 - d3) < ROOTVSMALL)
+        {
+            WarningIn
+            (
+                "pointHit triangle<Point, PointRef>::nearestPointClassify"
+                "("
+                    "const point& p,"
+                    "label& nearType,"
+                    "label& nearLabel"
+                ") const"
+            )
+                << "Degenerate triangle:" << nl << *this << nl
+                << "d1, d3: " << d1 << ", " << d3 << endl;
+
+            // For d1 = d3, a_ and b_ are likely coincident.
+
+            nearType = POINT;
+            nearLabel = 0;
+            return pointHit(false, a_, Foam::mag(a_ - p), true);
+        }
+
         // barycentric coordinates (1-v, v, 0)
         scalar v = d1/(d1 - d3);
 
@@ -556,6 +584,27 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
 
     if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0)
     {
+        if ((d2 - d6) < ROOTVSMALL)
+        {
+            WarningIn
+            (
+                "pointHit triangle<Point, PointRef>::nearestPointClassify"
+                "("
+                    "const point& p,"
+                    "label& nearType,"
+                    "label& nearLabel"
+                ") const"
+            )
+                << "Degenerate triangle:" << nl << *this << nl
+                << "d2, d6: " << d2 << ", " << d6 << endl;
+
+            // For d2 = d6, a_ and c_ are likely coincident.
+
+            nearType = POINT;
+            nearLabel = 0;
+            return pointHit(false, a_, Foam::mag(a_ - p), true);
+        }
+
         // barycentric coordinates (1-w, 0, w)
         scalar w = d2/(d2 - d6);
 
@@ -570,6 +619,28 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
 
     if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0)
     {
+        if (((d4 - d3) + (d5 - d6)) < ROOTVSMALL)
+        {
+            WarningIn
+            (
+                "pointHit triangle<Point, PointRef>::nearestPointClassify"
+                "("
+                    "const point& p,"
+                    "label& nearType,"
+                    "label& nearLabel"
+                ") const"
+            )
+                << "Degenerate triangle:" << nl << *this << nl
+                << "(d4 - d3), (d6 - d5): " << (d4 - d3) << ", " << (d6 - d5)
+                << endl;
+
+            // For (d4 - d3) = (d6 - d5), b_ and c_ are likely coincident.
+
+            nearType = POINT;
+            nearLabel = 1;
+            return pointHit(false, b_, Foam::mag(b_ - p), true);
+        }
+
         // barycentric coordinates (0, 1-w, w)
         scalar w = (d4 - d3)/((d4 - d3) + (d5 - d6));
 
@@ -581,6 +652,28 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
 
     // P inside face region. Compute Q through its barycentric
     // coordinates (u, v, w)
+
+    if ((va + vb + vc) < ROOTVSMALL)
+    {
+        WarningIn
+        (
+            "pointHit triangle<Point, PointRef>::nearestPointClassify"
+            "("
+                "const point& p,"
+                "label& nearType,"
+                "label& nearLabel"
+            ") const"
+        )
+            << "Degenerate triangle:" << nl << *this << nl
+            << "va, vb, vc: " << va << ", " << vb << ", " << vc
+            << endl;
+
+        point nearPt = centre();
+        nearType = NONE,
+        nearLabel = -1;
+        return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
+    }
+
     scalar denom = 1.0/(va + vb + vc);
     scalar v = vb * denom;
     scalar w = vc * denom;