diff --git a/applications/utilities/preProcessing/mapFields/mapFields.C b/applications/utilities/preProcessing/mapFields/mapFields.C
index 0fb8e40fff..5d4a6471e4 100644
--- a/applications/utilities/preProcessing/mapFields/mapFields.C
+++ b/applications/utilities/preProcessing/mapFields/mapFields.C
@@ -504,7 +504,7 @@ int main(int argc, char *argv[])
                     !bbsTargetSet[procITarget]
                   || (
                       bbsTargetSet[procITarget]
-                   && bbsTarget[procITarget].intersects(bbSource)
+                   && bbsTarget[procITarget].overlaps(bbSource)
                      )
                 )
                 {
@@ -533,7 +533,7 @@ int main(int argc, char *argv[])
                     bbsTarget[procITarget] = meshTarget.bounds();
                     bbsTargetSet[procITarget] = true;
 
-                    if (bbsTarget[procITarget].intersects(bbSource))
+                    if (bbsTarget[procITarget].overlaps(bbSource))
                     {
                         if (consistent)
                         {
diff --git a/src/OpenFOAM/containers/Lists/SortableList/SortableList.C b/src/OpenFOAM/containers/Lists/SortableList/SortableList.C
index 7cde863a63..a4264578d1 100644
--- a/src/OpenFOAM/containers/Lists/SortableList/SortableList.C
+++ b/src/OpenFOAM/containers/Lists/SortableList/SortableList.C
@@ -86,28 +86,7 @@ void Foam::SortableList<Type>::sort()
         indices_[i] = i;
     }
 
-    Foam::sort(indices_, less(*this));
-
-    List<Type> tmpValues(this->size());
-
-    forAll(indices_, i)
-    {
-        tmpValues[i] = this->operator[](indices_[i]);
-    }
-
-    List<Type>::transfer(tmpValues);
-}
-
-
-
-template <class Type>
-void Foam::SortableList<Type>::stableSort()
-{
-    forAll(indices_, i)
-    {
-        indices_[i] = i;
-    }
-
+    //Foam::sort(indices_, less(*this));
     Foam::stableSort(indices_, less(*this));
 
     List<Type> tmpValues(this->size());
diff --git a/src/OpenFOAM/containers/Lists/SortableList/SortableList.H b/src/OpenFOAM/containers/Lists/SortableList/SortableList.H
index a9ab2ae829..dd0f0af0d4 100644
--- a/src/OpenFOAM/containers/Lists/SortableList/SortableList.H
+++ b/src/OpenFOAM/containers/Lists/SortableList/SortableList.H
@@ -109,12 +109,9 @@ public:
         //- Size the list. If grow can cause undefined indices (until next sort)
         void setSize(const label);
 
-        //- Sort the list (if changed after construction time)
+        //- (stable) sort the list (if changed after construction time)
         void sort();
 
-        //- Sort the list (if changed after construction time)
-        void stableSort();
-
 
     // Member Operators
 
diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.H b/src/OpenFOAM/meshes/boundBox/boundBox.H
index a11aa68266..40287e53eb 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBox.H
+++ b/src/OpenFOAM/meshes/boundBox/boundBox.H
@@ -118,7 +118,7 @@ public:
         // Query
 
             //- Intersects other boundingbox?
-            bool intersects(const boundBox& bb) const
+            bool overlaps(const boundBox& bb) const
             {
                 if
                 (
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C
index 796cb0983a..5421b328d8 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C
@@ -31,7 +31,11 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::mapDistribute::calcSchedule() const
+Foam::List<Foam::labelPair> Foam::mapDistribute::schedule
+(
+    const labelListList& subMap,
+    const labelListList& constructMap
+)
 {
     // Communications: send and receive processor
     List<labelPair> allComms;
@@ -40,16 +44,16 @@ void Foam::mapDistribute::calcSchedule() const
         HashSet<labelPair, labelPair::Hash<> > commsSet(Pstream::nProcs());
 
         // Find what communication is required
-        forAll(subMap_, procI)
+        forAll(subMap, procI)
         {
             if (procI != Pstream::myProcNo())
             {
-                if (subMap_[procI].size() > 0)
+                if (subMap[procI].size() > 0)
                 {
                     // I need to send to procI
                     commsSet.insert(labelPair(Pstream::myProcNo(), procI));
                 }
-                if (constructMap_[procI].size() > 0)
+                if (constructMap[procI].size() > 0)
                 {
                     // I need to receive from procI
                     commsSet.insert(labelPair(procI, Pstream::myProcNo()));
@@ -120,13 +124,7 @@ void Foam::mapDistribute::calcSchedule() const
     );
 
     // Processors involved in my schedule
-    schedulePtr_.reset
-    (
-        new List<labelPair>
-        (
-            IndirectList<labelPair>(allComms, mySchedule)
-        )
-    );
+    return IndirectList<labelPair>(allComms, mySchedule);
 
 
     //if (debug)
@@ -152,6 +150,22 @@ void Foam::mapDistribute::calcSchedule() const
 }
 
 
+const Foam::List<Foam::labelPair>& Foam::mapDistribute::schedule() const
+{
+    if (!schedulePtr_.valid())
+    {
+        schedulePtr_.reset
+        (
+            new List<labelPair>
+            (
+                schedule(subMap_, constructMap_)
+            )
+        );
+    }
+    return schedulePtr_();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 //- Construct from components
@@ -257,13 +271,4 @@ Foam::mapDistribute::mapDistribute
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
index bdc71c286d..65888ea79e 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
@@ -36,6 +36,8 @@ Note:
     Schedule is a list of processor pairs (one send, one receive. One of
     them will be myself) which forms a scheduled (i.e. non-buffered) exchange.
     See distribute on how to use it.
+    Note2: number of items send on one processor have to equal the number
+    of items received on the other processor.
 
 
 SourceFiles
@@ -80,8 +82,6 @@ class mapDistribute
 
     // Private Member Functions
 
-        void calcSchedule() const;
-
         //- Disallow default bitwise copy construct
         mapDistribute(const mapDistribute&);
 
@@ -142,15 +142,15 @@ public:
                 return constructMap_;
             }
 
+            //- Calculate a schedule. See above.
+            static List<labelPair> schedule
+            (
+                const labelListList& subMap,
+                const labelListList& constructMap
+            );
+
             //- Return a schedule. Demand driven. See above.
-            const List<labelPair>& schedule() const
-            {
-                if (!schedulePtr_.valid())
-                {
-                    calcSchedule();
-                }
-                return schedulePtr_();
-            }
+            const List<labelPair>& schedule() const;
 
 
         // Other
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
index 86ee3d96cf..6ad85c7971 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
@@ -48,15 +48,28 @@ void Foam::mapDistribute::distribute
         // Send sub field to neighbour
         for (label domain = 0; domain < Pstream::nProcs(); domain++)
         {
-            if (domain != Pstream::myProcNo())
+            const labelList& map = subMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
             {
+                List<T> subField(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = field[map[i]];
+                }
                 OPstream toNbr(Pstream::blocking, domain);
-                toNbr << IndirectList<T>(field, subMap[domain])();
+                toNbr << subField;
             }
         }
 
         // Subset myself
-        List<T> subField(IndirectList<T>(field, subMap[Pstream::myProcNo()]));
+        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+        List<T> subField(mySubMap.size());
+        forAll(mySubMap, i)
+        {
+            subField[i] = field[mySubMap[i]];
+        }
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -71,7 +84,11 @@ void Foam::mapDistribute::distribute
         // Receive sub field from neighbour
         for (label domain = 0; domain < Pstream::nProcs(); domain++)
         {
-            if (domain != Pstream::myProcNo())
+            if
+            (
+                domain != Pstream::myProcNo()
+             && constructMap[domain].size() > 0
+            )
             {
                 IPstream fromNbr(Pstream::blocking, domain);
                 List<T> subField(fromNbr);
@@ -93,7 +110,13 @@ void Foam::mapDistribute::distribute
         List<T> newField(constructSize);
 
         // Subset myself
-        List<T> subField(IndirectList<T>(field, subMap[Pstream::myProcNo()]));
+        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+        List<T> subField(mySubMap.size());
+        forAll(mySubMap, i)
+        {
+            subField[i] = field[mySubMap[i]];
+        }
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -112,8 +135,16 @@ void Foam::mapDistribute::distribute
             if (Pstream::myProcNo() == sendProc)
             {
                 // I am sender. Send to recvProc.
+                const labelList& map = subMap[recvProc];
+
+                List<T> subField(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = field[map[i]];
+                }
+
                 OPstream toNbr(Pstream::scheduled, recvProc);
-                toNbr << IndirectList<T>(field, subMap[recvProc])();
+                toNbr << subField;
             }
             else
             {
@@ -136,7 +167,13 @@ void Foam::mapDistribute::distribute
         List<T> newField(constructSize);
 
         // Subset myself
-        List<T> subField(IndirectList<T>(field, subMap[Pstream::myProcNo()]));
+        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+        List<T> subField(mySubMap.size());
+        forAll(mySubMap, i)
+        {
+            subField[i] = field[mySubMap[i]];
+        }
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -149,10 +186,19 @@ void Foam::mapDistribute::distribute
         // Send sub field to neighbour
         for (label domain = 0; domain < Pstream::nProcs(); domain++)
         {
-            if (domain != Pstream::myProcNo())
+            const labelList& map = subMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
             {
+
+                List<T> subField(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = field[map[i]];
+                }
+
                 OPstream toNbr(Pstream::nonBlocking, domain);
-                toNbr << IndirectList<T>(field, subMap[domain])();
+                toNbr << subField;
             }
         }
 
@@ -160,13 +206,13 @@ void Foam::mapDistribute::distribute
         // Receive sub field from neighbour
         for (label domain = 0; domain < Pstream::nProcs(); domain++)
         {
-            if (domain != Pstream::myProcNo())
+            const labelList& map = constructMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
             {
                 IPstream fromNbr(Pstream::nonBlocking, domain);
                 List<T> subField(fromNbr);
 
-                const labelList& map = constructMap[domain];
-
                 forAll(map, i)
                 {
                     newField[map[i]] = subField[i];
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
index af86f37fed..f4e073d4e1 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "coupledPolyPatch.H"
-#include "SortableList.H"
 #include "ListOps.H"
 #include "transform.H"
 #include "OFstream.H"
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C
index a47c9c7d8c..835f884399 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C
@@ -68,6 +68,42 @@ PrimitivePatch<Face, FaceList, PointField, PointType>::PrimitivePatch
 {}
 
 
+// Construct from components
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+PrimitivePatch<Face, FaceList, PointField, PointType>::PrimitivePatch
+(
+    FaceList<Face>& faces,
+    Field<PointType>& points,
+    const bool reUse
+)
+:
+    FaceList<Face>(faces, reUse),
+    points_(points, reUse),
+    edgesPtr_(NULL),
+    nInternalEdges_(-1),
+    boundaryPointsPtr_(NULL),
+    faceFacesPtr_(NULL),
+    edgeFacesPtr_(NULL),
+    faceEdgesPtr_(NULL),
+    pointEdgesPtr_(NULL),
+    pointFacesPtr_(NULL),
+    localFacesPtr_(NULL),
+    meshPointsPtr_(NULL),
+    meshPointMapPtr_(NULL),
+    edgeLoopsPtr_(NULL),
+    localPointsPtr_(NULL),
+    localPointOrderPtr_(NULL),
+    faceNormalsPtr_(NULL),
+    pointNormalsPtr_(NULL)
+{}
+
+
 // Construct as copy
 template
 <
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
index c22514d25f..49b3be8260 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
@@ -235,6 +235,14 @@ public:
             const Field<PointType>& points
         );
 
+        //- Construct from components, reuse storage
+        PrimitivePatch
+        (
+            FaceList<Face>& faces,
+            Field<PointType>& points,
+            const bool reUse
+        );
+
         //- Construct as copy
         PrimitivePatch
         (
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C
index 8b5fa7bb03..6b3039ff30 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C
@@ -467,7 +467,8 @@ const edgeList& primitiveMesh::edges() const
 {
     if (!edgesPtr_)
     {
-        calcEdges(true);
+        //calcEdges(true);
+        calcEdges(false);
     }
 
     return *edgesPtr_;
@@ -477,10 +478,8 @@ const labelListList& primitiveMesh::pointEdges() const
 {
     if (!pePtr_)
     {
-        //// Invert edges
-        //pePtr_ = new labelListList(nPoints());
-        //invertManyToMany(nPoints(), edges(), *pePtr_);
-        calcEdges(true);
+        //calcEdges(true);
+        calcEdges(false);
     }
 
     return *pePtr_;
@@ -491,12 +490,53 @@ const labelListList& primitiveMesh::faceEdges() const
 {
     if (!fePtr_)
     {
-        calcEdges(true);
+        if (debug)
+        {
+            Pout<< "primitiveMesh::faceEdges() : "
+                << "calculating faceEdges" << endl;
+        }
+
+        //calcEdges(true);
+        const faceList& fcs = faces();
+        const labelListList& pe = pointEdges();
+        const edgeList& es = edges();
+
+        fePtr_ = new labelListList(fcs.size());
+        labelListList& faceEdges = *fePtr_;
+
+        forAll(fcs, faceI)
+        {
+            const face& f = fcs[faceI];
+
+            labelList& fEdges = faceEdges[faceI];
+            fEdges.setSize(f.size());
+
+            forAll(f, fp)
+            {
+                label pointI = f[fp];
+                label nextPointI = f[f.fcIndex(fp)];
+
+                // Find edge between pointI, nextPontI
+                const labelList& pEdges = pe[pointI];
+
+                forAll(pEdges, i)
+                {
+                    label edgeI = pEdges[i];
+
+                    if (es[edgeI].otherVertex(pointI) == nextPointI)
+                    {
+                        fEdges[fp] = edgeI;
+                        break;
+                    }
+                }
+            }
+        }
     }
 
     return *fePtr_;
 }
 
+
 void primitiveMesh::clearOutEdges()
 {
     deleteDemandDrivenData(edgesPtr_);
diff --git a/src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C b/src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C
index a5778bcdf0..8b3092f08f 100644
--- a/src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C
+++ b/src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C
@@ -41,8 +41,7 @@ extern "C"
 #   include "parmetis.h"
 }
 
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
@@ -57,6 +56,8 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
 //- Does prevention of 0 cell domains and calls parmetis.
 Foam::label Foam::parMetisDecomp::decompose
 (
@@ -76,6 +77,16 @@ Foam::label Foam::parMetisDecomp::decompose
     // Number of dimensions
     int nDims = 3;
 
+
+    if (cellCentres.size() != xadj.size()-1)
+    {
+        FatalErrorIn("parMetisDecomp::decompose(..)")
+            << "cellCentres:" << cellCentres.size()
+            << " xadj:" << xadj.size()
+            << abort(FatalError);
+    }
+
+
     // Get number of cells on all processors
     List<int> nLocalCells(Pstream::nProcs());
     nLocalCells[Pstream::myProcNo()] = xadj.size()-1;
@@ -106,12 +117,12 @@ Foam::label Foam::parMetisDecomp::decompose
     // Make sure every domain has at least one cell
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // (Metis falls over with zero sized domains)
-    // Trickle cells from processors that have them down to those that
+    // Trickle cells from processors that have them up to those that
     // don't.
 
 
-    // Number of cells to send down (is same as number of cells next processor
-    // has to receive)
+    // Number of cells to send to the next processor
+    // (is same as number of cells next processor has to receive)
     List<int> nSendCells(Pstream::nProcs(), 0);
 
     for (label procI = nLocalCells.size()-1; procI >=1; procI--)
@@ -135,6 +146,15 @@ Foam::label Foam::parMetisDecomp::decompose
         Field<int> prevCellWeights(fromPrevProc);
         Field<int> prevFaceWeights(fromPrevProc);
 
+        if (prevXadj.size() != nSendCells[Pstream::myProcNo()-1])
+        {
+            FatalErrorIn("parMetisDecomp::decompose(..)")
+                << "Expected from processor " << Pstream::myProcNo()-1
+                << " connectivity for " << nSendCells[Pstream::myProcNo()-1]
+                << " nCells but only received " << prevXadj.size()
+                << abort(FatalError);
+        }
+
         // Insert adjncy
         prepend(prevAdjncy, adjncy);
         // Adapt offsets and prepend xadj
@@ -222,6 +242,14 @@ Foam::label Foam::parMetisDecomp::decompose
     }
 
 
+    if (nLocalCells[Pstream::myProcNo()] != (xadj.size()-1))
+    {
+        FatalErrorIn("parMetisDecomp::decompose(..)")
+            << "Have connectivity for " << xadj.size()-1
+            << " cells but nLocalCells:" << nLocalCells[Pstream::myProcNo()]
+            << abort(FatalError);
+    }
+
     // Weight info
     int wgtFlag = 0;
     int* vwgtPtr = NULL;
@@ -292,6 +320,15 @@ Foam::label Foam::parMetisDecomp::decompose
 
         List<int> nextFinalDecomp(fromNextProc);
 
+        if (nextFinalDecomp.size() != nSendCells[Pstream::myProcNo()])
+        {
+            FatalErrorIn("parMetisDecomp::decompose(..)")
+                << "Expected from processor " << Pstream::myProcNo()+1
+                << " decomposition for " << nSendCells[Pstream::myProcNo()]
+                << " nCells but only received " << nextFinalDecomp.size()
+                << abort(FatalError);
+        }
+
         append(nextFinalDecomp, finalDecomp);
     }
 
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
index fedd050a78..da7e6368c1 100644
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
@@ -24,7 +24,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "SortableList.H"
 #include "dynamicRefineFvMesh.H"
 #include "addToRunTimeSelectionTable.H"
 #include "volFields.H"
@@ -32,7 +31,6 @@ License
 #include "surfaceFields.H"
 #include "fvCFD.H"
 #include "syncTools.H"
-#include "ListListOps.H"
 #include "pointFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C b/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C
index 2249bc006e..2ce870af9d 100644
--- a/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C
+++ b/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C
@@ -399,7 +399,7 @@ bool Foam::octreeDataFaceList::overlaps
     const treeBoundBox& sampleBb
 ) const
 {
-    return sampleBb.intersects(allBb_[index]);
+    return sampleBb.overlaps(allBb_[index]);
 }
 
 
diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
index 9c6a50177a..6d97285d73 100644
--- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
+++ b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
@@ -32,7 +32,6 @@ License
 #include "octreeDataFace.H"
 #include "octree.H"
 #include "OFstream.H"
-#include "SortableList.H"
 #include "IndirectList.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
index 83f621363e..e519da59d9 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
@@ -29,6 +29,8 @@ License
 #include "ListListOps.H"
 #include "meshSearch.H"
 #include "mapDistribute.H"
+#include "meshTools.H"
+#include "OFstream.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,16 +49,18 @@ void Foam::directMappedPolyPatch::collectSamples
 (
     pointField& samples,
     labelList& patchFaceProcs,
-    labelList& patchFaces
+    labelList& patchFaces,
+    pointField& patchFc
 ) const
 {
-    const vectorField::subField fc = this->faceCentres();
 
     // Collect all sample points and the faces they come from.
+    List<pointField> globalFc(Pstream::nProcs());
     List<pointField> globalSamples(Pstream::nProcs());
     labelListList globalFaces(Pstream::nProcs());
 
-    globalSamples[Pstream::myProcNo()] = fc+offset_;
+    globalFc[Pstream::myProcNo()] = this->faceCentres();
+    globalSamples[Pstream::myProcNo()] = globalFc[Pstream::myProcNo()]+offset_;
     globalFaces[Pstream::myProcNo()] = identity(size());
 
     // Distribute to all processors
@@ -64,6 +68,8 @@ void Foam::directMappedPolyPatch::collectSamples
     Pstream::scatterList(globalSamples);
     Pstream::gatherList(globalFaces);
     Pstream::scatterList(globalFaces);
+    Pstream::gatherList(globalFc);
+    Pstream::scatterList(globalFc);
 
     // Rework into straight list
     samples = ListListOps::combine<pointField>
@@ -76,6 +82,11 @@ void Foam::directMappedPolyPatch::collectSamples
         globalFaces,
         accessOp<labelList>()
     );
+    patchFc = ListListOps::combine<pointField>
+    (
+        globalFc,
+        accessOp<pointField>()
+    );
 
     patchFaceProcs.setSize(patchFaces.size());
     labelList nPerProc
@@ -103,11 +114,14 @@ void Foam::directMappedPolyPatch::findSamples
 (
     const pointField& samples,
     labelList& sampleCellProcs,
-    labelList& sampleCells
+    labelList& sampleCells,
+    pointField& sampleCc
 ) const
 {
     sampleCellProcs.setSize(samples.size());
     sampleCells.setSize(samples.size());
+    sampleCc.setSize(samples.size());
+    sampleCc = point(-GREAT, -GREAT, -GREAT);
 
     {
         // Octree based search engine
@@ -124,6 +138,8 @@ void Foam::directMappedPolyPatch::findSamples
             else
             {
                 sampleCellProcs[sampleI] = Pstream::myProcNo();
+                sampleCc[sampleI] =
+                    boundaryMesh().mesh().cellCentres()[sampleCells[sampleI]];
             }
         }
     }
@@ -136,6 +152,9 @@ void Foam::directMappedPolyPatch::findSamples
     Pstream::listCombineGather(sampleCellProcs, maxEqOp<label>());
     Pstream::listCombineScatter(sampleCellProcs);
 
+    Pstream::listCombineGather(sampleCc, maxEqOp<point>());
+    Pstream::listCombineScatter(sampleCc);
+
     if (debug)
     {
         Info<< "directMappedPolyPatch::findSamples : " << endl;
@@ -143,7 +162,8 @@ void Foam::directMappedPolyPatch::findSamples
         {
             Info<< "    " << sampleI << " coord:" << samples[sampleI]
                 << " found on processor:" << sampleCellProcs[sampleI]
-                << " in cell:" << sampleCells[sampleI] << endl;
+                << " in cell:" << sampleCells[sampleI]
+                << " with cc:" << sampleCc[sampleI] << endl;
         }
     }
 
@@ -213,12 +233,14 @@ void Foam::directMappedPolyPatch::calcMapping() const
     pointField samples;
     labelList patchFaceProcs;
     labelList patchFaces;
-    collectSamples(samples, patchFaceProcs, patchFaces);
+    pointField patchFc;
+    collectSamples(samples, patchFaceProcs, patchFaces, patchFc);
 
     // Find processor and cell samples are in
     labelList sampleCellProcs;
     labelList sampleCells;
-    findSamples(samples, sampleCellProcs, sampleCells);
+    pointField sampleCc;
+    findSamples(samples, sampleCellProcs, sampleCells, sampleCc);
 
 
     // Now we have all the data we need:
@@ -227,6 +249,60 @@ void Foam::directMappedPolyPatch::calcMapping() const
     // - cell sample is in (so source when mapping)
     //   sampleCells, sampleCellProcs.
 
+
+    if (debug && Pstream::master())
+    {
+        OFstream str
+        (
+            boundaryMesh().mesh().time().path()
+          / name()
+          + "_directMapped.obj"
+        );
+        Pout<< "Dumping mapping as lines from patch faceCentres to"
+            << " sampled cellCentres to file " << str.name() << endl;
+
+        label vertI = 0;
+
+        forAll(patchFc, i)
+        {
+            meshTools::writeOBJ(str, patchFc[i]);
+            vertI++;
+            meshTools::writeOBJ(str, sampleCc[i]);
+            vertI++;
+            str << "l " << vertI-1 << ' ' << vertI << nl;
+        }
+    }
+
+
+    // Check that actual offset vector (sampleCc - patchFc) is more or less
+    // constant.
+    if (Pstream::master())
+    {
+        const scalarField magOffset(mag(sampleCc - patchFc));
+        const scalar avgOffset(average(magOffset));
+
+        forAll(magOffset, sampleI)
+        {
+            if (mag(magOffset[sampleI]-avgOffset) > 0.001*avgOffset)
+            {
+                WarningIn("directMappedPolyPatch::calcMapping() const")
+                    << "The actual cell centres picked up using offset "
+                    << offset_ << " are not" << endl
+                    << "    on a single plane."
+                    << " This might give numerical problems." << endl
+                    << "    At patchface " << patchFc[sampleI]
+                    << " the sampled cell " << sampleCc[sampleI] << endl
+                    << "    is not on a plane " << avgOffset
+                    << " offset from the patch." << endl
+                    << "    You might want to shift your plane offset."
+                    << " Set the debug flag to get a dump of sampled cells."
+                    << endl;
+                break;
+            }
+        }
+    }
+
+
     // Determine schedule.
     mapDistribute distMap(sampleCellProcs, patchFaceProcs);
 
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
index 6439020f5d..06837c5f9b 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
@@ -82,7 +82,8 @@ class directMappedPolyPatch
         (
             pointField&,
             labelList& patchFaceProcs,
-            labelList& patchFaces
+            labelList& patchFaces,
+            pointField& patchFc
         ) const;
 
         //- Find cells containing samples
@@ -90,7 +91,8 @@ class directMappedPolyPatch
         (
             const pointField&,
             labelList& sampleCellProcs,
-            labelList& sampleCells
+            labelList& sampleCells,
+            pointField& sampleCc
         ) const;
 
         //- Calculate matching
diff --git a/src/meshTools/indexedOctree/indexedOctree.C b/src/meshTools/indexedOctree/indexedOctree.C
index f0f8bd94d0..1a55a88547 100644
--- a/src/meshTools/indexedOctree/indexedOctree.C
+++ b/src/meshTools/indexedOctree/indexedOctree.C
@@ -40,7 +40,7 @@ namespace Foam
 // Does bb intersect a sphere around sample? Or is any corner point of bb
 // closer than nearestDistSqr to sample.
 template <class Type>
-bool indexedOctree<Type>::intersects
+bool indexedOctree<Type>::overlaps
 (
     const point& p0,
     const point& p1,
@@ -84,7 +84,7 @@ bool indexedOctree<Type>::intersects
 // Does bb intersect a sphere around sample? Or is any corner point of bb
 // closer than nearestDistSqr to sample.
 template <class Type>
-bool indexedOctree<Type>::intersects
+bool indexedOctree<Type>::overlaps
 (
     const treeBoundBox& parentBb,
     const direction octant,
@@ -94,7 +94,7 @@ bool indexedOctree<Type>::intersects
 {
     //- Speeded up version of
     //     treeBoundBox subBb(parentBb.subBbox(mid, octant))
-    //     intersects
+    //     overlaps
     //     (
     //          subBb.min(),
     //          subBb.max(),
@@ -136,7 +136,7 @@ bool indexedOctree<Type>::intersects
 
     const point mid(0.5*(min+max));
 
-    return intersects(mid, other, nearestDistSqr, sample);
+    return overlaps(mid, other, nearestDistSqr, sample);
 }
 
 
@@ -567,7 +567,7 @@ void indexedOctree<Type>::findNearest
 
             const treeBoundBox& subBb = nodes_[subNodeI].bb_;
 
-            if (intersects(subBb.min(), subBb.max(), nearestDistSqr, sample))
+            if (overlaps(subBb.min(), subBb.max(), nearestDistSqr, sample))
             {
                 findNearest
                 (
@@ -584,7 +584,7 @@ void indexedOctree<Type>::findNearest
         {
             if
             (
-                intersects
+                overlaps
                 (
                     nod.bb_,
                     octant,
@@ -639,7 +639,7 @@ void indexedOctree<Type>::findNearest
         {
             const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
 
-            if (subBb.intersects(tightest))
+            if (subBb.overlaps(tightest))
             {
                 findNearest
                 (
@@ -657,7 +657,7 @@ void indexedOctree<Type>::findNearest
         {
             const treeBoundBox subBb(nodeBb.subBbox(octant));
 
-            if (subBb.intersects(tightest))
+            if (subBb.overlaps(tightest))
             {
                 shapes_.findNearest
                 (
@@ -1121,7 +1121,7 @@ void indexedOctree<Type>::findBox
         {
             const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
 
-            if (subBb.intersects(searchBox))
+            if (subBb.overlaps(searchBox))
             {
                 findBox(getNode(index), searchBox, elements);
             }
@@ -1130,7 +1130,7 @@ void indexedOctree<Type>::findBox
         {
             const treeBoundBox subBb(nodeBb.subBbox(octant));
 
-            if (subBb.intersects(searchBox))
+            if (subBb.overlaps(searchBox))
             {
                 const labelList& indices = contents_[getContent(index)];
 
diff --git a/src/meshTools/indexedOctree/indexedOctree.H b/src/meshTools/indexedOctree/indexedOctree.H
index 249eb1b5f1..97b6ab607c 100644
--- a/src/meshTools/indexedOctree/indexedOctree.H
+++ b/src/meshTools/indexedOctree/indexedOctree.H
@@ -36,9 +36,9 @@ SourceFiles
 #ifndef indexedOctree_H
 #define indexedOctree_H
 
+#include "treeBoundBox.H"
 #include "pointIndexHit.H"
 #include "FixedList.H"
-#include "treeBoundBox.H"
 #include "Ostream.H"
 #include "labelHashSet.H"
 #include "labelBits.H"
@@ -146,7 +146,7 @@ private:
 
         //- Like above but now bb is implicitly provided as parent bb + mid
         //  + octant
-        static bool intersects
+        static bool overlaps
         (
             const treeBoundBox& parentBb,
             const direction octant,
@@ -497,7 +497,7 @@ public:
 
             //- Helper: does bb intersect a sphere around sample? Or is any
             //  corner point of bb closer than nearestDistSqr to sample.
-            static bool intersects
+            static bool overlaps
             (
                 const point& bbMin,
                 const point& bbMax,
diff --git a/src/meshTools/indexedOctree/treeDataCell.C b/src/meshTools/indexedOctree/treeDataCell.C
index 2ade66780c..cb41899306 100644
--- a/src/meshTools/indexedOctree/treeDataCell.C
+++ b/src/meshTools/indexedOctree/treeDataCell.C
@@ -137,11 +137,11 @@ bool Foam::treeDataCell::overlaps
 {
     if (cacheBb_)
     {
-        return cubeBb.intersects(bbs_[index]);
+        return cubeBb.overlaps(bbs_[index]);
     }
     else
     {
-        return cubeBb.intersects(calcCellBb(cellLabels_[index]));
+        return cubeBb.overlaps(calcCellBb(cellLabels_[index]));
     }
 }
 
diff --git a/src/meshTools/indexedOctree/treeDataEdge.C b/src/meshTools/indexedOctree/treeDataEdge.C
index d918fab944..6b8d54399e 100644
--- a/src/meshTools/indexedOctree/treeDataEdge.C
+++ b/src/meshTools/indexedOctree/treeDataEdge.C
@@ -110,11 +110,11 @@ bool Foam::treeDataEdge::overlaps
 {
     if (cacheBb_)
     {
-        return cubeBb.intersects(bbs_[index]);
+        return cubeBb.overlaps(bbs_[index]);
     }
     else
     {
-        return cubeBb.intersects(calcBb(edgeLabels_[index]));
+        return cubeBb.overlaps(calcBb(edgeLabels_[index]));
     }
 }
 
diff --git a/src/meshTools/indexedOctree/treeDataFace.C b/src/meshTools/indexedOctree/treeDataFace.C
index e8d148ec23..5764841bbd 100644
--- a/src/meshTools/indexedOctree/treeDataFace.C
+++ b/src/meshTools/indexedOctree/treeDataFace.C
@@ -412,14 +412,14 @@ bool Foam::treeDataFace::overlaps
     // 1. Quick rejection: bb does not intersect face bb at all
     if (cacheBb_)
     {
-        if (!cubeBb.intersects(bbs_[index]))
+        if (!cubeBb.overlaps(bbs_[index]))
         {
             return false;
         }
     }
     else
     {
-        if (!cubeBb.intersects(calcBb(faceLabels_[index])))
+        if (!cubeBb.overlaps(calcBb(faceLabels_[index])))
         {
             return false;
         }
diff --git a/src/meshTools/indexedOctree/treeDataTriSurface.C b/src/meshTools/indexedOctree/treeDataTriSurface.C
index 23ba7f5116..c226de9499 100644
--- a/src/meshTools/indexedOctree/treeDataTriSurface.C
+++ b/src/meshTools/indexedOctree/treeDataTriSurface.C
@@ -284,13 +284,13 @@ bool Foam::treeDataTriSurface::overlaps
     triBb.max() = max(triBb.max(), p2);
 
     //- For testing: robust one
-    //return cubeBb.intersects(triBb);
+    //return cubeBb.overlaps(triBb);
 
     //- Exact test of triangle intersecting bb
 
     // Quick rejection. If whole bounding box of tri is outside cubeBb then
     // there will be no intersection.
-    if (!cubeBb.intersects(triBb))
+    if (!cubeBb.overlaps(triBb))
     {
         return false;
     }
diff --git a/src/meshTools/octree/PointIndexHit.H b/src/meshTools/octree/PointIndexHit.H
index 9b72d99605..f29965748c 100644
--- a/src/meshTools/octree/PointIndexHit.H
+++ b/src/meshTools/octree/PointIndexHit.H
@@ -94,6 +94,12 @@ public:
             index_(-1)
         {}
 
+        //- Construct from Istream
+        PointIndexHit(Istream& is)
+        {
+            is >> *this;
+        }
+
 
     // Member Functions
 
@@ -193,13 +199,45 @@ public:
 
         friend Ostream& operator<< (Ostream& os, const PointIndexHit& pHit)
         {
-            return os << pHit.hit_ << token::SPACE << pHit.hitPoint_
-                << token::SPACE << pHit.index_;
+            if (os.format() == IOstream::ASCII)
+            {
+                os  << pHit.hit_ << token::SPACE << pHit.hitPoint_
+                    << token::SPACE << pHit.index_;
+            }
+            else
+            {
+                os.write
+                (
+                    reinterpret_cast<const char*>(&pHit),
+                    sizeof(PointIndexHit)
+                );
+            }
+
+            // Check state of Ostream
+            os.check("Ostream& operator<<(Ostream&, const PointIndexHit&)");
+
+            return os;
         }
 
         friend Istream& operator>>(Istream& is, PointIndexHit& pHit)
         {
-            return is >> pHit.hit_ >> pHit.hitPoint_ >> pHit.index_;
+            if (is.format() == IOstream::ASCII)
+            {
+                return is >> pHit.hit_ >> pHit.hitPoint_ >> pHit.index_;
+            }
+            else
+            {
+                is.read
+                (
+                    reinterpret_cast<char*>(&pHit),
+                    sizeof(PointIndexHit)
+                );
+            }
+
+            // Check state of Istream
+            is.check("Istream& operator>>(Istream&, PointIndexHit&)");
+
+            return is;
         }
 
 };
diff --git a/src/meshTools/octree/octreeDataCell.C b/src/meshTools/octree/octreeDataCell.C
index 81a3600be9..17bba8b2c9 100644
--- a/src/meshTools/octree/octreeDataCell.C
+++ b/src/meshTools/octree/octreeDataCell.C
@@ -113,7 +113,7 @@ bool Foam::octreeDataCell::overlaps
     const treeBoundBox& cubeBb
 ) const
 {
-    return cubeBb.intersects(bbs_[index]);
+    return cubeBb.overlaps(bbs_[index]);
 }
 
 
diff --git a/src/meshTools/octree/octreeDataEdges.C b/src/meshTools/octree/octreeDataEdges.C
index 3974a0b91d..9828599e61 100644
--- a/src/meshTools/octree/octreeDataEdges.C
+++ b/src/meshTools/octree/octreeDataEdges.C
@@ -107,7 +107,7 @@ bool Foam::octreeDataEdges::overlaps
     const treeBoundBox& sampleBb
 ) const
 {
-    return sampleBb.intersects(allBb_[index]);
+    return sampleBb.overlaps(allBb_[index]);
 }
 
 
diff --git a/src/meshTools/octree/octreeDataFace.C b/src/meshTools/octree/octreeDataFace.C
index 46f03993cc..6a8c3f875a 100644
--- a/src/meshTools/octree/octreeDataFace.C
+++ b/src/meshTools/octree/octreeDataFace.C
@@ -507,12 +507,12 @@ bool Foam::octreeDataFace::overlaps
     const treeBoundBox& sampleBb
 ) const
 {
-    //return sampleBb.intersects(allBb_[index]);
+    //return sampleBb.overlaps(allBb_[index]);
 
     //- Exact test of face intersecting bb
 
     // 1. Quick rejection: bb does not intersect face bb at all
-    if (!sampleBb.intersects(allBb_[index]))
+    if (!sampleBb.overlaps(allBb_[index]))
     {
         return false;
     }
diff --git a/src/meshTools/octree/pointIndexHit.H b/src/meshTools/octree/pointIndexHit.H
index e41fd28499..b43f12d551 100644
--- a/src/meshTools/octree/pointIndexHit.H
+++ b/src/meshTools/octree/pointIndexHit.H
@@ -39,7 +39,14 @@ Description
 
 namespace Foam
 {
-    typedef PointIndexHit<point> pointIndexHit;
+
+typedef PointIndexHit<point> pointIndexHit;
+
+
+//- Specify data associated with pointIndexHit type is contiguous
+template<>
+inline bool contiguous<pointIndexHit>() {return true;}
+
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/meshTools/octree/treeBoundBox.C b/src/meshTools/octree/treeBoundBox.C
index 042209ec7c..26530b5a3a 100644
--- a/src/meshTools/octree/treeBoundBox.C
+++ b/src/meshTools/octree/treeBoundBox.C
@@ -297,6 +297,45 @@ Foam::treeBoundBox Foam::treeBoundBox::subBbox
 }
 
 
+bool Foam::treeBoundBox::overlaps
+(
+    const point& centre,
+    const scalar radiusSqr
+) const
+{
+    // Find out where centre is in relation to bb.
+    // Find nearest point on bb.
+    scalar distSqr = 0;
+
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        scalar d0 = min()[dir] - centre[dir];
+        scalar d1 = max()[dir] - centre[dir];
+
+        if ((d0 > 0) != (d1 > 0))
+        {
+            // centre inside both extrema. This component does not add any
+            // distance.
+        }
+        else if (Foam::mag(d0) < Foam::mag(d1))
+        {
+            distSqr += d0*d0;
+        }
+        else
+        {
+            distSqr += d1*d1;
+        }
+
+        if (distSqr > radiusSqr)
+        {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+
 // line intersection. Returns true if line (start to end) inside
 // bb or intersects bb. Sets pt to intersection.
 //
diff --git a/src/meshTools/octree/treeBoundBox.H b/src/meshTools/octree/treeBoundBox.H
index 01d7a9a6bd..02dbac38af 100644
--- a/src/meshTools/octree/treeBoundBox.H
+++ b/src/meshTools/octree/treeBoundBox.H
@@ -263,8 +263,11 @@ public:
                  FixedList<direction, 8>& octantOrder
             ) const;
 
-            //- Intersects other boundingbox?
-            inline bool intersects(const treeBoundBox&) const;
+            //- Overlaps other boundingbox?
+            inline bool overlaps(const treeBoundBox&) const;
+
+            //- Overlaps boundingSphere (centre + sqr(radius))?
+            bool overlaps(const point&, const scalar radiusSqr) const;
 
             //- Intersects segment; set point to intersection position,
             //  return true if intersection found.
diff --git a/src/meshTools/octree/treeBoundBoxI.H b/src/meshTools/octree/treeBoundBoxI.H
index 431eb8a04d..c444b2808a 100644
--- a/src/meshTools/octree/treeBoundBoxI.H
+++ b/src/meshTools/octree/treeBoundBoxI.H
@@ -24,9 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
 #include "treeBoundBox.H"
-#include "point.H"
 #include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -414,9 +412,9 @@ inline void treeBoundBox::searchOrder
 
 // true if bb's intersect or overlap.
 // Note: <= to make sure we catch all.
-inline bool treeBoundBox::intersects(const treeBoundBox& bb) const
+inline bool treeBoundBox::overlaps(const treeBoundBox& bb) const
 {
-    return boundBox::intersects(bb);
+    return boundBox::overlaps(bb);
 }
 
 
diff --git a/src/meshTools/octree/treeNode.C b/src/meshTools/octree/treeNode.C
index 678bcfb694..771ef2cdd9 100644
--- a/src/meshTools/octree/treeNode.C
+++ b/src/meshTools/octree/treeNode.C
@@ -798,7 +798,7 @@ bool treeNode<Type>::findTightest
                 // Node: recurse into subnodes
                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
 
-                if (subNodePtr->bb().intersects(tightest))
+                if (subNodePtr->bb().overlaps(tightest))
                 {
                     // there might be a better fit inside this subNode
                     changed |=
@@ -815,7 +815,7 @@ bool treeNode<Type>::findTightest
                 // Leaf: let leaf::find handle this
                 const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
 
-                if (subLeafPtr->bb().intersects(tightest))
+                if (subLeafPtr->bb().overlaps(tightest))
                 {
                     // there might be a better fit inside this subLeaf
                     changed |=
@@ -884,7 +884,7 @@ bool treeNode<Type>::findNearest
                 // Node
                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
 
-                if (subNodePtr->bb().intersects(tightest))
+                if (subNodePtr->bb().overlaps(tightest))
                 {
                     // there might be a better fit inside this subNode
                     changed |=
@@ -903,7 +903,7 @@ bool treeNode<Type>::findNearest
                 // Leaf: let leaf::find handle this
                 const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
 
-                if (subLeafPtr->bb().intersects(tightest))
+                if (subLeafPtr->bb().overlaps(tightest))
                 {
                     // there might be a better fit inside this subNode
                     changed |=
@@ -975,7 +975,7 @@ bool treeNode<Type>::findNearest
                 // Node
                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
 
-                if (subNodePtr->bb().intersects(tightest))
+                if (subNodePtr->bb().overlaps(tightest))
                 {
                     // there might be a better fit inside this subNode
                     changed |=
@@ -995,7 +995,7 @@ bool treeNode<Type>::findNearest
                 // Leaf: let leaf::find handle this
                 const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
 
-                if (subLeafPtr->bb().intersects(tightest))
+                if (subLeafPtr->bb().overlaps(tightest))
                 {
                     // there might be a better fit inside this subNode
                     changed |=
@@ -1060,7 +1060,7 @@ bool treeNode<Type>::findBox
                 // Node
                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
 
-                if (subNodePtr->bb().intersects(box))
+                if (subNodePtr->bb().overlaps(box))
                 {
                     // Visit sub node.
                     changed |= subNodePtr->findBox(shapes, box, elements);
@@ -1071,7 +1071,7 @@ bool treeNode<Type>::findBox
                 // Leaf: let leaf::find handle this
                 const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
 
-                if (subLeafPtr->bb().intersects(box))
+                if (subLeafPtr->bb().overlaps(box))
                 {
                     // Visit sub leaf.
                     changed |= subLeafPtr->findBox(shapes, box, elements);
diff --git a/src/meshTools/searchableSurface/searchableSphere.C b/src/meshTools/searchableSurface/searchableSphere.C
index 47abfbfff5..04f3f2f2b5 100644
--- a/src/meshTools/searchableSurface/searchableSphere.C
+++ b/src/meshTools/searchableSurface/searchableSphere.C
@@ -26,7 +26,6 @@ License
 
 #include "searchableSphere.H"
 #include "addToRunTimeSelectionTable.H"
-#include "SortableList.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
diff --git a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C
index f31b30602b..7e716c1679 100644
--- a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C
+++ b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C
@@ -27,8 +27,6 @@ License
 #include "intersectedSurface.H"
 #include "surfaceIntersection.H"
 #include "faceList.H"
-#include "SortableList.H"
-#include "triSurfaceTools.H"
 #include "faceTriangulation.H"
 #include "treeBoundBox.H"
 #include "OFstream.H"
diff --git a/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C b/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C
index 6a21879f35..f3dd4f3e03 100644
--- a/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C
+++ b/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C
@@ -385,12 +385,12 @@ bool Foam::octreeDataTriSurface::overlaps
     const treeBoundBox& cubeBb
 ) const
 {
-    //return cubeBb.intersects(allBb_[index]);
+    //return cubeBb.overlaps(allBb_[index]);
 
     //- Exact test of triangle intersecting bb
 
     // Quick rejection.
-    if (!cubeBb.intersects(allBb_[index]))
+    if (!cubeBb.overlaps(allBb_[index]))
     {
         return false;
     }
diff --git a/src/meshTools/triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.C b/src/meshTools/triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.C
index 6582cfd239..ca4a34975c 100644
--- a/src/meshTools/triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.C
+++ b/src/meshTools/triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.C
@@ -54,7 +54,7 @@ bool Foam::treeLeaf<Foam::octreeDataTriSurface>::findNearest
         label faceI = indices_[i];
 
         // Quick rejection test.
-        if (tightest.intersects(allBb[faceI]))
+        if (tightest.overlaps(allBb[faceI]))
         {
             // Full calculation
             scalar dist = shapes.calcNearest(faceI, sample, nearest);
diff --git a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
index 56f8799aec..e8c9a29552 100644
--- a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
+++ b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
@@ -101,6 +101,7 @@ bool Foam::triangleFuncs::intersectAxesBundle
     // Since direction is coordinate axis there is no need to do projection,
     // we can directly check u,v components for inclusion in triangle.
 
+    scalar localScale = max(max(magSqr(V10), magSqr(V20)), 1.0);
 
     // Get other components
     label i1 = (i0 + 1) % 3;
@@ -114,8 +115,11 @@ bool Foam::triangleFuncs::intersectAxesBundle
 
     scalar det = v2*u1 - u2*v1;
 
-    // Fix for V0:(-31.71428 0 -15.10714) V10:(-1.285715 8.99165e-16 -1.142858) V20:(0 0 -1.678573) i0:0
-    if (Foam::mag(det)/max(max(mag(V10),mag(V20)),1) < SMALL)
+    // Fix for  V0:(-31.71428 0 -15.10714)
+    //          V10:(-1.285715 8.99165e-16 -1.142858)
+    //          V20:(0 0 -1.678573)
+    //          i0:0
+    if (Foam::mag(det)/localScale < SMALL)
     {
         // Triangle parallel to dir
         return false;
@@ -132,7 +136,7 @@ bool Foam::triangleFuncs::intersectAxesBundle
         scalar beta = 0;
         bool inter = false;
 
-        if (Foam::mag(u1) < SMALL)
+        if (Foam::mag(u1)/localScale < SMALL)
         {
             beta = u0/u2;
             if ((beta >= 0) && (beta <= 1))
diff --git a/src/triSurface/triSurface/triSurface.C b/src/triSurface/triSurface/triSurface.C
index adc28176c5..89a48a30d6 100644
--- a/src/triSurface/triSurface/triSurface.C
+++ b/src/triSurface/triSurface/triSurface.C
@@ -31,6 +31,7 @@ License
 #include "Time.H"
 #include "boundBox.H"
 #include "SortableList.H"
+#include "PackedList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -626,7 +627,7 @@ surfacePatchList triSurface::calcPatches(labelList& faceMap) const
     {
         sortedRegion[faceI] = operator[](faceI).region();
     }
-    sortedRegion.stableSort();
+    sortedRegion.sort();
 
     faceMap = sortedRegion.indices();
 
@@ -743,6 +744,26 @@ triSurface::triSurface
 {}
 
 
+triSurface::triSurface
+(
+    List<labelledTri>& triangles,
+    const geometricSurfacePatchList& patches,
+    pointField& points,
+    const bool reUse
+)
+:
+    PrimitivePatch<labelledTri, ::Foam::List, pointField>
+    (
+        triangles,
+        points,
+        reUse
+    ),
+    patches_(patches),
+    sortedEdgeFacesPtr_(NULL),
+    edgeOwnerPtr_(NULL)
+{}
+
+
 triSurface::triSurface
 (
     const List<labelledTri>& triangles,
@@ -1148,9 +1169,7 @@ triSurface triSurface::subsetMesh
     }
 
     // Construct subsurface
-    triSurface subSurface(newTriangles, patches(), newPoints);
-
-    return subSurface;
+    return triSurface(newTriangles, patches(), newPoints, true);
 }
 
 
@@ -1187,30 +1206,36 @@ void triSurface::write(const Time& d) const
 
 void triSurface::writeStats(Ostream& os) const
 {
-    // Calculate bounding box without any additional addressing
-    // Copy of treeBoundBox code. Cannot use meshTools from triSurface...
+    // Unfortunately nPoints constructs meshPoints() so do compact version
+    // ourselves.
+    PackedList<1> pointIsUsed(points().size());
+    pointIsUsed = 0U;
+
+    label nPoints = 0;
     boundBox bb
     (
         point(VGREAT, VGREAT, VGREAT),
         point(-VGREAT, -VGREAT, -VGREAT)
     );
+
     forAll(*this, triI)
     {
         const labelledTri& f = operator[](triI);
 
         forAll(f, fp)
         {
-            const point& pt = points()[f[fp]];
-            bb.min() = ::Foam::min(bb.min(), pt);
-            bb.max() = ::Foam::max(bb.max(), pt);
+            label pointI = f[fp];
+            if (pointIsUsed.set(pointI, 1))
+            {
+                bb.min() = ::Foam::min(bb.min(), points()[pointI]);
+                bb.max() = ::Foam::max(bb.max(), points()[pointI]);
+                nPoints++;
+            }
         }
     }
 
-    // Unfortunately nPoints constructs meshPoints() ...
-
     os  << "Triangles    : " << size() << endl
-        //<< "Edges        : " << nEdges() << endl
-        << "Vertices     : " << nPoints() << endl
+        << "Vertices     : " << nPoints << endl
         << "Bounding Box : " << bb << endl;
 }
 
diff --git a/src/triSurface/triSurface/triSurface.H b/src/triSurface/triSurface/triSurface.H
index 56b1abc849..624b793117 100644
--- a/src/triSurface/triSurface/triSurface.H
+++ b/src/triSurface/triSurface/triSurface.H
@@ -36,9 +36,9 @@ SourceFiles
 #ifndef triSurface_H
 #define triSurface_H
 
+#include "PrimitivePatch.H"
 #include "pointField.H"
 #include "labelledTri.H"
-#include "PrimitivePatch.H"
 #include "boolList.H"
 #include "geometricSurfacePatchList.H"
 #include "surfacePatchList.H"
@@ -215,6 +215,15 @@ public:
             const pointField&
         );
 
+        //- Construct from triangles, patches, points. Reuse storage.
+        triSurface
+        (
+            List<labelledTri>&,
+            const geometricSurfacePatchList&,
+            pointField&,
+            const bool reUse
+        );
+
         //- Construct from triangles, points. Set patchnames to default.
         triSurface(const List<labelledTri>&, const pointField&);
 
diff --git a/tutorials/oodles/pitzDailyDirectMapped/system/changeDictionaryDict b/tutorials/oodles/pitzDailyDirectMapped/system/changeDictionaryDict
index efb0bc6379..282852ed87 100644
--- a/tutorials/oodles/pitzDailyDirectMapped/system/changeDictionaryDict
+++ b/tutorials/oodles/pitzDailyDirectMapped/system/changeDictionaryDict
@@ -22,7 +22,7 @@ dictionaryReplacement
         inlet
         {
             type            directMappedPatch;
-            offset          (0.05 0 0);
+            offset          (0.0495 0 0);
         }
     }
 }