diff --git a/applications/utilities/preProcessing/setWaves/setWaves.C b/applications/utilities/preProcessing/setWaves/setWaves.C
index f462135776..314218ed80 100644
--- a/applications/utilities/preProcessing/setWaves/setWaves.C
+++ b/applications/utilities/preProcessing/setWaves/setWaves.C
@@ -102,15 +102,15 @@ int main(int argc, char *argv[])
         );
 
         // Create modelled fields on both cells and points
-        volScalarField height
+        volScalarField h
         (
-            IOobject("height", runTime.timeName(), mesh),
+            IOobject("h", runTime.timeName(), mesh),
             mesh,
             dimensionedScalar("0", dimLength, 0)
         );
-        pointScalarField heightp
+        pointScalarField hp
         (
-            IOobject("heightp", runTime.timeName(), mesh),
+            IOobject("hp", runTime.timeName(), mesh),
             pMesh,
             dimensionedScalar("0", dimLength, 0)
         );
@@ -126,15 +126,15 @@ int main(int argc, char *argv[])
             pMesh,
             dimensionedVector("0", dimLength, vector::zero)
         );
-        volVectorField uLiquid
+        volVectorField uLiq
         (
-            IOobject("uLiquid", runTime.timeName(), mesh),
+            IOobject("uLiq", runTime.timeName(), mesh),
             mesh,
             dimensionedVector("0", dimVelocity, vector::zero)
         );
-        pointVectorField uLiquidp
+        pointVectorField uLiqp
         (
-            IOobject("uLiquidp", runTime.timeName(), mesh),
+            IOobject("uLiqp", runTime.timeName(), mesh),
             pMesh,
             dimensionedVector("0", dimLength, vector::zero)
         );
@@ -145,9 +145,6 @@ int main(int argc, char *argv[])
         // Whether the alpha conditions refer to the liquid phase
         bool liquid = false;
 
-        // The mean velocity of one of the wave patches
-        vector UMeanp = vector::zero;
-
         // Loop the patches, averaging and superimposing wave model data
         forAll(mesh.boundary(), patchi)
         {
@@ -179,8 +176,6 @@ int main(int argc, char *argv[])
             const waveSuperposition& waves =
                 refCast<waveVelocityFvPatchVectorField>(Up).waves();
 
-            UMeanp = waves.UMean();
-
             const bool liquidp =
                 refCast<waveAlphaFvPatchScalarField>(alphap).liquid();
             if (nWaves > 0 && liquidp != liquid)
@@ -198,47 +193,66 @@ int main(int argc, char *argv[])
             const pointField& pts = mesh.points();
 
             // Internal field superposition
-            height.primitiveFieldRef() += waves.height(t, ccs);
-            heightp.primitiveFieldRef() += waves.height(t, pts);
-            uGas.primitiveFieldRef() += waves.UGas(t, ccs) - UMeanp;
-            uGasp.primitiveFieldRef() += waves.UGas(t, pts) - UMeanp;
-            uLiquid.primitiveFieldRef() += waves.ULiquid(t, ccs) - UMeanp;
-            uLiquidp.primitiveFieldRef() += waves.ULiquid(t, pts) - UMeanp;
+            h.primitiveFieldRef() += waves.height(t, ccs);
+            hp.primitiveFieldRef() += waves.height(t, pts);
+            uGas.primitiveFieldRef() += waves.UGas(t, ccs) - waves.UMean();
+            uGasp.primitiveFieldRef() += waves.UGas(t, pts) - waves.UMean();
+            uLiq.primitiveFieldRef() += waves.ULiquid(t, ccs) - waves.UMean();
+            uLiqp.primitiveFieldRef() += waves.ULiquid(t, pts) - waves.UMean();
 
             // Boundary field superposition
             forAll(mesh.boundary(), patchj)
             {
                 const pointField& fcs = mesh.boundary()[patchj].Cf();
-                height.boundaryFieldRef()[patchj] += waves.height(t, fcs);
-                uGas.boundaryFieldRef()[patchj] += waves.UGas(t, fcs) - UMeanp;
-                uLiquid.boundaryFieldRef()[patchj] +=
-                    waves.ULiquid(t, fcs) - UMeanp;
+                h.boundaryFieldRef()[patchj] += waves.height(t, fcs);
+                uGas.boundaryFieldRef()[patchj] +=
+                    waves.UGas(t, fcs) - waves.UMean();
+                uLiq.boundaryFieldRef()[patchj] +=
+                    waves.ULiquid(t, fcs) - waves.UMean();
             }
 
             ++ nWaves;
         }
 
