diff --git a/applications/test/primitivePatch/Test-PrimitivePatch.C b/applications/test/primitivePatch/Test-PrimitivePatch.C
index 9a7e577676..b7c5ee617f 100644
--- a/applications/test/primitivePatch/Test-PrimitivePatch.C
+++ b/applications/test/primitivePatch/Test-PrimitivePatch.C
@@ -232,6 +232,7 @@ int main(int argc, char *argv[])
         const labelListList& edgeFaces = pp.edgeFaces();
         const labelListList& faceEdges = pp.faceEdges();
 
+        Pout<< "box: " << pp.box() << endl;
 
         checkFaceEdges(localFaces, edges, faceEdges);
 
diff --git a/src/OpenFOAM/meshes/meshShapes/cell/cell.C b/src/OpenFOAM/meshes/meshShapes/cell/cell.C
index 2ec1846256..a020f24e74 100644
--- a/src/OpenFOAM/meshes/meshShapes/cell/cell.C
+++ b/src/OpenFOAM/meshes/meshShapes/cell/cell.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +28,7 @@ License
 
 #include "cell.H"
 #include "pyramid.H"
+#include "primitiveMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -254,6 +255,36 @@ Foam::scalar Foam::cell::mag
 }
 
 
+Foam::Pair<Foam::point>
+Foam::cell::box
+(
+    const UList<point>& meshPoints,
+    const faceUList& meshFaces
+) const
+{
+    Pair<point> bb(point::rootMax, point::rootMin);
+
+    for (const label facei : *this)
+    {
+        for (const label pointi : meshFaces[facei])
+        {
+            const point& p = meshPoints[pointi];
+
+            bb.first()  = min(bb.first(), p);
+            bb.second() = max(bb.second(), p);
+        }
+    }
+
+    return bb;
+}
+
+
+Foam::Pair<Foam::point> Foam::cell::box(const primitiveMesh& mesh) const
+{
+    return cell::box(mesh.points(), mesh.faces());
+}
+
+
 // * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * * //
 
 bool Foam::operator==(const cell& a, const cell& b)
diff --git a/src/OpenFOAM/meshes/meshShapes/cell/cell.H b/src/OpenFOAM/meshes/meshShapes/cell/cell.H
index d9656d66ed..779aa4fecc 100644
--- a/src/OpenFOAM/meshes/meshShapes/cell/cell.H
+++ b/src/OpenFOAM/meshes/meshShapes/cell/cell.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,8 +37,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef cell_H
-#define cell_H
+#ifndef Foam_cell_H
+#define Foam_cell_H
 
 #include "faceList.H"
 #include "oppositeFace.H"
@@ -48,6 +48,9 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
+class primitiveMesh;
+
 /*---------------------------------------------------------------------------*\
                            Class cell Declaration
 \*---------------------------------------------------------------------------*/
@@ -146,6 +149,16 @@ public:
             const UList<point>& meshPoints,
             const faceUList& meshFaces
         ) const;
+
+        //- The bounding box for the cell
+        Pair<point> box
+        (
+            const UList<point>& meshPoints,
+            const faceUList& meshFaces
+        ) const;
+
+        //- The bounding box for the cell associated with given mesh
+        Pair<point> box(const primitiveMesh& mesh) const;
 };
 
 
diff --git a/src/OpenFOAM/meshes/meshShapes/edge/edge.H b/src/OpenFOAM/meshes/meshShapes/edge/edge.H
index 5d299d08bb..60feefd475 100644
--- a/src/OpenFOAM/meshes/meshShapes/edge/edge.H
+++ b/src/OpenFOAM/meshes/meshShapes/edge/edge.H
@@ -272,6 +272,9 @@ public:
         //  No special handling of negative point labels.
         inline scalar mag(const UList<point>& pts) const;
 
+        //- The enclosing (bounding) box for the edge
+        inline Pair<point> box(const UList<point>& pts) const;
+
         //- Return edge line
         //  No special handling of negative point labels.
         inline linePointRef line(const UList<point>& pts) const;
diff --git a/src/OpenFOAM/meshes/meshShapes/edge/edgeI.H b/src/OpenFOAM/meshes/meshShapes/edge/edgeI.H
index 31b7a65e0f..448bd43346 100644
--- a/src/OpenFOAM/meshes/meshShapes/edge/edgeI.H
+++ b/src/OpenFOAM/meshes/meshShapes/edge/edgeI.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -430,6 +430,22 @@ inline Foam::scalar Foam::edge::mag(const UList<point>& pts) const
 }
 
 
