ENH: Expand triSurfaceSearch and add triSurfaceRegionSearch

+ Move most of the octree searching functionality from triSurfaceMesh into
  triSurfaceSearch. Much of it was duplicated anyway.
+ Add triSurfaceRegionSearch, which constructs an octree for each region in
  a surface and performs searches on specified regions only.

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H

@@ -54,6 +54,7 @@ namespace Foam
 
 class polyMesh;
 class PackedBoolList;
+class boundBox;
 
 /*---------------------------------------------------------------------------*\
                          Class PatchTools Declaration
@@ -140,6 +141,21 @@ public:
         labelList& faceMap
     );
 
+    //-
+    template
+    <
+        class Face,
+        template<class> class FaceList,
+        class PointField,
+        class PointType
+    >
+    static void calcBounds
+    (
+        const PrimitivePatch<Face, FaceList, PointField, PointType>& p,
+        boundBox& bb,
+        label& nPoints
+    );
+
     //- Return edge-face addressing sorted by angle around the edge.
     //  Orientation is anticlockwise looking from edge.vec(localPoints())
     template
@@ -271,7 +287,6 @@ public:
         List<Face>& mergedFaces,
         labelList& pointMergeMap
     );
-
 };
 
 

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSearch.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "PatchTools.H"
+#include "PackedBoolList.H"
+#include "boundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -222,4 +224,45 @@ Foam::PatchTools::subsetMap
 }
 
 
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+void Foam::PatchTools::calcBounds
+(
+    const PrimitivePatch<Face, FaceList, PointField, PointType>& p,
+    boundBox& bb,
+    label& nPoints
+)
+{
+    // Unfortunately nPoints constructs meshPoints() so do compact version
+    // ourselves
+    const PointField& points = p.points();
+
+    PackedBoolList pointIsUsed(points.size());
+
+    nPoints = 0;
+    bb = boundBox::invertedBox;
+
+    forAll(p, faceI)
+    {
+        const Face& f = p[faceI];
+
+        forAll(f, fp)
+        {
+            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++;
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveShapes/plane/plane.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,7 @@ Class
 
 Description
     Geometric class that creates a 2D plane and can return the intersection
-     point between a line and the plane.
+    point between a line and the plane.
 
 SourceFiles
     plane.C
@@ -62,6 +62,14 @@ class plane
 {
 public:
 
+    //- Side of the plane
+    enum side
+    {
+        NORMAL,
+        FLIP
+    };
+
+
         //- A direction and a reference point
         class ray
         {
@@ -182,6 +190,10 @@ public:
         //- Return the cutting point between this plane and two other planes
         point planePlaneIntersect(const plane&, const plane&) const;
 
+        //- Return the side of the plane that the point is on.
+        //  If the point is on the plane, then returns NORMAL.
+        side sideOfPlane(const point& p) const;
+
         //- Write to dictionary
         void writeDict(Ostream&) const;
 

src/meshTools/Make/files

@@ -160,6 +160,7 @@ $(intersectedSurface)/intersectedSurface.C
 $(intersectedSurface)/edgeSurface.C
 
 triSurface/triSurfaceSearch/triSurfaceSearch.C
+triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
 triSurface/triangleFuncs/triangleFuncs.C
 triSurface/surfaceFeatures/surfaceFeatures.C
 triSurface/triSurfaceTools/triSurfaceTools.C

src/meshTools/searchableSurface/triSurfaceMesh.C

@@ -29,7 +29,7 @@ License
 #include "EdgeMap.H"
 #include "triSurfaceFields.H"
 #include "Time.H"
-#include "PackedBoolList.H"
+#include "PatchTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -218,48 +218,6 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
 }
 
 
-void Foam::triSurfaceMesh::calcBounds(boundBox& bb, label& nPoints) const
-{
-    const triSurface& s = static_cast<const triSurface&>(*this);
-
-    calcBounds(s, bb, nPoints);
-}
-
-
-template <class PatchType>
-void Foam::triSurfaceMesh::calcBounds
-(
-    const PatchType& patch,
-    boundBox& bb,
-    label& nPoints
-) const
-{
-    // Unfortunately nPoints constructs meshPoints() so do compact version
-    // ourselves
-
-    PackedBoolList pointIsUsed(points()().size());
-
-    nPoints = 0;
-    bb = boundBox::invertedBox;
-
-    forAll(patch, faceI)
-    {
-        const typename PatchType::FaceType& f = patch[faceI];
-
-        forAll(f, fp)
-        {
-            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++;
-            }
-        }
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
@@ -279,9 +237,8 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
         )
     ),
     triSurface(s),
-    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
+    triSurfaceRegionSearch(s),
     minQuality_(-1),
-    maxTreeDepth_(10),
     surfaceClosed_(-1)
 {
     bounds() = boundBox(points());
@@ -326,9 +283,8 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
             searchableSurface::objectPath()
         )
     ),
-    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
+    triSurfaceRegionSearch(static_cast<const triSurface&>(*this)),
     minQuality_(-1),
-    maxTreeDepth_(10),
     surfaceClosed_(-1)
 {
     bounds() = boundBox(points());
@@ -376,9 +332,8 @@ Foam::triSurfaceMesh::triSurfaceMesh
             searchableSurface::objectPath()
         )
     ),
-    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
+    triSurfaceRegionSearch(static_cast<const triSurface&>(*this), dict),
     minQuality_(-1),
-    maxTreeDepth_(10),
     surfaceClosed_(-1)
 {
     scalar scaleFactor = 0;
@@ -394,13 +349,6 @@ Foam::triSurfaceMesh::triSurfaceMesh
 
     bounds() = boundBox(points());
 
-    // Have optional non-standard search tolerance for gappy surfaces.
-    if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0)
-    {
-        Info<< searchableSurface::name() << " : using intersection tolerance "
-            << tolerance_ << endl;
-    }
-
     // Have optional minimum quality for normal calculation
     if (dict.readIfPresent("minQuality", minQuality_) && minQuality_ > 0)
     {
@@ -408,13 +356,6 @@ Foam::triSurfaceMesh::triSurfaceMesh
             << " : ignoring triangles with quality < "
             << minQuality_ << " for normals calculation." << endl;
     }
-
-    // Have optional non-standard tree-depth to limit storage.
-    if (dict.readIfPresent("maxTreeDepth", maxTreeDepth_) && maxTreeDepth_ > 0)
-    {
-        Info<< searchableSurface::name() << " : using maximum tree depth "
-            << maxTreeDepth_ << endl;
-    }
 }
 
 
@@ -428,7 +369,7 @@ Foam::triSurfaceMesh::~triSurfaceMesh()
 
 void Foam::triSurfaceMesh::clearOut()
 {
-    tree_.clear();
+    triSurfaceRegionSearch::clearOut();
     edgeTree_.clear();
     triSurface::clearOut();
 }
@@ -470,176 +411,12 @@ bool Foam::triSurfaceMesh::overlaps(const boundBox& bb) const
 
 void Foam::triSurfaceMesh::movePoints(const pointField& newPoints)
 {
-    tree_.clear();
+    triSurfaceRegionSearch::clearOut();
     edgeTree_.clear();
     triSurface::movePoints(newPoints);
 }
 
 
-const Foam::indexedOctree<Foam::treeDataTriSurface>&
-Foam::triSurfaceMesh::tree() const
-{
-    if (tree_.empty())
-    {
-        // Calculate bb without constructing local point numbering.
-        treeBoundBox bb;
-        label nPoints;
-        calcBounds(bb, nPoints);
-
-        if (nPoints != points()().size())
-        {
-            WarningIn("triSurfaceMesh::tree() const")
-                << "Surface " << searchableSurface::name()
-                << " does not have compact point numbering."
-                << " Of " << points()().size() << " only " << nPoints
-                << " are used. This might give problems in some routines."
-                << endl;
-        }
-
-
-        // Random number generator. Bit dodgy since not exactly random ;-)
-        Random rndGen(65431);
-
-        // Slightly extended bb. Slightly off-centred just so on symmetric
-        // geometry there are less face/edge aligned items.
-        bb = bb.extend(rndGen, 1e-4);
-        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-        scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
-        indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
-
-        tree_.reset
-        (
-            new indexedOctree<treeDataTriSurface>
-            (
-                treeDataTriSurface(true, *this, tolerance_),
-                bb,
-                maxTreeDepth_,  // maxLevel
-                10,             // leafsize
-                3.0             // duplicity
-            )
-        );
-
-        indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
-    }
-
-    return tree_();
-}
-
-
-const Foam::PtrList
-<
-    Foam::indexedOctree
-    <
-        Foam::triSurfaceMesh::treeDataTriSurfacePrimitivePatch
-    >
->&
-Foam::triSurfaceMesh::treeByRegion() const
-{
-    if (treeByRegion_.empty())
-    {
-        Map<label> regionSizes;
-        forAll(*this, fI)
-        {
-            const label regionI = this->triSurface::operator[](fI).region();
-
-            regionSizes(regionI)++;
-        }
-
-        label nRegions = regionSizes.size();
-
-        indirectRegionPatches_.setSize(nRegions);
-        treeByRegion_.setSize(nRegions);
-
-        labelListList regionsAddressing(nRegions);
-
-        forAll(regionsAddressing, regionI)
-        {
-            regionsAddressing[regionI] = labelList(regionSizes[regionI], -1);
-        }
-
-        labelList nFacesInRegions(nRegions, 0);
-
-        forAll(*this, fI)
-        {
-            const label regionI = this->triSurface::operator[](fI).region();
-
-            regionsAddressing[regionI][nFacesInRegions[regionI]++] = fI;
-        }
-
-        forAll(regionsAddressing, regionI)
-        {
-            scalar oldTol =
-                indexedOctree<treeDataTriSurfacePrimitivePatch>::perturbTol();
-            indexedOctree<treeDataTriSurfacePrimitivePatch>::perturbTol() =
-                tolerance_;
-
-            indirectRegionPatches_.set
-            (
-                regionI,
-                new indirectTriSurfacePrimitivePatch
-                (
-                    IndirectList<labelledTri>
-                    (
-                        *this,
-                        regionsAddressing[regionI]
-                    ),
-                    points()
-                )
-            );
-
-            // Calculate bb without constructing local point numbering.
-            treeBoundBox bb;
-            label nPoints;
-            calcBounds(indirectRegionPatches_[regionI], bb, nPoints);
-
-            if (nPoints != points()().size())
-            {
-                WarningIn("triSurfaceMesh::treeByRegion() const")
-                    << "Surface " << searchableSurface::name()
-                    << " does not have compact point numbering."
-                    << " Of " << points()().size() << " only " << nPoints
-                    << " are used. This might give problems in some routines."
-                    << endl;
-            }
-
-            // Random number generator. Bit dodgy since not exactly random ;-)
-            Random rndGen(65431);
-
-            // Slightly extended bb. Slightly off-centred just so on symmetric
-            // geometry there are fewer face/edge aligned items.
-            bb = bb.extend(rndGen, 1e-4);
-            bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-            bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-            treeByRegion_.set
-            (
-                regionI,
-                new indexedOctree<treeDataTriSurfacePrimitivePatch>
-                (
-                    treeDataTriSurfacePrimitivePatch
-                    (
-                        true,
-                        indirectRegionPatches_[regionI],
-                        tolerance_
-                    ),
-                    bb,
-                    maxTreeDepth_,  // maxLevel
-                    10,             // leafsize
-                    3.0             // duplicity
-                )
-            );
-
-            indexedOctree<treeDataTriSurfacePrimitivePatch>::perturbTol() =
-                oldTol;
-        }
-    }
-
-    return treeByRegion_;
-}
-
-
 const Foam::indexedOctree<Foam::treeDataEdge>&
 Foam::triSurfaceMesh::edgeTree() const
 {
@@ -658,7 +435,12 @@ Foam::triSurfaceMesh::edgeTree() const
 
         treeBoundBox bb;
         label nPoints;
-        calcBounds(bb, nPoints);
+        PatchTools::calcBounds
+        (
+            static_cast<const triSurface&>(*this),
+            bb,
+            nPoints
+        );
 
         // Random number generator. Bit dodgy since not exactly random ;-)
         Random rndGen(65431);
@@ -671,7 +453,7 @@ Foam::triSurfaceMesh::edgeTree() const
         bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 
         scalar oldTol = indexedOctree<treeDataEdge>::perturbTol();
-        indexedOctree<treeDataEdge>::perturbTol() = tolerance_;
+        indexedOctree<treeDataEdge>::perturbTol() = tolerance();
 
         edgeTree_.reset
         (
@@ -685,7 +467,7 @@ Foam::triSurfaceMesh::edgeTree() const
                     bEdges          // selected edges
                 ),
                 bb,                 // bb
-                maxTreeDepth_,      // maxLevel
+                maxTreeDepth(),      // maxLevel
                 10,                 // leafsize
                 3.0                 // duplicity
             )
@@ -697,12 +479,6 @@ Foam::triSurfaceMesh::edgeTree() const
 }
 
 
-Foam::scalar Foam::triSurfaceMesh::tolerance() const
-{
-    return tolerance_;
-}
-
-
 const Foam::wordList& Foam::triSurfaceMesh::regions() const
 {
     if (regions_.empty())
@@ -743,24 +519,7 @@ void Foam::triSurfaceMesh::findNearest
     List<pointIndexHit>& info
 ) const
 {
-    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    triSurfaceSearch::findNearest(samples, nearestDistSqr, info);
-    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
-
-    const indexedOctree<treeDataTriSurface>& octree = tree();
-
-    info.setSize(samples.size());
-
-    forAll(samples, i)
-    {
-        static_cast<pointIndexHit&>(info[i]) = octree.findNearest
-        (
-            samples[i],
-            nearestDistSqr[i],
-            treeDataTriSurface::findNearestOp(octree)
-        );
-    }
-
-    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 
@@ -772,67 +531,13 @@ void Foam::triSurfaceMesh::findNearest
     List<pointIndexHit>& info
 ) const
 {
-    if (regionIndices.empty())
+    triSurfaceRegionSearch::findNearest
-    {
+    (
-        findNearest(samples, nearestDistSqr, info);
+        samples,
-    }
+        nearestDistSqr,
-    else
+        regionIndices,
-    {
+        info
-        scalar oldTol =
+    );
-            indexedOctree<treeDataTriSurfacePrimitivePatch>::perturbTol();
-        indexedOctree<treeDataTriSurfacePrimitivePatch>::perturbTol() =
-            tolerance_;
-
-        const PtrList<indexedOctree<treeDataTriSurfacePrimitivePatch> >&
-            octrees = treeByRegion();
-
-        info.setSize(samples.size());
-
-        forAll(octrees, treeI)
-        {
-            if (findIndex(regionIndices, treeI) == -1)
-            {
-                continue;
-            }
-
-            const indexedOctree<treeDataTriSurfacePrimitivePatch>& octree =
-                octrees[treeI];
-
-            forAll(samples, i)
-            {
-//                if (!octree.bb().contains(samples[i]))
-//                {
-//                    continue;
-//                }
-
-                pointIndexHit currentRegionHit = octree.findNearest
-                (
-                    samples[i],
-                    nearestDistSqr[i],
-                    treeDataTriSurfacePrimitivePatch::findNearestOp(octree)
-                );
-
-                if
-                (
-                    currentRegionHit.hit()
-                 &&
-                    (
-                        !info[i].hit()
-                     ||
-                        (
-                            magSqr(currentRegionHit.hitPoint() - samples[i])
-                          < magSqr(info[i].hitPoint() - samples[i])
-                        )
-                    )
-                )
-                {
-                    info[i] = currentRegionHit;
-                }
-            }
-        }
-
-        indexedOctree<treeDataTriSurfacePrimitivePatch>::perturbTol() = oldTol;
-    }
 }
 
 
@@ -843,23 +548,7 @@ void Foam::triSurfaceMesh::findLine
     List<pointIndexHit>& info
 ) const
 {
-    const indexedOctree<treeDataTriSurface>& octree = tree();
+    triSurfaceSearch::findLine(start, end, info);
-
-    info.setSize(start.size());
-
-    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
-    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
-
-    forAll(start, i)
-    {
-        static_cast<pointIndexHit&>(info[i]) = octree.findLine
-        (
-            start[i],
-            end[i]
-        );
-    }
-
-    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 
@@ -870,23 +559,7 @@ void Foam::triSurfaceMesh::findLineAny
     List<pointIndexHit>& info
 ) const
 {
-    const indexedOctree<treeDataTriSurface>& octree = tree();
+    triSurfaceSearch::findLineAny(start, end, info);
-
-    info.setSize(start.size());
-
-    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
-    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
-
-    forAll(start, i)
-    {
-        static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
-        (
-            start[i],
-            end[i]
-        );
-    }
-
-    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 
@@ -897,51 +570,7 @@ void Foam::triSurfaceMesh::findLineAll
     List<List<pointIndexHit> >& info
 ) const
 {
-    const indexedOctree<treeDataTriSurface>& octree = tree();
+    triSurfaceSearch::findLineAll(start, end, info);
-
-    info.setSize(start.size());
-
-    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
-    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
-
-    // Work array
-    DynamicList<pointIndexHit, 1, 1> hits;
-
-    DynamicList<label> shapeMask;
-
-    treeDataTriSurface::findAllIntersectOp allIntersectOp(octree, shapeMask);
-
-    forAll(start, pointI)
-    {
-        hits.clear();
-        shapeMask.clear();
-
-        while (true)
-        {
-            // See if any intersection between pt and end
-            pointIndexHit inter = octree.findLine
-            (
-                start[pointI],
-                end[pointI],
-                allIntersectOp
-            );
-
-            if (inter.hit())
-            {
-                hits.append(inter);
-
-                shapeMask.append(inter.index());
-            }
-            else
-            {
-                break;
-            }
-        }
-
-        info[pointI].transfer(hits);
-    }
-
-    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 
@@ -1103,7 +732,7 @@ void Foam::triSurfaceMesh::getVolumeType
     volType.setSize(points.size());
 
     scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
-    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
 
     forAll(points, pointI)
     {

src/meshTools/searchableSurface/triSurfaceMesh.H

@@ -50,6 +50,7 @@ SourceFiles
 #include "treeDataEdge.H"
 #include "EdgeMap.H"
 #include "triSurface.H"
+#include "triSurfaceRegionSearch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,44 +65,17 @@ class triSurfaceMesh
 :
     public searchableSurface,
     public objectRegistry,      // so we can store fields
-    public triSurface
+    public triSurface,
+    public triSurfaceRegionSearch
 {
 private:
 
-    // Private typedefs
-
-        typedef PrimitivePatch
-        <
-            labelledTri,
-            IndirectList,
-            const pointField&
-        > indirectTriSurfacePrimitivePatch;
-
-        typedef treeDataPrimitivePatch<indirectTriSurfacePrimitivePatch>
-            treeDataTriSurfacePrimitivePatch;
-
     // Private member data
 
-        //- Optional tolerance to use in searches
-        scalar tolerance_;
-
         //- Optional min triangle quality. Triangles below this get
         //  ignored for normal calculation
         scalar minQuality_;
 
-        //- Optional max tree depth of octree
-        label maxTreeDepth_;
-
-        //- Search tree (triangles)
-        mutable autoPtr<indexedOctree<treeDataTriSurface> > tree_;
-
-        mutable PtrList<indirectTriSurfacePrimitivePatch>
-            indirectRegionPatches_;
-
-        //- Search tree for each region
-        mutable PtrList<indexedOctree<treeDataTriSurfacePrimitivePatch> >
-            treeByRegion_;
-
         //- Search tree for boundary edges.
         mutable autoPtr<indexedOctree<treeDataEdge> > edgeTree_;
 
@@ -157,20 +131,6 @@ private:
         void operator=(const triSurfaceMesh&);
 
 
-protected:
-
-        //- Calculate (number of) used points and their bounding box
-        void calcBounds(boundBox& bb, label& nPoints) const;
-
-        template <class PatchType>
-        void calcBounds
-        (
-            const PatchType& patch,
-            boundBox& bb,
-            label& nPoints
-        ) const;
-
-
 public:
 
     //- Runtime type information
@@ -206,17 +166,9 @@ public:
         //- Move points
         virtual void movePoints(const pointField&);
 
-        //- Demand driven construction of octree
-        const indexedOctree<treeDataTriSurface>& tree() const;
-
-        //- Demand driven construction of octree by region
-        const PtrList<indexedOctree<treeDataTriSurfacePrimitivePatch> >&
-        treeByRegion() const;
-
         //- Demand driven contruction of octree for boundary edges
         const indexedOctree<treeDataEdge>& edgeTree() const;
 
-        scalar tolerance() const;
 
         // searchableSurface implementation
 

src/meshTools/triSurface/orientedSurface/orientedSurface.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -258,7 +258,8 @@ void Foam::orientedSurface::findZoneSide
     zoneFaceI = -1;
     isOutside = false;
 
-    List<pointIndexHit> hits;
+    pointField start(1, outsidePoint);
+    List<List<pointIndexHit> > hits(1, List<pointIndexHit>());
 
     forAll(faceZone, faceI)
     {
@@ -273,7 +274,7 @@ void Foam::orientedSurface::findZoneSide
             // Check if normal different enough to decide upon
             if (magD > SMALL && (mag(n & d/magD) > 1e-6))
             {
-                point end = fc + d;
+                pointField end(1, fc + d);
 
                 //Info<< "Zone " << zoneI << " : Shooting ray"
                 //    << " from " << outsidePoint
@@ -281,12 +282,13 @@ void Foam::orientedSurface::findZoneSide
                 //    << " to pierce triangle " << faceI
                 //    << " with centre " << fc << endl;
 
-                surfSearches.findLineAll(outsidePoint, end, hits);
+
+                surfSearches.findLineAll(start, end, hits);
 
                 label zoneIndex = -1;
-                forAll(hits, i)
+                forAll(hits[0], i)
                 {
-                    if (hits[i].index() == faceI)
+                    if (hits[0][i].index() == faceI)
                     {
                         zoneIndex = i;
                         break;
@@ -573,4 +575,19 @@ bool Foam::orientedSurface::orient
 }
 
 
+bool Foam::orientedSurface::orientAll(triSurface& s)
+{
+    if (orientConsistent(s))
+    {
+        Info<< "Some triangles are not oriented consistently. Aborting."
+            << endl;
+        return false;
+    }
+
+    labelList flipState(s.size(), FLIP);
+
+    return flipSurface(s, flipState);
+}
+
+
 // ************************************************************************* //

src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C

@@ -0,0 +1,246 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "triSurfaceRegionSearch.H"
+#include "indexedOctree.H"
+#include "triSurface.H"
+#include "PatchTools.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::triSurfaceRegionSearch::triSurfaceRegionSearch(const triSurface& surface)
+:
+    triSurfaceSearch(surface),
+    indirectRegionPatches_(),
+    treeByRegion_()
+{}
+
+
+Foam::triSurfaceRegionSearch::triSurfaceRegionSearch
+(
+    const triSurface& surface,
+    const dictionary& dict
+)
+:
+    triSurfaceSearch(surface, dict),
+    indirectRegionPatches_(),
+    treeByRegion_()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::triSurfaceRegionSearch::~triSurfaceRegionSearch()
+{
+    clearOut();
+}
+
+
+void Foam::triSurfaceRegionSearch::clearOut()
+{
+    triSurfaceSearch::clearOut();
+    treeByRegion_.clear();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::PtrList<Foam::triSurfaceRegionSearch::treeType>&
+Foam::triSurfaceRegionSearch::treeByRegion() const
+{
+    if (treeByRegion_.empty())
+    {
+        Map<label> regionSizes;
+        forAll(surface(), fI)
+        {
+            const label regionI = surface()[fI].region();
+
+            regionSizes(regionI)++;
+        }
+
+        label nRegions = regionSizes.size();
+
+        indirectRegionPatches_.setSize(nRegions);
+        treeByRegion_.setSize(nRegions);
+
+        labelListList regionsAddressing(nRegions);
+
+        forAll(regionsAddressing, regionI)
+        {
+            regionsAddressing[regionI] = labelList(regionSizes[regionI], -1);
+        }
+
+        labelList nFacesInRegions(nRegions, 0);
+
+        forAll(surface(), fI)
+        {
+            const label regionI = surface()[fI].region();
+
+            regionsAddressing[regionI][nFacesInRegions[regionI]++] = fI;
+        }
+
+        forAll(regionsAddressing, regionI)
+        {
+            scalar oldTol = treeType::perturbTol();
+            treeType::perturbTol() = tolerance();
+
+            indirectRegionPatches_.set
+            (
+                regionI,
+                new indirectTriSurface
+                (
+                    IndirectList<labelledTri>
+                    (
+                        surface(),
+                        regionsAddressing[regionI]
+                    ),
+                    surface().points()
+                )
+            );
+
+            // Calculate bb without constructing local point numbering.
+            treeBoundBox bb;
+            label nPoints;
+            PatchTools::calcBounds
+            (
+                indirectRegionPatches_[regionI],
+                bb,
+                nPoints
+            );
+
+//            if (nPoints != surface().points().size())
+//            {
+//                WarningIn("triSurfaceRegionSearch::treeByRegion() const")
+//                    << "Surface does not have compact point numbering. "
+//                    << "Of " << surface().points().size()
+//                    << " only " << nPoints
+//                    << " are used. This might give problems in some routines."
+//                    << endl;
+//            }
+
+            // Random number generator. Bit dodgy since not exactly random ;-)
+            Random rndGen(65431);
+
+            // Slightly extended bb. Slightly off-centred just so on symmetric
+            // geometry there are fewer face/edge aligned items.
+            bb = bb.extend(rndGen, 1e-4);
+            bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+            bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+            treeByRegion_.set
+            (
+                regionI,
+                new treeType
+                (
+                    treeDataIndirectTriSurface
+                    (
+                        true,
+                        indirectRegionPatches_[regionI],
+                        tolerance()
+                    ),
+                    bb,
+                    maxTreeDepth(),  // maxLevel
+                    10,              // leafsize
+                    3.0              // duplicity
+                )
+            );
+
+            treeType::perturbTol() = oldTol;
+        }
+    }
+
+    return treeByRegion_;
+}
+
+
+void Foam::triSurfaceRegionSearch::findNearest
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    const labelList& regionIndices,
+    List<pointIndexHit>& info
+) const
+{
+    if (regionIndices.empty())
+    {
+        triSurfaceSearch::findNearest(samples, nearestDistSqr, info);
+    }
+    else
+    {
+        scalar oldTol = treeType::perturbTol();
+        treeType::perturbTol() = tolerance();
+
+        const PtrList<treeType>& octrees = treeByRegion();
+
+        info.setSize(samples.size());
+
+        forAll(octrees, treeI)
+        {
+            if (findIndex(regionIndices, treeI) == -1)
+            {
+                continue;
+            }
+
+            const treeType& octree = octrees[treeI];
+
+            forAll(samples, i)
+            {
+//                if (!octree.bb().contains(samples[i]))
+//                {
+//                    continue;
+//                }
+
+                pointIndexHit currentRegionHit = octree.findNearest
+                (
+                    samples[i],
+                    nearestDistSqr[i],
+                    treeDataIndirectTriSurface::findNearestOp(octree)
+                );
+
+                if
+                (
+                    currentRegionHit.hit()
+                 &&
+                    (
+                        !info[i].hit()
+                     ||
+                        (
+                            magSqr(currentRegionHit.hitPoint() - samples[i])
+                          < magSqr(info[i].hitPoint() - samples[i])
+                        )
+                    )
+                )
+                {
+                    info[i] = currentRegionHit;
+                }
+            }
+        }
+
+        treeType::perturbTol() = oldTol;
+    }
+}
+
+
+// ************************************************************************* //

src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.H

@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::triSurfaceRegionSearch
+
+Description
+    Helper class to search on triSurface. Creates an octree for each region of
+    the surface and only searches on the specified regions.
+
+SourceFiles
+    triSurfaceRegionSearch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef triSurfaceRegionSearch_H
+#define triSurfaceRegionSearch_H
+
+#include "pointField.H"
+#include "pointIndexHit.H"
+#include "triSurfaceSearch.H"
+#include "labelledTri.H"
+#include "IndirectList.H"
+#include "PtrList.H"
+#include "indexedOctree.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class triSurfaceRegionSearch Declaration
+\*---------------------------------------------------------------------------*/
+
+class triSurfaceRegionSearch
+:
+    public triSurfaceSearch
+{
+    // Private typedefs
+
+        typedef PrimitivePatch
+        <
+            labelledTri,
+            IndirectList,
+            const pointField&
+        > indirectTriSurface;
+
+        typedef treeDataPrimitivePatch<indirectTriSurface>
+            treeDataIndirectTriSurface;
+
+        typedef indexedOctree<treeDataIndirectTriSurface> treeType;
+
+
+    // Private data
+
+        //- Surface is split into patches by region
+        mutable PtrList<indirectTriSurface> indirectRegionPatches_;
+
+        //- Search tree for each region
+        mutable PtrList<treeType> treeByRegion_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        triSurfaceRegionSearch(const triSurfaceRegionSearch&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const triSurfaceRegionSearch&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from surface. Holds reference to surface!
+        explicit triSurfaceRegionSearch(const triSurface&);
+
+        //- Construct from surface and dictionary. Holds reference to surface!
+        triSurfaceRegionSearch(const triSurface&, const dictionary& dict);
+
+
+    //- Destructor
+    ~triSurfaceRegionSearch();
+
+        //- Clear storage
+        void clearOut();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Demand driven construction of octree for each region.
+            //  @todo Currently creates a tree for each region; could optimise
+            //        by only constructing trees when they are in regionIndices
+            const PtrList<treeType>& treeByRegion() const;
+
+
+        // Query
+
+            //- Find the nearest point on the surface out of the regions
+            //  supplied in the list regionIndices. Ignores regions that are
+            //  not specified
+            void findNearest
+            (
+                const pointField& samples,
+                const scalarField& nearestDistSqr,
+                const labelList& regionIndices,
+                List<pointIndexHit>& info
+            ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C

@@ -24,63 +24,126 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "triSurfaceSearch.H"
-#include "indexedOctree.H"
-#include "boolList.H"
-#include "treeDataTriSurface.H"
 #include "triSurface.H"
-#include "line.H"
+#include "PatchTools.H"
-#include "cpuTime.H"
+#include "volumeType.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-const Foam::point Foam::triSurfaceSearch::greatPoint(GREAT, GREAT, GREAT);
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from surface. Holds reference!
 Foam::triSurfaceSearch::triSurfaceSearch(const triSurface& surface)
 :
     surface_(surface),
+    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
+    maxTreeDepth_(10),
+    treePtr_(NULL)
+{}
+
+
+Foam::triSurfaceSearch::triSurfaceSearch
+(
+    const triSurface& surface,
+    const dictionary& dict
+)
+:
+    surface_(surface),
+    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
+    maxTreeDepth_(10),
     treePtr_(NULL)
 {
-    // Random number generator. Bit dodgy since not exactly random ;-)
+    // Have optional non-standard search tolerance for gappy surfaces.
-    Random rndGen(65431);
+    if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0)
+    {
+        Info<< "    using intersection tolerance " << tolerance_ << endl;
+    }
 
-    // Slightly extended bb. Slightly off-centred just so on symmetric
+    // Have optional non-standard tree-depth to limit storage.
-    // geometry there are less face/edge aligned items.
+    if (dict.readIfPresent("maxTreeDepth", maxTreeDepth_) && maxTreeDepth_ > 0)
-    treeBoundBox treeBb
+    {
-    (
+        Info<< "    using maximum tree depth " << maxTreeDepth_ << endl;
-        treeBoundBox(surface_.points(), surface_.meshPoints()).extend
+    }
-        (
+}
-            rndGen,
-            1e-4
-        )
-    );
-    treeBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-    treeBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 
-    treePtr_.reset
+
-    (
+Foam::triSurfaceSearch::triSurfaceSearch
-        new indexedOctree<treeDataTriSurface>
+(
-        (
+    const triSurface& surface,
-            treeDataTriSurface
+    const scalar tolerance,
-            (
+    const label maxTreeDepth
-                true,
+)
-                surface_,
+:
-                indexedOctree<treeDataTriSurface>::perturbTol()
+    surface_(surface),
-            ),
+    tolerance_(tolerance),
-            treeBb,
+    maxTreeDepth_(maxTreeDepth),
-            8,      // maxLevel
+    treePtr_(NULL)
-            10,     // leafsize
+{}
-            3.0     // duplicity
+
-        )
+
-    );
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::triSurfaceSearch::~triSurfaceSearch()
+{
+    clearOut();
+}
+
+
+void Foam::triSurfaceSearch::clearOut()
+{
+    treePtr_.clear();
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+const Foam::indexedOctree<Foam::treeDataTriSurface>&
+Foam::triSurfaceSearch::tree() const
+{
+    if (treePtr_.empty())
+    {
+        // Calculate bb without constructing local point numbering.
+        treeBoundBox bb;
+        label nPoints;
+        PatchTools::calcBounds(surface(), bb, nPoints);
+
+        if (nPoints != surface().points().size())
+        {
+            WarningIn("triSurfaceSearch::tree() const")
+                << "Surface does not have compact point numbering."
+                << " Of " << surface().points().size() << " only " << nPoints
+                << " are used. This might give problems in some routines."
+                << endl;
+        }
+
+        // Random number generator. Bit dodgy since not exactly random ;-)
+        Random rndGen(65431);
+
+        // Slightly extended bb. Slightly off-centred just so on symmetric
+        // geometry there are less face/edge aligned items.
+        bb = bb.extend(rndGen, 1e-4);
+        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+        scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+        indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
+        treePtr_.reset
+        (
+            new indexedOctree<treeDataTriSurface>
+            (
+                treeDataTriSurface(true, surface_, tolerance_),
+                bb,
+                maxTreeDepth_,  // maxLevel
+                10,             // leafsize
+                3.0             // duplicity
+            )
+        );
+
+        indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
+    }
+
+    return treePtr_();
+}
+
+
 // Determine inside/outside for samples
 Foam::boolList Foam::triSurfaceSearch::calcInside
 (
@@ -97,11 +160,7 @@ Foam::boolList Foam::triSurfaceSearch::calcInside
         {
             inside[sampleI] = false;
         }
-        else if
+        else if (tree().getVolumeType(sample) == volumeType::INSIDE)
-        (
-            tree().getVolumeType(sample)
-         == indexedOctree<treeDataTriSurface>::INSIDE
-        )
         {
             inside[sampleI] = true;
         }
@@ -114,65 +173,31 @@ Foam::boolList Foam::triSurfaceSearch::calcInside
 }
 
 
-Foam::labelList Foam::triSurfaceSearch::calcNearestTri
+void Foam::triSurfaceSearch::findNearest
 (
     const pointField& samples,
-    const vector& span
+    const scalarField& nearestDistSqr,
+    List<pointIndexHit>& info
 ) const
 {
-    labelList nearest(samples.size());
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
 
-    const scalar nearestDistSqr = 0.25*magSqr(span);
+    const indexedOctree<treeDataTriSurface>& octree = tree();
 
-    pointIndexHit hitInfo;
+    info.setSize(samples.size());
 
-    forAll(samples, sampleI)
+    forAll(samples, i)
     {
-        hitInfo = tree().findNearest(samples[sampleI], nearestDistSqr);
+        static_cast<pointIndexHit&>(info[i]) = octree.findNearest
-
+        (
-        if (hitInfo.hit())
+            samples[i],
-        {
+            nearestDistSqr[i],
-            nearest[sampleI] = hitInfo.index();
+            treeDataTriSurface::findNearestOp(octree)
-        }
+        );
-        else
-        {
-            nearest[sampleI] = -1;
-        }
     }
 
-    return nearest;
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
-}
-
-
-// Nearest point on surface
-Foam::tmp<Foam::pointField> Foam::triSurfaceSearch::calcNearest
-(
-    const pointField& samples,
-    const vector& span
-) const
-{
-    const scalar nearestDistSqr = 0.25*magSqr(span);
-
-    tmp<pointField> tnearest(new pointField(samples.size()));
-    pointField& nearest = tnearest();
-
-    pointIndexHit hitInfo;
-
-    forAll(samples, sampleI)
-    {
-        hitInfo = tree().findNearest(samples[sampleI], nearestDistSqr);
-
-        if (hitInfo.hit())
-        {
-            nearest[sampleI] = hitInfo.hitPoint();
-        }
-        else
-        {
-            nearest[sampleI] = greatPoint;
-        }
-    }
-
-    return tnearest;
 }
 
 
@@ -189,84 +214,112 @@ const
 }
 
 
+void Foam::triSurfaceSearch::findLine
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    const indexedOctree<treeDataTriSurface>& octree = tree();
+
+    info.setSize(start.size());
+
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
+
+    forAll(start, i)
+    {
+        static_cast<pointIndexHit&>(info[i]) = octree.findLine
+        (
+            start[i],
+            end[i]
+        );
+    }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
+}
+
+
+void Foam::triSurfaceSearch::findLineAny
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    const indexedOctree<treeDataTriSurface>& octree = tree();
+
+    info.setSize(start.size());
+
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
+
+    forAll(start, i)
+    {
+        static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
+        (
+            start[i],
+            end[i]
+        );
+    }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
+}
+
+
 void Foam::triSurfaceSearch::findLineAll
 (
-    const point& start,
+    const pointField& start,
-    const point& end,
+    const pointField& end,
-    List<pointIndexHit>& hits
+    List<List<pointIndexHit> >& info
-)
+) const
-const
 {
-    // See if any intersection between pt and end
+    const indexedOctree<treeDataTriSurface>& octree = tree();
-    pointIndexHit inter = tree().findLine(start, end);
 
-    if (inter.hit())
+    info.setSize(start.size());
+
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
+
+    // Work array
+    DynamicList<pointIndexHit, 1, 1> hits;
+
+    DynamicList<label> shapeMask;
+
+    treeDataTriSurface::findAllIntersectOp allIntersectOp(octree, shapeMask);
+
+    forAll(start, pointI)
     {
-        hits.setSize(1);
+        hits.clear();
-        hits[0] = inter;
+        shapeMask.clear();
-
-        const vector dirVec(end-start);
-        const scalar magSqrDirVec(magSqr(dirVec));
-        const vector smallVec
-        (
-            indexedOctree<treeDataTriSurface>::perturbTol()*dirVec
-          + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
-        );
-
-
-        // Initial perturbation amount
-        vector perturbVec(smallVec);
 
         while (true)
         {
-            // Start tracking from last hit.
-            point pt = hits.last().hitPoint() + perturbVec;
-
-            if (((pt-start)&dirVec) > magSqrDirVec)
-            {
-                return;
-            }
-
             // See if any intersection between pt and end
-            pointIndexHit inter = tree().findLine(pt, end);
+            pointIndexHit inter = octree.findLine
+            (
+                start[pointI],
+                end[pointI],
+                allIntersectOp
+            );
 
-            if (!inter.hit())
+            if (inter.hit())
             {
-                return;
+                hits.append(inter);
-            }
 
-            // Check if already found this intersection
+                shapeMask.append(inter.index());
-            bool duplicateHit = false;
-            forAllReverse(hits, i)
-            {
-                if (hits[i].index() == inter.index())
-                {
-                    duplicateHit = true;
-                    break;
-                }
-            }
-
-
-            if (duplicateHit)
-            {
-                // Hit same triangle again. Increase perturbVec and try again.
-                perturbVec *= 2;
             }
             else
             {
-                // Proper hit
+                break;
-                label sz = hits.size();
-                hits.setSize(sz+1);
-                hits[sz] = inter;
-                // Restore perturbVec
-                perturbVec = smallVec;
             }
         }
+
+        info[pointI].transfer(hits);
     }
-    else
+
-    {
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
-        hits.clear();
-    }
 }
 
 

