diff --git a/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C b/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C
index 264eb28dea..8d109f03b3 100644
--- a/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C
+++ b/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C
@@ -85,6 +85,8 @@ int main(int argc, char *argv[])
 
         #include "continuityErrs.H"
 
+        p = pd + rho*gh;
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/utilities/mesh/advanced/refineWallLayer/refineWallLayer.C b/applications/utilities/mesh/advanced/refineWallLayer/refineWallLayer.C
index 61a2f4d854..f4336bce92 100644
--- a/applications/utilities/mesh/advanced/refineWallLayer/refineWallLayer.C
+++ b/applications/utilities/mesh/advanced/refineWallLayer/refineWallLayer.C
@@ -184,10 +184,6 @@ int main(int argc, char *argv[])
         << endl;
 
     // Transfer DynamicLists to straight ones.
-    labelList cutEdges;
-    cutEdges.transfer(allCutEdges);
-    allCutEdges.clear();
-
     scalarField cutEdgeWeights;
     cutEdgeWeights.transfer(allCutEdgeWeights);
     allCutEdgeWeights.clear();
@@ -199,7 +195,7 @@ int main(int argc, char *argv[])
         mesh,
         cutCells.toc(),     // cells candidate for cutting
         labelList(0),       // cut vertices
-        cutEdges,           // cut edges
+        allCutEdges,        // cut edges
         cutEdgeWeights      // weight on cut edges
     );
 
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 4737f38e53..edcfbf2aee 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -433,7 +433,7 @@ int main(int argc, char *argv[])
     Info<< "Finished meshing in = "
         << runTime.elapsedCpuTime() << " s." << endl;
 
-    Pout<< "End\n" << endl;
+    Info<< "End\n" << endl;
 
     return(0);
 }
diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposeMesh.C b/applications/utilities/parallelProcessing/decomposePar/decomposeMesh.C
index 0d81b1aedd..55f5c96579 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposeMesh.C
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposeMesh.C
@@ -35,11 +35,8 @@ Description
 #include "IOstreams.H"
 #include "SLPtrList.H"
 #include "boolList.H"
-#include "cellList.H"
 #include "primitiveMesh.H"
-#include "cyclicFvPatch.H"
-#include "fvPatchList.H"
-#include "DynamicList.H"
+#include "cyclicPolyPatch.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
index bee2d17b25..dd107c2f2f 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
@@ -33,7 +33,7 @@ License
 #include "OSspecific.H"
 #include "Map.H"
 #include "globalMeshData.H"
-
+#include "DynamicList.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
diff --git a/applications/utilities/postProcessing/dataConversion/foamToFieldview9/foamToFieldview9.C b/applications/utilities/postProcessing/dataConversion/foamToFieldview9/foamToFieldview9.C
index 3030a30316..7cf137cd0e 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToFieldview9/foamToFieldview9.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToFieldview9/foamToFieldview9.C
@@ -45,7 +45,6 @@ Description
 #include "IOobjectList.H"
 #include "boolList.H"
 #include "stringList.H"
-#include "DynamicList.H"
 #include "cellModeller.H"
 
 #include "floatScalar.H"
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/internalWriter.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/internalWriter.C
index bf1fb3b775..1ac0d8b729 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/internalWriter.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/internalWriter.C
@@ -27,9 +27,6 @@ License
 #include "internalWriter.H"
 #include "writeFuns.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriter.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriter.C
index bbf24b1b17..ea7da50a40 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriter.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriter.C
@@ -29,9 +29,6 @@ License
 #include "Cloud.H"
 #include "passiveParticle.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writePatchGeom.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/writePatchGeom.C
index 62b8661810..6ad993e0aa 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/writePatchGeom.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/writePatchGeom.C
@@ -78,6 +78,7 @@ void writePatchGeom
     writeFuns::write(pStream, binary, vertLabels);
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writePointSet.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/writePointSet.C
index ea7dbf4301..9e495e073d 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/writePointSet.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/writePointSet.C
@@ -103,6 +103,7 @@ void writePointSet
     writeFuns::write(pStream, binary, pointIDs);
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.C b/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.C
index 8297002b82..24d1edea53 100644
--- a/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.C
+++ b/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.C
@@ -33,7 +33,7 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList(Istream& is)
 :
     List<T>(is),
-    nextFree_(List<T>::size())
+    allocSize_(List<T>::size())
 {}
 
 
@@ -44,9 +44,6 @@ Foam::Ostream& Foam::operator<<
     const Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>& DL
 )
 {
-    const_cast<DynamicList<T, SizeInc, SizeMult, SizeDiv>&>(DL)
-        .setSize(DL.nextFree_);
-
     os << static_cast<const List<T>&>(DL);
     return os;
 }
@@ -60,7 +57,7 @@ Foam::Istream& Foam::operator>>
 )
 {
     is >> static_cast<List<T>&>(DL);
-    DL.nextFree_ = DL.List<T>::size();
+    DL.allocSize_ = DL.List<T>::size();
 
     return is;
 }
diff --git a/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H b/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H
index a98be44fef..684d393745 100644
--- a/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H
+++ b/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H
@@ -81,24 +81,29 @@ class DynamicList
 {
     // Private data
 
-        //- Number of next free element
-        label nextFree_;
+        //- Allocated size for underlying List.
+        label allocSize_;
 
 
 public:
 
+    // Related types
+
+        //- Declare friendship with the List class
+        friend class List<T>;
+
     // Constructors
 
         //- Construct null
         inline DynamicList();
 
-        //- Construct given size
+        //- Construct given size.
         explicit inline DynamicList(const label);
 
-        //- Construct from UList. nextFree_ set to size().
+        //- Construct from UList. Size set to UList size.
         explicit inline DynamicList(const UList<T>&);
 
-        //- Construct from Istream. nextFree_ set to size().
+        //- Construct from Istream. Size set to size of read list.
         explicit DynamicList(Istream&);
 
 
@@ -106,22 +111,24 @@ public:
 
         // Access
 
-            //- Size of the active part of the list.
-            //  Direct over-ride of list size member function
-            inline label size() const;
+            //- Size of the underlying storage.
+            inline label allocSize() const;
 
 
         // Edit
 
             //- Reset size of List.
-            void setSize(const label);
+            inline void setSize(const label);
 
             //- Reset size of List and value for new elements.
-            void setSize(const label, const T&);
+            inline void setSize(const label, const T&);
 
-            //- Clear the list, i.e. set next free to zero.
+            //- Clear the list, i.e. set the size to zero.
             //  Allocated size does not change
-            void clear();
+            inline void clear();
+
+            //- Clear the list and delete storage.
+            inline void clearStorage();
 
             //- Shrink the List<T> to the number of elements used
             inline DynamicList<T, SizeInc, SizeMult, SizeDiv>& shrink();
@@ -130,11 +137,11 @@ public:
             //  and annull the argument list. Is same as List::transfer except
             //  checks that you're not changing the underlying list to something
             //  smaller than nextFree_.
-            void transfer(List<T>&);
+            inline void transfer(List<T>&);
 
             //- Transfer the contents of the argument DynamicList into this
             //  DynamicList and annull the argument list.
-            void transfer(DynamicList<T, SizeInc, SizeMult, SizeDiv>&);
+            inline void transfer(DynamicList<T, SizeInc, SizeMult, SizeDiv>&);
 
 
     // Member Operators
@@ -152,9 +159,15 @@ public:
         //- Assignment of all entries to the given value
         inline void operator=(const T&);
 
-        //- Assignment to List<T>
+        //- Assignment from List<T>
         inline void operator=(const List<T>&);
 
+        //- Assignment from DynamicList<T>
+        inline void operator=
+        (
+            const DynamicList<T, SizeInc, SizeMult, SizeDiv>&
+        );
+
 
     // IOstream operators
 
diff --git a/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H b/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
index 8d451ccd7a..d9a14a1111 100644
--- a/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
+++ b/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
@@ -31,8 +31,10 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList()
 :
     List<T>(SizeInc),
-    nextFree_(0)
-{}
+    allocSize_(SizeInc)
+{
+    List<T>::size() = 0;
+}
 
 
 //- Construct given size
@@ -43,8 +45,10 @@ inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList
 )
 :
     List<T>(s),
-    nextFree_(0)
-{}
+    allocSize_(s)
+{
+    List<T>::size() = 0;
+}
 
 
 //- Construct given size
@@ -55,17 +59,17 @@ inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList
 )
 :
     List<T>(s),
-    nextFree_(s.size())
+    allocSize_(s.size())
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
-inline Foam::label Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::size()
+inline Foam::label Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::allocSize()
 const
 {
-    return nextFree_;
+    return allocSize_;
 }
 
 
@@ -75,13 +79,14 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
     const label s
 )
 {
-    if (s < nextFree_)
+    if (s < List<T>::size())
     {
-        nextFree_ = s;
+        List<T>::size() = s;
     }
     else
     {
-        List<T>::setSize(s);
+        allocSize_ = s;
+        List<T>::setSize(allocSize_);
     }
 }
 
@@ -93,14 +98,14 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
     const T& t
 )
 {
-    if (s < nextFree_)
+    if (s < List<T>::size())
     {
-        nextFree_ = s;
+        List<T>::size() = s;
     }
     else
     {
-        List<T>::setSize(s, t);
-        nextFree_ = s;
+        allocSize_ = s;
+        List<T>::setSize(allocSize_, t);
     }
 }
 
@@ -108,7 +113,16 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::clear()
 {
-    nextFree_ = 0;
+    List<T>::size() = 0;
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::clearStorage()
+{
+    List<T>::size() = allocSize_;    // make List<T> consistent
+    List<T>::clear();
+    allocSize_ = 0;
 }
 
 
@@ -116,7 +130,8 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>&
 Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::shrink()
 {
-    List<T>::setSize(nextFree_);
+    allocSize_ = List<T>::size();
+    List<T>::setSize(allocSize_);
     return *this;
 }
 
@@ -125,20 +140,20 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline void
 Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::transfer(List<T>& l)
 {
-    if (l.size() < nextFree_)
+    if (l.size() < List<T>::size())
     {
         FatalErrorIn
         (
             "void DynamicList<T, SizeInc, SizeMult"
             ", SizeDiv>::transfer(List<T>&)"
         )   << "Cannot replace the underlying storage of this DynamicList"
-            << " of which " << nextFree_ << " elements are used" << nl
+            << " of which " << List<T>::size() << " elements are used" << nl
             << "with a List of size " << l.size() << abort(FatalError);
     }
     else
     {
+        allocSize_ = l.size();
         List<T>::transfer(l);   // take over storage
-        l.clear();              // set nextFree of l to 0
     }
 }
 
@@ -150,40 +165,40 @@ Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::transfer
     DynamicList<T, SizeInc, SizeMult, SizeDiv>& l
 )
 {
+    allocSize_ = l.allocSize();
     List<T>::transfer(l);       // take over storage
-    nextFree_ = l.size();       // take over used size
-    l.clear();                  // set nextFree of l to 0
+    l.allocSize_ = 0;
 }
 
 
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
-inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::append
-(
-    const T& e
-)
+inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::append(const T& e)
 {
-    nextFree_++;
+    // Work on copy free index since gets overwritten by setSize
+    label nextFree = List<T>::size();
 
-    if (nextFree_ > List<T>::size())
+    nextFree++;
+
+    if (nextFree > allocSize_)
     {
-        List<T>::setSize
+        allocSize_ = max
         (
-            max
-            (
-                nextFree_,
-                label(SizeMult*List<T>::size()/SizeDiv + SizeInc)
-            )
+            nextFree,
+            label(SizeMult*allocSize_/SizeDiv + SizeInc)
         );
+        List<T>::setSize(allocSize_);
     }
 
-    this->operator[](nextFree_ - 1) = e;
+    this->operator[](nextFree - 1) = e;
+
+    List<T>::size() = nextFree;
 }
 
 
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline T Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::remove()
 {
-    if (nextFree_ == 0)
+    if (List<T>::size() == 0)
     {
         FatalErrorIn
         (
@@ -191,7 +206,7 @@ inline T Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::remove()
         )   << "List is empty" << abort(FatalError);
     }
 
-    return List<T>::operator[](--nextFree_);
+    return List<T>::operator[](--List<T>::size());
 }
 
 
@@ -203,20 +218,22 @@ inline T& Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator()
     const label i
 )
 {
-    nextFree_ = max(nextFree_, i + 1);
+    label nextFree = List<T>::size();
 
-    if (nextFree_ > List<T>::size())
+    nextFree = max(nextFree, i + 1);
+
+    if (nextFree > allocSize_)
     {
-        List<T>::setSize
+       allocSize_ = max
         (
-            max
-            (
-                nextFree_,
-                label(SizeMult*List<T>::size()/SizeDiv + SizeInc)
-            )
+            nextFree,
+            label(SizeMult*allocSize_/SizeDiv + SizeInc)
         );
+        List<T>::setSize(allocSize_);
     }
 
+    List<T>::size() = nextFree;
+
     return this->operator[](i);
 }
 
@@ -228,7 +245,7 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator=
 )
 {
     List<T>::operator=(t);
-    nextFree_ = List<T>::size();
+    allocSize_ = List<T>::size();
 }
 
 
@@ -239,7 +256,18 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator=
 )
 {
     List<T>::operator=(l);
-    nextFree_ = l.size();
+    allocSize_ = List<T>::size();
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator=
+(
+    const DynamicList<T, SizeInc, SizeMult, SizeDiv>& l
+)
+{
+    List<T>::operator=(l);
+    allocSize_ = l.allocSize();
 }
 
 
diff --git a/src/OpenFOAM/containers/Lists/List/List.C b/src/OpenFOAM/containers/Lists/List/List.C
index 249d1fcf2b..21cc042050 100644
--- a/src/OpenFOAM/containers/Lists/List/List.C
+++ b/src/OpenFOAM/containers/Lists/List/List.C
@@ -420,6 +420,23 @@ void List<T>::transfer(List<T>& a)
 }
 
 