+inline Foam::Pair<Foam::point>
+Foam::edge::box(const UList<point>& pts) const
+{
+    #ifdef FULLDEBUG
+    if (first() < 0 || second() < 0)
+    {
+        FatalErrorInFunction
+            << "negative point index on edge " << *this
+            << abort(FatalError);
+    }
+    #endif
+
+    return linePointRef::box(pts[first()], pts[second()]);
+}
+
+
 inline Foam::linePointRef Foam::edge::line(const UList<point>& pts) const
 {
     #ifdef FULLDEBUG
diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.H b/src/OpenFOAM/meshes/meshShapes/face/face.H
index 1e64cd1881..b90c8d61e8 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.H
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -207,6 +207,9 @@ public:
         //- Magnitude of face area
         inline scalar mag(const UList<point>& p) const;
 
+        //- The enclosing (bounding) box for the face
+        inline Pair<point> box(const UList<point>& pts) const;
+
         //- Return face with reverse direction
         //  The starting points of the original and reverse face are identical.
         face reverseFace() const;
diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceI.H b/src/OpenFOAM/meshes/meshShapes/face/faceI.H
index f9379903c2..bc4ca27003 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/faceI.H
+++ b/src/OpenFOAM/meshes/meshShapes/face/faceI.H
@@ -115,6 +115,21 @@ inline Foam::scalar Foam::face::mag(const UList<point>& p) const
 }
 
 
+inline Foam::Pair<Foam::point>
+Foam::face::box(const UList<point>& pts) const
+{
+    Pair<point> bb(point::rootMax, point::rootMin);
+
+    for (const label pointi : *this)
+    {
+        bb.first()  = min(bb.first(),  pts[pointi]);
+        bb.second() = max(bb.second(), pts[pointi]);
+    }
+
+    return bb;
+}
+
+
 inline Foam::label Foam::face::nEdges() const noexcept
 {
     // for a closed polygon a number of edges is the same as number of points
diff --git a/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H b/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H
index 2a7489e5a0..957808e845 100644
--- a/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H
+++ b/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H
@@ -185,6 +185,9 @@ public:
         //- Magnitude of face area
         inline scalar mag(const UList<point>& points) const;
 
+        //- The enclosing (bounding) box for the face
+        inline Pair<point> box(const UList<point>& points) const;
+
         //- Number of triangles after splitting == 1
         inline label nTriangles() const noexcept;
 
diff --git a/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H b/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H
index 69ae3f889f..6a2d568609 100644
--- a/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H
+++ b/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H
@@ -217,6 +217,18 @@ inline Foam::scalar Foam::triFace::mag(const UList<point>& points) const
 }
 
 
