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/meshTools/searchableSurface/distributedTriSurfaceMesh.C b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
new file mode 100644
index 0000000000..b0ac235797
--- /dev/null
+++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
@@ -0,0 +1,2217 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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<Pair<point> >& 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(Pair<point>(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(Pair<point>(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(Pair<point>(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(Pair<point>(clipPt0, clipPt1));
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+Foam::autoPtr<Foam::mapDistribute>
+Foam::distributedTriSurfaceMesh::constructSegments
+(
+    const pointField& start,
+    const pointField& end,
+
+    List<Pair<point> >& 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<Pair<point> > 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]);
+            dynSendMap[procI].clear();
+        }
+
+        allSegments.transfer(dynAllSegments.shrink());
+        dynAllSegments.clear();
+        allSegmentMap.transfer(dynAllSegmentMap.shrink());
+        dynAllSegmentMap.clear();
+    }
+
+
+    // 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<Pair<point> > 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]);
+            dynSendMap[procI].clear();
+        }
+
+        allCentres.transfer(dynAllCentres.shrink());
+        dynAllCentres.clear();
+        allRadiusSqr.transfer(dynAllRadiusSqr.shrink());
+        dynAllRadiusSqr.clear();
+        allSegmentMap.transfer(dynAllSegmentMap.shrink());
+        dynAllSegmentMap.clear();
+    }
+
+
+    // 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..3e7dd3b57c
--- /dev/null
+++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H
@@ -0,0 +1,453 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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;
+
+/*---------------------------------------------------------------------------*\
+                           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<Pair<point> >&,
+                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<Pair<point> >& 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;
+        }
+    }
+}
+
+
+// ************************************************************************* //