diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
index 2e90e4d62c..30fbae0a1a 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
+++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
@@ -128,7 +128,20 @@ Foam::plane::plane(const vector& normalVector)
 :
     unitVector_(normalVector),
     basePoint_(vector::zero)
-{}
+{
+    scalar magUnitVector(mag(unitVector_));
+
+    if (magUnitVector > VSMALL)
+    {
+        unitVector_ /= magUnitVector;
+    }
+    else
+    {
+        FatalErrorIn("plane::plane(const vector&)")
+        << "plane normal has got zero length"
+            << abort(FatalError);
+    }
+}
 
 
 // Construct from point and normal vector
diff --git a/src/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C b/src/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
index 8cf27498a3..71361d826f 100644
--- a/src/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
+++ b/src/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
@@ -83,14 +83,20 @@ std::vector<Vb::Point> uniformGrid::initialPoints() const
     Random rndGen(1735621);
     scalar pert = randomPerturbationCoeff_*cmptMin(delta);
 
-    pointField points(ni*nj*nk);
-
-    label pI = 0;
+    std::vector<Vb::Point> initialPoints;
 
     for (int i = 0; i < ni; i++)
     {
         for (int j = 0; j < nj; j++)
         {
+            // Generating, testing and adding points one line at a time to
+            // reduce the memory requirement for cases with bounding boxes that
+            // are very large in comparison to the volume to be filled
+
+            label pI = 0;
+
+            pointField points(nk);
+
             for (int k = 0; k < nk; k++)
             {
                 point p
@@ -109,24 +115,22 @@ std::vector<Vb::Point> uniformGrid::initialPoints() const
 
                 points[pI++] = p;
             }
-        }
-    }
 
-    std::vector<Vb::Point> initialPoints;
+            Field<bool> insidePoints = cvMesh_.geometryToConformTo().wellInside
+            (
+                points,
+                minimumSurfaceDistance_*minimumSurfaceDistance_
+            );
 
-    Field<bool> insidePoints = cvMesh_.geometryToConformTo().wellInside
-    (
-        points,
-        minimumSurfaceDistance_*minimumSurfaceDistance_
-    );
+            forAll(insidePoints, i)
+            {
+                if (insidePoints[i])
+                {
+                    const point& p(points[i]);
 
-    forAll(insidePoints, i)
-    {
-        if (insidePoints[i])
-        {
-            const point& p(points[i]);
-
-            initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
+                    initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
+                }
+            }
         }
     }
 
diff --git a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
index 7f86c0f658..aceef22a5e 100644
--- a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
@@ -79,6 +79,23 @@ bool Foam::distributedTriSurfaceMesh::read()
     // Merge distance
     mergeDist_ = readScalar(dict_.lookup("mergeDistance"));
 
+    // Calculate the overall boundBox of the undistributed surface
+    point overallMin(VGREAT, VGREAT, VGREAT);
+    point overallMax(-VGREAT, -VGREAT, -VGREAT);
+
+    forAll(procBb_, procI)
+    {
+        const List<treeBoundBox>& bbs = procBb_[procI];
+
+        forAll(bbs, bbI)
+        {
+            overallMin = ::Foam::min(overallMin, bbs[bbI].min());
+            overallMax = ::Foam::max(overallMax, bbs[bbI].max());
+        }
+    }
+
+    bounds() = boundBox(overallMin, overallMax);
+
     return true;
 }
 
@@ -876,7 +893,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
         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][0].max() = point(-VGREAT, -VGREAT, -VGREAT);
     }
 
     forAll (s, triI)
@@ -2172,6 +2189,23 @@ void Foam::distributedTriSurfaceMesh::distribute
         {
             procBb_.transfer(newProcBb);
             dict_.set("bounds", procBb_[Pstream::myProcNo()]);
+
+            // Calculate the overall boundBox of the undistributed surface
+            point overallMin(VGREAT, VGREAT, VGREAT);
+            point overallMax(-VGREAT, -VGREAT, -VGREAT);
+
+            forAll(procBb_, procI)
+            {
+                const List<treeBoundBox>& bbs = procBb_[procI];
+
+                forAll(bbs, bbI)
+                {
+                    overallMin = ::Foam::min(overallMin, bbs[bbI].min());
+                    overallMax = ::Foam::max(overallMax, bbs[bbI].max());
+                }
+            }
+
+            bounds() = boundBox(overallMin, overallMax);
         }
     }
 