+inline Foam::Pair<Foam::point>
+Foam::triFace::box(const UList<point>& points) const
+{
+    return triPointRef::box
+    (
+        points[operator[](0)],
+        points[operator[](1)],
+        points[operator[](2)]
+    );
+}
+
+
 inline Foam::label Foam::triFace::nTriangles() const noexcept
 {
     return 1;
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
index 4a5d635310..0ad2f0c0d9 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
@@ -461,6 +461,9 @@ public:
             //- Return point normals for patch
             const Field<point_type>& pointNormals() const;
 
+            //- The enclosing (bounding) box for the patch points
+            Pair<point_type> box() const;
+
 
         // Storage Management
 
@@ -482,6 +485,7 @@ public:
             bool hasPointEdges() const { return bool(pointEdgesPtr_); }
             bool hasPointFaces() const { return bool(pointFacesPtr_); }
 
+            bool hasMeshPoints() const { return bool(meshPointsPtr_); }
             bool hasMeshPointMap() const { return bool(meshPointMapPtr_); }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C
index f2fdbea1c3..2148b78fd9 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -122,4 +122,42 @@ Foam::PrimitivePatch<FaceList, PointField>::calcBdryPoints() const
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class FaceList, class PointField>
+Foam::Pair<typename Foam::PrimitivePatch<FaceList, PointField>::point_type>
+Foam::PrimitivePatch<FaceList, PointField>::box() const
+{
+    Pair<point_type> bb
+    (
+        point_type::uniform(ROOTVGREAT),
+        point_type::uniform(-ROOTVGREAT)
+    );
+
+    if (hasMeshPoints())
+    {
+        // Less looping if meshPoints() are already available
+        for (const label pointi : meshPoints())
+        {
+            bb.first()  = min(bb.first(),  points_[pointi]);
+            bb.second() = max(bb.second(), points_[pointi]);
+        }
+    }
+    else
+    {
+        // Walk the points on each face
+        for (const face_type& f : *this)
+        {
+            for (const label pointi : f)
+            {
+                bb.first()  = min(bb.first(),  points_[pointi]);
+                bb.second() = max(bb.second(), points_[pointi]);
+            }
+        }
+    }
+
+    return bb;
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveShapes/line/line.C b/src/OpenFOAM/meshes/primitiveShapes/line/line.C
index 826c9dd512..fcaad9fcaa 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/line/line.C
+++ b/src/OpenFOAM/meshes/primitiveShapes/line/line.C
@@ -114,7 +114,7 @@ scalar line<point2D, const point2D&>::nearestDist
         {
             // maxEdge inside interval of *this
             edgePt = maxEdgePt;
-            thisPt = nearestDist(edgePt).rawPoint();
+            thisPt = nearestDist(edgePt).point();
         }
         else
         {
@@ -124,13 +124,13 @@ scalar line<point2D, const point2D&>::nearestDist
                 // Edge completely envelops this. Take any this point and
                 // determine nearest on edge.
                 thisPt = minThisPt;
-                edgePt = e.nearestDist(thisPt).rawPoint();
+                edgePt = e.nearestDist(thisPt).point();
             }
             else if (minEdge < maxThis)
             {
                 // minEdge inside this interval.
                 edgePt = minEdgePt;
-                thisPt = nearestDist(edgePt).rawPoint();
+                thisPt = nearestDist(edgePt).point();
             }
             else
             {
diff --git a/src/OpenFOAM/meshes/primitiveShapes/line/line.H b/src/OpenFOAM/meshes/primitiveShapes/line/line.H
index e00dfcdfca..d1c0d436fd 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/line/line.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/line/line.H
@@ -43,8 +43,9 @@ SourceFiles
 #include "vector.H"
 #include "PointHit.H"
 #include "FixedList.H"
-#include "Pair.H"
 #include "UList.H"
+#include "Pair.H"
+#include "Tuple2.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -114,6 +115,27 @@ public:
 
         //- The second vertex
         point& b() { return Pair<point>::second(); }
+
+        //- Return as line reference
+        inline linePointRef ln() const;
+
+
+    // Properties
+
+        //- Return centre (centroid)
+        inline point centre() const;
+
+        //- The magnitude (length) of the line
+        inline scalar mag() const;
+
+        //- Return start-to-end vector
+        inline vector vec() const;
+
+        //- Return the unit vector (start-to-end)
+        inline vector unitVec() const;
+
+        //- The enclosing (bounding) box for the line
+        inline Pair<point> box() const;
 };
 
 
@@ -132,7 +154,6 @@ class line
         //- Second point
         PointRef b_;
 
-
 public:
 
     // Constructors
@@ -178,12 +199,18 @@ public:
         PointRef end() const noexcept { return b_; }
 
 
+    // Line properties (static calculations)
+
+        //- The enclosing (bounding) box for two points
+        inline static Pair<Point> box(const Point& p0, const Point& p1);
+
+
     // Properties
 
         //- Return centre (centroid)
         inline Point centre() const;
 
-        //- Return scalar magnitude
+        //- The magnitude (length) of the line
         inline scalar mag() const;
 
         //- Return start-to-end vector
@@ -192,6 +219,9 @@ public:
         //- Return the unit vector (start-to-end)
         inline Point unitVec() const;
 
+        //- The enclosing (bounding) box for the line
+        inline Pair<Point> box() const;
+
         //- Return nearest distance to line from a given point
         //  If the nearest point is on the line, return a hit
         PointHit<Point> nearestDist(const Point& p) const;
diff --git a/src/OpenFOAM/meshes/primitiveShapes/line/lineI.H b/src/OpenFOAM/meshes/primitiveShapes/line/lineI.H
index 217d0dee31..d4bb63aada 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/line/lineI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/line/lineI.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,7 +48,11 @@ inline Foam::linePoints::linePoints
 
 
 template<class Point, class PointRef>
-inline Foam::line<Point, PointRef>::line(const Point& from, const Point& to)
+inline Foam::line<Point, PointRef>::line
+(
+    const Point& from,
+    const Point& to
+)
 :
     a_(from),
     b_(to)
@@ -76,6 +80,12 @@ inline Foam::line<Point, PointRef>::line(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline Foam::linePointRef Foam::linePoints::ln() const
+{
+    return linePointRef(a(), b());
+}
+
+
 template<class Point, class PointRef>
 inline Point Foam::line<Point, PointRef>::centre() const
 {
@@ -83,17 +93,35 @@ inline Point Foam::line<Point, PointRef>::centre() const
 }
 
 
+inline Foam::point Foam::linePoints::centre() const
+{
+    return 0.5*(a() + b());
+}
+
+
 template<class Point, class PointRef>
 inline Foam::scalar Foam::line<Point, PointRef>::mag() const
 {
-    return ::Foam::mag(vec());
+    return ::Foam::mag(b() - a());
+}
+
+
+inline Foam::scalar Foam::linePoints::mag() const
+{
+    return ::Foam::mag(b() - a());
 }
 
 
 template<class Point, class PointRef>
 inline Point Foam::line<Point, PointRef>::vec() const
 {
-    return b_ - a_;
+    return (b() - a());
+}
+
+
+inline Foam::vector Foam::linePoints::vec() const
+{
+    return (b() - a());
 }
 
 
@@ -107,6 +135,36 @@ inline Point Foam::line<Point, PointRef>::unitVec() const
 }
 
 
+inline Foam::vector Foam::linePoints::unitVec() const
+{
+    return normalised(b() - a());
+}
+
+
+template<class Point, class PointRef>
+inline Foam::Pair<Point> Foam::line<Point, PointRef>::box
+(
+    const Point& p0,
+    const Point& p1
+)
+{
+    return Pair<Point>(min(p0, p1), max(p0, p1));
+}
+
+
+template<class Point, class PointRef>
+inline Foam::Pair<Point> Foam::line<Point, PointRef>::box() const
+{
+    return line<Point, PointRef>::box(a_, b_);
+}
+
+
+inline Foam::Pair<Foam::point> Foam::linePoints::box() const
+{
+    return linePointRef::box(a(), b());
+}
+
+
 template<class Point, class PointRef>
 Foam::PointHit<Point> Foam::line<Point, PointRef>::nearestDist
 (
@@ -186,12 +244,12 @@ Foam::scalar Foam::line<Point, PointRef>::nearestDist
             {
                 if (this0.distance() < this1.distance())
                 {
-                    thisPt = this0.rawPoint();
+                    thisPt = this0.point();
                     edgePt = edge.start();
                 }
                 else
                 {
-                    thisPt = this1.rawPoint();
+                    thisPt = this1.point();
                     edgePt = edge.end();
                 }
             }
@@ -200,12 +258,12 @@ Foam::scalar Foam::line<Point, PointRef>::nearestDist
                 if (edge0.distance() < edge1.distance())
                 {
                     thisPt = start();
-                    edgePt = edge0.rawPoint();
+                    edgePt = edge0.point();
                 }
                 else
                 {
                     thisPt = end();
-                    edgePt = edge1.rawPoint();
+                    edgePt = edge1.point();
                 }
             }
         }
@@ -243,7 +301,7 @@ Foam::scalar Foam::line<Point, PointRef>::nearestDist
         {
             // maxEdge inside interval of *this
             edgePt = maxEdgePt;
-            thisPt = nearestDist(edgePt).rawPoint();
+            thisPt = nearestDist(edgePt).point();
         }
         else
         {
@@ -253,13 +311,13 @@ Foam::scalar Foam::line<Point, PointRef>::nearestDist
                 // Edge completely envelops this. Take any this point and
                 // determine nearest on edge.
                 thisPt = minThisPt;
-                edgePt = edge.nearestDist(thisPt).rawPoint();
+                edgePt = edge.nearestDist(thisPt).point();
             }
             else if (minEdge < maxThis)
             {
                 // minEdge inside this interval.
                 edgePt = minEdgePt;
-                thisPt = nearestDist(edgePt).rawPoint();
+                thisPt = nearestDist(edgePt).point();
             }
             else
             {
diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
index 4373becba3..f158b5a2d3 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
@@ -96,40 +96,41 @@ class plane
 {
 public:
 
-    //- Side of the plane
-    enum side
-    {
-        FRONT = 1,  //!< The front (positive normal) side of the plane
-        BACK  = -1, //!< The back (negative normal) side of the plane
-        NORMAL = 1, //!< Same as FRONT
-        FLIP  = -1  //!< Same as BACK
-    };
+    // Public Data Types
 
-
-    //- A reference point and direction
-    class ray
-    {
-        point  pt_;
-        vector dir_;
-
-    public:
-
-        ray(const point& pt, const vector& dir)
-        :
-            pt_(pt),
-            dir_(dir)
-        {}
-
-        const point& refPoint() const noexcept
+        //- Side of the plane
+        enum side
         {
-            return pt_;
-        }
+            FRONT = 1,  //!< The front (positive normal) side of the plane
+            BACK  = -1, //!< The back (negative normal) side of the plane
+            NORMAL = FRONT,   //!< Alias for FRONT
+            FLIP   = BACK     //!< Alias for BACK
+        };
 
-        const vector& dir() const noexcept
+
+    // Public Classes
+
+        //- A reference point and direction
+        class ray
         {
-            return dir_;
-        }
-    };
+            point  point_;
+            vector dir_;
+
+        public:
+
+            //- Construct from reference point and direction
+            ray(const point& p, const vector& dir)
+            :
+                point_(p),
+                dir_(dir)
+            {}
+
+            //- Return the reference point
+            const point& refPoint() const noexcept { return point_; }
+
+            //- Return the direction
+            const vector& dir() const noexcept { return dir_; }
+        };
 
 
 private:
@@ -255,9 +256,9 @@ public:
         }
 
         //- Return the cutting point between the plane and
