Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-dev-OpenCFD into develop

Conflicts:
	src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

src/meshTools/AABBTree/AABBTree.C

@@ -0,0 +1,531 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "AABBTree.H"
+#include "meshTools.H"
+#include "PackedBoolList.H"
+//#include "OFstream.H"
+
+template<class Type>
+Foam::scalar Foam::AABBTree<Type>::tolerance_ = 1e-4;
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class Type>
+void Foam::AABBTree<Type>::writeOBJ
+(
+    const bool writeLinesOnly,
+    const treeBoundBox& bb,
+    label& vertI,
+    Ostream& os
+) const
+{
+    const pointField pts(bb.points());
+    forAll(pts, i)
+    {
+        meshTools::writeOBJ(os, pts[i]);
+    }
+
+    if (writeLinesOnly)
+    {
+        forAll(bb.edges, i)
+        {
+            const edge& e = bb.edges[i];
+            os  << "l " << e[0] + vertI + 1 << ' ' << e[1] + vertI + 1 << nl;
+        }
+    }
+    else
+    {
+        forAll(bb.faces, i)
+        {
+            const face& f = bb.faces[i];
+
+            os  << 'f';
+            forAll(f, fp)
+            {
+                os  << ' ' << f[fp] + vertI + 1;
+            }
+            os  << nl;
+        }
+    }
+
+    vertI += pts.size();
+}
+
+
+template<class Type>
+void Foam::AABBTree<Type>::writeOBJ
+(
+    const bool leavesOnly,
+    const bool writeLinesOnly,
+    const treeBoundBox& bb,
+    const label nodeI,
+    const List<Pair<treeBoundBox> >& bbs,
+    const List<Pair<label> >& nodes,
+    label& vertI,
+    Ostream& os
+) const
+{
+    if (!leavesOnly || nodeI < 0)
+    {
+        writeOBJ(writeLinesOnly, bb, vertI, os);
+    }
+
+    // recurse to find leaves
+    if (nodeI >= 0)
+    {
+        writeOBJ
+        (
+            leavesOnly,
+            writeLinesOnly,
+            bbs[nodeI].first(),
+            nodes[nodeI].first(),
+            bbs,
+            nodes,
+            vertI,
+            os
+        );
+        writeOBJ
+        (
+            leavesOnly,
+            writeLinesOnly,
+            bbs[nodeI].second(),
+            nodes[nodeI].second(),
+            bbs,
+            nodes,
+            vertI,
+            os
+        );
+    }
+}
+
+
+template<class Type>
+void Foam::AABBTree<Type>::createBoxes
+(
+    const bool equalBinSize,
+    const label level,
+    const List<Type>& objects,
+    const pointField& points,
+    const DynamicList<label>& objectIDs,
+    const treeBoundBox& bb,
+    const label nodeI,
+
+    DynamicList<Pair<treeBoundBox> >& bbs,
+    DynamicList<labelPair>& nodes,
+    DynamicList<labelList>& addressing
+) const
+{
+    const vector span = bb.span();
+
+    // Determine which direction to divide the box
+
+    direction maxDir = 0;
+    scalar maxSpan = span[maxDir];
+    for (label dirI = 1; dirI < 3; dirI++)
+    {
+        if (span[dirI] > maxSpan)
+        {
+            maxSpan = span[dirI];
+            maxDir = dirI;
+        }
+    }
+
+
+    scalar divide;
+
+    if (equalBinSize)
+    {
+        // Pick up points used by this set of objects
+
+        PackedBoolList isUsedPoint(points.size());
+        DynamicList<scalar> component(points.size());
+
+        forAll(objectIDs, i)
+        {
+            const label objI = objectIDs[i];
+            const Type& obj = objects[objI];
+
+            forAll(obj, pI)
+            {
+                const label pointI = obj[pI];
+                if (isUsedPoint.set(pointI))
+                {
+                    component.append(points[pointI][maxDir]);
+                }
+            }
+        }
+
+        // Determine the median
+
+        Foam::sort(component);
+
+        divide = component[component.size()/2];
+    }
+    else
+    {
+        // Geometric middle
+        divide = bb.min()[maxDir] + 0.5*maxSpan;
+    }
+
+
+    scalar divMin = divide + tolerance_*maxSpan;
+    scalar divMax = divide - tolerance_*maxSpan;
+
+
+    // Assign the objects to min or max bin
+
+    DynamicList<label> minBinObjectIDs(objectIDs.size());
+    treeBoundBox minBb(point::max, point::min);
+
+    DynamicList<label> maxBinObjectIDs(objectIDs.size());
+    treeBoundBox maxBb(point::max, point::min);
+
+    forAll(objectIDs, i)
+    {
+        const label objI = objectIDs[i];
+        const Type& obj = objects[objI];
+
+        bool intoMin = false;
+        bool intoMax = false;
+
+        forAll(obj, pI)
+        {
+            const label pointI = obj[pI];
+            const point& pt = points[pointI];
+            if (pt[maxDir] < divMin)
+            {
+                intoMin = true;
+            }
+            if (pt[maxDir] > divMax)
+            {
+                intoMax = true;
+            }
+        }
+
+
+        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());
+        }
+        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());
+        }
+    }
+
+    // inflate box in case geometry reduces to 2-D
+    if (minBinObjectIDs.size())
+    {
+        minBb.inflate(0.01);
+    }
+    if (maxBinObjectIDs.size())
+    {
+        maxBb.inflate(0.01);
+    }
+
+    minBinObjectIDs.shrink();
+    maxBinObjectIDs.shrink();
+
+
+    label minI;
+    if (minBinObjectIDs.size() > minLeafSize_ && level < maxLevel_)
+    {
+        // new leaf
+        minI = nodes.size();
+        nodes.append(labelPair(-1, -1));
+    }
+    else
+    {
+        // update existing leaf
+        minI = -addressing.size() - 1;
+        addressing.append(minBinObjectIDs);
+    }
+
+    label maxI;
+    if (maxBinObjectIDs.size() > minLeafSize_ && level < maxLevel_)
+    {
+        // new leaf
+        maxI = nodes.size();
+        nodes.append(labelPair(-1, -1));
+    }
+    else
+    {
+        // update existing leaf
+        maxI = -addressing.size() - 1;
+        addressing.append(maxBinObjectIDs);
+    }
+
+    nodes(nodeI) = labelPair(minI, maxI);
+    bbs(nodeI) = Pair<treeBoundBox>(minBb, maxBb);
+
+    // recurse
+    if (minI >= 0)
+    {
+        createBoxes
+        (
+            equalBinSize,
+            level + 1,
+            objects,
+            points,
+            minBinObjectIDs,
+            minBb,
+            minI,
+            bbs,
+            nodes,
+            addressing
+        );
+    }
+    if (maxI >= 0)
+    {
+        createBoxes
+        (
+            equalBinSize,
+            level + 1,
+            objects,
+            points,
+            maxBinObjectIDs,
+            maxBb,
+            maxI,
+            bbs,
+            nodes,
+            addressing
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::AABBTree<Type>::AABBTree()
+:
+    maxLevel_(0),
+    minLeafSize_(0),
+    boundBoxes_(),
+    addressing_()
+{}
+
+
+template<class Type>
+Foam::AABBTree<Type>::AABBTree
+(
+    const UList<Type>& objects,
+    const pointField& points,
+    const bool equalBinSize,
+    const label maxLevel,
+    const label minLeafSize
+)
+:
+    maxLevel_(maxLevel),
+    minLeafSize_(minLeafSize),
+    boundBoxes_(),
+    addressing_()
+{
+    if (objects.empty())
+    {
+        return;
+    }
+
+
+    DynamicList<Pair<treeBoundBox> > bbs(maxLevel);
+    DynamicList<labelPair> nodes(maxLevel);
+    DynamicList<labelList> addr(maxLevel);
+
+    nodes.append(labelPair(-1, -1));
+    treeBoundBox topBb(points);
+    topBb.inflate(0.01);
+
+    DynamicList<label> objectIDs(identity(objects.size()));
+
+    createBoxes
+    (
+        equalBinSize,
+        0,          // starting at top level
+        objects,
+        points,
+        objectIDs,
+        topBb,
+        0,          // starting node
+
+        bbs,
+        nodes,
+        addr
+    );
+
+
+    //{
+    //    OFstream os("tree.obj");
+    //    label vertI = 0;
+    //    writeOBJ
+    //    (
+    //        true,       // leavesOnly
+    //        false,      // writeLinesOnly
+    //
+    //        topBb,
+    //        0,
+    //        bbs,
+    //        nodes,
+    //        vertI,
+    //        os
+    //    );
+    //}
+
+
+    // transfer flattened tree to persistent storage
+    DynamicList<treeBoundBox> boundBoxes(2*bbs.size());
+    DynamicList<labelList> addressing(2*addr.size());
+    forAll(nodes, nodeI)
+    {
+        if (nodes[nodeI].first() < 0)
+        {
+            boundBoxes.append(bbs[nodeI].first());
+            addressing.append(addr[nodeI + 1]);
+        }
+        if (nodes[nodeI].second() < 0)
+        {
+            boundBoxes.append(bbs[nodeI].second());
+            addressing.append(addr[nodeI + 1]);
+        }
+    }
+
+    boundBoxes_.transfer(boundBoxes);
+    addressing_.transfer(addressing);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+const Foam::List<Foam::treeBoundBox>& Foam::AABBTree<Type>::boundBoxes() const
+{
+    return boundBoxes_;
+}
+
+
+template<class Type>
+const Foam::List<Foam::labelList>& Foam::AABBTree<Type>::addressing() const
+{
+    return addressing_;
+}
+
+
+template<class Type>
+bool Foam::AABBTree<Type>::pointInside(const point& pt) const
+{
+    forAll(boundBoxes_, i)
+    {
+        const treeBoundBox& bb = boundBoxes_[i];
+
+        if (bb.contains(pt))
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+template<class Type>
+bool Foam::AABBTree<Type>::overlaps(const boundBox& bbIn) const
+{
+    forAll(boundBoxes_, i)
+    {
+        const treeBoundBox& bb = boundBoxes_[i];
+
+        if (bb.overlaps(bbIn))
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+// * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Ostream& Foam::operator<<(Ostream& os, const AABBTree<Type>& tree)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << tree.maxLevel_ << token::SPACE
+            << tree.minLeafSize_ << token::SPACE
+            << tree.boundBoxes_ << token::SPACE
+            << tree.addressing_ << token::SPACE;
+    }
+    else
+    {
+        os.write
+        (
+            reinterpret_cast<const char*>(&tree.maxLevel_),
+            sizeof(tree.maxLevel_)
+          + sizeof(tree.minLeafSize_)
+        );
+        os  << tree.boundBoxes_
+            << tree.addressing_;
+    }
+
+    os.check("Ostream& operator<<(Ostream&, const AABBTree<Type>&)");
+    return os;
+}
+
+
+template<class Type>
+Foam::Istream& Foam::operator>>(Istream& is, AABBTree<Type>& tree)
+{
+    if (is.format() == IOstream::ASCII)
+    {
+        is  >> tree.maxLevel_
+            >> tree.minLeafSize_
+            >> tree.boundBoxes_
+            >> tree.addressing_;
+    }
+    else
+    {
+        is.read
+        (
+            reinterpret_cast<char*>(&tree.maxLevel_),
+            sizeof(tree.maxLevel_)
+          + sizeof(tree.minLeafSize_)
+        );
+        is  >> tree.boundBoxes_
+            >> tree.addressing_;
+    }
+
+    is.check("Istream& operator>>(Istream&, AABBTree<Type>&)");
+    return is;
+}
+
+
+// ************************************************************************* //

src/meshTools/AABBTree/AABBTree.H

@@ -0,0 +1,197 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::AABBTree
+
+Description
+    Templated tree of axis-aligned bounding boxes (AABB)
+
+    Designed to be templated on either faces or cells, the AABBTree will
+    decompose the input into a tree of AABB's.  The maximum number of tree
+    levels and minimum number of objects per leaf are provided on construction,
+    and the contents (addressing) is stored.
+
+SourceFiles
+    AABBTree.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef AABBTree_H
+#define AABBTree_H
+
+#include "labelList.H"
+#include "labelPair.H"
+#include "DynamicList.H"
+#include "pointField.H"
+#include "treeBoundBox.H"
+#include "Ostream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class Type>
+class AABBTree;
+
+template<class Type>
+Istream& operator>>(Istream&, AABBTree<Type>&);
+
+template<class Type>
+Ostream& operator<<(Ostream&, const AABBTree<Type>&);
+
+/*---------------------------------------------------------------------------*\
+                          Class AABBTree Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class AABBTree
+{
+
+protected:
+
+    // Protected Data
+
+        //- Tolerance
+        static scalar tolerance_;
+
+        //- Maximum tree level
+        label maxLevel_;
+
+        //- Minimum points per leaf
+        label minLeafSize_;
+
+        //- Bounding boxes making up the tree
+        List<treeBoundBox> boundBoxes_;
+
+        //- Leaf adressing
+        List<labelList> addressing_;
+
+
+    // Protected Member Functions
+
+        //- Write OBJ file of bounding box
+        void writeOBJ
+        (
+            const bool writeLinesOnly,
+            const treeBoundBox& bb,
+            label& vertI,
+            Ostream& os
+        ) const;
+
+        //- Write OBJ for all bounding boxes
+        void writeOBJ
+        (
+            const bool leavesOnly,
+            const bool writeLinesOnly,
+            const treeBoundBox& bb,
+            const label nodeI,
+            const List<Pair<treeBoundBox> >& bbs,
+            const List<Pair<label> >& nodes,
+            label& vertI,
+            Ostream& os
+        ) const;
+
+        //- Create the bounding boxes by interrogating points
+        void createBoxes
+        (
+            const bool equalBinSize,
+            const label level,
+            const List<Type>& objects,
+            const pointField& points,
+            const DynamicList<label>& objectIDs,
+            const treeBoundBox& bb,
+            const label nodeI,
+
+            DynamicList<Pair<treeBoundBox> >& bbs,
+            DynamicList<labelPair>& nodes,
+            DynamicList<labelList>& addressing
+        ) const;
+
+
+public:
+
+    // Constructors
+
+        //- Null constructor
+        AABBTree();
+
+        //- Construct from components
+        //  equalBinSize: divide into equal number of elements or
+        //                equal span
+        AABBTree
+        (
+            const UList<Type>& objects,
+            const pointField& points,
+            const bool equalBinSize = true,
+            const label maxLevel = 3,
+            const label minBinSize = 100
+        );
+
+
+    // Public Member Functions
+
+        // Access
+
+            //- Return the bounding boxes making up the tree
+            const List<treeBoundBox>& boundBoxes() const;
+
+            //- Return the contents addressing
+            const List<labelList>& addressing() const;
+
+
+        // Evaluation
+
+            //- Determine whether a point is inside the bounding boxes
+            bool pointInside(const point& pt) const;
+
+            //- Determine whether a bounding box overlaps the tree bounding
+            //  boxes
+            bool overlaps(const boundBox& bbIn) const;
+
+
+    // IOstream operators
+
+        friend Istream& operator>> <Type>(Istream&, AABBTree&);
+        friend Ostream& operator<< <Type>(Ostream&, const AABBTree&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "AABBTree.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,6 +27,7 @@ License
 #include "AMIMethod.H"
 #include "meshTools.H"
 #include "mapDistribute.H"
+#include "flipOp.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1050,27 +1051,33 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
         // send data back to originating procs. Note that contributions
         // from different processors get added (ListAppendEqOp)
 
-        mapDistribute::distribute
+        mapDistributeBase::distribute
         (
             Pstream::nonBlocking,
             List<labelPair>(),
             tgtPatch.size(),
             map.constructMap(),
+            false,                      // has flip
             map.subMap(),
+            false,                      // has flip
             tgtAddress_,
             ListAppendEqOp<label>(),
+            flipOp(),                   // flip operation
             labelList()
         );
 
-        mapDistribute::distribute
+        mapDistributeBase::distribute
         (
             Pstream::nonBlocking,
             List<labelPair>(),
             tgtPatch.size(),
             map.constructMap(),
+            false,
             map.subMap(),
+            false,
             tgtWeights_,
             ListAppendEqOp<scalar>(),
+            flipOp(),
             scalarList()
         );
 

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationParallelOps.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,6 +26,7 @@ License
 #include "AMIInterpolation.H"
 #include "mergePoints.H"
 #include "mapDistribute.H"
+#include "AABBTree.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -97,11 +98,11 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcOverlappingProcs
 
     forAll(procBb, procI)
     {
-        const List<treeBoundBox>& bbs = procBb[procI];
+        const treeBoundBoxList& bbp = procBb[procI];
 
-        forAll(bbs, bbI)
+        forAll(bbp, bbI)
         {
-            if (bbs[bbI].overlaps(bb))
+            if (bbp[bbI].overlaps(bb))
             {
                 overlaps[procI] = true;
                 nOverlaps++;
@@ -109,6 +110,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcOverlappingProcs
             }
         }
     }
+
     return nOverlaps;
 }
 
@@ -338,33 +340,22 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcProcMap
 
     if (srcPatch.size())
     {
-        procBb[Pstream::myProcNo()] = treeBoundBoxList
-        (
-            1,  // For now single bounding box per proc
-            treeBoundBox
+        procBb[Pstream::myProcNo()] =
+            AABBTree<face>
             (
-                srcPatch.points(),
-                srcPatch.meshPoints()
-            )
-        );
+                srcPatch.localFaces(),
+                srcPatch.localPoints(),
+                false
+            ).boundBoxes();
     }
     else
     {
         procBb[Pstream::myProcNo()] = treeBoundBoxList();
     }
 
-    // slightly increase size of bounding boxes to allow for cases where
-    // bounding boxes are perfectly alligned
-    forAll(procBb[Pstream::myProcNo()], bbI)
-    {
-        treeBoundBox& bb = procBb[Pstream::myProcNo()][bbI];
-        bb.inflate(0.01);
-    }
-
     Pstream::gatherList(procBb);
     Pstream::scatterList(procBb);
 
-
     if (debug)
     {
         Info<< "Determining extent of srcPatch per processor:" << nl
@@ -375,7 +366,6 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcProcMap
         }
     }
 
-
     // Determine which faces of tgtPatch overlaps srcPatch per proc
     const faceList& faces = tgtPatch.localFaces();
     const pointField& points = tgtPatch.localPoints();

src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -221,7 +221,7 @@ void Foam::AMIMethod<SourcePatch, TargetPatch>::resetTree()
     // Clear the old octree
     treePtr_.clear();
 
-    treeBoundBox bb(tgtPatch_.points());
+    treeBoundBox bb(tgtPatch_.points(), tgtPatch_.meshPoints());
     bb.inflate(0.01);
 
     if (!treePtr_.valid())

src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -310,16 +310,14 @@ void Foam::directAMI<SourcePatch, TargetPatch>::calculate
     forAll(srcAddr, i)
     {
         scalar magSf = this->srcMagSf_[i];
-//        srcWeights[i] = scalarList(srcAddr[i].size(), magSf);
-        srcWeights[i] = scalarList(1, magSf);
         srcAddress[i].transfer(srcAddr[i]);
+        srcWeights[i] = scalarList(1, magSf);
     }
     forAll(tgtAddr, i)
     {
         scalar magSf = this->tgtMagSf_[i];
-//        tgtWeights[i] = scalarList(tgtAddr[i].size(), magSf);
-        tgtWeights[i] = scalarList(1, magSf);
         tgtAddress[i].transfer(tgtAddr[i]);
+        tgtWeights[i] = scalarList(1, magSf);
     }
 }
 

src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,6 +41,7 @@ SourceFiles
 #include "FixedList.H"
 #include "plane.H"
 #include "face.H"
+#include "triPoints.H"
 #include "NamedEnum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -56,8 +57,6 @@ class faceAreaIntersect
 {
 public:
 
-    typedef FixedList<point, 3> triPoints;
-
     enum triangulationMode
     {
         tmFan,

src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.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) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,7 +41,7 @@ inline void Foam::faceAreaIntersect::setTriPoints
 }
 
 
-inline Foam::faceAreaIntersect::triPoints Foam::faceAreaIntersect::getTriPoints
+inline Foam::triPoints Foam::faceAreaIntersect::getTriPoints
 (
     const pointField& points,
     const face& f,

src/meshTools/Make/files

@@ -63,7 +63,9 @@ indexedOctree/treeDataTriSurface.C
 
 searchableSurface = searchableSurface
 $(searchableSurface)/searchableBox.C
+$(searchableSurface)/searchableRotatedBox.C
 $(searchableSurface)/searchableCylinder.C
+$(searchableSurface)/searchableCone.C
 $(searchableSurface)/searchableDisk.C
 $(searchableSurface)/searchablePlane.C
 $(searchableSurface)/searchablePlate.C
@@ -75,6 +77,7 @@ $(searchableSurface)/searchableSurfacesQueries.C
 $(searchableSurface)/searchableSurfaceWithGaps.C
 $(searchableSurface)/triSurfaceMesh.C
 $(searchableSurface)/closedTriSurfaceMesh.C
+$(searchableSurface)/subTriSurfaceMesh.C
 
 topoSets = sets/topoSets
 $(topoSets)/cellSet.C

src/meshTools/algorithms/MeshWave/FaceCellWave.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -810,6 +810,79 @@ void Foam::FaceCellWave<Type, TrackingData>::handleAMICyclicPatches()
 }
 
 
+template<class Type, class TrackingData>
+void Foam::FaceCellWave<Type, TrackingData>::handleExplicitConnections()
+{
+    // Collect changed information
+
+    DynamicList<label> f0Baffle(explicitConnections_.size());
+    DynamicList<Type> f0Info(explicitConnections_.size());
+
+    DynamicList<label> f1Baffle(explicitConnections_.size());
+    DynamicList<Type> f1Info(explicitConnections_.size());
+
+    forAll(explicitConnections_, connI)
+    {
+        const labelPair& baffle = explicitConnections_[connI];
+
+        label f0 = baffle[0];
+        if (changedFace_[f0])
+        {
+            f0Baffle.append(connI);
+            f0Info.append(allFaceInfo_[f0]);
+        }
+
+        label f1 = baffle[1];
+        if (changedFace_[f1])
+        {
+            f1Baffle.append(connI);
+            f1Info.append(allFaceInfo_[f1]);
+        }
+    }
+
+
+    // Update other side with changed information
+
+    forAll(f1Info, index)
+    {
+        const labelPair& baffle = explicitConnections_[f1Baffle[index]];
+
+        label f0 = baffle[0];
+        Type& currentWallInfo = allFaceInfo_[f0];
+
+        if (!currentWallInfo.equal(f1Info[index], td_))
+        {
+            updateFace
+            (
+                f0,
+                f1Info[index],
+                propagationTol_,
+                currentWallInfo
+            );
+        }
+    }
+
+    forAll(f0Info, index)
+    {
+        const labelPair& baffle = explicitConnections_[f0Baffle[index]];
+
+        label f1 = baffle[1];
+        Type& currentWallInfo = allFaceInfo_[f1];
+
+        if (!currentWallInfo.equal(f0Info[index], td_))
+        {
+            updateFace
+            (
+                f1,
+                f0Info[index],
+                propagationTol_,
+                currentWallInfo
+            );
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Set up only. Use setFaceInfo and iterate() to do actual calculation.
@@ -823,6 +896,7 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
 )
 :
     mesh_(mesh),
+    explicitConnections_(0),
     allFaceInfo_(allFaceInfo),
     allCellInfo_(allCellInfo),
     td_(td),
@@ -878,6 +952,7 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
 )
 :
     mesh_(mesh),
+    explicitConnections_(0),
     allFaceInfo_(allFaceInfo),
     allCellInfo_(allCellInfo),
     td_(td),
@@ -939,6 +1014,87 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
 }
 
 
+// Iterate, propagating changedFacesInfo across mesh, until no change (or
+// maxIter reached). Initial cell values specified.
+template<class Type, class TrackingData>
+Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
+(
+    const polyMesh& mesh,
+    const List<labelPair>& explicitConnections,
+    const bool handleCyclicAMI,
+    const labelList& changedFaces,
+    const List<Type>& changedFacesInfo,
+    UList<Type>& allFaceInfo,
+    UList<Type>& allCellInfo,
+    const label maxIter,
+    TrackingData& td
+)
+:
+    mesh_(mesh),
+    explicitConnections_(explicitConnections),
+    allFaceInfo_(allFaceInfo),
+    allCellInfo_(allCellInfo),
+    td_(td),
+    changedFace_(mesh_.nFaces(), false),
+    changedFaces_(mesh_.nFaces()),
+    nChangedFaces_(0),
+    changedCell_(mesh_.nCells(), false),
+    changedCells_(mesh_.nCells()),
+    nChangedCells_(0),
+    hasCyclicPatches_(hasPatch<cyclicPolyPatch>()),
+    hasCyclicAMIPatches_
+    (
+        handleCyclicAMI
+     && returnReduce(hasPatch<cyclicAMIPolyPatch>(), orOp<bool>())
+    ),
+    nEvals_(0),
+    nUnvisitedCells_(mesh_.nCells()),
+    nUnvisitedFaces_(mesh_.nFaces())
+{
+    if
+    (
+        allFaceInfo.size() != mesh_.nFaces()
+     || allCellInfo.size() != mesh_.nCells()
+    )
+    {
+        FatalErrorIn
+        (
+            "FaceCellWave<Type, TrackingData>::FaceCellWave"
+            "(const polyMesh&, const List<labelPair>&, const labelList&"
+            ", const List<Type>,"
+            " UList<Type>&, UList<Type>&, const label maxIter)"
+        )   << "face and cell storage not the size of mesh faces, cells:"
+            << endl
+            << "    allFaceInfo   :" << allFaceInfo.size() << endl
+            << "    mesh_.nFaces():" << mesh_.nFaces() << endl
+            << "    allCellInfo   :" << allCellInfo.size() << endl
+            << "    mesh_.nCells():" << mesh_.nCells()
+            << exit(FatalError);
+    }
+
+    // Copy initial changed faces data
+    setFaceInfo(changedFaces, changedFacesInfo);
+
+    // Iterate until nothing changes
+    label iter = iterate(maxIter);
+
+    if ((maxIter > 0) && (iter >= maxIter))
+    {
+        FatalErrorIn
+        (
+            "FaceCellWave<Type, TrackingData>::FaceCellWave"
+            "(const polyMesh&, const List<labelPair>&, const labelList&"
+            ", const List<Type>, UList<Type>&, UList<Type>&"
+            ", const label maxIter)"
+        )   << "Maximum number of iterations reached. Increase maxIter." << endl
+            << "    maxIter:" << maxIter << endl
+            << "    nChangedCells:" << nChangedCells_ << endl
+            << "    nChangedFaces:" << nChangedFaces_ << endl
+            << exit(FatalError);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
@@ -1093,6 +1249,10 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::cellToFace()
     // Handled all changed cells by now
     nChangedCells_ = 0;
 
+
+    // Transfer across any explicitly provided internal connections
+    handleExplicitConnections();
+
     if (hasCyclicPatches_)
     {
         // Transfer changed faces across cyclic halves

src/meshTools/algorithms/MeshWave/FaceCellWave.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,6 +50,7 @@ SourceFiles
 #include "boolList.H"
 #include "labelList.H"
 #include "primitiveFieldsFwd.H"
+#include "labelPair.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -81,6 +82,9 @@ class FaceCellWave
         //- Reference to mesh
         const polyMesh& mesh_;
 
+        //- Optional boundary faces that information should travel through
+        const List<labelPair> explicitConnections_;
+
         //- Information for all faces
         UList<Type>& allFaceInfo_;
 
@@ -232,6 +236,10 @@ class FaceCellWave
             //- Merge data from across AMI cyclics
             void handleAMICyclicPatches();
 
+            //- Merge data across explicitly provided local connections (usually
+            //  baffles)
+            void handleExplicitConnections();
+
 
       // Private static data
 
@@ -286,6 +294,23 @@ public:
             TrackingData& td = dummyTrackData_
         );
 
+        //- Construct from mesh and explicitly connected boundary faces
+        //  and list of changed faces with the Type
+        //  for these faces. Iterates until nothing changes or maxIter reached.
+        //  (maxIter can be 0)
+        FaceCellWave
+        (
+            const polyMesh&,
+            const List<labelPair>& explicitConnections,
+            const bool handleCyclicAMI,
+            const labelList& initialChangedFaces,
+            const List<Type>& changedFacesInfo,
+            UList<Type>& allFaceInfo,
+            UList<Type>& allCellInfo,
+            const label maxIter,
+            TrackingData& td = dummyTrackData_
+        );
+
 
     // Member Functions
 

src/meshTools/algorithms/PointEdgeWave/PointData.C

@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+-------------------------------------------------------------------------------
+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 "PointData.H"
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+template<class DataType>
+Foam::Ostream& Foam::operator<<(Ostream& os, const PointData<DataType>& pd)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        return os
+            << static_cast<const pointEdgePoint&>(pd)
+            << token::SPACE << pd.data();
+    }
+    else
+    {
+        return os
+            << static_cast<const pointEdgePoint&>(pd)
+            << pd.data();
+    }
+}
+
+
+template<class DataType>
+Foam::Istream& Foam::operator>>(Istream& is, PointData<DataType>& pd)
+{
+    return is >> static_cast<pointEdgePoint&>(pd) >> pd.data_;
+}
+
+
+// ************************************************************************* //

src/meshTools/algorithms/PointEdgeWave/PointData.H

@@ -0,0 +1,211 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+-------------------------------------------------------------------------------
+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::PointData
+
+Description
+    Variant of pointEdgePoint with some transported additional data. Templated
+    on the transported data type.
+
+SourceFiles
+    PointDataI.H
+    PointData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PointData_H
+#define PointData_H
+
+#include "pointEdgePoint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+template<class DataType>
+class PointData;
+
+// Forward declaration of friend functions and operators
+template<class DataType>
+Ostream& operator<<(Ostream&, const PointData<DataType>&);
+template<class DataType>
+Istream& operator>>(Istream&, PointData<DataType>&);
+
+/*---------------------------------------------------------------------------*\
+                           Class PointData Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class DataType>
+class PointData
+:
+    public pointEdgePoint
+{
+private:
+
+    // Private data
+
+        //- Additional transported data
+        DataType data_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline PointData();
+
+        //- Construct from origin, distance and data
+        inline PointData
+        (
+            const point& origin,
+            const scalar distSqr,
+            const DataType& data
+        );
+
+        //- Construct as copy
+        inline PointData(const PointData&);
+
+
+    // Member Functions
+
+        // Access
+
+            //- Const access the data
+            inline const DataType& data() const;
+
+
+        // Needed by meshWave
+
+            //- Apply rotation matrix to the data
+            template<class TrackingData>
+            inline void transform
+            (
+                const tensor& rotTensor,
+                TrackingData& td
+            );
+
+            //- Influence of edge on point
+            template<class TrackingData>
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const label edgeI,
+                const PointData<DataType>& edgeInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of different value on same point.
+            //  Merge new and old info.
+            template<class TrackingData>
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const PointData<DataType>& newPointInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of different value on same point.
+            //  No information about current position whatsoever.
+            template<class TrackingData>
+            inline bool updatePoint
+            (
+                const PointData<DataType>& newPointInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of point on edge.
+            template<class TrackingData>
+            inline bool updateEdge
+            (
+                const polyMesh& mesh,
+                const label edgeI,
+                const label pointI,
+                const PointData<DataType>& pointInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+    // Member Operators
+
+        inline bool operator==(const PointData<DataType>&) const;
+        inline bool operator!=(const PointData<DataType>&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<< <DataType>
+        (
+            Ostream&,
+            const PointData<DataType>&
+        );
+        friend Istream& operator>> <DataType>
+        (
+            Istream&,
+            PointData<DataType>&
+        );
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Data associated with PointData types is contiguous
+
+template<>
+inline bool contiguous<PointData<scalar> >()
+{
+    return true;
+}
+
+template<>
+inline bool contiguous<PointData<vector> >()
+{
+    return true;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "PointDataI.H"
+
+#ifdef NoRepository
+#   include "PointData.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/algorithms/PointEdgeWave/PointDataI.H

@@ -0,0 +1,234 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+-------------------------------------------------------------------------------
+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 "polyMesh.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class DataType>
+inline Foam::PointData<DataType>::PointData()
+:
+    pointEdgePoint()
+{}
+
+
+template<class DataType>
+inline Foam::PointData<DataType>::PointData
+(
+    const point& origin,
+    const scalar distSqr,
+    const DataType& data
+)
+:
+    pointEdgePoint(origin, distSqr),
+    data_(data)
+{}
+
+
+template<class DataType>
+inline Foam::PointData<DataType>::PointData(const PointData<DataType>& wpt)
+:
+    pointEdgePoint(wpt),
+    data_(wpt.data())
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class DataType>
+inline const DataType& Foam::PointData<DataType>::data() const
+{
+    return data_;
+}
+
+
+template<class DataType>
+template<class TrackingData>
+inline void Foam::PointData<DataType>::transform
+(
+    const tensor& rotTensor,
+    TrackingData& td
+)
+{
+    pointEdgePoint::transform(rotTensor, td);
+    data_ = Foam::transform(rotTensor, data_);
+}
+
+
+template<class DataType>
+template<class TrackingData>
+inline bool Foam::PointData<DataType>::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const label edgeI,
+    const PointData<DataType>& edgeInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if
+    (
+        pointEdgePoint::updatePoint
+        (
+            mesh,
+            pointI,
+            edgeI,
+            edgeInfo,
+            tol,
+            td
+        )
+    )
+    {
+        data_ = edgeInfo.data_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class DataType>
+template<class TrackingData>
+inline bool Foam::PointData<DataType>::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const PointData<DataType>& newPointInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if
+    (
+        pointEdgePoint::updatePoint
+        (
+            mesh,
+            pointI,
+            newPointInfo,
+            tol,
+            td
+        )
+    )
+    {
+        data_ = newPointInfo.data_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class DataType>
+template<class TrackingData>
+inline bool Foam::PointData<DataType>::updatePoint
+(
+    const PointData<DataType>& newPointInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (pointEdgePoint::updatePoint(newPointInfo, tol, td))
+    {
+        data_ = newPointInfo.data_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class DataType>
+template<class TrackingData>
+inline bool Foam::PointData<DataType>::updateEdge
+(
+    const polyMesh& mesh,
+    const label edgeI,
+    const label pointI,
+    const PointData<DataType>& pointInfo,
+    const scalar tol,
+    TrackingData& td
+
+)
+{
+    if
+    (
+        pointEdgePoint::updateEdge
+        (
+            mesh,
+            edgeI,
+            pointI,
+            pointInfo,
+            tol,
+            td
+        )
+    )
+    {
+        data_ = pointInfo.data_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class DataType>
+inline bool Foam::PointData<DataType>::operator==
+(
+    const Foam::PointData<DataType>& rhs
+)
+const
+{
+    return pointEdgePoint::operator==(rhs) && (data() == rhs.data());
+}
+
+
+template<class DataType>
+inline bool Foam::PointData<DataType>::operator!=
+(
+    const Foam::PointData<DataType>& rhs
+)
+const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //

src/meshTools/indexedOctree/treeDataPrimitivePatch.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -188,8 +188,10 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
             << " nearest face:" << faceI;
     }
 
-    const pointField& points = patch_.localPoints();
-    const typename PatchType::FaceType& f = patch_.localFaces()[faceI];
+    const typename PatchType::FaceType& localF = patch_.localFaces()[faceI];
+    const typename PatchType::FaceType& f = patch_[faceI];
+    const pointField& points = patch_.points();
+    const labelList& mp = patch_.meshPoints();
 
     // Retest to classify where on face info is. Note: could be improved. We
     // already have point.
@@ -242,7 +244,7 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
 
             return indexedOctree<treeDataPrimitivePatch>::getSide
             (
-                patch_.pointNormals()[f[fp]],
+                patch_.pointNormals()[localF[fp]],
                 sample - curPt
             );
         }
@@ -280,8 +282,8 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
     {
         label edgeI = fEdges[fEdgeI];
         const edge& e = patch_.edges()[edgeI];
-
-        pointHit edgeHit = e.line(points).nearestDist(sample);
+        const linePointRef ln(points[mp[e.start()]], points[mp[e.end()]]);
+        pointHit edgeHit = ln.nearestDist(sample);
 
         if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < planarTol_)
         {
@@ -322,11 +324,7 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
 
     forAll(f, fp)
     {
-        pointHit edgeHit = linePointRef
-        (
-            points[f[fp]],
-            fc
-        ).nearestDist(sample);
+        pointHit edgeHit = linePointRef(points[f[fp]], fc).nearestDist(sample);
 
         if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < planarTol_)
         {
@@ -369,7 +367,8 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
 
         forAll(f, fp)
         {
-            Pout<< "    vertex:" << f[fp] << "  coord:" << points[f[fp]]
+            Pout<< "    vertex:" << f[fp]
+                << "  coord:" << points[f[fp]]
                 << endl;
         }
     }

src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.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) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -837,26 +837,28 @@ void Foam::mappedPatchBase::calcAMI() const
 
     AMIPtr_.clear();
 
+
+    const polyPatch& nbr = samplePolyPatch();
+
+    // Transform neighbour patch to local system
+    pointField nbrPoints(samplePoints(nbr.localPoints()));
+
+    primitivePatch nbrPatch0
+    (
+        SubList<face>
+        (
+            nbr.localFaces(),
+            nbr.size()
+        ),
+        nbrPoints
+    );
+
+
     if (debug)
     {
-        const polyPatch& nbr = samplePolyPatch();
-
-        pointField nbrPoints(nbr.localPoints());
-
         OFstream os(patch_.name() + "_neighbourPatch-org.obj");
         meshTools::writeOBJ(os, samplePolyPatch().localFaces(), nbrPoints);
 
-        // transform neighbour patch to local system
-        primitivePatch nbrPatch0
-        (
-            SubList<face>
-            (
-                nbr.localFaces(),
-                nbr.size()
-            ),
-            nbrPoints
-        );
-
         OFstream osN(patch_.name() + "_neighbourPatch-trans.obj");
         meshTools::writeOBJ(osN, nbrPatch0, nbrPoints);
 
@@ -870,7 +872,7 @@ void Foam::mappedPatchBase::calcAMI() const
         new AMIPatchToPatchInterpolation
         (
             patch_,
-            samplePolyPatch(), // nbrPatch0,
+            nbrPatch0,
             surfPtr(),
             faceAreaIntersect::tmMesh,
             true,

src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBaseTemplates.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,15 +61,18 @@ void Foam::mappedPatchBase::distribute
         }
         default:
         {
-            map().distribute
+            mapDistributeBase::distribute
             (
                 Pstream::defaultCommsType,
                 map().schedule(),
                 map().constructSize(),
                 map().subMap(),
+                false,
                 map().constructMap(),
+                false,
                 lst,
                 cop,
+                flipOp(),
                 pTraits<Type>::zero
             );
         }
@@ -117,15 +120,18 @@ void Foam::mappedPatchBase::reverseDistribute
         default:
         {
             label cSize = sampleSize();
-            map().distribute
+            mapDistributeBase::distribute
             (
                 Pstream::defaultCommsType,
                 map().schedule(),
                 cSize,
                 map().constructMap(),
+                false,
                 map().subMap(),
+                false,
                 lst,
                 cop,
+                flipOp(),
                 pTraits<Type>::zero
             );
             break;

src/meshTools/meshTools/meshTools.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,7 +29,6 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// Check if n is in same direction as normals of all faceLabels
 bool Foam::meshTools::visNormal
 (
     const vector& n,
@@ -309,7 +308,6 @@ bool Foam::meshTools::edgeOnFace
 }
 
 
-// Return true if faceI part of cellI
 bool Foam::meshTools::faceOnCell
 (
     const primitiveMesh& mesh,
@@ -465,7 +463,6 @@ Foam::label Foam::meshTools::getSharedFace
 }
 
 
-// Get the two faces on cellI using edgeI.
 void Foam::meshTools::getEdgeFaces
 (
     const primitiveMesh& mesh,
@@ -511,7 +508,6 @@ void Foam::meshTools::getEdgeFaces
 }
 
 
-// Return label of other edge connected to vertex
 Foam::label Foam::meshTools::otherEdge
 (
     const primitiveMesh& mesh,
@@ -549,7 +545,6 @@ Foam::label Foam::meshTools::otherEdge
 }
 
 
-// Return face on other side of edgeI
 Foam::label Foam::meshTools::otherFace
 (
     const primitiveMesh& mesh,
@@ -574,7 +569,6 @@ Foam::label Foam::meshTools::otherFace
 }
 
 
-// Return face on other side of edgeI
 Foam::label Foam::meshTools::otherCell
 (
     const primitiveMesh& mesh,
@@ -602,8 +596,6 @@ Foam::label Foam::meshTools::otherCell
 }
 
 
-// Returns label of edge nEdges away from startEdge (in the direction of
-// startVertI)
 Foam::label Foam::meshTools::walkFace
 (
     const primitiveMesh& mesh,
@@ -688,7 +680,6 @@ void Foam::meshTools::constrainToMeshCentre
 }
 
 
-//- Set the constrained components of directions/velocity to zero
 void Foam::meshTools::constrainDirection
 (
     const polyMesh& mesh,
@@ -814,7 +805,6 @@ Foam::vector Foam::meshTools::edgeToCutDir
 }
 
 
-// Find edges most aligned with cutDir
 Foam::label Foam::meshTools::cutDirToEdge
 (
     const primitiveMesh& mesh,

src/meshTools/meshTools/meshTools.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -200,7 +200,7 @@ namespace meshTools
             const label v1
         );
 
-        //- Return edge between two vertices. Returns -1 if no edge.
+        //- Return edge between two mesh vertices. Returns -1 if no edge.
         label findEdge
         (
             const primitiveMesh&,

src/meshTools/momentOfInertia/momentOfInertia.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -192,7 +192,8 @@ void Foam::momentOfInertia::massPropertiesShell
     scalar density,
     scalar& mass,
     vector& cM,
-    tensor& J
+    tensor& J,
+    bool doReduce
 )
 {
     // Reset properties for accumulation
@@ -221,6 +222,12 @@ void Foam::momentOfInertia::massPropertiesShell
         mass += triMag;
     }
 
+    if (doReduce)
+    {
+        reduce(cM, sumOp<vector>());
+        reduce(mass, sumOp<scalar>());
+    }
+
     cM /= mass;
 
     mass *= density;
@@ -238,6 +245,11 @@ void Foam::momentOfInertia::massPropertiesShell
             pts[tri[2]]
         ).inertia(cM, density);
     }
+
+    if (doReduce)
+    {
+        reduce(J, sumOp<tensor>());
+    }
 }
 
 
@@ -267,7 +279,8 @@ void Foam::momentOfInertia::massPropertiesShell
     scalar density,
     scalar& mass,
     vector& cM,
-    tensor& J
+    tensor& J,
+    bool doReduce
 )
 {
     triFaceList faces(surf.size());
@@ -277,7 +290,41 @@ void Foam::momentOfInertia::massPropertiesShell
         faces[i] = triFace(surf[i]);
     }
 
-    massPropertiesShell(surf.points(), faces, density, mass, cM, J);
+    massPropertiesShell(surf.points(), faces, density, mass, cM, J, doReduce);
+}
+
+
+void Foam::momentOfInertia::massPropertiesPatch
+(
+    const polyPatch& pp,
+    scalar density,
+    scalar& mass,
+    vector& cM,
+    tensor& J,
+    bool doReduce
+)
+{
+    DynamicList<triFace> faces(3*pp.size());
+
+    // decompose patch faces using triangle fan
+    forAll(pp, faceI)
+    {
+        const face& f = pp[faceI];
+
+        if (f.size() > 2)
+        {
+            const label v0 = 0;
+
+            for (label i = 1; i < f.size() - 1; i++)
+            {
+                faces.append(triFace(f[v0], f[i],f[i + 1]));
+            }
+        }
+    }
+
+    triFaceList triFaces;
+    triFaces.transfer(faces);
+    massPropertiesShell(pp.points(), triFaces, density, mass, cM, J, doReduce);
 }
 
 

src/meshTools/momentOfInertia/momentOfInertia.H

@@ -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) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -72,7 +72,8 @@ public:
         scalar density,
         scalar& mass,
         vector& cM,
-        tensor& J
+        tensor& J,
+        bool doReduce = false
     );
 
     static void massPropertiesSolid
@@ -90,7 +91,18 @@ public:
         scalar density,
         scalar& mass,
         vector& cM,
-        tensor& J
+        tensor& J,
+        bool doReduce = false
+    );
+
+    static void massPropertiesPatch
+    (
+        const polyPatch& pp,
+        scalar density,
+        scalar& mass,
+        vector& cM,
+        tensor& J,
+        bool doReduce = false
     );
 
     static tensor applyParallelAxisTheorem

src/meshTools/regionSplit/localPointRegion.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -242,6 +242,7 @@ void Foam::localPointRegion::countPointRegions
 void Foam::localPointRegion::calcPointRegions
 (
     const polyMesh& mesh,
+    const labelPairList& baffles,
     boolList& candidatePoint
 )
 {
@@ -423,6 +424,13 @@ void Foam::localPointRegion::calcPointRegions
             minEqOpFace(),
             Foam::dummyTransform()  // dummy transformation
         );
+        forAll(baffles, i)
+        {
+            label f0 = baffles[i].first();
+            label f1 = baffles[i].second();
+            minEqOpFace()(minRegion[f0], minRegion[f1]);
+            minRegion[f1] = minRegion[f0];
+        }
     }
 
 
@@ -469,7 +477,7 @@ Foam::localPointRegion::localPointRegion(const polyMesh& mesh)
         }
     }
 
-    calcPointRegions(mesh, candidatePoint);
+    calcPointRegions(mesh, labelPairList(0), candidatePoint);
 }
 
 
@@ -492,7 +500,31 @@ Foam::localPointRegion::localPointRegion
         candidatePoint[candidatePoints[i]] = true;
     }
 
-    calcPointRegions(mesh, candidatePoint);
+    calcPointRegions(mesh, labelPairList(0), candidatePoint);
+}
+
+
+Foam::localPointRegion::localPointRegion
+(
+    const polyMesh& mesh,
+    const labelPairList& baffles,
+    const labelList& candidatePoints
+)
+:
+    //nRegions_(0),
+    meshPointMap_(0),
+    pointRegions_(0),
+    meshFaceMap_(0),
+    faceRegions_(0)
+{
+    boolList candidatePoint(mesh.nPoints(), false);
+
+    forAll(candidatePoints, i)
+    {
+        candidatePoint[candidatePoints[i]] = true;
+    }
+
+    calcPointRegions(mesh, baffles, candidatePoint);
 }
 
 
@@ -630,14 +662,18 @@ Foam::List<Foam::labelPair> Foam::localPointRegion::findDuplicateFacePairs
                     << " processorPolyPatch."
                     << "This is not allowed." << nl
                     << "Face:" << meshFace0
+                    << " fc:" << mesh.faceCentres()[meshFace0]
                     << " is on patch:" << patches[patch0].name()
                     << nl
                     << "Face:" << meshFace1
+                    << " fc:" << mesh.faceCentres()[meshFace1]
                     << " is on patch:" << patches[patch1].name()
                     << abort(FatalError);
             }
-
-            baffles.append(labelPair(meshFace0, meshFace1));
+            else
+            {
+                baffles.append(labelPair(meshFace0, meshFace1));
+            }
         }
     }
     return baffles.shrink();

src/meshTools/regionSplit/localPointRegion.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -101,10 +101,10 @@ class localPointRegion
         void calcPointRegions
         (
             const polyMesh& mesh,
+            const labelPairList& baffles,
             boolList& candidatePoint
         );
 
-
         //- Check if two faces are equal. If forward = false checks f1 in
         //  reverse order.
         static bool isDuplicate
@@ -133,6 +133,14 @@ public:
             const labelList& candidatePoints
         );
 
+        //- Construct from mesh and candidate points for duplication
+        localPointRegion
+        (
+            const polyMesh& mesh,
+            const labelPairList& baffles,
+            const labelList& candidatePoints
+        );
+
 
     // Member Functions
 

src/meshTools/regionSplit/minData.H

@@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::minData
+
+Description
+    For use with FaceCellWave. Transports minimum passive data
+
+SourceFiles
+    minDataI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef minData_H
+#define minData_H
+
+#include "point.H"
+#include "tensor.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyPatch;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class minData Declaration
+\*---------------------------------------------------------------------------*/
+
+class minData
+{
+    // Private data
+
+        //- Starting data
+        label data_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline minData();
+
+        //- Construct from count
+        inline minData(const label data);
+
+
+    // Member Functions
+
+        // Access
+
+            inline label data() const
+            {
+                return data_;
+            }
+
+
+        // Needed by FaceCellWave
+
+            //- Check whether origin has been changed at all or
+            //  still contains original (invalid) value.
+            template<class TrackingData>
+            inline bool valid(TrackingData& td) const;
+
+            //- Check for identical geometrical data. Used for cyclics checking.
+            template<class TrackingData>
+            inline bool sameGeometry
+            (
+                const polyMesh&,
+                const minData&,
+                const scalar,
+                TrackingData& td
+            ) const;
+
+            //- Convert any absolute coordinates into relative to (patch)face
+            //  centre
+            template<class TrackingData>
+            inline void leaveDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFaceI,
+                const point& faceCentre,
+                TrackingData& td
+            );
+
+            //- Reverse of leaveDomain
+            template<class TrackingData>
+            inline void enterDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFaceI,
+                const point& faceCentre,
+                TrackingData& td
+            );
+
+            //- Apply rotation matrix to any coordinates
+            template<class TrackingData>
+            inline void transform
+            (
+                const polyMesh&,
+                const tensor&,
+                TrackingData& td
+            );
+
+            //- Influence of neighbouring face.
+            template<class TrackingData>
+            inline bool updateCell
+            (
+                const polyMesh&,
+                const label thisCellI,
+                const label neighbourFaceI,
+                const minData& neighbourInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of neighbouring cell.
+            template<class TrackingData>
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFaceI,
+                const label neighbourCellI,
+                const minData& neighbourInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of different value on same face.
+            template<class TrackingData>
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFaceI,
+                const minData& neighbourInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Same (like operator==)
+            template<class TrackingData>
+            inline bool equal(const minData&, TrackingData& td) const;
+
+    // Member Operators
+
+        // Needed for List IO
+        inline bool operator==(const minData&) const;
+
+        inline bool operator!=(const minData&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const minData&);
+        friend Istream& operator>>(Istream&, minData&);
+};
+
+
+//- Data associated with minData type are contiguous
+template<>
+inline bool contiguous<minData>()
+{
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "minDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/regionSplit/minDataI.H

@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "polyMesh.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::minData::minData()
+:
+    data_(labelMax)
+{}
+
+
+inline Foam::minData::minData(const label data)
+:
+    data_(data)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class TrackingData>
+inline bool Foam::minData::valid(TrackingData& td) const
+{
+    return data_ != labelMax;
+}
+
+
+template<class TrackingData>
+inline bool Foam::minData::sameGeometry
+(
+    const polyMesh&,
+    const minData&,
+    const scalar,
+    TrackingData&
+) const
+{
+    return true;
+}
+
+
+template<class TrackingData>
+inline void Foam::minData::leaveDomain
+(
+    const polyMesh&,
+    const polyPatch& patch,
+    const label patchFaceI,
+    const point& faceCentre,
+    TrackingData&
+)
+{}
+
+
+template<class TrackingData>
+inline void Foam::minData::transform
+(
+    const polyMesh&,
+    const tensor& rotTensor,
+    TrackingData&
+)
+{}
+
+
+template<class TrackingData>
+inline void Foam::minData::enterDomain
+(
+    const polyMesh&,
+    const polyPatch& patch,
+    const label patchFaceI,
+    const point& faceCentre,
+    TrackingData&
+)
+{}
+
+
+template<class TrackingData>
+inline bool Foam::minData::updateCell
+(
+    const polyMesh&,
+    const label thisCellI,
+    const label neighbourFaceI,
+    const minData& neighbourInfo,
+    const scalar tol,
+    TrackingData&
+)
+{
+    if (neighbourInfo.data_ < data_)
+    {
+        operator=(neighbourInfo);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class TrackingData>
+inline bool Foam::minData::updateFace
+(
+    const polyMesh& mesh,
+    const label thisFaceI,
+    const label neighbourCellI,
+    const minData& neighbourInfo,
+    const scalar tol,
+    TrackingData&
+)
+{
+    // From cell to its faces.
+
+    if (neighbourInfo.data_ < data_)
+    {
+        operator=(neighbourInfo);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class TrackingData>
+inline bool Foam::minData::updateFace
+(
+    const polyMesh&,
+    const label thisFaceI,
+    const minData& neighbourInfo,
+    const scalar tol,
+    TrackingData&
+)
+{
+    // From face to face (e.g. coupled faces)
+    if (neighbourInfo.data_ < data_)
+    {
+        operator=(neighbourInfo);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class TrackingData>
+inline bool Foam::minData::equal
+(
+    const minData& rhs,
+    TrackingData& td
+) const
+{
+    return operator==(rhs);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::minData::operator==
+(
+    const Foam::minData& rhs
+) const
+{
+    return data() == rhs.data();
+}
+
+
+inline bool Foam::minData::operator!=
+(
+    const Foam::minData& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Foam::Ostream& os,
+    const Foam::minData& wDist
+)
+{
+    return os << wDist.data_;
+}
+
+
+Foam::Istream& Foam::operator>>
+(
+    Foam::Istream& is,
+    Foam::minData& wDist
+)
+{
+    return is >> wDist.data_;
+}
+
+
+// ************************************************************************* //

src/meshTools/searchableSurface/searchableCone.H

@@ -0,0 +1,296 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     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::searchableCone
+
+Description
+    Searching on (optionally hollow) cone.
+
+    \heading Function object usage
+    \table
+        Property     | Description             | Required    | Default value
+        point1       | coordinate of endpoint  | yes         |
+        radius1      | radius at point1        | yes         | yes
+        innerRadius1 | inner radius at point1  | no          |
+        point2       | coordinate of endpoint  | yes         |
+        radius2      | radius at point2        | yes         | yes
+        innerRadius2 | inner radius at point2  | no          |
+    \endtable
+
+Note
+    Initial implementation, might suffer from robustness (e.g. radius1==radius2)
+
+SourceFiles
+    searchableCone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableCone_H
+#define searchableCone_H
+
+#include "searchableSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class searchableCone Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableCone
+:
+    public searchableSurface
+{
+    // Private Member Data
+
+        //- 'Left' point
+        const point point1_;
+
+        //- Outer radius at point1
+        const scalar radius1_;
+
+        //- Inner radius at point1
+        const scalar innerRadius1_;
+
+
+        //- 'Right' point
+        const point point2_;
+
+        //- Outer radius at point2
+        const scalar radius2_;
+
+        //- Inner radius at point2
+        const scalar innerRadius2_;
+
+
+        //- Length of vector point2-point1
+        const scalar magDir_;
+
+        //- Normalised vector point2-point1
+        const vector unitDir_;
+
+
+        //- Names of regions
+        mutable wordList regions_;
+
+
+   // Private Member Functions
+
+        //- Find nearest point on cylinder.
+        void findNearestAndNormal
+        (
+            const point& sample,
+            const scalar nearestDistSqr,
+            pointIndexHit & nearInfo,
+            vector& normal
+        ) const;
+
+        //- Determine radial coordinate (squared)
+        static scalar radius2(const searchableCone& cone, const point& pt);
+
+        //- Find both intersections with cone. innerRadii supplied externally.
+        void findLineAll
+        (
+            const searchableCone& cone,
+            const scalar innerRadius1,
+            const scalar innerRadius2,
+            const point& start,
+            const point& end,
+            pointIndexHit& near,
+            pointIndexHit& far
+        ) const;
+
+        //- Insert a hit if it differs (by a tolerance) from the existing ones
+        void insertHit
+        (
+            const point& start,
+            const point& end,
+            List<pointIndexHit>& info,
+            const pointIndexHit& hit
+        ) const;
+
+        //- Return the boundBox of the cylinder
+        boundBox calcBounds() const;
+
+        //- Disallow default bitwise copy construct
+        searchableCone(const searchableCone&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const searchableCone&);
+
+
+public:
+
+    //- Runtime type information
+   TypeName("searchableCone");
+
+
+    // Constructors
+
+        //- Construct from components
+       searchableCone
+       (
+           const IOobject& io,
+           const point& point1,
+           const scalar radius1,
+           const scalar innerRadius1,
+           const point& point2,
+           const scalar radius2,
+           const scalar innerRadius2
+       );
+
+        //- Construct from dictionary (used by searchableSurface)
+        searchableCone
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+
+        virtual ~searchableCone();
+
+
+    // Member Functions
+
+        virtual const wordList& regions() const;
+
+        //- Whether supports volume type below
+        virtual bool hasVolumeType() const
+        {
+            return true;
+        }
+
+        //- Range of local indices that can be returned.
+        virtual label size() const
+        {
+            return 1;
+        }
+
+        //- Get representative set of element coordinates
+        //  Usually the element centres (should be of length size()).
+        virtual tmp<pointField> coordinates() const;
+
+        //- Get bounding spheres (centre and radius squared), one per element.
+        //  Any point on element is guaranteed to be inside.
+        virtual void boundingSpheres
+        (
+            pointField& centres,
+            scalarField& radiusSqr
+        ) const;
+
+        //- Get the points that define the surface.
+        virtual tmp<pointField> points() const;
+
+        //- Does any part of the surface overlap the supplied bound box?
+        virtual bool overlaps(const boundBox& bb) const
+        {
+            notImplemented
+            (
+                "searchableCone::overlaps(const boundBox&) const"
+            );
+
+            return false;
+        }
+
+
+        // Multiple point queries.
+
+            //- Find nearest point on cylinder
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>&
+            ) const;
+
+            //- Find nearest intersection on line from start to end
+            virtual void findLine
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            //- Find all intersections in order from start to end
+            virtual void findLineAll
+            (
+                const pointField& start,
+                const pointField& end,
+                List<List<pointIndexHit> >&
+            ) const;
+
+            //- Find any intersection on line from start to end
+            virtual void findLineAny
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            //- From a set of points and indices get the region
+            virtual void getRegion
+            (
+                const List<pointIndexHit>&,
+                labelList& region
+            ) const;
+
+            //- From a set of points and indices get the normal
+            virtual void getNormal
+            (
+                const List<pointIndexHit>&,
+                vectorField& normal
+            ) const;
+
+            //- Determine type (inside/outside/mixed) for point. unknown if
+            //  cannot be determined (e.g. non-manifold surface)
+            virtual void getVolumeType
+            (
+                const pointField&,
+                List<volumeType>&
+            ) const;
+
+
+        // regIOobject implementation
+
+            virtual bool writeData(Ostream&) const
+            {
+                notImplemented("searchableCone::writeData(Ostream&) const");
+                return false;
+            }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/searchableSurface/searchableCylinder.H

@@ -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) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,8 +43,6 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
-
 /*---------------------------------------------------------------------------*\
                            Class searchableCylinder Declaration
 \*---------------------------------------------------------------------------*/

src/meshTools/searchableSurface/searchableRotatedBox.C

@@ -0,0 +1,400 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "searchableRotatedBox.H"
+#include "addToRunTimeSelectionTable.H"
+#include "axesRotation.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(searchableRotatedBox, 0);
+    addToRunTimeSelectionTable(searchableSurface, searchableRotatedBox, dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchableRotatedBox::searchableRotatedBox
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    searchableSurface(io),
+    box_
+    (
+        IOobject
+        (
+            io.name() + "_box",
+            io.instance(),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            false      //io.registerObject(),
+        ),
+        treeBoundBox(point::zero, dict.lookup("span"))
+    ),
+    transform_
+    (
+        "rotation",
+        dict.lookup("origin"),
+        dict.lookup("e3"),
+        dict.lookup("e1")
+    )
+{
+    points_ = transform_.globalPosition(box_.points());
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchableRotatedBox::~searchableRotatedBox()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::wordList& Foam::searchableRotatedBox::regions() const
+{
+    return box_.regions();
+}
+
+
+Foam::tmp<Foam::pointField> Foam::searchableRotatedBox::coordinates() const
+{
+    return transform_.globalPosition(box_.coordinates());
+}
+
+
+void Foam::searchableRotatedBox::boundingSpheres
+(
+    pointField& centres,
+    scalarField& radiusSqr
+) const
+{
+    box_.boundingSpheres(centres, radiusSqr);
+    centres = transform_.globalPosition(centres);
+}
+
+
+Foam::tmp<Foam::pointField> Foam::searchableRotatedBox::points() const
+{
+    return points_;
+}
+
+
+bool Foam::searchableRotatedBox::overlaps(const boundBox& bb) const
+{
+    // (from treeDataPrimitivePatch.C)
+
+    // 1. bounding box
+    if (!bb.overlaps(bounds()))
+    {
+        return false;
+    }
+
+    // 2. Check if one or more face points inside
+    const faceList& fcs = treeBoundBox::faces;
+    forAll(fcs, faceI)
+    {
+        if (bb.containsAny(points_))
+        {
+            return true;
+        }
+    }
+
+    // 3. Difficult case: all points are outside but connecting edges might
+    // go through cube.
+
+    const treeBoundBox treeBb(bb);
+
+    // 3a. my edges through bb faces
+    const edgeList& edges = treeBoundBox::edges;
+    forAll(edges, edgeI)
+    {
+        const edge& e = edges[edgeI];
+
+        point inter;
+        if (treeBb.intersects(points_[e[0]], points_[e[1]], inter))
+        {
+            return true;
+        }
+    }
+
+    // 3b. bb edges through my faces
+
+    const pointField bbPoints(bb.points());
+
+    forAll(fcs, faceI)
+    {
+        const face& f = fcs[faceI];
+        point fc = f.centre(points_);
+
+        forAll(edges, edgeI)
+        {
+            const edge& e = edges[edgeI];
+
+            pointHit inter = f.intersection
+            (
+                bbPoints[e[0]],
+                bbPoints[e[1]],
+                fc,
+                points_,
+                intersection::HALF_RAY
+            );
+
+            if (inter.hit() && inter.distance() <= 1)
+            {
+                return true;
+            }
+        }
+    }
+
+    return false;
+}
+
+
+Foam::pointIndexHit Foam::searchableRotatedBox::findNearest
+(
+    const point& sample,
+    const scalar nearestDistSqr
+) const
+{
+    pointIndexHit boxNearest
+    (
+        box_.findNearest
+        (
+            transform_.localPosition(sample),
+            nearestDistSqr
+        )
+    );
+
+    boxNearest.rawPoint() = transform_.globalPosition(boxNearest.rawPoint());
+
+    return boxNearest;
+}
+
+
+Foam::pointIndexHit Foam::searchableRotatedBox::findNearest
+(
+    const linePointRef& ln,
+    treeBoundBox& tightest,
+    point& linePoint
+) const
+{
+    notImplemented
+    (
+        "searchableRotatedBox::findNearest"
+        "(const linePointRef&, treeBoundBox&, point&)"
+    );
+    return pointIndexHit();
+}
+
+
+Foam::pointIndexHit Foam::searchableRotatedBox::findLine
+(
+    const point& start,
+    const point& end
+) const
+{
+    pointIndexHit boxHit
+    (
+        box_.findLine
+        (
+            transform_.localPosition(start),
+            transform_.localPosition(end)
+        )
+    );
+
+    boxHit.rawPoint() = transform_.globalPosition(boxHit.rawPoint());
+
+    return boxHit;
+}
+
+
+Foam::pointIndexHit Foam::searchableRotatedBox::findLineAny
+(
+    const point& start,
+    const point& end
+) const
+{
+    return findLine(start, end);
+}
+
+
+void Foam::searchableRotatedBox::findNearest
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(samples.size());
+
+    forAll(samples, i)
+    {
+        info[i] = findNearest(samples[i], nearestDistSqr[i]);
+    }
+}
+
+
+void Foam::searchableRotatedBox::findLine
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(start.size());
+
+    forAll(start, i)
+    {
+        info[i] = findLine(start[i], end[i]);
+    }
+}
+
+
+void Foam::searchableRotatedBox::findLineAny
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(start.size());
+
+    forAll(start, i)
+    {
+        info[i] = findLineAny(start[i], end[i]);
+    }
+}
+
+
+void Foam::searchableRotatedBox::findLineAll
+(
+    const pointField& start,
+    const pointField& end,
+    List<List<pointIndexHit> >& info
+) const
+{
+    info.setSize(start.size());
+
+    // Work array
+    DynamicList<pointIndexHit, 1, 1> hits;
+
+    // Tolerances:
+    // To find all intersections we add a small vector to the last intersection
+    // This is chosen such that
+    // - it is significant (SMALL is smallest representative relative tolerance;
+    //   we need something bigger since we're doing calculations)
+    // - if the start-end vector is zero we still progress
+    const vectorField dirVec(end-start);
+    const scalarField magSqrDirVec(magSqr(dirVec));
+    const vectorField smallVec
+    (
+        Foam::sqrt(SMALL)*dirVec
+      + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
+    );
+
+    forAll(start, pointI)
+    {
+        // See if any intersection between pt and end
+        pointIndexHit inter = findLine(start[pointI], end[pointI]);
+
+        if (inter.hit())
+        {
+            hits.clear();
+            hits.append(inter);
+
+            point pt = inter.hitPoint() + smallVec[pointI];
+
+            while (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
+            {
+                // See if any intersection between pt and end
+                pointIndexHit inter = findLine(pt, end[pointI]);
+
+                // Check for not hit or hit same face as before (can happen
+                // if vector along surface of face)
+                if
+                (
+                    !inter.hit()
+                 || (inter.index() == hits.last().index())
+                )
+                {
+                    break;
+                }
+                hits.append(inter);
+
+                pt = inter.hitPoint() + smallVec[pointI];
+            }
+
+            info[pointI].transfer(hits);
+        }
+        else
+        {
+            info[pointI].clear();
+        }
+    }
+}
+
+
+void Foam::searchableRotatedBox::getRegion
+(
+    const List<pointIndexHit>& info,
+    labelList& region
+) const
+{
+    region.setSize(info.size());
+    region = 0;
+}
+
+
+void Foam::searchableRotatedBox::getNormal
+(
+    const List<pointIndexHit>& info,
+    vectorField& normal
+) const
+{
+    // searchableBox does not use hitPoints so no need to transform
+    box_.getNormal(info, normal);
+
+    normal = transform_.globalVector(normal);
+}
+
+
+void Foam::searchableRotatedBox::getVolumeType
+(
+    const pointField& points,
+    List<volumeType>& volType
+) const
+{
+    box_.getVolumeType(transform_.localPosition(points), volType);
+}
+
+
+// ************************************************************************* //

src/meshTools/searchableSurface/searchableRotatedBox.H

@@ -0,0 +1,265 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::searchableRotatedBox
+
+Description
+    Searching on a rotated box
+
+    Box defined as min and max coordinate. Rotation by coordinate system
+    at box middle.
+
+    E.g. box with sides 1 1 1 rotated 45 degrees around z-axis at
+         origin (0.5 0.5 0.5)
+    \verbatim
+        span     (1 1 1);
+        origin  (0.5 0.5 0.5);
+        e1      (1 1 0);
+        e3      (0 0 1);
+    \endverbatim
+
+SourceFiles
+    searchableRotatedBox.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableRotatedBox_H
+#define searchableRotatedBox_H
+
+#include "searchableSurface.H"
+#include "coordinateSystem.H"
+#include "searchableBox.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class searchableRotatedBox Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableRotatedBox
+:
+    public searchableSurface
+{
+private:
+
+    // Private Member Data
+
+        //- box in local coordinate system
+        searchableBox box_;
+
+        //- transformation from local to global coordinates
+        coordinateSystem transform_;
+
+        //- (global) corner points (in treeBoundBox order)
+        pointField points_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        searchableRotatedBox(const searchableRotatedBox&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const searchableRotatedBox&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("searchableRotatedBox");
+
+
+    // Constructors
+
+        //- Construct from dictionary (used by searchableSurface)
+        searchableRotatedBox
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~searchableRotatedBox();
+
+
+    // Member Functions
+
+        virtual const wordList& regions() const;
+
+        //- Whether supports volume type below
+        virtual bool hasVolumeType() const
+        {
+            return true;
+        }
+
+        //- Range of local indices that can be returned.
+        virtual label size() const
+        {
+            return 6;
+        }
+
+        //- Get representative set of element coordinates
+        //  Usually the element centres (should be of length size()).
+        virtual tmp<pointField> coordinates() const;
+
+        //- Get bounding spheres (centre and radius squared), one per element.
+        //  Any point on element is guaranteed to be inside.
+        virtual void boundingSpheres
+        (
+            pointField& centres,
+            scalarField& radiusSqr
+        ) const;
+
+        //- Get the points that define the surface.
+        virtual tmp<pointField> points() const;
+
+        // Does any part of the surface overlap the supplied bound box?
+        virtual bool overlaps(const boundBox& bb) const;
+
+        // Single point queries.
+
+            //- Inherit findNearest from searchableSurface
+            using searchableSurface::findNearest;
+
+            //- Calculate nearest point on surface.
+            //  Returns
+            //  - bool : any point found nearer than nearestDistSqr
+            //  - label: relevant index in surface (=face 0..5)
+            //  - point: actual nearest point found
+            pointIndexHit findNearest
+            (
+                const point& sample,
+                const scalar nearestDistSqr
+            ) const;
+
+            //- Find nearest to segment.
+            //  Returns
+            //  - bool : any point found?
+            //  - label: relevant index in shapes (=face 0..5)
+            //  - point: actual nearest point found
+            //  sets:
+            //  - tightest  : bounding box
+            //  - linePoint : corresponding nearest point on line
+            pointIndexHit findNearest
+            (
+                const linePointRef& ln,
+                treeBoundBox& tightest,
+                point& linePoint
+            ) const;
+
+            //- Find nearest intersection of line between start and end.
+            pointIndexHit findLine
+            (
+                const point& start,
+                const point& end
+            ) const;
+
+            //- Find any intersection of line between start and end.
+            pointIndexHit findLineAny
+            (
+                const point& start,
+                const point& end
+            ) const;
+
+
+        // Multiple point queries.
+
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLine
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLineAny
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            //- Get all intersections in order from start to end.
+            virtual void findLineAll
+            (
+                const pointField& start,
+                const pointField& end,
+                List<List<pointIndexHit> >&
+            ) const;
+
+            //- From a set of points and indices get the region
+            virtual void getRegion
+            (
+                const List<pointIndexHit>&,
+                labelList& region
+            ) const;
+
+            //- From a set of points and indices get the normal
+            virtual void getNormal
+            (
+                const List<pointIndexHit>&,
+                vectorField& normal
+            ) const;
+
+            //- Determine type (inside/outside/mixed) for point. unknown if
+            //  cannot be determined (e.g. non-manifold surface)
+            virtual void getVolumeType
+            (
+                const pointField&,
+                List<volumeType>&
+            ) const;
+
+
+        // regIOobject implementation
+
+            bool writeData(Ostream&) const
+            {
+                notImplemented
+                (
+                    "searchableRotatedBox::writeData(Ostream&) const"
+                );
+                return false;
+            }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/searchableSurface/searchableSurfaceCollection.H

@@ -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) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,8 +25,8 @@ Class
     Foam::searchableSurfaceCollection
 
 Description
-    Set of transformed searchableSurfaces. Does not do boolean operations.
-    So when meshing might find parts 'inside'.
+    Set of transformed searchableSurfaces. Does not do boolean operations
+    so when meshing might find parts 'inside'.
 
 SourceFiles
     searchableSurfaceCollection.C
@@ -47,8 +47,6 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
-
 /*---------------------------------------------------------------------------*\
                  Class searchableSurfaceCollection Declaration
 \*---------------------------------------------------------------------------*/

src/meshTools/searchableSurface/searchableSurfaceWithGaps.H

@@ -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) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,6 +34,22 @@ Description
     - extend the test vector slightly (with SMALL) to account for numerical
       inaccuracies.
 
+    \verbatim
+    sphere.stl
+    {
+        type triSurfaceMesh;
+    }
+
+    sphere
+    {
+        type searchableSurfaceWithGaps;
+        // Underlying surface
+        surface sphere.stl;
+        // Perturb distance
+        gap 1e-3;
+    }
+    \endverbatim
+
 SourceFiles
     searchableSurfaceWithGaps.C
 

src/meshTools/searchableSurface/searchableSurfaces.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,9 +27,7 @@ License
 #include "searchableSurfacesQueries.H"
 #include "ListOps.H"
 #include "Time.H"
-//#include "vtkSetWriter.H"
 #include "DynamicField.H"
-//#include "OBJstream.H"
 #include "PatchTools.H"
 #include "triSurfaceMesh.H"
 
@@ -405,9 +403,9 @@ void Foam::searchableSurfaces::findNearest
 // Find nearest. Return -1 or nearest point
 void Foam::searchableSurfaces::findNearest
 (
+    const labelListList& regionIndices,
     const pointField& samples,
     const scalarField& nearestDistSqr,
-    const labelList& regionIndices,
     labelList& nearestSurfaces,
     List<pointIndexHit>& nearestInfo
 ) const
@@ -416,9 +414,11 @@ void Foam::searchableSurfaces::findNearest
     (
         *this,
         allSurfaces_,
+        regionIndices,
+
         samples,
         nearestDistSqr,
-        regionIndices,
+
         nearestSurfaces,
         nearestInfo
     );

src/meshTools/searchableSurface/searchableSurfaces.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -200,9 +200,9 @@ public:
 
             void findNearest
             (
+                const labelListList& regionIndices,
                 const pointField& samples,
                 const scalarField& nearestDistSqr,
-                const labelList& regionIndices,
                 labelList& nearestSurfaces,
                 List<pointIndexHit>& nearestInfo
             ) const;

src/meshTools/searchableSurface/searchableSurfacesQueries.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -672,9 +672,11 @@ void Foam::searchableSurfacesQueries::findNearest
 (
     const PtrList<searchableSurface>& allSurfaces,
     const labelList& surfacesToTest,
+    const labelListList& regionIndices,
+
     const pointField& samples,
     const scalarField& nearestDistSqr,
-    const labelList& regionIndices,
+
     labelList& nearestSurfaces,
     List<pointIndexHit>& nearestInfo
 )
@@ -707,7 +709,7 @@ void Foam::searchableSurfacesQueries::findNearest
         (
             samples,
             minDistSqr,
-            regionIndices,
+            regionIndices[testI],
             hitInfo
         );
 

src/meshTools/searchableSurface/searchableSurfacesQueries.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -183,14 +183,14 @@ public:
                 List<pointIndexHit>&
             );
 
-            //- Find nearest points to a specific region of the surface
+            //- Find nearest points to a specific regions of the surface
             static void findNearest
             (
                 const PtrList<searchableSurface>& allSurfaces,
                 const labelList& surfacesToTest,
+                const labelListList& regionIndices,
                 const pointField& samples,
                 const scalarField& nearestDistSqr,
-                const labelList& regionIndices,
                 labelList& nearestSurfaces,
                 List<pointIndexHit>& nearestInfo
             );

src/meshTools/searchableSurface/subTriSurfaceMesh.C

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     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 "subTriSurfaceMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(subTriSurfaceMesh, 0);
+addToRunTimeSelectionTable(searchableSurface, subTriSurfaceMesh, dict);
+
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::wordList Foam::subTriSurfaceMesh::patchNames(const triSurface& s)
+{
+    const geometricSurfacePatchList& patches = s.patches();
+
+    wordList names(patches.size());
+    forAll(patches, patchI)
+    {
+        names[patchI] = patches[patchI].name();
+    }
+    return names;
+}
+
+
+Foam::labelList Foam::subTriSurfaceMesh::selectedRegions
+(
+    const triSurface& s,
+    const wordReList& regionNames
+)
+{
+    const wordList names(patchNames(s));
+
+    labelList regions(names.size());
+
+    label compactI = 0;
+
+    forAll(names, regionI)
+    {
+        const word& name = names[regionI];
+
+        forAll(regionNames, i)
+        {
+            if (regionNames[i].match(name))
+            {
+                regions[compactI++] = regionI;
+            }
+        }
+    }
+    regions.setSize(compactI);
+
+    return regions;
+}
+
+
+Foam::triSurface Foam::subTriSurfaceMesh::subset
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+{
+    const word subGeomName(dict.lookup("surface"));
+
+    const triSurfaceMesh& s =
+        io.db().lookupObject<triSurfaceMesh>(subGeomName);
+
+    const wordReList regionNames(dict.lookup("patches"));
+
+    labelList regionMap(selectedRegions(s, regionNames));
+
+    if (regionMap.size() == 0)
+    {
+        FatalIOErrorIn
+        (
+            "subTriSurfaceMesh::subset"
+            "(\n"
+            "    const IOobject&,\n"
+            "    const dictionary&\n"
+            ")",
+            dict
+        )   << "Found no regions in triSurface matching " << regionNames
+            << ". Valid regions are " << patchNames(s)
+            << exit(FatalIOError);
+    }
+
+    labelList reverseRegionMap(s.patches().size(), -1);
+    forAll(regionMap, i)
+    {
+        reverseRegionMap[regionMap[i]] = i;
+    }
+
+    boolList isSelected(s.size(), false);
+    forAll(s, triI)
+    {
+        if (reverseRegionMap[s.triSurface::operator[](triI).region()] != -1)
+        {
+            isSelected[triI] = true;
+        }
+    }
+
+    labelList pointMap;
+    labelList faceMap;
+    return s.subsetMesh(isSelected, pointMap, faceMap);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::subTriSurfaceMesh::subTriSurfaceMesh
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    triSurfaceMesh(io, subset(io, dict))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::subTriSurfaceMesh::~subTriSurfaceMesh()
+{}
+
+
+// ************************************************************************* //

src/meshTools/searchableSurface/subTriSurfaceMesh.H

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     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::subTriSurfaceMesh
+
+Description
+    On-the-fly subset of triSurfaceMesh (based on selected patches)
+
+Note
+
+    Example usage:
+    \verbatim
+    spheres.stl
+    {
+        type subTriSurfaceMesh;
+
+        // Surface to operate on. Avoid duplicate loading
+        surface spheres.stl;
+
+        // Regions to operate on (regular expressions allowed)
+        patches (solid1);
+    }
+    \endverbatim
+
+
+SourceFiles
+    subTriSurfaceMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef subTriSurfaceMesh_H
+#define subTriSurfaceMesh_H
+
+#include "triSurfaceMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class subTriSurfaceMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class subTriSurfaceMesh
+:
+    public triSurfaceMesh
+{
+    // Private member functions
+
+        //- Extract patch names of triSurface
+        static wordList patchNames(const triSurface& s);
+
+        //- Select regions by name
+        static labelList selectedRegions
+        (
+            const triSurface& s,
+            const wordReList& regionNames
+        );
+
+        //- Subset triSurface based on regions
+        static triSurface subset(const IOobject&, const dictionary&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("subTriSurfaceMesh");
+
+
+    // Constructors
+
+        //- Construct from IO and dictionary.
+        subTriSurfaceMesh
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~subTriSurfaceMesh();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/searchableSurface/triSurfaceMesh.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,6 +39,8 @@ namespace Foam
 defineTypeNameAndDebug(triSurfaceMesh, 0);
 addToRunTimeSelectionTable(searchableSurface, triSurfaceMesh, dict);
 
+word triSurfaceMesh::meshSubDir = "triSurface";
+
 }
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -231,7 +233,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
         (
             io.name(),
             io.instance(),
-            io.local(),
+            io.local(),  //"triSurfaceFields",
             io.db(),
             io.readOpt(),
             io.writeOpt(),
@@ -239,9 +241,10 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
         )
     ),
     triSurface(s),
-    triSurfaceRegionSearch(s),
+    triSurfaceRegionSearch(static_cast<const triSurface&>(*this)),
     minQuality_(-1),
-    surfaceClosed_(-1)
+    surfaceClosed_(-1),
+    outsideVolType_(volumeType::UNKNOWN)
 {
     const pointField& pts = triSurface::points();
 
@@ -272,7 +275,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
         (
             io.name(),
             static_cast<const searchableSurface&>(*this).instance(),
-            io.local(),
+            io.local(), //"triSurfaceFields",
             io.db(),
             io.readOpt(),
             io.writeOpt(),
@@ -289,7 +292,8 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
     ),
     triSurfaceRegionSearch(static_cast<const triSurface&>(*this)),
     minQuality_(-1),
-    surfaceClosed_(-1)
+    surfaceClosed_(-1),
+    outsideVolType_(volumeType::UNKNOWN)
 {
     const pointField& pts = triSurface::points();
 
@@ -323,7 +327,7 @@ Foam::triSurfaceMesh::triSurfaceMesh
         (
             io.name(),
             static_cast<const searchableSurface&>(*this).instance(),
-            io.local(),
+            io.local(), //"triSurfaceFields",
             io.db(),
             io.readOpt(),
             io.writeOpt(),
@@ -340,7 +344,8 @@ Foam::triSurfaceMesh::triSurfaceMesh
     ),
     triSurfaceRegionSearch(static_cast<const triSurface&>(*this), dict),
     minQuality_(-1),
-    surfaceClosed_(-1)
+    surfaceClosed_(-1),
+    outsideVolType_(volumeType::UNKNOWN)
 {
     scalar scaleFactor = 0;
 
@@ -377,6 +382,7 @@ Foam::triSurfaceMesh::~triSurfaceMesh()
 
 void Foam::triSurfaceMesh::clearOut()
 {
+    outsideVolType_ = volumeType::UNKNOWN;
     triSurfaceRegionSearch::clearOut();
     edgeTree_.clear();
     triSurface::clearOut();
@@ -447,9 +453,22 @@ bool Foam::triSurfaceMesh::overlaps(const boundBox& bb) const
 
 void Foam::triSurfaceMesh::movePoints(const pointField& newPoints)
 {
+    outsideVolType_ = volumeType::UNKNOWN;
+
+    // Update local information (instance, event number)
+    searchableSurface::instance() = objectRegistry::time().timeName();
+    objectRegistry::instance() = searchableSurface::instance();
+
+    label event = getEvent();
+    searchableSurface::eventNo() = event;
+    objectRegistry::eventNo() = searchableSurface::eventNo();
+
+    // Clear additional addressing
     triSurfaceRegionSearch::clearOut();
     edgeTree_.clear();
     triSurface::movePoints(newPoints);
+
+    bounds() = boundBox(triSurface::points());
 }
 
 
@@ -714,27 +733,41 @@ void Foam::triSurfaceMesh::getNormal
 
 void Foam::triSurfaceMesh::setField(const labelList& values)
 {
-    autoPtr<triSurfaceLabelField> fldPtr
-    (
-        new triSurfaceLabelField
+    if (foundObject<triSurfaceLabelField>("values"))
+    {
+        triSurfaceLabelField& fld = const_cast<triSurfaceLabelField&>
         (
-            IOobject
+            lookupObject<triSurfaceLabelField>
             (
-                "values",
-                objectRegistry::time().timeName(),  // instance
-                "triSurface",                       // local
+                "values"
+            )
+        );
+        fld.field() = values;
+    }
+    else
+    {
+        autoPtr<triSurfaceLabelField> fldPtr
+        (
+            new triSurfaceLabelField
+            (
+                IOobject
+                (
+                    "values",
+                    objectRegistry::time().timeName(),  // instance
+                    meshSubDir,                         // local
+                    *this,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
                 *this,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            *this,
-            dimless,
-            labelField(values)
-        )
-    );
+                dimless,
+                labelField(values)
+            )
+        );
 
-    // Store field on triMesh
-    fldPtr.ptr()->store();
+        // Store field on triMesh
+        fldPtr.ptr()->store();
+    }
 }
 
 
@@ -781,17 +814,26 @@ void Foam::triSurfaceMesh::getVolumeType
 
         if (!tree().bb().contains(pt))
         {
-            // Have to calculate directly as outside the octree
-            volType[pointI] = tree().shapes().getVolumeType(tree(), pt);
+            if (hasVolumeType())
+            {
+                // Precalculate and cache value for this outside point
+                if (outsideVolType_ == volumeType::UNKNOWN)
+                {
+                    outsideVolType_ = tree().shapes().getVolumeType(tree(), pt);
+                }
+                volType[pointI] = outsideVolType_;
+            }
+            else
+            {
+                // Have to calculate directly as outside the octree
+                volType[pointI] = tree().shapes().getVolumeType(tree(), pt);
+            }
         }
         else
         {
             // - use cached volume type per each tree node
             volType[pointI] = tree().getVolumeType(pt);
         }
-
-//        Info<< "octree : " << pt << " = "
-//            << volumeType::names[volType[pointI]] << endl;
     }
 
     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
@@ -806,6 +848,24 @@ bool Foam::triSurfaceMesh::writeObject
     IOstream::compressionType cmp
 ) const
 {
+    const Time& runTime = searchableSurface::time();
+    const fileName& instance = searchableSurface::instance();
+
+    if
+    (
+        instance != runTime.timeName()
+     && instance != runTime.system()
+     && instance != runTime.caseSystem()
+     && instance != runTime.constant()
+     && instance != runTime.caseConstant()
+    )
+    {
+        const_cast<triSurfaceMesh&>(*this).searchableSurface::instance() =
+            runTime.timeName();
+        const_cast<triSurfaceMesh&>(*this).objectRegistry::instance() =
+            runTime.timeName();
+    }
+
     fileName fullPath(searchableSurface::objectPath());
 
     if (!mkDir(fullPath.path()))

src/meshTools/searchableSurface/triSurfaceMesh.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -85,6 +85,9 @@ private:
         //- Is surface closed
         mutable label surfaceClosed_;
 
+        //- If surface is closed, what is type of outside points
+        mutable volumeType outsideVolType_;
+
 
     // Private Member Functions
 
@@ -137,6 +140,9 @@ public:
     //- Runtime type information
     TypeName("triSurfaceMesh");
 
+    //- Return the mesh sub-directory name (usually "triSurface")
+    static word meshSubDir;
+
 
     // Constructors
 

src/meshTools/sets/cellSources/regionToCell/regionToCell.C

@@ -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) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,7 +47,7 @@ addToRunTimeSelectionTable(topoSetSource, regionToCell, istream);
 Foam::topoSetSource::addToUsageTable Foam::regionToCell::usage_
 (
     regionToCell::typeName,
-    "\n    Usage: regionToCell subCellSet (pt0 .. ptn)\n\n"
+    "\n    Usage: regionToCell subCellSet (pt0 .. ptn) nErode\n\n"
     "    Select all cells in the connected region containing"
     " points (pt0..ptn).\n"
 );
@@ -433,7 +433,7 @@ Foam::regionToCell::regionToCell
     topoSetSource(mesh),
     setName_(checkIs(is)),
     insidePoints_(checkIs(is)),
-    nErode_(0)
+    nErode_(readLabel(checkIs(is)))
 {}
 
 

src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.H

@@ -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) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -65,16 +65,16 @@ class surfaceToPoint
         static addToUsageTable usage_;
 
         //- Name of surface file
-        fileName surfName_;
+        const fileName surfName_;
 
         //- If > 0 : include points with distance to surface less than nearDist.
-        scalar nearDist_;
+        const scalar nearDist_;
 
         //- Include all points on opposite sign of normal on nearest surface.
-        bool includeInside_;
+        const bool includeInside_;
 
         //- Include all points on this sign of normal on nearest surface.
-        bool includeOutside_;
+        const bool includeOutside_;
 
 
     // Private Member Functions

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) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,6 +39,12 @@ namespace Foam
     defineTypeNameAndDebug(tetOverlapVolume, 0);
 }
 
+// When to consider a tet to be zero volume. We want to avoid doing clipping
+// against negative volume tets. Tet volume can be calculated incorrectly
+// due to truncation errors. The value below works for single and double
+// precision but could probably be tighter for double precision.
+Foam::scalar Foam::tetOverlapVolume::minTetVolume_ = SMALL*SMALL;
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -48,79 +54,12 @@ Foam::tetOverlapVolume::tetOverlapVolume()
 
 // * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol
-(
-    const tetPoints& tetA,
-    const tetPoints& tetB
-) const
-{
-    static tetPointRef::tetIntersectionList insideTets;
-    label nInside = 0;
-    static tetPointRef::tetIntersectionList cutInsideTets;
-    label nCutInside = 0;
-
-    tetPointRef::storeOp inside(insideTets, nInside);
-    tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
-    tetPointRef::sumVolOp volInside;
-    tetPointRef::dummyOp outside;
-
-    if ((tetA.tet().mag() < SMALL*SMALL) || (tetB.tet().mag() < SMALL*SMALL))
-    {
-        return 0.0;
-    }
-
-    // face0
-    plane pl0(tetB[1], tetB[3], tetB[2]);
-    tetA.tet().sliceWithPlane(pl0, cutInside, outside);
-    if (nCutInside == 0)
-    {
-        return 0.0;
-    }
-
-    // face1
-    plane pl1(tetB[0], tetB[2], tetB[3]);
-    nInside = 0;
-    for (label i = 0; i < nCutInside; i++)
-    {
-        const tetPointRef t = cutInsideTets[i].tet();
-        t.sliceWithPlane(pl1, inside, outside);
-    }
-    if (nInside == 0)
-    {
-        return 0.0;
-    }
-
-    // face2
-    plane pl2(tetB[0], tetB[3], tetB[1]);
-    nCutInside = 0;
-    for (label i = 0; i < nInside; i++)
-    {
-        const tetPointRef t = insideTets[i].tet();
-        t.sliceWithPlane(pl2, cutInside, outside);
-    }
-    if (nCutInside == 0)
-    {
-        return 0.0;
-    }
-
-    // face3
-    plane pl3(tetB[0], tetB[1], tetB[2]);
-    for (label i = 0; i < nCutInside; i++)
-    {
-        const tetPointRef t = cutInsideTets[i].tet();
-        t.sliceWithPlane(pl3, volInside, outside);
-    }
-
-    return volInside.vol_;
-}
-
-
 Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
 (
     const pointField& points,
     const face& f,
     const point& fc
-) const
+)
 {
     treeBoundBox bb(fc, fc);
     forAll(f, fp)
@@ -145,122 +84,18 @@ bool Foam::tetOverlapVolume::cellCellOverlapMinDecomp
     const scalar threshold
 ) const
 {
-    const cell& cFacesA = meshA.cells()[cellAI];
-    const point& ccA = meshA.cellCentres()[cellAI];
+    hasOverlapOp overlapCheckOp(threshold);
+    cellCellOverlapMinDecomp<hasOverlapOp>
+    (
+        meshA,
+        cellAI,
+        meshB,
+        cellBI,
+        cellBbB,
+        overlapCheckOp
+    );
 
-    const cell& cFacesB = meshB.cells()[cellBI];
-    const point& ccB = meshB.cellCentres()[cellBI];
-
-    scalar vol = 0.0;
-
-    forAll(cFacesA, cFA)
-    {
-        label faceAI = cFacesA[cFA];
-
-        const face& fA = meshA.faces()[faceAI];
-        const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
-        if (!pyrA.overlaps(cellBbB))
-        {
-            continue;
-        }
-
-        bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
-
-        label tetBasePtAI = 0;
-
-        const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
-
-        for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
-        {
-            label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
-            label otherFacePtAI = fA.fcIndex(facePtAI);
-
-            label pt0I = -1;
-            label pt1I = -1;
-
-            if (ownA)
-            {
-                pt0I = fA[facePtAI];
-                pt1I = fA[otherFacePtAI];
-            }
-            else
-            {
-                pt0I = fA[otherFacePtAI];
-                pt1I = fA[facePtAI];
-            }
-
-            const tetPoints tetA
-            (
-                ccA,
-                tetBasePtA,
-                meshA.points()[pt0I],
-                meshA.points()[pt1I]
-            );
-            const treeBoundBox tetABb(tetA.bounds());
-
-
-            // Loop over tets of cellB
-            forAll(cFacesB, cFB)
-            {
-                label faceBI = cFacesB[cFB];
-
-                const face& fB = meshB.faces()[faceBI];
-                const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
-                if (!pyrB.overlaps(pyrA))
-                {
-                    continue;
-                }
-
-                bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
-
-                label tetBasePtBI = 0;
-
-                const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
-
-                for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
-                {
-                    label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
-                    label otherFacePtBI = fB.fcIndex(facePtBI);
-
-                    label pt0I = -1;
-                    label pt1I = -1;
-
-                    if (ownB)
-                    {
-                        pt0I = fB[facePtBI];
-                        pt1I = fB[otherFacePtBI];
-                    }
-                    else
-                    {
-                        pt0I = fB[otherFacePtBI];
-                        pt1I = fB[facePtBI];
-                    }
-
-                    const tetPoints tetB
-                    (
-                        ccB,
-                        tetBasePtB,
-                        meshB.points()[pt0I],
-                        meshB.points()[pt1I]
-                    );
-
-                    if (!tetB.bounds().overlaps(tetABb))
-                    {
-                        continue;
-                    }
-
-                    vol += tetTetOverlapVol(tetA, tetB);
-
-                    if (vol > threshold)
-                    {
-                        return true;
-                    }
-                }
-            }
-        }
-    }
-
-    return false;
+    return overlapCheckOp.ok_;
 }
 
 
@@ -268,121 +103,50 @@ Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
 (
     const primitiveMesh& meshA,
     const label cellAI,
+
     const primitiveMesh& meshB,
     const label cellBI,
     const treeBoundBox& cellBbB
 ) const
 {
-    const cell& cFacesA = meshA.cells()[cellAI];
-    const point& ccA = meshA.cellCentres()[cellAI];
+    sumOverlapOp overlapSumOp;
+    cellCellOverlapMinDecomp<sumOverlapOp>
+    (
+        meshA,
+        cellAI,
+        meshB,
+        cellBI,
+        cellBbB,
+        overlapSumOp
+    );
 
-    const cell& cFacesB = meshB.cells()[cellBI];
-    const point& ccB = meshB.cellCentres()[cellBI];
-
-    scalar vol = 0.0;
-
-    forAll(cFacesA, cFA)
-    {
-        label faceAI = cFacesA[cFA];
-
-        const face& fA = meshA.faces()[faceAI];
-        const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
-        if (!pyrA.overlaps(cellBbB))
-        {
-            continue;
-        }
-
-        bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
-
-        label tetBasePtAI = 0;
-
-        const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
-
-        for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
-        {
-            label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
-            label otherFacePtAI = fA.fcIndex(facePtAI);
-
-            label pt0I = -1;
-            label pt1I = -1;
-
-            if (ownA)
-            {
-                pt0I = fA[facePtAI];
-                pt1I = fA[otherFacePtAI];
-            }
-            else
-            {
-                pt0I = fA[otherFacePtAI];
-                pt1I = fA[facePtAI];
-            }
-
-            const tetPoints tetA
-            (
-                ccA,
-                tetBasePtA,
-                meshA.points()[pt0I],
-                meshA.points()[pt1I]
-            );
-            const treeBoundBox tetABb(tetA.bounds());
+    return overlapSumOp.iop_.vol_;
+}
 
 
-            // Loop over tets of cellB
-            forAll(cFacesB, cFB)
-            {
-                label faceBI = cFacesB[cFB];
+Foam::Tuple2<Foam::scalar, Foam::point>
+Foam::tetOverlapVolume::cellCellOverlapMomentMinDecomp
+(
+    const primitiveMesh& meshA,
+    const label cellAI,
 
-                const face& fB = meshB.faces()[faceBI];
-                const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
-                if (!pyrB.overlaps(pyrA))
-                {
-                    continue;
-                }
+    const primitiveMesh& meshB,
+    const label cellBI,
+    const treeBoundBox& cellBbB
+) const
+{
+    sumOverlapMomentOp overlapSumOp;
+    cellCellOverlapMinDecomp<sumOverlapMomentOp>
+    (
+        meshA,
+        cellAI,
+        meshB,
+        cellBI,
+        cellBbB,
+        overlapSumOp
+    );
 
-                bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
-
-                label tetBasePtBI = 0;
-
-                const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
-
-                for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
-                {
-                    label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
-                    label otherFacePtBI = fB.fcIndex(facePtBI);
-
-                    label pt0I = -1;
-                    label pt1I = -1;
-
-                    if (ownB)
-                    {
-                        pt0I = fB[facePtBI];
-                        pt1I = fB[otherFacePtBI];
-                    }
-                    else
-                    {
-                        pt0I = fB[otherFacePtBI];
-                        pt1I = fB[facePtBI];
-                    }
-
-                    const tetPoints tetB
-                    (
-                        ccB,
-                        tetBasePtB,
-                        meshB.points()[pt0I],
-                        meshB.points()[pt1I]
-                    );
-                    if (!tetB.bounds().overlaps(tetABb))
-                    {
-                        continue;
-                    }
-
-                    vol += tetTetOverlapVol(tetA, tetB);
-                }
-            }
-        }
-    }
-
-    return vol;
+    return overlapSumOp.iop_.vol_;
 }
 
 

src/meshTools/tetOverlapVolume/tetOverlapVolume.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,6 +29,7 @@ Description
 
 SourceFiles
     tetOverlapVolume.C
+    tetOverlapVolumeTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -38,13 +39,14 @@ SourceFiles
 #include "FixedList.H"
 #include "labelList.H"
 #include "treeBoundBox.H"
+#include "Tuple2.H"
+#include "tetrahedron.H"
 
 namespace Foam
 {
 
 class primitiveMesh;
 class polyMesh;
-class tetPoints;
 
 /*---------------------------------------------------------------------------*\
                       Class tetOverlapVolume Declaration
@@ -52,22 +54,130 @@ class tetPoints;
 
 class tetOverlapVolume
 {
+    // Private data
+
+        //- Minimum tet volume to skip test
+        static scalar minTetVolume_;
+
+
+    // Private classes
+
+        //- tetPoints handling : sum resulting volumes
+        class sumMomentOp
+        {
+        public:
+            Tuple2<scalar, point> vol_;
+
+            inline sumMomentOp()
+            :
+                vol_(0.0, vector::zero)
+            {}
+
+            inline void operator()(const tetPoints& tet)
+            {
+                const tetPointRef t(tet.tet());
+                scalar tetVol = t.mag();
+                vol_.first() += tetVol;
+                vol_.second() += (tetVol*t.centre());
+            }
+        };
+
+        //- tetPoints combining : check for overlap
+        class hasOverlapOp
+        {
+        public:
+
+            const scalar threshold_;
+            tetPointRef::sumVolOp iop_;
+            bool ok_;
+
+            inline hasOverlapOp(const scalar threshold)
+            :
+                threshold_(threshold),
+                iop_(),
+                ok_(false)
+            {}
+
+            //- Overlap two tets
+            inline bool operator()(const tetPoints& A, const tetPoints& B)
+            {
+                tetTetOverlap<tetPointRef::sumVolOp>(A, B, iop_);
+                ok_ = (iop_.vol_ > threshold_);
+                return ok_;
+            }
+        };
+
+        //- tetPoints combining : sum overlap volume
+        class sumOverlapOp
+        {
+        public:
+
+            tetPointRef::sumVolOp iop_;
+
+            inline sumOverlapOp()
+            :
+                iop_()
+            {}
+
+            //- Overlap two tets
+            inline bool operator()(const tetPoints& A, const tetPoints& B)
+            {
+                tetTetOverlap<tetPointRef::sumVolOp>(A, B, iop_);
+                return false;
+            }
+        };
+
+        //- tetPoints combining : sum overlap volume
+        class sumOverlapMomentOp
+        {
+        public:
+
+            sumMomentOp iop_;
+
+            inline sumOverlapMomentOp()
+            :
+                iop_()
+            {}
+
+            //- Overlap two tets
+            inline bool operator()(const tetPoints& A, const tetPoints& B)
+            {
+                tetTetOverlap<sumMomentOp>(A, B, iop_);
+                return false;
+            }
+        };
+
+
     // Private member functions
 
-        //- Tet overlap volume
-        scalar tetTetOverlapVol
+        //- Tet overlap calculation
+        template<class tetPointsOp>
+        static void tetTetOverlap
         (
             const tetPoints& tetA,
-            const tetPoints& tetB
-        ) const;
+            const tetPoints& tetB,
+            tetPointsOp& insideOp
+        );
+
+        //- Cell overlap calculation
+        template<class tetsOp>
+        static void cellCellOverlapMinDecomp
+        (
+            const primitiveMesh& meshA,
+            const label cellAI,
+            const primitiveMesh& meshB,
+            const label cellBI,
+            const treeBoundBox& cellBbB,
+            tetsOp& combineTetsOp
+        );
 
         //- Return a const treeBoundBox
-        treeBoundBox pyrBb
+        static treeBoundBox pyrBb
         (
             const pointField& points,
             const face& f,
             const point& fc
-        ) const;
+        );
 
 
 public:
@@ -114,6 +224,17 @@ public:
             const label cellBI,
             const treeBoundBox& cellBbB
         ) const;
+
+        //- Calculates the overlap volume and moment
+        Tuple2<scalar, point> cellCellOverlapMomentMinDecomp
+        (
+            const primitiveMesh& meshA,
+            const label cellAI,
+
+            const primitiveMesh& meshB,
+            const label cellBI,
+            const treeBoundBox& cellBbB
+        ) const;
 };
 
 
@@ -123,6 +244,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "tetOverlapVolumeTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C

@@ -0,0 +1,220 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "tetOverlapVolume.H"
+#include "primitiveMesh.H"
+
+// * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
+
+template<class tetPointsOp>
+void Foam::tetOverlapVolume::tetTetOverlap
+(
+    const tetPoints& tetA,
+    const tetPoints& tetB,
+    tetPointsOp& insideOp
+)
+{
+    static tetPointRef::tetIntersectionList insideTets;
+    label nInside = 0;
+    static tetPointRef::tetIntersectionList cutInsideTets;
+    label nCutInside = 0;
+
+    tetPointRef::storeOp inside(insideTets, nInside);
+    tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
+    tetPointRef::dummyOp outside;
+
+    if (tetA.tet().mag() < minTetVolume_ || tetB.tet().mag() < minTetVolume_)
+    {
+        return;
+    }
+
+    // face0
+    plane pl0(tetB[1], tetB[3], tetB[2]);
+    tetA.tet().sliceWithPlane(pl0, cutInside, outside);
+    if (nCutInside == 0)
+    {
+        return;
+    }
+
+    // face1
+    plane pl1(tetB[0], tetB[2], tetB[3]);
+    nInside = 0;
+    for (label i = 0; i < nCutInside; i++)
+    {
+        const tetPointRef t = cutInsideTets[i].tet();
+        t.sliceWithPlane(pl1, inside, outside);
+    }
+    if (nInside == 0)
+    {
+        return;
+    }
+
+    // face2
+    plane pl2(tetB[0], tetB[3], tetB[1]);
+    nCutInside = 0;
+    for (label i = 0; i < nInside; i++)
+    {
+        const tetPointRef t = insideTets[i].tet();
+        t.sliceWithPlane(pl2, cutInside, outside);
+    }
+    if (nCutInside == 0)
+    {
+        return;
+    }
+
+    // face3
+    plane pl3(tetB[0], tetB[1], tetB[2]);
+    for (label i = 0; i < nCutInside; i++)
+    {
+        const tetPointRef t = cutInsideTets[i].tet();
+        t.sliceWithPlane(pl3, insideOp, outside);
+    }
+}
+
+
+template<class tetsOp>
+void Foam::tetOverlapVolume::cellCellOverlapMinDecomp
+(
+    const primitiveMesh& meshA,
+    const label cellAI,
+
+    const primitiveMesh& meshB,
+    const label cellBI,
+    const treeBoundBox& cellBbB,
+    tetsOp& combineTetsOp
+)
+{
+    const cell& cFacesA = meshA.cells()[cellAI];
+    const point& ccA = meshA.cellCentres()[cellAI];
+
+    const cell& cFacesB = meshB.cells()[cellBI];
+    const point& ccB = meshB.cellCentres()[cellBI];
+
+    forAll(cFacesA, cFA)
+    {
+        label faceAI = cFacesA[cFA];
+
+        const face& fA = meshA.faces()[faceAI];
+        const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
+        if (!pyrA.overlaps(cellBbB))
+        {
+            continue;
+        }
+
+        bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
+
+        label tetBasePtAI = 0;
+
+        const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
+
+        for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
+        {
+            label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
+            label otherFacePtAI = fA.fcIndex(facePtAI);
+
+            label pt0I = -1;
+            label pt1I = -1;
+
+            if (ownA)
+            {
+                pt0I = fA[facePtAI];
+                pt1I = fA[otherFacePtAI];
+            }
+            else
+            {
+                pt0I = fA[otherFacePtAI];
+                pt1I = fA[facePtAI];
+            }
+
+            const tetPoints tetA
+            (
+                ccA,
+                tetBasePtA,
+                meshA.points()[pt0I],
+                meshA.points()[pt1I]
+            );
+            const treeBoundBox tetABb(tetA.bounds());
+
+            // Loop over tets of cellB
+            forAll(cFacesB, cFB)
+            {
+                label faceBI = cFacesB[cFB];
+
+                const face& fB = meshB.faces()[faceBI];
+                const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
+                if (!pyrB.overlaps(pyrA))
+                {
+                    continue;
+                }
+
+                bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
+
+                label tetBasePtBI = 0;
+
+                const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
+
+                for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
+                {
+                    label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
+                    label otherFacePtBI = fB.fcIndex(facePtBI);
+
+                    label pt0I = -1;
+                    label pt1I = -1;
+
+                    if (ownB)
+                    {
+                        pt0I = fB[facePtBI];
+                        pt1I = fB[otherFacePtBI];
+                    }
+                    else
+                    {
+                        pt0I = fB[otherFacePtBI];
+                        pt1I = fB[facePtBI];
+                    }
+
+                    const tetPoints tetB
+                    (
+                        ccB,
+                        tetBasePtB,
+                        meshB.points()[pt0I],
+                        meshB.points()[pt1I]
+                    );
+                    if (!tetB.bounds().overlaps(tetABb))
+                    {
+                        continue;
+                    }
+
+                    if (combineTetsOp(tetA, tetB))
+                    {
+                        return;
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/meshTools/triSurface/booleanOps/booleanSurface/booleanSurface.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,8 +36,24 @@ License
 namespace Foam
 {
 defineTypeNameAndDebug(booleanSurface, 0);
+
+template<>
+const char* Foam::NamedEnum
+<
+    Foam::booleanSurface::booleanOpType,
+    4
+>::names[] =
+{
+    "union",
+    "intersection",
+    "difference",
+    "all"
+};
 }
 
+const Foam::NamedEnum<Foam::booleanSurface::booleanOpType, 4>
+Foam::booleanSurface::booleanOpTypeNames;
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 

src/meshTools/triSurface/booleanOps/booleanSurface/booleanSurface.H

@@ -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) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -170,6 +170,11 @@ public:
                           // surface. (Produces multiply connected surface)
         };
 
+    // Static data
+
+        static const NamedEnum<booleanOpType, 4> booleanOpTypeNames;
+
+
 
     // Constructors
 

src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -716,66 +716,35 @@ Foam::faceList Foam::intersectedSurface::resplitFace
     edgeSurface& eSurf
 )
 {
-    {
-        // Dump face for debugging.
-        Pout<< "Writing face:" << faceI << " to face.obj" << endl;
-        OFstream str("face.obj");
-        writeOBJ(eSurf.points(), eSurf.edges(), eSurf.faceEdges()[faceI], str);
-    }
-
-
     // Count the number of times point has been visited so we
     // can compare number to facePointEdges.
     Map<label> pointVisited(2*facePointEdges.size());
 
-
-    {
-        OFstream str0("visitedNone.obj");
-        OFstream str1("visitedOnce.obj");
-        OFstream str2("visitedTwice.obj");
-        forAll(eSurf.points(), pointI)
-        {
-            const point& pt = eSurf.points()[pointI];
-
-            str0 << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
-            str1 << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
-            str2 << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
-        }
-
-
     forAllConstIter(Map<label>, visited, iter)
     {
         label edgeI = iter.key();
-
         const edge& e = eSurf.edges()[edgeI];
-
         label stat = iter();
 
         if (stat == STARTTOEND || stat == ENDTOSTART)
         {
             incCount(pointVisited, e[0], 1);
             incCount(pointVisited, e[1], 1);
-
-            str1 << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
         }
         else if (stat == BOTH)
         {
             incCount(pointVisited, e[0], 2);
             incCount(pointVisited, e[1], 2);
-
-            str2 << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
         }
         else if (stat == UNVISITED)
         {
             incCount(pointVisited, e[0], 0);
             incCount(pointVisited, e[1], 0);
-
-            str0 << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
         }
     }
-    }
 
 
+    if (debug)
     {
         forAllConstIter(Map<label>, pointVisited, iter)
         {
@@ -838,8 +807,6 @@ Foam::faceList Foam::intersectedSurface::resplitFace
 
 
     // Find second intersection.
-    label visitedVert1 = -1;
-    label unvisitedVert1 = -1;
     {
         scalar minDist = GREAT;
 
@@ -876,8 +843,6 @@ Foam::faceList Foam::intersectedSurface::resplitFace
                     if (nearDist < minDist)
                     {
                         minDist = nearDist;
-                        visitedVert1 = nearVertI;
-                        unvisitedVert1 = pointI;
                     }
                 }
             }
@@ -885,32 +850,24 @@ Foam::faceList Foam::intersectedSurface::resplitFace
     }
 
 
-    Pout<< "resplitFace : adding intersection from " << visitedVert0
-        << " to " << unvisitedVert0 << endl
-        << " and from " << visitedVert1
-        << " to " << unvisitedVert1 << endl;
-
-
-//    // Add the new intersection edges to the edgeSurface
-//    edgeList additionalEdges(2);
-//    additionalEdges[0] = edge(visitedVert0, unvisitedVert0);
-//    additionalEdges[1] = edge(visitedVert1, unvisitedVert1);
-
     // Add the new intersection edges to the edgeSurface
     edgeList additionalEdges(1);
     additionalEdges[0] = edge(visitedVert0, unvisitedVert0);
 
     eSurf.addIntersectionEdges(faceI, additionalEdges);
 
-    fileName newFName("face_" + Foam::name(faceI) + "_newEdges.obj");
-    Pout<< "Dumping face:" << faceI << " to " << newFName << endl;
-    writeLocalOBJ
-    (
-        eSurf.points(),
-        eSurf.edges(),
-        eSurf.faceEdges()[faceI],
-        newFName
-    );
+    if (debug)
+    {
+        fileName newFName("face_" + Foam::name(faceI) + "_newEdges.obj");
+        Pout<< "Dumping face:" << faceI << " to " << newFName << endl;
+        writeLocalOBJ
+        (
+            eSurf.points(),
+            eSurf.edges(),
+            eSurf.faceEdges()[faceI],
+            newFName
+        );
+    }
 
     // Retry splitFace. Use recursion since is rare situation.
     return splitFace(surf, faceI, eSurf);
@@ -1085,20 +1042,6 @@ Foam::faceList Foam::intersectedSurface::splitFace
         }
         else if (stat != BOTH)
         {
-            //{
-            //    Pout<< "Dumping faces so far to faces.obj" << nl
-            //        << faces << endl;
-            //
-            //    OFstream str("faces.obj");
-            //
-            //    forAll(faces, i)
-            //    {
-            //        writeOBJ(points, faces[i], str);
-            //    }
-            //}
-
-            Pout<< "** Resplitting **" << endl;
-
             // Redo face after introducing extra edge. Edge introduced
             // should be one nearest to any fully visited edge.
             return resplitFace
@@ -1220,34 +1163,6 @@ Foam::intersectedSurface::intersectedSurface
             forAll(newFaces, newFaceI)
             {
                 const face& newF = newFaces[newFaceI];
-
-//                {
-//                    fileName fName
-//                    (
-//                        "face_"
-//                      + Foam::name(faceI)
-//                      + "_subFace_"
-//                      + Foam::name(newFaceI)
-//                      + ".obj"
-//                    );
-//                    Pout<< "Writing original face:" << faceI << " subFace:"
-//                        << newFaceI << " to " << fName << endl;
-//
-//                    OFstream str(fName);
-//
-//                    forAll(newF, fp)
-//                    {
-//                        meshTools::writeOBJ(str, eSurf.points()[newF[fp]]);
-//                    }
-//                    str << 'l';
-//                    forAll(newF, fp)
-//                    {
-//                        str << ' ' << fp+1;
-//                    }
-//                    str<< " 1" << nl;
-//                }
-
-
                 const vector& n = surf.faceNormals()[faceI];
                 const label region = surf[faceI].region();
 

src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,7 +54,6 @@ void Foam::edgeIntersections::checkEdges(const triSurface& surf)
 {
     const pointField& localPoints = surf.localPoints();
     const edgeList& edges = surf.edges();
-    const labelListList& edgeFaces = surf.edgeFaces();
 
     treeBoundBox bb(localPoints);
 
@@ -78,19 +77,6 @@ void Foam::edgeIntersections::checkEdges(const triSurface& surf)
                 << "This might lead to problems in intersection"
                 << endl;
         }
-
-        if (edgeFaces[edgeI].size() == 1)
-        {
-            WarningIn
-            (
-                "Foam::edgeIntersections::checkEdges(const triSurface& surf)"
-            )   << "Edge " << edgeI << " vertices " << e
-                << " coords:" << localPoints[e[0]] << ' '
-                << localPoints[e[1]] << " has only one face connected to it:"
-                << edgeFaces[edgeI] << endl
-                << "This might lead to problems in intersection"
-                << endl;
-        }
     }
 }
 
@@ -466,7 +452,6 @@ Foam::edgeIntersections::edgeIntersections()
 {}
 
 
-// Construct from surface and tolerance
 Foam::edgeIntersections::edgeIntersections
 (
     const triSurface& surf1,
@@ -478,16 +463,9 @@ Foam::edgeIntersections::edgeIntersections
     classification_(surf1.nEdges())
 {
     checkEdges(surf1);
-    checkEdges(query2.surface());
 
     // Current set of edges to test
-    labelList edgesToTest(surf1.nEdges());
-
-    // Start off with all edges
-    forAll(edgesToTest, i)
-    {
-        edgesToTest[i] = i;
-    }
+    labelList edgesToTest(identity(surf1.nEdges()));
 
     // Determine intersections for edgesToTest
     intersectEdges

src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -242,9 +242,9 @@ void Foam::surfaceIntersection::storeIntersection
                 WarningIn
                 (
                     "Foam::surfaceIntersection::storeIntersection"
-                    "(const bool isFirstSurf, const labelList& facesA,"
-                    "const label faceB, DynamicList<edge>& allCutEdges,"
-                    "DynamicList<point>& allCutPoints)"
+                    "(const bool, const labelList&,"
+                    "const label, DynamicList<edge>&,"
+                    "DynamicList<point>&)"
                 )   << "Encountered degenerate edge between face "
                     << twoFaces[0] << " on first surface"
                     << " and face " << twoFaces[1] << " on second surface"
@@ -815,7 +815,7 @@ Foam::surfaceIntersection::surfaceIntersection
 }
 
 
-// Construct from full intersection Poutrmation
+// Construct from full intersection information
 Foam::surfaceIntersection::surfaceIntersection
 (
     const triSurface& surf1,

src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -163,7 +163,7 @@ Foam::triSurfaceRegionSearch::treeByRegion() const
                 (
                     treeDataIndirectTriSurface
                     (
-                        true,
+                        false,              //true,
                         indirectRegionPatches_[regionI],
                         tolerance()
                     ),
@@ -211,6 +211,7 @@ void Foam::triSurfaceRegionSearch::findNearest
             }
 
             const treeType& octree = octrees[treeI];
+            const treeDataIndirectTriSurface::findNearestOp nearOp(octree);
 
             forAll(samples, i)
             {
@@ -223,7 +224,7 @@ void Foam::triSurfaceRegionSearch::findNearest
                 (
                     samples[i],
                     nearestDistSqr[i],
-                    treeDataIndirectTriSurface::findNearestOp(octree)
+                    nearOp
                 );
 
                 if

src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -234,7 +234,7 @@ Foam::triSurfaceSearch::tree() const
         (
             new indexedOctree<treeDataTriSurface>
             (
-                treeDataTriSurface(true, surface_, tolerance_),
+                treeDataTriSurface(false, surface_, tolerance_),
                 bb,
                 maxTreeDepth_,  // maxLevel
                 10,             // leafsize
@@ -290,15 +290,17 @@ void Foam::triSurfaceSearch::findNearest
 
     const indexedOctree<treeDataTriSurface>& octree = tree();
 
+    const treeDataTriSurface::findNearestOp fOp(octree);
+
     info.setSize(samples.size());
 
     forAll(samples, i)
     {
-        static_cast<pointIndexHit&>(info[i]) = octree.findNearest
+        info[i] = octree.findNearest
         (
             samples[i],
             nearestDistSqr[i],
-            treeDataTriSurface::findNearestOp(octree)
+            fOp
         );
     }
 
@@ -335,11 +337,7 @@ void Foam::triSurfaceSearch::findLine
 
     forAll(start, i)
     {
-        static_cast<pointIndexHit&>(info[i]) = octree.findLine
-        (
-            start[i],
-            end[i]
-        );
+        info[i] = octree.findLine(start[i], end[i]);
     }
 
     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
@@ -362,11 +360,7 @@ void Foam::triSurfaceSearch::findLineAny
 
     forAll(start, i)
     {
-        static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
-        (
-            start[i],
-            end[i]
-        );
+        info[i] = octree.findLineAny(start[i], end[i]);
     }
 
     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;