diff --git a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H
index 4d5f7dfd17..fbdcfc83d2 100644
--- a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H
+++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H
@@ -69,7 +69,7 @@ inline bool contiguous<segment>() {return contiguous<point>();}
 
 
 /*---------------------------------------------------------------------------*\
-                           Class distributedTriSurfaceMesh Declaration
+                   Class distributedTriSurfaceMesh Declaration
 \*---------------------------------------------------------------------------*/
 
 class distributedTriSurfaceMesh
diff --git a/src/meshTools/searchableSurface/searchableCylinder.C b/src/meshTools/searchableSurface/searchableCylinder.C
index 2aafe100ee..e495ea479f 100644
--- a/src/meshTools/searchableSurface/searchableCylinder.C
+++ b/src/meshTools/searchableSurface/searchableCylinder.C
@@ -354,17 +354,47 @@ void Foam::searchableCylinder::findLineAll
 
 Foam::boundBox Foam::searchableCylinder::calcBounds() const
 {
-    // Approximating the boundBox by the extents of spheres, centred at the
-    // endpoints, of the same radius as the cylinder
 
-    pointField bbPoints(4);
+    // Adapted from http://www.gamedev.net/community/forums/topic.asp?topic_id=338522&forum_id=20&gforum_id=0
 
-    bbPoints[0] = point1_ + radius_*vector::one;
-    bbPoints[1] = point1_ - radius_*vector::one;
-    bbPoints[2] = point2_ + radius_*vector::one;
-    bbPoints[3] = point2_ - radius_*vector::one;
+    // Let cylinder have end points A,B and radius r,
 
-    return boundBox(bbPoints);
+    // Bounds in direction X (same for Y and Z) can be found as:
+    // Let A.X<B.X (otherwise swap points)
+    // Good approximate lowest bound is A.X-r and highest is B.X+r (precise for
+    // capsule). At worst, in one direction it can be larger than needed by 2*r.
+
+    // Accurate bounds for cylinder is
+    // A.X-kx*r, B.X+kx*r
+    // where
+    // kx=sqrt(((A.Y-B.Y)^2+(A.Z-B.Z)^2)/((A.X-B.X)^2+(A.Y-B.Y)^2+(A.Z-B.Z)^2))
+
+    // similar thing for Y and Z
+    // (i.e.
+    // ky=sqrt(((A.X-B.X)^2+(A.Z-B.Z)^2)/((A.X-B.X)^2+(A.Y-B.Y)^2+(A.Z-B.Z)^2))
+    // kz=sqrt(((A.X-B.X)^2+(A.Y-B.Y)^2)/((A.X-B.X)^2+(A.Y-B.Y)^2+(A.Z-B.Z)^2))
+    // )
+
+    // How derived: geometric reasoning. Bounds of cylinder is same as for 2
+    // circles centered on A and B. This sqrt thingy gives sine of angle between
+    // axis and direction, used to find projection of radius.
+
+    vector kr
+    (
+        sqrt(sqr(unitDir_.y()) + sqr(unitDir_.z())),
+        sqrt(sqr(unitDir_.x()) + sqr(unitDir_.z())),
+        sqrt(sqr(unitDir_.x()) + sqr(unitDir_.y()))
+    );
+
+    kr *= radius_;
+
+    point min = point1_ - kr;
+    point max = point1_ + kr;
+
+    min = ::Foam::min(min, point2_ - kr);
+    max = ::Foam::max(max, point2_ + kr);
+
+    return boundBox(min, max);
 }
 
 
diff --git a/src/meshTools/searchableSurface/searchablePlane.C b/src/meshTools/searchableSurface/searchablePlane.C
index cb2a4e2fb9..3fb9c458ed 100644
--- a/src/meshTools/searchableSurface/searchablePlane.C
+++ b/src/meshTools/searchableSurface/searchablePlane.C
@@ -67,6 +67,31 @@ Foam::pointIndexHit Foam::searchablePlane::findLine
 }
 
 