-        // a line passing through the supplied points
-        template<class Point, class PointRef>
-        scalar lineIntersect(const line<Point, PointRef>& l) const
+        //- a line passing through the supplied points
+        template<class PointType, class PointRef>
+        scalar lineIntersect(const line<PointType, PointRef>& l) const
         {
             return normalIntersect(l.start(), l.vec());
         }
@@ -278,15 +279,18 @@ public:
         point somePointInPlane(const scalar dist = 1e-3) const;
 
         //- Return the side of the plane that the point is on.
-        //  If the point is on the plane, then returns NORMAL.
-        inline side sideOfPlane(const point& p) const;
+        //  If the point is on the plane, then returns FRONT.
+        //  \return
+        //  - +1 (FRONT): above or on plane (front-side)
+        //  - -1 (BACK):  below plane (back-side)
+        inline side whichSide(const point& p) const;
 
         //- The sign for the side of the plane that the point is on.
         //  Uses the supplied tolerance for rounding around zero.
         //  \return
         //  -  0: on plane
-        //  - +1: front-side
-        //  - -1: back-side
+        //  - +1 (FRONT): above plane (front-side)
+        //  - -1 (BACK):  below plane (back-side)
         inline int sign(const point& p, const scalar tol = SMALL) const;
 
         //- Mirror the supplied point in the plane. Return the mirrored point.
