Added proper bounds() to searchableCylinder and searchablePlane. Fixed error in searchableSphere. Fixed constructor in plane. Added overall boundBox calculation to distributedTriSurfaceMesh. Reinstated findNearestIntersection in searchableSurfacesQueries.

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

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);
+    std::vector<Vb::Point> initialPoints;
-
-    label pI = 0;
 
     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,10 +115,6 @@ std::vector<Vb::Point> uniformGrid::initialPoints() const
 
                 points[pI++] = p;
             }
-        }
-    }
-
-    std::vector<Vb::Point> initialPoints;
 
             Field<bool> insidePoints = cvMesh_.geometryToConformTo().wellInside
             (
@@ -129,6 +131,8 @@ std::vector<Vb::Point> uniformGrid::initialPoints() const
                     initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
                 }
             }
+        }
+    }
 
     return initialPoints;
 }

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;
 }
 
@@ -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);
         }
     }
 

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;
+    // Let cylinder have end points A,B and radius r,
-    bbPoints[1] = point1_ - radius_*vector::one;
-    bbPoints[2] = point2_ + radius_*vector::one;
-    bbPoints[3] = point2_ - radius_*vector::one;
 
-    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);
 }
 
 

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  * * * * * * * * * * * * * * * //

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&);

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
     );
 }
 

src/meshTools/searchableSurface/searchableSurfacesQueries.C

@@ -552,93 +552,92 @@ void Foam::searchableSurfacesQueries::findAllIntersections
 }
 
 
-//// Find intersections of edge nearest to both endpoints.
+//Find intersections of edge nearest to both endpoints.
-//void Foam::searchableSurfacesQueries::findNearestIntersection
+void Foam::searchableSurfacesQueries::findNearestIntersection
-//(
+(
-//    const PtrList<searchableSurface>& allSurfaces,
+   const PtrList<searchableSurface>& allSurfaces,
-//    const labelList& surfacesToTest,
+   const labelList& surfacesToTest,
-//    const pointField& start,
+   const pointField& start,
-//    const pointField& end,
+   const pointField& end,
-//
+   labelList& surface1,
-//    labelList& surface1,
+   List<pointIndexHit>& hit1,
-//    List<pointIndexHit>& hit1,
+   labelList& surface2,
-//    labelList& surface2,
+   List<pointIndexHit>& hit2
-//    List<pointIndexHit>& hit2
+)
-//)
+{
-//{
+   // 1. intersection from start to end
-//    // 1. intersection from start to end
+   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
-//
+   // Initialize arguments
-//    // Initialize arguments
+   surface1.setSize(start.size());
-//    surface1.setSize(start.size());
+   surface1 = -1;
-//    surface1 = -1;
+   hit1.setSize(start.size());
-//    hit1.setSize(start.size());
+
-//
+   // Current end of segment to test.
-//    // Current end of segment to test.
+   pointField nearest(end);
-//    pointField nearest(end);
+   // Work array
-//    // Work array
+   List<pointIndexHit> nearestInfo(start.size());
-//    List<pointIndexHit> nearestInfo(start.size());
+
-//
+   forAll(surfacesToTest, testI)
-//    forAll(surfacesToTest, testI)
+   {
-//    {
+       // See if any intersection between start and current nearest
-//        // See if any intersection between start and current nearest
+       allSurfaces[surfacesToTest[testI]].findLine
-//        allSurfaces[surfacesToTest[testI]].findLine
+       (
-//        (
+           start,
-//            start,
+           nearest,
-//            nearest,
+           nearestInfo
-//            nearestInfo
+       );
-//        );
+
-//
+       forAll(nearestInfo, pointI)
-//        forAll(nearestInfo, pointI)
+       {
-//        {
+           if (nearestInfo[pointI].hit())
-//            if (nearestInfo[pointI].hit())
+           {
-//            {
+               hit1[pointI] = nearestInfo[pointI];
-//                hit1[pointI] = nearestInfo[pointI];
+               surface1[pointI] = testI;
-//                surface1[pointI] = testI;
+               nearest[pointI] = hit1[pointI].hitPoint();
-//                nearest[pointI] = hit1[pointI].hitPoint();
+           }
-//            }
+       }
-//        }
+   }
-//    }
+
-//
+
-//
+   // 2. intersection from end to last intersection
-//    // 2. intersection from end to last intersection
+   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
-//
+   // Find the nearest intersection from end to start. Note that we
-//    // Find the nearest intersection from end to start. Note that we
+   // initialize to the first intersection (if any).
-//    // initialize to the first intersection (if any).
+   surface2 = surface1;
-//    surface2 = surface1;
+   hit2 = hit1;
-//    hit2 = hit1;
+
-//
+   // Set current end of segment to test.
-//    // Set current end of segment to test.
+   forAll(nearest, pointI)
-//    forAll(nearest, pointI)
+   {
-//    {
+       if (hit1[pointI].hit())
-//        if (hit1[pointI].hit())
+       {
-//        {
+           nearest[pointI] = hit1[pointI].hitPoint();
-//            nearest[pointI] = hit1[pointI].hitPoint();
+       }
-//        }
+       else
-//        else
+       {
-//        {
+           // Disable testing by setting to end.
-//            // Disable testing by setting to end.
+           nearest[pointI] = end[pointI];
-//            nearest[pointI] = end[pointI];
+       }
-//        }
+   }
-//    }
+
-//
+   forAll(surfacesToTest, testI)
-//    forAll(surfacesToTest, testI)
+   {
-//    {
+       // See if any intersection between end and current nearest
-//        // See if any intersection between end and current nearest
+       allSurfaces[surfacesToTest[testI]].findLine(end, nearest, nearestInfo);
-//        allSurfaces[surfacesToTest[i]].findLine(end, nearest, nearestInfo);
+
-//
+       forAll(nearestInfo, pointI)
-//        forAll(nearestInfo, pointI)
+       {
-//        {
+           if (nearestInfo[pointI].hit())
-//            if (nearestInfo[pointI].hit())
+           {
-//            {
+               hit2[pointI] = nearestInfo[pointI];
-//                hit2[pointI] = nearestInfo[pointI];
+               surface2[pointI] = testI;
-//                surface2[pointI] = testI;
+               nearest[pointI] = hit2[pointI].hitPoint();
-//                nearest[pointI] = hit2[pointI].hitPoint();
+           }
-//            }
+       }
-//        }
+   }
-//    }
+}
-//}
 
 
 // Find nearest. Return -1 or nearest point

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