diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
index 888b51b9c6..92e8694874 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
@@ -284,6 +284,7 @@ Foam::distanceSurface::distanceSurface
     distance_(readScalar(dict.lookup("distance"))),
     signed_(readBool(dict.lookup("signed"))),
     regularise_(dict.lookupOrDefault("regularise", true)),
+    average_(dict.lookupOrDefault("average", false)),
     zoneName_(word::null),
     needsUpdate_(true),
     isoSurfPtr_(NULL),
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
index b46546abbf..9fd2ffd9f0 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
@@ -68,6 +68,9 @@ class distanceSurface
         //- Whether to coarsen
         const Switch regularise_;
 
+        //- Whether to recalculate cell values as average of point values
+        const Switch average_;
+
         //- zone name (if restricted to zones)
         word zoneName_;
 
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
index 1a66028f95..ea54214fb6 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
@@ -62,7 +62,15 @@ Foam::distanceSurface::interpolateField
     );
 
     // Sample.
-    return surface().interpolate(volFld, pointFld());
+    return surface().interpolate
+    (
+        (
+            average_
+          ? pointAverage(pointFld())()
+          : volFld
+        ),
+        pointFld()
+    );
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 235af8607a..11a0d74cae 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -164,7 +164,10 @@ void Foam::sampledIsoSurface::getIsoFields() const
         // point field.
         if (average_)
         {
-            storedVolFieldPtr_.reset(average(fvm, *pointFieldPtr_).ptr());
+            storedVolFieldPtr_.reset
+            (
+                pointAverage(*pointFieldPtr_).ptr()
+            );
             volFieldPtr_ = storedVolFieldPtr_.operator->();
         }
 
@@ -265,7 +268,7 @@ void Foam::sampledIsoSurface::getIsoFields() const
         {
             storedVolSubFieldPtr_.reset
             (
-                average(subFvm, *pointSubFieldPtr_).ptr()
+                pointAverage(*pointSubFieldPtr_).ptr()
             );
             volSubFieldPtr_ = storedVolSubFieldPtr_.operator->();
         }
@@ -286,99 +289,6 @@ void Foam::sampledIsoSurface::getIsoFields() const
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::sampledIsoSurface::average
-(
-    const fvMesh& mesh,
-    const pointScalarField& pfld
-) const
-{
-    tmp<volScalarField> tcellAvg
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "cellAvg",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh,
-            dimensionedScalar("zero", dimless, scalar(0.0))
-        )
-    );
-    volScalarField& cellAvg = tcellAvg();
-
-    labelField nPointCells(mesh.nCells(), 0);
-    {
-        for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
-        {
-            const labelList& pCells = mesh.pointCells(pointI);
-
-            forAll(pCells, i)
-            {
-                label cellI = pCells[i];
-
-                cellAvg[cellI] += pfld[pointI];
-                nPointCells[cellI]++;
-            }
-        }
-    }
-    forAll(cellAvg, cellI)
-    {
-        cellAvg[cellI] /= nPointCells[cellI];
-    }
-    // Give value to calculatedFvPatchFields
-    cellAvg.correctBoundaryConditions();
-
-    return tcellAvg;
-}
-
-
-Foam::tmp<Foam::pointScalarField> Foam::sampledIsoSurface::average
-(
-    const pointMesh& pMesh,
-    const volScalarField& fld
-) const
-{
-    tmp<pointScalarField> tpointAvg
-    (
-        new pointScalarField
-        (
-            IOobject
-            (
-                "pointAvg",
-                fld.time().timeName(),
-                fld.db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            pMesh,
-            dimensionedScalar("zero", dimless, scalar(0.0))
-        )
-    );
-    pointScalarField& pointAvg = tpointAvg();
-
-    for (label pointI = 0; pointI < fld.mesh().nPoints(); pointI++)
-    {
-        const labelList& pCells = fld.mesh().pointCells(pointI);
-
-        forAll(pCells, i)
-        {
-            pointAvg[pointI] += fld[pCells[i]];
-        }
-        pointAvg[pointI] /= pCells.size();
-    }
-    // Give value to calculatedFvPatchFields
-    pointAvg.correctBoundaryConditions();
-
-    return tpointAvg;
-}
-
-
 bool Foam::sampledIsoSurface::updateGeometry() const
 {
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
index ccaa4d8750..990ee4c476 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -118,18 +118,6 @@ class sampledIsoSurface
         //- Get fields needed to recreate iso surface.
         void getIsoFields() const;
 
-        tmp<volScalarField> average
-        (
-            const fvMesh&,
-            const pointScalarField&
-        ) const;
-
-        tmp<pointScalarField> average
-        (
-            const pointMesh&,
-            const volScalarField& fld
-        ) const;
-
         //- Create iso surface (if time has changed)
         //  Do nothing (and return false) if no update was needed
         bool updateGeometry() const;
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
index a8dfa78506..40f05bd206 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
@@ -71,7 +71,15 @@ Foam::sampledIsoSurface::interpolateField
             volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
 
         // Sample.
-        return surface().interpolate(volSubFld, tpointSubFld());
+        return surface().interpolate
+        (
+            (
+                average_
+              ? pointAverage(tpointSubFld())()
+              : volSubFld
+            ),
+            tpointSubFld()
+        );
     }
     else
     {
@@ -79,7 +87,15 @@ Foam::sampledIsoSurface::interpolateField
             volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
 
         // Sample.
-        return surface().interpolate(volFld, tpointFld());
+        return surface().interpolate
+        (
+            (
+                average_
+              ? pointAverage(tpointFld())()
+              : volFld
+            ),
+            tpointFld()
+        );
     }
 }
 
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
index 2bcf3784b6..5e369ed533 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -254,6 +254,7 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
     plane_(dict),
     mergeTol_(dict.lookupOrDefault("mergeTol", 1E-6)),
     regularise_(dict.lookupOrDefault("regularise", true)),