@@ -300,6 +304,9 @@ public:
 
         //- The plane base point (same as origin)
         const point& refPoint() const noexcept { return origin_; }
+
+        //- Same as whichSide()
+        side sideOfPlane(const point& p) const { return whichSide(p); }
 };
 
 
diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/planeI.H b/src/OpenFOAM/meshes/primitiveShapes/plane/planeI.H
index 695f191afb..f47c2c7cdc 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/plane/planeI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/plane/planeI.H
@@ -78,11 +78,11 @@ inline Foam::scalar Foam::plane::signedDistance(const point& p) const
 }
 
 
-inline Foam::plane::side Foam::plane::sideOfPlane(const point& p) const
+inline Foam::plane::side Foam::plane::whichSide(const point& p) const
 {
     const scalar dist = signedDistance(p);
 
-    return (dist < 0 ? BACK : FRONT);
+    return (dist < 0 ? side::BACK : side::FRONT);
 }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C
index 0db8ddd168..1e28580e79 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C
+++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C
@@ -53,7 +53,7 @@ Foam::pointHit Foam::tetrahedron<Point, PointRef>::containmentSphere
 
     pointHit pHit(circumCentre());
     pHit.setHit();
-    scalar minRadiusSqr = magSqr(pHit.rawPoint() - a_);
+    scalar minRadiusSqr = magSqr(pHit.point() - a_);
 
 
     // 2. Try circumcentre of tet triangles. Create circumcircle for triFace and
