diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
index a147425d9d..1cb8bf418e 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
@@ -61,36 +61,13 @@ Foam::sampledIsoSurfaceCell::sampleOnPoints
 {
     updateGeometry();  // Recreate geometry if time has changed
 
-    const labelList& elements = meshCells_;
-
-    // One value per point
-    auto tvalues = tmp<Field<Type>>::New(points().size());
-    auto& values = tvalues.ref();
-
-    const faceList& fcs = faces();
-    const pointField& pts = points();
-
-    bitSet pointDone(points().size());
-
-    forAll(faces(), cutFacei)
-    {
-        const face& f = fcs[cutFacei];
-        const label celli = elements[cutFacei];
-
-        for (const label pointi : f)
-        {
-            if (pointDone.set(pointi))
-            {
-                values[pointi] = interpolator.interpolate
-                (
-                    pts[pointi],
-                    celli
-                );
-            }
-        }
-    }
-
-    return tvalues;
+    return sampledSurface::sampleOnPoints
+    (
+        interpolator,
+        meshCells_,
+        faces(),
+        points()
+    );
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
index 81943de7a5..8c3b67c55e 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
@@ -61,36 +61,13 @@ Foam::sampledIsoSurfaceTopo::sampleOnPoints
 {
     updateGeometry();  // Recreate geometry if time has changed
 
-    const labelList& elements = meshCells_;
-
-    // One value per point
-    auto tvalues = tmp<Field<Type>>::New(points().size());
-    auto& values = tvalues.ref();
-
-    const faceList& fcs = faces();
-    const pointField& pts = points();
-
-    bitSet pointDone(points().size());
-
-    forAll(faces(), cutFacei)
-    {
-        const face& f = fcs[cutFacei];
-        const label celli = elements[cutFacei];
-
-        for (const label pointi : f)
-        {
-            if (pointDone.set(pointi))
-            {
-                values[pointi] = interpolator.interpolate
-                (
-                    pts[pointi],
-                    celli
-                );
-            }
-        }
-    }
-
-    return tvalues;
+    return sampledSurface::sampleOnPoints
+    (
+        interpolator,
+        meshCells_,
+        faces(),
+        points()
+    );
 }
 
 
diff --git a/src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurfaceTemplates.C b/src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurfaceTemplates.C
index 31b34783ac..4a0442e9bc 100644
--- a/src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurfaceTemplates.C
@@ -53,38 +53,13 @@ Foam::sampledCuttingSurface::sampleOnPoints
     const interpolation<Type>& interpolator
 ) const
 {
-    // elements to sample
-    const labelList& elements = meshCells();
-
-    // One value per point.
-    // Initialize with Zero to handle missed/degenerate faces
-
-    auto tvalues = tmp<Field<Type>>::New(points().size(), Zero);
-    auto& values = tvalues.ref();
-
-    bitSet pointDone(points().size());
-
-    const faceList& fcs = faces();
-
-    forAll(fcs, facei)
-    {
-        const face& f = fcs[facei];
-        const label celli = elements[facei];
-
-        for (const label pointi : f)
-        {
-            if (pointDone.set(pointi))
-            {
-                values[pointi] = interpolator.interpolate
-                (
-                    points()[pointi],
-                    celli
-                );
-            }
-        }
-    }
-
-    return tvalues;
+    return sampledSurface::sampleOnPoints
+    (
+        interpolator,
+        meshCells(),
+        faces(),
+        points()
+    );
 }
 
 
diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C b/src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C
index 027cefb5e8..552c9baf98 100644
--- a/src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C
+++ b/src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,38 +54,13 @@ Foam::sampledPlane::sampleOnPoints
     const interpolation<Type>& interpolator
 ) const
 {
-    // elements to sample
-    const labelList& elements = meshCells();
-
-    // One value per point.
-    // Initialize with Zero to handle missed/degenerate faces
-
-    auto tvalues = tmp<Field<Type>>::New(points().size(), Zero);
-    auto& values = tvalues.ref();
-
-    bitSet pointDone(points().size());
-
-    const faceList& fcs = faces();
-
-    forAll(fcs, facei)
-    {
-        const face& f = fcs[facei];
-        const label celli = elements[facei];
-
-        for (const label pointi : f)
-        {
-            if (pointDone.set(pointi))
-            {
-                values[pointi] = interpolator.interpolate
-                (
-                    points()[pointi],
-                    celli
-                );
-            }
-        }
-    }
-
-    return tvalues;
+    return sampledSurface::sampleOnPoints
+    (
+        interpolator,
+        meshCells(),
+        faces(),
+        points()
+    );
 }
 
 
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
index d135ab0693..d6490746d8 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
@@ -132,7 +132,7 @@ protected:
 
     // Protected Member Functions
 
