diff --git a/applications/test/triangleIntersection/Make/files b/applications/test/triangleIntersection/Make/files
new file mode 100644
index 0000000000..21bf7bf1b6
--- /dev/null
+++ b/applications/test/triangleIntersection/Make/files
@@ -0,0 +1,3 @@
+Test-triangleIntersection.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-triangleIntersection
diff --git a/applications/test/triangleIntersection/Make/options b/applications/test/triangleIntersection/Make/options
new file mode 100644
index 0000000000..54c035b8f5
--- /dev/null
+++ b/applications/test/triangleIntersection/Make/options
@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools
diff --git a/applications/test/triangleIntersection/Test-triangleIntersection.C b/applications/test/triangleIntersection/Test-triangleIntersection.C
new file mode 100644
index 0000000000..1aca11601c
--- /dev/null
+++ b/applications/test/triangleIntersection/Test-triangleIntersection.C
@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Test bounding box / triangle intersection
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "line.H"
+#include "Random.H"
+#include "triangle.H"
+#include "triangleFuncs.H"
+#include "treeBoundBox.H"
+#include "ListOps.H"
+
+using namespace Foam;
+
+//- simple helper to create a cube
+boundBox cube(scalar start, scalar width)
+{
+    return boundBox
+    (
+        point::uniform(start),
+        point::uniform(start + width)
+    );
+}
+
+
+void printEdges(const treeBoundBox& bb)
+{
+    pointField pts(bb.points());
+
+    for (const edge& e : treeBoundBox::edges)
+    {
+        Info<< pts[e.first()] << " -> " << pts[e.second()] << nl;
+    }
+}
+
+
+void testIntersect(const treeBoundBox& bb, const triPoints& tri)
+{
+    int ninside = 0;
+
+    if (bb.contains(tri.a())) ++ninside;
+    if (bb.contains(tri.b())) ++ninside;
+    if (bb.contains(tri.c())) ++ninside;
+
+    Info<< "box: " << bb << endl;
+    Info<< "tri: " << tri << endl;
+
+    Info<< "num inside: " << ninside << nl;
+    Info<< "intersects: " << bb.intersects(tri.tri())
+        << ' ' << triangleFuncs::intersectBb(tri.tri(), bb) << nl << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    argList args(argc, argv);
+
+    // Info<<"tree-bb faces: " << treeBoundBox::faces << nl
+    //     <<"tree-bb edges: " << treeBoundBox::edges << endl;
+
+    treeBoundBox bb;
+    bb = cube(0, 1);
+
+    triPoints tri;
+    tri.a() = point(-0.1, 0.5, 0.5);
+    tri.b() = point(0.1, 0.6, 0.5);
+    tri.c() = point(0.2, 0.4, 0.5);
+
+    testIntersect(bb, tri);
+
+    tri.a().z() = 1.1;
+    tri.b().z() = 1.1;
+    tri.c().z() = 1.1;
+
+    testIntersect(bb, tri);
+
+    tri.a().z() = 1 + 1e-15;
+    tri.b().z() = 1 + 1e-15;
+    tri.c().z() = 1 + 1e-15;
+
+    testIntersect(bb, tri);
+
+    tri.a() = point(-1, -1, -1);
+    tri.b() = point(2, 2, -2);
+    tri.c() = point(0, 0, 2);
+
+    testIntersect(bb, tri);
+
+    tri.a() = point(0.9, 1.1, 0);
+    tri.b() = point(1.1, 0.9, 0);
+    tri.c() = point(1, 1, 1.1);
+
+    testIntersect(bb, tri);
+
+
+    tri.a() = point(-1e-3, 1e-3,  -1);
+    tri.b() = point( 1e-3, -1e-3, -1);
+    tri.c() = point(0, 0, 2);
+
+    testIntersect(bb, tri);
+
+    tri.a() = point(-1e-3, 1e-3,  -1);
+    tri.b() = point( 1e-3, -1e-3, -1);
+    tri.c() = point(-1e-4, -1e-4, 2);
+
+    testIntersect(bb, tri);
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.C b/src/OpenFOAM/meshes/boundBox/boundBox.C
index 965a3b7f84..6bc66967bc 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBox.C
+++ b/src/OpenFOAM/meshes/boundBox/boundBox.C
@@ -231,6 +231,145 @@ bool Foam::boundBox::intersects(const plane& pln) const
 }
 
 