+Foam::boundBox Foam::searchablePlane::calcBounds() const
+{
+    point max(VGREAT, VGREAT, VGREAT);
+
+    if (mag(normal() & vector(1,0,0)) > 1 - SMALL)
+    {
+        max.x() = 0;
+    }
+
+    if (mag(normal() & vector(0,1,0)) > 1 - SMALL)
+    {
+        max.y() = 0;
+    }
+
+    if (mag(normal() & vector(0,0,1)) > 1 - SMALL)
+    {
+        max.z() = 0;
+    }
+
+    point min = -max;
+
+    return boundBox(min, max);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::searchablePlane::searchablePlane
@@ -78,7 +103,9 @@ Foam::searchablePlane::searchablePlane
 :
     searchableSurface(io),
     plane(basePoint, normal)
-{}
+{
+    bounds() = calcBounds();
+}
 
 
 Foam::searchablePlane::searchablePlane
@@ -89,7 +116,9 @@ Foam::searchablePlane::searchablePlane
 :
     searchableSurface(io),
     plane(dict)
-{}
+{
+    bounds() = calcBounds();
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
diff --git a/src/meshTools/searchableSurface/searchablePlane.H b/src/meshTools/searchableSurface/searchablePlane.H
index 87777e3209..b18bb739f1 100644
--- a/src/meshTools/searchableSurface/searchablePlane.H
+++ b/src/meshTools/searchableSurface/searchablePlane.H
@@ -70,6 +70,8 @@ private:
             const point& end
         ) const;
 
+        //- Return the boundBox of the plane
+        boundBox calcBounds() const;
 
         //- Disallow default bitwise copy construct
         searchablePlane(const searchablePlane&);
diff --git a/src/meshTools/searchableSurface/searchableSphere.C b/src/meshTools/searchableSurface/searchableSphere.C
index 1de89fa019..a75fced3a4 100644
--- a/src/meshTools/searchableSurface/searchableSphere.C
+++ b/src/meshTools/searchableSurface/searchableSphere.C
@@ -136,8 +136,8 @@ Foam::searchableSphere::searchableSphere
 {
     bounds() = boundBox
     (
-        centre_ + radius_*vector::one,
-        centre_ - radius_*vector::one
+        centre_ - radius_*vector::one,
+        centre_ + radius_*vector::one
     );
 }
 
@@ -154,8 +154,8 @@ Foam::searchableSphere::searchableSphere
 {
     bounds() = boundBox
     (
-        centre_ + radius_*vector::one,
-        centre_ - radius_*vector::one
+        centre_ - radius_*vector::one,
+        centre_ + radius_*vector::one
     );
 }
 
diff --git a/src/meshTools/searchableSurface/searchableSurfacesQueries.C b/src/meshTools/searchableSurface/searchableSurfacesQueries.C
index 8ae7e0b1cc..8f0eb8bdac 100644
--- a/src/meshTools/searchableSurface/searchableSurfacesQueries.C
+++ b/src/meshTools/searchableSurface/searchableSurfacesQueries.C
@@ -552,93 +552,92 @@ void Foam::searchableSurfacesQueries::findAllIntersections
 }
 
 