@@ -337,16 +337,4 @@ void Foam::tetrahedron<Point, PointRef>::gradNiGradNj
 }
 
 
-template<class Point, class PointRef>
-Foam::boundBox Foam::tetrahedron<Point, PointRef>::bounds() const
-{
-    return
-        boundBox
-        (
-            min(a(), min(b(), min(c(), d()))),
-            max(a(), max(b(), max(c(), d())))
-        );
-}
-
-
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
index ae8323f243..c578722a53 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
@@ -149,6 +149,9 @@ public:
         //- Invert tetrahedron by swapping third and fourth vertices
         inline void flip();
 
+        //- The enclosing (bounding) box for the tetrahedron
+        inline Pair<point> box() const;
+
         //- The bounding box for the tetrahedron
         inline treeBoundBox bounds() const;
 };
@@ -296,6 +299,18 @@ public:
         inline triPointRef tri(const label facei) const;
 
 
+    // Tet properties (static calculations)
+
+        //- The enclosing (bounding) box for four points
+        inline static Pair<Point> box
+        (
+            const Point& p0,
+            const Point& p1,
+            const Point& p2,
+            const Point& p3
+        );
+
+
         // Properties
 
             //- Face area normal for side a
@@ -316,6 +331,12 @@ public:
             //- Return volume
             inline scalar mag() const;
 
+            //- The enclosing (bounding) box for the tetrahedron
+            inline Pair<Point> box() const;
+
+            //- The bounding box for the tetrahedron
+            inline treeBoundBox bounds() const;
+
             //- Return circum-centre
             inline Point circumCentre() const;
 
@@ -392,9 +413,6 @@ public:
 
             void gradNiGradNj(tensorField& buffer) const;
 
-            //- Calculate the bounding box
-            boundBox bounds() const;
-
 
     // IOstream Operators
 
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
index 4dc069a134..196a60d459 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
@@ -139,11 +139,50 @@ inline void Foam::tetPoints::flip()
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Point, class PointRef>
+inline Foam::Pair<Point> Foam::tetrahedron<Point, PointRef>::box
+(
+    const Point& p0,
+    const Point& p1,
+    const Point& p2,
+    const Point& p3
+)
+{
+    return Pair<Point>
+    (
+        min(p0, min(p1, min(p2, p3))),
+        max(p0, max(p1, max(p2, p3)))
+    );
+}
+
+
+template<class Point, class PointRef>
+inline Foam::Pair<Point> Foam::tetrahedron<Point, PointRef>::box() const
+{
+    return tetrahedron<Point, PointRef>::box(a_, b_, c_, d_);
+}
+
+
+inline Foam::Pair<Foam::point> Foam::tetPoints::box() const
+{
+    return tetPointRef::box(a(), b(), c(), d());
+}
+
+
+template<class Point, class PointRef>
+inline Foam::treeBoundBox Foam::tetrahedron<Point, PointRef>::bounds() const
+{
+    Pair<point> bb(tetrahedron<Point, PointRef>::box());
+    return treeBoundBox(bb.first(), bb.second());
+}
+
+
 inline Foam::treeBoundBox Foam::tetPoints::bounds() const
 {
-    treeBoundBox bb;
-    bb.add(static_cast<const FixedList<point, 4>&>(*this));
-    return bb;
+    Pair<point> bb(tetPoints::box());
+    return treeBoundBox(bb.first(), bb.second());
 }
 
 
@@ -285,7 +324,7 @@ inline Foam::scalar Foam::tetrahedron<Point, PointRef>::quality() const
         mag()
        /(
             8.0/(9.0*sqrt(3.0))
-           *pow3(min(circumRadius(), GREAT))
+           *pow3(Foam::min(circumRadius(), GREAT))
           + ROOTVSMALL
         );
 }
@@ -397,7 +436,7 @@ inline Foam::pointHit Foam::tetrahedron<Point, PointRef>::nearestPoint
 
             if (info.distance() < minOutsideDistance)
             {
-                closestPt = info.rawPoint();
+                closestPt = info.point();
 
                 minOutsideDistance = info.distance();
             }
@@ -417,7 +456,7 @@ inline Foam::pointHit Foam::tetrahedron<Point, PointRef>::nearestPoint
 
             if (info.distance() < minOutsideDistance)
             {
-                closestPt = info.rawPoint();
+                closestPt = info.point();
 
                 minOutsideDistance = info.distance();
             }
