diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index 7e9274616e..9e35ee1edc 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -25,11 +25,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sampledTriSurfaceMesh.H"
-#include "dictionary.H"
-#include "polyMesh.H"
-#include "polyPatch.H"
-#include "volFields.H"
+#include "treeDataPoint.H"
 #include "meshSearch.H"
+#include "Tuple2.H"
+#include "globalIndex.H"
 
 #include "addToRunTimeSelectionTable.H"
 
@@ -44,6 +43,25 @@ namespace Foam
         sampledTriSurfaceMesh,
         word
     );
+
+    //- Private class for finding nearest
+    //  - global index
+    //  - sqr(distance)
+    typedef Tuple2<scalar, label> nearInfo;
+
+    class nearestEqOp
+    {
+
+    public:
+
+        void operator()(nearInfo& x, const nearInfo& y) const
+        {
+            if (y.first() < x.first())
+            {
+                x = y;
+            }
+        }
+    };
 }
 
 
@@ -63,16 +81,16 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
         (
             name,
             mesh.time().constant(), // instance
-            "triSurface",       // local
-            mesh,               // registry
+            "triSurface",           // local
+            mesh,                   // registry
             IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::NO_WRITE,
+            false
         )
     ),
     needsUpdate_(true),
     cellLabels_(0),
-    pointMap_(0),
-    faceMap_(0)
+    pointToFace_(0)
 {}
 
 
@@ -90,16 +108,16 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
         (
             dict.lookup("surface"),
             mesh.time().constant(), // instance
-            "triSurface",       // local
-            mesh,               // registry
+            "triSurface",           // local
+            mesh,                   // registry
             IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::NO_WRITE,
+            false
         )
     ),
     needsUpdate_(true),
     cellLabels_(0),
-    pointMap_(0),
-    faceMap_(0)
+    pointToFace_(0)
 {}
 
 
@@ -128,8 +146,7 @@ bool Foam::sampledTriSurfaceMesh::expire()
     sampledSurface::clearGeom();
     MeshStorage::clear();
     cellLabels_.clear();
-    pointMap_.clear();
-    faceMap_.clear();
+    pointToFace_.clear();
 
     needsUpdate_ = true;
     return true;
@@ -148,93 +165,137 @@ bool Foam::sampledTriSurfaceMesh::update()
     // Does approximation by looking at the face centres only
     const pointField& fc = surface_.faceCentres();
 
-    cellLabels_.setSize(fc.size());
-
     meshSearch meshSearcher(mesh(), false);
 
+    const indexedOctree<treeDataPoint>& cellCentreTree =
+        meshSearcher.cellCentreTree();
+
+
+    // Global numbering for cells - only used to uniquely identify local cells.
+    globalIndex globalCells(mesh().nCells());
+    List<nearInfo> nearest(fc.size());
+    forAll(nearest, i)
+    {
+        nearest[i].first() = GREAT;
+        nearest[i].second() = labelMax;
+    }
+
+    // Search triangles using nearest on local mesh
     forAll(fc, triI)
     {
-        cellLabels_[triI] = meshSearcher.findCell
+        pointIndexHit nearInfo = cellCentreTree.findNearest
         (
             fc[triI],
-            -1,                 // seed
-            true                // use tree
+            sqr(GREAT)
         );
-    }
-
-
-    boolList include(surface_.size(), true);
-
-    label nNotFound = 0;
-
-    forAll(cellLabels_, triI)
-    {
-        if (cellLabels_[triI] == -1)
+        if (nearInfo.hit())
         {
-            include[triI] = false;
-            nNotFound++;
+            nearest[triI].first() = magSqr(nearInfo.hitPoint()-fc[triI]);
+            nearest[triI].second() = globalCells.toGlobal(nearInfo.index());
         }
     }
-    label nTotalNotFound = returnReduce(nNotFound, sumOp<label>());
+
+    // See which processor has the nearest.
+    Pstream::listCombineGather(nearest, nearestEqOp());
+    Pstream::listCombineScatter(nearest);
+
+    boolList include(surface_.size(), false);
+
+    cellLabels_.setSize(fc.size());
+    cellLabels_ = -1;
+
+    label nFound = 0;
+    forAll(nearest, triI)
+    {
+        if (nearest[triI].second() == labelMax)
+        {
+            // Not found on any processor. How to map?
+        }
+        else if (globalCells.isLocal(nearest[triI].second()))
+        {
+            cellLabels_[triI] = globalCells.toLocal(nearest[triI].second());
+
+            include[triI] = true;
+            nFound++;
+        }
+    }
+
 
     if (debug)
     {
-        Pout<< "surface:" << surface_.size()
-            << "  included:" << surface_.size()-nNotFound
-            << "  total not found" << nTotalNotFound << endl;
+        Pout<< "Local out of faces:" << cellLabels_.size()
+            << " keeping:" << nFound << endl;
     }
 
+    // Now subset the surface. Do not use triSurface::subsetMesh since requires
+    // original surface to be in compact numbering.
+
+    const triSurface& s = surface_;
+
+    // Compact to original triangle
+    labelList faceMap(s.size());
+    // Compact to original points
+    labelList pointMap(s.points().size());
+    // From original point to compact points
+    labelList reversePointMap(s.points().size(), -1);
 
