diff --git a/applications/test/boundBox/Test-boundBox.C b/applications/test/boundBox/Test-boundBox.C
index 8d99361a84..9c06ccf711 100644
--- a/applications/test/boundBox/Test-boundBox.C
+++ b/applications/test/boundBox/Test-boundBox.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,8 +40,8 @@ boundBox cube(scalar start, scalar width)
 {
     return boundBox
     (
-        point(start, start, start),
-        point(start + width, start + width, start + width)
+        point::uniform(start),
+        point::uniform(start + width)
     );
 }
 
@@ -59,10 +59,19 @@ int main(int argc, char *argv[])
 
     Info<<"boundBox faces: " << boundBox::faces << endl;
     Info<<"hex faces: " << hex.modelFaces() << endl;
+    Info<<"tree-bb faces: " << treeBoundBox::faces << endl;
+    Info<<"tree-bb edges: " << treeBoundBox::edges << endl;
 
     boundBox bb = boundBox::greatBox;
     Info<<"great box: " << bb << endl;
 
+    // bb.clear();
+    // Info<<"zero box: " << bb << endl;
+
+    bb = boundBox::invertedBox;
+    Info<<"invalid box: " << bb << endl;
+    Info<< nl << endl;
+
     if (Pstream::parRun())
     {
         bb = cube(Pstream::myProcNo(), 1.1);
@@ -71,6 +80,48 @@ int main(int argc, char *argv[])
         bb.reduce();
         Pout<<"reduced: " << bb << endl;
     }
+    else
+    {
+        bb = cube(0, 1);
+        Info<<"starting box: " << bb << endl;
+
+        point pt(Zero);
+        bb.add(pt);
+        Info<<"enclose point " << pt << " -> " << bb << endl;
+
+        pt = point(0,1.5,0.5);
+        bb.add(pt);
+        Info<<"enclose point " << pt << " -> " << bb << endl;
+
+        pt = point(5,2,-2);
+        bb.add(pt);
+        Info<<"enclose point " << pt << " -> " << bb << endl;
+
+        // restart with same points
+        bb = boundBox::invertedBox;
+        bb.add(point(1,1,1));
+        bb.add(point::zero);
+        bb.add(point(0,1.5,0.5));
+        bb.add(point(5,2,-2));
+
+        Info<<"repeated " << bb << endl;
+
+        boundBox box1 = cube(0, 1);
+        boundBox box2 = cube(0, 0.75);
+        boundBox box3 = cube(0.5, 1);
+        boundBox box4 = cube(-1, 0.5);
+
+        Info<<"union of " << box1 << " and " << box2 << " => ";
+
+        box1.add(box2);
+        Info<< box1 << endl;
+
+        box1.add(box3);
+        Info<<"union with " << box3 << " => " << box1 << endl;
+
+        box1.add(box4);
+        Info<<"union with " << box4 << " => " << box1 << endl;
+    }
 
     return 0;
 }
diff --git a/applications/test/wallDist/Test-wallDist.C b/applications/test/wallDist/Test-wallDist.C
index 7a2d4f5cdf..d25bffc006 100644
--- a/applications/test/wallDist/Test-wallDist.C
+++ b/applications/test/wallDist/Test-wallDist.C
@@ -66,14 +66,14 @@ int main(int argc, char *argv[])
 
     pointField newPoints(mesh.points());
 
-    const point half(0.5*(meshBb.min() + meshBb.max()));
+    const point half = meshBb.midpoint();
 
     forAll(newPoints, pointi)
     {
         point& pt = newPoints[pointi];
 
         // expand around half
-        pt.y() +=  pt.y() - half.y();
+        pt.y() += pt.y() - half.y();
     }
 
     mesh.movePoints(newPoints);
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index 2974598fa4..53c2324aed 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -719,17 +719,13 @@ void Foam::backgroundMeshDecomposition::buildPatchAndTree()
     Pstream::gatherList(allBackgroundMeshBounds_);
     Pstream::scatterList(allBackgroundMeshBounds_);
 
-    point bbMin(GREAT, GREAT, GREAT);
-    point bbMax(-GREAT, -GREAT, -GREAT);
-
+    // find global bounding box
+    globalBackgroundBounds_ = treeBoundBox::invertedBox;
     forAll(allBackgroundMeshBounds_, proci)
     {
-        bbMin = min(bbMin, allBackgroundMeshBounds_[proci].min());
-        bbMax = max(bbMax, allBackgroundMeshBounds_[proci].max());
+        globalBackgroundBounds_.add(allBackgroundMeshBounds_[proci]);
     }
 
-    globalBackgroundBounds_ = treeBoundBox(bbMin, bbMax);
-
     if (false)
     {
         OFstream fStr
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index 2b58f79335..cad195d013 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -558,18 +558,11 @@ Foam::label Foam::checkTopology
 
                 if (allGeometry)
                 {
-                    const pointField& pts = pp.points();
                     const labelList& mp = pp.meshPoints();
 
                     if (returnReduce(mp.size(), sumOp<label>()) > 0)
                     {
-                        boundBox bb(point::max, point::min);
-                        forAll(mp, i)
-                        {
-                            bb.min() = min(bb.min(), pts[mp[i]]);
-                            bb.max() = max(bb.max(), pts[mp[i]]);
-                        }
-                        bb.reduce();
+                        boundBox bb(pp.points(), mp, true); // reduce
                         Info<< ' ' << bb;
                     }
                 }
diff --git a/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C b/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
index 251b6523f4..ef68688bf7 100644
--- a/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
+++ b/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -336,10 +336,7 @@ boundBox procBounds
             )
         );
 
-        boundBox domainBb(points, false);
-
-        bb.min() = min(bb.min(), domainBb.min());
-        bb.max() = max(bb.max(), domainBb.max());
+        bb.add(points);
     }
 
     return bb;
diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.C b/src/OpenFOAM/meshes/boundBox/boundBox.C
index 982fe6bf53..1a4baf4af9 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBox.C
+++ b/src/OpenFOAM/meshes/boundBox/boundBox.C
@@ -49,33 +49,23 @@ const Foam::faceList Foam::boundBox::faces
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::boundBox::calculate(const UList<point>& points, const bool doReduce)
+void Foam::boundBox::calculate(const UList<point>& points, bool doReduce)
 {
     if (points.empty())
     {
-        min_ = Zero;
-        max_ = Zero;
-
         if (doReduce && Pstream::parRun())
         {
-            // Use values that get overwritten by reduce minOp, maxOp below
-            min_ = point(VGREAT, VGREAT, VGREAT);
-            max_ = point(-VGREAT, -VGREAT, -VGREAT);
+            // Values that get overwritten by subsequent reduce operation
+            operator=(invertedBox);
         }
     }
     else
     {
-        min_ = points[0];
-        max_ = points[0];
-
-        for (label i = 1; i < points.size(); i++)
-        {
-            min_ = ::Foam::min(min_, points[i]);
-            max_ = ::Foam::max(max_, points[i]);
-        }
+        operator=(invertedBox);
+        add(points);
     }
 
-    // Reduce parallel information
+    // Parallel reduction
     if (doReduce)
     {
         reduce();
@@ -85,7 +75,7 @@ void Foam::boundBox::calculate(const UList<point>& points, const bool doReduce)
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::boundBox::boundBox(const UList<point>& points, const bool doReduce)
+Foam::boundBox::boundBox(const UList<point>& points, bool doReduce)
 :
     min_(Zero),
     max_(Zero)
@@ -94,13 +84,13 @@ Foam::boundBox::boundBox(const UList<point>& points, const bool doReduce)
 }
 
 
-Foam::boundBox::boundBox(const tmp<pointField>& points, const bool doReduce)
+Foam::boundBox::boundBox(const tmp<pointField>& tpoints, bool doReduce)
 :
     min_(Zero),
     max_(Zero)
 {
-    calculate(points(), doReduce);
-    points.clear();
+    calculate(tpoints(), doReduce);
+    tpoints.clear();
 }
 
 
@@ -108,7 +98,7 @@ Foam::boundBox::boundBox
 (
     const UList<point>& points,
     const labelUList& indices,
-    const bool doReduce
+    bool doReduce
 )
 :
     min_(Zero),
@@ -118,24 +108,17 @@ Foam::boundBox::boundBox
     {
         if (doReduce && Pstream::parRun())
         {
-            // Use values that get overwritten by reduce minOp, maxOp below
-            min_ = point(VGREAT, VGREAT, VGREAT);
-            max_ = point(-VGREAT, -VGREAT, -VGREAT);
+            // Values that get overwritten by subsequent reduce operation
+            operator=(invertedBox);
         }
     }
     else
     {
-        min_ = points[indices[0]];
-        max_ = points[indices[0]];
-
-        for (label i=1; i < indices.size(); ++i)
-        {
-            min_ = ::Foam::min(min_, points[indices[i]]);
-            max_ = ::Foam::max(max_, points[indices[i]]);
-        }
+        operator=(invertedBox);
+        add(points, indices);
     }
 
-    // Reduce parallel information
+    // Parallel reduction
     if (doReduce)
     {
         reduce();
@@ -147,8 +130,8 @@ Foam::boundBox::boundBox
 
 Foam::tmp<Foam::pointField> Foam::boundBox::points() const
 {
-    tmp<pointField> tPts = tmp<pointField>(new pointField(8));
-    pointField& pt = tPts.ref();
+    tmp<pointField> tpoints = tmp<pointField>(new pointField(8));
+    pointField& pt = tpoints.ref();
 
     pt[0] = min_;                                   // min-x, min-y, min-z
     pt[1] = point(max_.x(), min_.y(), min_.z());    // max-x, min-y, min-z
@@ -159,13 +142,13 @@ Foam::tmp<Foam::pointField> Foam::boundBox::points() const
     pt[6] = max_;                                   // max-x, max-y, max-z
     pt[7] = point(min_.x(), max_.y(), max_.z());    // min-x, max-y, max-z
 
-    return tPts;
+    return tpoints;
 }
 
 
 void Foam::boundBox::inflate(const scalar s)
 {
-    vector ext = vector::one*s*mag();
+    const vector ext = vector::one*s*mag();
 
     min_ -= ext;
     max_ += ext;
@@ -179,6 +162,15 @@ void Foam::boundBox::reduce()
 }
 
 
+bool Foam::boundBox::intersect(const boundBox& bb)
+{
+    min_ = ::Foam::max(min_, bb.min_);
+    max_ = ::Foam::min(max_, bb.max_);
+
+    return !empty();
+}
+
+
 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 51ae4f52f2..a7e6e569f4 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBox.H
+++ b/src/OpenFOAM/meshes/boundBox/boundBox.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -60,22 +60,19 @@ class boundBox
 {
     // Private data
 
-        //- Minimum and maximum describing the bounding box
+        //- Minimum and maximum points describing the bounding box
         point min_, max_;
 
     // Private Member Functions
 
         //- Calculate the bounding box from the given points.
         //  Does parallel communication (doReduce = true)
-        void calculate(const UList<point>&, const bool doReduce = true);
+        void calculate(const UList<point>& points, bool doReduce = true);
 
 public:
 
     // Static data members
 
-        //- The great value used for greatBox and invertedBox
-        static const scalar great;
-
         //- A very large boundBox: min/max == -/+ VGREAT
         static const boundBox greatBox;
 
@@ -91,46 +88,52 @@ public:
         //- Construct null, setting points to zero
         inline boundBox();
 
+        //- Construct a bounding box containing a single initial point
+        explicit inline boundBox(const point& pt);
+
         //- Construct from components
         inline boundBox(const point& min, const point& max);
 
         //- Construct as the bounding box of the given points
         //  Does parallel communication (doReduce = true)
-        boundBox(const UList<point>&, const bool doReduce = true);
+        explicit boundBox(const UList<point>& points, bool doReduce = true);
 
         //- Construct as the bounding box of the given temporary pointField.
         //  Does parallel communication (doReduce = true)
-        boundBox(const tmp<pointField>&, const bool doReduce = true);
+        explicit boundBox(const tmp<pointField>& tpoints, bool doReduce = true);
 
-        //- Construct bounding box as subset of the pointField.
+        //- Construct bounding box as an indirect subset of the points.
         //  The indices could be from cell/face etc.
         //  Does parallel communication (doReduce = true)
         boundBox
         (
-            const UList<point>&,
+            const UList<point>& points,
             const labelUList& indices,
-            const bool doReduce = true
+            bool doReduce = true
         );
 
-        //- Construct bounding box as subset of the pointField.
+        //- Construct bounding box as an indirect subset of the points.
         //  The indices could be from edge/triFace etc.
         //  Does parallel communication (doReduce = true)
         template<unsigned Size>
         boundBox
         (
-            const UList<point>&,
+            const UList<point>& points,
             const FixedList<label, Size>& indices,
-            const bool doReduce = true
+            bool doReduce = true
         );
 
         //- Construct from Istream
-        inline boundBox(Istream&);
+        inline boundBox(Istream& is);
 
 
     // Member functions
 
         // Access
 
+            //- Bounding box is inverted, contains no points.
+            inline bool empty() const;
+
             //- Minimum describing the bounding box
             inline const point& min() const;
 
@@ -170,83 +173,130 @@ public:
 
         // Manipulate
 
+            //- Extend to include the second box.
+            inline void add(const boundBox& bb);
+
+            //- Extend to include the point.
+            inline void add(const point& pt);
+
+            //- Extend to include the points.
+            inline void add(const UList<point>& points);
+
+            //- Extend to include the points from the temporary point field.
+            inline void add(const tmp<pointField>& tpoints);
+
+            //- Extend to include the subset of the point field.
+            //  The indices could be from cell/face etc.
+            inline void add
+            (
+                const UList<point>& points,
+                const labelUList& indices
+            );
+
+            //- Extend to include the points.
+            template<unsigned Size>
+            void add(const FixedList<point, Size>& points);
+
+            //- Extend to include the subset of the point field.
+            //  The indices could be from edge/triFace etc.
+            template<unsigned Size>
+            void add
+            (
+                const UList<point>& points,
+                const FixedList<label, Size>& indices
+            );
+
             //- Inflate box by factor*mag(span) in all dimensions
             void inflate(const scalar s);
 
             //- Parallel reduction of min/max values
             void reduce();
 
+            //- Intersection (union) with the second box.
+            //  The return value is true if the intersection is non-empty.
+            bool intersect(const boundBox& bb);
+
 
         // Query
 
             //- Overlaps/touches boundingBox?
-            inline bool overlaps(const boundBox&) const;
+            inline bool overlaps(const boundBox& bb) const;
 
             //- Overlaps boundingSphere (centre + sqr(radius))?
-            inline bool overlaps(const point&, const scalar radiusSqr) const;
+            inline bool overlaps
+            (
+                const point& centre,
+                const scalar radiusSqr
+            ) const;
 
             //- Contains point? (inside or on edge)
-            inline bool contains(const point&) const;
+            inline bool contains(const point& pt) const;
 
             //- Fully contains other boundingBox?
-            inline bool contains(const boundBox&) const;
+            inline bool contains(const boundBox& bb) const;
 
             //- Contains point? (inside only)
-            inline bool containsInside(const point&) const;
+            inline bool containsInside(const point& pt) const;
 
             //- Contains all of the points? (inside or on edge)
-            bool contains(const UList<point>&) const;
+            bool contains(const UList<point>& points) const;
 
-            //- Contains all of the points? (inside or on edge)
+            //- Contains all of the subset of points? (inside or on edge)
             bool contains
             (
-                const UList<point>&,
+                const UList<point>& points,
                 const labelUList& indices
             ) const;
 
-            //- Contains all of the points? (inside or on edge)
+            //- Contains all of the subset of points? (inside or on edge)
             template<unsigned Size>
             bool contains
             (
-                const UList<point>&,
+                const UList<point>& points,
                 const FixedList<label, Size>& indices
             ) const;
 
 
             //- Contains any of the points? (inside or on edge)
-            bool containsAny(const UList<point>&) const;
+            bool containsAny(const UList<point>& points) const;
 
-            //- Contains any of the points? (inside or on edge)
+            //- Contains any of the subset of points? (inside or on edge)
             bool containsAny
             (
-                const UList<point>&,
+                const UList<point>& points,
                 const labelUList& indices
             ) const;
 
-            //- Contains any of the points? (inside or on edge)
+            //- Contains any of the subset of points? (inside or on edge)
             template<unsigned Size>
             bool containsAny
             (
-                const UList<point>&,
+                const UList<point>& points,
                 const FixedList<label, Size>& indices
             ) const;
 
             //- Return the nearest point on the boundBox to the supplied point.
             //  If point is inside the boundBox then the point is returned
             //  unchanged.
-            point nearest(const point&) const;
+            point nearest(const point& pt) const;
+
+
+    // Member Operators
+
+        //- Extend box to include the second box, as per the add() method.
+        inline void operator+=(const boundBox& bb);
 
 
     // Friend Operators
 
-        inline friend bool operator==(const boundBox&, const boundBox&);
-        inline friend bool operator!=(const boundBox&, const boundBox&);
+        inline friend bool operator==(const boundBox& a, const boundBox& b);
+        inline friend bool operator!=(const boundBox& a, const boundBox& b);
 
 
     // IOstream operator
 
-        friend Istream& operator>>(Istream&, boundBox&);
-        friend Ostream& operator<<(Ostream&, const boundBox&);
+        friend Istream& operator>>(Istream& is, boundBox& bb);
+        friend Ostream& operator<<(Ostream& os, const boundBox& bb);
 };
 
 
diff --git a/src/OpenFOAM/meshes/boundBox/boundBoxI.H b/src/OpenFOAM/meshes/boundBox/boundBoxI.H
index 13af23e99c..611130477a 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBoxI.H
+++ b/src/OpenFOAM/meshes/boundBox/boundBoxI.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "boundBox.H"
-#include "pointField.H"
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -36,6 +35,13 @@ inline Foam::boundBox::boundBox()
 {}
 
 
+inline Foam::boundBox::boundBox(const point& pt)
+:
+    min_(pt),
+    max_(pt)
+{}
+
+
 inline Foam::boundBox::boundBox(const point& min, const point& max)
 :
     min_(min),
@@ -51,6 +57,12 @@ inline Foam::boundBox::boundBox(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline bool Foam::boundBox::empty() const
+{
+    return (min_ > max_);
+}
+
+
 inline const Foam::point& Foam::boundBox::min() const
 {
     return min_;
@@ -117,6 +129,52 @@ inline Foam::scalar Foam::boundBox::avgDim() const
 }
 
 
+inline void Foam::boundBox::add(const boundBox& bb)
+{
+    min_ = ::Foam::min(min_, bb.min_);
+    max_ = ::Foam::max(max_, bb.max_);
+}
+
+
+inline void Foam::boundBox::add(const point& pt)
+{
+    min_ = ::Foam::min(min_, pt);
+    max_ = ::Foam::max(max_, pt);
+}
+
+
+inline void Foam::boundBox::add(const UList<point>& points)
+{
+    forAll(points, i)
+    {
+        add(points[i]);
+    }
+}
+
+
+inline void Foam::boundBox::add(const tmp<pointField>& tpoints)
+{
+    add(tpoints());
+    tpoints.clear();
+}
+
+
+inline void Foam::boundBox::add
+(
+    const UList<point>& points,
+    const labelUList& indices
+)
+{
+    if (!points.empty())
+    {
+        forAll(indices, i)
+        {
+            add(points[indices[i]]);
+        }
+    }
+}
+
+
 inline bool Foam::boundBox::overlaps(const boundBox& bb) const
 {
     return
@@ -138,7 +196,7 @@ inline bool Foam::boundBox::overlaps
     // Find nearest point on bb.
     scalar distSqr = 0;
 
-    for (direction dir = 0; dir < vector::nComponents; dir++)
+    for (direction dir = 0; dir < vector::nComponents; ++dir)
     {
         scalar d0 = min_[dir] - centre[dir];
         scalar d1 = max_[dir] - centre[dir];
@@ -178,7 +236,6 @@ inline bool Foam::boundBox::contains(const point& pt) const
 }
 
 
-// this.bb fully contains bb
 inline bool Foam::boundBox::contains(const boundBox& bb) const
 {
     return contains(bb.min()) && contains(bb.max());
@@ -196,6 +253,14 @@ inline bool Foam::boundBox::containsInside(const point& pt) const
 }
 
 
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline void Foam::boundBox::operator+=(const boundBox& bb)
+{
+    add(bb);
+}
+
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::operator==(const boundBox& a, const boundBox& b)
diff --git a/src/OpenFOAM/meshes/boundBox/boundBoxTemplates.C b/src/OpenFOAM/meshes/boundBox/boundBoxTemplates.C
index 1ce09cb3fa..0bb02715a8 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBoxTemplates.C
+++ b/src/OpenFOAM/meshes/boundBox/boundBoxTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,7 +26,6 @@ License
 #include "boundBox.H"
 #include "FixedList.H"
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<unsigned Size>
@@ -34,7 +33,7 @@ Foam::boundBox::boundBox
 (
     const UList<point>& points,
     const FixedList<label, Size>& indices,
-    const bool doReduce
+    bool doReduce
 )
 :
     min_(Zero),
@@ -45,24 +44,17 @@ Foam::boundBox::boundBox
     {
         if (doReduce && Pstream::parRun())
         {
-            // Use values that get overwritten by reduce minOp, maxOp below
-            min_ = point(VGREAT, VGREAT, VGREAT);
-            max_ = point(-VGREAT, -VGREAT, -VGREAT);
+            // Values that get overwritten by subsequent reduce operation
+            operator=(invertedBox);
         }
     }
     else
     {
-        min_ = points[indices[0]];
-        max_ = points[indices[0]];
-
-        for (unsigned i=1; i < Size; ++i)
-        {
-            min_ = ::Foam::min(min_, points[indices[i]]);
-            max_ = ::Foam::max(max_, points[indices[i]]);
-        }
+        operator=(invertedBox);
+        add(points, indices);
     }
 
-    // Reduce parallel information
+    // Parallel reduction
     if (doReduce)
     {
         reduce();
@@ -72,6 +64,38 @@ Foam::boundBox::boundBox
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+template<unsigned Size>
+void Foam::boundBox::add
+(
+    const FixedList<point, Size>& points
+)
+{
+    // a FixedList is never empty
+    for (unsigned i=0; i < Size; ++i)
+    {
+        add(points[i]);
+    }
+}
+
+
+template<unsigned Size>
+void Foam::boundBox::add
+(
+    const UList<point>& points,
+    const FixedList<label, Size>& indices
+)
+{
+    // points may be empty, but a FixedList is never empty
+    if (!points.empty())
+    {
+        for (unsigned i=0; i < Size; ++i)
+        {
+            add(points[indices[i]]);
+        }
+    }
+}
+
+
 template<unsigned Size>
 bool Foam::boundBox::contains
 (
@@ -79,7 +103,7 @@ bool Foam::boundBox::contains
     const FixedList<label, Size>& indices
 ) const
 {
-    // a FixedList is never empty
+    // points may be empty, but a FixedList is never empty
     if (points.empty())
     {
         return false;
@@ -104,7 +128,7 @@ bool Foam::boundBox::containsAny
     const FixedList<label, Size>& indices
 ) const
 {
-    // a FixedList is never empty
+    // points may be empty, but a FixedList is never empty
     if (points.empty())
     {
         return false;
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSearch.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSearch.C
index d621df5ca3..66fde47413 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSearch.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSearch.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -236,9 +236,8 @@ void Foam::PatchTools::calcBounds
             label pointi = f[fp];
             if (pointIsUsed.set(pointi, 1u))
             {
-                bb.min() = ::Foam::min(bb.min(), points[pointi]);
-                bb.max() = ::Foam::max(bb.max(), points[pointi]);
-                nPoints++;
+                bb.add(points[pointi]);
+                ++nPoints;
             }
         }
     }
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C
index 31ce6616fd..c26745d96f 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,8 +48,7 @@ bool Foam::primitiveMesh::pointInCellBB
 
     if (inflationFraction > SMALL)
     {
-        vector inflation = inflationFraction*vector::one*mag(bb.span());
-        bb = boundBox(bb.min() - inflation, bb.max() + inflation);
+        bb.inflate(inflationFraction);
     }
 
     return bb.contains(p);
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H
index cb42d33d4b..8ac906ca62 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -92,11 +92,10 @@ public:
         //- Calculate the bounding box
         inline treeBoundBox bounds() const
         {
-            treeBoundBox bb(operator[](0), operator[](0));
-            for (label i = 1; i < size(); i++)
+            treeBoundBox bb(operator[](0));
+            for (label i = 1; i < size(); ++i)
             {
-                bb.min() = min(bb.min(), operator[](i));
-                bb.max() = max(bb.max(), operator[](i));
+                bb.add(operator[](i));
             }
             return bb;
         }
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triPoints.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triPoints.H
index 67c7efb902..7f1dc8d42a 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triPoints.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triPoints.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -78,11 +78,10 @@ public:
         //- Calculate the bounding box
         inline treeBoundBox bounds() const
         {
-            treeBoundBox bb(operator[](0), operator[](0));
-            for (label i = 1; i < size(); i++)
+            treeBoundBox bb(operator[](0));
+            for (label i = 1; i < size(); ++i)
             {
-                bb.min() = min(bb.min(), operator[](i));
-                bb.max() = max(bb.max(), operator[](i));
+                bb.add(operator[](i));
             }
             return bb;
         }
diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C
index cd804fcf94..00af118874 100644
--- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C
+++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,50 +24,33 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "treeBoundBox.H"
-#include "ListOps.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const Foam::scalar Foam::treeBoundBox::great(GREAT);
-
 const Foam::treeBoundBox Foam::treeBoundBox::greatBox
 (
-    vector(-GREAT, -GREAT, -GREAT),
-    vector(GREAT, GREAT, GREAT)
+    point::uniform(-GREAT),
+    point::uniform(GREAT)
 );
 
-
 const Foam::treeBoundBox Foam::treeBoundBox::invertedBox
 (
-    vector(GREAT, GREAT, GREAT),
-    vector(-GREAT, -GREAT, -GREAT)
+    point::uniform(GREAT),
+    point::uniform(-GREAT)
 );
 
-
-//! \cond ignoreDocumentation
-//- Skip documentation : local scope only
-const Foam::label facesArray[6][4] =
-{
-    {0, 4, 6, 2}, // left
-    {1, 3, 7, 5}, // right
-    {0, 1, 5, 4}, // bottom
-    {2, 6, 7, 3}, // top
-    {0, 2, 3, 1}, // back
-    {4, 5, 7, 6}  // front
-};
-//! \endcond
-
-
 const Foam::faceList Foam::treeBoundBox::faces
-(
-    initListList<face, label, 6, 4>(facesArray)
-);
+({
+    face{0, 4, 6, 2}, // left
+    face{1, 3, 7, 5}, // right
+    face{0, 1, 5, 4}, // bottom
+    face{2, 6, 7, 3}, // top
+    face{0, 2, 3, 1}, // back
+    face{4, 5, 7, 6}  // front
+});
 
-
-//! \cond  ignoreDocumentation
-//- Skip documentation : local scope only
-const Foam::label edgesArray[12][2] =
-{
+const Foam::edgeList Foam::treeBoundBox::edges
+({
     {0, 1}, // 0
     {1, 3},
     {2, 3}, // 2
@@ -80,48 +63,17 @@ const Foam::label edgesArray[12][2] =
     {1, 5},
     {3, 7}, // 10
     {2, 6}
-};
-//! \endcond
-
-
-const Foam::edgeList Foam::treeBoundBox::edges
-(
-    //initListList<edge, label, 12, 2>(edgesArray)
-    calcEdges(edgesArray)
-);
-
+});
 
 const Foam::FixedList<Foam::vector, 6> Foam::treeBoundBox::faceNormals
-(
-    calcFaceNormals()
-);
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::edgeList Foam::treeBoundBox::calcEdges(const label edgesArray[12][2])
-{
-    edgeList edges(12);
-    forAll(edges, edgeI)
-    {
-        edges[edgeI][0] = edgesArray[edgeI][0];
-        edges[edgeI][1] = edgesArray[edgeI][1];
-    }
-    return edges;
-}
-
-
-Foam::FixedList<Foam::vector, 6> Foam::treeBoundBox::calcFaceNormals()
-{
-    FixedList<vector, 6> normals;
-    normals[LEFT]   = vector(-1,  0,  0);
-    normals[RIGHT]  = vector( 1,  0,  0);
-    normals[BOTTOM] = vector( 0, -1,  0);
-    normals[TOP]    = vector( 0,  1,  0);
-    normals[BACK]   = vector( 0,  0, -1);
-    normals[FRONT]  = vector( 0,  0,  1);
-    return normals;
-}
+({
+    vector(-1,  0,  0), // left
+    vector( 1,  0,  0), // right
+    vector( 0, -1,  0), // bottom
+    vector( 0,  1,  0), // top
+    vector( 0,  0, -1), // back
+    vector( 0,  0,  1)  // front
+});
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -164,15 +116,15 @@ Foam::treeBoundBox::treeBoundBox
 
 Foam::tmp<Foam::pointField> Foam::treeBoundBox::points() const
 {
-    tmp<pointField> tPts = tmp<pointField>(new pointField(8));
-    pointField& points = tPts.ref();
+    tmp<pointField> tpoints = tmp<pointField>(new pointField(8));
+    pointField& pts = tpoints.ref();
 
-    forAll(points, octant)
+    forAll(pts, octant)
     {
-        points[octant] = corner(octant);
+        pts[octant] = corner(octant);
     }
 
-    return tPts;
+    return tpoints;
 }
 
 
@@ -394,7 +346,7 @@ bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const
 {
     // Compare all components against min and max of bb
 
-    for (direction cmpt=0; cmpt<3; cmpt++)
+    for (direction cmpt=0; cmpt<3; ++cmpt)
     {
         if (pt[cmpt] < min()[cmpt])
         {
diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H
index 7632b15fa1..4a18d906ab 100644
--- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H
+++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,7 +46,6 @@ Description
 
     For the front plane add 4 to the point labels.
 
-
 SourceFiles
     treeBoundBoxI.H
     treeBoundBox.C
@@ -74,11 +73,11 @@ class Random;
 
 class treeBoundBox;
 
-bool operator==(const treeBoundBox&, const treeBoundBox&);
-bool operator!=(const treeBoundBox&, const treeBoundBox&);
+bool operator==(const treeBoundBox& a, const treeBoundBox& b);
+bool operator!=(const treeBoundBox& a, const treeBoundBox& b);
 
-Istream& operator>>(Istream& is, treeBoundBox&);
-Ostream& operator<<(Ostream& os, const treeBoundBox&);
+Istream& operator>>(Istream& is, treeBoundBox& bb);
+Ostream& operator<<(Ostream& os, const treeBoundBox& bb);
 
 
 /*---------------------------------------------------------------------------*\
@@ -90,21 +89,10 @@ class treeBoundBox
     public boundBox
 {
 
-private:
-
-        //- To initialise edges.
-        static edgeList calcEdges(const label[12][2]);
-
-        //- To initialise faceNormals.
-        static FixedList<vector, 6> calcFaceNormals();
-
 public:
 
     // Static data members
 
-        //- The great value used for greatBox and invertedBox
-        static const scalar great;
-
         //- As per boundBox::greatBox, but with GREAT instead of VGREAT
         static const treeBoundBox greatBox;
 
@@ -169,7 +157,7 @@ public:
         //- Edge to point addressing
         static const edgeList edges;
 
-        //- Per face the unit normal
+        //- The unit normal per face
         static const FixedList<vector, 6> faceNormals;
 
 
@@ -181,16 +169,19 @@ public:
         //- Construct from a boundBox
         explicit inline treeBoundBox(const boundBox& bb);
 
+        //- Construct a bounding box containing a single initial point
+        explicit inline treeBoundBox(const point& pt);
+
         //- Construct from components
         inline treeBoundBox(const point& min, const point& max);
 
         //- Construct as the bounding box of the given pointField.
         //  Local processor domain only (no reduce as in boundBox)
-        explicit treeBoundBox(const UList<point>&);
+        explicit treeBoundBox(const UList<point>& points);
 
         //- Construct as subset of points
         //  Local processor domain only (no reduce as in boundBox)
-        treeBoundBox(const UList<point>&, const labelUList& indices);
+        treeBoundBox(const UList<point>& points, const labelUList& indices);
 
         //- Construct as subset of points
         //  The indices could be from edge/triFace etc.
@@ -198,13 +189,13 @@ public:
         template<unsigned Size>
         treeBoundBox
         (
-            const UList<point>&,
+            const UList<point>& points,
             const FixedList<label, Size>& indices
         );
 
 
         //- Construct from Istream
-        inline treeBoundBox(Istream&);
+        inline treeBoundBox(Istream& is);
 
 
     // Member functions
@@ -220,13 +211,13 @@ public:
 
         // Check
 
-            //- Corner point given octant
-            inline point corner(const direction) const;
+            //- Corner point of given octant
+            inline point corner(const direction octant) const;
 
-            //- Sub box given by octant number. Midpoint calculated.
-            treeBoundBox subBbox(const direction) const;
+            //- Sub-box of given octant. Midpoint calculated.
+            treeBoundBox subBbox(const direction octant) const;
 
-            //- Sub box given by octant number. Midpoint provided.
+            //- Sub-box given by octant number. Midpoint provided.
             treeBoundBox subBbox(const point& mid, const direction) const;
 
             //- Returns octant number given point and the calculated midpoint.
@@ -318,10 +309,10 @@ public:
             bool contains(const vector& dir, const point&) const;
 
             //- Code position of point on bounding box faces
-            direction faceBits(const point&) const;
+            direction faceBits(const point& pt) const;
 
             //- Position of point relative to bounding box
-            direction posBits(const point&) const;
+            direction posBits(const point& pt) const;
 
             //- Calculate nearest and furthest (to point) vertex coords of
             //  bounding box
@@ -333,7 +324,7 @@ public:
             ) const;
 
             //- Returns distance point to furthest away corner.
-            scalar maxDist(const point&) const;
+            scalar maxDist(const point& pt) const;
 
             //- Compare distance to point with other bounding box
             //  return:
@@ -342,22 +333,24 @@ public:
             //  +1 : all vertices of my bounding box are further away than
             //       any of other
             //   0 : none of the above.
-            label distanceCmp(const point&, const treeBoundBox& other) const;
+            label distanceCmp(const point& pt, const treeBoundBox& other) const;
 
             //- Return slightly wider bounding box
             //  Extends all dimensions with s*span*Random::scalar01()
             //  and guarantees in any direction s*mag(span) minimum width
-            inline treeBoundBox extend(Random&, const scalar s) const;
+            inline treeBoundBox extend(Random& rndGen, const scalar s) const;
+
 
     // Friend Operators
 
-        friend bool operator==(const treeBoundBox&, const treeBoundBox&);
-        friend bool operator!=(const treeBoundBox&, const treeBoundBox&);
+        friend bool operator==(const treeBoundBox& a, const treeBoundBox& b);
+        friend bool operator!=(const treeBoundBox& a, const treeBoundBox& b);
+
 
     // IOstream operator
 
-        friend Istream& operator>>(Istream& is, treeBoundBox&);
-        friend Ostream& operator<<(Ostream& os, const treeBoundBox&);
+        friend Istream& operator>>(Istream& is, treeBoundBox& bb);
+        friend Ostream& operator<<(Ostream& os, const treeBoundBox& bb);
 };
 
 
diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H
index 3480eac00d..4b073fd975 100644
--- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H
+++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,18 +34,24 @@ inline Foam::treeBoundBox::treeBoundBox()
 {}
 
 
-inline Foam::treeBoundBox::treeBoundBox(const point& min, const point& max)
-:
-    boundBox(min, max)
-{}
-
-
 inline Foam::treeBoundBox::treeBoundBox(const boundBox& bb)
 :
     boundBox(bb)
 {}
 
 
+inline Foam::treeBoundBox::treeBoundBox(const point& pt)
+:
+    boundBox(pt)
+{}
+
+
+inline Foam::treeBoundBox::treeBoundBox(const point& min, const point& max)
+:
+    boundBox(min, max)
+{}
+
+
 inline Foam::treeBoundBox::treeBoundBox(Istream& is)
 :
     boundBox(is)
@@ -312,7 +318,6 @@ inline void Foam::treeBoundBox::searchOrder
 }
 
 
-//- Return slightly wider bounding box
 inline Foam::treeBoundBox Foam::treeBoundBox::extend
 (
     Random& rndGen,
@@ -326,7 +331,7 @@ inline Foam::treeBoundBox Foam::treeBoundBox::extend
     // Make 3D
     scalar minSpan = s * Foam::mag(newSpan);
 
-    for (direction dir = 0; dir < vector::nComponents; dir++)
+    for (direction dir = 0; dir < vector::nComponents; ++dir)
     {
         newSpan[dir] = Foam::max(newSpan[dir], minSpan);
     }
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
index fab34ba8cd..b7cd414123 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
+++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -125,8 +125,7 @@ void Foam::polyMeshGeometry::updateCellCentresAndVols
         {
             const point& fc = faceCentres_[cFaces[cFaceI]];
             cEst += fc;
-            bb.max() = max(bb.max(), fc);
-            bb.min() = min(bb.min(), fc);
+            bb.add(fc);
         }
         cEst /= cFaces.size();
 
diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
index c6d79d45b1..768618bf5e 100644
--- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
+++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -199,21 +199,9 @@ void Foam::shellSurfaces::orient()
 
             if (shell.triSurface::size())
             {
-                const pointField& points = shell.points();
-
                 hasSurface = true;
-
-                boundBox shellBb(points[0], points[0]);
                 // Assume surface is compact!
-                forAll(points, i)
-                {
-                    const point& pt = points[i];
-                    shellBb.min() = min(shellBb.min(), pt);
-                    shellBb.max() = max(shellBb.max(), pt);
-                }
-
-                overallBb.min() = min(overallBb.min(), shellBb.min());
-                overallBb.max() = max(overallBb.max(), shellBb.max());
+                overallBb.add(shell.points());
             }
         }
     }
diff --git a/src/meshTools/AABBTree/AABBTree.C b/src/meshTools/AABBTree/AABBTree.C
index adbe0ab5ad..19c11f33a9 100644
--- a/src/meshTools/AABBTree/AABBTree.C
+++ b/src/meshTools/AABBTree/AABBTree.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -198,10 +198,10 @@ void Foam::AABBTree<Type>::createBoxes
     // Assign the objects to min or max bin
 
     DynamicList<label> minBinObjectIDs(objectIDs.size());
-    treeBoundBox minBb(point::max, point::min);
+    treeBoundBox minBb(boundBox::invertedBox);
 
     DynamicList<label> maxBinObjectIDs(objectIDs.size());
-    treeBoundBox maxBb(point::max, point::min);
+    treeBoundBox maxBb(boundBox::invertedBox);
 
     forAll(objectIDs, i)
     {
@@ -229,16 +229,12 @@ void Foam::AABBTree<Type>::createBoxes
         if (intoMin)
         {
             minBinObjectIDs.append(objI);
-            const boundBox objBb(points, obj, false);
-            minBb.min() = min(minBb.min(), objBb.min());
-            minBb.max() = max(minBb.max(), objBb.max());
+            minBb.add(points, obj);
         }
         if (intoMax)
         {
             maxBinObjectIDs.append(objI);
-            const boundBox objBb(points, obj, false);
-            maxBb.min() = min(maxBb.min(), objBb.min());
-            maxBb.max() = max(maxBb.max(), objBb.max());
+            maxBb.add(points, obj);
         }
     }
 
diff --git a/src/meshTools/indexedOctree/treeDataFace.C b/src/meshTools/indexedOctree/treeDataFace.C
index 471740a924..73983df95d 100644
--- a/src/meshTools/indexedOctree/treeDataFace.C
+++ b/src/meshTools/indexedOctree/treeDataFace.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,20 +41,7 @@ namespace Foam
 
 Foam::treeBoundBox Foam::treeDataFace::calcBb(const label facei) const
 {
-    const pointField& points = mesh_.points();
-
-    const face& f = mesh_.faces()[facei];
-
-    treeBoundBox bb(points[f[0]], points[f[0]]);
-
-    for (label fp = 1; fp < f.size(); fp++)
-    {
-        const point& p = points[f[fp]];
-
-        bb.min() = min(bb.min(), p);
-        bb.max() = max(bb.max(), p);
-    }
-    return bb;
+    return treeBoundBox(mesh_.points(), mesh_.faces()[facei]);
 }
 
 
diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
index 00223a009e..f1c969762b 100644
--- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
+++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,26 +31,6 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class PatchType>
-Foam::treeBoundBox Foam::treeDataPrimitivePatch<PatchType>::calcBb
-(
-    const pointField& points,
-    const face& f
-)
-{
-    treeBoundBox bb(points[f[0]], points[f[0]]);
-
-    for (label fp = 1; fp < f.size(); fp++)
-    {
-        const point& p = points[f[fp]];
-
-        bb.min() = min(bb.min(), p);
-        bb.max() = max(bb.max(), p);
-    }
-    return bb;
-}
-
-
 template<class PatchType>
 void Foam::treeDataPrimitivePatch<PatchType>::update()
 {
@@ -60,7 +40,7 @@ void Foam::treeDataPrimitivePatch<PatchType>::update()
 
         forAll(patch_, i)
         {
-            bbs_[i] = calcBb(patch_.points(), patch_[i]);
+            bbs_[i] = treeBoundBox(patch_.points(), patch_[i]);
         }
     }
 }
@@ -386,19 +366,14 @@ bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
 ) const
 {
     // 1. Quick rejection: bb does not intersect face bb at all
-    if (cacheBb_)
+    if
+    (
+        cacheBb_
+      ? !cubeBb.overlaps(bbs_[index])
+      : !cubeBb.overlaps(treeBoundBox(patch_.points(), patch_[index]))
+    )
     {
-        if (!cubeBb.overlaps(bbs_[index]))
-        {
-            return false;
-        }
-    }
-    else
-    {
-        if (!cubeBb.overlaps(calcBb(patch_.points(), patch_[index])))
-        {
-            return false;
-        }
+        return false;
     }
 
 
@@ -459,19 +434,14 @@ bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
 ) const
 {
     // 1. Quick rejection: sphere does not intersect face bb at all
-    if (cacheBb_)
+    if
+    (
+        cacheBb_
+      ? !bbs_[index].overlaps(centre, radiusSqr)
+      : !treeBoundBox(patch_.points(),patch_[index]).overlaps(centre, radiusSqr)
+    )
     {
-        if (!bbs_[index].overlaps(centre, radiusSqr))
-        {
-            return false;
-        }
-    }
-    else
-    {
-        if (!calcBb(patch_.points(), patch_[index]).overlaps(centre, radiusSqr))
-        {
-            return false;
-        }
+        return false;
     }
 
     const pointField& points = patch_.points();
diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.H b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H
index b1dd87b517..ad0a00b1ff 100644
--- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.H
+++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H
@@ -80,9 +80,6 @@ class treeDataPrimitivePatch
 
     // Private Member Functions
 
-        //- Calculate face bounding box
-        static treeBoundBox calcBb(const pointField&, const face&);
-
         //- Initialise all member data
         void update();
 
diff --git a/src/meshTools/searchableSurface/searchableSurfaceCollection.C b/src/meshTools/searchableSurface/searchableSurfaceCollection.C
index e85200a8b0..e61b690e9e 100644
--- a/src/meshTools/searchableSurface/searchableSurfaceCollection.C
+++ b/src/meshTools/searchableSurface/searchableSurfaceCollection.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -244,8 +244,8 @@ Foam::searchableSurfaceCollection::searchableSurfaceCollection
     subGeom_.setSize(surfI);
     indexOffset_.setSize(surfI+1);
 
-    // Bounds is the overall bounds
-    bounds() = boundBox(point::max, point::min);
+    // Bounds is the overall bounds - prepare for min/max ops
+    bounds() = boundBox::invertedBox;
 
     forAll(subGeom_, surfI)
     {
diff --git a/src/meshTools/searchableSurface/searchableSurfacesQueries.C b/src/meshTools/searchableSurface/searchableSurfacesQueries.C
index 12b88dbbc9..9b8788c1a4 100644
--- a/src/meshTools/searchableSurface/searchableSurfacesQueries.C
+++ b/src/meshTools/searchableSurface/searchableSurfacesQueries.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -659,18 +659,14 @@ Foam::boundBox Foam::searchableSurfacesQueries::bounds
     const labelList& surfacesToTest
 )
 {
-    pointField bbPoints(2*surfacesToTest.size());
+    boundBox bb(boundBox::invertedBox);
 
-    forAll(surfacesToTest, testI)
+    forAll(surfacesToTest, testi)
     {
-        const searchableSurface& surface(allSurfaces[surfacesToTest[testI]]);
-
-        bbPoints[2*testI] = surface.bounds().min();
-
-        bbPoints[2*testI + 1] = surface.bounds().max();
+        bb.add(allSurfaces[surfacesToTest[testi]].bounds());
     }
 
-    return boundBox(bbPoints);
+    return bb;
 }
 
 
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
index cc54720dd3..73a5b67fbd 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,13 +55,9 @@ Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
     const point& fc
 )
 {
-    treeBoundBox bb(fc, fc);
-    forAll(f, fp)
-    {
-        const point& pt = points[f[fp]];
-        bb.min() = min(bb.min(), pt);
-        bb.max() = max(bb.max(), pt);
-    }
+    treeBoundBox bb(fc);
+    bb.add(points, f);
+
     return bb;
 }
 
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
index 830fc436e7..68fdd723cb 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -891,28 +891,19 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 
     // Find bounding box for all triangles on new distribution.
 
-    // Initialise to inverted box (VGREAT, -VGREAT)
+    // Initialise to inverted box (GREAT, -GREAT)
     List<List<treeBoundBox>> bbs(Pstream::nProcs());
     forAll(bbs, procI)
     {
-        bbs[procI].setSize(1);
-        //bbs[procI][0] = boundBox::invertedBox;
-        bbs[procI][0].min() = point( VGREAT,  VGREAT,  VGREAT);
-        bbs[procI][0].max() = point(-VGREAT, -VGREAT, -VGREAT);
+        bbs[procI].setSize(1, treeBoundBox::invertedBox);
     }
 
     forAll(s, triI)
     {
-        point& bbMin = bbs[distribution[triI]][0].min();
-        point& bbMax = bbs[distribution[triI]][0].max();
-
         const triSurface::FaceType& f = s[triI];
-        forAll(f, fp)
-        {
-            const point& pt = s.points()[f[fp]];
-            bbMin = ::Foam::min(bbMin, pt);
-            bbMax = ::Foam::max(bbMax, pt);
-        }
+
+        treeBoundBox& bb = bbs[distribution[triI]][0];
+        bb.add(s.points(), f);
     }
 
     // Now combine for all processors and convert to correct format.
@@ -936,17 +927,14 @@ bool Foam::distributedTriSurfaceMesh::overlaps
     const point& p2
 )
 {
+    treeBoundBox triBb(p0);
+    triBb.add(p1);
+    triBb.add(p2);
+
     forAll(bbs, bbI)
     {
         const treeBoundBox& bb = bbs[bbI];
 
-        treeBoundBox triBb(p0, p0);
-        triBb.min() = min(triBb.min(), p1);
-        triBb.min() = min(triBb.min(), p2);
-
-        triBb.max() = max(triBb.max(), p1);
-        triBb.max() = max(triBb.max(), p2);
-
         // Exact test of triangle intersecting bb
 
         // Quick rejection. If whole bounding box of tri is outside cubeBb then
@@ -2022,7 +2010,7 @@ Foam::triSurface Foam::distributedTriSurfaceMesh::overlappingSurface
     const scalar eps = 1.0e-4;
     forAll(bbs, i)
     {
-        const point mid = 0.5*(bbs[i].min() + bbs[i].max());
+        const point mid = bbs[i].midpoint();
         const vector halfSpan = (1.0+eps)*(bbs[i].max() - mid);
 
         bbsX[i].min() = mid - halfSpan;
diff --git a/src/sampling/meshToMesh/meshToMeshParallelOps.C b/src/sampling/meshToMesh/meshToMeshParallelOps.C
index 3d2173693f..b50341c42b 100644
--- a/src/sampling/meshToMesh/meshToMeshParallelOps.C
+++ b/src/sampling/meshToMesh/meshToMeshParallelOps.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -177,14 +177,13 @@ Foam::autoPtr<Foam::mapDistribute> Foam::meshToMesh::calcProcMap
             const cell& c = cells[celli];
 
             // determine bounding box of tgt cell
-            boundBox cellBb(point::max, point::min);
+            boundBox cellBb(boundBox::invertedBox);
             forAll(c, facei)
             {
                 const face& f = faces[c[facei]];
                 forAll(f, fp)
                 {
-                    cellBb.min() = min(cellBb.min(), points[f[fp]]);
-                    cellBb.max() = max(cellBb.max(), points[f[fp]]);
+                    cellBb.add(points, f);
                 }
             }
 
diff --git a/src/sampling/probes/patchProbes.C b/src/sampling/probes/patchProbes.C
index 80694beb77..e550f43533 100644
--- a/src/sampling/probes/patchProbes.C
+++ b/src/sampling/probes/patchProbes.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -79,12 +79,9 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
             {
                 bndFaces[nFaces++] = pp.start()+i;
                 const face& f = pp[i];
-                forAll(f, fp)
-                {
-                    const point& pt = pp.points()[f[fp]];
-                    overallBb.min() = min(overallBb.min(), pt);
-                    overallBb.max() = max(overallBb.max(), pt);
-                }
+
+                // Without reduction.
+                overallBb.add(pp.points(), f);
             }
         }
 
diff --git a/src/sampling/sampledSet/patchCloud/patchCloudSet.C b/src/sampling/sampledSet/patchCloud/patchCloudSet.C
index 93b1f9c6d0..c5833eeac7 100644
--- a/src/sampling/sampledSet/patchCloud/patchCloudSet.C
+++ b/src/sampling/sampledSet/patchCloud/patchCloudSet.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -74,7 +74,7 @@ void Foam::patchCloudSet::calcSamples
 
     labelList patchFaces(sz);
     sz = 0;
-    treeBoundBox bb(point::max, point::min);
+    treeBoundBox bb(treeBoundBox::invertedBox);
     forAllConstIter(labelHashSet, patchSet_, iter)
     {
         const polyPatch& pp = mesh().boundaryMesh()[iter.key()];
@@ -84,11 +84,8 @@ void Foam::patchCloudSet::calcSamples
             patchFaces[sz++] = pp.start()+i;
         }
 
-        // Do not do reduction.
-        const boundBox patchBb(pp.points(), pp.meshPoints(), false);
-
-        bb.min() = min(bb.min(), patchBb.min());
-        bb.max() = max(bb.max(), patchBb.max());
+        // Without reduction.
+        bb.add(pp.points(), pp.meshPoints());
     }
 
     // Not very random
@@ -97,8 +94,8 @@ void Foam::patchCloudSet::calcSamples
     bb = bb.extend(rndGen, 1e-4);
 
     // Make sure bb is 3D.
-    bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-    bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    bb.min() -= point::uniform(ROOTVSMALL);
+    bb.max() += point::uniform(ROOTVSMALL);
 
 
     indexedOctree<treeDataFace> patchTree
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index 3b225adc6a..f79bdf3a68 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -778,8 +778,8 @@ bool Foam::sampledTriSurfaceMesh::update()
         surface_.triSurface::points(),
         surface_.triSurface::meshPoints()
     );
-    bb.min() = max(bb.min(), mesh().bounds().min());
-    bb.max() = min(bb.max(), mesh().bounds().max());
+
+    bb.intersect(mesh().bounds());
 
     // Extend a bit
     const vector span(bb.span());
diff --git a/src/triSurface/triSurface/triSurface.C b/src/triSurface/triSurface/triSurface.C
index eb262d6cea..2ea0e6eea8 100644
--- a/src/triSurface/triSurface/triSurface.C
+++ b/src/triSurface/triSurface/triSurface.C
@@ -1106,7 +1106,7 @@ void Foam::triSurface::writeStats(Ostream& os) const
     PackedBoolList pointIsUsed(points().size());
 
     label nPoints = 0;
-    boundBox bb = boundBox::invertedBox;
+    boundBox bb(boundBox::invertedBox);
 
     forAll(*this, facei)
     {
@@ -1117,9 +1117,8 @@ void Foam::triSurface::writeStats(Ostream& os) const
             label pointi = f[fp];
             if (pointIsUsed.set(pointi, 1))
             {
-                bb.min() = ::Foam::min(bb.min(), points()[pointi]);
-                bb.max() = ::Foam::max(bb.max(), points()[pointi]);
-                nPoints++;
+                bb.add(points()[pointi]);
+                ++nPoints;
             }
         }
     }