diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index 859db2f9af..3e3d56c6ed 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -76,6 +76,63 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+const Foam::indexedOctree<Foam::treeDataFace>&
+Foam::sampledTriSurfaceMesh::nonCoupledboundaryTree() const
+{
+    // Variant of meshSearch::boundaryTree() that only does non-coupled
+    // boundary faces.
+
+    if (!boundaryTreePtr_.valid())
+    {
+        // all non-coupled boundary faces (not just walls)
+        const polyBoundaryMesh& patches = mesh().boundaryMesh();
+
+        labelList bndFaces(mesh().nFaces()-mesh().nInternalFaces());
+        label bndI = 0;
+        forAll(patches, patchI)
+        {
+            const polyPatch& pp = patches[patchI];
+            if (!pp.coupled())
+            {
+                forAll(pp, i)
+                {
+                    bndFaces[bndI++] = pp.start()+i;
+                }
+            }
+        }
+        bndFaces.setSize(bndI);
+
+
+        treeBoundBox overallBb(mesh().points());
+        Random rndGen(123456);
+        overallBb = overallBb.extend(rndGen, 1E-4);
+        overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+        boundaryTreePtr_.reset
+        (
+            new indexedOctree<treeDataFace>
+            (
+                treeDataFace    // all information needed to search faces
+                (
+                    false,                      // do not cache bb
+                    mesh(),
+                    bndFaces                    // boundary faces only
+                ),
+                overallBb,                      // overall search domain
+                8,                              // maxLevel
+                10,                             // leafsize
+                3.0                             // duplicity
+            )
+        );
+    }
+
+    return boundaryTreePtr_();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
@@ -159,6 +216,8 @@ bool Foam::sampledTriSurfaceMesh::expire()
 
     sampledSurface::clearGeom();
     MeshStorage::clear();
+
+    boundaryTreePtr_.clear();
     sampleElements_.clear();
     samplePoints_.clear();
 
@@ -224,7 +283,10 @@ bool Foam::sampledTriSurfaceMesh::update()
     {
         // Search for nearest boundaryFace
 
-        const indexedOctree<treeDataFace>& bTree = meshSearcher.boundaryTree();
+        ////- Search on all (including coupled) boundary faces
+        //const indexedOctree<treeDataFace>& bTree = meshSearcher.boundaryTree()
+        //- Search on all non-coupled boundary faces
+        const indexedOctree<treeDataFace>& bTree = nonCoupledboundaryTree();
 
         forAll(fc, triI)
         {
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
index 20426e576b..c54bf28dc0 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
@@ -28,12 +28,31 @@ 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 ...
+    - it either samples cells or (non-coupled) boundary faces
 
-    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.
+    - 4 different modes:
+        - source=cells, interpolate=false:
+            finds per triangle centre the nearest cell centre and uses its value
+        - source=cells, interpolate=true
+            finds per triangle centre the nearest cell centre.
+            Per surface point checks if this nearest cell is the one containing
+            point; otherwise projects the point onto the nearest point on
+            the boundary of the cell (to make sure interpolateCellPoint
+            gets a valid location)
+
+        - source=boundaryFaces, interpolate=false:
+            finds per triangle centre the nearest point on the boundary
+            (uncoupled faces only) and uses the value (or 0 if the nearest
+            is on an empty boundary)
+        - source=boundaryFaces, interpolate=true:
+            finds per triangle centre the nearest point on the boundary
+            (uncoupled faces only).
+            Per surface point projects the point onto this boundary face
+            (to make sure interpolateCellPoint gets a valid location)
+
+    - since it finds a nearest per triangle each triangle is guaranteed
+    to be on one processor only. So after stitching (by sampledSurfaces)
+    the original surface should be complete.
 
 SourceFiles
     sampledTriSurfaceMesh.C
@@ -52,6 +71,8 @@ SourceFiles
 namespace Foam
 {
 
+class treeDataFace;
+
 /*---------------------------------------------------------------------------*\
                        Class sampledTriSurfaceMesh Declaration
 \*---------------------------------------------------------------------------*/
@@ -88,6 +109,9 @@ private:
         //- Track if the surface needs an update
         mutable bool needsUpdate_;
 
+        //- Search tree for all non-coupled boundary faces
+        mutable autoPtr<indexedOctree<treeDataFace> > boundaryTreePtr_;
+
         //- From local surface triangle to mesh cell/face.
         labelList sampleElements_;
 
@@ -97,6 +121,9 @@ private:
 
     // Private Member Functions
 
+        //- Get tree of all non-coupled boundary faces
+        const indexedOctree<treeDataFace>& nonCoupledboundaryTree() const;
+
         //- sample field on faces
         template <class Type>
         tmp<Field<Type> > sampleField