-    // Add to master all triangles are outside all meshes.
     {
-        boolList onAnyProc(surface_.size(), false);
-        Pstream::listCombineGather(onAnyProc, orEqOp<bool>());
+        label newPointI = 0;
+        label newTriI = 0;
 
-        if (Pstream::master())
+        forAll(s, triI)
         {
-            label nAdded = 0;
-            forAll(onAnyProc, triI)
+            if (include[triI])
             {
-                if (!onAnyProc[triI])
+                faceMap[newTriI++] = triI;
+
+                const labelledTri& f = s[triI];
+                forAll(f, fp)
                 {
-                    include[triI] = true;
-                    nAdded++;
+                    if (reversePointMap[f[fp]] == -1)
+                    {
+                        pointMap[newPointI] = f[fp];
+                        reversePointMap[f[fp]] = newPointI++;
+                    }
                 }
             }
+        }
+        faceMap.setSize(newTriI);
+        pointMap.setSize(newPointI);
+    }
 
-            if (debug)
-            {
-                Pout<< "nAdded to master:" << nAdded << endl;
-            }
+    // Subset cellLabels
+    cellLabels_ = UIndirectList<label>(cellLabels_, faceMap)();
+
+    // Store any face per point
+    pointToFace_.setSize(pointMap.size());
+
+    // Create faces and points for subsetted surface
+    faceList& faces = this->storedFaces();
+    faces.setSize(faceMap.size());
+    forAll(faceMap, i)
+    {
+        const triFace& f = s[faceMap[i]];
+        triFace newF
+        (
+            reversePointMap[f[0]],
+            reversePointMap[f[1]],
+            reversePointMap[f[2]]
+        );
+        faces[i] = newF.triFaceFace();
+
+        forAll(newF, fp)
+        {
+            pointToFace_[newF[fp]] = i;
         }
     }
 
-
-    // Now subset the surface
-    triSurface localSurface
-    (
-        surface_.subsetMesh
-        (
-            include,
-            pointMap_,
-            faceMap_
-        )
-    );
-
-    // And convert into faces.
-    faceList& faces = this->storedFaces();
-    faces.setSize(localSurface.size());
-    forAll(localSurface, i)
-    {
-        faces[i] = localSurface[i].triFaceFace();
-    }
-
-    this->storedPoints() = localSurface.points();
+    this->storedPoints() = pointField(s.points(), pointMap);
 
     if (debug)
     {
         print(Pout);
-        Pout << endl;
+        Pout<< endl;
     }
 
     needsUpdate_ = false;
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
index eef8a16224..10ab6edb9d 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
@@ -29,15 +29,13 @@ Description
     A sampledSurface from a triSurfaceMesh. It samples on the points/triangles
     of the triSurface.
 
+    It samples using the cell nearest to the triangle centre so does
+    not check the cell the centre is actually in ...
+
     In parallel every processor just operates on the part of the surface
     where the face centres are inside the mesh. It is then up to the
     caller to stitch the partial surfaces together.
 
-    No check is done for triangle centres being on multiple processors
-    (should never happen but truncation errors ...)
-
-    Any value on a triangle outside the mesh will be set to 0.
-
 SourceFiles
     sampledTriSurfaceMesh.C
 
@@ -75,14 +73,11 @@ class sampledTriSurfaceMesh
         //- Track if the surface needs an update
         mutable bool needsUpdate_;
 
-        //- Local cell labels
+        //- From local surface triangle to mesh cell.
         labelList cellLabels_;
 
         //- From local surface back to surface_
-        labelList pointMap_;
-
-        //- From local surface back to surface_
-        labelList faceMap_;
+        labelList pointToFace_;
 
 
     // Private Member Functions
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
index 0dd08f8f26..b2a8ee956c 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
@@ -36,20 +36,12 @@ Foam::sampledTriSurfaceMesh::sampleField
 ) const
 {
     // One value per face
-    tmp<Field<Type> > tvalues(new Field<Type>(faceMap_.size()));
+    tmp<Field<Type> > tvalues(new Field<Type>(cellLabels_.size()));
     Field<Type>& values = tvalues();
 
-    forAll(faceMap_, i)
+    forAll(cellLabels_, triI)
     {
-        label cellI = cellLabels_[faceMap_[i]];
-        if (cellI != -1)
-        {
-            values[i] = vField[cellI];
-        }
-        else
-        {
-            values[i] = pTraits<Type>::zero;
-        }
+        values[triI] = vField[cellLabels_[triI]];
     }
 
     return tvalues;
@@ -64,24 +56,15 @@ Foam::sampledTriSurfaceMesh::interpolateField
 ) const
 {
     // One value per vertex
-    tmp<Field<Type> > tvalues(new Field<Type>(pointMap_.size()));
+    tmp<Field<Type> > tvalues(new Field<Type>(pointToFace_.size()));
     Field<Type>& values = tvalues();
 
-    forAll(pointMap_, i)
+    forAll(pointToFace_, pointI)
     {
-        label origPointI = pointMap_[i];
-        label origTriI = surface_.pointFaces()[origPointI][0];
+        label triI = pointToFace_[pointI];
+        label cellI = cellLabels_[triI];
 
-        label cellI = cellLabels_[origTriI];
-
-        if (cellI != -1)
-        {
-            values[i] = interpolator.interpolate(points()[i], cellI);
-        }
-        else
-        {
-            values[i] = pTraits<Type>::zero;
-        }
+        values[pointI] = interpolator.interpolate(points()[pointI], cellI);
     }
 
     return tvalues;