-        // Warn and skip to the next time if no wave patches were found
-        if (nWaves == 0)
-        {
-            WarningInFunction
-                << "No " << waveAlphaFvPatchScalarField::typeName << " or "
-                << waveVelocityFvPatchVectorField::typeName << " patch fields "
-                << "were found. No waves have been set." << endl;
-
-            continue;
-        }
-
         // Create the mean velocity field
         volVectorField UMean
         (
             IOobject("UMean", runTime.timeName(), mesh),
             mesh,
-            dimensionedVector("UMean", dimVelocity, UMeanp)
+            dimensionedVector("UMean", dimVelocity, Zero)
         );
 
-        if (nWaves > 1)
+        if (nWaves == 0)
         {
+            // Warn and skip to the next time if there are no wave patches
+            WarningInFunction
+                << "No " << waveAlphaFvPatchScalarField::typeName << " or "
+                << waveVelocityFvPatchVectorField::typeName << " patch fields "
+                << "were found. No waves have been set." << endl;
+
+            continue;
+        }
+        else if (nWaves == 1)
+        {
+            // Set a mean velocity equal to that on the only wave patch
+            forAll(mesh.boundary(), patchi)
+            {
+                const fvPatchVectorField& Up = U.boundaryField()[patchi];
+                if (!isA<waveVelocityFvPatchVectorField>(Up))
+                {
+                    continue;
+                }
+
+                const waveSuperposition& waves =
+                    refCast<const waveVelocityFvPatchVectorField>(Up).waves();
+
+                UMean == dimensionedVector("UMean", dimVelocity, waves.UMean());
+            }
+        }
+        else if (nWaves > 1)
+        {
+            // Set the mean velocity by distance weighting from the wave patches
+
             // Create weighted average fields for the mean velocity
             volScalarField weight
             (
@@ -256,20 +270,14 @@ int main(int argc, char *argv[])
             // Loop the patches, inverse-distance weighting the mean velocities
             forAll(mesh.boundary(), patchi)
             {
-                fvPatchScalarField& alphap = alpha.boundaryFieldRef()[patchi];
-                fvPatchVectorField& Up = U.boundaryFieldRef()[patchi];
-
-                const bool isWave = isA<waveAlphaFvPatchScalarField>(alphap);
-
-                if (!isWave)
+                const fvPatchVectorField& Up = U.boundaryField()[patchi];
+                if (!isA<waveVelocityFvPatchVectorField>(Up))
                 {
                     continue;
                 }
 
                 const waveSuperposition& waves =
-                    refCast<waveVelocityFvPatchVectorField>(Up).waves();
-
-                UMeanp = waves.UMean();
+                    refCast<const waveVelocityFvPatchVectorField>(Up).waves();
 
                 const volScalarField w
                 (
@@ -279,62 +287,26 @@ int main(int argc, char *argv[])
                      + dimensionedScalar("ySmall", dimLength, SMALL)
                     )
                 );
+
                 weight += w;
                 weightUMean +=
-                    w*dimensionedVector("UMeanp", dimVelocity, UMeanp);
+                    w*dimensionedVector("wUMean", dimVelocity, waves.UMean());
             }
 
             // Complete the average for the mean velocity
             UMean = weightUMean/weight;
         }
 