+// Transfer the contents of the argument DynamicList into this List
+// and anull the argument list
+template<class T>
+template<unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+void List<T>::transfer(DynamicList<T, SizeInc, SizeMult, SizeDiv>& a)
+{
+    if (this->v_) delete[] this->v_;
+
+    this->size_ = a.size_;
+    this->v_ = a.v_;
+
+    a.size_ = 0;
+    a.allocSize_ = 0;
+    a.v_ = 0;
+}
+
+
 template<class T>
 void sort(List<T>& a)
 {
diff --git a/src/OpenFOAM/containers/Lists/List/List.H b/src/OpenFOAM/containers/Lists/List/List.H
index 3ee814a0f6..7269159ef7 100644
--- a/src/OpenFOAM/containers/Lists/List/List.H
+++ b/src/OpenFOAM/containers/Lists/List/List.H
@@ -61,6 +61,8 @@ template<class T> Istream& operator>>(Istream&, List<T>&);
 template<class T, label Size> class FixedList;
 template<class T> class PtrList;
 template<class T> class SLList;
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+ class DynamicList;
 template<class T> class IndirectList;
 template<class T> class BiIndirectList;
 
@@ -156,6 +158,11 @@ public:
             //  and annull the argument list.
             void transfer(List<T>&);
 
+            //- Transfer the contents of the argument List into this List
+            //  and annull the argument list.
+            template<unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+            void transfer(DynamicList<T, SizeInc, SizeMult, SizeDiv>&);
+
             //- Return subscript-checked element of UList.
             inline T& newElmt(const label);
 
diff --git a/src/OpenFOAM/containers/Lists/PackedList/PackedList.C b/src/OpenFOAM/containers/Lists/PackedList/PackedList.C
index ab98111778..7c20b200b3 100644
--- a/src/OpenFOAM/containers/Lists/PackedList/PackedList.C
+++ b/src/OpenFOAM/containers/Lists/PackedList/PackedList.C
@@ -40,10 +40,7 @@ PackedList<nBits>::PackedList(const label size, const unsigned int val)
     List<unsigned int>(intSize(size)),
     size_(size)
 {
-    for (label i = 0; i < size; i++)
-    {
-        set(i, val);
-    }
+    operator=(val);
 }
 
 
diff --git a/src/OpenFOAM/containers/Lists/PackedList/PackedList.H b/src/OpenFOAM/containers/Lists/PackedList/PackedList.H
index 8eb84b9749..5a966c6e68 100644
--- a/src/OpenFOAM/containers/Lists/PackedList/PackedList.H
+++ b/src/OpenFOAM/containers/Lists/PackedList/PackedList.H
@@ -160,13 +160,15 @@ public:
             //- Number of packed elements
             inline label size() const;
 
-
             //- Get value at index I
             inline unsigned int get(const label i) const;
 
             //- Set value at index I. Return true if value changed.
             inline bool set(const label i, const unsigned int val);
 
+            //- Underlying storage
+            inline List<unsigned int>& storage();
+
 
     // Member operators
 
@@ -179,7 +181,8 @@ public:
             //- Assignment operator. Takes linear time.
             void operator=(const PackedList<nBits>&);
 
-            //- Assignment of all entries to the given value
+            //- Assignment of all entries to the given value. Does set on all
+            //  elements.
             inline void operator=(const unsigned int val);
 
             //- Return as labelList
diff --git a/src/OpenFOAM/containers/Lists/PackedList/PackedListI.H b/src/OpenFOAM/containers/Lists/PackedList/PackedListI.H
index 42b7afb81a..8e24a291bc 100644
--- a/src/OpenFOAM/containers/Lists/PackedList/PackedListI.H
+++ b/src/OpenFOAM/containers/Lists/PackedList/PackedListI.H
@@ -83,7 +83,7 @@ inline void PackedList<nBits>::checkValue(const unsigned int val) const
 {
     if (val>=(1u << nBits))
     {
-        FatalErrorIn("PackedList<T>::set(const unsigned int)")
+        FatalErrorIn("PackedList<T>::checkValue(const unsigned int)")
             << "value " << label(val) << " out of range 0 ... "
             << label((1u << nBits)-1)
             << " representable by " << nBits << " bits"
@@ -183,6 +183,13 @@ inline bool PackedList<nBits>::set(const label i, const unsigned int val)
 }
 
 
+template<int nBits>
+inline List<unsigned int>& PackedList<nBits>::storage()
+{
+    return static_cast<List<unsigned int>&>(*this);
+}
+
+
 template<int nBits>
 inline ::Foam::reference PackedList<nBits>::operator[](const label i)
 {
@@ -213,9 +220,21 @@ inline void PackedList<nBits>::operator=(const unsigned int val)
 #   ifdef DEBUGList
     checkValue(val);
 #   endif
-    List<unsigned int>::operator=(val);
+
+    if (val == 0)
+    {
+        List<unsigned int>::operator=(val);
+    }
+    else
+    {
+        for (label i = 0; i < size_; i++)
+        {
+            set(i, val);
+        }
+    }
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/containers/Lists/PackedList/PackedListIO.C b/src/OpenFOAM/containers/Lists/PackedList/PackedListIO.C
deleted file mode 100644
index dfa961a873..0000000000
--- a/src/OpenFOAM/containers/Lists/PackedList/PackedListIO.C
+++ /dev/null
@@ -1,152 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "PackedList.H"
-#include "Ostream.H"
-#include "token.H"
-#include "contiguous.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * Ostream Operator *  * * * * * * * * * * * * //
-
-template<class T>
-void PackedList<T>::writeEntry(Ostream& os) const
-{
-    if
-    (
-        size()
-     && token::compound::isCompound
-        (
-            "List<" + word(pTraits<T>::typeName) + '>'
-        )
-    )
-    {
-        os  << word("List<" + word(pTraits<T>::typeName) + '>') << " ";
-    }
-    
-    os << *this;
-}
-
-
-template<class T>
-void PackedList<T>::writeEntry(const word& keyword, Ostream& os) const
-{
-    os.writeKeyword(keyword);
-    writeEntry(os);
-    os << token::END_STATEMENT << endl;
-}
-
-
-template<class T>
-Ostream& operator<<(Ostream& os, const PackedList<T>& L)
-{
-    // Write list contents depending on data format
-    if (os.format() == IOstream::ASCII || !contiguous<T>())
-    {
-        bool uniform = false;
-
-        if (L.size() > 1 && contiguous<T>())
-        {
-            uniform = true;
-
-            unsigned int L0 = L.get(0);
-
-            forAll(L, i)
-            {
-                if (L.get(i) != L0)
-                {
-                    uniform = false;
-                    break;
-                }
-            }
-        }
-
-        if (uniform)
-        {
-            // Write size of list and start contents delimiter
-            os << L.size() << token::BEGIN_BLOCK;
-
-            // Write list contents
-            os << L.get(0);
-
-            // Write end of contents delimiter
-            os << token::END_BLOCK;
-        }
-        else if (L.size() < 11 && contiguous<T>())
-        {
-            // Write size of list and start contents delimiter
-            os << L.size() << token::BEGIN_LIST;
-
-            // Write list contents
-            forAll(L, i)
-            {
-                if (i > 0) os << token::SPACE;
-                os << L.get(i);
-            }
-
-            // Write end of contents delimiter
-            os << token::END_LIST;
-        }
-        else
-        {
-            // Write size of list and start contents delimiter
-            os << nl << L.size() << nl << token::BEGIN_LIST;
-
-            // Write list contents
-            forAll(L, i)
-            {
-                os << nl << L.get(i);
-            }
-
-            // Write end of contents delimiter
-            os << nl << token::END_LIST << nl;
-        }
-    }
-    else
-    {
-        os << nl << L.size() << nl;
-        if (L.size())
-        {
-            os.write((const char*)L.v_, L.size()*sizeof(T));
-        }
-    }
-
-    // Check state of IOstream
-    os.check("Ostream& operator<<(Ostream&, const PackedList&)");
-
-    return os;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.C b/src/OpenFOAM/meshes/boundBox/boundBox.C
index c87e552c85..a4b82f9dc4 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBox.C
+++ b/src/OpenFOAM/meshes/boundBox/boundBox.C
@@ -79,10 +79,9 @@ boundBox::boundBox(const pointField& points, const bool doReduce)
 
 
 boundBox::boundBox(Istream& is)
-:
-    min_(is),
-    max_(is)
-{}
+{
+    operator>>(is, *this);
+}
 
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
@@ -113,7 +112,7 @@ Istream& operator>>(Istream& is, boundBox& bb)
 {
     if (is.format() == IOstream::ASCII)
     {
-        return is >> bb.min_ >>  bb.max_;
+        return is >> bb.min_ >> bb.max_;
     }
     else
     {
diff --git a/src/OpenFOAM/meshes/meshShapes/edge/edge.H b/src/OpenFOAM/meshes/meshShapes/edge/edge.H
index 4e481c24f1..8b37a9f3be 100644
--- a/src/OpenFOAM/meshes/meshShapes/edge/edge.H
+++ b/src/OpenFOAM/meshes/meshShapes/edge/edge.H
@@ -141,6 +141,9 @@ inline label Hash<edge>::operator()(const edge& e) const
     return e[0]*e[1] + e[0]+e[1];
 }
 
+template<>
+inline bool contiguous<edge>()  {return true;}
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H b/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H
index 38301fec58..48b79eebe4 100644
--- a/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H
+++ b/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H
@@ -166,6 +166,9 @@ inline label Hash<triFace>::operator()(const triFace& t) const
     return (t[0]*t[1]*t[2] + t[0]+t[1]+t[2]);
 }
 
+template<>
+inline bool contiguous<triFace>()  {return true;}
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C
index 2760104238..eda9946856 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C
@@ -26,12 +26,6 @@ License
 
 #include "globalIndex.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::globalIndex::globalIndex(const label localSize)
@@ -68,15 +62,6 @@ Foam::globalIndex::globalIndex(Istream& is)
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 Foam::Istream& Foam::operator>>(Istream& is, globalIndex& gi)
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
index df8d61fa21..8fdee370a0 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
@@ -397,13 +397,11 @@ void Foam::globalMeshData::calcSharedEdges() const
     sharedEdgeLabelsPtr_ = new labelList();
     labelList& sharedEdgeLabels = *sharedEdgeLabelsPtr_;
     sharedEdgeLabels.transfer(dynSharedEdgeLabels);
-    dynSharedEdgeLabels.clear();
 
     dynSharedEdgeAddr.shrink();
     sharedEdgeAddrPtr_ = new labelList();
     labelList& sharedEdgeAddr = *sharedEdgeAddrPtr_;
     sharedEdgeAddr.transfer(dynSharedEdgeAddr);
-    dynSharedEdgeAddr.clear();
 
     if (debug)
     {
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
index 6ad85c7971..86a0a618c1 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
@@ -84,16 +84,32 @@ void Foam::mapDistribute::distribute
         // Receive sub field from neighbour
         for (label domain = 0; domain < Pstream::nProcs(); domain++)
         {
-            if
-            (
-                domain != Pstream::myProcNo()
-             && constructMap[domain].size() > 0
-            )
+            const labelList& map = constructMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
             {
                 IPstream fromNbr(Pstream::blocking, domain);
                 List<T> subField(fromNbr);
 
-                const labelList& map = constructMap[domain];
+                if (subField.size() != map.size())
+                {
+                    FatalErrorIn
+                    (
+                        "template<class T>\n"
+                        "void mapDistribute::distribute\n"
+                        "(\n"
+                        "    const Pstream::commsTypes commsType,\n"
+                        "    const List<labelPair>& schedule,\n"
+                        "    const label constructSize,\n"
+                        "    const labelListList& subMap,\n"
+                        "    const labelListList& constructMap,\n"
+                        "    List<T>& field\n"
+                        ")\n"
+                    )   << "Expected from processor " << domain
+                        << " " << map.size() << " but received "
+                        << subField.size() << " elements."
+                        << abort(FatalError);
+                }
 
                 forAll(map, i)
                 {
@@ -126,6 +142,7 @@ void Foam::mapDistribute::distribute
             newField[map[i]] = subField[i];
         }
 
+        // Schedule will already have pruned 0-sized comms
         forAll(schedule, i)
         {
             const labelPair& twoProcs = schedule[i];
@@ -154,6 +171,26 @@ void Foam::mapDistribute::distribute
 
                 const labelList& map = constructMap[sendProc];
 
+                if (subField.size() != map.size())
+                {
+                    FatalErrorIn
+                    (
+                        "template<class T>\n"
+                        "void mapDistribute::distribute\n"
+                        "(\n"
+                        "    const Pstream::commsTypes commsType,\n"
+                        "    const List<labelPair>& schedule,\n"
+                        "    const label constructSize,\n"
+                        "    const labelListList& subMap,\n"
+                        "    const labelListList& constructMap,\n"
+                        "    List<T>& field\n"
+                        ")\n"
+                    )   << "Expected from processor " << sendProc
+                        << " " << map.size() << " but received "
+                        << subField.size() << " elements."
+                        << abort(FatalError);
+                }
+
                 forAll(map, i)
                 {
                     newField[map[i]] = subField[i];
@@ -164,7 +201,74 @@ void Foam::mapDistribute::distribute
     }
     else if (commsType == Pstream::nonBlocking)
     {
-        List<T> newField(constructSize);
+        if (!contiguous<T>())
+        {
+            FatalErrorIn
+            (
+                "template<class T>\n"
+                "void mapDistribute::distribute\n"
+                "(\n"
+                "    const Pstream::commsTypes commsType,\n"
+                "    const List<labelPair>& schedule,\n"
+                "    const label constructSize,\n"
+                "    const labelListList& subMap,\n"
+                "    const labelListList& constructMap,\n"
+                "    List<T>& field\n"
+                ")\n"
+            )   << "Non-blocking only supported for contiguous data."
+                << exit(FatalError);
+        }
+
+        // Set up sends to neighbours
+
+        List<List<T > > sendFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                List<T>& subField = sendFields[domain];
+                subField.setSize(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = field[map[i]];
+                }
+
+                OPstream::write
+                (
+                    Pstream::nonBlocking,
+                    domain,
+                    reinterpret_cast<const char*>(subField.begin()),
+                    subField.size()
+                );
+            }
+        }
+
+        // Set up receives from neighbours
+
+        List<List<T > > recvFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = constructMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                recvFields[domain].setSize(map.size());
+                IPstream::read
+                (
+                    Pstream::nonBlocking,
+                    domain,
+                    reinterpret_cast<char*>(recvFields[domain].begin()),
+                    recvFields[domain].size()
+                );
+            }
+        }
+
+
+        // Combine bits. Note that can reuse field storage
 
         // Subset myself
         const labelList& mySubMap = subMap[Pstream::myProcNo()];
@@ -175,54 +279,56 @@ void Foam::mapDistribute::distribute
             subField[i] = field[mySubMap[i]];
         }
 
+        field.setSize(constructSize);
+
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
 
         forAll(map, i)
         {
-            newField[map[i]] = subField[i];
-        }
-
-        // Send sub field to neighbour
-        for (label domain = 0; domain < Pstream::nProcs(); domain++)
-        {
-            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 << subField;
-            }
+            field[map[i]] = subField[i];
         }
 
 
-        // Receive sub field from neighbour
+        // Wait for all to finish
+
+        OPstream::waitRequests();
+        IPstream::waitRequests();
+        
+        // Collect neighbour fields
+
         for (label domain = 0; domain < Pstream::nProcs(); domain++)
         {
             const labelList& map = constructMap[domain];
 
             if (domain != Pstream::myProcNo() && map.size() > 0)
             {
-                IPstream fromNbr(Pstream::nonBlocking, domain);
-                List<T> subField(fromNbr);
+                if (recvFields[domain].size() != map.size())
+                {
+                    FatalErrorIn
+                    (
+                        "template<class T>\n"
+                        "void mapDistribute::distribute\n"
+                        "(\n"
+                        "    const Pstream::commsTypes commsType,\n"
+                        "    const List<labelPair>& schedule,\n"
+                        "    const label constructSize,\n"
+                        "    const labelListList& subMap,\n"
+                        "    const labelListList& constructMap,\n"
+                        "    List<T>& field\n"
+                        ")\n"
+                    )   << "Expected from processor " << domain
+                        << " " << map.size() << " but received "
+                        << recvFields[domain].size() << " elements."
+                        << abort(FatalError);
+                }
 
                 forAll(map, i)
                 {
-                    newField[map[i]] = subField[i];
+                    field[map[i]] = recvFields[domain][i];
                 }
             }
         }
-        OPstream::waitRequests();
-        IPstream::waitRequests();
-
-        field.transfer(newField);
     }
     else
     {
diff --git a/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H b/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H
index 7312e08e55..48c0c566fc 100644
--- a/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H
+++ b/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H
@@ -221,6 +221,41 @@ public:
                 const bool applySeparation
             );
 
+        // PackedList versions
+
+            template <int nBits, class CombineOp>
+            static void syncFaceList
+            (
+                const polyMesh& mesh,
+                PackedList<nBits>& faceValues,
+                const CombineOp& cop
+            );
+
+            template <int nBits>
+            static void swapFaceList
+            (
+                const polyMesh& mesh,
+                PackedList<nBits>& faceValues
+            );
+
+            template <int nBits, class CombineOp>
+            static void syncPointList
+            (
+                const polyMesh& mesh,
+                PackedList<nBits>& pointValues,
+                const CombineOp& cop,
+                const unsigned int nullValue
+            );
+
+            template <int nBits, class CombineOp>
+            static void syncEdgeList
+            (
+                const polyMesh& mesh,
+                PackedList<nBits>& edgeValues,
+                const CombineOp& cop,
+                const unsigned int nullValue
+            );
+
         // Other
 
             //- Get per point whether is it master (of a coupled set of points)
diff --git a/src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C b/src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C
index 9715ba65b2..3ee897d97d 100644
--- a/src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C
@@ -940,17 +940,7 @@ void Foam::syncTools::syncPointList
                     fromNbr >> nbrPatchInfo;
                 }
                 // Null any value which is not on neighbouring processor
-                label nbrSize = nbrPatchInfo.size();
-                nbrPatchInfo.setSize(procPatch.nPoints());
-                for
-                (
-                    label i = nbrSize;
-                    i < procPatch.nPoints();
-                    i++
-                )
-                {
-                    nbrPatchInfo[i] = nullValue;
-                }
+                nbrPatchInfo.setSize(procPatch.nPoints(), nullValue);
 
                 if (!procPatch.parallel())
                 {
@@ -1225,18 +1215,7 @@ void Foam::syncTools::syncEdgeList
                     fromNeighb >> nbrPatchInfo;
                 }
                 // Null any value which is not on neighbouring processor
-                label nbrSize = nbrPatchInfo.size();
-                nbrPatchInfo.setSize(procPatch.nEdges());
-                for
-                (
-                    label i = nbrSize;
-                    i < procPatch.nEdges();
-                    i++
-                )
-                {
-                    nbrPatchInfo[i] = nullValue;
-                }
-
+                nbrPatchInfo.setSize(procPatch.nEdges(), nullValue);
 
                 if (!procPatch.parallel())
                 {
@@ -1591,4 +1570,469 @@ void Foam::syncTools::swapFaceList
 }
 
 
+template <int nBits, class CombineOp>
+void Foam::syncTools::syncFaceList
+(
+    const polyMesh& mesh,
+    PackedList<nBits>& faceValues,
+    const CombineOp& cop
+)
+{
+    if (faceValues.size() != mesh.nFaces())
+    {
+        FatalErrorIn
+        (
+            "syncTools<int nBits, class CombineOp>::syncFaceList"
+            "(const polyMesh&, PackedList<nBits>&, const CombineOp&)"
+        )   << "Number of values " << faceValues.size()
+            << " is not equal to the number of faces in the mesh "
+            << mesh.nFaces() << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    if (!hasCouples(patches))
+    {
+        return;
+    }
+
+    // Patch data (proc patches only).
+    List<List<unsigned int> > patchValues(patches.size());
+
+    if (Pstream::parRun())
+    {
+        // Send
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].size() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                patchValues[patchI].setSize(procPatch.size());
+                forAll(procPatch, i)
+                {
+                    patchValues[patchI][i] =
+                        faceValues.get(procPatch.start()+i);
+                }
+
+                OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
+                toNbr << patchValues[patchI];
+            }
+        }
+
+
+        // Receive and combine.
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].size() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                {
+                    IPstream fromNbr
+                    (
+                        Pstream::blocking,
+                        procPatch.neighbProcNo()
+                    );
+                    fromNbr >> patchValues[patchI];
+                }
+
+                // Combine (bitwise)
+                forAll(procPatch, i)
+                {
+                    unsigned int patchVal = patchValues[patchI][i];
+                    label meshFaceI = procPatch.start()+i;
+                    unsigned int faceVal = faceValues.get(meshFaceI);
+                    cop(faceVal, patchVal);
+                    faceValues.set(meshFaceI, faceVal);
+                }
+            }
+        }
+    }
+
+    // Do the cyclics.
+    forAll(patches, patchI)
+    {
+        if (isA<cyclicPolyPatch>(patches[patchI]))
+        {
+            const cyclicPolyPatch& cycPatch =
+                refCast<const cyclicPolyPatch>(patches[patchI]);
+
+            label half = cycPatch.size()/2;
+
+            for (label i = 0; i < half; i++)
+            {
+                label meshFace0 = cycPatch.start()+i;
+                unsigned int val0 = faceValues.get(meshFace0);
+                label meshFace1 = meshFace0 + half;
+                unsigned int val1 = faceValues.get(meshFace1);
+
+                unsigned int t = val0;
+                cop(t, val1);
+                faceValues.set(meshFace0, t);
+
+                cop(val1, val0);                
+                faceValues.set(meshFace1, val1);
+            }
+        }
+    }
+}
+
+
+template <int nBits>
+void Foam::syncTools::swapFaceList
+(
+    const polyMesh& mesh,
+    PackedList<nBits>& faceValues
+)
+{
+    syncFaceList(mesh, faceValues, eqOp<unsigned int>());
+}
+
+
+template <int nBits, class CombineOp>
+void Foam::syncTools::syncPointList
+(
+    const polyMesh& mesh,
+    PackedList<nBits>& pointValues,
+    const CombineOp& cop,
+    const unsigned int nullValue
+)
+{
+    if (pointValues.size() != mesh.nPoints())
+    {
+        FatalErrorIn
+        (
+            "syncTools<int nBits, class CombineOp>::syncPointList"
+            "(const polyMesh&, PackedList<nBits>&, const CombineOp&"
+            ", const unsigned int&)"
+        )   << "Number of values " << pointValues.size()
+            << " is not equal to the number of points in the mesh "
+            << mesh.nPoints() << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    if (!hasCouples(patches))
+    {
+        return;
+    }
+
+    // Patch data (proc patches only).
+    List<List<unsigned int> > patchValues(patches.size());
+
+    if (Pstream::parRun())
+    {
+        // Send
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nPoints() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                patchValues[patchI].setSize(procPatch.nPoints());
+                patchValues[patchI] = nullValue;
+
+                const labelList& meshPts = procPatch.meshPoints();
+                const labelList& nbrPts = procPatch.neighbPoints();
+
+                forAll(nbrPts, pointI)
+                {
+                    label nbrPointI = nbrPts[pointI];
+                    if (nbrPointI >= 0 && nbrPointI < procPatch.nPoints())
+                    {
+                        patchValues[patchI][nbrPointI] =
+                            pointValues.get(meshPts[pointI]);
+                    }
+                }
+
+                OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
+                toNbr << patchValues[patchI];
+            }
+        }
+
+
+        // Receive and combine.
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nPoints() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                {
+                    // We do not know the number of points on the other side
+                    // so cannot use Pstream::read.
+                    IPstream fromNbr
+                    (
+                        Pstream::blocking,
+                        procPatch.neighbProcNo()
+                    );
+                    fromNbr >> patchValues[patchI];
+                }
+
+                // Null any value which is not on neighbouring processor
+                patchValues[patchI].setSize(procPatch.nPoints(), nullValue);
+
+                const labelList& meshPts = procPatch.meshPoints();
+
+                forAll(meshPts, pointI)
+                {
+                    label meshPointI = meshPts[pointI];
+                    unsigned int pointVal = pointValues.get(meshPointI);
+                    cop(pointVal, patchValues[patchI][pointI]);
+                    pointValues.set(meshPointI, pointVal);
+                }
+            }
+        }
+    }
+
+    // Do the cyclics.
+    forAll(patches, patchI)
+    {
+        if (isA<cyclicPolyPatch>(patches[patchI]))
+        {
+            const cyclicPolyPatch& cycPatch =
+                refCast<const cyclicPolyPatch>(patches[patchI]);
+
+            const edgeList& coupledPoints = cycPatch.coupledPoints();
+            const labelList& meshPts = cycPatch.meshPoints();
+
+            forAll(coupledPoints, i)
+            {
+                const edge& e = coupledPoints[i];
+
+                label point0 = meshPts[e[0]];
+                label point1 = meshPts[e[1]];
+
+                unsigned int val0 = pointValues.get(point0);
+                unsigned int t = val0;
+                unsigned int val1 = pointValues.get(point1);
+
+                cop(t, val1);
+                pointValues.set(point0, t);
+                cop(val1, val0);
+                pointValues.set(point1, val1);
+            }
+        }
+    }
+
+    // Synchronize multiple shared points.
+    const globalMeshData& pd = mesh.globalData();
+
+    if (pd.nGlobalPoints() > 0)
+    {
+        // Values on shared points. Use unpacked storage for ease!
+        List<unsigned int> sharedPts(pd.nGlobalPoints(), nullValue);
+
+        forAll(pd.sharedPointLabels(), i)
+        {
+            label meshPointI = pd.sharedPointLabels()[i];
+            // Fill my entries in the shared points
+            sharedPts[pd.sharedPointAddr()[i]] = pointValues.get(meshPointI);
+        }
+
+        // Combine on master.
+        Pstream::listCombineGather(sharedPts, cop);
+        Pstream::listCombineScatter(sharedPts);
+
+        // Now we will all have the same information. Merge it back with
+        // my local information.
+        forAll(pd.sharedPointLabels(), i)
+        {
+            label meshPointI = pd.sharedPointLabels()[i];
+            pointValues.set(meshPointI, sharedPts[pd.sharedPointAddr()[i]]);
+        }
+    }
+}
+
+
+template <int nBits, class CombineOp>
+void Foam::syncTools::syncEdgeList
+(
+    const polyMesh& mesh,
+    PackedList<nBits>& edgeValues,
+    const CombineOp& cop,
+    const unsigned int nullValue
+)
+{
+    if (edgeValues.size() != mesh.nEdges())
+    {
+        FatalErrorIn
+        (
+            "syncTools<int nBits, class CombineOp>::syncEdgeList"
+            "(const polyMesh&, PackedList<nBits>&, const CombineOp&"
+            ", const unsigned int&)"
+        )   << "Number of values " << edgeValues.size()
+            << " is not equal to the number of edges in the mesh "
+            << mesh.nEdges() << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    if (!hasCouples(patches))
+    {
+        return;
+    }
+
+    // Patch data (proc patches only).
+    List<List<unsigned int> > patchValues(patches.size());
+
+    if (Pstream::parRun())
+    {
+        // Send
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nEdges() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                patchValues[patchI].setSize(procPatch.nEdges(), nullValue);
+
+                const labelList& meshEdges = procPatch.meshEdges();
+                const labelList& neighbEdges = procPatch.neighbEdges();
+
+                forAll(neighbEdges, edgeI)
+                {
+                    label nbrEdgeI = neighbEdges[edgeI];
+                    if (nbrEdgeI >= 0 && nbrEdgeI < procPatch.nEdges())
+                    {
+                        patchValues[patchI][nbrEdgeI] =
+                            edgeValues.get(meshEdges[edgeI]);
+                    }
+                }
+
+                OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
+                toNbr << patchValues[patchI];
+            }
+        }
+
+
+        // Receive and combine.
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nEdges() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                {
+                    IPstream fromNeighb
+                    (
+                        Pstream::blocking,
+                        procPatch.neighbProcNo()
+                    );
+                    fromNeighb >> patchValues[patchI];
+                }
+
+                patchValues[patchI].setSize(procPatch.nEdges(), nullValue);
+
+                const labelList& meshEdges = procPatch.meshEdges();
+
+                forAll(meshEdges, edgeI)
+                {
+                    unsigned int patchVal = patchValues[patchI][edgeI];
+                    label meshEdgeI = meshEdges[edgeI];
+                    unsigned int edgeVal = edgeValues.get(meshEdgeI);
+                    cop(edgeVal, patchVal);
+                    edgeValues.set(meshEdgeI, edgeVal);
+                }
+            }
+        }
+    }
+
+    // Do the cyclics.
+    forAll(patches, patchI)
+    {
+        if (isA<cyclicPolyPatch>(patches[patchI]))
+        {
+            const cyclicPolyPatch& cycPatch =
+                refCast<const cyclicPolyPatch>(patches[patchI]);
+
+            const edgeList& coupledEdges = cycPatch.coupledEdges();
+            const labelList& meshEdges = cycPatch.meshEdges();
+
+            forAll(coupledEdges, i)
+            {
+                const edge& e = coupledEdges[i];
+
+                label edge0 = meshEdges[e[0]];
+                label edge1 = meshEdges[e[1]];
+
+                unsigned int val0 = edgeValues.get(edge0);
+                unsigned int t = val0;
+                unsigned int val1 = edgeValues.get(edge1);
+
+                cop(t, val1);
+                edgeValues.set(edge0, t);
+                cop(val1, val0);
+                edgeValues.set(edge1, val1);
+            }
+        }
+    }
+
+    // Synchronize multiple shared edges.
+    const globalMeshData& pd = mesh.globalData();
+
+    if (pd.nGlobalEdges() > 0)
+    {
+        // Values on shared edges. Use unpacked storage for ease!
+        List<unsigned int> sharedPts(pd.nGlobalEdges(), nullValue);
+
+        forAll(pd.sharedEdgeLabels(), i)
+        {
+            label meshEdgeI = pd.sharedEdgeLabels()[i];
+            // Fill my entries in the shared edges
+            sharedPts[pd.sharedEdgeAddr()[i]] = edgeValues.get(meshEdgeI);
+        }
+
+        // Combine on master.
+        Pstream::listCombineGather(sharedPts, cop);
+        Pstream::listCombineScatter(sharedPts);
+
+        // Now we will all have the same information. Merge it back with
+        // my local information.
+        forAll(pd.sharedEdgeLabels(), i)
+        {
+            label meshEdgeI = pd.sharedEdgeLabels()[i];
+            edgeValues.set(meshEdgeI, sharedPts[pd.sharedEdgeAddr()[i]]);
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatchExtra/PrimitivePatchExtraSearch.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatchExtra/PrimitivePatchExtraSearch.C
index 453cdd63c9..ffad4f005f 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatchExtra/PrimitivePatchExtraSearch.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatchExtra/PrimitivePatchExtraSearch.C
@@ -102,8 +102,9 @@ void Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::markZone
             break;
         }
 
-        changedFaces.transfer(newChangedFaces.shrink());
-        newChangedFaces.clear();
+        // New dynamiclist: can leave dynamicList unshrunk
+        //changedFaces.transfer(newChangedFaces.shrink());
+        changedFaces.transfer(newChangedFaces);
     }
 }
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C
index b60b366079..b64eda7266 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C
@@ -441,11 +441,10 @@ void primitiveMesh::calcEdges(const bool doFaceEdges) const
         forAll(pe, pointI)
         {
             DynamicList<label>& pEdges = pe[pointI];
-            inplaceRenumber(oldToNew, pEdges);
             pEdges.shrink();
+            inplaceRenumber(oldToNew, pEdges);
             pointEdges[pointI].transfer(pEdges);
             Foam::sort(pointEdges[pointI]);
-            pEdges.clear();
         }
 
         // faceEdges
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitivePatch/primitivePatch.H b/src/OpenFOAM/meshes/primitiveMesh/primitivePatch/primitivePatch.H
index af035e9a03..a51ed319ad 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitivePatch/primitivePatch.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitivePatch/primitivePatch.H
@@ -26,7 +26,7 @@ Typedef
     Foam::primitivePatch
 
 Description
