triangle, tetrahedron: Consolidate circumCircle/Sphere functions

applications/test/triTet/Test-triTet.C

@@ -7,17 +7,20 @@ using namespace Foam;
 
 int main()
 {
-    triangle<point, point> tri
+    const triangle<point, point> tri
     (
         vector(0, 0, 0),
         vector(1, 0, 0),
         vector(1, 1, 0)
     );
 
-    Info<< "tri circumCentre = " << tri.circumCentre() << endl;
+    const Tuple2<point, scalar> triCircle = tri.circumCircle();
-    Info<< "tri circumRadius = " << tri.circumRadius() << endl;
 
-    tetrahedron<point, point> tet
+    Info<< "tri circumCentre = " << triCircle.first() << endl
+        << "tri circumRadius = " << triCircle.second() << endl
+        << "     tri quality = " << tri.quality() << endl;
+
+    const tetrahedron<point, point> tet
     (
         vector(1, 0, 0),
         vector(0, 1, 0),
@@ -25,16 +28,11 @@ int main()
         vector(0.5773502, 0.5773502, 0.5773502)
     );
 
-    Info<< "tet circumCentre = " << tet.circumCentre() << endl;
+    const Tuple2<point, scalar> tetSphere = tet.circumSphere();
-    Info<< "tet circumRadius = " << tet.circumRadius() << endl;
 
-    vector a(Sin);
+    Info<< "tet circumCentre = " << tetSphere.first() << endl
-    vector b(Sin);
+        << "tet circumRadius = " << tetSphere.second() << endl
-    vector c(Sin);
+        << "     tet quality = " << tet.quality() << endl;
-    vector d(Sin);
-
-    Info<< "tet circumRadius = "
-        << tetrahedron<point, point>(a, b, c, d).circumRadius() << endl;
 
     return 0;
 }