src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H

@@ -48,7 +48,6 @@ namespace Foam
 
 // Forward declaration of classes
 class triSurface;
-template<class Type> class indexedOctree;
 
 /*---------------------------------------------------------------------------*\
                            Class triSurfaceSearch Declaration
@@ -61,64 +60,82 @@ class triSurfaceSearch
         //- Reference to surface to work on
         const triSurface& surface_;
 
+        //- Optional tolerance to use in searches
+        scalar tolerance_;
+
+        //- Optional max tree depth of octree
+        label maxTreeDepth_;
+
         //- Octree for searches
-        autoPtr<indexedOctree<treeDataTriSurface> > treePtr_;
+        mutable autoPtr<indexedOctree<treeDataTriSurface> > treePtr_;
 
 
     // Private Member Functions
 
-
         //- Disallow default bitwise copy construct
         triSurfaceSearch(const triSurfaceSearch&);
 
         //- Disallow default bitwise assignment
         void operator=(const triSurfaceSearch&);
 
+
 public:
 
-    // Static data members
-
-        //- Point far away; used for illegal finds
-        static const point greatPoint;
-
-
     // Constructors
 
         //- Construct from surface. Holds reference to surface!
-        triSurfaceSearch(const triSurface&);
+        explicit triSurfaceSearch(const triSurface&);
 
+        //- Construct from surface and dictionary.
+        triSurfaceSearch(const triSurface&, const dictionary& dict);
+
+        //- Construct from components
+        triSurfaceSearch
+        (
+            const triSurface& surface,
+            const scalar tolerance,
+            const label maxTreeDepth
+        );
+
+
+    //- Destructor
+    ~triSurfaceSearch();
+
+        //- Clear storage
+        void clearOut();
 
 
     // Member Functions
 
-        const indexedOctree<treeDataTriSurface>& tree() const
+        //- Demand driven construction of the octree
-        {
+        const indexedOctree<treeDataTriSurface>& tree() const;
-            return treePtr_();
-        }
 
+        //- Return reference to the surface.
         const triSurface& surface() const
         {
             return surface_;
         }
 
+        //- Return tolerance to use in searches
+        scalar tolerance() const
+        {
+            return tolerance_;
+        }
+
+        //- Return max tree depth of octree
+        label maxTreeDepth() const
+        {
+            return maxTreeDepth_;
+        }
+
         //- Calculate for each searchPoint inside/outside status.
         boolList calcInside(const pointField& searchPoints) const;
 
-        //- Calculate index of nearest triangle (or -1) for each sample.
+        void findNearest
-        //  Looks only in box of size 2*span around sample.
-        labelList calcNearestTri
         (
             const pointField& samples,
-            const vector& span
+            const scalarField& nearestDistSqr,
-        ) const;
+            List<pointIndexHit>& info
-
-        //- Calculate nearest points (to searchPoints) on surface.
-        //  Looks only in box of size 2*span around sample. Returns greatPoint
-        //  if not found.
-        tmp<pointField> calcNearest
-        (
-            const pointField& samples,
-            const vector& span
         ) const;
 
         //- Calculate nearest point on surface for single searchPoint. Returns
@@ -128,14 +145,27 @@ public:
         //  - index()    : surface triangle label
         pointIndexHit nearest(const point&, const vector& span) const;
 
+        void findLine
+        (
+            const pointField& start,
+            const pointField& end,
+            List<pointIndexHit>& info
+        ) const;
+
+        void findLineAny
+        (
+            const pointField& start,
+            const pointField& end,
+            List<pointIndexHit>& info
+        ) const;
+
         //- Calculate all intersections from start to end
         void findLineAll
         (
-            const point& start,
+            const pointField& start,
-            const point& end,
+            const pointField& end,
-            List<pointIndexHit>&
+            List<List<pointIndexHit> >& info
         ) const;
-
 };
 
 

src/meshTools/twoDPointCorrector/twoDPointCorrector.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,7 +31,7 @@ Description
     and thus introduce a third dimension into a 2-D problem.
 
     The operation is performed by looping through all edges approximately
-    normal to the plane and enforcing their orthoginality onto the plane by
+    normal to the plane and enforcing their orthogonality onto the plane by
     adjusting points on their ends.
 
 SourceFiles

src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,6 +37,7 @@ License
 #include "geomDecomp.H"
 #include "vectorList.H"
 #include "PackedBoolList.H"
+#include "PatchTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -2356,7 +2357,7 @@ void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
 {
     boundBox bb;
     label nPoints;
-    calcBounds(bb, nPoints);
+    PatchTools::calcBounds(static_cast<const triSurface&>(*this), bb, nPoints);
     reduce(bb.min(), minOp<point>());
     reduce(bb.max(), maxOp<point>());