-    Foam::primitivePatch
+    Addressing for a faceList slice.
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/src/OpenFOAM/meshes/primitiveShapes/line/linePointRef.H b/src/OpenFOAM/meshes/primitiveShapes/line/linePointRef.H
index b0e84c7aa9..bc2a25d381 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/line/linePointRef.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/line/linePointRef.H
@@ -26,6 +26,7 @@ Typedef
     Foam::linePointRef
 
 Description
+    Line using referred points
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/src/OpenFOAM/primitives/contiguous/contiguous.H b/src/OpenFOAM/primitives/contiguous/contiguous.H
index d634bb0b60..59c206dd22 100644
--- a/src/OpenFOAM/primitives/contiguous/contiguous.H
+++ b/src/OpenFOAM/primitives/contiguous/contiguous.H
@@ -38,6 +38,8 @@ Description
 #ifndef contiguous_H
 #define contiguous_H
 
+#include "label.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -45,54 +47,117 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Forward declaration of friend functions and operators
+template<class T, label Size> class FixedList;
+template<class T> class Pair;
+
+
 // Assume the data associated with type T is not contiguous
 template<class T>
 inline bool contiguous()                        {return false;}
 
 
-// Specify data associated with primitive types is contiguous
+
+// Specify data associated with primitive types (and simple fixed size
+// containers - only size 2 defined here) is contiguous
 
 template<>
-inline bool contiguous<bool>()                  {return true;}
+inline bool contiguous<bool>()                              {return true;}
+template<>
+inline bool contiguous<FixedList<bool, 2> >()               {return true;}
+template<>
+inline bool contiguous<Pair<bool> >()                       {return true;}
 
 template<>
-inline bool contiguous<char>()                  {return true;}
+inline bool contiguous<char>()                              {return true;}
+template<>
+inline bool contiguous<FixedList<char, 2> >()               {return true;}
+template<>
+inline bool contiguous<Pair<char> >()                       {return true;}
 
 template<>
