ENH: face, triFace sign() method analogous to plane

- a quick test for which side of the face a point is located

applications/test/faces/Test-faces.C

@@ -31,19 +31,81 @@ Description
 
 #include "argList.H"
 #include "labelledTri.H"
+#include "pointList.H"
+#include "ListOps.H"
 
 using namespace Foam;
 
+
+template<class Face>
+void faceInfo(const Face& f, const UList<point>& points)
+{
+    Info<< f
+        << " points:" << f.points(points)
+        << " normal:" << f.unitNormal(points);
+}
+
+
+template<class Face>
+void testSign
+(
+    const Face& f,
+    const UList<point>& points,
+    const UList<point>& testPoints
+)
+{
+    for (const point& p : testPoints)
+    {
+        Info<< "  point:" << p << " sign=" << f.sign(p, points) << nl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 //  Main program:
 
 int main(int argc, char *argv[])
 {
+    pointList points1
+    ({
+        { 0, 0, 0},
+        { -1, -1, 1},
+        { 1, -1, -1},
+        { 1, 1, -1},
+        { -1, 1, 1}
+    });
+
+    pointList points2 = ListOps::create<point>
+    (
+        points1,
+        [](const point& p){ return point(p.x(), p.y(), -p.z()); }
+    );
+
+    pointList testPoints
+    ({
+        { -2, -2, -2},
+        { -2, -2, 2},
+        { 0, 0, 0},
+        { 2, 2, -2},
+        { 2, 2, 2}
+    });
+
     face f1{ 1, 2, 3, 4 };
-    Info<< "face:" << f1 << nl;
+    Info<< "face:"; faceInfo(f1, points1); Info << nl;
+    testSign(f1, points1, testPoints);
+
+    Info<< "face:"; faceInfo(f1, points2); Info << nl;
+    testSign(f1, points2, testPoints);
+    Info<< nl;
 
     triFace t1{ 1, 2, 3 };
-    Info<< "triFace:" << t1 << nl;
+    Info<< "triFace:"; faceInfo(t1, points1); Info << nl;
+    testSign(t1, points1, testPoints);
+
+    Info<< "triFace:"; faceInfo(t1, points2); Info << nl;
+    testSign(t1, points2, testPoints);
+    Info<< nl;
+
 
     f1 = t1;
     Info<< "face:" << f1 << nl;

src/OpenFOAM/meshes/meshShapes/face/face.H

@@ -299,6 +299,20 @@ public:
             label& nearLabel
         ) const;
 
+        //- The sign for the side of the face plane the point is on,
+        //- using three evenly distributed face points for the estimated normal.
+        //  Uses the supplied tolerance for rounding around zero.
+        //  \return
+        //  -  0: on plane
+        //  - +1: front-side
+        //  - -1: back-side
+        int sign
+        (
+            const point& p,
+            const UList<point>& points,
+            const scalar tol = SMALL
+        ) const;
+
         //- Return contact sphere diameter
         scalar contactSphereDiameter
         (

src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C

@@ -309,4 +309,23 @@ Foam::pointHit Foam::face::nearestPointClassify
 }
 
 
+int Foam::face::sign
+(
+    const point& p,
+    const UList<point>& points,
+    const scalar tol
+) const
+{
+    // Take three points [0, 1/3, 2/3] from the face
+    // - assumes the face is not severely warped
+
+    return triPointRef
+    (
+        points[operator[](0)],
+        points[operator[](size()/3)],
+        points[operator[]((2*size())/3)]
+    ).sign(p, tol);
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/meshShapes/triFace/triFace.H

@@ -221,6 +221,19 @@ public:
             label& nearLabel
         ) const;
 
+        //- The sign for which side of the face plane the point is on.
+        //  Uses the supplied tolerance for rounding around zero.
+        //  \return
+        //  -  0: on plane
+        //  - +1: front-side
+        //  - -1: back-side
+        inline int sign
+        (
+            const point& p,
+            const UList<point>& points,
+            const scalar tol = SMALL
+        ) const;
+
         //- Return number of edges
         inline label nEdges() const;
 

src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H

@@ -332,6 +332,17 @@ inline Foam::pointHit Foam::triFace::nearestPointClassify
 }
 
 
+inline int Foam::triFace::sign
+(
+    const point& p,
+    const UList<point>& points,
+    const scalar tol
+) const
+{
+    return this->tri(points).sign(p, tol);
+}
+
+
 inline Foam::label Foam::triFace::nEdges() const
 {
     return 3;

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

@@ -323,6 +323,14 @@ public:
                 pointHit& edgePoint
             ) const;
 
+            //- The sign for which side of the face plane the point is on.
+            //  Uses the supplied tolerance for rounding around zero.
+            //  \return
+            //  -  0: on plane
+            //  - +1: front-side
+            //  - -1: back-side
+            inline int sign(const point& p, const scalar tol = SMALL) const;
+
             //- Decompose triangle into triangles above and below plane
             template<class AboveOp, class BelowOp>
             inline void sliceWithPlane

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

@@ -815,6 +815,19 @@ inline Foam::pointHit Foam::triangle<Point, PointRef>::nearestPoint
 }
 
 
+template<class Point, class PointRef>
+inline int Foam::triangle<Point, PointRef>::sign
+(
+    const point& p,
+    const scalar tol
+) const
+{
+    const scalar dist = ((p - a_) & unitNormal());
+
+    return ((dist < -tol) ? -1 : (dist > tol) ? +1 : 0);
+}
+
+
 template<class Point, class PointRef>
 inline void Foam::triangle<Point, PointRef>::dummyOp::operator()
 (