applications/utilities/surface/surfaceFeatures/surfaceFeatures.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -616,24 +616,22 @@ namespace Foam
             forAll(surf, fi)
             {
                 const triPointRef& tri = surf[fi].tri(surf.points());
-                const point& triCentre = tri.circumCentre();
 
-                const scalar radiusSqr = min
+                const Tuple2<point, scalar> circle = tri.circumCircle();
-                (
+                const point& c = circle.first();
-                    sqr(4*tri.circumRadius()),
+                const scalar rSqr =
-                    sqr(searchDistance)
+                    min(sqr(4*circle.second()), sqr(searchDistance));
-                );
 
                 pointIndexHitList hitList;
 
-                feMesh.allNearestFeatureEdges(triCentre, radiusSqr, hitList);
+                feMesh.allNearestFeatureEdges(c, rSqr, hitList);
                 featureProximity[fi] = min
                 (
                     feMesh.minDisconnectedDist(hitList),
                     featureProximity[fi]
                 );
 
-                feMesh.allNearestFeaturePoints(triCentre, radiusSqr, hitList);
+                feMesh.allNearestFeaturePoints(c, rSqr, hitList);
                 featureProximity[fi] = min
                 (
                     minDist(hitList),

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

@@ -152,11 +152,8 @@ public:
             //- Return volume
             inline scalar mag() const;
 
-            //- Return circum-centre
+            //- Return the circum centre and radius
-            inline Point circumCentre() const;
+            inline Tuple2<Point, scalar> circumSphere() const;
-
-            //- Return circum-radius
-            inline scalar circumRadius() const;
 
             //- Return quality: Ratio of tetrahedron and circum-sphere
             //  volume, scaled so that a regular tetrahedron has a

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

@@ -174,68 +174,42 @@ inline Foam::scalar Foam::tetrahedron<Point, PointRef>::mag() const
 
 
 template<class Point, class PointRef>
-inline Point Foam::tetrahedron<Point, PointRef>::circumCentre() const
+inline Foam::Tuple2<Point, Foam::scalar>
+Foam::tetrahedron<Point, PointRef>::circumSphere() const
 {
-    vector a = b_ - a_;
+    const vector a = b_ - a_;
-    vector b = c_ - a_;
+    const vector b = c_ - a_;
-    vector c = d_ - a_;
+    const vector c = d_ - a_;
 
-    scalar lambda = magSqr(c) - (a & c);
+    const vector ba = b ^ a;
-    scalar mu = magSqr(b) - (a & b);
+    const vector ca = c ^ a;
 
-    vector ba = b ^ a;
+    const scalar lambda = magSqr(c) - (a & c);
-    vector ca = c ^ a;
+    const scalar mu = magSqr(b) - (a & b);
 
-    vector num = lambda*ba - mu*ca;
+    const vector num = lambda*ba - mu*ca;
-    scalar denom = (c & ba);
+    const scalar denom = (c & ba);
 
     if (Foam::mag(denom) < rootVSmall)
     {
-        // Degenerate tetrahedron, returning centre instead of circumCentre.
+        // Degenerate. Return a point far away.
-
+        static const scalar sqrt3 = sqrt(scalar(3));
-        return centre();
+        return Tuple2<Point, scalar>(point::uniform(great), sqrt3*great);
     }
-
+    else
-    return a_ + 0.5*(a + num/denom);
-}
-
-
-template<class Point, class PointRef>
-inline Foam::scalar Foam::tetrahedron<Point, PointRef>::circumRadius() const
-{
-    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)
     {
-        // Degenerate tetrahedron, returning great for circumRadius.
+        const vector v = (a + num/denom)/2;
-        return great;
+        return Tuple2<Point, scalar>(a_ + v, Foam::mag(v));
     }
-
-    return Foam::mag(0.5*(a + num/denom));
 }
 
 
 template<class Point, class PointRef>
 inline Foam::scalar Foam::tetrahedron<Point, PointRef>::quality() const
 {
-    return
+    const scalar r = circumSphere().second();
-        mag()
+    static const scalar sqrt3 = sqrt(scalar(3));
-       /(
+    return mag()/((8.0/27.0)*sqrt3*pow3(min(r, great)) + rootVSmall);
-            8.0/(9.0*sqrt(3.0))
-           *pow3(min(circumRadius(), great))
-          + rootVSmall
-        );
 }
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -140,11 +140,8 @@ public:
             //- Return unit normal
             inline vector normal() const;
 
-            //- Return circum-centre
+            //- Return the circum centre and radius
-            inline Point circumCentre() const;
+            inline Tuple2<Point, scalar> circumCircle() const;
-
-            //- Return circum-radius
-            inline scalar circumRadius() const;
 
             //- Return quality: Ratio of triangle and circum-circle
             //  area, scaled so that an equilateral triangle has a

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -116,52 +116,33 @@ inline Foam::vector Foam::triangle<Point, PointRef>::normal() const
 
 
 template<class Point, class PointRef>
-inline Point Foam::triangle<Point, PointRef>::circumCentre() const
+inline Foam::Tuple2<Point, Foam::scalar>
+Foam::triangle<Point, PointRef>::circumCircle() const
 {
-    scalar d1 =  (c_ - a_) & (b_ - a_);
+    const scalar d1 =  (c_ - a_) & (b_ - a_);
-    scalar d2 = -(c_ - b_) & (b_ - a_);
+    const scalar d2 = -(c_ - b_) & (b_ - a_);
-    scalar d3 =  (c_ - a_) & (c_ - b_);
+    const scalar d3 =  (c_ - a_) & (c_ - b_);
 
-    scalar c1 = d2*d3;
+    const scalar c1 = d2*d3;
-    scalar c2 = d3*d1;
+    const scalar c2 = d3*d1;
-    scalar c3 = d1*d2;
+    const scalar c3 = d1*d2;
 
-    scalar c = c1 + c2 + c3;
+    const scalar denom = c1 + c2 + c3;
 
-    if (Foam::mag(c) < rootVSmall)
+    if (Foam::mag(denom) < rootVSmall)
     {
-        // Degenerate triangle, returning centre instead of circumCentre.
+        // Degenerate. Return a point far away.
-
+        static const scalar sqrt3 = sqrt(scalar(3));
-        return centre();
+        return Tuple2<Point, scalar>(point::uniform(great), sqrt3*great);
-    }
-
-    return
-    (
-        ((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*c)
-    );
-}
-
-
-template<class Point, class PointRef>
-inline Foam::scalar Foam::triangle<Point, PointRef>::circumRadius() const
-{
-    scalar d1 =  (c_ - a_) & (b_ - a_);
-    scalar d2 = -(c_ - b_) & (b_ - a_);
-    scalar d3 =  (c_ - a_) & (c_ - b_);
-
-    scalar denom = d2*d3 + d3*d1 + d1*d2;
-
-    if (Foam::mag(denom) < vSmall)
-    {
-        // Degenerate triangle, returning great for circumRadius.
-
-        return great;
     }
     else
     {
-        scalar a = (d1 + d2)*(d2 + d3)*(d3 + d1) / denom;
+        return
-
+            Tuple2<Point, scalar>
-        return 0.5*Foam::sqrt(min(great, max(0, a)));
+            (
+                ((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*denom),
+                sqrt(max(0, (d1 + d2)*(d2 + d3)*(d3 + d1)/(4*denom)))
+            );
     }
 }
 
@@ -169,15 +150,9 @@ inline Foam::scalar Foam::triangle<Point, PointRef>::circumRadius() const
 template<class Point, class PointRef>
 inline Foam::scalar Foam::triangle<Point, PointRef>::quality() const
 {
-    scalar c = circumRadius();
+    const scalar r = circumCircle().second();
-
+    static const scalar sqrt3 = sqrt(scalar(3));
-    if (c < rootVSmall)
+    return mag()/(0.75*sqrt3*sqr(min(r, great)) + rootVSmall);
-    {
-        // zero circumRadius, something has gone wrong.
-        return small;
-    }
-
-    return mag()/(Foam::sqr(c)*3.0*sqrt(3.0)/4.0);
 }