@@ -437,7 +476,7 @@ inline Foam::pointHit Foam::tetrahedron<Point, PointRef>::nearestPoint
 
             if (info.distance() < minOutsideDistance)
             {
-                closestPt = info.rawPoint();
+                closestPt = info.point();
 
                 minOutsideDistance = info.distance();
             }
@@ -457,7 +496,7 @@ inline Foam::pointHit Foam::tetrahedron<Point, PointRef>::nearestPoint
 
             if (info.distance() < minOutsideDistance)
             {
-                closestPt = info.rawPoint();
+                closestPt = info.point();
 
                 minOutsideDistance = info.distance();
             }
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
index 83d61172e4..5bdd6979be 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
@@ -29,9 +29,11 @@ Class
 
 Description
     A triangle primitive used to calculate face normals and swept volumes.
+    Uses referred points.
 
 SourceFiles
     triangleI.H
+    triangle.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -47,6 +49,8 @@ SourceFiles
 #include "FixedList.H"
 #include "UList.H"
 #include "line.H"
+#include "Pair.H"
+#include "Tuple2.H"
 #include "barycentric2D.H"
 #include "treeBoundBox.H"
 
@@ -150,8 +154,8 @@ public:
         //- The magnitude of the triangle area
         inline scalar mag() const;
 
-        //- The bounding box for the tetrahedron
-        inline treeBoundBox bounds() const;
+        //- The enclosing (bounding) box for the triangle
+        inline Pair<point> box() const;
 };
 
 
@@ -311,6 +315,14 @@ public:
             const Point& p2
         );
 
+        //- The enclosing (bounding) box for three points
+        inline static Pair<Point> box
+        (
+            const Point& p0,
+            const Point& p1,
+            const Point& p2
+        );
+
 
         // Properties
 
@@ -334,6 +346,9 @@ public:
             //- The magnitude of the triangle area
             inline scalar mag() const;
 
+            //- The enclosing (bounding) box for the triangle
+            inline Pair<Point> box() const;
+
             //- Return circum-centre
             inline Point circumCentre() const;
 
@@ -374,6 +389,22 @@ public:
                 barycentric2D& bary
             ) const;
 
+            //- Fast intersection detection with a plane.
+            inline bool intersects
+            (
+                const point& origin,
+                const vector& normal
+            ) const;
+
+            //- Fast intersection detection with an \b axis plane.
+            inline bool intersects
+            (
+                //! Origin of the plane
+                const point& origin,
+                //! Normal of the plane (vector::X, vector::Y, vector::Z)
+                const vector::components axis
+            ) const;
+
             //- Return point intersection with a ray.
             //  For a hit, the distance is signed. Positive number
             //  represents the point in front of triangle.
@@ -439,8 +470,8 @@ public:
             //  Uses the supplied tolerance for rounding around zero.
             //  \return
             //  -  0: on plane
-            //  - +1: front-side
-            //  - -1: back-side
+            //  - +1: above plane
+            //  - -1: below plane
             inline int sign(const point& p, const scalar tol = SMALL) const;
 
             //- Decompose triangle into triangles above and below plane
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
index 3be9f01eba..45adc5fbc2 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
@@ -184,7 +184,7 @@ inline Foam::vector Foam::triangle<Point, PointRef>::unitNormal
 template<class Point, class PointRef>
 inline Foam::vector Foam::triangle<Point, PointRef>::unitNormal() const
 {
-    return triPointRef::unitNormal(a_, b_, c_);
+    return triangle<Point, PointRef>::unitNormal(a_, b_, c_);
 }
 
 
@@ -194,6 +194,35 @@ inline Foam::vector Foam::triPoints::unitNormal() const
 }
 
 
+template<class Point, class PointRef>
+inline Foam::Pair<Point> Foam::triangle<Point, PointRef>::box
+(
+    const Point& p0,
+    const Point& p1,
+    const Point& p2
+)
+{
+    return Pair<Point>
+    (
+        min(p0, min(p1, p2)),
+        max(p0, max(p1, p2))
+    );
+}
+
+
+template<class Point, class PointRef>
+inline Foam::Pair<Point> Foam::triangle<Point, PointRef>::box() const
+{
+    return triangle<Point, PointRef>::box(a_, b_, c_);
+}
+
+
+inline Foam::Pair<Foam::point> Foam::triPoints::box() const
+{
+    return triPointRef::box(a(), b(), c());
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline Foam::triPointRef Foam::triPoints::tri() const
@@ -211,14 +240,6 @@ inline void Foam::triPoints::flip()
 }
 
 