+bool Foam::boundBox::intersects(const triPointRef& tri) const
+{
+    // Require a full 3D box
+    if (nDim() != 3)
+    {
+        return false;
+    }
+
+    // Simplest check - if any points are inside
+    if (contains(tri.a()) || contains(tri.b()) || contains(tri.c()))
+    {
+        return true;
+    }
+
+
+    // Extent of box points projected onto axis
+    const auto project_box = []
+    (
+        const boundBox& bb,
+        const vector& axis,
+        scalarMinMax& extent
+    ) -> void
+    {
+        extent.reset(axis & bb.hexCorner<0>());
+        extent.add(axis & bb.hexCorner<1>());
+        extent.add(axis & bb.hexCorner<2>());
+        extent.add(axis & bb.hexCorner<3>());
+        extent.add(axis & bb.hexCorner<4>());
+        extent.add(axis & bb.hexCorner<5>());
+        extent.add(axis & bb.hexCorner<6>());
+        extent.add(axis & bb.hexCorner<7>());
+    };
+
+
+    // Use separating axis theorem to determine if triangle and
+    // (axis-aligned) bounding box intersect.
+
+    scalarMinMax tri_extent(0);
+    scalarMinMax box_extent(0);
+    const boundBox& bb = *this;
+
+    // 1.
+    // Test separating axis defined by the box normals
+    // (project triangle points)
+    // - do first (largely corresponds to normal bound box rejection test)
+    //
+    // No intersection if extent of projected triangle points are outside
+    // of the box range
+
+    {
+        // vector::X
+        tri_extent.reset(tri.a().x());
+        tri_extent.add(tri.b().x());
+        tri_extent.add(tri.c().x());
+
+        box_extent.reset(bb.min().x(), bb.max().x());
+
+        if (!tri_extent.overlaps(box_extent))
+        {
+            return false;
+        }
+
+        // vector::Y
+        tri_extent.reset(tri.a().y());
+        tri_extent.add(tri.b().y());
+        tri_extent.add(tri.c().y());
+
+        box_extent.reset(bb.min().y(), bb.max().y());
+
+        if (!tri_extent.overlaps(box_extent))
+        {
+            return false;
+        }
+
+        // vector::Z
+        tri_extent.reset(tri.a().z());
+        tri_extent.add(tri.b().z());
+        tri_extent.add(tri.c().z());
+
+        box_extent.reset(bb.min().z(), bb.max().z());
+
+        if (!tri_extent.overlaps(box_extent))
+        {
+            return false;
+        }
+    }
+
+
+    // 2.
+    // Test separating axis defined by the triangle normal
+    // (project box points)
+    // - can use area or unit normal since any scaling is applied to both
+    //   sides of the comparison.
+    // - by definition all triangle points lie in the plane defined by
+    //   the normal. It doesn't matter which of the points we use to define
+    //   the triangle offset (extent) when projected onto the triangle normal
+
+    vector axis = tri.areaNormal();
+
+    tri_extent.reset(axis & tri.a());
+    project_box(bb, axis, box_extent);
+
+    if (!tri_extent.overlaps(box_extent))
+    {
+        return false;
+    }
+
+
+    // 3.
+    // Test separating axes defined by the triangle edges, which are the
+    // cross product of the edge vectors and the box face normals
+
+    for (const vector& edgeVec : { tri.vecA(), tri.vecB(), tri.vecC() })
+    {
+        for (direction faceDir = 0; faceDir < vector::nComponents; ++faceDir)
+        {
+            axis = Zero;
+            axis[faceDir] = 1;
+
+            axis = (edgeVec ^ axis);
+
+            // project tri
+            tri_extent.reset(axis & tri.a());
+            tri_extent.add(axis & tri.b());
+            tri_extent.add(axis & tri.c());
+
+            project_box(bb, axis, box_extent);
+
+            if (!tri_extent.overlaps(box_extent))
+            {
+                return false;
+            }
+        }
+    }
+
+    return true;
+}
+
+
 bool Foam::boundBox::contains(const UList<point>& points) const
 {
     if (points.empty())
diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.H b/src/OpenFOAM/meshes/boundBox/boundBox.H
index f3eb5639c6..1373c5d754 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBox.H
+++ b/src/OpenFOAM/meshes/boundBox/boundBox.H
@@ -46,6 +46,7 @@ SeeAlso
 #include "pointField.H"
 #include "faceList.H"
 #include "Pair.H"
+#include "triangleFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -378,6 +379,12 @@ public:
             //  exactly with a box face
             bool intersects(const plane& pln) const;
 
+            //- Does triangle intersect this bounding box
+            //- or is contained within this bounding box.
+            //  \note results may be unreliable when it coincides almost
+            //  exactly with a box face
+            bool intersects(const triPointRef& tri) const;
+
             //- Overlaps/touches boundingBox?
             inline bool overlaps(const boundBox& bb) const;