-        // Set the internal fields
-        alpha.ref() = levelSetFraction(mesh, height, heightp, !liquid);
-        U.ref() =
-            UMean
-          + levelSetAverage
-            (
-                mesh,
-                height,
-                heightp,
-                uGas,
-                uGasp,
-                uLiquid,
-                uLiquidp
-            );
+        // Set the fields
+        alpha == levelSetFraction(h, hp, !liquid);
+        U == UMean + levelSetAverage(mesh, h, hp, uGas, uGasp, uLiq, uLiqp);
 
         // Set the boundary fields
         forAll(mesh.boundary(), patchi)
         {
             fvPatchScalarField& alphap = alpha.boundaryFieldRef()[patchi];
             fvPatchVectorField& Up = U.boundaryFieldRef()[patchi];
-
-            const bool isWave = isA<waveAlphaFvPatchScalarField>(alphap);
-
-            if (!isWave)
-            {
-                alphap ==
-                    levelSetFraction
-                    (
-                        mesh.boundary()[patchi],
-                        height.boundaryField()[patchi],
-                        heightp.boundaryField()[patchi].patchInternalField(),
-                        !liquid
-                    );
-                Up ==
-                    UMean.boundaryField()[patchi]
-                  + levelSetAverage
-                    (
-                        mesh.boundary()[patchi],
-                        height.boundaryField()[patchi],
-                        heightp.boundaryField()[patchi].patchInternalField()(),
-                        uGas.boundaryField()[patchi],
-                        uGasp.boundaryField()[patchi].patchInternalField()(),
-                        uLiquid.boundaryField()[patchi],
-                        uLiquidp.boundaryField()[patchi].patchInternalField()()
-                    );
-            }
-            else
+            if (isA<waveAlphaFvPatchScalarField>(alphap))
             {
                 alphap == refCast<waveAlphaFvPatchScalarField>(alphap).alpha();
                 Up == refCast<waveVelocityFvPatchVectorField>(Up).U();
diff --git a/src/finiteVolume/cfdTools/general/levelSet/levelSet.C b/src/finiteVolume/cfdTools/general/levelSet/levelSet.C
index e5a3fc9489..9442cb2c29 100644
--- a/src/finiteVolume/cfdTools/general/levelSet/levelSet.C
+++ b/src/finiteVolume/cfdTools/general/levelSet/levelSet.C
@@ -30,8 +30,7 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
-Foam::levelSetFraction
+Foam::tmp<Foam::scalarField> Foam::levelSetFraction
 (
     const fvMesh& mesh,
     const scalarField& levelC,
@@ -39,21 +38,8 @@ Foam::levelSetFraction
     const bool above
 )
 {
-    tmp<DimensionedField<scalar, volMesh>> tResult
-    (
-        new DimensionedField<scalar, volMesh>
-        (
-            IOobject
-            (
-                "levelSetFraction",
-                mesh.time().timeName(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar("0", dimless, 0)
-        )
-    );
-    DimensionedField<scalar, volMesh>& result = tResult.ref();
+    tmp<scalarField> tResult(new scalarField(mesh.nCells(), Zero));
+    scalarField& result = tResult.ref();
 
     forAll(result, cI)
     {
@@ -156,4 +142,55 @@ Foam::tmp<Foam::scalarField> Foam::levelSetFraction
     return tResult;
 }
 
+
+Foam::tmp<Foam::volScalarField> Foam::levelSetFraction
+(
+    const volScalarField& levelC,
+    const pointScalarField& levelP,
+    const bool above
+)
+{
+    const fvMesh& mesh = levelC.mesh();
+
+    tmp<volScalarField> tResult
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "levelSetFraction",
+                mesh.time().timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar("0", dimless, 0)
+        )
+    );
+    volScalarField& result = tResult.ref();
+
+    result.primitiveFieldRef() =
+        levelSetFraction
+        (
+            mesh,
+            levelC.primitiveField(),
+            levelP.primitiveField(),
+            above
+        );
+
+    forAll(mesh.boundary(), patchi)
+    {
+        result.boundaryFieldRef()[patchi] =
+            levelSetFraction
+            (
+                mesh.boundary()[patchi],
+                levelC.boundaryField()[patchi],
+                levelP.boundaryField()[patchi].patchInternalField()(),
+                above
+            );
+    }
+
+    return tResult;
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/levelSet/levelSet.H b/src/finiteVolume/cfdTools/general/levelSet/levelSet.H
index f7b7004a4f..3c885f249b 100644
--- a/src/finiteVolume/cfdTools/general/levelSet/levelSet.H
+++ b/src/finiteVolume/cfdTools/general/levelSet/levelSet.H
@@ -34,10 +34,8 @@ SourceFiles:
 #ifndef levelSet_H
 #define levelSet_H
 
-#include "DimensionedField.H"
-#include "fvMesh.H"
-#include "pointMesh.H"
-#include "volMesh.H"
+#include "volFields.H"
+#include "pointFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,15 +47,15 @@ namespace Foam
 //- Calculate the average value of two fields, one on each side of a level set
 //  The level set and the fields are given both on the points and cell-centres.
 template<class Type>
-tmp<DimensionedField<Type, volMesh>> levelSetAverage
+tmp<Field<Type>> levelSetAverage
 (
     const fvMesh& mesh,
     const scalarField& levelC,
     const scalarField& levelP,
-    const DimensionedField<Type, volMesh>& positiveC,
-    const DimensionedField<Type, pointMesh>& positiveP,
-    const DimensionedField<Type, volMesh>& negativeC,
-    const DimensionedField<Type, pointMesh>& negativeP
+    const Field<Type>& positiveC,
+    const Field<Type>& positiveP,
+    const Field<Type>& negativeC,
+    const Field<Type>& negativeP
 );
 
 //- As the above overload, but on the faces of a patch
@@ -73,10 +71,22 @@ tmp<Field<Type>> levelSetAverage
     const Field<Type>& negativeP
 );
 
+//- As the above oveload, but both in cells and on patch faces
+template<class Type>
+tmp<GeometricField<Type, fvPatchField, volMesh>> levelSetAverage
+(
+    const volScalarField& levelC,
+    const pointScalarField& levelP,
+    const GeometricField<Type, fvPatchField, volMesh>& positiveC,
+    const GeometricField<Type, pointPatchField, pointMesh>& positiveP,
+    const GeometricField<Type, fvPatchField, volMesh>& negativeC,
+    const GeometricField<Type, pointPatchField, pointMesh>& negativeP
+);
+
 //- Calculate the volume-fraction that a level set occupies. This gives the the
 //  same result as levelSetAverage if the fields passed to the latter are
 //  uniformly 0 and 1. The above flag flips the direction.
-tmp<DimensionedField<scalar, volMesh>> levelSetFraction
+tmp<scalarField> levelSetFraction
 (
     const fvMesh& mesh,
     const scalarField& levelC,
@@ -93,6 +103,14 @@ tmp<scalarField> levelSetFraction
     const bool above
 );
 
+//- As the above oveload, but both in cells and on patch faces
+tmp<volScalarField> levelSetFraction
+(
+    const volScalarField& levelC,
+    const pointScalarField& levelP,
+    const bool above
+);
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C b/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C
index de89bbb76d..9d901f5d9e 100644
--- a/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C
+++ b/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C
@@ -31,32 +31,19 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type>
-Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh>> Foam::levelSetAverage
+Foam::tmp<Foam::Field<Type>> Foam::levelSetAverage
 (
     const fvMesh& mesh,
     const scalarField& levelC,
     const scalarField& levelP,
-    const DimensionedField<Type, volMesh>& positiveC,
-    const DimensionedField<Type, pointMesh>& positiveP,
-    const DimensionedField<Type, volMesh>& negativeC,
-    const DimensionedField<Type, pointMesh>& negativeP
+    const Field<Type>& positiveC,
+    const Field<Type>& positiveP,
+    const Field<Type>& negativeC,
+    const Field<Type>& negativeP
 )
 {
-    tmp<DimensionedField<Type, volMesh>> tResult
-    (
-        new DimensionedField<Type, volMesh>
-        (
-            IOobject
-            (
-                positiveC.name() + ":levelSetAverage",
-                mesh.time().timeName(),
-                mesh
-            ),
-            mesh,
-            dimensioned<Type>("0", positiveC.dimensions(), Zero)
-        )
-    );
-    DimensionedField<Type, volMesh>& result = tResult.ref();
+    tmp<Field<Type>> tResult(new Field<Type>(mesh.nCells(), Zero));
+    Field<Type>& result = tResult.ref();
 
     forAll(result, cI)
     {
@@ -184,4 +171,65 @@ Foam::tmp<Foam::Field<Type>> Foam::levelSetAverage
 }
 
 
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
+Foam::levelSetAverage
+(
+    const volScalarField& levelC,
+    const pointScalarField& levelP,
+    const GeometricField<Type, fvPatchField, volMesh>& positiveC,
+    const GeometricField<Type, pointPatchField, pointMesh>& positiveP,
+    const GeometricField<Type, fvPatchField, volMesh>& negativeC,
+    const GeometricField<Type, pointPatchField, pointMesh>& negativeP
+)
+{
+    const fvMesh& mesh = levelC.mesh();
+
+    tmp<GeometricField<Type, fvPatchField, volMesh>> tResult
+    (
+        new GeometricField<Type, fvPatchField, volMesh>
+        (
+            IOobject
+            (
+                positiveC.name() + ":levelSetAverage",
+                mesh.time().timeName(),
+                mesh
+            ),
+            mesh,
+            dimensioned<Type>("0", positiveC.dimensions(), Zero)
+        )
+    );
+    GeometricField<Type, fvPatchField, volMesh>& result = tResult.ref();
+
+    result.primitiveFieldRef() =
+        levelSetAverage
+        (
+            mesh,
+            levelC.primitiveField(),
+            levelP.primitiveField(),
+            positiveC.primitiveField(),
+            positiveP.primitiveField(),
+            negativeC.primitiveField(),
+            negativeP.primitiveField()
+        );
+
+    forAll(mesh.boundary(), patchi)
+    {
+        result.boundaryField()[patchi] =
+            levelSetAverage
+            (
+                mesh.boundary()[patchi],
+                levelC.boundaryField()[patchi],
+                levelP.boundaryField()[patchi].patchInternalField()(),
+                positiveC.boundaryField()[patchi],
+                negativeP.boundaryField()[patchi].patchInternalField()(),
+                positiveC.boundaryField()[patchi],
+                negativeP.boundaryField()[patchi].patchInternalField()()
+            );
+    }
+
+    return tResult;
+}
+
+
 // ************************************************************************* //