From e5006a62d73e6ded3158a3a3a24dd47d7dde471a Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 1 Nov 2022 12:15:08 +0100
Subject: [PATCH] ENH: use simpler boundBox handling

- use default initialize boundBox instead of invertedBox
- reset() instead of assigning from invertedBox
- extend (three parameter version) and grow method
- inflate(Random) instead of extend + re-assigning
---
 .../findCell-octree/Test-findCell-octree.C    |  9 ++---
 .../extrude/extrudeMesh/extrudeMesh.C         |  3 +-
 .../backgroundMeshDecomposition.C             |  4 +-
 .../conformalVoronoiMeshConformToSurface.C    | 10 +----
 .../conformationSurfaces.C                    |  7 +---
 .../foamyHexMeshSurfaceSimplify_non_octree.C  |  6 +--
 .../manipulation/checkMesh/checkTopology.C    | 15 ++++----
 .../viewFactorsGen/searchingEngine.H          |  5 +--
 .../viewFactorsGen/searchingEngine_CGAL.H     |  5 +--
 .../surface/surfaceHookUp/surfaceHookUp.C     |  5 +--
 .../surfaceRedistributePar.C                  |  7 +---
 src/OpenFOAM/algorithms/AABBTree/AABBTree.C   |  4 +-
 src/OpenFOAM/meshes/polyMesh/polyMesh.C       |  7 +---
 .../PatchTools/PatchToolsSearch.C             |  2 +-
 .../polyMeshGeometry/polyMeshGeometry.C       |  2 +-
 .../faMesh/faMeshDemandDrivenData.C           | 22 +++--------
 .../turbulentDFSEMInletFvPatchVectorField.C   |  4 +-
 .../fieldExtents/fieldExtentsTemplates.C      |  6 +--
 .../field/streamLine/streamLineBase.C         |  2 +-
 .../energySpectrum/energySpectrum.C           | 10 ++---
 ...surfaceDisplacementPointPatchVectorField.C |  5 +--
 ...aceSlipDisplacementPointPatchVectorField.C |  7 ++--
 .../basic/InteractionLists/InteractionLists.C | 38 +++++--------------
 src/mesh/blockMesh/PDRblockMesh/PDRblock.C    |  2 +-
 .../projectVertex/projectVertex.C             |  2 +-
 .../meshRefinement/meshRefinement.C           |  9 +++--
 .../refinementFeatures/refinementFeatures.C   |  8 +---
 .../shellSurfaces/shellSurfaces.C             |  2 +-
 .../snappySnapDriverFeature.C                 |  4 +-
 .../advancingFrontAMI/advancingFrontAMI.C     |  2 +-
 .../faceAreaWeightAMI2D/faceAreaWeightAMI2D.C |  6 +--
 .../extendedEdgeMesh/extendedEdgeMesh.C       | 15 ++------
 .../mappedPolyPatch/mappedPatchBase.C         | 22 +++++------
 src/meshTools/meshSearch/meshSearch.C         | 38 +++++++------------
 src/meshTools/processorLOD/cellBox/cellBox.C  | 14 +++----
 .../searchableExtrudedCircle.C                | 25 +++++-------
 .../searchableSurfaceCollection.C             |  4 +-
 .../searchableSurfacesQueries.C               |  2 +-
 .../triSurfaceMesh/triSurfaceMesh.C           | 27 ++++---------
 .../triSurfaceSearch/triSurfaceRegionSearch.C |  6 +--
 .../triSurfaceSearch/triSurfaceSearch.C       |  8 ++--
 .../inverseDistanceCellCellStencil.C          | 17 +++------
 .../trackingInverseDistanceCellCellStencil.C  | 18 ++++-----
 .../distributedTriSurfaceMesh.C               | 19 +++++-----
 .../exact/exactPatchDistMethod.C              |  4 +-
 .../meshToMeshMethod/meshToMeshMethod.C       |  8 +---
 .../calculateMeshToMesh0Addressing.C          | 25 +++++-------
 src/sampling/probes/patchProbes.C             |  6 +--
 .../sampledSet/cellCentre/cellCentreSet.C     |  2 +-
 .../sampledSet/cellCentre/cellCentreSet.H     |  2 +-
 .../sampledSet/patchCloud/patchCloudSet.C     |  7 +---
 .../sampledSet/patchSeed/patchSeedSet.C       | 21 ++--------
 .../sampledCuttingSurface.C                   |  2 +-
 .../sampledCuttingSurface.H                   |  2 +-
 .../sampledPlane/sampledPlane.C               |  2 +-
 .../surface/isoSurface/isoSurfaceParams.C     |  2 +-
 .../surface/isoSurface/isoSurfacePoint.C      |  2 +-
 src/surfMesh/triSurface/triSurfaceIO.C        |  2 +-
 .../solarLoad/faceReflecting/faceReflecting.C | 10 ++---
 .../solarLoad/faceShading/faceShading.C       | 10 ++---
 60 files changed, 189 insertions(+), 353 deletions(-)

diff --git a/applications/test/findCell-octree/Test-findCell-octree.C b/applications/test/findCell-octree/Test-findCell-octree.C
index aa8db650aa..2210753eda 100644
--- a/applications/test/findCell-octree/Test-findCell-octree.C
+++ b/applications/test/findCell-octree/Test-findCell-octree.C
@@ -53,20 +53,17 @@ int main(int argc, char *argv[])
     const polyMesh::cellDecomposition decompMode = polyMesh::CELL_TETS;
 
     treeBoundBox meshBb(mesh.bounds());
+    treeBoundBox shiftedBb(meshBb);
 
     // Calculate typical cell related size to shift bb by.
     scalar typDim = meshBb.avgDim()/(2.0*Foam::cbrt(scalar(mesh.nCells())));
 
-    treeBoundBox shiftedBb
-    (
-        meshBb.min(),
-        meshBb.max() + vector(typDim, typDim, typDim)
-    );
+    shiftedBb.max() += vector::uniform(typDim);
 
     Info<< "Mesh" << endl;
     Info<< "   bounding box           : " << meshBb << endl;
     Info<< "   bounding box (shifted) : " << shiftedBb << endl;
-    Info<< "   typical dimension      : " << shiftedBb.typDim() << endl;
+    Info<< "   typical dimension      : " << shiftedBb.avgDim() << endl;
 
     Info<< "Initialised mesh in "
         << runTime.cpuTimeIncrement() << " s" << endl;
diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
index cd84e849cc..7ea12e1f2a 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
@@ -890,11 +890,10 @@ int main(int argc, char *argv[])
 
 
     const boundBox& bb = mesh.bounds();
-    const vector span = bb.span();
     const scalar mergeDim = mergeTol * bb.minDim();
 
     Info<< "Mesh bounding box : " << bb << nl
-        << "        with span : " << span << nl
+        << "        with span : " << bb.span() << nl
         << "Merge distance    : " << mergeDim << nl
         << endl;
 
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index ad3db870f4..a224a42ebd 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -577,7 +577,7 @@ bool Foam::backgroundMeshDecomposition::refineCell
 //            // weightEstimate += sampleVol/pow3(s);
 //        }
 //
-//        if (sqr(spanScale_)*sqr(minCellSize) < magSqr(cellBb.span()))
+//        if (sqr(spanScale_)*sqr(minCellSize) < cellBb.magSqr())
 //        {
 //            return true;
 //        }
@@ -695,7 +695,7 @@ void Foam::backgroundMeshDecomposition::buildPatchAndTree()
     Pstream::allGatherList(allBackgroundMeshBounds_);
 
     // find global bounding box
-    globalBackgroundBounds_ = treeBoundBox(boundBox::invertedBox);
+    globalBackgroundBounds_.reset();
     forAll(allBackgroundMeshBounds_, proci)
     {
         globalBackgroundBounds_.add(allBackgroundMeshBounds_[proci]);
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index 8591f6bb0c..57bcc5f761 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -1965,12 +1965,9 @@ void Foam::conformalVoronoiMesh::buildEdgeLocationTree
 {
     treeBoundBox overallBb
     (
-        geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4)
+        geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4, ROOTVSMALL)
     );
 
-    overallBb.min() -= Foam::point::uniform(ROOTVSMALL);
-    overallBb.max() += Foam::point::uniform(ROOTVSMALL);
-
     edgeLocationTreePtr_.reset
     (
         new dynamicIndexedOctree<dynamicTreeDataPoint>
@@ -1992,12 +1989,9 @@ void Foam::conformalVoronoiMesh::buildSurfacePtLocationTree
 {
     treeBoundBox overallBb
     (
-        geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4)
+        geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4, ROOTVSMALL)
     );
 