-//// Find intersections of edge nearest to both endpoints.
-//void Foam::searchableSurfacesQueries::findNearestIntersection
-//(
-//    const PtrList<searchableSurface>& allSurfaces,
-//    const labelList& surfacesToTest,
-//    const pointField& start,
-//    const pointField& end,
-//
-//    labelList& surface1,
-//    List<pointIndexHit>& hit1,
-//    labelList& surface2,
-//    List<pointIndexHit>& hit2
-//)
-//{
-//    // 1. intersection from start to end
-//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-//
-//    // Initialize arguments
-//    surface1.setSize(start.size());
-//    surface1 = -1;
-//    hit1.setSize(start.size());
-//
-//    // Current end of segment to test.
-//    pointField nearest(end);
-//    // Work array
-//    List<pointIndexHit> nearestInfo(start.size());
-//
-//    forAll(surfacesToTest, testI)
-//    {
-//        // See if any intersection between start and current nearest
-//        allSurfaces[surfacesToTest[testI]].findLine
-//        (
-//            start,
-//            nearest,
-//            nearestInfo
-//        );
-//
-//        forAll(nearestInfo, pointI)
-//        {
-//            if (nearestInfo[pointI].hit())
-//            {
-//                hit1[pointI] = nearestInfo[pointI];
-//                surface1[pointI] = testI;
-//                nearest[pointI] = hit1[pointI].hitPoint();
-//            }
-//        }
-//    }
-//
-//
-//    // 2. intersection from end to last intersection
-//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-//
-//    // Find the nearest intersection from end to start. Note that we
-//    // initialize to the first intersection (if any).
-//    surface2 = surface1;
-//    hit2 = hit1;
-//
-//    // Set current end of segment to test.
-//    forAll(nearest, pointI)
-//    {
-//        if (hit1[pointI].hit())
-//        {
-//            nearest[pointI] = hit1[pointI].hitPoint();
-//        }
-//        else
-//        {
-//            // Disable testing by setting to end.
-//            nearest[pointI] = end[pointI];
-//        }
-//    }
-//
-//    forAll(surfacesToTest, testI)
-//    {
-//        // See if any intersection between end and current nearest
-//        allSurfaces[surfacesToTest[i]].findLine(end, nearest, nearestInfo);
-//
-//        forAll(nearestInfo, pointI)
-//        {
-//            if (nearestInfo[pointI].hit())
-//            {
-//                hit2[pointI] = nearestInfo[pointI];
-//                surface2[pointI] = testI;
-//                nearest[pointI] = hit2[pointI].hitPoint();
-//            }
-//        }
-//    }
-//}
+//Find intersections of edge nearest to both endpoints.
+void Foam::searchableSurfacesQueries::findNearestIntersection
+(
+   const PtrList<searchableSurface>& allSurfaces,
+   const labelList& surfacesToTest,
+   const pointField& start,
+   const pointField& end,
+   labelList& surface1,
+   List<pointIndexHit>& hit1,
+   labelList& surface2,
+   List<pointIndexHit>& hit2
+)
+{
+   // 1. intersection from start to end
+   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+   // Initialize arguments
+   surface1.setSize(start.size());
+   surface1 = -1;
+   hit1.setSize(start.size());
+
+   // Current end of segment to test.
+   pointField nearest(end);
+   // Work array
+   List<pointIndexHit> nearestInfo(start.size());
+
+   forAll(surfacesToTest, testI)
+   {
+       // See if any intersection between start and current nearest
+       allSurfaces[surfacesToTest[testI]].findLine
+       (
+           start,
+           nearest,
+           nearestInfo
+       );
+
+       forAll(nearestInfo, pointI)
+       {
+           if (nearestInfo[pointI].hit())
+           {
+               hit1[pointI] = nearestInfo[pointI];
+               surface1[pointI] = testI;
+               nearest[pointI] = hit1[pointI].hitPoint();
+           }
+       }
+   }
+
+
+   // 2. intersection from end to last intersection
+   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+   // Find the nearest intersection from end to start. Note that we
+   // initialize to the first intersection (if any).
+   surface2 = surface1;
+   hit2 = hit1;
+
+   // Set current end of segment to test.
+   forAll(nearest, pointI)
+   {
+       if (hit1[pointI].hit())
+       {
+           nearest[pointI] = hit1[pointI].hitPoint();
+       }
+       else
+       {
+           // Disable testing by setting to end.
+           nearest[pointI] = end[pointI];
+       }
+   }
+
+   forAll(surfacesToTest, testI)
+   {
+       // See if any intersection between end and current nearest
+       allSurfaces[surfacesToTest[testI]].findLine(end, nearest, nearestInfo);
+
+       forAll(nearestInfo, pointI)
+       {
+           if (nearestInfo[pointI].hit())
+           {
+               hit2[pointI] = nearestInfo[pointI];
+               surface2[pointI] = testI;
+               nearest[pointI] = hit2[pointI].hitPoint();
+           }
+       }
+   }
+}
 
 
 // Find nearest. Return -1 or nearest point
diff --git a/src/meshTools/searchableSurface/searchableSurfacesQueries.H b/src/meshTools/searchableSurface/searchableSurfacesQueries.H
index 7e336fc68f..ff6e0695bf 100644
--- a/src/meshTools/searchableSurface/searchableSurfacesQueries.H
+++ b/src/meshTools/searchableSurface/searchableSurfacesQueries.H
@@ -160,6 +160,19 @@ public:
                 List<List<pointIndexHit> >& surfaceHits
             );
 
+            //Find intersections of edge nearest to both endpoints.
+            static void findNearestIntersection
+            (
+                const PtrList<searchableSurface>& allSurfaces,
+                const labelList& surfacesToTest,
+                const pointField& start,
+                const pointField& end,
+                labelList& surface1,
+                List<pointIndexHit>& hit1,
+                labelList& surface2,
+                List<pointIndexHit>& hit2
+            );
+
             //- Find nearest. Return -1 (and a miss()) or surface and nearest
             //  point.
             static void findNearest