-        //- General loop for sampling elements to faces.
+        //- Loop for sampling volume elements to faces.
         //  The defaultValue is used for invalid (negative) elements
         template<class Type>
         static tmp<Field<Type>> sampleOnFaces
@@ -145,6 +145,17 @@ protected:
         );
 
 
+        //- Loop for interpolating volume elements to face points.
+        template<class Type>
+        static tmp<Field<Type>> sampleOnPoints
+        (
+            const interpolation<Type>& interpolator,
+            const labelUList& elements,
+            const faceList& fcs,
+            const pointField& pts
+        );
+
+
         //- Create cell values by averaging the point values
         template<class Type>
         static tmp<GeometricField<Type, fvPatchField, volMesh>> pointAverage
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C b/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
index 2b494d5027..f283b7eea9 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
@@ -74,6 +74,56 @@ Foam::sampledSurface::sampleOnFaces
 }
 
 
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::sampledSurface::sampleOnPoints
+(
+    const interpolation<Type>& interpolator,
+    const labelUList& elements,
+    const faceList& fcs,
+    const pointField& pts
+)
+{
+    const label len = elements.size();
+
+    if (len != fcs.size())
+    {
+        FatalErrorInFunction
+            << "size mismatch: "
+            << "sampled elements (" << len
+            << ") != faces (" << fcs.size() << ')'
+            << exit(FatalError);
+    }
+
+    // One value per point
+    // Initialize with Zero to handle missed/degenerate faces
+    auto tvalues = tmp<Field<Type>>::New(pts.size(), Zero);
+    auto& values = tvalues.ref();
+
+    bitSet pointDone(pts.size());
+
+    forAll(fcs, facei)
+    {
+        const face& f = fcs[facei];
+        const label celli = elements[facei];
+
+        for (const label pointi : f)
+        {
+            if (pointDone.set(pointi))
+            {
+                values[pointi] = interpolator.interpolate
+                (
+                    pts[pointi],
+                    celli
+                );
+            }
+        }
+    }
+
+    return tvalues;
+}
+
+
 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
 Foam::sampledSurface::pointAverage
diff --git a/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFacesTemplates.C b/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFacesTemplates.C
index f8590657e9..8e3e8f04f8 100644
--- a/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFacesTemplates.C
+++ b/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFacesTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,27 +42,15 @@ Foam::sampledThresholdCellFaces::sampleOnFaces
     const interpolation<Type>& sampler
 ) const
 {
-    updateGeometry(); // Recreate geometry if time has changed
+    updateGeometry();  // Recreate geometry if time has changed
 
-    const labelList& elements = meshCells_;
-
-    const label len = faces().size();
-
-    auto tvalues = tmp<Field<Type>>::New(len);
-    auto& values = tvalues.ref();
-
-    const faceList& fcs = faces();
-    const pointField& pts = points();
-
-    for (label i=0; i < len; ++i)
-    {
-        const label celli = elements[i];
-        const point pt = fcs[i].centre(pts);
-
-        values[i] = sampler.interpolate(pt, celli);
-    }
-
-    return tvalues;
+    return sampledSurface::sampleOnFaces
+    (
+        sampler,
+        meshCells_,
+        faces(),
+        points()
+    );
 }
 
 
@@ -73,38 +61,15 @@ Foam::sampledThresholdCellFaces::sampleOnPoints
     const interpolation<Type>& interpolator
 ) const
 {
-    updateGeometry(); // Recreate geometry if time has changed
+    updateGeometry();  // Recreate geometry if time has changed
 
-    const labelList& elements = meshCells_;
-
-    // One value per point
-    auto tvalues = tmp<Field<Type>>::New(points().size(), Zero);
-    auto& values = tvalues.ref();
-
-    bitSet pointDone(points().size());
-
-    const faceList& fcs = faces();
-    const pointField& pts = points();
-
-    forAll(fcs, i)
-    {
-        const face& f = fcs[i];
-        const label celli = elements[i];
-
-        for (const label pointi : f)
-        {
-            if (pointDone.set(pointi))
-            {
-                values[pointi] = interpolator.interpolate
-                (
-                    pts[pointi],
-                    celli
-                );
-            }
-        }
-    }
-
-    return tvalues;
+    return sampledSurface::sampleOnPoints
+    (
+        interpolator,
+        meshCells_,
+        faces(),
+        points()
+    );
 }