-    overallBb.min() -= Foam::point::uniform(ROOTVSMALL);
-    overallBb.max() += Foam::point::uniform(ROOTVSMALL);
-
     surfacePtLocationTreePtr_.reset
     (
         new dynamicIndexedOctree<dynamicTreeDataPoint>
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index a52d0c6dcc..57cb21f66e 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -564,12 +564,9 @@ Foam::conformationSurfaces::conformationSurfaces
 
     // Extend the global bounds to stop the bound box sitting on the surfaces
     // to be conformed to
-    //globalBounds_ = globalBounds_.extend(rndGen, 1e-4);
+    //globalBounds_.inflate(rndGen, 1e-4);
 
-    vector newSpan = 1e-4*globalBounds_.span();
-
-    globalBounds_.min() -= newSpan;
-    globalBounds_.max() += newSpan;
+    globalBounds_.grow(1e-4*globalBounds_.span());
 
     // Look at all surfaces at determine whether the locationInMesh point is
     // inside or outside each, to establish a signature for the domain to be
diff --git a/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshSurfaceSimplify/foamyHexMeshSurfaceSimplify_non_octree.C b/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshSurfaceSimplify/foamyHexMeshSurfaceSimplify_non_octree.C
index 1cfd57f9e8..442c83c700 100644
--- a/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshSurfaceSimplify/foamyHexMeshSurfaceSimplify_non_octree.C
+++ b/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshSurfaceSimplify/foamyHexMeshSurfaceSimplify_non_octree.C
@@ -122,11 +122,7 @@ int main(int argc, char *argv[])
 
     // Extend
     treeBoundBox bb = geometryToConformTo.globalBounds();
-    {
-        const vector smallVec = 0.1*bb.span();
-        bb.min() -= smallVec;
-        bb.max() += smallVec;
-    }
+    bb.grow(0.1*bb.span());
 
     Info<< "Meshing to bounding box " << bb << nl << endl;
 
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index 8aa4e36325..f461e45c43 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -871,17 +871,16 @@ Foam::label Foam::checkTopology
             Info<< "    "
                 << setw(20) << "PointZone"
                 << setw(8) << "Points"
-                << "BoundingBox" << endl;
+                << "BoundingBox" << nl;
 
             for (const auto& zone : pointZones)
             {
-                boundBox bb;
-                for (const label pointi : zone)
-                {
-                    bb.add(mesh.points()[pointi]);
-                }
-
-                bb.reduce();  // Global min/max
+                boundBox bb
+                (
+                    mesh.points(),
+                    static_cast<const labelUList&>(zone),
+                    true  // Reduce (global min/max)
+                );
 
                 Info<< "    "
                     << setw(20) << zone.name()
diff --git a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
index df38326ded..39664a65ad 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
@@ -4,10 +4,7 @@ Random rndGen(653213);
 List<treeBoundBox> meshBb
 (
     1,
-    treeBoundBox
-    (
-        boundBox(coarseMesh.points(), false)
-    ).extend(rndGen, 1e-3)
+    treeBoundBox(coarseMesh.points()).extend(rndGen, 1e-3)
 );
 
 // Dummy bounds dictionary
diff --git a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H
index 796e2a3dc3..66375be8b3 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H
@@ -4,10 +4,7 @@ Random rndGen(653213);
 List<treeBoundBox> meshBb
 (
     1,
-    treeBoundBox
-    (
-        boundBox(coarseMesh.points(), false)
-    ).extend(rndGen, 1e-3)
+    treeBoundBox(coarseMesh.points()).extend(rndGen, 1e-3)
 );
 
 // Dummy bounds dictionary
diff --git a/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C b/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C
index bc99d0e172..955680be35 100644
--- a/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C
+++ b/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C
@@ -178,12 +178,9 @@ void createBoundaryEdgeTrees
         // geometry there are less face/edge aligned items.
         treeBoundBox bb
         (
-            treeBoundBox(UList<point>(surf.localPoints())).extend(rndGen, 1e-4)
+            treeBoundBox(surf.localPoints()).extend(rndGen, 1e-4, ROOTVSMALL)
         );
 
-        bb.min() -= point::uniform(ROOTVSMALL);
-        bb.max() += point::uniform(ROOTVSMALL);
-
         bEdgeTrees.set
         (
             surfI,
diff --git a/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C b/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C
index 3b8552ea55..d5d887d973 100644
--- a/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C
+++ b/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C
@@ -193,10 +193,7 @@ int main(int argc, char *argv[])
         meshBb[Pstream::myProcNo()] = List<treeBoundBox>
         (
             1,
-            treeBoundBox
-            (
-                boundBox(mesh.points(), false)
-            ).extend(rndGen, 1e-3)
+            treeBoundBox(mesh.points()).extend(rndGen, 1e-3)
         );
         Pstream::allGatherList(meshBb);
     }
@@ -243,7 +240,7 @@ int main(int argc, char *argv[])
         }
         else
         {
-            bbs = List<treeBoundBox>(1, treeBoundBox(boundBox::invertedBox));
+            bbs = List<treeBoundBox>(1, treeBoundBox::null());
         }
 
         dictionary dict;
diff --git a/src/OpenFOAM/algorithms/AABBTree/AABBTree.C b/src/OpenFOAM/algorithms/AABBTree/AABBTree.C
index 4bd0231819..df4a4c48e4 100644
--- a/src/OpenFOAM/algorithms/AABBTree/AABBTree.C
+++ b/src/OpenFOAM/algorithms/AABBTree/AABBTree.C
@@ -195,10 +195,10 @@ void Foam::AABBTree<Type>::createBoxes
     // Assign the objects to min or max bin
 
     DynamicList<label> minBinObjectIDs(objectIDs.size());
-    treeBoundBox minBb(boundBox::invertedBox);
+    treeBoundBox minBb;
 
     DynamicList<label> maxBinObjectIDs(objectIDs.size());
-    treeBoundBox maxBb(boundBox::invertedBox);
+    treeBoundBox maxBb;
 
     for (const label objI : objectIDs)
     {
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index ddf6111a33..5a8a820d65 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -935,13 +935,10 @@ Foam::polyMesh::cellTree() const
 {
     if (!cellTreePtr_)
     {
-        treeBoundBox overallBb(points());
-
         Random rndGen(261782);
 
-        overallBb = overallBb.extend(rndGen, 1e-4);
-        overallBb.min() -= point::uniform(ROOTVSMALL);
-        overallBb.max() += point::uniform(ROOTVSMALL);
+        treeBoundBox overallBb(points());
+        overallBb.inflate(rndGen, 1e-4, ROOTVSMALL);
 
         cellTreePtr_.reset
         (
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSearch.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSearch.C
index fe01c0e8db..c76f1cb5d2 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSearch.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSearch.C
@@ -188,7 +188,7 @@ void Foam::PatchTools::calcBounds
     bitSet pointUsed(points.size());
 
     nPoints = 0;
-    bb = boundBox::invertedBox;
+    bb.reset();
 
     for (const auto& f : p)
     {
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
index 8113c2262d..966f278e07 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
+++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
@@ -129,7 +129,7 @@ void Foam::polyMeshGeometry::updateCellCentresAndVols
 
         // Estimate the cell centre and bounding box using the face centres
         vector cEst(Zero);
-        boundBox bb(boundBox::invertedBox);
+        boundBox bb;
 
         for (const label facei : cFaces)
         {
diff --git a/src/finiteArea/faMesh/faMeshDemandDrivenData.C b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
index 8be7175779..1d98ea822d 100644
--- a/src/finiteArea/faMesh/faMeshDemandDrivenData.C
+++ b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
@@ -1924,24 +1924,15 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(vectorField& result) const
 
                 // Remove duplicate points
                 vectorField allPoints(nAllPoints, Zero);
-                boundBox bb(allPointsExt, false);
-                scalar tol = 0.001*mag(bb.max() - bb.min());
+                const scalar tol = 0.001*treeBoundBox(allPointsExt).mag();
 
                 nAllPoints = 0;
-                forAll(allPointsExt, pI)
+                for (const point& pt : allPointsExt)
                 {
                     bool duplicate = false;
-                    for (label i=0; i<nAllPoints; ++i)
+                    for (label i = 0; i < nAllPoints; ++i)
                     {
-                        if
-                        (
-                            mag
-                            (
-                                allPoints[i]
-                              - allPointsExt[pI]
-                            )
-                          < tol
-                        )
+                        if (pt.dist(allPoints[i]) < tol)
                         {
                             duplicate = true;
                             break;
@@ -1950,7 +1941,7 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(vectorField& result) const
 
                     if (!duplicate)
                     {
-                        allPoints[nAllPoints++] = allPointsExt[pI];
+                        allPoints[nAllPoints++] = pt;
                     }
                 }
 
@@ -2069,8 +2060,7 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(vectorField& result) const
 
                 procLsPoints[Pstream::myProcNo()] = locPoints;
 
-                const boundBox bb(locPoints, false);
-                tol = 0.001*mag(bb.max() - bb.min());
+                tol = 0.001*treeBoundBox(locPoints).mag();
             }
 
             Pstream::allGatherList(procLsPoints);
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
index d618ec946b..8d39187889 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
@@ -626,7 +626,7 @@ turbulentDFSEMInletFvPatchVectorField
     triCumulativeMagSf_(),
     sumTriMagSf_(Pstream::nProcs() + 1, Zero),
     patchNormal_(Zero),
-    patchBounds_(boundBox::invertedBox),
+    patchBounds_(),
 
     eddies_(Zero),
     v0_(Zero),
@@ -709,7 +709,7 @@ turbulentDFSEMInletFvPatchVectorField
     triCumulativeMagSf_(),
     sumTriMagSf_(Pstream::nProcs() + 1, Zero),
     patchNormal_(Zero),
-    patchBounds_(boundBox::invertedBox),
+    patchBounds_(),
 
     eddies_(),
     v0_(Zero),
diff --git a/src/functionObjects/field/fieldExtents/fieldExtentsTemplates.C b/src/functionObjects/field/fieldExtents/fieldExtentsTemplates.C
index c568fc59be..90761a51bd 100644
--- a/src/functionObjects/field/fieldExtents/fieldExtentsTemplates.C
+++ b/src/functionObjects/field/fieldExtents/fieldExtentsTemplates.C
@@ -64,7 +64,7 @@ void Foam::functionObjects::fieldExtents::calcFieldExtents
 
     auto extents = [this](const scalarField& mask, const vectorField& C)
     {
-        boundBox extents(boundBox::invertedBox);
+        boundBox extents;
         forAll(mask, i)
         {
             if (mask[i] > 0.5)
@@ -77,7 +77,7 @@ void Foam::functionObjects::fieldExtents::calcFieldExtents
 
         if (extents.empty())
         {
-            extents.add(point::zero);
+            extents.reset(point::zero);
         }
 
         return extents;
@@ -97,7 +97,7 @@ void Foam::functionObjects::fieldExtents::calcFieldExtents
         Log << "    internal field: " << bb << nl;
         file() << bb;
 
-        this->setResult(fieldName + "_internal_min" , bb.min());
+        this->setResult(fieldName + "_internal_min", bb.min());
         this->setResult(fieldName + "_internal_max", bb.max());
     }
 
diff --git a/src/functionObjects/field/streamLine/streamLineBase.C b/src/functionObjects/field/streamLine/streamLineBase.C
index 71a476e83a..512fff95e1 100644
--- a/src/functionObjects/field/streamLine/streamLineBase.C
+++ b/src/functionObjects/field/streamLine/streamLineBase.C
@@ -927,7 +927,7 @@ bool Foam::functionObjects::streamLineBase::read(const dictionary& dict)
     }
 
 
-    bounds_ = boundBox::invertedBox;
+    bounds_.reset();
     if (dict.readIfPresent("bounds", bounds_) && !bounds_.empty())
     {
         Info<< "    clipping all segments to " << bounds_ << nl << endl;
diff --git a/src/functionObjects/randomProcesses/energySpectrum/energySpectrum.C b/src/functionObjects/randomProcesses/energySpectrum/energySpectrum.C
index 045fce032d..8e1d24af70 100644
--- a/src/functionObjects/randomProcesses/energySpectrum/energySpectrum.C
+++ b/src/functionObjects/randomProcesses/energySpectrum/energySpectrum.C
@@ -152,8 +152,8 @@ bool Foam::functionObjects::energySpectrum::read(const dictionary& dict)
     // Assume all cells are the same size...
     boundBox cellBb(mesh_.cellBb(0));
 
-    const vector L(meshBb.max() - meshBb.min());
-    const vector nCellXYZ(cmptDivide(L, cellBb.max() - cellBb.min()));
+    const vector L(meshBb.span());
+    const vector nCellXYZ(cmptDivide(L, cellBb.span()));
 
     N_ = Vector<int>
     (
@@ -163,7 +163,7 @@ bool Foam::functionObjects::energySpectrum::read(const dictionary& dict)
     );
 
     // Check that the mesh is a structured box
-    vector cellDx(cellBb.max() - cellBb.min());
+    vector cellDx(cellBb.span());
     vector expectedMax(N_.x()*cellDx.x(), N_.y()*cellDx.y(), N_.z()*cellDx.z());
     vector relativeSize(cmptDivide(L, expectedMax));
     for (direction i = 0; i < 3; ++i)
@@ -183,8 +183,8 @@ bool Foam::functionObjects::energySpectrum::read(const dictionary& dict)
     // Map into i-j-k co-ordinates
     const vectorField& C = mesh_.C();
     c0_ = returnReduce(min(C), minOp<vector>());
-    const vector cMax = returnReduce(max(C), maxOp<vector>());
-    deltaC_ = cMax - c0_;
+    deltaC_ = returnReduce(max(C), maxOp<vector>());
+    deltaC_ -= c0_;
 
     forAll(C, celli)
     {
diff --git a/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C
index 66fe7f4e43..3d24412823 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C
@@ -64,14 +64,13 @@ void Foam::surfaceDisplacementPointPatchVectorField::calcProjection
     // we're guaranteed to hit something.
 
     //- Per point projection vector:
-    const scalar projectLen = mag(mesh.bounds().max()-mesh.bounds().min());
+    const scalar projectLen = mesh.bounds().mag();
 
     // For case of fixed projection vector:
     vector projectVec(Zero);
     if (projectMode_ == FIXEDNORMAL)
     {
-        vector n = projectDir_/mag(projectDir_);
-        projectVec = projectLen*n;
+        projectVec = projectLen * normalised(projectDir_);
     }
 
 
diff --git a/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
index b560085117..645313a746 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
@@ -64,14 +64,13 @@ void Foam::surfaceSlipDisplacementPointPatchVectorField::calcProjection
     // we're guaranteed to hit something.
 
     //- Per point projection vector:
-    const scalar projectLen = mag(mesh.bounds().max()-mesh.bounds().min());
+    const scalar projectLen = mesh.bounds().mag();
 
     // For case of fixed projection vector:
-    vector projectVec(0, 0, 0);
+    vector projectVec(Zero);
     if (projectMode_ == FIXEDNORMAL)
     {
-        vector n = projectDir_/mag(projectDir_);
-        projectVec = projectLen*n;
+        projectVec = projectLen * normalised(projectDir_);
     }
 
 
diff --git a/src/lagrangian/basic/InteractionLists/InteractionLists.C b/src/lagrangian/basic/InteractionLists/InteractionLists.C
index 36ccf7184d..352f51224c 100644
--- a/src/lagrangian/basic/InteractionLists/InteractionLists.C
+++ b/src/lagrangian/basic/InteractionLists/InteractionLists.C
@@ -44,19 +44,12 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
 
     Random rndGen(419715);
 
-    const vector interactionVec = maxDistance_*vector::one;
-
-    treeBoundBox procBb(treeBoundBox(mesh_.points()));
-
-    treeBoundBox extendedProcBb
-    (
-        procBb.min() - interactionVec,
-        procBb.max() + interactionVec
-    );
+    treeBoundBox procBb(mesh_.points());
 
     treeBoundBoxList allExtendedProcBbs(Pstream::nProcs());
 
-    allExtendedProcBbs[Pstream::myProcNo()] = extendedProcBb;
+    allExtendedProcBbs[Pstream::myProcNo()] = procBb;
+    allExtendedProcBbs[Pstream::myProcNo()].grow(maxDistance_);
 
     Pstream::allGatherList(allExtendedProcBbs);
 
@@ -194,16 +187,11 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
             globalTransforms.transformIndex(ciat)
         );
 
-        treeBoundBox tempTransformedBb
+        treeBoundBox extendedBb
         (
             transform.invTransformPosition(cellBbsToExchange[bbI].points())
         );
-
-        treeBoundBox extendedBb
-        (
-            tempTransformedBb.min() - interactionVec,
-            tempTransformedBb.max() + interactionVec
-        );
+        extendedBb.grow(maxDistance_);
 
         // Find all elements intersecting box.
         labelList interactingElems
@@ -419,16 +407,11 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
             globalTransforms.transformIndex(wfiat)
         );
 
-        treeBoundBox tempTransformedBb
+        treeBoundBox extendedBb
         (
             transform.invTransformPosition(wallFaceBbsToExchange[bbI].points())
         );
-
-        treeBoundBox extendedBb
-        (
-            tempTransformedBb.min() - interactionVec,
-            tempTransformedBb.max() + interactionVec
-        );
+        extendedBb.grow(maxDistance_);
 
         // Find all elements intersecting box.
         labelList interactingElems
@@ -592,11 +575,8 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
     {
         const treeBoundBox& cellBb = cellBbs[celli];
 
-        treeBoundBox extendedBb
-        (
-            cellBb.min() - interactionVec,
-            cellBb.max() + interactionVec
-        );
+        treeBoundBox extendedBb(cellBb);
+        extendedBb.grow(maxDistance_);
 
         // Find all cells intersecting extendedBb
         labelList interactingElems(allCellsTree.findBox(extendedBb));
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblock.C b/src/mesh/blockMesh/PDRblockMesh/PDRblock.C
index 9559073305..c1b6f3b6f0 100644
--- a/src/mesh/blockMesh/PDRblockMesh/PDRblock.C
+++ b/src/mesh/blockMesh/PDRblockMesh/PDRblock.C
@@ -130,7 +130,7 @@ void Foam::PDRblock::adjustSizes()
         grid_.y().clear();
         grid_.z().clear();
 
-        bounds_ = boundBox::invertedBox;
+        bounds_.reset();
         edgeLimits_.min() = 0;
         edgeLimits_.max() = 0;
         return;
diff --git a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C
index 99f4db3011..f8edf82b95 100644
--- a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C
+++ b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C
@@ -94,7 +94,7 @@ Foam::blockVertices::projectVertex::operator point() const
         geometry_,
         surfaces_,
         start,
-        scalarField(start.size(), magSqr(bb.span())),
+        scalarField(start.size(), bb.magSqr()),
         boundaryNear,
         boundaryConstraint
     );
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
index 6ac13188b0..397fb18255 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
@@ -3228,10 +3228,11 @@ void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
         Random rndGen(653213);
 
         // Get local mesh bounding box. Single box for now.
-        List<treeBoundBox> meshBb(1);
-        treeBoundBox& bb = meshBb[0];
-        bb = treeBoundBox(mesh_.points());
-        bb = bb.extend(rndGen, 1e-4);
+        List<treeBoundBox> meshBb
+        (
+            1,
+            treeBoundBox(mesh_.points()).extend(rndGen, 1e-4)
+        );
 
         // Distribute all geometry (so refinementSurfaces and shellSurfaces)
         searchableSurfaces& geometry =
diff --git a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
index afef9a7de0..763d488eb9 100644
--- a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
+++ b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
@@ -297,9 +297,7 @@ void Foam::refinementFeatures::buildTrees(const label featI)
 
     // Slightly extended bb. Slightly off-centred just so on symmetric
     // geometry there are less face/edge aligned items.
-    bb = bb.extend(rndGen, 1e-4);
-    bb.min() -= point::uniform(ROOTVSMALL);
-    bb.max() += point::uniform(ROOTVSMALL);
+    bb.inflate(rndGen, 1e-4, ROOTVSMALL);
 
     edgeTrees_.set
     (
@@ -438,9 +436,7 @@ Foam::refinementFeatures::regionEdgeTrees() const
 
             // Slightly extended bb. Slightly off-centred just so on symmetric
             // geometry there are less face/edge aligned items.
-            bb = bb.extend(rndGen, 1e-4);
-            bb.min() -= point::uniform(ROOTVSMALL);
-            bb.max() += point::uniform(ROOTVSMALL);
+            bb.inflate(rndGen, 1e-4, ROOTVSMALL);
 
             trees.set
             (
diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
index 9340842d8c..198be01bcf 100644
--- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
+++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
@@ -160,7 +160,7 @@ void Foam::shellSurfaces::setAndCheckLevels
 void Foam::shellSurfaces::orient()
 {
     // Determine outside point.
-    boundBox overallBb = boundBox::invertedBox;
+    boundBox overallBb;
 
     bool hasSurface = false;
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C
index 6c3ee03d07..1e8d1f287c 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C
@@ -3172,9 +3172,7 @@ void Foam::snappySnapDriver::reverseAttractMeshPoints
 
         // Slightly extended bb. Slightly off-centred just so on symmetric
         // geometry there are less face/edge aligned items.
-        bb = bb.extend(rndGen, 1e-4);
-        bb.min() -= point::uniform(ROOTVSMALL);
-        bb.max() += point::uniform(ROOTVSMALL);
+        bb.inflate(rndGen, 1e-4, ROOTVSMALL);
     }
 
     // Collect candidate points for attraction
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C
index 6ff8f94379..bc27474b66 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C
@@ -300,7 +300,7 @@ Foam::label Foam::advancingFrontAMI::findTargetFace
         srcFacePti == -1 ? bb.centre() : srcPts[srcFace[srcFacePti]];
 
     const pointIndexHit sample =
-        treePtr_->findNearest(srcPt, magSqr(bb.max() - bb.centre()), fnOp);
+        treePtr_->findNearest(srcPt, 0.25*bb.magSqr(), fnOp);
 
     if (sample.hit() && isCandidate(srcFacei, sample.index()))
     {
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C b/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C
index da6e87b50c..0281437c53 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C
@@ -337,9 +337,6 @@ bool Foam::faceAreaWeightAMI2D::calculate
     // Find an approximate face length scale
     const scalar lf(Cbb_*Foam::sqrt(gAverage(srcMagSf_)));
 
-    // Expansion to apply to source face bounding box
-    const vector d(lf*vector::one);
-
     if (validSize)
     {
         // Create the tgt tree
@@ -371,8 +368,7 @@ bool Foam::faceAreaWeightAMI2D::calculate
             const face& srcFace = src[srcFacei];
 
             treeBoundBox srcFaceBb(srcPoints, srcFace);
-            srcFaceBb.min() -= d;
-            srcFaceBb.max() += d;
+            srcFaceBb.grow(lf);
 
             const labelList tgtFaces
             (
diff --git a/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C b/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
index ec030926ae..82a018e2d0 100644
--- a/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
+++ b/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
@@ -833,12 +833,9 @@ Foam::extendedEdgeMesh::pointTree() const
         // geometry there are less face/edge aligned items.
         treeBoundBox bb
         (
-            treeBoundBox(points()).extend(rndGen, 1e-4)
+            treeBoundBox(points()).extend(rndGen, 1e-4, ROOTVSMALL)
         );
 
-        bb.min() -= point::uniform(ROOTVSMALL);
-        bb.max() += point::uniform(ROOTVSMALL);
-
         const labelList featurePointLabels = identity(nonFeatureStart_);
 
         pointTree_.reset
@@ -873,12 +870,9 @@ Foam::extendedEdgeMesh::edgeTree() const
         // geometry there are less face/edge aligned items.
         treeBoundBox bb
         (
-            treeBoundBox(points()).extend(rndGen, 1e-4)
+            treeBoundBox(points()).extend(rndGen, 1e-4, ROOTVSMALL)
         );
 
-        bb.min() -= point::uniform(ROOTVSMALL);
-        bb.max() += point::uniform(ROOTVSMALL);
-
         labelList allEdges(identity(edges().size()));
 
         edgeTree_.reset
@@ -915,12 +909,9 @@ Foam::extendedEdgeMesh::edgeTreesByType() const
         // geometry there are less face/edge aligned items.
         treeBoundBox bb
         (
-            treeBoundBox(points()).extend(rndGen, 1e-4)
+            treeBoundBox(points()).extend(rndGen, 1e-4, ROOTVSMALL)
         );
 
-        bb.min() -= point::uniform(ROOTVSMALL);
-        bb.max() += point::uniform(ROOTVSMALL);
-
         labelListList sliceEdges(nEdgeTypes);
 
         // External edges
diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
index b3a1a45a1c..383f6f7932 100644
--- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
+++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
@@ -337,7 +337,7 @@ void Foam::mappedPatchBase::findLocalSamples
                         cc,
                         celli
                     );
-                    near.first().second().first() = magSqr(cc-sample);
+                    near.first().second().first() = sample.distSqr(cc);
                     near.first().second().second() = myRank;
                     near.second() = mySampleWorld;
                 }
@@ -396,11 +396,10 @@ void Foam::mappedPatchBase::findLocalSamples
                     treeBoundBox(pp.points(), pp.meshPoints()).extend
                     (
                         rndGen,
-                        1e-4
+                        1e-4,
+                        ROOTVSMALL
                     )
                 );
-                patchBb.min() -= point::uniform(ROOTVSMALL);
-                patchBb.max() += point::uniform(ROOTVSMALL);
 
                 indexedOctree<treeDataFace> boundaryTree
                 (
@@ -425,7 +424,7 @@ void Foam::mappedPatchBase::findLocalSamples
                     nearInfo = boundaryTree.findNearest
                     (
                         sample,
-                        magSqr(patchBb.span())
+                        patchBb.magSqr()
                     );
 
                     if (!nearInfo.hit())
@@ -438,7 +437,7 @@ void Foam::mappedPatchBase::findLocalSamples
                     {
                         point fc(pp[nearInfo.index()].centre(pp.points()));
                         nearInfo.setPoint(fc);
-                        near.first().second().first() = magSqr(fc-sample);
+                        near.first().second().first() = sample.distSqr(fc);
                         near.first().second().second() = myRank;
                         near.second() = mySampleWorld;
                     }
@@ -471,11 +470,10 @@ void Foam::mappedPatchBase::findLocalSamples
                     treeBoundBox(pp.points(), pp.meshPoints()).extend
                     (
                         rndGen,
-                        1e-4
+                        1e-4,
+                        ROOTVSMALL
                     )
                 );
-                patchBb.min() -= point::uniform(ROOTVSMALL);
-                patchBb.max() += point::uniform(ROOTVSMALL);
 
                 indexedOctree<treeDataPoint> boundaryTree
                 (
@@ -499,7 +497,7 @@ void Foam::mappedPatchBase::findLocalSamples
                     nearInfo = boundaryTree.findNearest
                     (
                         sample,
-                        magSqr(patchBb.span())
+                        patchBb.magSqr()
                     );
 
                     if (!nearInfo.hit())
@@ -512,7 +510,7 @@ void Foam::mappedPatchBase::findLocalSamples
                     {
                         const point& pt = nearInfo.point();
 
-                        near.first().second().first() = magSqr(pt-sample);
+                        near.first().second().first() = sample.distSqr(pt);
                         near.first().second().second() = myRank;
                         near.second() = mySampleWorld;
                     }
@@ -552,7 +550,7 @@ void Foam::mappedPatchBase::findLocalSamples
                     const point& fc = mesh.faceCentres()[facei];
 
                     near.first().first() = pointIndexHit(true, fc, facei);
-                    near.first().second().first() = magSqr(fc-sample);
+                    near.first().second().first() = sample.distSqr(fc);
                     near.first().second().second() = myRank;
                     near.second() = mySampleWorld;
                 }
diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C
index 2f840ac26d..8a2240d06e 100644
--- a/src/meshTools/meshSearch/meshSearch.C
+++ b/src/meshTools/meshSearch/meshSearch.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -136,7 +136,7 @@ bool Foam::meshSearch::findNearer
 
     forAll(points, pointi)
     {
-        scalar distSqr = magSqr(points[pointi] - sample);
+        scalar distSqr = sample.distSqr(points[pointi]);
 
         if (distSqr < nearestDistSqr)
         {
@@ -165,7 +165,7 @@ bool Foam::meshSearch::findNearer
     {
         label pointi = indices[i];
 
-        scalar distSqr = magSqr(points[pointi] - sample);
+        scalar distSqr = sample.distSqr(points[pointi]);
 
         if (distSqr < nearestDistSqr)
         {
@@ -184,11 +184,7 @@ Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const
 {
     const indexedOctree<treeDataCell>& tree = cellTree();
 
-    pointIndexHit info = tree.findNearest
-    (
-        location,
-        magSqr(tree.bb().max()-tree.bb().min())
-    );
+    pointIndexHit info = tree.findNearest(location, tree.bb().magSqr());
 
     if (!info.hit())
     {
@@ -203,7 +199,7 @@ Foam::label Foam::meshSearch::findNearestCellLinear(const point& location) const
     const vectorField& centres = mesh_.cellCentres();
 
     label nearestIndex = 0;
-    scalar minProximity = magSqr(centres[nearestIndex] - location);
+    scalar minProximity = location.distSqr(centres[nearestIndex]);
 
     findNearer
     (
@@ -232,7 +228,7 @@ Foam::label Foam::meshSearch::findNearestCellWalk
     // Walk in direction of face that decreases distance
 
     label curCelli = seedCelli;
-    scalar distanceSqr = magSqr(mesh_.cellCentres()[curCelli] - location);
+    scalar distanceSqr = location.distSqr(mesh_.cellCentres()[curCelli]);
 
     bool closer;
 
@@ -259,11 +255,7 @@ Foam::label Foam::meshSearch::findNearestFaceTree(const point& location) const
     const indexedOctree<treeDataCell>& tree = cellTree();
 
     // Search with decent span
-    pointIndexHit info = tree.findNearest
-    (
-        location,
-        magSqr(tree.bb().max()-tree.bb().min())
-    );
+    pointIndexHit info = tree.findNearest(location, tree.bb().magSqr());
 
     if (!info.hit())
     {
@@ -277,7 +269,7 @@ Foam::label Foam::meshSearch::findNearestFaceTree(const point& location) const
     const cell& ownFaces = mesh_.cells()[info.index()];
 
     label nearestFacei = ownFaces[0];
-    scalar minProximity = magSqr(centres[nearestFacei] - location);
+    scalar minProximity = location.distSqr(centres[nearestFacei]);
 
     findNearer
     (
@@ -297,7 +289,7 @@ Foam::label Foam::meshSearch::findNearestFaceLinear(const point& location) const
     const vectorField& centres = mesh_.faceCentres();
 
     label nearestFacei = 0;
-    scalar minProximity = magSqr(centres[nearestFacei] - location);
+    scalar minProximity = location.distSqr(centres[nearestFacei]);
 
     findNearer
     (
@@ -329,7 +321,7 @@ Foam::label Foam::meshSearch::findNearestFaceWalk
     // Walk in direction of face that decreases distance
 
     label curFacei = seedFacei;
-    scalar distanceSqr = magSqr(centres[curFacei] - location);
+    scalar distanceSqr = location.distSqr(centres[curFacei]);
 
     while (true)
     {
@@ -415,7 +407,7 @@ Foam::label Foam::meshSearch::findCellWalk
 
     // Walk in direction of face that decreases distance
     label curCelli = seedCelli;
-    scalar nearestDistSqr = magSqr(mesh_.cellCentres()[curCelli] - location);
+    scalar nearestDistSqr = location.distSqr(mesh_.cellCentres()[curCelli]);
 
     while(true)
     {
@@ -444,7 +436,7 @@ Foam::label Foam::meshSearch::findCellWalk
                 }
 
                 // Also calculate the nearest cell
-                scalar distSqr = magSqr(mesh_.cellCentres()[celli] - location);
+                scalar distSqr = location.distSqr(mesh_.cellCentres()[celli]);
 
                 if (distSqr < nearestDistSqr)
                 {
@@ -615,9 +607,7 @@ const Foam::treeBoundBox& Foam::meshSearch::dataBoundBox() const
         treeBoundBox& overallBb = overallBbPtr_();
 
         // Extend slightly and make 3D
-        overallBb = overallBb.extend(rndGen, 1e-4);
-        overallBb.min() -= point::uniform(ROOTVSMALL);
-        overallBb.max() += point::uniform(ROOTVSMALL);
+        overallBb.inflate(rndGen, 1e-4, ROOTVSMALL);
     }
 
     return *overallBbPtr_;
@@ -816,7 +806,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFace
             pointIndexHit info = boundaryTree().findNearest
             (
                 location,
-                magSqr(tree.bb().max()-tree.bb().min())
+                tree.bb().magSqr()
             );
 
             if (!info.hit())
diff --git a/src/meshTools/processorLOD/cellBox/cellBox.C b/src/meshTools/processorLOD/cellBox/cellBox.C
index 56900afa5b..df9c0e190c 100644
--- a/src/meshTools/processorLOD/cellBox/cellBox.C
+++ b/src/meshTools/processorLOD/cellBox/cellBox.C
@@ -44,12 +44,11 @@ Foam::boundBox Foam::processorLODs::cellBox::calcSrcBox
     const label srcObji
 ) const
 {
-    const UList<label>& cellFaces = srcCells_[srcObji];
+    boundBox bb;
 
-    boundBox bb(srcPoints_, srcFaces_[cellFaces[0]], false);
-    for (label i = 1; i < cellFaces.size(); ++i)
+    for (const label facei : srcCells_[srcObji])
     {
-        bb.add(srcPoints_, srcFaces_[cellFaces[i]]);
+        bb.add(srcPoints_, srcFaces_[facei]);
     }
 
     return bb;
@@ -61,12 +60,11 @@ Foam::boundBox Foam::processorLODs::cellBox::calcTgtBox
     const label tgtObji
 ) const
 {
-    const UList<label>& cellFaces = tgtCells_[tgtObji];
+    boundBox bb;
 
-    boundBox bb(tgtPoints_, tgtFaces_[cellFaces[0]], false);
-    for (label i = 1; i < cellFaces.size(); ++i)
+    for (const label facei : tgtCells_[tgtObji])
     {
-        bb.add(tgtPoints_, tgtFaces_[cellFaces[i]]);
+        bb.add(tgtPoints_, tgtFaces_[facei]);
     }
 
     return bb;
diff --git a/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C b/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C
index ad5e8d3c9c..c58a76fdf0 100644
--- a/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C
+++ b/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C
@@ -89,19 +89,16 @@ Foam::searchableExtrudedCircle::searchableExtrudedCircle
     const edgeList& edges = eMesh.edges();
     bounds() = boundBox(points, false);
 
-    vector halfSpan(0.5*bounds().span());
-    point ctr(bounds().centre());
+    // Make the boundBox into a perfect cube around its centre
+    const scalar halfWidth = mag(0.5*bounds().span());
 
-    bounds().min() = ctr - mag(halfSpan) * vector::one;
-    bounds().max() = ctr + mag(halfSpan) * vector::one;
-
-    // Calculate bb of all points
-    treeBoundBox bb(bounds());
+    bounds().reset(bounds().centre());
+    bounds().grow(halfWidth);
 
     // Slightly extended bb. Slightly off-centred just so on symmetric
     // geometry there are less face/edge aligned items.
-    bb.min() -= point::uniform(ROOTVSMALL);
-    bb.max() += point::uniform(ROOTVSMALL);
+    treeBoundBox bb(bounds());
+    bb.grow(ROOTVSMALL);
 
     edgeTree_.reset
     (
@@ -230,7 +227,7 @@ void Foam::searchableExtrudedCircle::findParametricNearest
     const pointField& points = mesh.points();
     const labelListList& pointEdges = mesh.pointEdges();
 
-    const scalar maxDistSqr(Foam::magSqr(bounds().span()));
+    const scalar maxDistSqr = bounds().magSqr();
 
     // Normalise lambdas
     const scalarField lambdas
@@ -446,16 +443,14 @@ void Foam::searchableExtrudedCircle::getNormal
     normal.setSize(info.size());
     normal = Zero;
 
+    const scalar distSqr = bounds().magSqr();
+
     forAll(info, i)
     {
         if (info[i].hit())
         {
             // Find nearest on curve
-            pointIndexHit curvePt = tree.findNearest
-            (
-                info[i].point(),
-                Foam::magSqr(bounds().span())
-            );
+            pointIndexHit curvePt = tree.findNearest(info[i].point(), distSqr);
 
             normal[i] = info[i].hitPoint()-curvePt.hitPoint();
 
diff --git a/src/meshTools/searchableSurfaces/searchableSurfaceCollection/searchableSurfaceCollection.C b/src/meshTools/searchableSurfaces/searchableSurfaceCollection/searchableSurfaceCollection.C
index f3e8ce0b5e..9d47c4cd41 100644
--- a/src/meshTools/searchableSurfaces/searchableSurfaceCollection/searchableSurfaceCollection.C
+++ b/src/meshTools/searchableSurfaces/searchableSurfaceCollection/searchableSurfaceCollection.C
@@ -100,7 +100,7 @@ void Foam::searchableSurfaceCollection::findNearest
                     )
                 );
 
-                scalar distSqr = magSqr(globalPt - samples[pointi]);
+                scalar distSqr = globalPt.distSqr(samples[pointi]);
 
                 if (distSqr < minDistSqr[pointi])
                 {
@@ -249,7 +249,7 @@ Foam::searchableSurfaceCollection::searchableSurfaceCollection
     indexOffset_.setSize(surfI+1);
 
     // Bounds is the overall bounds - prepare for min/max ops
-    bounds() = boundBox::invertedBox;
+    bounds().reset();
 
     forAll(subGeom_, surfI)
     {
diff --git a/src/meshTools/searchableSurfaces/searchableSurfacesQueries/searchableSurfacesQueries.C b/src/meshTools/searchableSurfaces/searchableSurfacesQueries/searchableSurfacesQueries.C
index 34c17f8d5a..03266ce8e4 100644
--- a/src/meshTools/searchableSurfaces/searchableSurfacesQueries/searchableSurfacesQueries.C
+++ b/src/meshTools/searchableSurfaces/searchableSurfacesQueries/searchableSurfacesQueries.C
@@ -695,7 +695,7 @@ Foam::boundBox Foam::searchableSurfacesQueries::bounds
     const labelUList& surfacesToTest
 )
 {
-    boundBox bb(boundBox::invertedBox);
+    boundBox bb;
 
     for (const label surfi : surfacesToTest)
     {
diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
index 1acd5127d0..4714300105 100644
--- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
@@ -304,9 +304,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
     surfaceClosed_(-1),
     outsideVolType_(volumeType::UNKNOWN)
 {
-    const pointField& pts = triSurface::points();
-
-    bounds() = boundBox(pts, false);
+    bounds() = boundBox(triSurface::points(), false);
 }
 
 
@@ -334,9 +332,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
     surfaceClosed_(-1),
     outsideVolType_(volumeType::UNKNOWN)
 {
-    const pointField& pts = triSurface::points();
-
-    bounds() = boundBox(pts, false);
+    bounds() = boundBox(triSurface::points(), false);
 }
 
 
@@ -387,9 +383,7 @@ Foam::triSurfaceMesh::triSurfaceMesh
             << " : using scale " << scaleFactor << endl;
     }
 
-    const pointField& pts = triSurface::points();
-
-    bounds() = boundBox(pts, false);
+    bounds() = boundBox(triSurface::points(), false);
 
     // Have optional minimum quality for normal calculation
     if (dict.readIfPresent("minQuality", minQuality_) && minQuality_ > 0)
@@ -494,8 +488,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const readAction r)
         }
     }
 
-    const pointField& pts = triSurface::points();
-    bounds() = boundBox(pts, false);
+    bounds() = boundBox(triSurface::points(), false);
 }
 
 
@@ -629,8 +622,7 @@ Foam::triSurfaceMesh::triSurfaceMesh
     }
 
 
-    const pointField& pts = triSurface::points();
-    bounds() = boundBox(pts, false);
+    bounds() = boundBox(triSurface::points(), false);
 
     // Have optional minimum quality for normal calculation
     if (dict.readIfPresent("minQuality", minQuality_) && minQuality_ > 0)
@@ -706,8 +698,7 @@ void Foam::triSurfaceMesh::boundingSpheres
         const point& fc = centres[facei];
         for (const label pointi : f)
         {
-            const point& pt = pts[pointi];
-            radiusSqr[facei] = max(radiusSqr[facei], Foam::magSqr(fc-pt));
+            radiusSqr[facei] = max(radiusSqr[facei], fc.distSqr(pts[pointi]));
         }
     }
 
@@ -782,7 +773,7 @@ Foam::triSurfaceMesh::edgeTree() const
             identity(nEdges() - nInternalEdges(), nInternalEdges())
         );
 
-        treeBoundBox bb(Zero, Zero);
+        treeBoundBox bb(point::zero);
 
         if (bEdges.size())
         {
@@ -800,9 +791,7 @@ Foam::triSurfaceMesh::edgeTree() const
             // Slightly extended bb. Slightly off-centred just so on symmetric
             // geometry there are less face/edge aligned items.
 
-            bb = bb.extend(rndGen, 1e-4);
-            bb.min() -= point::uniform(ROOTVSMALL);
-            bb.max() += point::uniform(ROOTVSMALL);
+            bb.inflate(rndGen, 1e-4, ROOTVSMALL);
         }
 
 
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
index d9dfb971ca..54e5c35245 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
@@ -126,7 +126,7 @@ Foam::triSurfaceRegionSearch::treeByRegion() const
             );
 
             // Calculate bb without constructing local point numbering.
-            treeBoundBox bb(Zero, Zero);
+            treeBoundBox bb(point::zero);
 
             if (indirectRegionPatches_[regionI].size())
             {
@@ -156,9 +156,7 @@ Foam::triSurfaceRegionSearch::treeByRegion() const
                 // Slightly extended bb. Slightly off-centred just so
                 // on symmetric geometry there are fewer face/edge
                 // aligned items.
-                bb = bb.extend(rndGen, 1e-4);
-                bb.min() -= point::uniform(ROOTVSMALL);
-                bb.max() += point::uniform(ROOTVSMALL);
+                bb.inflate(rndGen, 1e-4, ROOTVSMALL);
             }
 
             treeByRegion_.set
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
index c281014c8b..63853a2d3b 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -208,7 +208,7 @@ Foam::triSurfaceSearch::tree() const
     if (!treePtr_)
     {
         // Calculate bb without constructing local point numbering.
-        treeBoundBox bb(Zero, Zero);
+        treeBoundBox bb(point::zero);
 
         if (surface().size())
         {
@@ -230,9 +230,7 @@ Foam::triSurfaceSearch::tree() const
 
             // Slightly extended bb. Slightly off-centred just so on symmetric
             // geometry there are less face/edge aligned items.
-            bb = bb.extend(rndGen, 1e-4);
-            bb.min() -= point::uniform(ROOTVSMALL);
-            bb.max() += point::uniform(ROOTVSMALL);
+            bb.inflate(rndGen, 1e-4, ROOTVSMALL);
         }
 
         const scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
diff --git a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
index 37fea38605..29f5d55102 100644
--- a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
+++ b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
@@ -195,8 +195,7 @@ void Foam::cellCellStencils::inverseDistance::markBoundaries
 
                 // Mark in voxel mesh
                 boundBox faceBb(pp.points(), pp[i], false);
-                faceBb.min() -= smallVec;
-                faceBb.max() += smallVec;
+                faceBb.grow(smallVec);
 
                 if (bb.overlaps(faceBb))
                 {
@@ -223,8 +222,7 @@ void Foam::cellCellStencils::inverseDistance::markBoundaries
 
                 // Mark in voxel mesh
                 boundBox faceBb(pp.points(), pp[i], false);
-                faceBb.min() -= smallVec;
-                faceBb.max() += smallVec;
+                faceBb.grow(smallVec);
 
                 if (bb.overlaps(faceBb))
                 {
@@ -313,8 +311,7 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
             {
                 label celli = tgtCellMap[tgtCelli];
                 treeBoundBox cBb(mesh_.cellBb(celli));
-                cBb.min() -= smallVec_;
-                cBb.max() += smallVec_;
+                cBb.grow(smallVec_);
 
                 if
                 (
@@ -383,8 +380,7 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
                 {
                     label celli = tgtCellMap[tgtCelli];
                     treeBoundBox cBb(mesh_.cellBb(celli));
-                    cBb.min() -= smallVec_;
-                    cBb.max() += smallVec_;
+                    cBb.grow(smallVec_);
 
                     if
                     (
@@ -542,8 +538,7 @@ void Foam::cellCellStencils::inverseDistance::markDonors
         if (srcOverlapProcs.size())
         {
             treeBoundBox subBb(mesh_.cellBb(celli));
-            subBb.min() -= smallVec_;
-            subBb.max() += smallVec_;
+            subBb.grow(smallVec_);
 
             forAll(srcOverlapProcs, i)
             {
@@ -1818,7 +1813,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
     }
     Pstream::listCombineReduce(nCellsPerZone, plusEqOp<label>());
 
-    const boundBox& allBb(mesh_.bounds());
+    const boundBox& allBb = mesh_.bounds();
 
     PtrList<fvMeshSubset> meshParts(nZones);
     List<treeBoundBoxList> meshBb(nZones);
diff --git a/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C b/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
index 515744366d..3ec9093fd7 100644
--- a/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
+++ b/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
@@ -91,8 +91,7 @@ bool Foam::cellCellStencils::trackingInverseDistance::markBoundaries
 
                 // Mark in voxel mesh
                 boundBox faceBb(pp.points(), pp[i], false);
-                faceBb.min() -= smallVec;
-                faceBb.max() += smallVec;
+                faceBb.grow(smallVec);
 
                 if (bb.overlaps(faceBb))
                 {
@@ -127,8 +126,8 @@ bool Foam::cellCellStencils::trackingInverseDistance::markBoundaries
 
                 // Mark in voxel mesh
                 boundBox faceBb(pp.points(), pp[i], false);
-                faceBb.min() -= smallVec;
-                faceBb.max() += smallVec;
+                faceBb.grow(smallVec);
+
                 if (!bb.contains(faceCentres[i]))
                 {
                     if (!hasWarned)
@@ -224,8 +223,7 @@ void Foam::cellCellStencils::trackingInverseDistance::markPatchesAsHoles
             {
                 label celli = tgtCellMap[tgtCelli];
                 boundBox cBb(mesh_.cellBb(celli));
-                cBb.min() -= smallVec_;
-                cBb.max() += smallVec_;
+                cBb.grow(smallVec_);
 
                 if
                 (
@@ -293,8 +291,7 @@ void Foam::cellCellStencils::trackingInverseDistance::markPatchesAsHoles
                 {
                     label celli = tgtCellMap[tgtCelli];
                     boundBox cBb(mesh_.cellBb(celli));
-                    cBb.min() -= smallVec_;
-                    cBb.max() += smallVec_;
+                    cBb.grow(smallVec_);
 
                     if
                     (
@@ -410,8 +407,7 @@ void Foam::cellCellStencils::trackingInverseDistance::markDonors
         if (srcOverlapProcs.size())
         {
             treeBoundBox subBb(mesh_.cellBb(celli));
-            subBb.min() -= smallVec_;
-            subBb.max() += smallVec_;
+            subBb.grow(smallVec_);
 
             forAll(srcOverlapProcs, i)
             {
@@ -600,7 +596,7 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
             << exit(FatalIOError);
     }
 
-    const boundBox& allBb(mesh_.bounds());
+    const boundBox& allBb = mesh_.bounds();
 
     List<treeBoundBoxList> meshBb(nZones);
 
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
index 80c7cb5332..107be52c0b 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
@@ -950,7 +950,7 @@ Foam::distributedTriSurfaceMesh::findBestProcs
                 // Minimum search distance to find the triangle
                 point near, far;
                 bbs[bbi].calcExtremities(centre, near, far);
-                minDistSqr = min(minDistSqr, magSqr(centre-far));
+                minDistSqr = min(minDistSqr, centre.distSqr(far));
             }
         }
     }
@@ -977,11 +977,11 @@ Foam::distributedTriSurfaceMesh::findBestProcs
                     point near, far;
                     bbs[bbi].calcExtremities(centre, near, far);
 
-                    scalar d2 = magSqr(centre-near);
+                    scalar d2 = centre.distSqr(near);
                     if (d2 < minDistSqr)
                     {
                         minDistSqr = d2;
-                        maxDistSqr = min(radiusSqr, magSqr(centre-far));
+                        maxDistSqr = min(radiusSqr, centre.distSqr(far));
                         minProci = proci;
                     }
                 }
@@ -1748,12 +1748,11 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 
 
     // Initialise to inverted box
-    List<List<treeBoundBox>> bbs(Pstream::nProcs());
-    forAll(bbs, proci)
-    {
-        bbs[proci].setSize(1, treeBoundBox(boundBox::invertedBox));
-    }
-
+    List<List<treeBoundBox>> bbs
+    (
+        Pstream::nProcs(),
+        List<treeBoundBox>(1, treeBoundBox::null())
+    );
 
     const globalIndex& triIndexer = globalTris();
 
@@ -3735,7 +3734,7 @@ void Foam::distributedTriSurfaceMesh::findLineAll
     const vectorField smallVec
     (
         ROOTSMALL*dirVec
-      + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
+      + vector::uniform(ROOTVSMALL)
     );
 
     // Copy to input and compact any hits
diff --git a/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.C b/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.C
index 3fc1ac3f25..991a43b162 100644
--- a/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.C
+++ b/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.C
@@ -52,13 +52,13 @@ Foam::patchDistMethods::exact::patchSurface() const
 
         Random rndGen(0);
 
-        boundBox localBb(mesh_.points(), false);
+        treeBoundBox localBb(mesh_.points());
 
         // Determine mesh bounding boxes:
         List<treeBoundBox> meshBb
         (
             1,
-            treeBoundBox(localBb).extend(rndGen, 1E-3)
+            localBb.extend(rndGen, 1E-3)
         );
 
         // Dummy bounds dictionary
diff --git a/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
index 3d0034a8c0..513ba5f1ff 100644
--- a/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
+++ b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
@@ -46,12 +46,8 @@ Foam::scalar Foam::meshToMeshMethod::tolerance_ = 1e-6;
 
 Foam::labelList Foam::meshToMeshMethod::maskCells() const
 {
-    boundBox intersectBb
-    (
-        max(src_.bounds().min(), tgt_.bounds().min()),
-        min(src_.bounds().max(), tgt_.bounds().max())
-    );
-
+    boundBox intersectBb(src_.bounds());
+    intersectBb &= tgt_.bounds();
     intersectBb.inflate(0.01);
 
     DynamicList<label> cells(src_.nCells());
diff --git a/src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C b/src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
index 6ca92bc00f..b644db13a6 100644
--- a/src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
+++ b/src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
@@ -77,20 +77,16 @@ void Foam::meshToMesh0::calcAddressing()
     }
 
     treeBoundBox meshBb(fromPoints);
+    treeBoundBox shiftedBb(meshBb);
 
     scalar typDim = meshBb.avgDim()/(2.0*cbrt(scalar(fromCells.size())));
-
-    treeBoundBox shiftedBb
-    (
-        meshBb.min(),
-        meshBb.max() + vector(typDim, typDim, typDim)
-    );
+    shiftedBb.max() += vector::uniform(typDim);
 
     DebugInfo
         << "\nMesh" << nl
         << "   bounding box           : " << meshBb << nl
         << "   bounding box (shifted) : " << shiftedBb << nl
-        << "   typical dimension      : " << shiftedBb.typDim() << endl;
+        << "   typical dimension      : " << shiftedBb.avgDim() << endl;
 
     indexedOctree<treeDataCell> oc
     (
@@ -155,14 +151,12 @@ void Foam::meshToMesh0::calcAddressing()
             else
             {
                 treeBoundBox wallBb(fromPatch.localPoints());
+                treeBoundBox shiftedBb(wallBb);
+
                 scalar typDim =
                     wallBb.avgDim()/(2.0*sqrt(scalar(fromPatch.size())));
 
-                treeBoundBox shiftedBb
-                (
-                    wallBb.min(),
-                    wallBb.max() + vector(typDim, typDim, typDim)
-                );
+                shiftedBb.max() += vector::uniform(typDim);
 
                 // Note: allow more levels than in meshSearch. Assume patch
                 // is not as big as all boundary faces
@@ -180,7 +174,7 @@ void Foam::meshToMesh0::calcAddressing()
 
                 boundaryAddressing_[patchi].setSize(toPatch.size());
 
-                scalar distSqr = sqr(wallBb.mag());
+                scalar distSqr = wallBb.magSqr();
 
                 forAll(toPatch, toi)
                 {
@@ -229,7 +223,7 @@ void Foam::meshToMesh0::cellAddresses
         const vector& p = points[toI];
 
         // set the sqr-distance
-        scalar distSqr = magSqr(p - centresFrom[curCell]);
+        scalar distSqr = p.distSqr(centresFrom[curCell]);
 
         bool closer;
 
@@ -242,8 +236,7 @@ void Foam::meshToMesh0::cellAddresses
 
             forAll(neighbours, nI)
             {
-                scalar curDistSqr =
-                    magSqr(p - centresFrom[neighbours[nI]]);
+                scalar curDistSqr = p.distSqr(centresFrom[neighbours[nI]]);
 
                 // search through all the neighbours.
                 // If the cell is closer, reset current cell and distance
diff --git a/src/sampling/probes/patchProbes.C b/src/sampling/probes/patchProbes.C
index 5f50142d24..9e325f2901 100644
--- a/src/sampling/probes/patchProbes.C
+++ b/src/sampling/probes/patchProbes.C
@@ -71,7 +71,7 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
     {
         // Collect mesh faces and bounding box
         labelList bndFaces(nFaces);
-        treeBoundBox overallBb(boundBox::invertedBox);
+        treeBoundBox overallBb;
 
         nFaces = 0;
         forAll(patchIDs, i)
@@ -88,9 +88,7 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
         }
 
         Random rndGen(123456);
-        overallBb = overallBb.extend(rndGen, 1e-4);
-        overallBb.min() -= point::uniform(ROOTVSMALL);
-        overallBb.max() += point::uniform(ROOTVSMALL);
+        overallBb.inflate(rndGen, 1e-4, ROOTVSMALL);
 
 
         const indexedOctree<treeDataFace> boundaryTree
diff --git a/src/sampling/sampledSet/cellCentre/cellCentreSet.C b/src/sampling/sampledSet/cellCentre/cellCentreSet.C
index 2f523ef323..45b60315c6 100644
--- a/src/sampling/sampledSet/cellCentre/cellCentreSet.C
+++ b/src/sampling/sampledSet/cellCentre/cellCentreSet.C
@@ -133,7 +133,7 @@ Foam::cellCentreSet::cellCentreSet
         searchEngine,
         dict.getOrDefault<word>("axis", "xyz")
     ),
-    bounds_(dict.getOrDefault("bounds", boundBox::invertedBox))
+    bounds_(dict.getOrDefault("bounds", boundBox::null()))
 {
     genSamples();
 }
diff --git a/src/sampling/sampledSet/cellCentre/cellCentreSet.H b/src/sampling/sampledSet/cellCentre/cellCentreSet.H
index 83ae7c9423..27efa412ad 100644
--- a/src/sampling/sampledSet/cellCentre/cellCentreSet.H
+++ b/src/sampling/sampledSet/cellCentre/cellCentreSet.H
@@ -88,7 +88,7 @@ public:
             const polyMesh& mesh,
             const meshSearch& searchEngine,
             const word& axis,
-            const boundBox& bbox = boundBox::invertedBox
+            const boundBox& bbox = boundBox::null()
         );
 
         //- Construct from dictionary
diff --git a/src/sampling/sampledSet/patchCloud/patchCloudSet.C b/src/sampling/sampledSet/patchCloud/patchCloudSet.C
index 4358b1af1b..b1701cf28c 100644
--- a/src/sampling/sampledSet/patchCloud/patchCloudSet.C
+++ b/src/sampling/sampledSet/patchCloud/patchCloudSet.C
@@ -71,7 +71,7 @@ void Foam::patchCloudSet::calcSamples
 
     labelList patchFaces(sz);
     sz = 0;
-    treeBoundBox bb(boundBox::invertedBox);
+    treeBoundBox bb;
     for (const label patchi : patchSet_)
     {
         const polyPatch& pp = mesh().boundaryMesh()[patchi];
@@ -88,11 +88,8 @@ void Foam::patchCloudSet::calcSamples
     // Not very random
     Random rndGen(123456);
     // Make bb asymmetric just to avoid problems on symmetric meshes
-    bb = bb.extend(rndGen, 1e-4);
-
     // Make sure bb is 3D.
-    bb.min() -= point::uniform(ROOTVSMALL);
-    bb.max() += point::uniform(ROOTVSMALL);
+    bb.inflate(rndGen, 1e-4, ROOTVSMALL);
 
 
     indexedOctree<treeDataFace> patchTree
diff --git a/src/sampling/sampledSet/patchSeed/patchSeedSet.C b/src/sampling/sampledSet/patchSeed/patchSeedSet.C
index fcc633fb58..f600dc8cf2 100644
--- a/src/sampling/sampledSet/patchSeed/patchSeedSet.C
+++ b/src/sampling/sampledSet/patchSeed/patchSeedSet.C
@@ -108,14 +108,9 @@ void Foam::patchSeedSet::calcSamples
 
             treeBoundBox patchBb
             (
-                treeBoundBox(pp.points(), pp.meshPoints()).extend
-                (
-                    rndGen,
-                    1e-4
-                )
+                treeBoundBox(pp.points(), pp.meshPoints())
+                    .extend(rndGen, 1e-4, ROOTVSMALL)
             );
-            patchBb.min() -= point::uniform(ROOTVSMALL);
-            patchBb.max() += point::uniform(ROOTVSMALL);
 
             indexedOctree<treeDataFace> boundaryTree
             (
@@ -135,15 +130,7 @@ void Foam::patchSeedSet::calcSamples
             // to be found
             const scalar globalDistSqr
             (
-                //magSqr
-                //(
-                //    boundBox
-                //    (
-                //        pp.points(),
-                //        pp.meshPoints(),
-                //        true
-                //    ).span()
-                //)
+                //boundBox(pp.points(), pp.meshPoints(), true).magSqr()
                 GREAT
             );
 
@@ -168,7 +155,7 @@ void Foam::patchSeedSet::calcSamples
                 {
                     point fc(pp[nearInfo.index()].centre(pp.points()));
                     nearInfo.setPoint(fc);
-                    nearest[sampleI].second().first() = magSqr(fc-sample);
+                    nearest[sampleI].second().first() = sample.magSqr(fc);
                     nearest[sampleI].second().second() =
                         Pstream::myProcNo();
                 }
diff --git a/src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurface.C b/src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurface.C
index 14b5151da1..74ab0bc4f9 100644
--- a/src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurface.C
+++ b/src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurface.C
@@ -94,7 +94,7 @@ Foam::sampledCuttingSurface::sampledCuttingSurface
     sampledSurface(defaultSurfaceName, mesh, dict),
     cuttingSurface(defaultSurfaceName, mesh, dict),
     zoneNames_(),
-    bounds_(dict.getOrDefault("bounds", boundBox::invertedBox)),
+    bounds_(dict.getOrDefault("bounds", boundBox::null())),
     triangulate_(dict.getOrDefault("triangulate", true)),
     needsUpdate_(true)
 {
diff --git a/src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurface.H b/src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurface.H
index 20776dfb94..4815dfa5d3 100644
--- a/src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurface.H
+++ b/src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurface.H
@@ -125,7 +125,7 @@ public:
             const word& surfaceType,
             const word& surfaceName,
             const bool triangulate = true,
-            const boundBox& bounds = boundBox::invertedBox
+            const boundBox& bounds = boundBox::null()
         );
 
         //- Construct from dictionary
diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C
index 0de30644c3..238299e02c 100644
--- a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C
+++ b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C
@@ -150,7 +150,7 @@ Foam::sampledPlane::sampledPlane
     sampledSurface(name, mesh, dict),
     cuttingPlane(definePlane(mesh, dict)),
     zoneNames_(),
-    bounds_(dict.getOrDefault("bounds", boundBox::invertedBox)),
+    bounds_(dict.getOrDefault("bounds", boundBox::null())),
     triangulate_(dict.getOrDefault("triangulate", true)),
     needsUpdate_(true)
 {
diff --git a/src/sampling/surface/isoSurface/isoSurfaceParams.C b/src/sampling/surface/isoSurface/isoSurfaceParams.C
index 4a79d01c82..ac3aca3bd6 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceParams.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceParams.C
@@ -141,7 +141,7 @@ Foam::isoSurfaceParams::isoSurfaceParams
     filter_(filter),
     snap_(true),
     mergeTol_(1e-6),
-    clipBounds_(boundBox::invertedBox)
+    clipBounds_()
 {}
 
 
diff --git a/src/sampling/surface/isoSurface/isoSurfacePoint.C b/src/sampling/surface/isoSurface/isoSurfacePoint.C
index 348a2bce34..0eaa1c1e71 100644
--- a/src/sampling/surface/isoSurface/isoSurfacePoint.C
+++ b/src/sampling/surface/isoSurface/isoSurfacePoint.C
@@ -1115,7 +1115,7 @@ void Foam::isoSurfacePoint::trimToBox
     }
 
     // Generate inwards pointing planes
-    PtrList<plane> planes(treeBoundBox::faceNormals.size());
+    PtrList<plane> planes(boundBox::nFaces());
     forAll(treeBoundBox::faceNormals, faceI)
     {
         const vector& n = treeBoundBox::faceNormals[faceI];
diff --git a/src/surfMesh/triSurface/triSurfaceIO.C b/src/surfMesh/triSurface/triSurfaceIO.C
index e4da39fc5b..7a58c7bf7a 100644
--- a/src/surfMesh/triSurface/triSurfaceIO.C
+++ b/src/surfMesh/triSurface/triSurfaceIO.C
@@ -359,7 +359,7 @@ void Foam::triSurface::writeStats(Ostream& os) const
 
     bitSet pointIsUsed(points().size());
 
-    boundBox bb(boundBox::invertedBox);
+    boundBox bb;
     labelHashSet regionsUsed;
 
     for (const auto& f : *this)
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceReflecting/faceReflecting.C b/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceReflecting/faceReflecting.C
index 49beeee788..0980bc0b72 100644
--- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceReflecting/faceReflecting.C
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceReflecting/faceReflecting.C
@@ -197,10 +197,7 @@ void Foam::faceReflecting::initialise(const dictionary& coeffs)
     List<treeBoundBox> meshBb
     (
         1,
-        treeBoundBox
-        (
-            boundBox(mesh_.points(), false)
-        ).extend(rndGen, 1e-3)
+        treeBoundBox(mesh_.points()).extend(rndGen, 1e-3)
     );
 
     // Dummy bounds dictionary
@@ -332,9 +329,8 @@ void Foam::faceReflecting::calculate()
     Pstream::listCombineReduce(refDisDirsIndex, maxEqOp<label>());
     Pstream::mapCombineReduce(refFacesDirIndex, minEqOp<label>());
 
-    scalar maxBounding = 5.0*mag(mesh_.bounds().max() - mesh_.bounds().min());
-
-    reduce(maxBounding, maxOp<scalar>());
+    const scalar maxBounding =
+        returnReduce(5.0*mesh_.bounds().mag(), maxOp<scalar>());
 
     // Shoot Rays
     // From faces t = 0, r = 0 and a > 0 to all 'used' discrete reflected
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.C b/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.C
index 6395e426dc..ee2cf3a330 100644
--- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.C
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.C
@@ -194,10 +194,7 @@ void Foam::faceShading::calculate()
     List<treeBoundBox> meshBb
     (
         1,
-        treeBoundBox
-        (
-            boundBox(mesh_.points(), false)
-        ).extend(rndGen, 1e-3)
+        treeBoundBox(mesh_.points()).extend(rndGen, 1e-3)
     );
 
     // Dummy bounds dictionary
@@ -263,9 +260,8 @@ void Foam::faceShading::calculate()
         surfacesMesh.searchableSurface::write();
     }
 
-    scalar maxBounding = 5.0*mag(mesh_.bounds().max() - mesh_.bounds().min());
-
-    reduce(maxBounding, maxOp<scalar>());
+    const scalar maxBounding =
+        returnReduce(5.0*mesh_.bounds().mag(), maxOp<scalar>());
 
     // Calculate index of faces which have a direct hit (local)
     DynamicList<label> rayStartFace(nFaces + 0.01*nFaces);