-inline Foam::treeBoundBox Foam::triPoints::bounds() const
-{
-    treeBoundBox bb;
-    bb.add(static_cast<const FixedList<point, 3>&>(*this));
-    return bb;
-}
-
-
 inline Foam::scalar Foam::triPoints::mag() const
 {
     return ::Foam::mag(areaNormal());
@@ -419,6 +440,46 @@ inline Foam::scalar Foam::triangle<Point, PointRef>::pointToBarycentric
 }
 
 
+template<class Point, class PointRef>
+inline bool Foam::triangle<Point, PointRef>::intersects
+(
+    const point& origin,
+    const vector& normal
+) const
+{
+    // Check points cut below(1) or above(2). Mix of above/below == 3
+    // Cf. plane::whichSide()
+    return
+    (
+        (
+            (((a_ - origin) & normal) < 0 ? 1 : 2)
+          | (((b_ - origin) & normal) < 0 ? 1 : 2)
+          | (((c_ - origin) & normal) < 0 ? 1 : 2)
+        ) == 3
+    );
+}
+
+
+template<class Point, class PointRef>
+inline bool Foam::triangle<Point, PointRef>::intersects
+(
+    const point& origin,
+    const vector::components axis
+) const
+{
+    // Direct check of points cut below(1) or above(2)
+    // Cf. plane::whichSide()
+    return
+    (
+        (
+            (a_[axis] < origin[axis] ? 1 : 2)
+          | (b_[axis] < origin[axis] ? 1 : 2)
+          | (c_[axis] < origin[axis] ? 1 : 2)
+        ) == 3
+    );
+}
+
+
 template<class Point, class PointRef>
 inline Foam::pointHit Foam::triangle<Point, PointRef>::ray
 (
@@ -476,7 +537,7 @@ inline Foam::pointHit Foam::triangle<Point, PointRef>::ray
 
         if (hit)
         {
-            pInter = fastInter.rawPoint();
+            pInter = fastInter.point();
         }
         else
         {
@@ -521,7 +582,7 @@ inline Foam::pointHit Foam::triangle<Point, PointRef>::ray
         inter.setMiss(eligible);
 
         // The miss point is the nearest point on the triangle
-        inter.setPoint(nearestPoint(p).rawPoint());
+        inter.setPoint(nearestPoint(p).point());
 
         // The distance to the miss is the distance between the
         // original point and plane of intersection
diff --git a/src/OpenFOAM/primitives/ranges/MinMax/MinMax.H b/src/OpenFOAM/primitives/ranges/MinMax/MinMax.H
index 7792a242dd..ee2e6c479c 100644
--- a/src/OpenFOAM/primitives/ranges/MinMax/MinMax.H
+++ b/src/OpenFOAM/primitives/ranges/MinMax/MinMax.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -95,8 +95,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef MinMax_H
-#define MinMax_H
+#ifndef Foam_MinMax_H
+#define Foam_MinMax_H
 
 #include "scalar.H"
 #include "Pair.H"
@@ -109,7 +109,7 @@ Description
 namespace Foam
 {
 
-// Forward declarations
+// Forward Declarations
 template<class T> class MinMax;
 class zero;
 
@@ -205,8 +205,11 @@ public:
         //- Range is valid if it is not inverted
         inline bool valid() const;
 
-        //- Reset to an invalid, inverted range
-        inline void clear();
+        //- Reset to an inverted (invalid) range
+        inline void reset();
+
+        //- Same as reset() - reset to an inverted (invalid) range
+        void clear() { reset(); }
 
 
     // Testing / Query
diff --git a/src/OpenFOAM/primitives/ranges/MinMax/MinMaxI.H b/src/OpenFOAM/primitives/ranges/MinMax/MinMaxI.H
index 02d167f050..c656b0328e 100644
--- a/src/OpenFOAM/primitives/ranges/MinMax/MinMaxI.H
+++ b/src/OpenFOAM/primitives/ranges/MinMax/MinMaxI.H
@@ -168,7 +168,7 @@ inline bool Foam::MinMax<T>::valid() const
 
 
 template<class T>
-inline void Foam::MinMax<T>::clear()
+inline void Foam::MinMax<T>::reset()
 {
     min() = pTraits<T>::max;
     max() = pTraits<T>::min;