-inline bool contiguous<unsigned char>()         {return true;}
+inline bool contiguous<unsigned char>()                     {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned char, 2> >()      {return true;}
+template<>
+inline bool contiguous<Pair<unsigned char> >()              {return true;}
 
 template<>
-inline bool contiguous<short>()                 {return true;}
+inline bool contiguous<short>()                             {return true;}
+template<>
+inline bool contiguous<FixedList<short, 2> >()              {return true;}
+template<>
+inline bool contiguous<Pair<short> >()                      {return true;}
 
 template<>
-inline bool contiguous<unsigned short>()        {return true;}
+inline bool contiguous<unsigned short>()                    {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned short, 2> >()     {return true;}
+template<>
+inline bool contiguous<Pair<unsigned short> >()             {return true;}
 
 template<>
-inline bool contiguous<int>()                   {return true;}
+inline bool contiguous<int>()                               {return true;}
+template<>
+inline bool contiguous<FixedList<int, 2> >()                {return true;}
+template<>
+inline bool contiguous<Pair<int> >()                        {return true;}
 
 template<>
-inline bool contiguous<unsigned int>()          {return true;}
+inline bool contiguous<unsigned int>()                      {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned int, 2> >()       {return true;}
+template<>
+inline bool contiguous<Pair<unsigned int> >()               {return true;}
 
 template<>
-inline bool contiguous<long>()                  {return true;}
+inline bool contiguous<long>()                              {return true;}
+template<>
+inline bool contiguous<FixedList<long, 2> >()               {return true;}
+template<>
+inline bool contiguous<Pair<long> >()                       {return true;}
 
 template<>
-inline bool contiguous<unsigned long>()         {return true;}
+inline bool contiguous<unsigned long>()                     {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned long, 2> >()      {return true;}
+template<>
+inline bool contiguous<Pair<unsigned long> >()              {return true;}
 
 template<>
-inline bool contiguous<long long>()             {return true;}
+inline bool contiguous<long long>()                         {return true;}
+template<>
+inline bool contiguous<FixedList<long long, 2> >()          {return true;}
+template<>
+inline bool contiguous<Pair<long long> >()                  {return true;}
 
 template<>
-inline bool contiguous<unsigned long long>()    {return true;}
+inline bool contiguous<unsigned long long>()                {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned long long, 2> >() {return true;}
+template<>
+inline bool contiguous<Pair<unsigned long long> >()         {return true;}
 
 template<>
-inline bool contiguous<float>()                 {return true;}
+inline bool contiguous<float>()                             {return true;}
+template<>
+inline bool contiguous<FixedList<float, 2> >()              {return true;}
+template<>
+inline bool contiguous<Pair<float> >()                       {return true;}
 
 template<>
-inline bool contiguous<double>()                {return true;}
+inline bool contiguous<double>()                            {return true;}
+template<>
+inline bool contiguous<FixedList<double, 2> >()             {return true;}
+template<>
+inline bool contiguous<Pair<double> >()                     {return true;}
 
 template<>
-inline bool contiguous<long double>()           {return true;}
+inline bool contiguous<long double>()                       {return true;}
+template<>
+inline bool contiguous<FixedList<long double, 2> >()        {return true;}
+template<>
+inline bool contiguous<Pair<long double> >()                {return true;}
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index ca0c0cbddc..57a4cb5b64 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -25,6 +25,7 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "autoRefineDriver.H"
+#include "meshRefinement.H"
 #include "fvMesh.H"
 #include "Time.H"
 #include "boundBox.H"
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
index f1747a29be..ce788e4b1e 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef autoRefineDriver_H
 #define autoRefineDriver_H
 
-#include "meshRefinement.H"
+#include "treeBoundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +45,9 @@ namespace Foam
 // Forward declaration of classes
 class featureEdgeMesh;
 class refinementParameters;
+class meshRefinement;
+class decompositionMethod;
+class fvMeshDistribute;
 
 /*---------------------------------------------------------------------------*\
                            Class autoRefineDriver Declaration
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 7f5bc29a79..626fb2fe25 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -51,6 +51,9 @@ License
 #include "meshTools.H"
 #include "OFstream.H"
 #include "geomDecomp.H"
+#include "Random.H"
+#include "searchableSurfaces.H"
+#include "treeBoundBox.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -122,13 +125,36 @@ void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
 {
     const pointField& cellCentres = mesh_.cellCentres();
 
-    label nTotEdges = returnReduce(surfaceIndex_.size(), sumOp<label>());
-    label nChangedFaces = returnReduce(changedFaces.size(), sumOp<label>());
+    // Stats on edges to test. Count proc faces only once.
+    PackedList<1> isMasterFace(syncTools::getMasterFaces(mesh_));
+
+    {
+        label nMasterFaces = 0;
+        forAll(isMasterFace, faceI)
+        {
+            if (isMasterFace.get(faceI) == 1)
+            {
+                nMasterFaces++;
+            }
+        }
+        reduce(nMasterFaces, sumOp<label>());
+
+        label nChangedFaces = 0;
+        forAll(changedFaces, i)
+        {
+            if (isMasterFace.get(changedFaces[i]) == 1)
+            {
+                nChangedFaces++;
+            }
+        }
+        reduce(nChangedFaces, sumOp<label>());
+
+        Info<< "Edge intersection testing:" << nl
+            << "    Number of edges             : " << nMasterFaces << nl
+            << "    Number of edges to retest   : " << nChangedFaces
+            << endl;
+    }
 
-    Info<< "Edge intersection testing:" << nl
-        << "    Number of edges             : " << nTotEdges << nl
-        << "    Number of edges to retest   : " << nChangedFaces
-        << endl;
 
     // Get boundary face centre and level. Coupled aware.
     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
@@ -838,11 +864,14 @@ Foam::meshRefinement::meshRefinement
 
 Foam::label Foam::meshRefinement::countHits() const
 {
+    // Stats on edges to test. Count proc faces only once.
+    PackedList<1> isMasterFace(syncTools::getMasterFaces(mesh_));
+
     label nHits = 0;
 
     forAll(surfaceIndex_, faceI)
     {
-        if (surfaceIndex_[faceI] >= 0)
+        if (surfaceIndex_[faceI] >= 0 && isMasterFace.get(faceI) == 1)
         {
             nHits++;
         }
@@ -996,11 +1025,6 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 
     if (Pstream::parRun())
     {
-        //Info<< nl
-        //    << "Doing final balancing" << nl
-        //    << "---------------------" << nl
-        //    << endl;
-        //
         //if (debug_)
         //{
         //    const_cast<Time&>(mesh_.time())++;
@@ -1014,17 +1038,12 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
             // Faces where owner and neighbour are not 'connected' so can
             // go to different processors.
             boolList blockedFace(mesh_.nFaces(), true);
+            label nUnblocked = 0;
             // Pairs of baffles
             List<labelPair> couples;
 
             if (keepZoneFaces)
             {
-                label nNamed = surfaces().getNamedSurfaces().size();
-
-                Info<< "Found " << nNamed << " surfaces with faceZones."
-                    << " Applying special decomposition to keep those together."
-                    << endl;
-
                 // Determine decomposition to keep/move surface zones
                 // on one processor. The reason is that snapping will make these
                 // into baffles, move and convert them back so if they were
@@ -1039,8 +1058,6 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                 // Get faces whose owner and neighbour should stay together,
                 // i.e. they are not 'blocked'.
 
-                label nZoned = 0;
-
                 forAll(fzNames, surfI)
                 {
                     if (fzNames[surfI].size() > 0)
@@ -1055,14 +1072,18 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                             if (blockedFace[fZone[i]])
                             {
                                 blockedFace[fZone[i]] = false;
-                                nZoned++;
+                                nUnblocked++;
                             }
                         }
                     }
                 }
-                Info<< "Found " << returnReduce(nZoned, sumOp<label>())
-                    << " zoned faces to keep together."
-                    << endl;
+            }
+            reduce(nUnblocked, sumOp<label>());
+
+            if (keepZoneFaces)
+            {
+                Info<< "Found " << nUnblocked
+                    << " zoned faces to keep together." << endl;
             }
 
             if (keepBaffles)
@@ -1073,29 +1094,43 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                     identity(mesh_.nFaces()-mesh_.nInternalFaces())
                    +mesh_.nInternalFaces()
                 );
+            }
+            label nCouples = returnReduce(couples.size(), sumOp<label>());
 
-                Info<< "Found " << returnReduce(couples.size(), sumOp<label>())
-                    << " baffles to keep together."
+            if (keepBaffles)
+            {
+                Info<< "Found " << nCouples << " baffles to keep together."
                     << endl;
             }
 
-            distribution = decomposeCombineRegions
-            (
-                blockedFace,
-                couples,
-                decomposer
-            );
-
-            labelList nProcCells(distributor.countCells(distribution));
-            Pstream::listCombineGather(nProcCells, plusEqOp<label>());
-            Pstream::listCombineScatter(nProcCells);
-
-            Info<< "Calculated decomposition:" << endl;
-            forAll(nProcCells, procI)
+            if (nUnblocked > 0 || nCouples > 0)
             {
-                Info<< "    " << procI << '\t' << nProcCells[procI] << endl;
+                Info<< "Applying special decomposition to keep baffles"
+                    << " and zoned faces together." << endl;
+
+                distribution = decomposeCombineRegions
+                (
+                    blockedFace,
+                    couples,
+                    decomposer
+                );
+
+                labelList nProcCells(distributor.countCells(distribution));
+                Pstream::listCombineGather(nProcCells, plusEqOp<label>());
+                Pstream::listCombineScatter(nProcCells);
+
+                Info<< "Calculated decomposition:" << endl;
+                forAll(nProcCells, procI)
+                {
+                    Info<< "    " << procI << '\t' << nProcCells[procI] << endl;
+                }
+                Info<< endl;
+            }
+            else
+            {
+                // Normal decomposition
+                distribution = decomposer.decompose(mesh_.cellCentres());
             }
-            Info<< endl;
         }
         else
         {
@@ -1715,6 +1750,36 @@ void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
     {
         map.distributeFaceData(userFaceData_[i].second());
     }
+
+    // Redistribute surface and any fields on it.
+    {
+        Random rndGen(653213);
+
+        // Get local mesh bounding box. Single box for now.
+        List<treeBoundBox> meshBb(1);
+        treeBoundBox& bb = meshBb[0];
+        bb = boundBox(mesh_.points(), false);
+        bb = bb.extend(rndGen, 1E-4);
+
+        // Distribute all geometry (so refinementSurfaces and shellSurfaces)
+        searchableSurfaces& geometry =
+            const_cast<searchableSurfaces&>(surfaces_.geometry());
+
+        forAll(geometry, i)
+        {
+            autoPtr<mapDistribute> faceMap;
+            autoPtr<mapDistribute> pointMap;
+            geometry[i].distribute
+            (
+                meshBb,
+                false,          // do not keep outside triangles
+                faceMap,
+                pointMap
+            );
+            faceMap.clear();
+            pointMap.clear();
+        }
+    }
 }
 
 
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index e6bd37f670..58c74d3fb4 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -386,7 +386,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
         FatalErrorIn
         (
             "meshRefinement::createBaffles"
-            "(const label, const labelList&, const labelList&)"
+            "(const labelList&, const labelList&)"
         )   << "Illegal size :"
             << " ownPatch:" << ownPatch.size()
             << " neiPatch:" << neiPatch.size()
@@ -412,7 +412,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
                 FatalErrorIn
                 (
                     "meshRefinement::createBaffles"
-                    "(const label, const labelList&, const labelList&)"
+                    "(const labelList&, const labelList&)"
                 )   << "Non synchronised at face:" << faceI
                     << " on patch:" << mesh_.boundaryMesh().whichPatch(faceI)
                     << " fc:" << mesh_.faceCentres()[faceI] << endl
@@ -461,7 +461,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
 
     //- Redo the intersections on the newly create baffle faces. Note that
     //  this changes also the cell centre positions.
-    labelHashSet baffledFacesSet(2*nBaffles);
+    faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
 
     const labelList& reverseFaceMap = map().reverseFaceMap();
     const labelList& faceMap = map().faceMap();
@@ -496,6 +496,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
             }
         }
     }
+    baffledFacesSet.sync(mesh_);
 
     updateMesh(map, baffledFacesSet.toc());
 
@@ -1886,15 +1887,15 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
         labelList facePatch
         (
-            //markFacesOnProblemCells
-            //(
-            //    globalToPatch
-            //)
-            markFacesOnProblemCellsGeometric
+            markFacesOnProblemCells
             (
-                motionDict,
                 globalToPatch
             )
+            //markFacesOnProblemCellsGeometric
+            //(
+            //    motionDict,
+            //    globalToPatch
+            //)
         );
         Info<< "Analyzed problem cells in = "
             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index e9eb2cb445..b457f0f5c7 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -1091,8 +1091,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
         //    << "    local allowable refinement:" << nAllowRefine << nl
         //    << endl;
 
-        //- Disable refinement shortcut
-        label nAllowRefine = labelMax;
+        //- Disable refinement shortcut. nAllowRefine is per processor limit.
+        label nAllowRefine = labelMax / Pstream::nProcs();
 
         // Marked for refinement (>= 0) or not (-1). Actual value is the
         // index of the surface it intersects.
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index 1654e3cd7f..9629b53bdc 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -25,19 +25,13 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "refinementSurfaces.H"
-#include "orientedSurface.H"
 #include "Time.H"
 #include "searchableSurfaces.H"
 #include "shellSurfaces.H"
 #include "triSurfaceMesh.H"
 #include "labelPair.H"
 #include "searchableSurfacesQueries.H"
-
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+#include "UPtrList.H"
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -438,31 +432,6 @@ Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const
 }
 
 
-// Orient surface (if they're closed) before any searching is done.
-void Foam::refinementSurfaces::orientSurface
-(
-    const point& keepPoint,
-    triSurfaceMesh& s
-)
-{
-    // Flip surface so keepPoint is outside.
-    bool anyFlipped = orientedSurface::orient(s, keepPoint, true);
-
-    if (anyFlipped)
-    {
-        // orientedSurface will have done a clearOut of the surface.
-        // we could do a clearout of the triSurfaceMeshes::trees()
-        // but these aren't affected by orientation (except for cached
-        // sideness which should not be set at this point. !!Should
-        // check!)
-
-        Info<< "orientSurfaces : Flipped orientation of surface "
-            << s.searchableSurface::name()
-            << " so point " << keepPoint << " is outside." << endl;
-    }
-}
-
-
 // Count number of triangles per surface region
 Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s)
 {
@@ -485,7 +454,7 @@ void Foam::refinementSurfaces::setMinLevelFields
     const shellSurfaces& shells
 )
 {
-    minLevelFields_.setSize(surfaces_.size());
+    //minLevelFields_.setSize(surfaces_.size());
 
     forAll(surfaces_, surfI)
     {
@@ -495,26 +464,24 @@ void Foam::refinementSurfaces::setMinLevelFields
         {
             const triSurfaceMesh& triMesh = refCast<const triSurfaceMesh>(geom);
 
-            minLevelFields_.set
+            autoPtr<triSurfaceLabelField> minLevelFieldPtr
             (
-                surfI,
                 new triSurfaceLabelField
                 (
                     IOobject
                     (
-                        triMesh.searchableSurface::name(),
+                        "minLevel",
                         triMesh.objectRegistry::time().constant(),// directory
                         "triSurface",               // instance
                         triMesh,
                         IOobject::NO_READ,
-                        IOobject::AUTO_WRITE,
-                        false
+                        IOobject::AUTO_WRITE
                     ),
                     triMesh,
                     dimless
                 )
             );
-            triSurfaceLabelField& minLevelField = minLevelFields_[surfI];
+            triSurfaceLabelField& minLevelField = minLevelFieldPtr();
 
             const triSurface& s = static_cast<const triSurface&>(triMesh);
 
@@ -542,6 +509,9 @@ void Foam::refinementSurfaces::setMinLevelFields
                     shellLevel[triI]
                 );
             }
+
+            // Store field on triMesh
+            minLevelFieldPtr.ptr()->store();
         }
     }
 }
@@ -569,17 +539,38 @@ void Foam::refinementSurfaces::findHigherIntersection
         return;
     }
 
+
+    // Precalculate per surface whether it has a minlevelfield
+    UPtrList<triSurfaceLabelField> minLevelFields(surfaces_.size());
+    forAll(surfaces_, surfI)
+    {
+        const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+
+        if (isA<triSurfaceMesh>(geom))
+        {
+            const triSurfaceMesh& triMesh = refCast<const triSurfaceMesh>(geom);
+            minLevelFields.set
+            (
+                surfI,
+               &const_cast<triSurfaceLabelField&>
+                (
+                    triMesh.lookupObject<triSurfaceLabelField>("minLevel")
+                )
+            );
+        }
+    }
+
     // Work arrays
     labelList hitMap(identity(start.size()));
     pointField p0(start);
     pointField p1(end);
     List<pointIndexHit> hitInfo(start.size());
 
-
     forAll(surfaces_, surfI)
     {
         allGeometry_[surfaces_[surfI]].findLineAny(p0, p1, hitInfo);
 
+
         // Copy all hits into arguments, continue with misses
         label newI = 0;
         forAll(hitInfo, hitI)
@@ -592,12 +583,12 @@ void Foam::refinementSurfaces::findHigherIntersection
                 // Check if minLevelField for this surface.
                 if
                 (
-                    minLevelFields_.set(surfI)
-                 && minLevelFields_[surfI].size() > 0
+                    minLevelFields.set(surfI)
+                 && minLevelFields[surfI].size() > 0
                 )
                 {
                     minLocalLevel =
-                        minLevelFields_[surfI][hitInfo[hitI].index()];
+                        minLevelFields[surfI][hitInfo[hitI].index()];
                 }
                 else
                 {
@@ -668,30 +659,57 @@ void Foam::refinementSurfaces::findAllHigherIntersections
     {
         allGeometry_[surfaces_[surfI]].findLineAll(start, end, hitInfo);
 
+        // Repack hits for surface into flat list
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // To avoid overhead of calling getRegion for every point
+
+        label n = 0;
+        forAll(hitInfo, pointI)
+        {
+            n += hitInfo[pointI].size();
+        }
+
+        List<pointIndexHit> surfInfo(n);
+        labelList pointMap(n);
+        n = 0;
+
         forAll(hitInfo, pointI)
         {
             const List<pointIndexHit>& pHits = hitInfo[pointI];
-            allGeometry_[surfaces_[surfI]].getRegion(pHits, pRegions);
-            allGeometry_[surfaces_[surfI]].getNormal(pHits, pNormals);
 
-            // Extract those hits that are on higher levelled surfaces.
-            // Note: should move extraction of region, normal outside of loop
-            // below for if getRegion/getNormal have high overhead.
-
-            forAll(pHits, pHitI)
+            forAll(pHits, i)
             {
-                label region = globalRegion(surfI, pRegions[pHitI]);
+                surfInfo[n] = pHits[i];
+                pointMap[n] = pointI;
+                n++;
+            }
+        }
 
-                if (maxLevel_[region] > currentLevel[pointI])
-                {
-                    // Append to pointI info
-                    label sz = surfaceNormal[pointI].size();
-                    surfaceNormal[pointI].setSize(sz+1);
-                    surfaceNormal[pointI][sz] = pNormals[pHitI];
+        labelList surfRegion(n);
+        vectorField surfNormal(n);
+        allGeometry_[surfaces_[surfI]].getRegion(surfInfo, surfRegion);
+        allGeometry_[surfaces_[surfI]].getNormal(surfInfo, surfNormal);
 
-                    surfaceLevel[pointI].setSize(sz+1);
-                    surfaceLevel[pointI][sz] = maxLevel_[region];
-                }
+        surfInfo.clear();
+
+
+        // Extract back into pointwise
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        forAll(surfRegion, i)
+        {
+            label region = globalRegion(surfI, surfRegion[i]);
+            label pointI = pointMap[i];
+
+            if (maxLevel_[region] > currentLevel[pointI])
+            {
+                // Append to pointI info
+                label sz = surfaceNormal[pointI].size();
+                surfaceNormal[pointI].setSize(sz+1);
+                surfaceNormal[pointI][sz] = surfNormal[i];
+
+                surfaceLevel[pointI].setSize(sz+1);
+                surfaceLevel[pointI][sz] = maxLevel_[region];
             }
         }
     }
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index ee5dd219a5..a20df19f29 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -90,9 +90,6 @@ class refinementSurfaces
         //- From global region number to refinement level
         labelList maxLevel_;
 
-        //- Per surface refinement level adapted for shells.
-        PtrList<triSurfaceLabelField> minLevelFields_;
-
 
     // Private Member Functions
 
@@ -213,13 +210,6 @@ public:
                 const shellSurfaces& shells
             );
 
-            //- Helper: orient (closed only) surfaces so keepPoint is outside.
-            static void orientSurface
-            (
-                const point& keepPoint,
-                triSurfaceMesh& surface
-            );
-
             //- Helper: count number of triangles per region
             static labelList countRegions(const triSurface&);
 
diff --git a/src/dynamicMesh/meshCut/cellLooper/hexCellLooper.C b/src/dynamicMesh/meshCut/cellLooper/hexCellLooper.C
index dda1f3510d..62a90ae797 100644
--- a/src/dynamicMesh/meshCut/cellLooper/hexCellLooper.C
+++ b/src/dynamicMesh/meshCut/cellLooper/hexCellLooper.C
@@ -29,7 +29,6 @@ License
 #include "polyMesh.H"
 #include "cellModeller.H"
 #include "mathematicalConstants.H"
-#include "DynamicList.H"
 #include "plane.H"
 #include "ListOps.H"
 #include "meshTools.H"
diff --git a/src/dynamicMesh/meshCut/meshModifiers/meshCutAndRemove/meshCutAndRemove.H b/src/dynamicMesh/meshCut/meshModifiers/meshCutAndRemove/meshCutAndRemove.H
index ad8cdb46ed..4cc2e550b1 100644
--- a/src/dynamicMesh/meshCut/meshModifiers/meshCutAndRemove/meshCutAndRemove.H
+++ b/src/dynamicMesh/meshCut/meshModifiers/meshCutAndRemove/meshCutAndRemove.H
@@ -40,7 +40,6 @@ SourceFiles
 #include "boolList.H"
 #include "labelList.H"
 #include "typeInfo.H"
-#include "DynamicList.H"
 #include "Map.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.H b/src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.H
index a4a4b04bbc..4a33ea2ccc 100644
--- a/src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.H
+++ b/src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.H
@@ -113,10 +113,8 @@ SourceFiles
 #define meshCutter_H
 
 #include "edgeVertex.H"
-#include "boolList.H"
 #include "labelList.H"
 #include "typeInfo.H"
-#include "DynamicList.H"
 #include "Map.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/dynamicMesh/meshCut/meshModifiers/refinementIterator/refinementIterator.C b/src/dynamicMesh/meshCut/meshModifiers/refinementIterator/refinementIterator.C
index b7d7655dfb..f4dd69823c 100644
--- a/src/dynamicMesh/meshCut/meshModifiers/refinementIterator/refinementIterator.C
+++ b/src/dynamicMesh/meshCut/meshModifiers/refinementIterator/refinementIterator.C
@@ -26,12 +26,10 @@ License
 
 #include "refinementIterator.H"
 #include "polyMesh.H"
-#include "polyTopoChanger.H"
 #include "Time.H"
 #include "refineCell.H"
 #include "undoableMeshCutter.H"
 #include "polyTopoChange.H"
-#include "DynamicList.H"
 #include "mapPolyMesh.H"
 #include "cellCuts.H"
 #include "OFstream.H"
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
index da3d062fc1..ef711f432f 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
@@ -595,6 +595,7 @@ void Foam::edgeCollapser::updateMesh(const mapPolyMesh& map)
 {
     pointRegion_.setSize(mesh_.nPoints());
     pointRegion_ = -1;
+    // Reset count, do not remove underlying storage
     pointRegionMaster_.clear();
     freeRegions_.clear();
 }
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/faceCollapser.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/faceCollapser.H
index d31a4054dd..bd80ae1e77 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/faceCollapser.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/faceCollapser.H
@@ -54,7 +54,6 @@ SourceFiles
 #define faceCollapser_H
 
 #include "labelList.H"
-#include "DynamicList.H"
 #include "point.H"
 #include "Map.H"
 #include "labelHashSet.H"
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
index 0c567ee312..2a38afad1a 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
@@ -3021,7 +3021,9 @@ Foam::labelListList Foam::hexRef8::setRefinement
             << " Checking initial mesh just to make sure" << endl;
 
         checkMesh();
-        checkRefinementLevels(-1, labelList(0));
+        // Cannot call checkRefinementlevels since hanging points might
+        // get triggered by the mesher after subsetting.
+        //checkRefinementLevels(-1, labelList(0));
     }
 
     // Clear any saved point/cell data.
@@ -4813,7 +4815,7 @@ void Foam::hexRef8::checkRefinementLevels
             false                   // no separation
         );
 
-        OFstream str(mesh_.time().path()/"hangingPoints.obj");
+        //OFstream str(mesh_.time().path()/"hangingPoints.obj");
 
         label nHanging = 0;
 
@@ -4826,8 +4828,7 @@ void Foam::hexRef8::checkRefinementLevels
                 Pout<< "Hanging boundary point " << pointI
                     << " at " << mesh_.points()[pointI]
                     << endl;
-
-                meshTools::writeOBJ(str, mesh_.points()[pointI]);
+                //meshTools::writeOBJ(str, mesh_.points()[pointI]);
             }
         }
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChangeTemplates.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChangeTemplates.C
index 682161c0cc..f911802a44 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChangeTemplates.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChangeTemplates.C
@@ -96,13 +96,5 @@ void Foam::polyTopoChange::renumberKey
     elems.transfer(newElems);
 }
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
 
 // ************************************************************************* //
diff --git a/src/edgeMesh/edgeMesh.C b/src/edgeMesh/edgeMesh.C
index 27cb417b3e..ffc32de954 100644
--- a/src/edgeMesh/edgeMesh.C
+++ b/src/edgeMesh/edgeMesh.C
@@ -156,7 +156,7 @@ label edgeMesh::regions(labelList& edgeRegion) const
                 }
             }
 
-            edgesToVisit.transfer(newEdgesToVisit.shrink());
+            edgesToVisit.transfer(newEdgesToVisit);
         }
 
         currentRegion++;
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C
index f591376bf5..10299f6175 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C
@@ -1027,8 +1027,7 @@ CrankNicholsonDdtScheme<Type>::fvcDdtPhiCorr
             {
                 dUdt0 = rDtCoef0_(dUdt0)*
                 (
-                    rho.oldTime()*U.oldTime()
-                  - rho.oldTime().oldTime()*U.oldTime().oldTime()
+                    U.oldTime() - U.oldTime().oldTime()
                 ) - offCentre_(dUdt0());
             }
 
@@ -1054,11 +1053,8 @@ CrankNicholsonDdtScheme<Type>::fvcDdtPhiCorr
                       - (
                             fvc::interpolate
                             (
-                                rA*
-                                (
-                                    rDtCoef*rho.oldTime()*U.oldTime()
-                                    + offCentre_(dUdt0())
-                                )
+                                rA*rho.oldTime()
+                               *(rDtCoef*U.oldTime() + offCentre_(dUdt0()))
                             ) & mesh().Sf()
                         )
                     )
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 90f9de250e..64437bdda8 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -54,12 +54,15 @@ indexedOctree/treeDataPoint.C
 indexedOctree/treeDataTriSurface.C
 
 searchableSurface = searchableSurface
+$(searchableSurface)/distributedTriSurfaceMesh.C
 $(searchableSurface)/searchableBox.C
+$(searchableSurface)/searchablePlate.C
 $(searchableSurface)/searchableSphere.C
 $(searchableSurface)/searchableSurface.C
 $(searchableSurface)/searchableSurfaces.C
 $(searchableSurface)/searchableSurfacesQueries.C
 $(searchableSurface)/triSurfaceMesh.C
+$(searchableSurface)/searchableSurfaceWithGaps.C
 
 topoSets = sets/topoSets
 $(topoSets)/cellSet.C
diff --git a/src/meshTools/cellFeatures/cellFeatures.C b/src/meshTools/cellFeatures/cellFeatures.C
index 247aa0cfd4..498e154840 100644
--- a/src/meshTools/cellFeatures/cellFeatures.C
+++ b/src/meshTools/cellFeatures/cellFeatures.C
@@ -371,7 +371,7 @@ void Foam::cellFeatures::calcSuperFaces() const
                 {
                     superFace.shrink();
 
-                    faces[superFaceI] = face(superFace);
+                    faces[superFaceI].transfer(superFace);
                 }
             }
         }
diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C
index 8657f5c242..68b6026f53 100644
--- a/src/meshTools/meshSearch/meshSearch.C
+++ b/src/meshTools/meshSearch/meshSearch.C
@@ -921,14 +921,7 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
 
     hits.shrink();
 
-    // Copy into straight list
-    List<pointIndexHit> allHits(hits.size());
-
-    forAll(hits, hitI)
-    {
-        allHits[hitI] = hits[hitI];
-    }
-    return allHits;
+    return hits;
 }
 
 
diff --git a/src/meshTools/octree/pointIndexHit.H b/src/meshTools/octree/pointIndexHit.H
index b43f12d551..4a468b60e1 100644
--- a/src/meshTools/octree/pointIndexHit.H
+++ b/src/meshTools/octree/pointIndexHit.H
@@ -45,7 +45,7 @@ typedef PointIndexHit<point> pointIndexHit;
 
 //- Specify data associated with pointIndexHit type is contiguous
 template<>
-inline bool contiguous<pointIndexHit>() {return true;}
+inline bool contiguous<pointIndexHit>() {return contiguous<point>();}
 
 }
 
diff --git a/src/meshTools/octree/treeBoundBox.C b/src/meshTools/octree/treeBoundBox.C
index 26530b5a3a..1c67e72a7e 100644
--- a/src/meshTools/octree/treeBoundBox.C
+++ b/src/meshTools/octree/treeBoundBox.C
@@ -679,7 +679,11 @@ Foam::label Foam::treeBoundBox::distanceCmp
 
 bool Foam::operator==(const treeBoundBox& a, const treeBoundBox& b)
 {
-    return (a.min() == b.min()) && (a.max() == b.max());
+    return operator==
+    (
+        static_cast<const boundBox&>(a),
+        static_cast<const boundBox&>(b)
+    );
 }
 
 
@@ -689,12 +693,17 @@ bool Foam::operator!=(const treeBoundBox& a, const treeBoundBox& b)
 }
 
 
-// * * * * * * * * * * * * * * * IOstream Operator  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const treeBoundBox& bb)
+{
+    return os << static_cast<const boundBox&>(bb);
+}
+
 
 Foam::Istream& Foam::operator>>(Istream& is, treeBoundBox& bb)
 {
-    is >> bb.min() >> bb.max();
-    return is;
+    return is >> static_cast<boundBox&>(bb);
 }
 
 
diff --git a/src/meshTools/octree/treeBoundBox.H b/src/meshTools/octree/treeBoundBox.H
index 1805d88b76..926b2f3336 100644
--- a/src/meshTools/octree/treeBoundBox.H
+++ b/src/meshTools/octree/treeBoundBox.H
@@ -322,13 +322,19 @@ public:
         friend bool operator==(const treeBoundBox&, const treeBoundBox&);
         friend bool operator!=(const treeBoundBox&, const treeBoundBox&);
 
-
     // IOstream operator
 
-        friend Istream& operator>>(Istream&, treeBoundBox&);
+        friend Istream& operator>>(Istream& is, treeBoundBox&);
+        friend Ostream& operator<<(Ostream& os, const treeBoundBox&);
+
 };
 
 
+//- Specify data associated with boundBox type is contiguous
+template<>
+inline bool contiguous<treeBoundBox>() {return contiguous<boundBox>();}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
new file mode 100644
index 0000000000..5abc758ad3
--- /dev/null
+++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
@@ -0,0 +1,2210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "distributedTriSurfaceMesh.H"
+#include "mapDistribute.H"
+#include "Random.H"
+#include "addToRunTimeSelectionTable.H"
+#include "triangleFuncs.H"
+#include "matchPoints.H"
+#include "globalIndex.H"
+#include "Time.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(distributedTriSurfaceMesh, 0);
+addToRunTimeSelectionTable(searchableSurface, distributedTriSurfaceMesh, dict);
+
+scalar distributedTriSurfaceMesh::mergeDist_ = SMALL;
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Read my additional data from the dictionary
+bool Foam::distributedTriSurfaceMesh::read()
+{
+    // Get bb of all domains
+    procBb_.setSize(Pstream::nProcs());
+
+    procBb_[Pstream::myProcNo()] = List<treeBoundBox>(dict_.lookup("bounds"));
+    Pstream::gatherList(procBb_);
+    Pstream::scatterList(procBb_);
+
+    return true;
+}
+
+
+// Is segment fully local?
+bool Foam::distributedTriSurfaceMesh::isLocal
+(
+    const List<treeBoundBox>& myBbs,
+    const point& start,
+    const point& end
+)
+{
+    forAll(myBbs, bbI)
+    {
+        if (myBbs[bbI].contains(start) && myBbs[bbI].contains(end))
+        {
+            return true;
+        }
+    }
+    return false;
+}
+
+
+void Foam::distributedTriSurfaceMesh::splitSegment
+(
+    const label segmentI,
+    const point& start,
+    const point& end,
+
+    DynamicList<segment>& allSegments,
+    DynamicList<label>& allSegmentMap,
+    List<DynamicList<label> >& sendMap
+) const
+{
+    // Work points
+    point clipPt0, clipPt1;
+
+
+    // 1. Fully local already handled outside
+
+
+    // 2. Check if fully inside other processor. Rare occurrence
+    // but cheap to test.
+
+    forAll(procBb_, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            const List<treeBoundBox>& bbs = procBb_[procI];
+
+            forAll(bbs, bbI)
+            {
+                if (bbs[bbI].contains(start) && bbs[bbI].contains(end))
+                {
+                    //Pout<< "    Completely remote segment:"
+                    //    << start << end << " on proc:" << procI << endl;
+                    sendMap[procI].append(allSegments.size());
+                    allSegmentMap.append(segmentI);
+                    allSegments.append(segment(start, end));
+                    return;
+                }
+            }
+        }
+    }
+
+
+    // 3. Not contained in single processor. Clip against proc bbs.
+    forAll(procBb_, procI)
+    {
+        const List<treeBoundBox>& bbs = procBb_[procI];
+
+        forAll(bbs, bbI)
+        {
+            const treeBoundBox& bb = bbs[bbI];
+
+            if (bb.contains(start))
+            {
+                // start within, trim end to bb
+                bool clipped = bb.intersects(end, start, clipPt0);
+
+                if (clipped)
+                {
+                    //Pout<< "    Start of segment:"
+                    //    << start << end << " clips proc:" << procI
+                    //    << " at " << clipPt0 << endl;
+                    // segment from start to clippedStart passes
+                    // through proc.
+                    sendMap[procI].append(allSegments.size());
+                    allSegmentMap.append(segmentI);
+                    allSegments.append(segment(start, clipPt0));
+                }
+            }
+            else if (bb.contains(end))
+            {
+                // end within, trim start to bb
+                bool clipped = bb.intersects(start, end, clipPt0);
+
+                if (clipped)
+                {
+                    //Pout<< "    End of segment:"
+                    //    << start << end << " clips proc:" << procI
+                    //    << " at " << clipPt0 << endl;
+                    sendMap[procI].append(allSegments.size());
+                    allSegmentMap.append(segmentI);
+                    allSegments.append(segment(clipPt0, end));
+                }
+            }
+            else
+            {
+                // trim both
+                bool clippedStart = bb.intersects(start, end, clipPt0);
+
+                if (clippedStart)
+                {
+                    bool clippedEnd = bb.intersects(end, clipPt0, clipPt1);
+
+                    if (clippedEnd)
+                    {
+                        //Pout<< "    Middle of segment:"
+                        //    << start << end << " clips proc:" << procI
+                        //    << " at " << clipPt0 << clipPt1 << endl;
+                        // middle part of segment passes through
+                        // proc.
+                        sendMap[procI].append(allSegments.size());
+                        allSegmentMap.append(segmentI);
+                        allSegments.append(segment(clipPt0, clipPt1));
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+Foam::autoPtr<Foam::mapDistribute>
+Foam::distributedTriSurfaceMesh::constructSegments
+(
+    const pointField& start,
+    const pointField& end,
+
+    List<segment>& allSegments,
+    labelList& allSegmentMap
+) const
+{
+    // Determine send map
+    // ~~~~~~~~~~~~~~~~~~
+
+    labelListList sendMap(Pstream::nProcs());
+
+    {
+        // Since intersection test is quite expensive compared to memory
+        // allocation we use DynamicList to immediately store the segment
+        // in the correct bin.
+
+        // Segments to test
+        DynamicList<segment> dynAllSegments(start.size());
+        // Original index of segment
+        DynamicList<label> dynAllSegmentMap(start.size());
+        // Per processor indices into allSegments to send
+        List<DynamicList<label> > dynSendMap(Pstream::nProcs());
+
+        forAll(start, segmentI)
+        {
+            splitSegment
+            (
+                segmentI,
+                start[segmentI],
+                end[segmentI],
+
+                dynAllSegments,
+                dynAllSegmentMap,
+                dynSendMap
+            );
+        }
+
+        // Convert dynamicList to labelList
+        sendMap.setSize(Pstream::nProcs());
+        forAll(sendMap, procI)
+        {
+            dynSendMap[procI].shrink();
+            sendMap[procI].transfer(dynSendMap[procI]);
+        }
+
+        allSegments.transfer(dynAllSegments.shrink());
+        allSegmentMap.transfer(dynAllSegmentMap.shrink());
+    }
+
+
+    // Send over how many I need to receive.
+    labelListList sendSizes(Pstream::nProcs());
+    sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
+    forAll(sendMap, procI)
+    {
+        sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
+    }
+    Pstream::gatherList(sendSizes);
+    Pstream::scatterList(sendSizes);
+
+
+    // Determine order of receiving
+    labelListList constructMap(Pstream::nProcs());
+
+    // My local segments first
+    constructMap[Pstream::myProcNo()] = identity
+    (
+        sendMap[Pstream::myProcNo()].size()
+    );
+
+    label segmentI = constructMap[Pstream::myProcNo()].size();
+    forAll(constructMap, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            // What I need to receive is what other processor is sending to me.
+            label nRecv = sendSizes[procI][Pstream::myProcNo()];
+            constructMap[procI].setSize(nRecv);
+
+            for (label i = 0; i < nRecv; i++)
+            {
+                constructMap[procI][i] = segmentI++;
+            }
+        }
+    }
+
+    return autoPtr<mapDistribute>
+    (
+        new mapDistribute
+        (
+            segmentI,       // size after construction
+            sendMap,
+            constructMap,
+            true            // reuse storage
+        )
+    );
+}
+
+
+void Foam::distributedTriSurfaceMesh::findLine
+(
+    const bool nearestIntersection,
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    const indexedOctree<treeDataTriSurface>& octree = tree();
+
+    // Important:force synchronised construction of indexing
+    const globalIndex& triIndexer = globalTris();
+
+    // Initialise
+    info.setSize(start.size());
+    forAll(info, i)
+    {
+        info[i].setMiss();
+    }
+
+
+    // Do any local queries
+    // ~~~~~~~~~~~~~~~~~~~~
+
+    label nLocal = 0;
+
+    forAll(start, i)
+    {
+        if (isLocal(procBb_[Pstream::myProcNo()], start[i], end[i]))
+        {
+            if (nearestIntersection)
+            {
+                info[i] = octree.findLine(start[i], end[i]);
+            }
+            else
+            {
+                info[i] = octree.findLineAny(start[i], end[i]);
+            }
+
+            if (info[i].hit())
+            {
+                info[i].setIndex(triIndexer.toGlobal(info[i].index()));
+            }
+            nLocal++;
+        }
+    }
+
+
+    if
+    (
+        Pstream::parRun()
+     && (
+            returnReduce(nLocal, sumOp<label>())
+          < returnReduce(start.size(), sumOp<label>())
+        )
+    )
+    {
+        // Not all can be resolved locally. Build segments and map, send over
+        // segments, do intersections, send back and merge.
+
+
+        // Construct queries (segments)
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        // Segments to test
+        List<segment> allSegments(start.size());
+        // Original index of segment
+        labelList allSegmentMap(start.size());
+
+        const autoPtr<mapDistribute> mapPtr
+        (
+            constructSegments
+            (
+                start,
+                end,
+                allSegments,
+                allSegmentMap
+            )
+        );
+        const mapDistribute& map = mapPtr();
+
+        label nOldAllSegments = allSegments.size();
+
+
+        // Exchange the segments
+        // ~~~~~~~~~~~~~~~~~~~~~
+
+        map.distribute
+        (
+            Pstream::nonBlocking,   //Pstream::scheduled,
+            List<labelPair>(0),     //map.schedule(),
+            map.constructSize(),
+            map.subMap(),           // what to send
+            map.constructMap(),     // what to receive
+            allSegments
+        );
+
+
+        // Do tests I need to do
+        // ~~~~~~~~~~~~~~~~~~~~~
+
+        // Intersections
+        List<pointIndexHit> intersections(allSegments.size());
+
+        forAll(allSegments, i)
+        {
+            if (nearestIntersection)
+            {
+                intersections[i] = octree.findLine
+                (
+                    allSegments[i].first(),
+                    allSegments[i].second()
+                );
+            }
+            else
+            {
+                intersections[i] = octree.findLineAny
+                (
+                    allSegments[i].first(),
+                    allSegments[i].second()
+                );
+            }
+
+            // Convert triangle index to global numbering
+            if (intersections[i].hit())
+            {
+                intersections[i].setIndex
+                (
+                    triIndexer.toGlobal(intersections[i].index())
+                );
+            }
+        }
+
+
+        // Exchange the intersections (opposite to segments)
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        map.distribute
+        (
+            //Pstream::scheduled,
+            //map.schedule            // Note reverse schedule
+            //(
+            //    map.constructMap(),
+            //    map.subMap()
+            //),
+            Pstream::nonBlocking,
+            List<labelPair>(0),
+            nOldAllSegments,
+            map.constructMap(),     // what to send
+            map.subMap(),           // what to receive
+            intersections
+        );
+
+
+        // Extract the hits
+        // ~~~~~~~~~~~~~~~~
+
+        forAll(intersections, i)
+        {
+            const pointIndexHit& allInfo = intersections[i];
+            label segmentI = allSegmentMap[i];
+            pointIndexHit& hitInfo = info[segmentI];
+
+            if (allInfo.hit())
+            {
+                if (!hitInfo.hit())
+                {
+                    // No intersection yet so take this one
+                    hitInfo = allInfo;
+                }
+                else if (nearestIntersection)
+                {
+                    // Nearest intersection
+                    if
+                    (
+                        magSqr(allInfo.hitPoint()-start[segmentI])
+                      < magSqr(hitInfo.hitPoint()-start[segmentI])
+                    )
+                    {
+                        hitInfo = allInfo;
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+// Exchanges indices to the processor they come from.
+// - calculates exchange map
+// - uses map to calculate local triangle index
+Foam::autoPtr<Foam::mapDistribute>
+Foam::distributedTriSurfaceMesh::calcLocalQueries
+(
+    const List<pointIndexHit>& info,
+    labelList& triangleIndex
+) const
+{
+    triangleIndex.setSize(info.size());
+
+    const globalIndex& triIndexer = globalTris();
+
+
+    // Determine send map
+    // ~~~~~~~~~~~~~~~~~~
+
+    // Since determining which processor the query should go to is
+    // cheap we do a multi-pass algorithm to save some memory temporarily.
+
+    // 1. Count
+    labelList nSend(Pstream::nProcs(), 0);
+
+    forAll(info, i)
+    {
+        if (info[i].hit())
+        {
+            label procI = triIndexer.whichProcID(info[i].index());
+            nSend[procI]++;
+        }
+    }
+
+    // 2. Size sendMap
+    labelListList sendMap(Pstream::nProcs());
+    forAll(nSend, procI)
+    {
+        sendMap[procI].setSize(nSend[procI]);
+        nSend[procI] = 0;
+    }
+
+    // 3. Fill sendMap
+    forAll(info, i)
+    {
+        if (info[i].hit())
+        {
+            label procI = triIndexer.whichProcID(info[i].index());
+            triangleIndex[i] = triIndexer.toLocal(procI, info[i].index());
+            sendMap[procI][nSend[procI]++] = i;
+        }
+        else
+        {
+            triangleIndex[i] = -1;
+        }
+    }
+
+
+    // Send over how many I need to receive
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelListList sendSizes(Pstream::nProcs());
+    sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
+    forAll(sendMap, procI)
+    {
+        sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
+    }
+    Pstream::gatherList(sendSizes);
+    Pstream::scatterList(sendSizes);
+
+
+    // Determine receive map
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    labelListList constructMap(Pstream::nProcs());
+
+    // My local segments first
+    constructMap[Pstream::myProcNo()] = identity
+    (
+        sendMap[Pstream::myProcNo()].size()
+    );
+
+    label segmentI = constructMap[Pstream::myProcNo()].size();
+    forAll(constructMap, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            // What I need to receive is what other processor is sending to me.
+            label nRecv = sendSizes[procI][Pstream::myProcNo()];
+            constructMap[procI].setSize(nRecv);
+
+            for (label i = 0; i < nRecv; i++)
+            {
+                constructMap[procI][i] = segmentI++;
+            }
+        }
+    }
+
+
+    // Pack into distribution map
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    autoPtr<mapDistribute> mapPtr
+    (
+        new mapDistribute
+        (
+            segmentI,       // size after construction
+            sendMap,
+            constructMap,
+            true            // reuse storage
+        )
+    );
+    const mapDistribute& map = mapPtr();
+
+
+    // Send over queries
+    // ~~~~~~~~~~~~~~~~~
+
+    map.distribute
+    (
+        //Pstream::scheduled,
+        //map.schedule(),
+        Pstream::nonBlocking,
+        List<labelPair>(0),
+        map.constructSize(),
+        map.subMap(),           // what to send
+        map.constructMap(),     // what to receive
+        triangleIndex
+    );
+
+
+    return mapPtr;
+}
+
+
+Foam::label Foam::distributedTriSurfaceMesh::calcOverlappingProcs
+(
+    const point& centre,
+    const scalar radiusSqr,
+    boolList& overlaps
+) const
+{
+    overlaps = false;
+    label nOverlaps = 0;
+
+    forAll(procBb_, procI)
+    {
+        const List<treeBoundBox>& bbs = procBb_[procI];
+
+        forAll(bbs, bbI)
+        {
+            if (bbs[bbI].overlaps(centre, radiusSqr))
+            {
+                overlaps[procI] = true;
+                nOverlaps++;
+                break;
+            }
+        }
+    }
+    return nOverlaps;
+}
+
+
+// Generate queries for parallel distance calculation
+// - calculates exchange map
+// - uses map to exchange points and radius
+Foam::autoPtr<Foam::mapDistribute>
+Foam::distributedTriSurfaceMesh::calcLocalQueries
+(
+    const pointField& centres,
+    const scalarField& radiusSqr,
+
+    pointField& allCentres,
+    scalarField& allRadiusSqr,
+    labelList& allSegmentMap
+) const
+{
+    // Determine queries
+    // ~~~~~~~~~~~~~~~~~
+
+    labelListList sendMap(Pstream::nProcs());
+
+    {
+        // Queries
+        DynamicList<point> dynAllCentres(centres.size());
+        DynamicList<scalar> dynAllRadiusSqr(centres.size());
+        // Original index of segment
+        DynamicList<label> dynAllSegmentMap(centres.size());
+        // Per processor indices into allSegments to send
+        List<DynamicList<label> > dynSendMap(Pstream::nProcs());
+
+        // Work array - whether processor bb overlaps the bounding sphere.
+        boolList procBbOverlaps(Pstream::nProcs());
+
+        forAll(centres, centreI)
+        {
+            // Find the processor this sample+radius overlaps.
+            calcOverlappingProcs
+            (
+                centres[centreI],
+                radiusSqr[centreI],
+                procBbOverlaps
+            );
+
+            forAll(procBbOverlaps, procI)
+            {
+                if (procI != Pstream::myProcNo() && procBbOverlaps[procI])
+                {
+                    dynSendMap[procI].append(dynAllCentres.size());
+                    dynAllSegmentMap.append(centreI);
+                    dynAllCentres.append(centres[centreI]);
+                    dynAllRadiusSqr.append(radiusSqr[centreI]);
+                }
+            }
+        }
+
+        // Convert dynamicList to labelList
+        sendMap.setSize(Pstream::nProcs());
+        forAll(sendMap, procI)
+        {
+            dynSendMap[procI].shrink();
+            sendMap[procI].transfer(dynSendMap[procI]);
+        }
+
+        allCentres.transfer(dynAllCentres.shrink());
+        allRadiusSqr.transfer(dynAllRadiusSqr.shrink());
+        allSegmentMap.transfer(dynAllSegmentMap.shrink());
+    }
+
+
+    // Send over how many I need to receive.
+    labelListList sendSizes(Pstream::nProcs());
+    sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
+    forAll(sendMap, procI)
+    {
+        sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
+    }
+    Pstream::gatherList(sendSizes);
+    Pstream::scatterList(sendSizes);
+
+
+    // Determine order of receiving
+    labelListList constructMap(Pstream::nProcs());
+
+    // My local segments first
+    constructMap[Pstream::myProcNo()] = identity
+    (
+        sendMap[Pstream::myProcNo()].size()
+    );
+
+    label segmentI = constructMap[Pstream::myProcNo()].size();
+    forAll(constructMap, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            // What I need to receive is what other processor is sending to me.
+            label nRecv = sendSizes[procI][Pstream::myProcNo()];
+            constructMap[procI].setSize(nRecv);
+
+            for (label i = 0; i < nRecv; i++)
+            {
+                constructMap[procI][i] = segmentI++;
+            }
+        }
+    }
+
+    autoPtr<mapDistribute> mapPtr
+    (
+        new mapDistribute
+        (
+            segmentI,       // size after construction
+            sendMap,
+            constructMap,
+            true            // reuse storage
+        )
+    );
+    return mapPtr;
+}
+
+
+void Foam::distributedTriSurfaceMesh::calcBounds
+(
+    boundBox& bb,
+    label& nPoints
+) const
+{
+    // Unfortunately nPoints constructs meshPoints() so do compact version
+    // ourselves
+
+    PackedList<1> pointIsUsed(points().size());
+    pointIsUsed = 0U;
+
+    nPoints = 0;
+    bb.min() = point(VGREAT, VGREAT, VGREAT);
+    bb.max() = point(-VGREAT, -VGREAT, -VGREAT);
+
+    const triSurface& s = static_cast<const triSurface&>(*this);
+
+    forAll(s, triI)
+    {
+        const labelledTri& f = s[triI];
+
+        forAll(f, fp)
+        {
+            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++;
+            }
+        }
+    }
+}
+
+
+// Does any part of triangle overlap bb.
+bool Foam::distributedTriSurfaceMesh::overlaps
+(
+    const List<treeBoundBox>& bbs,
+    const point& p0,
+    const point& p1,
+    const point& p2
+)
+{
+    forAll(bbs, bbI)
+    {
+        const treeBoundBox& bb = bbs[bbI];
+
+        boundBox triBb(p0, p0);
+        triBb.min() = min(triBb.min(), p1);
+        triBb.min() = min(triBb.min(), p2);
+
+        triBb.max() = max(triBb.max(), p1);
+        triBb.max() = max(triBb.max(), p2);
+
+        //- Exact test of triangle intersecting bb
+
+        // Quick rejection. If whole bounding box of tri is outside cubeBb then
+        // there will be no intersection.
+        if (bb.overlaps(triBb))
+        {
+            // Check if one or more triangle point inside
+            if (bb.contains(p0) || bb.contains(p1) || bb.contains(p2))
+            {
+                // One or more points inside
+                return true;
+            }
+
+            // Now we have the difficult case: all points are outside but
+            // connecting edges might go through cube. Use fast intersection
+            // of bounding box.
+            bool intersect = triangleFuncs::intersectBb(p0, p1, p2, bb);
+
+            if (intersect)
+            {
+                return true;
+            }
+        }
+    }
+    return false;
+}
+
+
+void Foam::distributedTriSurfaceMesh::subsetMeshMap
+(
+    const triSurface& s,
+    const boolList& include,
+    const label nIncluded,
+    labelList& newToOldPoints,
+    labelList& oldToNewPoints,
+    labelList& newToOldFaces
+)
+{
+    newToOldFaces.setSize(nIncluded);
+    newToOldPoints.setSize(s.points().size());
+    oldToNewPoints.setSize(s.points().size());
+    oldToNewPoints = -1;
+    {
+        label faceI = 0;
+        label pointI = 0;
+
+        forAll(include, oldFacei)
+        {
+            if (include[oldFacei])
+            {
+                // Store new faces compact
+                newToOldFaces[faceI++] = oldFacei;
+
+                // Renumber labels for triangle
+                const labelledTri& tri = s[oldFacei];
+
+                forAll(tri, fp)
+                {
+                    label oldPointI = tri[fp];
+
+                    if (oldToNewPoints[oldPointI] == -1)
+                    {
+                        oldToNewPoints[oldPointI] = pointI;
+                        newToOldPoints[pointI++] = oldPointI;
+                    }
+                }
+            }
+        }
+        newToOldPoints.setSize(pointI);
+    }
+}
+
+
+Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
+(
+    const triSurface& s,
+    const labelList& newToOldPoints,
+    const labelList& oldToNewPoints,
+    const labelList& newToOldFaces
+)
+{
+    // Extract points
+    pointField newPoints(newToOldPoints.size());
+    forAll(newToOldPoints, i)
+    {
+        newPoints[i] = s.points()[newToOldPoints[i]];
+    }
+    // Extract faces
+    List<labelledTri> newTriangles(newToOldFaces.size());
+
+    forAll(newToOldFaces, i)
+    {
+        // Get old vertex labels
+        const labelledTri& tri = s[newToOldFaces[i]];
+
+        newTriangles[i][0] = oldToNewPoints[tri[0]];
+        newTriangles[i][1] = oldToNewPoints[tri[1]];
+        newTriangles[i][2] = oldToNewPoints[tri[2]];
+        newTriangles[i].region() = tri.region();
+    }
+
+    // Reuse storage.
+    return triSurface(newTriangles, s.patches(), newPoints, true);
+}
+
+
+Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
+(
+    const triSurface& s,
+    const boolList& include,
+    labelList& newToOldPoints,
+    labelList& newToOldFaces
+)
+{
+    label n = 0;
+
+    forAll(include, i)
+    {
+        if (include[i])
+        {
+            n++;
+        }
+    }
+
+    labelList oldToNewPoints;
+    subsetMeshMap
+    (
+        s,
+        include,
+        n,
+        newToOldPoints,
+        oldToNewPoints,
+        newToOldFaces
+    );
+
+    return subsetMesh
+    (
+        s,
+        newToOldPoints,
+        oldToNewPoints,
+        newToOldFaces
+    );
+}
+
+
+Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
+(
+    const triSurface& s,
+    const labelList& newToOldFaces,
+    labelList& newToOldPoints
+)
+{
+    const boolList include
+    (
+        createWithValues<boolList>
+        (
+            s.size(),
+            false,
+            newToOldFaces,
+            true
+        )
+    );
+
+    newToOldPoints.setSize(s.points().size());
+    labelList oldToNewPoints(s.points().size(), -1);
+    {
+        label pointI = 0;
+
+        forAll(include, oldFacei)
+        {
+            if (include[oldFacei])
+            {
+                // Renumber labels for triangle
+                const labelledTri& tri = s[oldFacei];
+
+                forAll(tri, fp)
+                {
+                    label oldPointI = tri[fp];
+
+                    if (oldToNewPoints[oldPointI] == -1)
+                    {
+                        oldToNewPoints[oldPointI] = pointI;
+                        newToOldPoints[pointI++] = oldPointI;
+                    }
+                }
+            }
+        }
+        newToOldPoints.setSize(pointI);
+    }
+
+    return subsetMesh
+    (
+        s,
+        newToOldPoints,
+        oldToNewPoints,
+        newToOldFaces
+    );
+}
+
+
+Foam::label Foam::distributedTriSurfaceMesh::findTriangle
+(
+    const List<labelledTri>& allFaces,
+    const labelListList& allPointFaces,
+    const labelledTri& otherF
+)
+{
+    // allFaces connected to otherF[0]
+    const labelList& pFaces = allPointFaces[otherF[0]];
+
+    forAll(pFaces, i)
+    {
+        const labelledTri& f = allFaces[pFaces[i]];
+
+        if (f.region() == otherF.region())
+        {
+            // Find index of otherF[0]
+            label fp0 = findIndex(f, otherF[0]);
+            // Check rest of triangle in same order
+            label fp1 = f.fcIndex(fp0);
+            label fp2 = f.fcIndex(fp1);
+
+            if (f[fp1] == otherF[1] && f[fp2] == otherF[2])
+            {
+                return pFaces[i];
+            }
+        }
+    }
+    return -1;
+}
+
+
+// Merge into allSurf.
+void Foam::distributedTriSurfaceMesh::merge
+(
+    const scalar mergeDist,
+    const List<labelledTri>& subTris,
+    const pointField& subPoints,
+
+    List<labelledTri>& allTris,
+    pointField& allPoints,
+
+    labelList& faceConstructMap,
+    labelList& pointConstructMap
+)
+{
+    labelList subToAll;
+    matchPoints
+    (
+        subPoints,
+        allPoints,
+        scalarField(subPoints.size(), mergeDist),   // match distance
+        false,                                      // verbose
+        pointConstructMap
+    );
+
+    label nOldAllPoints = allPoints.size();
+
+
+    // Add all unmatched points
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    label allPointI = nOldAllPoints;
+    forAll(pointConstructMap, pointI)
+    {
+        if (pointConstructMap[pointI] == -1)
+        {
+            pointConstructMap[pointI] = allPointI++;
+        }
+    }
+
+    if (allPointI > nOldAllPoints)
+    {
+        allPoints.setSize(allPointI);
+
+        forAll(pointConstructMap, pointI)
+        {
+            if (pointConstructMap[pointI] >= nOldAllPoints)
+            {
+                allPoints[pointConstructMap[pointI]] = subPoints[pointI];
+            }
+        }
+    }
+
+
+    // To check whether triangles are same we use pointFaces.
+    labelListList allPointFaces;
+    invertManyToMany(nOldAllPoints, allTris, allPointFaces);
+
+
+    // Add all unmatched triangles
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    label allTriI = allTris.size();
+    allTris.setSize(allTriI + subTris.size());
+
+    faceConstructMap.setSize(subTris.size());
+
+    forAll(subTris, triI)
+    {
+        const labelledTri& subTri = subTris[triI];
+
+        // Get triangle in new numbering
+        labelledTri mappedTri
+        (
+            pointConstructMap[subTri[0]],
+            pointConstructMap[subTri[1]],
+            pointConstructMap[subTri[2]],
+            subTri.region()
+        );
+
+
+        // Check if all points were already in surface
+        bool fullMatch = true;
+
+        forAll(mappedTri, fp)
+        {
+            if (mappedTri[fp] >= nOldAllPoints)
+            {
+                fullMatch = false;
+                break;
+            }
+        }
+
+        if (fullMatch)
+        {
+            // All three points are mapped to old points. See if same
+            // triangle.
+            label i = findTriangle
+            (
+                allTris,
+                allPointFaces,
+                mappedTri
+            );
+
+            if (i == -1)
+            {
+                // Add
+                faceConstructMap[triI] = allTriI;
+                allTris[allTriI] = mappedTri;
+                allTriI++;
+            }
+            else
+            {
+                faceConstructMap[triI] = i;
+            }
+        }
+        else
+        {
+            // Add
+            faceConstructMap[triI] = allTriI;
+            allTris[allTriI] = mappedTri;
+            allTriI++;
+        }
+    }
+    allTris.setSize(allTriI);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
+(
+    const IOobject& io,
+    const triSurface& s,
+    const dictionary& dict
+)
+:
+    triSurfaceMesh(io, s),
+    dict_(io, dict)
+{
+    read();
+}
+
+
+Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
+:
+    triSurfaceMesh(io),
+    dict_
+    (
+        IOobject
+        (
+            io.name() + "Dict",
+            io.instance(),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            io.registerObject()
+        )
+    )
+{
+    read();
+}
+
+
+Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    triSurfaceMesh(io, dict),
+    dict_
+    (
+        IOobject
+        (
+            io.name() + "Dict",
+            io.instance(),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            io.registerObject()
+        )
+    )
+{
+    read();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::distributedTriSurfaceMesh::~distributedTriSurfaceMesh()
+{
+    clearOut();
+}
+
+
+void Foam::distributedTriSurfaceMesh::clearOut()
+{
+    globalTris_.clear();
+    triSurfaceMesh::clearOut();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::globalIndex& Foam::distributedTriSurfaceMesh::globalTris() const
+{
+    if (!globalTris_.valid())
+    {
+        globalTris_.reset(new globalIndex(triSurface::size()));
+    }
+    return globalTris_;
+}
+
+
+void Foam::distributedTriSurfaceMesh::findNearest
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    List<pointIndexHit>& info
+) const
+{
+    const indexedOctree<treeDataTriSurface>& octree = tree();
+
+    // Important:force synchronised construction of indexing
+    const globalIndex& triIndexer = globalTris();
+
+
+    // Initialise
+    // ~~~~~~~~~~
+
+    info.setSize(samples.size());
+    forAll(info, i)
+    {
+        info[i].setMiss();
+    }
+
+
+
+    // Do any local queries
+    // ~~~~~~~~~~~~~~~~~~~~
+
+    label nLocal = 0;
+
+    {
+        // Work array - whether processor bb overlaps the bounding sphere.
+        boolList procBbOverlaps(Pstream::nProcs());
+
+        forAll(samples, i)
+        {
+            // Find the processor this sample+radius overlaps.
+            label nProcs = calcOverlappingProcs
+            (
+                samples[i],
+                nearestDistSqr[i],
+                procBbOverlaps
+            );
+
+            // Overlaps local processor?
+            if (procBbOverlaps[Pstream::myProcNo()])
+            {
+                info[i] = octree.findNearest(samples[i], nearestDistSqr[i]);
+                if (info[i].hit())
+                {
+                    info[i].setIndex(triIndexer.toGlobal(info[i].index()));
+                }
+                if (nProcs == 1)
+                {
+                    // Fully local
+                    nLocal++;
+                }
+            }
+        }
+    }
+
+
+    if
+    (
+        Pstream::parRun()
+     && (
+            returnReduce(nLocal, sumOp<label>())
+          < returnReduce(samples.size(), sumOp<label>())
+        )
+    )
+    {
+        // Not all can be resolved locally. Build queries and map, send over
+        // queries, do intersections, send back and merge.
+
+        // Calculate queries and exchange map
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        pointField allCentres;
+        scalarField allRadiusSqr;
+        labelList allSegmentMap;
+        autoPtr<mapDistribute> mapPtr
+        (
+            calcLocalQueries
+            (
+                samples,
+                nearestDistSqr,
+
+                allCentres,
+                allRadiusSqr,
+                allSegmentMap
+            )
+        );
+        const mapDistribute& map = mapPtr();
+
+
+        // swap samples to local processor
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        map.distribute
+        (
+            //Pstream::scheduled,
+            //map.schedule(),
+            Pstream::nonBlocking,
+            List<labelPair>(0),
+            map.constructSize(),
+            map.subMap(),           // what to send
+            map.constructMap(),     // what to receive
+            allCentres
+        );
+        map.distribute
+        (
+            //Pstream::scheduled,
+            //map.schedule(),
+            Pstream::nonBlocking,
+            List<labelPair>(0),
+            map.constructSize(),
+            map.subMap(),           // what to send
+            map.constructMap(),     // what to receive
+            allRadiusSqr
+        );
+
+
+        // Do my tests
+        // ~~~~~~~~~~~
+
+        List<pointIndexHit> allInfo(allCentres.size());
+        forAll(allInfo, i)
+        {
+            allInfo[i] = octree.findNearest
+            (
+                allCentres[i],
+                allRadiusSqr[i]
+            );
+            if (allInfo[i].hit())
+            {
+                allInfo[i].setIndex(triIndexer.toGlobal(allInfo[i].index()));
+            }
+        }
+
+
+        // Send back results
+        // ~~~~~~~~~~~~~~~~~
+
+        map.distribute
+        (
+            //Pstream::scheduled,
+            //map.schedule            // note reverse schedule
+            //(
+            //    map.constructMap(),
+            //    map.subMap()
+            //),
+            Pstream::nonBlocking,
+            List<labelPair>(0),
+            allSegmentMap.size(),
+            map.constructMap(),     // what to send
+            map.subMap(),           // what to receive
+            allInfo
+        );
+
+
+        // Extract information
+        // ~~~~~~~~~~~~~~~~~~~
+
+        forAll(allInfo, i)
+        {
+            if (allInfo[i].hit())
+            {
+                label pointI = allSegmentMap[i];
+
+                if (!info[pointI].hit())
+                {
+                    // No intersection yet so take this one
+                    info[pointI] = allInfo[i];
+                }
+                else
+                {
+                    // Nearest intersection
+                    if
+                    (
+                        magSqr(allInfo[i].hitPoint()-samples[pointI])
+                      < magSqr(info[pointI].hitPoint()-samples[pointI])
+                    )
+                    {
+                        info[pointI] = allInfo[i];
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::distributedTriSurfaceMesh::findLine
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    findLine
+    (
+        true,   // nearestIntersection
+        start,
+        end,
+        info
+    );
+}
+
+
+void Foam::distributedTriSurfaceMesh::findLineAny
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    findLine
+    (
+        true,   // nearestIntersection
+        start,
+        end,
+        info
+    );
+}
+
+
+void Foam::distributedTriSurfaceMesh::findLineAll
+(
+    const pointField& start,
+    const pointField& end,
+    List<List<pointIndexHit> >& info
+) const
+{
+    // Reuse fineLine. We could modify all of findLine to do multiple
+    // intersections but this would mean a lot of data transferred so
+    // for now we just find nearest intersection and retest from that
+    // intersection onwards.
+
+    // Work array.
+    List<pointIndexHit> hitInfo(start.size());
+
+    findLine
+    (
+        true,   // nearestIntersection
+        start,
+        end,
+        hitInfo
+    );
+
+    // 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)
+    );
+
+    // Copy to input and compact any hits
+    labelList pointMap(start.size());
+    pointField e0(start.size());
+    pointField e1(start.size());
+    label compactI = 0;
+
+    info.setSize(hitInfo.size());
+    forAll(hitInfo, pointI)
+    {
+        if (hitInfo[pointI].hit())
+        {
+            info[pointI].setSize(1);
+            info[pointI][0] = hitInfo[pointI];
+
+            point pt = hitInfo[pointI].hitPoint() + smallVec[pointI];
+
+            if (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
+            {
+                e0[compactI] = pt;
+                e1[compactI] = end[pointI];
+                pointMap[compactI] = pointI;
+                compactI++;
+            }
+        }
+        else
+        {
+            info[pointI].clear();
+        }
+    }
+
+    e0.setSize(compactI);
+    e1.setSize(compactI);
+    pointMap.setSize(compactI);
+
+    while (returnReduce(e0.size(), sumOp<label>()) > 0)
+    {
+        findLine
+        (
+            true,   // nearestIntersection
+            e0,
+            e1,
+            hitInfo
+        );
+
+
+        // Extract
+        label compactI = 0;
+        forAll(hitInfo, i)
+        {
+            if (hitInfo[i].hit())
+            {
+                label pointI = pointMap[i];
+
+                label sz = info[pointI].size();
+                info[pointI].setSize(sz+1);
+                info[pointI][sz] = hitInfo[i];
+
+                point pt = hitInfo[i].hitPoint() + smallVec[pointI];
+
+                if (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
+                {
+                    e0[compactI] = pt;
+                    e1[compactI] = end[pointI];
+                    pointMap[compactI] = pointI;
+                    compactI++;
+                }
+            }
+        }
+
+        // Trim
+        e0.setSize(compactI);
+        e1.setSize(compactI);
+        pointMap.setSize(compactI);
+    }
+}
+
+
+void Foam::distributedTriSurfaceMesh::getRegion
+(
+    const List<pointIndexHit>& info,
+    labelList& region
+) const
+{
+    if (!Pstream::parRun())
+    {
+        region.setSize(info.size());
+        forAll(info, i)
+        {
+            if (info[i].hit())
+            {
+                region[i] = triSurface::operator[](info[i].index()).region();
+            }
+            else
+            {
+                region[i] = -1;
+            }
+        }
+        return;
+    }
+
+
+    // Get query data (= local index of triangle)
+    // ~~~~~~~~~~~~~~
+
+    labelList triangleIndex(info.size());
+    autoPtr<mapDistribute> mapPtr
+    (
+        calcLocalQueries
+        (
+            info,
+            triangleIndex
+        )
+    );
+    const mapDistribute& map = mapPtr();
+
+
+    // Do my tests
+    // ~~~~~~~~~~~
+
+    const triSurface& s = static_cast<const triSurface&>(*this);
+
+    region.setSize(triangleIndex.size());
+
+    forAll(triangleIndex, i)
+    {
+        label triI = triangleIndex[i];
+        region[i] = s[triI].region();
+    }
+
+
+    // Send back results
+    // ~~~~~~~~~~~~~~~~~
+
+    map.distribute
+    (
+        //Pstream::scheduled,
+        //map.schedule            // note reverse schedule
+        //(
+        //    map.constructMap(),
+        //    map.subMap()
+        //),
+        Pstream::nonBlocking,
+        List<labelPair>(0),
+        info.size(),
+        map.constructMap(),     // what to send
+        map.subMap(),           // what to receive
+        region
+    );
+}
+
+
+void Foam::distributedTriSurfaceMesh::getNormal
+(
+    const List<pointIndexHit>& info,
+    vectorField& normal
+) const
+{
+    if (!Pstream::parRun())
+    {
+        normal.setSize(info.size());
+
+        forAll(info, i)
+        {
+            if (info[i].hit())
+            {
+                normal[i] = faceNormals()[info[i].index()];
+            }
+            else
+            {
+                // Set to what?
+                normal[i] = vector::zero;
+            }
+        }
+        return;
+    }
+
+
+    // Get query data (= local index of triangle)
+    // ~~~~~~~~~~~~~~
+
+    labelList triangleIndex(info.size());
+    autoPtr<mapDistribute> mapPtr
+    (
+        calcLocalQueries
+        (
+            info,
+            triangleIndex
+        )
+    );
+    const mapDistribute& map = mapPtr();
+
+
+    // Do my tests
+    // ~~~~~~~~~~~
+
+    const triSurface& s = static_cast<const triSurface&>(*this);
+
+    normal.setSize(triangleIndex.size());
+
+    forAll(triangleIndex, i)
+    {
+        label triI = triangleIndex[i];
+        normal[i] = s[triI].normal(s.points());
+        normal[i] /= mag(normal[i]) + VSMALL;
+    }
+
+
+    // Send back results
+    // ~~~~~~~~~~~~~~~~~
+
+    map.distribute
+    (
+        //Pstream::scheduled,
+        //map.schedule            // note reverse schedule
+        //(
+        //    map.constructMap(),
+        //    map.subMap()
+        //),
+        Pstream::nonBlocking,
+        List<labelPair>(0),
+        info.size(),
+        map.constructMap(),     // what to send
+        map.subMap(),           // what to receive
+        normal
+    );
+}
+
+
+void Foam::distributedTriSurfaceMesh::getField
+(
+    const word& fieldName,
+    const List<pointIndexHit>& info,
+    labelList& values
+) const
+{
+    const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
+    (
+        fieldName
+    );
+
+
+    if (!Pstream::parRun())
+    {
+        values.setSize(info.size());
+        forAll(info, i)
+        {
+            if (info[i].hit())
+            {
+                values[i] = fld[info[i].index()];
+            }
+        }
+        return;
+    }
+
+
+    // Get query data (= local index of triangle)
+    // ~~~~~~~~~~~~~~
+
+    labelList triangleIndex(info.size());
+    autoPtr<mapDistribute> mapPtr
+    (
+        calcLocalQueries
+        (
+            info,
+            triangleIndex
+        )
+    );
+    const mapDistribute& map = mapPtr();
+
+
+    // Do my tests
+    // ~~~~~~~~~~~
+
+    values.setSize(triangleIndex.size());
+
+    forAll(triangleIndex, i)
+    {
+        label triI = triangleIndex[i];
+        values[i] = fld[triI];
+    }
+
+
+    // Send back results
+    // ~~~~~~~~~~~~~~~~~
+
+    map.distribute
+    (
+        Pstream::nonBlocking,
+        List<labelPair>(0),
+        info.size(),
+        map.constructMap(),     // what to send
+        map.subMap(),           // what to receive
+        values
+    );
+}
+
+
+void Foam::distributedTriSurfaceMesh::getVolumeType
+(
+    const pointField& points,
+    List<volumeType>& volType
+) const
+{
+    FatalErrorIn
+    (
+        "distributedTriSurfaceMesh::getVolumeType"
+        "(const pointField&, List<volumeType>&) const"
+    )   << "Volume type not supported for distributed surfaces."
+        << exit(FatalError);
+}
+
+
+// Subset the part of surface that is overlapping bb.
+Foam::triSurface Foam::distributedTriSurfaceMesh::overlappingSurface
+(
+    const triSurface& s,
+    const List<treeBoundBox>& bbs,
+
+    labelList& subPointMap,
+    labelList& subFaceMap
+)
+{
+    // Determine what triangles to keep.
+    boolList includedFace(s.size(), false);
+
+    forAll(s, triI)
+    {
+        const labelledTri& f = s[triI];
+        const point& p0 = s.points()[f[0]];
+        const point& p1 = s.points()[f[1]];
+        const point& p2 = s.points()[f[2]];
+
+        if (overlaps(bbs, p0, p1, p2))
+        {
+            includedFace[triI] = true;
+        }
+    }
+
+    return subsetMesh(s, includedFace, subPointMap, subFaceMap);
+}
+
+
+void Foam::distributedTriSurfaceMesh::distribute
+(
+    const List<treeBoundBox>& bbs,
+    const bool keepNonLocal,
+    autoPtr<mapDistribute>& faceMap,
+    autoPtr<mapDistribute>& pointMap
+)
+{
+    // Get bb of all domains
+    {
+        List<List<treeBoundBox> > newProcBb(Pstream::nProcs());
+        newProcBb[Pstream::myProcNo()].setSize(bbs.size());
+        forAll(bbs, i)
+        {
+            newProcBb[Pstream::myProcNo()][i] = bbs[i];
+        }
+        Pstream::gatherList(newProcBb);
+        Pstream::scatterList(newProcBb);
+
+
+        //if (debug)
+        //{
+        //    Info<< "old bb:" << procBb_ << endl << endl;
+        //    Info<< "new bb:" << newProcBb << endl << endl;
+        //    Info<< "Same:" << (newProcBb == procBb_) << endl;
+        //}
+
+        if (newProcBb == procBb_)
+        {
+            return;
+        }
+        else
+        {
+            procBb_.transfer(newProcBb);
+        }
+    }
+
+
+    // Debug information
+    if (debug)
+    {
+        labelList nTris(Pstream::nProcs());
+        nTris[Pstream::myProcNo()] = triSurface::size();
+        Pstream::gatherList(nTris);
+        Pstream::scatterList(nTris);
+
+        Info<< "distributedTriSurfaceMesh::distribute : before distribution:"
+            << endl
+            << "\tproc\ttris" << endl;
+
+        forAll(nTris, procI)
+        {
+            Info<< '\t' << procI << '\t' << nTris[procI] << endl;
+        }
+        Info<< endl;
+    }
+
+
+    labelListList faceSendMap(Pstream::nProcs());
+    labelListList pointSendMap(Pstream::nProcs());
+
+    forAll(procBb_, procI)
+    {
+        overlappingSurface
+        (
+            *this,
+            procBb_[procI],
+            pointSendMap[procI],
+            faceSendMap[procI]
+        );
+
+        if (debug)
+        {
+            //Pout<< "Overlapping with proc " << procI
+            //    << " faces:" << faceSendMap[procI].size()
+            //    << " points:" << pointSendMap[procI].size() << endl << endl;
+        }
+    }
+
+    if (keepNonLocal)
+    {
+        // Include in faceSendMap/pointSendMap the triangles that are
+        // not mapped to any processor so they stay local.
+
+        const triSurface& s = static_cast<const triSurface&>(*this);
+
+        boolList includedFace(s.size(), true);
+
+        forAll(faceSendMap, procI)
+        {
+            if (procI != Pstream::myProcNo())
+            {
+                forAll(faceSendMap[procI], i)
+                {
+                    includedFace[faceSendMap[procI][i]] = false;
+                }
+            }
+        }
+
+        // Combine includedFace (all triangles that are not on any neighbour)
+        // with overlap.
+
+        forAll(faceSendMap[Pstream::myProcNo()], i)
+        {
+            includedFace[faceSendMap[Pstream::myProcNo()][i]] = true;
+        }
+
+        subsetMesh
+        (
+            s,
+            includedFace,
+            pointSendMap[Pstream::myProcNo()],
+            faceSendMap[Pstream::myProcNo()]
+        );
+    }
+
+
+    // Send over how many faces/points I need to receive
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelListList faceSendSizes(Pstream::nProcs());
+    faceSendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
+    forAll(faceSendMap, procI)
+    {
+        faceSendSizes[Pstream::myProcNo()][procI] = faceSendMap[procI].size();
+    }
+    Pstream::gatherList(faceSendSizes);
+    Pstream::scatterList(faceSendSizes);
+
+
+
+    // Exchange surfaces
+    // ~~~~~~~~~~~~~~~~~
+
+    // Storage for resulting surface
+    List<labelledTri> allTris;
+    pointField allPoints;
+
+    labelListList faceConstructMap(Pstream::nProcs());
+    labelListList pointConstructMap(Pstream::nProcs());
+
+
+    // My own surface first
+    // ~~~~~~~~~~~~~~~~~~~~
+
+    {
+        labelList pointMap;
+        triSurface subSurface
+        (
+            subsetMesh
+            (
+                *this,
+                faceSendMap[Pstream::myProcNo()],
+                pointMap
+            )
+        );
+
+        allTris = subSurface;
+        allPoints = subSurface.points();
+
+        faceConstructMap[Pstream::myProcNo()] = identity
+        (
+            faceSendMap[Pstream::myProcNo()].size()
+        );
+        pointConstructMap[Pstream::myProcNo()] = identity
+        (
+            pointSendMap[Pstream::myProcNo()].size()
+        );
+    }
+
+
+
+    // Send all
+    // ~~~~~~~~
+
+    forAll(faceSendSizes, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            if (faceSendSizes[Pstream::myProcNo()][procI] > 0)
+            {
+                OPstream str(Pstream::blocking, procI);
+
+                labelList pointMap;
+                triSurface subSurface
+                (
+                    subsetMesh
+                    (
+                        *this,
+                        faceSendMap[procI],
+                        pointMap
+                    )
+                );
+
+                //if (debug)
+                //{
+                //    Pout<< "Sending to " << procI
+                //        << " faces:" << faceSendMap[procI].size()
+                //        << " points:" << subSurface.points().size() << endl
+                //        << endl;
+                //}
+
+                str << subSurface;
+           }
+        }
+    }
+
+
+    // Receive and merge all
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    forAll(faceSendSizes, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            if (faceSendSizes[procI][Pstream::myProcNo()] > 0)
+            {
+                IPstream str(Pstream::blocking, procI);
+
+                // Receive
+                triSurface subSurface(str);
+
+                //if (debug)
+                //{
+                //    Pout<< "Received from " << procI
+                //        << " faces:" << subSurface.size()
+                //        << " points:" << subSurface.points().size() << endl
+                //        << endl;
+                //}
+
+                // Merge into allSurf
+                merge
+                (
+                    mergeDist_,
+                    subSurface,
+                    subSurface.points(),
+
+                    allTris,
+                    allPoints,
+                    faceConstructMap[procI],
+                    pointConstructMap[procI]
+                );
+
+                //if (debug)
+                //{
+                //    Pout<< "Current merged surface : faces:" << allTris.size()
+                //        << " points:" << allPoints.size() << endl << endl;
+                //}
+           }
+        }
+    }
+    
+
+    faceMap.reset
+    (
+        new mapDistribute
+        (
+            allTris.size(),
+            faceSendMap,
+            faceConstructMap,
+            true
+        )
+    );
+    pointMap.reset
+    (
+        new mapDistribute
+        (
+            allPoints.size(),
+            pointSendMap,
+            pointConstructMap,
+            true
+        )
+    );
+
+    // Construct triSurface. Reuse storage.
+    triSurface::operator=(triSurface(allTris, patches(), allPoints, true));
+
+    clearOut();
+
+    // Regions stays same
+    // Volume type stays same.
+
+    distributeFields<label>(faceMap());
+    distributeFields<scalar>(faceMap());
+    distributeFields<vector>(faceMap());
+    distributeFields<sphericalTensor>(faceMap());
+    distributeFields<symmTensor>(faceMap());
+    distributeFields<tensor>(faceMap());
+
+    if (debug)
+    {
+        labelList nTris(Pstream::nProcs());
+        nTris[Pstream::myProcNo()] = triSurface::size();
+        Pstream::gatherList(nTris);
+        Pstream::scatterList(nTris);
+
+        Info<< "distributedTriSurfaceMesh::distribute : after distribution:"
+            << endl
+            << "\tproc\ttris" << endl;
+
+        forAll(nTris, procI)
+        {
+            Info<< '\t' << procI << '\t' << nTris[procI] << endl;
+        }
+        Info<< endl;
+    }
+}
+
+
+//- Write using given format, version and compression
+bool Foam::distributedTriSurfaceMesh::writeObject
+(
+    IOstream::streamFormat fmt,
+    IOstream::versionNumber ver,
+    IOstream::compressionType cmp
+) const
+{
+    return
+        triSurfaceMesh::writeObject(fmt, ver, cmp)
+     && dict_.writeObject(fmt, ver, cmp);
+}
+
+
+void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
+{
+    boundBox bb;
+    label nPoints;
+    calcBounds(bb, nPoints);
+
+    os  << "Triangles    : " << returnReduce(triSurface::size(), sumOp<label>())
+        << endl
+        << "Vertices     : " << returnReduce(nPoints, sumOp<label>()) << endl
+        << "Bounding Box : " << bb << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H
new file mode 100644
index 0000000000..0fb5fa546e
--- /dev/null
+++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H
@@ -0,0 +1,459 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::distributedTriSurfaceMesh
+
+Description
+    IOoject and searching on triSurface
+
+SourceFiles
+    distributedTriSurfaceMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef distributedTriSurfaceMesh_H
+#define distributedTriSurfaceMesh_H
+
+#include "triSurfaceMesh.H"
+#include "IOdictionary.H"
+#include "Pair.H"
+#include "globalIndex.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class mapDistribute;
+
+// Typedefs
+typedef Pair<point> segment;
+template<>
+inline bool contiguous<segment>() {return contiguous<point>();}
+
+
+/*---------------------------------------------------------------------------*\
+                           Class distributedTriSurfaceMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class distributedTriSurfaceMesh
+:
+    public triSurfaceMesh
+{
+private:
+
+    // Static data
+
+        //- Merging distance
+        static scalar mergeDist_;
+
+
+    // Private member data
+
+        //- Bounding box settings
+        IOdictionary dict_;
+
+        //- bounding boxes of all processors
+        List<List<treeBoundBox> > procBb_;
+
+        //- Global triangle numbering
+        mutable autoPtr<globalIndex> globalTris_;
+
+
+    // Private Member Functions
+
+        // Read
+
+            //- Read my additional data
+            bool read();
+
+
+        // Line intersection
+
+            static bool isLocal
+            (
+                const List<treeBoundBox>& myBbs,
+                const point& start,
+                const point& end
+            );
+
+            //- Split segment into subsegments overlapping the processor
+            //  bounding boxes. Sort per processor.
+            void splitSegment
+            (
+                const label,
+                const point& start,
+                const point& end,
+
+                DynamicList<segment>&,
+                DynamicList<label>&,
+                List<DynamicList<label> >&
+            ) const;
+
+            //- Divide edges into local and remote segments. Construct map to
+            //  distribute and collect data.
+            autoPtr<mapDistribute> constructSegments
+            (
+                const pointField& start,
+                const pointField& end,
+
+                List<segment>& allSegments,
+                List<label>& allSegmentMap
+            ) const;
+
+            //- Split edges, distribute, test and collect.
+            void findLine
+            (
+                const bool nearestIntersection,
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>& info
+            ) const;
+
+
+        // Triangle index
+
+            //- Obtains global indices from pointIndexHit and swaps them back
+            //  to their original processor. Used to calculate local region
+            //  and normal.
+            autoPtr<mapDistribute> calcLocalQueries
+            (
+                const List<pointIndexHit>&,
+                labelList& triangleIndex
+            ) const;
+
+
+        // Nearest
+
+            label calcOverlappingProcs
+            (
+                const point& centre,
+                const scalar radiusSqr,
+                boolList& overlaps
+            ) const;
+
+            autoPtr<mapDistribute> calcLocalQueries
+            (
+                const pointField& centres,
+                const scalarField& radiusSqr,
+
+                pointField& allCentres,
+                scalarField& allRadiusSqr,
+                labelList& allSegmentMap
+            ) const;
+
+
+        // Surface redistribution
+
+            //- Calculate surface bounding box
+            void calcBounds(boundBox& bb, label& nPoints) const;
+
+            //- Does any part of triangle overlap bb.
+            static bool overlaps
+            (
+                const List<treeBoundBox>& bb,
+                const point& p0,
+                const point& p1,
+                const point& p2
+            );
+
+            //- Find points used in subset
+            static void subsetMeshMap
+            (
+                const triSurface& s,
+                const boolList& include,
+                const label nIncluded,
+                labelList& newToOldPoints,
+                labelList& oldToNewPoints,
+                labelList& newToOldFaces
+            );
+
+            //- Construct subsetted surface
+            static triSurface subsetMesh
+            (
+                const triSurface& s,
+                const labelList& newToOldPoints,
+                const labelList& oldToNewPoints,
+                const labelList& newToOldFaces
+            );
+
+            //- Subset given marked faces
+            static triSurface subsetMesh
+            (
+                const triSurface& s,
+                const boolList& include,
+                labelList& newToOldPoints,
+                labelList& newToOldFaces
+            );
+
+            //- Subset given marked faces
+            static triSurface subsetMesh
+            (
+                const triSurface& s,
+                const labelList& newToOldFaces,
+                labelList& newToOldPoints
+            );
+
+            //- Find triangle otherF in allFaces.
+            static label findTriangle
+            (
+                const List<labelledTri>& allFaces,
+                const labelListList& allPointFaces,
+                const labelledTri& otherF
+            );
+
+            //- Merge triSurface (subTris, subPoints) into allTris, allPoints.
+            static void merge
+            (
+                const scalar mergeDist,
+                const List<labelledTri>& subTris,
+                const pointField& subPoints,
+
+                List<labelledTri>& allTris,
+                pointField& allPoints,
+
+                labelList& faceConstructMap,
+                labelList& pointConstructMap
+            );
+
+            //- Distribute stored fields
+            template<class Type>
+            void distributeFields(const mapDistribute& map);
+
+
+        //- Disallow default bitwise copy construct
+        distributedTriSurfaceMesh(const distributedTriSurfaceMesh&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const distributedTriSurfaceMesh&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("distributedTriSurfaceMesh");
+
+
+    // Constructors
+
+        //- Construct from triSurface
+        distributedTriSurfaceMesh
+        (
+            const IOobject&,
+            const triSurface&,
+            const dictionary& dict
+        );
+
+        //- Construct read
+        distributedTriSurfaceMesh(const IOobject& io);
+
+        //- Construct from dictionary (used by searchableSurface)
+        distributedTriSurfaceMesh
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~distributedTriSurfaceMesh();
+
+        //- Clear storage
+        void clearOut();
+
+
+    // Member Functions
+
+        //- Return merge distance
+        static scalar mergeDist()
+        {
+            return mergeDist_;
+        }
+
+        //- Set the merge distance, returning the previous value
+        static scalar setMergeDist(const scalar t)
+        {
+            if (t < -VSMALL)
+            {
+                FatalErrorIn
+                (
+                    "scalar distributedTriSurfaceMesh::setMergeDist"
+                    "(const scalar t)"
+                )   << "Negative merge distance tolerance. This is not allowed."
+                    << abort(FatalError);
+            }
+
+            scalar oldTol = mergeDist_;
+            mergeDist_ = t;
+
+            return oldTol;
+        }
+
+
+        //- Triangle indexing (demand driven)
+        const globalIndex& globalTris() const;
+
+
+        // searchableSurface implementation
+
+            //- Whether supports volume type below. I.e. whether is closed.
+            //  Not supported.
+            virtual bool hasVolumeType() const
+            {
+                return false;
+            }
+
+            //- Range of global indices that can be returned.
+            virtual label globalSize() const
+            {
+                const labelList& offsets = globalTris().offsets();
+                return offsets[offsets.size()-1];
+            }
+
+            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;
+
+            //- Set bounds of surface. Bounds currently set as list of
+            //  bounding boxes. Will do redistribution of surface to locally
+            //  have all triangles overlapping bounds.
+            //  Larger bounds: more triangles (memory), more fully local tests
+            //  (quick).
+            //  keepNonLocal = true : keep triangles that do not overlap
+            //  any processor bounds.
+            //  Should really be split into a routine to determine decomposition
+            //  and one that does actual distribution but determining
+            //  decomposition with duplicate triangle merging requires
+            //  same amoun as work as actual distribution.
+            virtual void distribute
+            (
+                const List<treeBoundBox>&,
+                const bool keepNonLocal,
+                autoPtr<mapDistribute>& faceMap,
+                autoPtr<mapDistribute>& pointMap
+            );
+
+
+        // Other
+
+            //- Specific to triSurfaceMesh: from a set of hits (points and
+            //  indices) get the specified field. Misses do not get set.
+            virtual void getField
+            (
+                const word& fieldName,
+                const List<pointIndexHit>&,
+                labelList& values
+            ) const;
+
+            //- Subset the part of surface that is overlapping bounds.
+            static triSurface overlappingSurface
+            (
+                const triSurface&,
+                const List<treeBoundBox>&,
+                labelList& subPointMap,
+                labelList& subFaceMap
+            );
+
+            //- Print some stats. Parallel aware version of
+            //  triSurface::writeStats.
+            void writeStats(Ostream& os) const;
+
+
+        // regIOobject implementation
+
+            //- Write using given format, version and compression
+            virtual bool writeObject
+            (
+                IOstream::streamFormat fmt,
+                IOstream::versionNumber ver,
+                IOstream::compressionType cmp
+            ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "distributedTriSurfaceMeshTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/distributedTriSurfaceMeshTemplates.C b/src/meshTools/searchableSurface/distributedTriSurfaceMeshTemplates.C
new file mode 100644
index 0000000000..689e1a3f22
--- /dev/null
+++ b/src/meshTools/searchableSurface/distributedTriSurfaceMeshTemplates.C
@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "distributedTriSurfaceMesh.H"
+#include "triSurfaceFields.H"
+#include "mapDistribute.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//template<class Type>
+//void Foam::distributedTriSurfaceMesh::getField
+//(
+//    const word& fieldName,
+//    const List<pointIndexHit>& info,
+//    List<Type>& values
+//) const
+//{
+//    typedef DimensionedField<Type, triSurfaceGeoMesh> DimensionedSurfField;
+//
+//
+//    // Get query data (= local index of triangle)
+//    // ~~~~~~~~~~~~~~
+//
+//    labelList triangleIndex(info.size());
+//    autoPtr<mapDistribute> mapPtr
+//    (
+//        calcLocalQueries
+//        (
+//            info,
+//            triangleIndex
+//        )
+//    );
+//    const mapDistribute& map = mapPtr();
+//
+//
+//    // Do my tests
+//    // ~~~~~~~~~~~
+//
+//    const DimensionedSurfField& fld = lookupObject<DimensionedSurfField>
+//    (
+//        fieldName
+//    );
+//    const triSurface& s = static_cast<const triSurface&>(*this);
+//
+//    values.setSize(triangleIndex.size());
+//
+//    forAll(triangleIndex, i)
+//    {
+//        label triI = triangleIndex[i];
+//        values[i] = fld[triI];
+//    }
+//
+//
+//    // Send back results
+//    // ~~~~~~~~~~~~~~~~~
+//
+//    map.distribute
+//    (
+//        Pstream::nonBlocking,
+//        List<labelPair>(0),
+//        info.size(),
+//        map.constructMap(),     // what to send
+//        map.subMap(),           // what to receive
+//        values
+//    );
+//}
+
+
+template<class Type>
+void Foam::distributedTriSurfaceMesh::distributeFields
+(
+    const mapDistribute& map
+)
+{
+    typedef DimensionedField<Type, triSurfaceGeoMesh> DimensionedSurfField;
+
+    HashTable<const DimensionedSurfField*> fields
+    (
+        objectRegistry::lookupClass
+        <DimensionedSurfField >()
+    );
+
+    for
+    (
+        typename HashTable<const DimensionedSurfField*>::iterator fieldIter =
+            fields.begin();
+        fieldIter != fields.end();
+        ++fieldIter
+    )
+    {
+        DimensionedSurfField& field =
+            const_cast<DimensionedSurfField&>(*fieldIter());
+
+        label oldSize = field.size();
+
+        map.distribute
+        (
+            Pstream::nonBlocking,
+            List<labelPair>(0),
+            map.constructSize(),
+            map.subMap(),
+            map.constructMap(),
+            field
+        );
+
+        if (debug)
+        {
+            Info<< "Mapped " << field.typeName << ' ' << field.name()
+                << " from size " << oldSize << " to size " << field.size()
+                << endl;
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchableBox.C b/src/meshTools/searchableSurface/searchableBox.C
index 1858905d49..34e3faec99 100644
--- a/src/meshTools/searchableSurface/searchableBox.C
+++ b/src/meshTools/searchableSurface/searchableBox.C
@@ -484,7 +484,6 @@ void Foam::searchableBox::findLineAll
 
             hits.shrink();
             info[pointI].transfer(hits);
-            hits.clear();
         }
         else
         {
diff --git a/src/meshTools/searchableSurface/searchableBox.H b/src/meshTools/searchableSurface/searchableBox.H
index 1c51311511..a593d85065 100644
--- a/src/meshTools/searchableSurface/searchableBox.H
+++ b/src/meshTools/searchableSurface/searchableBox.H
@@ -121,6 +121,12 @@ public:
             return true;
         }
 
+        //- Range of local indices that can be returned.
+        virtual label size() const
+        {
+            return 6;
+        }
+
 
         // Single point queries.
 
diff --git a/src/meshTools/searchableSurface/searchablePlate.C b/src/meshTools/searchableSurface/searchablePlate.C
new file mode 100644
index 0000000000..9e8d5c061e
--- /dev/null
+++ b/src/meshTools/searchableSurface/searchablePlate.C
@@ -0,0 +1,366 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "searchablePlate.H"
+#include "addToRunTimeSelectionTable.H"
+#include "SortableList.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(searchablePlate, 0);
+addToRunTimeSelectionTable(searchableSurface, searchablePlate, dict);
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::direction Foam::searchablePlate::calcNormal(const point& span)
+{
+    direction normalDir = 3;
+
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        if (span[dir] < 0)
+        {
+            FatalErrorIn("searchablePlate::calcNormal()")
+                << "Span should have two positive and one zero entry. Now:"
+                << span << exit(FatalError);
+        }
+        else if (span[dir] < VSMALL)
+        {
+            if (normalDir == 3)
+            {
+                normalDir = dir;
+            }
+            else
+            {
+                // Multiple zero entries. Flag and exit.
+                normalDir = 3;
+                break;
+            }
+        }
+    }
+
+    if (normalDir == 3)
+    {
+        FatalErrorIn("searchablePlate::calcNormal()")
+            << "Span should have one and only zero entry. Now:" << span
+            << exit(FatalError);
+    }
+
+    return normalDir;
+}
+
+
+// Returns miss or hit with face (always 0)
+Foam::pointIndexHit Foam::searchablePlate::findNearest
+(
+    const point& sample,
+    const scalar nearestDistSqr
+) const
+{
+    // For every component direction can be
+    // left of min, right of max or inbetween.
+    // - outside points: project first one x plane (either min().x()
+    // or max().x()), then onto y plane and finally z. You should be left
+    // with intersection point
+    // - inside point: find nearest side (compare to mid point). Project onto
+    //   that.
+
+    // Project point on plane.
+    pointIndexHit info(true, sample, 0);
+    info.rawPoint()[normalDir_] = origin_[normalDir_];
+
+    // Clip to edges if outside
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        if (dir != normalDir_)
+        {
+            if (info.rawPoint()[dir] < origin_[dir])
+            {
+                info.rawPoint()[dir] = origin_[dir];
+            }
+            else if (info.rawPoint()[dir] > origin_[dir]+span_[dir])
+            {
+                info.rawPoint()[dir] = origin_[dir]+span_[dir];
+            }
+        }
+    }
+
+    // Check if outside. Optimisation: could do some checks on distance already
+    // on components above
+    if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
+    {
+        info.setMiss();
+        info.setIndex(-1);
+    }
+
+    return info;
+}
+
+
+Foam::pointIndexHit Foam::searchablePlate::findLine
+(
+    const point& start,
+    const point& end
+) const
+{
+    pointIndexHit info
+    (
+        true,
+        vector::zero,
+        0
+    );
+
+    const vector dir(end-start);
+
+    if (mag(dir[normalDir_]) < VSMALL)
+    {
+        info.setMiss();
+        info.setIndex(-1);
+    }
+    else
+    {
+        scalar t = (origin_[normalDir_]-start[normalDir_]) / dir[normalDir_];
+
+        if (t < 0 || t > 1)
+        {
+            info.setMiss();
+            info.setIndex(-1);
+        }
+        else
+        {
+            info.rawPoint() = start+t*dir;
+            info.rawPoint()[normalDir_] = origin_[normalDir_];
+
+            // Clip to edges
+            for (direction dir = 0; dir < vector::nComponents; dir++)
+            {
+                if (dir != normalDir_)
+                {
+                    if (info.rawPoint()[dir] < origin_[dir])
+                    {
+                        info.setMiss();
+                        info.setIndex(-1);
+                        break;
+                    }
+                    else if (info.rawPoint()[dir] > origin_[dir]+span_[dir])
+                    {
+                        info.setMiss();
+                        info.setIndex(-1);
+                        break;
+                    }
+                }
+            }
+        }
+    }
+
+    // Debug
+    if (info.hit())
+    {
+        treeBoundBox bb(origin_, origin_+span_);
+        bb.min()[normalDir_] -= 1E-6;
+        bb.max()[normalDir_] += 1E-6;
+
+        if (!bb.contains(info.hitPoint()))
+        {
+            FatalErrorIn("searchablePlate::findLine(..)")
+                << "bb:" << bb << endl
+                << "origin_:" << origin_ << endl
+                << "span_:" << span_ << endl
+                << "normalDir_:" << normalDir_ << endl
+                << "hitPoint:" << info.hitPoint()
+                << abort(FatalError);
+        }
+    }
+
+    return info;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchablePlate::searchablePlate
+(
+    const IOobject& io,
+    const point& origin,
+    const vector& span
+)
+:
+    searchableSurface(io),
+    origin_(origin),
+    span_(span),
+    normalDir_(calcNormal(span_))
+{}
+
+
+Foam::searchablePlate::searchablePlate
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    searchableSurface(io),
+    origin_(dict.lookup("origin")),
+    span_(dict.lookup("span")),
+    normalDir_(calcNormal(span_))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchablePlate::~searchablePlate()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::wordList& Foam::searchablePlate::regions() const
+{
+    if (regions_.size() == 0)
+    {
+        regions_.setSize(1);
+        regions_[0] = "region0";
+    }
+    return regions_;
+}
+
+
+void Foam::searchablePlate::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::searchablePlate::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::searchablePlate::findLineAny
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    findLine(start, end, info);
+}
+
+
+void Foam::searchablePlate::findLineAll
+(
+    const pointField& start,
+    const pointField& end,
+    List<List<pointIndexHit> >& info
+) const
+{
+    List<pointIndexHit> nearestInfo;
+    findLine(start, end, nearestInfo);
+
+    info.setSize(start.size());
+    forAll(info, pointI)
+    {
+        if (nearestInfo[pointI].hit())
+        {
+            info[pointI].setSize(1);
+            info[pointI][0] = nearestInfo[pointI];
+        }
+        else
+        {
+            info[pointI].clear();
+        }
+    }
+}
+
+
+void Foam::searchablePlate::getRegion
+(
+    const List<pointIndexHit>& info,
+    labelList& region
+) const
+{
+    region.setSize(info.size());
+    region = 0;
+}
+
+
+void Foam::searchablePlate::getNormal
+(
+    const List<pointIndexHit>& info,
+    vectorField& normal
+) const
+{
+    normal.setSize(info.size());
+    normal = vector::zero;
+    forAll(normal, i)
+    {
+        normal[i][normalDir_] = 1.0;
+    }
+}
+
+
+void Foam::searchablePlate::getVolumeType
+(
+    const pointField& points,
+    List<volumeType>& volType
+) const
+{
+    FatalErrorIn
+    (
+        "searchableCollection::getVolumeType(const pointField&"
+        ", List<volumeType>&) const"
+    )   << "Volume type not supported for plate."
+        << exit(FatalError);
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchablePlate.H b/src/meshTools/searchableSurface/searchablePlate.H
new file mode 100644
index 0000000000..0ad916fe79
--- /dev/null
+++ b/src/meshTools/searchableSurface/searchablePlate.H
@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::searchablePlate
+
+Description
+    Searching on finite plate. Plate has to be aligned with coordinate
+    axes!
+
+SourceFiles
+    searchablePlate.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchablePlate_H
+#define searchablePlate_H
+
+#include "searchableSurface.H"
+#include "treeBoundBox.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+/*---------------------------------------------------------------------------*\
+                           Class searchablePlate Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchablePlate
+:
+    public searchableSurface
+{
+private:
+
+    // Private Member Data
+
+        const point origin_;
+
+        const vector span_;
+
+        //- Coordinate direction which is normal
+        const direction normalDir_;
+
+        mutable wordList regions_;
+
+
+    // Private Member Functions
+
+        //- Calculate normal direction from span
+        static direction calcNormal(const point&);
+
+        pointIndexHit findNearest
+        (
+            const point& sample,
+            const scalar nearestDistSqr
+        ) const;
+
+        pointIndexHit findLine
+        (
+            const point& start,
+            const point& end
+        ) const;
+
+
+        //- Disallow default bitwise copy construct
+        searchablePlate(const searchablePlate&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const searchablePlate&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("searchablePlate");
+
+
+    // Constructors
+
+        //- Construct from components
+        searchablePlate
+        (
+            const IOobject& io,
+            const point& origin,
+            const point& span
+        );
+
+        //- Construct from dictionary (used by searchableSurface)
+        searchablePlate
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+    // Destructor
+
+        virtual ~searchablePlate();
+
+
+    // Member Functions
+
+        virtual const wordList& regions() const;
+
+        //- Whether supports volume type below
+        virtual bool hasVolumeType() const
+        {
+            return false;
+        }
+
+        //- Range of local indices that can be returned.
+        virtual label size() const
+        {
+            return 1;
+        }
+
+
+        // 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("searchablePlate::writeData(Ostream&) const");
+                return false;
+            }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchableSphere.H b/src/meshTools/searchableSurface/searchableSphere.H
index 61b114de4e..ecafdbf2e5 100644
--- a/src/meshTools/searchableSurface/searchableSphere.H
+++ b/src/meshTools/searchableSurface/searchableSphere.H
@@ -36,8 +36,8 @@ SourceFiles
 #ifndef searchableSphere_H
 #define searchableSphere_H
 
-#include "searchableSurface.H"
 #include "treeBoundBox.H"
+#include "searchableSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -127,6 +127,12 @@ public:
             return true;
         }
 
+        //- Range of local indices that can be returned.
+        virtual label size() const
+        {
+            return 1;
+        }
+
 
         // Multiple point queries.
 
diff --git a/src/meshTools/searchableSurface/searchableSurface.H b/src/meshTools/searchableSurface/searchableSurface.H
index f9dc8b1289..6c788448d6 100644
--- a/src/meshTools/searchableSurface/searchableSurface.H
+++ b/src/meshTools/searchableSurface/searchableSurface.H
@@ -56,6 +56,7 @@ namespace Foam
 
 // Forward declaration of classes
 class objectRegistry;
+class mapDistribute;
 class treeBoundBox;
 
 /*---------------------------------------------------------------------------*\
@@ -180,6 +181,15 @@ public:
         //- Whether supports volume type below.
         virtual bool hasVolumeType() const = 0;
 
+        //- Range of local indices that can be returned.
+        virtual label size() const = 0;
+
+        //- Range of global indices that can be returned.
+        virtual label globalSize() const
+        {
+            return size();
+        }
+
 
         // Single point queries.
 
@@ -296,16 +306,18 @@ public:
 
         // Other
 
-            ////- Get bounding box.
-            //const boundBox& bounds() const = 0;
-
-            ////- Set bounding box.
-            //void setBounds
-            //(
-            //    const boundBox&,
-            //    autoPtr<mapDistribute>& faceMap,
-            //    autoPtr<mapDistribute>& pointMap
-            //) = 0;
+            //- Set bounds of surface. Bounds currently set as list of
+            //  bounding boxes. The bounds are hints to the surface as for
+            //  the range of queries it can expect. faceMap/pointMap can be
+            //  set if the surface has done any redistribution.
+            virtual void distribute
+            (
+                const List<treeBoundBox>&,
+                const bool keepNonLocal,
+                autoPtr<mapDistribute>& faceMap,
+                autoPtr<mapDistribute>& pointMap
+            )
+            {}
 
 
 };
diff --git a/src/meshTools/searchableSurface/searchableSurfaceWithGaps.C b/src/meshTools/searchableSurface/searchableSurfaceWithGaps.C
new file mode 100644
index 0000000000..f138daa7a4
--- /dev/null
+++ b/src/meshTools/searchableSurface/searchableSurfaceWithGaps.C
@@ -0,0 +1,321 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "searchableSurfaceWithGaps.H"
+#include "addToRunTimeSelectionTable.H"
+#include "SortableList.H"
+#include "Time.H"
+#include "ListOps.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(searchableSurfaceWithGaps, 0);
+addToRunTimeSelectionTable(searchableSurface, searchableSurfaceWithGaps, dict);
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::Pair<Foam::vector> Foam::searchableSurfaceWithGaps::offsetVecs
+(
+    const point& start,
+    const point& end
+) const
+{
+    Pair<vector> offsets(vector::zero, vector::zero);
+
+    vector n(end-start);
+
+    scalar magN = mag(n);
+
+    if (magN > SMALL)
+    {
+        n /= magN;
+
+        // Do first offset vector. Is the coordinate axes with the smallest
+        // component along the vector n.
+        scalar minMag = GREAT;
+        direction minCmpt = 0;
+
+        for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
+        {
+            if (mag(n[cmpt]) < minMag)
+            {
+                minMag = mag(n[cmpt]);
+                minCmpt = cmpt;
+            }
+        }
+
+        offsets[0][minCmpt] = 1.0;
+        // Orthogonalise
+        offsets[0] -= n[minCmpt]*n;
+        // Scale
+        offsets[0] *= gap_/mag(offsets[0]);
+
+
+        // Do second offset vector perp to original edge and first offset vector
+        offsets[1] = n ^ offsets[0];
+        offsets[1] *= gap_/mag(offsets[1]);
+    }
+
+    return offsets;    
+}
+
+
+void Foam::searchableSurfaceWithGaps::offsetVecs
+(
+    const pointField& start,
+    const pointField& end,
+    pointField& offset0,
+    pointField& offset1
+) const
+{
+    offset0.setSize(start.size());
+    offset1.setSize(start.size());
+
+    forAll(start, i)
+    {
+        const Pair<vector> offsets(offsetVecs(start[i], end[i]));
+        offset0[i] = offsets[0];
+        offset1[i] = offsets[1];
+    }
+}
+
+
+Foam::label Foam::searchableSurfaceWithGaps::countMisses
+(
+    const List<pointIndexHit>& info,
+    labelList& missMap
+)
+{
+    label nMiss = 0;
+    forAll(info, i)
+    {
+        if (!info[i].hit())
+        {
+            nMiss++;
+        }
+    }
+
+    missMap.setSize(nMiss);
+    nMiss = 0;
+
+    forAll(info, i)
+    {
+        if (!info[i].hit())
+        {
+            missMap[nMiss++] = i;
+        }
+    }
+
+    return nMiss;
+}
+
+
+// Anything not a hit in both counts as a hit
+Foam::label Foam::searchableSurfaceWithGaps::countMisses
+(
+    const List<pointIndexHit>& plusInfo,
+    const List<pointIndexHit>& minInfo,
+    labelList& missMap
+)
+{
+    label nMiss = 0;
+    forAll(plusInfo, i)
+    {
+        if (!plusInfo[i].hit() || !minInfo[i].hit())
+        {
+            nMiss++;
+        }
+    }
+
+    missMap.setSize(nMiss);
+    nMiss = 0;
+
+    forAll(plusInfo, i)
+    {
+        if (!plusInfo[i].hit() || !minInfo[i].hit())
+        {
+            missMap[nMiss++] = i;
+        }
+    }
+
+    return nMiss;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceWithGaps::searchableSurfaceWithGaps
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    searchableSurface(io),
+    gap_(readScalar(dict.lookup("gap"))),
+    subGeom_(1)
+{
+    const word subGeomName(dict.lookup("surface"));
+
+    const searchableSurface& s =
+        io.db().lookupObject<searchableSurface>(subGeomName);
+
+    subGeom_.set(0, &const_cast<searchableSurface&>(s));
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceWithGaps::~searchableSurfaceWithGaps()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::searchableSurfaceWithGaps::findLine
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    surface().findLine(start, end, info);
+
+    // Count number of misses. Determine map
+    labelList compactMap;
+    label nMiss = countMisses(info, compactMap);
+
+    if (returnReduce(nMiss, sumOp<label>()) > 0)
+    {
+        pointField compactStart(start, compactMap);
+        pointField compactEnd(end, compactMap);
+
+        // Calculate offset vector
+        pointField offset0, offset1;
+        offsetVecs
+        (
+            compactStart,
+            compactEnd,
+            offset0,
+            offset1
+        );
+
+        // test in pairs: only if both perturbations hit something
+        // do we accept the hit.
+
+        List<pointIndexHit> plusInfo;
+        surface().findLine(start+offset0, end+offset0, plusInfo);
+        List<pointIndexHit> minInfo;
+        surface().findLine(start-offset0, end-offset0, minInfo);
+
+        // Extract any hits
+        forAll(plusInfo, i)
+        {
+            if (plusInfo[i].hit() && minInfo[i].hit())
+            {
+                info[compactMap[i]] = plusInfo[i].hitPoint()-offset0[i];
+            }
+        }
+
+        labelList plusMissMap;
+        nMiss = countMisses(plusInfo, minInfo, plusMissMap);
+
+        if (returnReduce(nMiss, sumOp<label>()) > 0)
+        {
+            // Extract (inplace possible because of order)
+            forAll(plusMissMap, i)
+            {
+                label mapI = plusMissMap[i];
+                compactStart[i] = compactStart[mapI];
+                compactEnd[i] = compactEnd[mapI];
+                compactMap[i] = compactMap[mapI];
+            }
+            compactStart.setSize(plusMissMap.size());
+            compactEnd.setSize(plusMissMap.size());
+            compactMap.setSize(plusMissMap.size());
+
+
+            surface().findLine(start+offset1, end+offset1, plusInfo);
+            surface().findLine(start-offset1, end-offset1, minInfo);
+
+            // Extract any hits
+            forAll(plusInfo, i)
+            {
+                if (plusInfo[i].hit() && minInfo[i].hit())
+                {
+                    info[compactMap[i]] = plusInfo[i].hitPoint()-offset1[i];
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::searchableSurfaceWithGaps::findLineAny
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    // To be done ...
+    findLine(start, end, info);
+}
+
+
+void Foam::searchableSurfaceWithGaps::findLineAll
+(
+    const pointField& start,
+    const pointField& end,
+    List<List<pointIndexHit> >& info
+) const
+{
+    // To be done. Assume for now only one intersection.
+    List<pointIndexHit> nearestInfo;
+    findLine(start, end, nearestInfo);
+
+    info.setSize(start.size());
+    forAll(info, pointI)
+    {
+        if (nearestInfo[pointI].hit())
+        {
+            info[pointI].setSize(1);
+            info[pointI][0] = nearestInfo[pointI];
+        }
+        else
+        {
+            info[pointI].clear();
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchableSurfaceWithGaps.H b/src/meshTools/searchableSurface/searchableSurfaceWithGaps.H
new file mode 100644
index 0000000000..1bdddb9706
--- /dev/null
+++ b/src/meshTools/searchableSurface/searchableSurfaceWithGaps.H
@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::searchableSurfaceWithGaps
+
+Description
+    searchableSurface using multiple slightly shifted underlying surfaces
+    to make sure pierces don't go through gaps. 
+
+SourceFiles
+    searchableSurfaceWithGaps.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableSurfaceWithGaps_H
+#define searchableSurfaceWithGaps_H
+
+#include "searchableSurface.H"
+#include "UPtrList.H"
+#include "Pair.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+/*---------------------------------------------------------------------------*\
+                           Class searchableSurfaceWithGaps Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableSurfaceWithGaps
+:
+    public searchableSurface
+{
+private:
+
+    // Private Member Data
+
+        //- gap size in meter
+        const scalar gap_;
+
+        //- Underlying geometry (size 1)
+        UPtrList<searchableSurface> subGeom_;
+
+
+    // Private Member Functions
+
+        Pair<vector> offsetVecs(const point&, const point&) const;
+
+        void offsetVecs
+        (
+            const pointField& start,
+            const pointField& end,
+            pointField& offset0,
+            pointField& offset1
+        ) const;
+
+        static label countMisses
+        (
+            const List<pointIndexHit>& info,
+            labelList& missMap
+        );
+
+        static label countMisses
+        (
+            const List<pointIndexHit>& plusInfo,
+            const List<pointIndexHit>& minInfo,
+            labelList& missMap
+        );
+
+
+        //- Disallow default bitwise copy construct
+        searchableSurfaceWithGaps(const searchableSurfaceWithGaps&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const searchableSurfaceWithGaps&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("searchableSurfaceWithGaps");
+
+
+    // Constructors
+
+        //- Construct from dictionary (used by searchableSurface)
+        searchableSurfaceWithGaps
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+    // Destructor
+
+        virtual ~searchableSurfaceWithGaps();
+
+
+    // Member Functions
+
+        const searchableSurface& surface() const
+        {
+            return subGeom_[0];
+        }
+
+
+        virtual const wordList& regions() const
+        {
+            return surface().regions();
+        }
+
+        //- Whether supports volume type below
+        virtual bool hasVolumeType() const
+        {
+            return surface().hasVolumeType();
+        }
+
+        //- Range of local indices that can be returned.
+        virtual label size() const
+        {
+            return surface().size();
+        }
+
+
+        // Multiple point queries.
+
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>& info
+            ) const
+            {
+                surface().findNearest
+                (
+                    sample,
+                    nearestDistSqr,
+                    info
+                );
+            }
+
+            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>& info,
+                labelList& region
+            ) const
+            {
+                surface().getRegion(info, region);
+            }
+
+            //- From a set of points and indices get the normal
+            virtual void getNormal
+            (
+                const List<pointIndexHit>& info,
+                vectorField& normal
+            ) const
+            {
+                surface().getNormal(info, normal);
+            }
+
+            //- Determine type (inside/outside/mixed) for point. unknown if
+            //  cannot be determined (e.g. non-manifold surface)
+            virtual void getVolumeType
+            (
+                const pointField& samples,
+                List<volumeType>& info
+            ) const
+            {
+                surface().getVolumeType(samples, info);
+            }
+
+
+        // regIOobject implementation
+
+            bool writeData(Ostream& os) const
+            {
+                return surface().writeData(os);
+            }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchableSurfaces.H b/src/meshTools/searchableSurface/searchableSurfaces.H
index 4e6c1bc2ed..0735c3d9d7 100644
--- a/src/meshTools/searchableSurface/searchableSurfaces.H
+++ b/src/meshTools/searchableSurface/searchableSurfaces.H
@@ -149,7 +149,7 @@ public:
                 const pointField& start,
                 const pointField& end,
                 labelListList& surfaces,
-                List<List<pointIndexHit> >& surfaceHits
+                List<List<pointIndexHit> >&
             ) const;
 
             //- Find nearest. Return -1 (and a miss()) or surface and nearest
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C
index a426741960..a51b3ff495 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.C
@@ -28,6 +28,8 @@ License
 #include "Random.H"
 #include "addToRunTimeSelectionTable.H"
 #include "EdgeMap.H"
+#include "triSurfaceFields.H"
+#include "Time.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -299,57 +301,6 @@ const Foam::wordList& Foam::triSurfaceMesh::regions() const
 }
 
 
-//Foam::pointIndexHit Foam::triSurfaceMesh::findNearest
-//(
-//    const point& sample,
-//    const scalar nearestDistSqr
-//) const
-//{
-//    return tree().findNearest(sample, nearestDistSqr);
-//}
-//
-//
-//Foam::pointIndexHit Foam::triSurfaceMesh::findNearestOnEdge
-//(
-//    const point& sample,
-//    const scalar nearestDistSqr
-//) const
-//{
-//    return = edgeTree().findNearest(sample, nearestDistSqr);
-//}
-//
-//
-//Foam::pointIndexHit Foam::triSurfaceMesh::findNearest
-//(
-//    const linePointRef& ln,
-//    treeBoundBox& tightest,
-//    point& linePoint
-//) const
-//{
-//    return tree().findNearest(ln, tightest, linePoint);
-//}
-//
-//
-//Foam::pointIndexHit Foam::triSurfaceMesh::findLine
-//(
-//    const point& start,
-//    const point& end
-//) const
-//{
-//    return tree().findLine(start, end);
-//}
-//
-//
-//Foam::pointIndexHit Foam::triSurfaceMesh::findLineAny
-//(
-//    const point& start,
-//    const point& end
-//) const
-//{
-//    return tree().findLineAny(start, end);
-//}
-
-
 // Find out if surface is closed.
 bool Foam::triSurfaceMesh::hasVolumeType() const
 {
@@ -491,7 +442,6 @@ void Foam::triSurfaceMesh::findLineAll
 
             hits.shrink();
             info[pointI].transfer(hits);
-            hits.clear();
         }
         else
         {
@@ -545,6 +495,29 @@ void Foam::triSurfaceMesh::getNormal
 }
 
 
+void Foam::triSurfaceMesh::getField
+(
+    const word& fieldName,
+    const List<pointIndexHit>& info,
+    labelList& values
+) const
+{
+    const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
+    (
+        fieldName
+    );
+
+    values.setSize(info.size());
+    forAll(info, i)
+    {
+        if (info[i].hit())
+        {
+            values[i] = fld[info[i].index()];
+        }
+    }
+}
+
+
 void Foam::triSurfaceMesh::getVolumeType
 (
     const pointField& points,
@@ -564,4 +537,31 @@ void Foam::triSurfaceMesh::getVolumeType
 }
 
 
+//- Write using given format, version and compression
+bool Foam::triSurfaceMesh::writeObject
+(
+    IOstream::streamFormat fmt,
+    IOstream::versionNumber ver,
+    IOstream::compressionType cmp
+) const
+{
+    fileName fullPath(searchableSurface::objectPath());
+
+    if (!mkDir(fullPath.path()))
+    {
+        return false;
+    }
+
+    triSurface::write(fullPath);
+
+    if (!file(fullPath))
+    {
+        return false;
+    }
+
+    //return objectRegistry::writeObject(fmt, ver, cmp);
+    return true;
+}
+
+
 // ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.H b/src/meshTools/searchableSurface/triSurfaceMesh.H
index 8a0709386c..12582c4e83 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.H
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.H
@@ -36,6 +36,7 @@ SourceFiles
 #ifndef triSurfaceMesh_H
 #define triSurfaceMesh_H
 
+#include "treeBoundBox.H"
 #include "searchableSurface.H"
 #include "objectRegistry.H"
 #include "indexedOctree.H"
@@ -142,8 +143,11 @@ public:
             //- Whether supports volume type below. I.e. whether is closed.
             virtual bool hasVolumeType() const;
 
-
-        // Multiple point queries.
+            //- Range of local indices that can be returned.
+            virtual label size() const
+            {
+                return triSurface::size();
+            }
 
             virtual void findNearest
             (
@@ -196,6 +200,30 @@ public:
                 List<volumeType>&
             ) const;
 
+            //- Set bounds of surface. Bounds currently set as list of
+            //  bounding boxes. The bounds are hints to the surface as for
+            //  the range of queries it can expect. faceMap/pointMap can be
+            //  set if the surface has done any redistribution.
+            virtual void distribute
+            (
+                const List<treeBoundBox>&,
+                const bool keepNonLocal,
+                autoPtr<mapDistribute>& faceMap,
+                autoPtr<mapDistribute>& pointMap
+            )
+            {}
+
+        // Other
+
+            //- Specific to triSurfaceMesh: from a set of hits (points and
+            //  indices) get the specified field. Misses do not get set.
+            virtual void getField
+            (
+                const word& fieldName,
+                const List<pointIndexHit>&,
+                labelList& values
+            ) const;
+
 
         // regIOobject implementation
 
@@ -205,6 +233,14 @@ public:
                 return false;
             }
 
+            //- Write using given format, version and compression
+            virtual bool writeObject
+            (
+                IOstream::streamFormat fmt,
+                IOstream::versionNumber ver,
+                IOstream::compressionType cmp
+            ) const;
+
 };
 
 
diff --git a/src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.C b/src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.C
index 7f588b8335..577d004982 100644
--- a/src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.C
+++ b/src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.C
@@ -288,11 +288,9 @@ Foam::edgeSurface::edgeSurface
 
     // Transfer.
     allEdges.shrink();
-    allEdges.clear();
     edges_.transfer(allEdges);
 
     allParentEdges.shrink();
-    allParentEdges.clear();
     parentEdges_.transfer(allParentEdges);
 
     forAll(allFaceEdges, faceI)
@@ -300,7 +298,6 @@ Foam::edgeSurface::edgeSurface
         DynamicList<label>& allFEdges = allFaceEdges[faceI];
 
         allFEdges.shrink();
-        allFEdges.clear();
         faceEdges_[faceI].transfer(allFEdges);
     }
 
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionTemplates.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionTemplates.C
index 7b349102b2..fd06704ab1 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionTemplates.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionTemplates.C
@@ -40,7 +40,6 @@ void Foam::surfaceIntersection::transfer
 {
     dList.shrink();
     lList.transfer(dList);
-    dList.clear();
 }
 
 
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
index 9f858b95cd..10662625f0 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
@@ -199,7 +199,6 @@ void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
     }
     featurePoints.shrink();
     featurePoints_.transfer(featurePoints);
-    featurePoints.clear();
 }
 
 
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H
index 5f4422c2d2..619632393b 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H
@@ -36,7 +36,6 @@ SourceFiles
 #ifndef triSurfaceSearch_H
 #define triSurfaceSearch_H
 
-#include "DynamicList.H"
 #include "pointField.H"
 #include "boolList.H"
 #include "pointIndexHit.H"
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
index 0083917df6..a8edd2d334 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
@@ -300,11 +300,7 @@ Foam::triSurface Foam::triSurfaceTools::doRefine
     allPoints.transfer(newPoints);
     newPoints.clear();
 
-    List<labelledTri> allFaces;
-    allFaces.transfer(newFaces);
-    newFaces.clear();
-
-    return triSurface(allFaces, surf.patches(), allPoints);
+    return triSurface(newFaces, surf.patches(), allPoints, true);
 }
 
 
@@ -1988,7 +1984,7 @@ Foam::triSurface Foam::triSurfaceTools::greenRefine
     newFaces.shrink();
     newPoints.setSize(newPointI);
 
-    return triSurface(newFaces, surf.patches(), newPoints);
+    return triSurface(newFaces, surf.patches(), newPoints, true);
 }
 
 
diff --git a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
index e8c9a29552..9fe49a1b1f 100644
--- a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
+++ b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
@@ -101,8 +101,6 @@ 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;
     label i2 = (i1 + 1) % 3;
@@ -113,6 +111,8 @@ bool Foam::triangleFuncs::intersectAxesBundle
     scalar u2 = V20[i1];
     scalar v2 = V20[i2];
 
+    scalar localScale = mag(u1)+mag(v1)+mag(u2)+mag(v2);
+
     scalar det = v2*u1 - u2*v1;
 
     // Fix for  V0:(-31.71428 0 -15.10714)
@@ -136,7 +136,7 @@ bool Foam::triangleFuncs::intersectAxesBundle
         scalar beta = 0;
         bool inter = false;
 
-        if (Foam::mag(u1)/localScale < SMALL)
+        if (Foam::mag(u1) < SMALL)
         {
             beta = u0/u2;
             if ((beta >= 0) && (beta <= 1))
diff --git a/src/sampling/cuttingPlane/cuttingPlane.C b/src/sampling/cuttingPlane/cuttingPlane.C
index 472d30e782..e1a0212ff0 100644
--- a/src/sampling/cuttingPlane/cuttingPlane.C
+++ b/src/sampling/cuttingPlane/cuttingPlane.C
@@ -154,7 +154,6 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
 
     dynCuttingPoints.shrink();
     cutPoints_.transfer(dynCuttingPoints);
-    dynCuttingPoints.clear();
 
     return edgePoint;
 }
diff --git a/src/sampling/sampledSet/midPoint/midPointSet.C b/src/sampling/sampledSet/midPoint/midPointSet.C
index 37d4d45fc7..e69a379bac 100644
--- a/src/sampling/sampledSet/midPoint/midPointSet.C
+++ b/src/sampling/sampledSet/midPoint/midPointSet.C
@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "midPointSet.H"
-#include "DynamicList.H"
 #include "polyMesh.H"
 #include "addToRunTimeSelectionTable.H"
 
diff --git a/src/sampling/sampledSet/midPointAndFace/midPointAndFaceSet.C b/src/sampling/sampledSet/midPointAndFace/midPointAndFaceSet.C
index b1399cc8ad..1d8a37176d 100644
--- a/src/sampling/sampledSet/midPointAndFace/midPointAndFaceSet.C
+++ b/src/sampling/sampledSet/midPointAndFace/midPointAndFaceSet.C
@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "midPointAndFaceSet.H"
-#include "DynamicList.H"
 #include "polyMesh.H"
 #include "addToRunTimeSelectionTable.H"
 
diff --git a/src/thermophysicalModels/combustion/chemistryReaders/foamChemistryReader/foamChemistryReader.H b/src/thermophysicalModels/combustion/chemistryReaders/foamChemistryReader/foamChemistryReader.H
index ac0d0651b1..6ec766546c 100644
--- a/src/thermophysicalModels/combustion/chemistryReaders/foamChemistryReader/foamChemistryReader.H
+++ b/src/thermophysicalModels/combustion/chemistryReaders/foamChemistryReader/foamChemistryReader.H
@@ -41,7 +41,6 @@ SourceFiles
 #include "typeInfo.H"
 #include "HashPtrTable.H"
 #include "SLPtrList.H"
-#include "DynamicList.H"
 #include "labelList.H"
 #include "speciesTable.H"
 #include "atomicWeights.H"
diff --git a/src/triSurface/tools/labelledTri/labelledTri.H b/src/triSurface/tools/labelledTri/labelledTri.H
index af61315d80..4470298e8e 100644
--- a/src/triSurface/tools/labelledTri/labelledTri.H
+++ b/src/triSurface/tools/labelledTri/labelledTri.H
@@ -105,6 +105,10 @@ public:
 };
 
 
+template<>
+inline bool contiguous<labelledTri>()  {return true;}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/triSurface/tools/labelledTri/labelledTriI.H b/src/triSurface/tools/labelledTri/labelledTriI.H
index 43195124ac..075b5b8205 100644
--- a/src/triSurface/tools/labelledTri/labelledTriI.H
+++ b/src/triSurface/tools/labelledTri/labelledTriI.H
@@ -58,16 +58,7 @@ inline labelledTri::labelledTri
 
 inline labelledTri::labelledTri(Istream& is)
 {
-    // Read beginning of labelledTri point pair
-    is.readBegin("labelledTri");
-
-    is >> static_cast<triFace&>(*this) >> region_;
-
-    // Read end of labelledTri point pair
-    is.readEnd("labelledTri");
-
-    // Check state of Istream
-    is.check("labelledTri::labelledTri(Istream& is)");
+    operator>>(is, *this);
 }
 
 
@@ -88,13 +79,20 @@ inline label& labelledTri::region()
 
 inline Istream& operator>>(Istream& is, labelledTri& t)
 {
-    // Read beginning of labelledTri point pair
-    is.readBegin("labelledTri");
+    if (is.format() == IOstream::ASCII)
+    {
+        // Read beginning of labelledTri point pair
+        is.readBegin("labelledTri");
 
-    is  >> static_cast<triFace&>(t) >> t.region_;
+        is  >> static_cast<triFace&>(t) >> t.region_;
 
-    // Read end of labelledTri point pair
-    is.readEnd("labelledTri");
+        // Read end of labelledTri point pair
+        is.readEnd("labelledTri");
+    }
+    else
+    {
+        is.read(reinterpret_cast<char*>(&t), sizeof(labelledTri));
+    }
 
     // Check state of Ostream
     is.check("Istream& operator>>(Istream&, labelledTri&)");
@@ -105,9 +103,24 @@ inline Istream& operator>>(Istream& is, labelledTri& t)
 
 inline Ostream& operator<<(Ostream& os, const labelledTri& t)
 {
-    os  << token::BEGIN_LIST
-        << static_cast<const triFace&>(t) << token::SPACE << t.region_
-        << token::END_LIST;
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << token::BEGIN_LIST
+            << static_cast<const triFace&>(t) << token::SPACE << t.region_
+            << token::END_LIST;
+    }
+    else
+    {
+        os.write
+        (
+            reinterpret_cast<const char*>(&t),
+            sizeof(labelledTri)
+        );
+    }
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const labelledTri&)");
+
 
     return os;
 }
diff --git a/src/triSurface/triSurface/interfaces/AC3D/readAC.C b/src/triSurface/triSurface/interfaces/AC3D/readAC.C
index 06afea1dba..85974d52e1 100644
--- a/src/triSurface/triSurface/interfaces/AC3D/readAC.C
+++ b/src/triSurface/triSurface/interfaces/AC3D/readAC.C
@@ -339,11 +339,7 @@ bool triSurface::readAC(const fileName& ACfileName)
     allPoints.transfer(points);
     points.clear();
 
-    List<labelledTri> allFaces;
-    allFaces.transfer(faces);
-    faces.clear();
-
-    *this = triSurface(allFaces, patches, allPoints);
+    *this = triSurface(faces, patches, allPoints, true);
 
     stitchTriangles(allPoints);
 
diff --git a/src/triSurface/triSurface/interfaces/NAS/readNAS.C b/src/triSurface/triSurface/interfaces/NAS/readNAS.C
index 1f4338ac13..2ec15c6706 100644
--- a/src/triSurface/triSurface/interfaces/NAS/readNAS.C
+++ b/src/triSurface/triSurface/interfaces/NAS/readNAS.C
@@ -350,7 +350,7 @@ bool triSurface::readNAS(const fileName& OBJfileName)
     points.clear();
 
     // Create triSurface
-    *this = triSurface(faces, patches, allPoints);
+    *this = triSurface(faces, patches, allPoints, true);
 
     return true;
 }
diff --git a/src/triSurface/triSurface/interfaces/OBJ/readOBJ.C b/src/triSurface/triSurface/interfaces/OBJ/readOBJ.C
index 018fd4932f..6b522c4f0e 100644
--- a/src/triSurface/triSurface/interfaces/OBJ/readOBJ.C
+++ b/src/triSurface/triSurface/interfaces/OBJ/readOBJ.C
@@ -195,14 +195,9 @@ bool triSurface::readOBJ(const fileName& OBJfileName)
     // Transfer DynamicLists to straight ones.
     pointField allPoints;
     allPoints.transfer(points);
-    points.clear();
-
-    List<labelledTri> allFaces;
-    allFaces.transfer(faces);
-    faces.clear();
 
     // Create triSurface
-    *this = triSurface(allFaces, patches, allPoints);
+    *this = triSurface(faces, patches, allPoints, true);
 
     return true;
 }
diff --git a/src/triSurface/triSurface/triSurface.C b/src/triSurface/triSurface/triSurface.C
index 89a48a30d6..a9f7b53e66 100644
--- a/src/triSurface/triSurface/triSurface.C
+++ b/src/triSurface/triSurface/triSurface.C
@@ -634,7 +634,8 @@ surfacePatchList triSurface::calcPatches(labelList& faceMap) const
     // Compact regions
 
     // Get last region
-    label maxRegion = 0;
+    //label maxRegion = 0;                  // for compacted regions
+    label maxRegion = patches_.size()-1;    // for non-compacted regions
 
     if (faceMap.size() > 0)
     {