diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
index 48bd03308e725c436bf0104a24bc5eccbf02b858..b4f83b373621b5448206ab761b8491d364b713c8 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
@@ -58,14 +58,14 @@ Foam::sampledTriSurfaceMesh::sampleOnFaces
     auto& values = tvalues.ref();
 
     const polyBoundaryMesh& pbm = mesh().boundaryMesh();
-
+    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
 
     // Create flat boundary field
-    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
 
     Field<Type> bVals(nBnd, Zero);
 
     const auto& bField = sampler.psi().boundaryField();
+
     forAll(bField, patchi)
     {
         const label bFacei = (pbm[patchi].start() - mesh().nInternalFaces());
diff --git a/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.C b/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.C
index 2e02fa1982d287cce2473cdd9cc772f3d475434a..dce44e0e63c5ec1b473edd66492b3bfa335d211c 100644
--- a/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.C
+++ b/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.H b/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.H
index 4dd47e2c8922d2a25321c8190ecfb3ef5400ed91..c7a9b90b9dd1df344594a798743a9218534b84ea 100644
--- a/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.H
+++ b/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,6 +91,13 @@ class surfMeshSampleDiscrete
         //- Transfer mesh content from SurfaceSource to surfMesh
         void transferContent();
 
+        //- Sample field (from volume field) onto surface faces
+        template<class Type>
+        tmp<Field<Type>> sampleOnFaces
+        (
+            const interpolation<Type>& sampler
+        ) const;
+
         //- Sample field on surface.
         template<class Type>
         bool sampleType
diff --git a/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscreteTemplates.C b/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscreteTemplates.C
index 4d56412539d5a676b0ed9ee993c0b80e87f53d94..78cf386c52147f43a567a4960f4d91f6254d3fb5 100644
--- a/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscreteTemplates.C
+++ b/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscreteTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,11 +24,32 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "surfMeshSampleDiscrete.H"
-#include "dimensionedType.H"
-#include "error.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::surfMeshSampleDiscrete::sampleOnFaces
+(
+    const interpolation<Type>& sampler
+) const
+{
+    if (onBoundary())
+    {
+        return SurfaceSource::sampleOnFaces(sampler);
+    }
+
+    // Sample cells
+    return surfMeshSample::sampleOnFaces
+    (
+        sampler,
+        sampleElements(),  //< sampleElements == meshCells
+        surface().faces(),
+        surface().points()
+    );
+}
+
+
 template<class Type>
 bool Foam::surfMeshSampleDiscrete::sampleType
 (
@@ -38,7 +59,7 @@ bool Foam::surfMeshSampleDiscrete::sampleType
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
 
-    const auto volFldPtr =
+    const auto* volFldPtr =
         SurfaceSource::mesh().lookupObjectPtr<VolFieldType>(fieldName);
 
     if (!volFldPtr)
@@ -48,52 +69,11 @@ bool Foam::surfMeshSampleDiscrete::sampleType
 
     auto samplerPtr = interpolation<Type>::New(sampleScheme, *volFldPtr);
 
-    getOrCreateSurfField<Type>(*volFldPtr).field()
-        = SurfaceSource::sampleOnFaces(*samplerPtr);
+    getOrCreateSurfField<Type>(*volFldPtr).field() =
+        sampleOnFaces(*samplerPtr);
 
     return true;
 }
 
 
-// template<class Type>
-// Foam::tmp<Foam::DimensionedField<Type, Foam::surfGeoMesh>>
-// Foam::surfMeshSampleDiscrete::sampleOnFaces
-// (
-//     const GeometricField<Type, fvPatchField, volMesh>& fld
-// ) const
-// {
-//     typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
-//
-//     tmp<Field<Type>> tfield = SurfaceSource::sampleOnFaces(fld);
-//     SurfFieldType& result = getOrCreateSurfField<Type>(fld);
-//
-//     // need to verify if this will be needed (in the future)
-//     const surfMesh& s = surface();
-//     if (result.size() != s.size())
-//     {
-//         // maybe resampling changed the surfMesh,
-//         // but the existing surfField wasn't updated
-//
-//         result.resize(s.size(), Zero);
-//     }
-//
-//     if (result.size() != sampleElements().size())
-//     {
-//         FatalErrorInFunction
-//             << "mismatch in field/mesh sizes "
-//             << result.name() << nl
-//             << " field has " << result.size()
-//             << " but sampleElements has " << sampleElements().size() << nl
-//             << endl
-//             << exit(FatalError);
-//
-//         Info<< "WARNING: "
-//             << endl;
-//     }
-//
-//     result = tfield;
-//     return result;
-// }
-
-
 // ************************************************************************* //
diff --git a/src/sampling/surface/triSurfaceMesh/discreteSurface.C b/src/sampling/surface/triSurfaceMesh/discreteSurface.C
index ad7ea4313199ad84d8ac295defa2c061818e4aa3..2c67b787c4b23f4b3d3778526df7c9559b58556d 100644
--- a/src/sampling/surface/triSurfaceMesh/discreteSurface.C
+++ b/src/sampling/surface/triSurfaceMesh/discreteSurface.C
@@ -744,6 +744,7 @@ bool Foam::discreteSurface::interpolate() const
     return interpolate_;
 }
 
+
 bool Foam::discreteSurface::interpolate(bool b)
 {
     if (interpolate_ == b)
diff --git a/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C b/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C
index a1a7b8a2c9a13805eb7e0047f63e44c22371e4d3..74256346f42d0a65645a07f77f5553649d9c7d93 100644
--- a/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C
+++ b/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C
@@ -175,24 +175,26 @@ Foam::discreteSurface::sampleOnFaces
 
         Field<Type> bVals(nBnd, Zero);
 
-        forAll(vField.boundaryField(), patchi)
+        const auto& bField = vField.boundaryField();
+
+        forAll(bField, patchi)
         {
-            label bFacei = pbm[patchi].start() - mesh().nInternalFaces();
+            const label bFacei = pbm[patchi].start() - mesh().nInternalFaces();
 
             SubList<Type>
             (
                 bVals,
-                vField.boundaryField()[patchi].size(),
+                bField[patchi].size(),
                 bFacei
-            ) = vField.boundaryField()[patchi];
+            ) = bField[patchi];
         }
 
         // Sample in flat boundary field
 
         for (label i=0; i < len; ++i)
         {
-            label facei = elements[i];
-            values[i] = bVals[facei-mesh().nInternalFaces()];
+            const label bFacei = (elements[i] - mesh().nInternalFaces());
+            values[i] = bVals[bFacei];
         }
     }
 
@@ -230,7 +232,7 @@ Foam::discreteSurface::sampleOnPoints
 
         forAll(samplePoints_, pointi)
         {
-            label facei = sampleElements_[pointi];
+            const label facei = sampleElements_[pointi];
 
             values[pointi] = interpolator.interpolate
             (