+    average_(dict.lookupOrDefault("average", false)),
     zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
     exposedPatchName_(word::null),
     needsUpdate_(true),
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
index fcf14db8bd..54e0f6627c 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
@@ -66,6 +66,9 @@ class sampledCuttingPlane
         //- Whether to coarsen
         const Switch regularise_;
 
+        //- Whether to recalculate cell values as average of point values
+        const Switch average_;
+
         //- zone name/index (if restricted to zones)
         mutable cellZoneID zoneID_;
 
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
index ad1d97ffe9..ea0ed6bf47 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
@@ -67,7 +67,15 @@ Foam::sampledCuttingPlane::interpolateField
             volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
 
         // Sample.
-        return surface().interpolate(volSubFld, tpointSubFld());
+        return surface().interpolate
+        (
+            (
+                average_
+              ? pointAverage(tpointSubFld())()
+              : volSubFld
+            ),
+            tpointSubFld()
+        );
     }
     else
     {
@@ -77,7 +85,15 @@ Foam::sampledCuttingPlane::interpolateField
         );
 
         // Sample.
-        return surface().interpolate(volFld, tpointFld());
+        return surface().interpolate
+        (
+            (
+                average_
+              ? pointAverage(tpointFld())()
+              : volFld
+            ),
+            tpointFld()
+        );
     }
 }
 
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
index 265b5c3db7..837e8aa159 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
@@ -312,6 +312,13 @@ public:
         //- Project field onto surface
         tmp<Field<vector> > project(const Field<tensor>&) const;
 
+        //- Interpolate from points to cell centre
+        template<class Type>
+        tmp<GeometricField<Type, fvPatchField, volMesh> > pointAverage
+        (
+            const GeometricField<Type, pointPatchField, pointMesh>& pfld
+        ) const;
+
         //- Sample field on surface
         virtual tmp<scalarField> sample
         (
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C b/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
index 3aff356cef..d943792bd2 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
@@ -155,4 +155,59 @@ Foam::sampledSurface::project
 }
 
 
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::sampledSurface::pointAverage
+(
+    const GeometricField<Type, pointPatchField, pointMesh>& pfld
+) const
+{
+    const fvMesh& mesh = dynamic_cast<const fvMesh&>(pfld.mesh()());
+
+    tmp<GeometricField<Type, fvPatchField, volMesh> > tcellAvg
+    (
+        new GeometricField<Type, fvPatchField, volMesh>
+        (
+            IOobject
+            (
+                "cellAvg",
+                mesh.time().timeName(),
+                pfld.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimensioned<Type>("zero", dimless, pTraits<Type>::zero)
+        )
+    );
+    GeometricField<Type, fvPatchField, volMesh>& cellAvg = tcellAvg();
+
+    labelField nPointCells(mesh.nCells(), 0);
+    {
+        for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
+        {
+            const labelList& pCells = mesh.pointCells(pointI);
+
+            forAll(pCells, i)
+            {
+                label cellI = pCells[i];
+
+                cellAvg[cellI] += pfld[pointI];
+                nPointCells[cellI]++;
+            }
+        }
+    }
+    forAll(cellAvg, cellI)
+    {
+        cellAvg[cellI] /= nPointCells[cellI];
+    }
+    // Give value to calculatedFvPatchFields
+    cellAvg.correctBoundaryConditions();
+
+    return tcellAvg;
+}
+
+
+
 // ************************************************************************* //