From 45f73bf64f3cfa5d7de9900d2d709fc42255a2ca Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Sat, 23 Apr 2016 23:07:28 +0100
Subject: [PATCH] GeometricField: New non-const access function
 boundaryFieldRef()

There is a need to specify const or non-const access to a non-const
object which is not currently possible with the "boundaryField()" access
function the const-ness of the return of which is defined by the
const-ness of the object for which it is called.  For consistency with
the latest "tmp" storage class in which non-const access is obtained
with the "ref()" function it is proposed to replace the non-const form
of "boundaryField()" with "boundaryFieldRef()".

Thanks to Mattijs Janssens for starting the process of migration to
"boundaryFieldRef()" and providing a patch for the OpenFOAM and
finiteVolume libraries.
---
 .../GeometricField/GeometricField.C           | 30 ++++-----
 .../GeometricField/GeometricField.H           | 10 ++-
 .../GeometricField/GeometricFieldFunctions.C  | 25 +++++---
 .../GeometricField/GeometricFieldFunctionsM.C | 23 ++++---
 .../GeometricField/MapGeometricFields.H       |  2 +-
 .../GeometricScalarField.C                    | 16 ++---
 .../transformGeometricField.C                 |  4 +-
 .../general/CorrectPhi/correctUphiBCs.C       | 42 ++++++------
 .../cfdTools/general/MRF/MRFZone.C            | 19 ++++--
 .../cfdTools/general/MRF/MRFZoneList.C        |  9 ++-
 .../cfdTools/general/MRF/MRFZoneTemplates.C   |  9 ++-
 .../general/SRF/SRFModel/SRFModel/SRFModel.C  |  9 ++-
 .../cfdTools/general/adjustPhi/adjustPhi.C    |  3 +-
 .../cfdTools/general/bound/bound.C            |  4 +-
 .../general/constrainHbyA/constrainHbyA.C     |  3 +-
 .../externalCoupledMixedFvPatchField.C        | 42 ++++++------
 ...mappedVelocityFluxFixedValueFvPatchField.C | 11 ++--
 .../waveSurfacePressureFvPatchScalarField.C   | 12 ++--
 .../CoEulerDdtScheme/CoEulerDdtScheme.C       | 12 ++--
 .../CrankNicolsonDdtScheme.C                  | 64 +++++++++++++------
 .../ddtSchemes/ddtScheme/ddtScheme.C          |  5 +-
 .../finiteVolume/fvc/fvcAverage.C             |  2 +-
 .../gradSchemes/fourthGrad/fourthGrad.C       |  2 +-
 .../gradSchemes/gaussGrad/gaussGrad.C         |  9 ++-
 .../leastSquaresGrad/leastSquaresGrad.C       |  2 +-
 .../leastSquaresGrad/leastSquaresVectors.C    | 10 +--
 .../faceCorrectedSnGrad/faceCorrectedSnGrad.C |  2 +-
 .../snGradSchemes/snGradScheme/snGradScheme.C |  8 ++-
 .../fvMatrices/fvMatrix/fvMatrix.C            | 10 +--
 .../fvMatrices/fvMatrix/fvMatrixSolve.C       |  2 +-
 .../fvScalarMatrix/fvScalarMatrix.C           |  2 +-
 .../solvers/MULES/CMULESTemplates.C           |  4 +-
 .../fvMatrices/solvers/MULES/MULESTemplates.C |  4 +-
 .../extendedCellToFaceStencilTemplates.C      |  2 +-
 ...extendedUpwindCellToFaceStencilTemplates.C |  2 +-
 src/finiteVolume/fvMesh/fvMesh.C              |  6 +-
 src/finiteVolume/fvMesh/fvMeshGeometry.C      |  7 +-
 .../fvMeshSubset/fvMeshSubsetInterpolate.C    |  6 +-
 .../Poisson/PoissonPatchDistMethod.C          |  2 +-
 .../advectionDiffusionPatchDistMethod.C       |  3 +-
 .../meshWave/meshWavePatchDistMethod.C        | 24 ++++---
 .../fvMesh/wallDist/wallDist/wallDist.C       |  4 +-
 .../LimitedScheme/LimitedScheme.C             |  2 +-
 .../limitedSchemes/PhiScheme/PhiScheme.C      |  2 +-
 .../limitedSurfaceInterpolationScheme.C       |  2 +-
 .../schemes/clippedLinear/clippedLinear.H     | 10 +--
 .../schemes/cubic/cubic.H                     |  9 ++-
 .../schemes/linearUpwind/linearUpwind.C       |  2 +-
 .../schemes/linearUpwind/linearUpwindV.C      |  2 +-
 .../schemes/localMax/localMax.H               |  7 +-
 .../schemes/localMin/localMin.H               |  7 +-
 .../schemes/midPoint/midPoint.H               | 10 +--
 .../schemes/pointLinear/pointLinear.C         |  2 +-
 .../schemes/reverseLinear/reverseLinear.H     | 12 ++--
 .../skewCorrected/skewCorrectionVectors.C     | 11 ++--
 .../surfaceInterpolation.C                    | 36 +++++++----
 .../surfaceInterpolationScheme.C              | 11 +++-
 .../pointConstraintsTemplates.C               |  7 +-
 .../volPointInterpolate.C                     | 15 +++--
 59 files changed, 369 insertions(+), 245 deletions(-)

diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
index 959a1182b1..b30ca6e5a7 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
@@ -712,7 +712,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::internalField()
 template<class Type, template<class> class PatchField, class GeoMesh>
 typename
 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricBoundaryField&
-Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryField()
+Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryFieldRef()
 {
     this->setUpToDate();
     storeOldTimes();
@@ -993,7 +993,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::T() const
     );
 
     Foam::T(result.ref().internalField(), internalField());
-    Foam::T(result.ref().boundaryField(), boundaryField());
+    Foam::T(result.ref().boundaryFieldRef(), boundaryField());
 
     return result;
 }
@@ -1030,7 +1030,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::component
     );
 
     Foam::component(Component.ref().internalField(), internalField(), d);
-    Foam::component(Component.ref().boundaryField(), boundaryField(), d);
+    Foam::component(Component.ref().boundaryFieldRef(), boundaryField(), d);
 
     return Component;
 }
@@ -1049,7 +1049,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::replace
 )
 {
     internalField().replace(d, gcf.internalField());
-    boundaryField().replace(d, gcf.boundaryField());
+    boundaryFieldRef().replace(d, gcf.boundaryField());
 }
 
 
@@ -1061,7 +1061,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::replace
 )
 {
     internalField().replace(d, ds.value());
-    boundaryField().replace(d, ds.value());
+    boundaryFieldRef().replace(d, ds.value());
 }
 
 
@@ -1072,7 +1072,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::max
 )
 {
     Foam::max(internalField(), internalField(), dt.value());
-    Foam::max(boundaryField(), boundaryField(), dt.value());
+    Foam::max(boundaryFieldRef(), boundaryField(), dt.value());
 }
 
 
@@ -1083,7 +1083,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::min
 )
 {
     Foam::min(internalField(), internalField(), dt.value());
-    Foam::min(boundaryField(), boundaryField(), dt.value());
+    Foam::min(boundaryFieldRef(), boundaryField(), dt.value());
 }
 
 
@@ -1091,7 +1091,7 @@ template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::negate()
 {
     internalField().negate();
-    boundaryField().negate();
+    boundaryFieldRef().negate();
 }
 
 
@@ -1115,7 +1115,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=
     // only equate field contents not ID
 
     dimensionedInternalField() = gf.dimensionedInternalField();
-    boundaryField() = gf.boundaryField();
+    boundaryFieldRef() = gf.boundaryField();
 }
 
 
@@ -1146,7 +1146,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=
         const_cast<Field<Type>&>(gf.internalField())
     );
 
-    boundaryField() = gf.boundaryField();
+    boundaryFieldRef() = gf.boundaryField();
 
     tgf.clear();
 }
@@ -1159,7 +1159,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=
 )
 {
     dimensionedInternalField() = dt;
-    boundaryField() = dt.value();
+    boundaryFieldRef() = dt.value();
 }
 
 
@@ -1176,7 +1176,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator==
     // only equate field contents not ID
 
     dimensionedInternalField() = gf.dimensionedInternalField();
-    boundaryField() == gf.boundaryField();
+    boundaryFieldRef() == gf.boundaryField();
 
     tgf.clear();
 }
@@ -1189,7 +1189,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator==
 )
 {
     dimensionedInternalField() = dt;
-    boundaryField() == dt.value();
+    boundaryFieldRef() == dt.value();
 }
 
 
@@ -1204,7 +1204,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator op              \
     checkField(*this, gf, #op);                                                \
                                                                                \
     dimensionedInternalField() op gf.dimensionedInternalField();               \
-    boundaryField() op gf.boundaryField();                                     \
+    boundaryFieldRef() op gf.boundaryField();                                  \
 }                                                                              \
                                                                                \
 template<class Type, template<class> class PatchField, class GeoMesh>          \
@@ -1224,7 +1224,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator op              \
 )                                                                              \
 {                                                                              \
     dimensionedInternalField() op dt;                                          \
-    boundaryField() op dt.value();                                             \
+    boundaryFieldRef() op dt.value();                                          \
 }
 
 COMPUTED_ASSIGNMENT(Type, +=)
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
index 4efc05d309..2f1b71f7f1 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
@@ -442,7 +442,15 @@ public:
         inline const InternalField& internalField() const;
 
         //- Return reference to GeometricBoundaryField
-        GeometricBoundaryField& boundaryField();
+        GeometricBoundaryField& boundaryFieldRef();
+
+        //- Return reference to GeometricBoundaryField
+        #ifndef BOUNDARY_FIELD_REF
+        GeometricBoundaryField& boundaryField()
+        {
+            return boundaryFieldRef();
+        }
+        #endif
 
         //- Return reference to GeometricBoundaryField for const field
         inline const GeometricBoundaryField& boundaryField() const;
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
index e04a1588b5..dda0738161 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
@@ -50,7 +50,7 @@ void component
 )
 {
     component(gcf.internalField(), gf.internalField(), d);
-    component(gcf.boundaryField(), gf.boundaryField(), d);
+    component(gcf.boundaryFieldRef(), gf.boundaryField(), d);
 }
 
 
@@ -62,7 +62,7 @@ void T
 )
 {
     T(gf.internalField(), gf1.internalField());
-    T(gf.boundaryField(), gf1.boundaryField());
+    T(gf.boundaryFieldRef(), gf1.boundaryField());
 }
 
 
@@ -80,7 +80,7 @@ void pow
 )
 {
     pow(gf.internalField(), gf1.internalField(), r);
-    pow(gf.boundaryField(), gf1.boundaryField(), r);
+    pow(gf.boundaryFieldRef(), gf1.boundaryField(), r);
 }
 
 template
@@ -174,7 +174,7 @@ void sqr
 )
 {
     sqr(gf.internalField(), gf1.internalField());
-    sqr(gf.boundaryField(), gf1.boundaryField());
+    sqr(gf.boundaryFieldRef(), gf1.boundaryField());
 }
 
 template<class Type, template<class> class PatchField, class GeoMesh>
@@ -262,7 +262,7 @@ void magSqr
 )
 {
     magSqr(gsf.internalField(), gf.internalField());
-    magSqr(gsf.boundaryField(), gf.boundaryField());
+    magSqr(gsf.boundaryFieldRef(), gf.boundaryField());
 }
 
 template<class Type, template<class> class PatchField, class GeoMesh>
@@ -334,7 +334,7 @@ void mag
 )
 {
     mag(gsf.internalField(), gf.internalField());
-    mag(gsf.boundaryField(), gf.boundaryField());
+    mag(gsf.boundaryFieldRef(), gf.boundaryField());
 }
 
 template<class Type, template<class> class PatchField, class GeoMesh>
@@ -411,7 +411,7 @@ void cmptAv
 )
 {
     cmptAv(gcf.internalField(), gf.internalField());
-    cmptAv(gcf.boundaryField(), gf.boundaryField());
+    cmptAv(gcf.boundaryFieldRef(), gf.boundaryField());
 }
 
 template<class Type, template<class> class PatchField, class GeoMesh>
@@ -600,7 +600,12 @@ void opFunc                                                                    \
 )                                                                              \
 {                                                                              \
     Foam::opFunc(gf.internalField(), gf1.internalField(), gf2.internalField());\
-    Foam::opFunc(gf.boundaryField(), gf1.boundaryField(), gf2.boundaryField());\
+    Foam::opFunc                                                               \
+    (                                                                          \
+        gf.boundaryFieldRef(),                                                 \
+        gf1.boundaryField(),                                                   \
+        gf2.boundaryField()                                                    \
+    );\
 }                                                                              \
                                                                                \
 template                                                                       \
@@ -746,7 +751,7 @@ void opFunc                                                                    \
 )                                                                              \
 {                                                                              \
     Foam::opFunc(gf.internalField(), gf1.internalField(), dvs.value());        \
-    Foam::opFunc(gf.boundaryField(), gf1.boundaryField(), dvs.value());        \
+    Foam::opFunc(gf.boundaryFieldRef(), gf1.boundaryField(), dvs.value());     \
 }                                                                              \
                                                                                \
 template                                                                       \
@@ -859,7 +864,7 @@ void opFunc                                                                    \
 )                                                                              \
 {                                                                              \
     Foam::opFunc(gf.internalField(), dvs.value(), gf1.internalField());        \
-    Foam::opFunc(gf.boundaryField(), dvs.value(), gf1.boundaryField());        \
+    Foam::opFunc(gf.boundaryFieldRef(), dvs.value(), gf1.boundaryField());     \
 }                                                                              \
                                                                                \
 template                                                                       \
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C
index 2d0b007c27..d044d677aa 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C
@@ -42,7 +42,7 @@ void Func                                                                      \
 )                                                                              \
 {                                                                              \
     Foam::Func(res.internalField(), gf1.internalField());                      \
-    Foam::Func(res.boundaryField(), gf1.boundaryField());                      \
+    Foam::Func(res.boundaryFieldRef(), gf1.boundaryField());                   \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -111,7 +111,7 @@ void OpFunc                                                                    \
 )                                                                              \
 {                                                                              \
     Foam::OpFunc(res.internalField(), gf1.internalField());                    \
-    Foam::OpFunc(res.boundaryField(), gf1.boundaryField());                    \
+    Foam::OpFunc(res.boundaryFieldRef(), gf1.boundaryField());                 \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -180,8 +180,13 @@ void Func                                                                      \
     const GeometricField<Type2, PatchField, GeoMesh>& gf2                      \
 )                                                                              \
 {                                                                              \
-    Foam::Func(res.internalField(), gf1.internalField(), gf2.internalField());\
-    Foam::Func(res.boundaryField(), gf1.boundaryField(), gf2.boundaryField());\
+    Foam::Func(res.internalField(), gf1.internalField(), gf2.internalField()); \
+    Foam::Func                                                                 \
+    (                                                                          \
+        res.boundaryFieldRef(),                                                \
+        gf1.boundaryField(),                                                   \
+        gf2.boundaryField()                                                    \
+    );\
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -309,7 +314,7 @@ void Func                                                                      \
 )                                                                              \
 {                                                                              \
     Foam::Func(res.internalField(), dt1.value(), gf2.internalField());         \
-    Foam::Func(res.boundaryField(), dt1.value(), gf2.boundaryField());         \
+    Foam::Func(res.boundaryFieldRef(), dt1.value(), gf2.boundaryField());      \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -400,7 +405,7 @@ void Func                                                                      \
 )                                                                              \
 {                                                                              \
     Foam::Func(res.internalField(), gf1.internalField(), dt2.value());         \
-    Foam::Func(res.boundaryField(), gf1.boundaryField(), dt2.value());         \
+    Foam::Func(res.boundaryFieldRef(), gf1.boundaryField(), dt2.value());      \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -500,7 +505,7 @@ void OpFunc                                                                    \
     Foam::OpFunc                                                               \
         (res.internalField(), gf1.internalField(), gf2.internalField());       \
     Foam::OpFunc                                                               \
-        (res.boundaryField(), gf1.boundaryField(), gf2.boundaryField());       \
+        (res.boundaryFieldRef(), gf1.boundaryField(), gf2.boundaryField());    \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -628,7 +633,7 @@ void OpFunc                                                                    \
 )                                                                              \
 {                                                                              \
     Foam::OpFunc(res.internalField(), dt1.value(), gf2.internalField());       \
-    Foam::OpFunc(res.boundaryField(), dt1.value(), gf2.boundaryField());       \
+    Foam::OpFunc(res.boundaryFieldRef(), dt1.value(), gf2.boundaryField());    \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -719,7 +724,7 @@ void OpFunc                                                                    \
 )                                                                              \
 {                                                                              \
     Foam::OpFunc(res.internalField(), gf1.internalField(), dt2.value());       \
-    Foam::OpFunc(res.boundaryField(), gf1.boundaryField(), dt2.value());       \
+    Foam::OpFunc(res.boundaryFieldRef(), gf1.boundaryField(), dt2.value());    \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H
index 637e4314f3..ca437acec4 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H
@@ -132,7 +132,7 @@ void MapGeometricFields
 
             // Map the patch fields
             typename GeometricField<Type, PatchField, GeoMesh>
-            ::GeometricBoundaryField& bfield = field.boundaryField();
+            ::GeometricBoundaryField& bfield = field.boundaryFieldRef();
             forAll(bfield, patchi)
             {
                 // Cannot check sizes for patch fields because of
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C b/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
index 41bbe8178f..6b58138dc7 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
@@ -44,7 +44,7 @@ void stabilise
 )
 {
     stabilise(result.internalField(), gsf.internalField(), ds.value());
-    stabilise(result.boundaryField(), gsf.boundaryField(), ds.value());
+    stabilise(result.boundaryFieldRef(), gsf.boundaryField(), ds.value());
 }
 
 
@@ -127,7 +127,7 @@ void pow
 )
 {
     pow(Pow.internalField(), gsf1.internalField(), gsf2.internalField());
-    pow(Pow.boundaryField(), gsf1.boundaryField(), gsf2.boundaryField());
+    pow(Pow.boundaryFieldRef(), gsf1.boundaryField(), gsf2.boundaryField());
 }
 
 
@@ -270,7 +270,7 @@ void pow
 )
 {
     pow(tPow.internalField(), gsf.internalField(), ds.value());
-    pow(tPow.boundaryField(), gsf.boundaryField(), ds.value());
+    pow(tPow.boundaryFieldRef(), gsf.boundaryField(), ds.value());
 }
 
 
@@ -359,7 +359,7 @@ void pow
 )
 {
     pow(tPow.internalField(), ds.value(), gsf.internalField());
-    pow(tPow.boundaryField(), ds.value(), gsf.boundaryField());
+    pow(tPow.boundaryFieldRef(), ds.value(), gsf.boundaryField());
 }
 
 
@@ -451,7 +451,7 @@ void atan2
 )
 {
     atan2(Atan2.internalField(), gsf1.internalField(), gsf2.internalField());
-    atan2(Atan2.boundaryField(), gsf1.boundaryField(), gsf2.boundaryField());
+    atan2(Atan2.boundaryFieldRef(), gsf1.boundaryField(), gsf2.boundaryField());
 }
 
 
@@ -578,7 +578,7 @@ void atan2
 )
 {
     atan2(tAtan2.internalField(), gsf.internalField(), ds.value());
-    atan2(tAtan2.boundaryField(), gsf.boundaryField(), ds.value());
+    atan2(tAtan2.boundaryFieldRef(), gsf.boundaryField(), ds.value());
 }
 
 
@@ -667,7 +667,7 @@ void atan2
 )
 {
     atan2(tAtan2.internalField(), ds.value(), gsf.internalField());
-    atan2(tAtan2.boundaryField(), ds.value(), gsf.boundaryField());
+    atan2(tAtan2.boundaryFieldRef(), ds.value(), gsf.boundaryField());
 }
 
 
@@ -800,7 +800,7 @@ void func                                                                      \
 )                                                                              \
 {                                                                              \
     func(gsf.internalField(), n, gsf1.internalField());                        \
-    func(gsf.boundaryField(), n, gsf1.boundaryField());                        \
+    func(gsf.boundaryFieldRef(), n, gsf1.boundaryField());                     \
 }                                                                              \
                                                                                \
 template<template<class> class PatchField, class GeoMesh>                      \
diff --git a/src/OpenFOAM/fields/GeometricFields/transformGeometricField/transformGeometricField.C b/src/OpenFOAM/fields/GeometricFields/transformGeometricField/transformGeometricField.C
index e3b75000b4..baa2f7041b 100644
--- a/src/OpenFOAM/fields/GeometricFields/transformGeometricField/transformGeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/transformGeometricField/transformGeometricField.C
@@ -46,7 +46,7 @@ void transform
 )
 {
     transform(rtf.internalField(), trf.internalField(), tf.internalField());
-    transform(rtf.boundaryField(), trf.boundaryField(), tf.boundaryField());
+    transform(rtf.boundaryFieldRef(), trf.boundaryField(), tf.boundaryField());
 }
 
 
@@ -132,7 +132,7 @@ void transform
 )
 {
     transform(rtf.internalField(), t.value(), tf.internalField());
-    transform(rtf.boundaryField(), t.value(), tf.boundaryField());
+    transform(rtf.boundaryFieldRef(), t.value(), tf.boundaryField());
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/CorrectPhi/correctUphiBCs.C b/src/finiteVolume/cfdTools/general/CorrectPhi/correctUphiBCs.C
index d9ab448b46..beff14e462 100644
--- a/src/finiteVolume/cfdTools/general/CorrectPhi/correctUphiBCs.C
+++ b/src/finiteVolume/cfdTools/general/CorrectPhi/correctUphiBCs.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,23 +37,25 @@ void Foam::correctUphiBCs
 
     if (mesh.changing())
     {
-        forAll(U.boundaryField(), patchi)
+        volVectorField::GeometricBoundaryField& Ubf = U.boundaryFieldRef();
+        surfaceScalarField::GeometricBoundaryField& phibf =
+            phi.boundaryFieldRef();
+
+        forAll(Ubf, patchi)
         {
-            if (U.boundaryField()[patchi].fixesValue())
+            if (Ubf[patchi].fixesValue())
             {
-                U.boundaryField()[patchi].initEvaluate();
+                Ubf[patchi].initEvaluate();
             }
         }
 
-        forAll(U.boundaryField(), patchi)
+        forAll(Ubf, patchi)
         {
-            if (U.boundaryField()[patchi].fixesValue())
+            if (Ubf[patchi].fixesValue())
             {
-                U.boundaryField()[patchi].evaluate();
+                Ubf[patchi].evaluate();
 
-                phi.boundaryField()[patchi] =
-                    U.boundaryField()[patchi]
-                  & mesh.Sf().boundaryField()[patchi];
+                phibf[patchi] = Ubf[patchi] & mesh.Sf().boundaryField()[patchi];
             }
         }
     }
@@ -71,24 +73,28 @@ void Foam::correctUphiBCs
 
     if (mesh.changing())
     {
-        forAll(U.boundaryField(), patchi)
+        volVectorField::GeometricBoundaryField& Ubf = U.boundaryFieldRef();
+        surfaceScalarField::GeometricBoundaryField& phibf =
+            phi.boundaryFieldRef();
+
+        forAll(Ubf, patchi)
         {
-            if (U.boundaryField()[patchi].fixesValue())
+            if (Ubf[patchi].fixesValue())
             {
-                U.boundaryField()[patchi].initEvaluate();
+                Ubf[patchi].initEvaluate();
             }
         }
 
-        forAll(U.boundaryField(), patchi)
+        forAll(Ubf, patchi)
         {
-            if (U.boundaryField()[patchi].fixesValue())
+            if (Ubf[patchi].fixesValue())
             {
-                U.boundaryField()[patchi].evaluate();
+                Ubf[patchi].evaluate();
 
-                phi.boundaryField()[patchi] =
+                phibf[patchi] =
                     rho.boundaryField()[patchi]
                    *(
-                        U.boundaryField()[patchi]
+                        Ubf[patchi]
                       & mesh.Sf().boundaryField()[patchi]
                     );
             }
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
index 90239b4984..00fb6ffe09 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
@@ -418,12 +418,15 @@ void Foam::MRFZone::makeRelative(volVectorField& U) const
     }
 
     // Included patches
+
+    volVectorField::GeometricBoundaryField& Ubf = U.boundaryFieldRef();
+
     forAll(includedFaces_, patchi)
     {
         forAll(includedFaces_[patchi], i)
         {
             label patchFacei = includedFaces_[patchi][i];
-            U.boundaryField()[patchi][patchFacei] = Zero;
+            Ubf[patchi][patchFacei] = Zero;
         }
     }
 
@@ -433,7 +436,7 @@ void Foam::MRFZone::makeRelative(volVectorField& U) const
         forAll(excludedFaces_[patchi], i)
         {
             label patchFacei = excludedFaces_[patchi][i];
-            U.boundaryField()[patchi][patchFacei] -=
+            Ubf[patchi][patchFacei] -=
                 (Omega
               ^ (C.boundaryField()[patchi][patchFacei] - origin_));
         }
@@ -484,12 +487,14 @@ void Foam::MRFZone::makeAbsolute(volVectorField& U) const
     }
 
     // Included patches
+    volVectorField::GeometricBoundaryField& Ubf = U.boundaryFieldRef();
+
     forAll(includedFaces_, patchi)
     {
         forAll(includedFaces_[patchi], i)
         {
             label patchFacei = includedFaces_[patchi][i];
-            U.boundaryField()[patchi][patchFacei] =
+            Ubf[patchi][patchFacei] =
                 (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin_));
         }
     }
@@ -500,7 +505,7 @@ void Foam::MRFZone::makeAbsolute(volVectorField& U) const
         forAll(excludedFaces_[patchi], i)
         {
             label patchFacei = excludedFaces_[patchi][i];
-            U.boundaryField()[patchi][patchFacei] +=
+            Ubf[patchi][patchFacei] +=
                 (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin_));
         }
     }
@@ -528,11 +533,13 @@ void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
     const vector Omega = this->Omega();
 
     // Included patches
+    volVectorField::GeometricBoundaryField& Ubf = U.boundaryFieldRef();
+
     forAll(includedFaces_, patchi)
     {
         const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
 
-        vectorField pfld(U.boundaryField()[patchi]);
+        vectorField pfld(Ubf[patchi]);
 
         forAll(includedFaces_[patchi], i)
         {
@@ -541,7 +548,7 @@ void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
             pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin_));
         }
 
-        U.boundaryField()[patchi] == pfld;
+        Ubf[patchi] == pfld;
     }
 }
 
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C
index 271157b9a6..761e7b7428 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C
@@ -393,19 +393,22 @@ void Foam::MRFZoneList::correctBoundaryFlux
     surfaceScalarField& phi
 ) const
 {
-    FieldField<fvsPatchField, scalar> phibf
+    FieldField<fvsPatchField, scalar> Uf
     (
         relative(mesh_.Sf().boundaryField() & U.boundaryField())
     );
 
+
+    surfaceScalarField::GeometricBoundaryField& phibf = phi.boundaryFieldRef();
+
     forAll(mesh_.boundary(), patchi)
     {
         if
         (
-            isA<fixedValueFvsPatchScalarField>(phi.boundaryField()[patchi])
+            isA<fixedValueFvsPatchScalarField>(phibf[patchi])
         )
         {
-            phi.boundaryField()[patchi] == phibf[patchi];
+            phibf[patchi] == Uf[patchi];
         }
     }
 }
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
index b4a59329df..effeed232c 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
@@ -54,7 +54,7 @@ void Foam::MRFZone::makeRelativeRhoFlux
         phii[facei] -= rho[facei]*(Omega ^ (Cfi[facei] - origin_)) & Sfi[facei];
     }
 
-    makeRelativeRhoFlux(rho.boundaryField(), phi.boundaryField());
+    makeRelativeRhoFlux(rho.boundaryField(), phi.boundaryFieldRef());
 }
 
 
@@ -154,6 +154,9 @@ void Foam::MRFZone::makeAbsoluteRhoFlux
         phii[facei] += rho[facei]*(Omega ^ (Cfi[facei] - origin_)) & Sfi[facei];
     }
 
+    surfaceScalarField::GeometricBoundaryField& phibf = phi.boundaryFieldRef();
+
+
     // Included patches
     forAll(includedFaces_, patchi)
     {
@@ -161,7 +164,7 @@ void Foam::MRFZone::makeAbsoluteRhoFlux
         {
             label patchFacei = includedFaces_[patchi][i];
 
-            phi.boundaryField()[patchi][patchFacei] +=
+            phibf[patchi][patchFacei] +=
                 rho.boundaryField()[patchi][patchFacei]
               * (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin_))
               & Sf.boundaryField()[patchi][patchFacei];
@@ -175,7 +178,7 @@ void Foam::MRFZone::makeAbsoluteRhoFlux
         {
             label patchFacei = excludedFaces_[patchi][i];
 
-            phi.boundaryField()[patchi][patchFacei] +=
+            phibf[patchi][patchFacei] +=
                 rho.boundaryField()[patchi][patchFacei]
               * (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin_))
               & Sf.boundaryField()[patchi][patchFacei];
diff --git a/src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C b/src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C
index 7077bc2f6b..a48b028b5e 100644
--- a/src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C
+++ b/src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C
@@ -225,10 +225,13 @@ Foam::tmp<Foam::volVectorField> Foam::SRF::SRFModel::Uabs() const
         )
     );
 
+    volVectorField& Uabs = tUabs.ref();
+
     // Add SRF contribution to internal field
-    tUabs.ref().internalField() += Urel_.internalField();
+    Uabs.internalField() += Urel_.internalField();
 
     // Add Urel boundary contributions
+    volVectorField::GeometricBoundaryField& Uabsbf = Uabs.boundaryFieldRef();
     const volVectorField::GeometricBoundaryField& bvf = Urel_.boundaryField();
 
     forAll(bvf, i)
@@ -241,12 +244,12 @@ Foam::tmp<Foam::volVectorField> Foam::SRF::SRFModel::Uabs() const
                 refCast<const SRFVelocityFvPatchVectorField>(bvf[i]);
             if (UrelPatch.relative())
             {
-                tUabs.ref().boundaryField()[i] += Urel_.boundaryField()[i];
+                Uabsbf[i] += Urel_.boundaryField()[i];
             }
         }
         else
         {
-            tUabs.ref().boundaryField()[i] += Urel_.boundaryField()[i];
+            Uabsbf[i] += Urel_.boundaryField()[i];
         }
     }
 
diff --git a/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C b/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
index 6f0742b26e..c34239a4c3 100644
--- a/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
+++ b/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
@@ -43,7 +43,8 @@ bool Foam::adjustPhi
         scalar fixedMassOut = 0.0;
         scalar adjustableMassOut = 0.0;
 
-        surfaceScalarField::GeometricBoundaryField& bphi = phi.boundaryField();
+        surfaceScalarField::GeometricBoundaryField& bphi =
+            phi.boundaryFieldRef();
 
         forAll(bphi, patchi)
         {
diff --git a/src/finiteVolume/cfdTools/general/bound/bound.C b/src/finiteVolume/cfdTools/general/bound/bound.C
index f47cfd33ff..009a14388d 100644
--- a/src/finiteVolume/cfdTools/general/bound/bound.C
+++ b/src/finiteVolume/cfdTools/general/bound/bound.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,7 +53,7 @@ Foam::bound(volScalarField& vsf, const dimensionedScalar& lowerBound)
             lowerBound.value()
         );
 
-        vsf.boundaryField() = max(vsf.boundaryField(), lowerBound.value());
+        vsf.boundaryFieldRef() = max(vsf.boundaryField(), lowerBound.value());
     }
 
     return vsf;
diff --git a/src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.C b/src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.C
index 38f84f96a7..bd51322e6e 100644
--- a/src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.C
+++ b/src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.C
@@ -49,6 +49,7 @@ Foam::tmp<Foam::volVectorField> Foam::constrainHbyA
     }
 
     volVectorField& HbyA = tHbyANew.ref();
+    volVectorField::GeometricBoundaryField& HbyAbf = HbyA.boundaryFieldRef();
 
     forAll(U.boundaryField(), patchi)
     {
@@ -61,7 +62,7 @@ Foam::tmp<Foam::volVectorField> Foam::constrainHbyA
             )
         )
         {
-            HbyA.boundaryField()[patchi] = U.boundaryField()[patchi];
+            HbyAbf[patchi] = U.boundaryField()[patchi];
         }
     }
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/externalCoupledMixed/externalCoupledMixedFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/externalCoupledMixed/externalCoupledMixedFvPatchField.C
index b3629e877f..568bc9a6f7 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/externalCoupledMixed/externalCoupledMixedFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/externalCoupledMixed/externalCoupledMixedFvPatchField.C
@@ -68,7 +68,7 @@ void Foam::externalCoupledMixedFvPatchField<Type>::setMaster
 
     volFieldType& vf = const_cast<volFieldType&>(cvf);
 
-    typename volFieldType::GeometricBoundaryField& bf = vf.boundaryField();
+    typename volFieldType::GeometricBoundaryField& bf = vf.boundaryFieldRef();
 
     // number of patches can be different in parallel...
     label nPatch = bf.size();
@@ -84,11 +84,11 @@ void Foam::externalCoupledMixedFvPatchField<Type>::setMaster
     // set the master patch
     forAll(patchIDs, i)
     {
-        label patchI = patchIDs[i];
+        label patchi = patchIDs[i];
 
-        patchType& pf = refCast<patchType>(bf[patchI]);
+        patchType& pf = refCast<patchType>(bf[patchi]);
 
-        offsets_[patchI][Pstream::myProcNo()] = pf.size();
+        offsets_[patchi][Pstream::myProcNo()] = pf.size();
 
         if (i == 0)
         {
@@ -109,10 +109,10 @@ void Foam::externalCoupledMixedFvPatchField<Type>::setMaster
     }
 
     label patchOffset = 0;
-    forAll(offsets_, patchI)
+    forAll(offsets_, patchi)
     {
         label sumOffset = 0;
-        List<label>& procOffsets = offsets_[patchI];
+        List<label>& procOffsets = offsets_[patchi];
 
         forAll(procOffsets, procI)
         {
@@ -249,9 +249,9 @@ void Foam::externalCoupledMixedFvPatchField<Type>::startWait() const
 
     forAll(coupledPatchIDs_, i)
     {
-        label patchI = coupledPatchIDs_[i];
+        label patchi = coupledPatchIDs_[i];
 
-        const patchType& pf = refCast<const patchType>(bf[patchI]);
+        const patchType& pf = refCast<const patchType>(bf[patchi]);
 
         if (pf.master())
         {
@@ -413,9 +413,9 @@ void Foam::externalCoupledMixedFvPatchField<Type>::writeData
 
     forAll(coupledPatchIDs_, i)
     {
-        label patchI = coupledPatchIDs_[i];
+        label patchi = coupledPatchIDs_[i];
 
-        const patchType& pf = refCast<const patchType>(bf[patchI]);
+        const patchType& pf = refCast<const patchType>(bf[patchi]);
 
         pf.transferData(os);
     }
@@ -612,16 +612,16 @@ void Foam::externalCoupledMixedFvPatchField<Type>::initialise
 
     volFieldType& vf = const_cast<volFieldType&>(cvf);
 
-    typename volFieldType::GeometricBoundaryField& bf = vf.boundaryField();
+    typename volFieldType::GeometricBoundaryField& bf = vf.boundaryFieldRef();
 
     // identify all coupled patches
     DynamicList<label> coupledPatchIDs(bf.size());
 
-    forAll(bf, patchI)
+    forAll(bf, patchi)
     {
-        if (isA<patchType>(bf[patchI]))
+        if (isA<patchType>(bf[patchi]))
         {
-            coupledPatchIDs.append(patchI);
+            coupledPatchIDs.append(patchi);
         }
     }
 
@@ -636,9 +636,9 @@ void Foam::externalCoupledMixedFvPatchField<Type>::initialise
 
         forAll(coupledPatchIDs_, i)
         {
-            label patchI = coupledPatchIDs_[i];
+            label patchi = coupledPatchIDs_[i];
 
-            patchType& pf = refCast<patchType>(bf[patchI]);
+            patchType& pf = refCast<patchType>(bf[patchi]);
 
             pf.setMaster(coupledPatchIDs_);
         }
@@ -652,9 +652,9 @@ void Foam::externalCoupledMixedFvPatchField<Type>::initialise
         {
             forAll(coupledPatchIDs_, i)
             {
-                label patchI = coupledPatchIDs_[i];
+                label patchi = coupledPatchIDs_[i];
 
-                patchType& pf = refCast<patchType>(bf[patchI]);
+                patchType& pf = refCast<patchType>(bf[patchi]);
 
                 pf.readData(transferFile);
             }
@@ -791,11 +791,11 @@ void Foam::externalCoupledMixedFvPatchField<Type>::writeGeometry() const
             << "writing collated patch faces to: " << osFaces.name() << endl;
     }
 
-    forAll(bf, patchI)
+    forAll(bf, patchi)
     {
-        if (isA<patchType>(bf[patchI]))
+        if (isA<patchType>(bf[patchi]))
         {
-            const patchType& pf = refCast<const patchType>(bf[patchI]);
+            const patchType& pf = refCast<const patchType>(bf[patchi]);
 
             pf.writeGeometry(osPoints, osFaces);
         }
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedVelocityFluxFixedValue/mappedVelocityFluxFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/mappedVelocityFluxFixedValue/mappedVelocityFluxFixedValueFvPatchField.C
index be3e6304b9..56a87b9805 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/mappedVelocityFluxFixedValue/mappedVelocityFluxFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedVelocityFluxFixedValue/mappedVelocityFluxFixedValueFvPatchField.C
@@ -156,10 +156,8 @@ void Foam::mappedVelocityFluxFixedValueFvPatchField::updateCoeffs()
     const volVectorField& UField =
         nbrMesh.lookupObject<volVectorField>(fieldName);
 
-    surfaceScalarField& phiField = const_cast<surfaceScalarField&>
-    (
-        nbrMesh.lookupObject<surfaceScalarField>(phiName_)
-    );
+    const surfaceScalarField& phiField =
+        nbrMesh.lookupObject<surfaceScalarField>(phiName_);
 
     vectorField newUValues;
     scalarField newPhiValues;
@@ -217,7 +215,10 @@ void Foam::mappedVelocityFluxFixedValueFvPatchField::updateCoeffs()
     }
 
     operator==(newUValues);
-    phiField.boundaryField()[patch().index()] == newPhiValues;
+    const_cast<surfaceScalarField&>
+    (
+        phiField
+    ).boundaryFieldRef()[patch().index()] == newPhiValues;
 
     // Restore tag
     UPstream::msgType() = oldTag;
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
index 04d46c8eed..4c9acd7f21 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -156,7 +156,7 @@ void Foam::waveSurfacePressureFvPatchScalarField::updateCoeffs()
         (
             db().lookupObject<volVectorField>(zetaName_)
         );
-    vectorField& zetap = zeta.boundaryField()[patchI];
+    vectorField& zetap = zeta.boundaryFieldRef()[patchI];
 
     // lookup d/dt scheme from database for zeta
     const word ddtSchemeName(zeta.mesh().ddtScheme(zeta.name()));
@@ -180,12 +180,14 @@ void Foam::waveSurfacePressureFvPatchScalarField::updateCoeffs()
         dZetap /= rhop;
     }
 
+    const volVectorField& zeta0 = zeta.oldTime();
+
     switch (ddtScheme)
     {
         case tsEuler:
         case tsCrankNicolson:
         {
-            zetap = zeta.oldTime().boundaryField()[patchI] + dZetap;
+            zetap = zeta0.boundaryField()[patchI] + dZetap;
 
             break;
         }
@@ -199,8 +201,8 @@ void Foam::waveSurfacePressureFvPatchScalarField::updateCoeffs()
 
             zetap =
                 (
-                    c0*zeta.oldTime().boundaryField()[patchI]
-                  - c00*zeta.oldTime().oldTime().boundaryField()[patchI]
+                    c0*zeta0.boundaryField()[patchI]
+                  - c00*zeta0.oldTime().boundaryField()[patchI]
                   + dZetap
                 )/c;
 
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
index 4165b62191..55e5c0c391 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
@@ -75,14 +75,12 @@ tmp<volScalarField> CoEulerDdtScheme<Type>::CorDeltaT() const
             max(corDeltaT[neighbour[faceI]], cofrDeltaT[faceI]);
     }
 
-    volScalarField::GeometricBoundaryField& bcorDeltaT =
-        corDeltaT.boundaryField();
+    const surfaceScalarField::GeometricBoundaryField& cofrDeltaTbf =
+        cofrDeltaT.boundaryField();
 
-    forAll(bcorDeltaT, patchi)
+    forAll(cofrDeltaTbf, patchi)
     {
-        const fvsPatchScalarField& pcofrDeltaT =
-            cofrDeltaT.boundaryField()[patchi];
-
+        const fvsPatchScalarField& pcofrDeltaT = cofrDeltaTbf[patchi];
         const fvPatch& p = pcofrDeltaT.patch();
         const labelUList& faceCells = p.patch().faceCells();
 
@@ -98,8 +96,6 @@ tmp<volScalarField> CoEulerDdtScheme<Type>::CorDeltaT() const
 
     corDeltaT.correctBoundaryConditions();
 
-    //corDeltaT = max(corDeltaT, max(corDeltaT)/100.0);
-
     return tcorDeltaT;
 }
 
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
index 67473222f5..b424943e34 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
@@ -366,6 +366,9 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
+        typename DDt0Field<GeometricField<Type, fvPatchField, volMesh>>::
+            GeometricBoundaryField& ddt0bf = ddt0.boundaryFieldRef();
+
         if (evaluate(ddt0))
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
@@ -379,13 +382,13 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                 ) - mesh().V00()*offCentre_(ddt0.internalField())
             )/mesh().V0();
 
-            ddt0.boundaryField() =
+            ddt0bf =
             (
                 rDtCoef0*
                 (
                     vf.oldTime().boundaryField()
                   - vf.oldTime().oldTime().boundaryField()
-                ) - offCentre_(ff(ddt0.boundaryField()))
+                ) - offCentre_(ff(ddt0bf))
             );
         }
 
@@ -406,7 +409,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                 rDtCoef.value()*
                 (
                     vf.boundaryField() - vf.oldTime().boundaryField()
-                ) - offCentre_(ff(ddt0.boundaryField()))
+                ) - offCentre_(ff(ddt0bf))
             )
         );
     }
@@ -456,6 +459,9 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
+        typename DDt0Field<GeometricField<Type, fvPatchField, volMesh>>::
+            GeometricBoundaryField& ddt0bf = ddt0.boundaryFieldRef();
+
         if (evaluate(ddt0))
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
@@ -469,13 +475,13 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                 ) - mesh().V00()*offCentre_(ddt0.internalField())
             )/mesh().V0();
 
-            ddt0.boundaryField() =
+            ddt0bf =
             (
                 rDtCoef0*rho.value()*
                 (
                     vf.oldTime().boundaryField()
                   - vf.oldTime().oldTime().boundaryField()
-                ) - offCentre_(ff(ddt0.boundaryField()))
+                ) - offCentre_(ff(ddt0bf))
             );
         }
 
@@ -496,7 +502,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                 rDtCoef.value()*rho.value()*
                 (
                     vf.boundaryField() - vf.oldTime().boundaryField()
-                ) - offCentre_(ff(ddt0.boundaryField()))
+                ) - offCentre_(ff(ddt0bf))
             )
         );
     }
@@ -546,6 +552,9 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
+        typename DDt0Field<GeometricField<Type, fvPatchField, volMesh>>::
+            GeometricBoundaryField& ddt0bf = ddt0.boundaryFieldRef();
+
         if (evaluate(ddt0))
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
@@ -561,7 +570,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                 ) - mesh().V00()*offCentre_(ddt0.internalField())
             )/mesh().V0();
 
-            ddt0.boundaryField() =
+            ddt0bf =
             (
                 rDtCoef0*
                 (
@@ -569,7 +578,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                    *vf.oldTime().boundaryField()
                   - rho.oldTime().oldTime().boundaryField()
                    *vf.oldTime().oldTime().boundaryField()
-                ) - offCentre_(ff(ddt0.boundaryField()))
+                ) - offCentre_(ff(ddt0bf))
             );
         }
 
@@ -592,7 +601,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                 (
                     rho.boundaryField()*vf.boundaryField()
                   - rho.oldTime().boundaryField()*vf.oldTime().boundaryField()
-                ) - offCentre_(ff(ddt0.boundaryField()))
+                ) - offCentre_(ff(ddt0bf))
             )
         );
     }
@@ -647,6 +656,9 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
+        typename DDt0Field<GeometricField<Type, fvPatchField, volMesh>>::
+            GeometricBoundaryField& ddt0bf = ddt0.boundaryFieldRef();
+
         if (evaluate(ddt0))
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
@@ -667,7 +679,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                 ) - mesh().V00()*offCentre_(ddt0.internalField())
             )/mesh().V0();
 
-            ddt0.boundaryField() =
+            ddt0bf =
             (
                 rDtCoef0*
                 (
@@ -678,7 +690,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                   - alpha.oldTime().oldTime().boundaryField()
                    *rho.oldTime().oldTime().boundaryField()
                    *vf.oldTime().oldTime().boundaryField()
-                ) - offCentre_(ff(ddt0.boundaryField()))
+                ) - offCentre_(ff(ddt0bf))
             );
         }
 
@@ -713,7 +725,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                   - alpha.oldTime().boundaryField()
                    *rho.oldTime().boundaryField()
                    *vf.oldTime().boundaryField()
-                ) - offCentre_(ff(ddt0.boundaryField()))
+                ) - offCentre_(ff(ddt0bf))
             )
         );
     }
@@ -782,6 +794,9 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
 
     if (mesh().moving())
     {
+        typename DDt0Field<GeometricField<Type, fvPatchField, volMesh>>::
+            GeometricBoundaryField& ddt0bf = ddt0.boundaryFieldRef();
+
         if (evaluate(ddt0))
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
@@ -796,14 +811,14 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
               - mesh().V00()*offCentre_(ddt0.internalField())
             )/mesh().V0();
 
-            ddt0.boundaryField() =
+            ddt0bf =
             (
                 rDtCoef0*
                 (
                     vf.oldTime().boundaryField()
                   - vf.oldTime().oldTime().boundaryField()
                 )
-              - offCentre_(ff(ddt0.boundaryField()))
+              - offCentre_(ff(ddt0bf))
             );
         }
 
@@ -864,6 +879,9 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
 
     if (mesh().moving())
     {
+        typename DDt0Field<GeometricField<Type, fvPatchField, volMesh>>::
+            GeometricBoundaryField& ddt0bf = ddt0.boundaryFieldRef();
+
         if (evaluate(ddt0))
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
@@ -878,14 +896,14 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
               - mesh().V00()*offCentre_(ddt0.internalField())
             )/mesh().V0();
 
-            ddt0.boundaryField() =
+            ddt0bf =
             (
                 rDtCoef0*rho.value()*
                 (
                     vf.oldTime().boundaryField()
                   - vf.oldTime().oldTime().boundaryField()
                 )
-              - offCentre_(ff(ddt0.boundaryField()))
+              - offCentre_(ff(ddt0bf))
             );
         }
 
@@ -947,6 +965,9 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
 
     if (mesh().moving())
     {
+        typename DDt0Field<GeometricField<Type, fvPatchField, volMesh>>::
+            GeometricBoundaryField& ddt0bf = ddt0.boundaryFieldRef();
+
         if (evaluate(ddt0))
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
@@ -963,7 +984,7 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
               - mesh().V00()*offCentre_(ddt0.internalField())
             )/mesh().V0();
 
-            ddt0.boundaryField() =
+            ddt0bf =
             (
                 rDtCoef0*
                 (
@@ -972,7 +993,7 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
                   - rho.oldTime().oldTime().boundaryField()
                    *vf.oldTime().oldTime().boundaryField()
                 )
-              - offCentre_(ff(ddt0.boundaryField()))
+              - offCentre_(ff(ddt0bf))
             );
         }
 
@@ -1039,6 +1060,9 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
 
     if (mesh().moving())
     {
+        typename DDt0Field<GeometricField<Type, fvPatchField, volMesh>>::
+            GeometricBoundaryField& ddt0bf = ddt0.boundaryFieldRef();
+
         if (evaluate(ddt0))
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
@@ -1060,7 +1084,7 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
               - mesh().V00()*offCentre_(ddt0.internalField())
             )/mesh().V0();
 
-            ddt0.boundaryField() =
+            ddt0bf =
             (
                 rDtCoef0*
                 (
@@ -1072,7 +1096,7 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
                    *rho.oldTime().oldTime().boundaryField()
                    *vf.oldTime().oldTime().boundaryField()
                 )
-              - offCentre_(ff(ddt0.boundaryField()))
+              - offCentre_(ff(ddt0bf))
             );
         }
 
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C
index cd0e1adc7e..d83c832a66 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C
@@ -150,6 +150,9 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
 
     surfaceScalarField& ddtCouplingCoeff = tddtCouplingCoeff.ref();
 
+    surfaceScalarField::GeometricBoundaryField& ccbf =
+        ddtCouplingCoeff.boundaryFieldRef();
+
     forAll(U.boundaryField(), patchi)
     {
         if
@@ -158,7 +161,7 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
          || isA<cyclicAMIFvPatch>(mesh().boundary()[patchi])
         )
         {
-            ddtCouplingCoeff.boundaryField()[patchi] = 0.0;
+            ccbf[patchi] = 0.0;
         }
     }
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcAverage.C b/src/finiteVolume/finiteVolume/fvc/fvcAverage.C
index 711a012f66..2ea76c1dc1 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcAverage.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcAverage.C
@@ -74,7 +74,7 @@ average
     );
 
     typename GeometricField<Type, fvPatchField, volMesh>::
-    GeometricBoundaryField& bav = av.boundaryField();
+    GeometricBoundaryField& bav = av.boundaryFieldRef();
 
     forAll(bav, patchi)
     {
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C
index 975ff9141d..367c33d673 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C
@@ -125,7 +125,7 @@ Foam::fv::fourthGrad<Type>::calcGrad
 
             const scalarField& lambdap = lambda.boundaryField()[patchi];
 
-            const fvPatch& p = fGrad.boundaryField()[patchi].patch();
+            const fvPatch& p = vsf.boundaryField()[patchi].patch();
 
             const labelUList& faceCells = p.faceCells();
 
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C
index 34b9170c62..1bba801db6 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C
@@ -150,6 +150,11 @@ void Foam::fv::gaussGrad<Type>::correctBoundaryConditions
     >& gGrad
 )
 {
+    typename GeometricField
+    <
+        typename outerProduct<vector, Type>::type, fvPatchField, volMesh
+    >::GeometricBoundaryField& gGradbf = gGrad.boundaryFieldRef();
+
     forAll(vsf.boundaryField(), patchi)
     {
         if (!vsf.boundaryField()[patchi].coupled())
@@ -160,10 +165,10 @@ void Foam::fv::gaussGrad<Type>::correctBoundaryConditions
               / vsf.mesh().magSf().boundaryField()[patchi]
             );
 
-            gGrad.boundaryField()[patchi] += n *
+            gGradbf[patchi] += n *
             (
                 vsf.boundaryField()[patchi].snGrad()
-              - (n & gGrad.boundaryField()[patchi])
+              - (n & gGradbf[patchi])
             );
         }
      }
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C
index 2a3839b5b2..d027065b4c 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C
@@ -104,7 +104,7 @@ Foam::fv::leastSquaresGrad<Type>::calcGrad
         const fvsPatchVectorField& patchOwnLs = ownLs.boundaryField()[patchi];
 
         const labelUList& faceCells =
-            lsGrad.boundaryField()[patchi].patch().faceCells();
+            vsf.boundaryField()[patchi].patch().faceCells();
 
         if (vsf.boundaryField()[patchi].coupled())
         {
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
index 3a6e8233fa..4bf8e9662a 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
@@ -114,10 +114,10 @@ void Foam::leastSquaresVectors::calcLeastSquaresVectors()
     }
 
 
-    surfaceVectorField::GeometricBoundaryField& blsP =
-        pVectors_.boundaryField();
+    surfaceVectorField::GeometricBoundaryField& pVectorsBf =
+        pVectors_.boundaryFieldRef();
 
-    forAll(blsP, patchi)
+    forAll(pVectorsBf, patchi)
     {
         const fvsPatchScalarField& pw = w.boundaryField()[patchi];
         const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
@@ -168,9 +168,9 @@ void Foam::leastSquaresVectors::calcLeastSquaresVectors()
         nVectors_[facei] = -w[facei]*magSfByMagSqrd*(invDd[nei] & d);
     }
 
-    forAll(blsP, patchi)
+    forAll(pVectorsBf, patchi)
     {
-        fvsPatchVectorField& patchLsP = blsP[patchi];
+        fvsPatchVectorField& patchLsP = pVectorsBf[patchi];
 
         const fvsPatchScalarField& pw = w.boundaryField()[patchi];
         const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C
index e67bc081e3..b0eb52a108 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C
@@ -116,7 +116,7 @@ Foam::fv::faceCorrectedSnGrad<Type>::fullGradCorrection
         sfCorr[facei] = dCorr&fgrad;
     }
 
-    tsfCorr.ref().boundaryField() = Zero;
+    tsfCorr.ref().boundaryFieldRef() = Zero;
 
     return tsfCorr;
 }
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C b/src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C
index f5c28f080c..7be3a297fd 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C
@@ -138,18 +138,20 @@ snGradScheme<Type>::snGrad
             deltaCoeffs[facei]*(vf[neighbour[facei]] - vf[owner[facei]]);
     }
 
+    typename GeometricField<Type, fvsPatchField, surfaceMesh>::
+        GeometricBoundaryField& ssfbf = ssf.boundaryFieldRef();
+
     forAll(vf.boundaryField(), patchi)
     {
         const fvPatchField<Type>& pvf = vf.boundaryField()[patchi];
 
         if (pvf.coupled())
         {
-            ssf.boundaryField()[patchi] =
-                pvf.snGrad(tdeltaCoeffs().boundaryField()[patchi]);
+            ssfbf[patchi] = pvf.snGrad(tdeltaCoeffs().boundaryField()[patchi]);
         }
         else
         {
-            ssf.boundaryField()[patchi] = pvf.snGrad();
+            ssfbf[patchi] = pvf.snGrad();
         }
     }
 
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 1d1e9ee529..afb76abede 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -313,7 +313,7 @@ Foam::fvMatrix<Type>::fvMatrix
        const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);
 
     label currentStatePsi = psiRef.eventNo();
-    psiRef.boundaryField().updateCoeffs();
+    psiRef.boundaryFieldRef().updateCoeffs();
     psiRef.eventNo() = currentStatePsi;
 }
 
@@ -928,10 +928,12 @@ flux() const
         }
     }
 
-    forAll(fieldFlux.boundaryField(), patchI)
+    typename GeometricField<Type, fvsPatchField, surfaceMesh>::
+        GeometricBoundaryField& ffbf = fieldFlux.boundaryFieldRef();
+
+    forAll(ffbf, patchI)
     {
-        fieldFlux.boundaryField()[patchI] =
-            InternalContrib[patchI] - NeighbourContrib[patchI];
+        ffbf[patchI] = InternalContrib[patchI] - NeighbourContrib[patchI];
     }
 
     if (faceFluxCorrectionPtr_)
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
index 26dab20b75..991d899ad8 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
@@ -241,7 +241,7 @@ Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::solveCoupled
     addBoundaryDiag(coupledMatrix.diag(), 0);
     addBoundarySource(coupledMatrix.source(), false);
 
-    coupledMatrix.interfaces() = psi.boundaryField().interfaces();
+    coupledMatrix.interfaces() = psi.boundaryFieldRef().interfaces();
     coupledMatrix.interfacesUpper() = boundaryCoeffs().component(0);
     coupledMatrix.interfacesLower() = internalCoeffs().component(0);
 
diff --git a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
index 46fa555289..252ac0bb79 100644
--- a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
@@ -164,7 +164,7 @@ Foam::solverPerformance Foam::fvMatrix<Foam::scalar>::solveSegregated
         *this,
         boundaryCoeffs_,
         internalCoeffs_,
-        psi.boundaryField().scalarInterfaces(),
+        psi_.boundaryField().scalarInterfaces(),
         solverControls
     )->solve(psi.internalField(), totalSource);
 
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
index 00480d48c0..e3f40d1638 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -197,7 +197,7 @@ void Foam::MULES::limiterCorr
 
     scalarField& lambdaIf = lambda;
     surfaceScalarField::GeometricBoundaryField& lambdaBf =
-        lambda.boundaryField();
+        lambda.boundaryFieldRef();
 
     scalarField psiMaxn(psiIf.size(), psiMin);
     scalarField psiMinn(psiIf.size(), psiMax);
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index 61e0d88550..6673e5d983 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -229,7 +229,7 @@ void Foam::MULES::limiter
 
     scalarField& lambdaIf = lambda;
     surfaceScalarField::GeometricBoundaryField& lambdaBf =
-        lambda.boundaryField();
+        lambda.boundaryFieldRef();
 
     scalarField psiMaxn(psiIf.size(), psiMin);
     scalarField psiMinn(psiIf.size(), psiMax);
diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencilTemplates.C b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencilTemplates.C
index d6aac47191..0870d3bf7f 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencilTemplates.C
+++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedCellToFaceStencilTemplates.C
@@ -135,7 +135,7 @@ Foam::extendedCellToFaceStencil::weightedSum
     // Boundaries. Either constrained or calculated so assign value
     // directly (instead of nicely using operator==)
     typename GeometricField<Type, fvsPatchField, surfaceMesh>::
-        GeometricBoundaryField& bSfCorr = sf.boundaryField();
+        GeometricBoundaryField& bSfCorr = sf.boundaryFieldRef();
 
     forAll(bSfCorr, patchi)
     {
diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencilTemplates.C b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencilTemplates.C
index e47683cb48..460b482b07 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencilTemplates.C
+++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencilTemplates.C
@@ -98,7 +98,7 @@ Foam::extendedUpwindCellToFaceStencil::weightedSum
     // Boundaries. Either constrained or calculated so assign value
     // directly (instead of nicely using operator==)
     typename GeometricField<Type, fvsPatchField, surfaceMesh>::
-        GeometricBoundaryField& bSfCorr = sf.boundaryField();
+        GeometricBoundaryField& bSfCorr = sf.boundaryFieldRef();
 
     forAll(bSfCorr, patchi)
     {
diff --git a/src/finiteVolume/fvMesh/fvMesh.C b/src/finiteVolume/fvMesh/fvMesh.C
index 1b26f9a55d..09c81f5845 100644
--- a/src/finiteVolume/fvMesh/fvMesh.C
+++ b/src/finiteVolume/fvMesh/fvMesh.C
@@ -768,10 +768,12 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
 
     const fvPatchList& patches = boundary();
 
+    surfaceScalarField::GeometricBoundaryField& phibf = phi.boundaryFieldRef();
+
     forAll(patches, patchI)
     {
-        phi.boundaryField()[patchI] = patches[patchI].patchSlice(sweptVols);
-        phi.boundaryField()[patchI] *= rDeltaT;
+        phibf[patchI] = patches[patchI].patchSlice(sweptVols);
+        phibf[patchI] *= rDeltaT;
     }
 
     // Update or delete the local geometric properties as early as possible so
diff --git a/src/finiteVolume/fvMesh/fvMeshGeometry.C b/src/finiteVolume/fvMesh/fvMeshGeometry.C
index 8cdeca38a4..1a2c77c5de 100644
--- a/src/finiteVolume/fvMesh/fvMeshGeometry.C
+++ b/src/finiteVolume/fvMesh/fvMeshGeometry.C
@@ -410,9 +410,12 @@ Foam::tmp<Foam::surfaceVectorField> Foam::fvMesh::delta() const
         delta[facei] = C[neighbour[facei]] - C[owner[facei]];
     }
 
-    forAll(delta.boundaryField(), patchi)
+    surfaceVectorField::GeometricBoundaryField& deltabf =
+        delta.boundaryFieldRef();
+
+    forAll(deltabf, patchi)
     {
-        delta.boundaryField()[patchi] = boundary()[patchi].delta();
+        deltabf[patchi] = boundary()[patchi].delta();
     }
 
     return tdelta;
diff --git a/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C b/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
index 7f3923aef8..eea4d29bbf 100644
--- a/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
+++ b/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
@@ -107,7 +107,7 @@ tmp<GeometricField<Type, fvPatchField, volMesh>> fvMeshSubset::interpolate
     //  constructor (with reference to the now correct internal field)
 
     typename GeometricField<Type, fvPatchField, volMesh>::
-        GeometricBoundaryField& bf = resF.boundaryField();
+        GeometricBoundaryField& bf = resF.boundaryFieldRef();
 
     forAll(bf, patchI)
     {
@@ -250,7 +250,7 @@ tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> fvMeshSubset::interpolate
     //  constructor (with reference to the now correct internal field)
 
     typename GeometricField<Type, fvsPatchField, surfaceMesh>::
-        GeometricBoundaryField& bf = resF.boundaryField();
+        GeometricBoundaryField& bf = resF.boundaryFieldRef();
 
     forAll(bf, patchI)
     {
@@ -413,7 +413,7 @@ fvMeshSubset::interpolate
     //  constructor (with reference to the now correct internal field)
 
     typename GeometricField<Type, pointPatchField, pointMesh>::
-        GeometricBoundaryField& bf = resF.boundaryField();
+        GeometricBoundaryField& bf = resF.boundaryFieldRef();
 
     forAll(bf, patchI)
     {
diff --git a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/Poisson/PoissonPatchDistMethod.C b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/Poisson/PoissonPatchDistMethod.C
index 03e03c6dfe..fa3d382d66 100644
--- a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/Poisson/PoissonPatchDistMethod.C
+++ b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/Poisson/PoissonPatchDistMethod.C
@@ -90,7 +90,7 @@ bool Foam::patchDistMethods::Poisson::correct
                 ),
                 mesh_,
                 dimensionedScalar("yPsi", sqr(dimLength), 0.0),
-                y.boundaryField().types()
+                y.boundaryFieldRef().types()
             )
         );
     }
diff --git a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.C b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.C
index 87ca3336d1..e2201a2408 100644
--- a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.C
+++ b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.C
@@ -110,11 +110,12 @@ bool Foam::patchDistMethods::advectionDiffusion::correct
     );
 
     const fvPatchList& patches = mesh_.boundary();
+    volVectorField::GeometricBoundaryField& nybf = ny.boundaryFieldRef();
 
     forAllConstIter(labelHashSet, patchIDs_, iter)
     {
         label patchi = iter.key();
-        ny.boundaryField()[patchi] == -patches[patchi].nf();
+        nybf[patchi] == -patches[patchi].nf();
     }
 
     int iter = 0;
diff --git a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/meshWave/meshWavePatchDistMethod.C b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/meshWave/meshWavePatchDistMethod.C
index 7893341b27..3bc87d3911 100644
--- a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/meshWave/meshWavePatchDistMethod.C
+++ b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/meshWave/meshWavePatchDistMethod.C
@@ -84,13 +84,15 @@ bool Foam::patchDistMethods::meshWave::correct(volScalarField& y)
     y.transfer(wave.distance());
 
     // Transfer values on patches into boundaryField of y
-    forAll(y.boundaryField(), patchI)
+    volScalarField::GeometricBoundaryField& ybf = y.boundaryFieldRef();
+
+    forAll(ybf, patchI)
     {
-        if (!isA<emptyFvPatchScalarField>(y.boundaryField()[patchI]))
+        if (!isA<emptyFvPatchScalarField>(ybf[patchI]))
         {
             scalarField& waveFld = wave.patchDistance()[patchI];
 
-            y.boundaryField()[patchI].transfer(waveFld);
+            ybf[patchI].transfer(waveFld);
         }
     }
 
@@ -112,9 +114,11 @@ bool Foam::patchDistMethods::meshWave::correct
     // Collect pointers to data on patches
     UPtrList<vectorField> patchData(mesh_.boundaryMesh().size());
 
-    forAll(n.boundaryField(), patchI)
+    volVectorField::GeometricBoundaryField& nbf = n.boundaryFieldRef();
+
+    forAll(nbf, patchI)
     {
-        patchData.set(patchI, &n.boundaryField()[patchI]);
+        patchData.set(patchI, &nbf[patchI]);
     }
 
     // Do mesh wave
@@ -132,17 +136,19 @@ bool Foam::patchDistMethods::meshWave::correct
     n.transfer(wave.cellData());
 
     // Transfer values on patches into boundaryField of y and n
-    forAll(y.boundaryField(), patchI)
+    volScalarField::GeometricBoundaryField& ybf = y.boundaryFieldRef();
+
+    forAll(ybf, patchI)
     {
         scalarField& waveFld = wave.patchDistance()[patchI];
 
-        if (!isA<emptyFvPatchScalarField>(y.boundaryField()[patchI]))
+        if (!isA<emptyFvPatchScalarField>(ybf[patchI]))
         {
-            y.boundaryField()[patchI].transfer(waveFld);
+            ybf[patchI].transfer(waveFld);
 
             vectorField& wavePatchData = wave.patchData()[patchI];
 
-            n.boundaryField()[patchI].transfer(wavePatchData);
+            nbf[patchI].transfer(wavePatchData);
         }
     }
 
diff --git a/src/finiteVolume/fvMesh/wallDist/wallDist/wallDist.C b/src/finiteVolume/fvMesh/wallDist/wallDist/wallDist.C
index ae5be722b1..9a3ba0974d 100644
--- a/src/finiteVolume/fvMesh/wallDist/wallDist/wallDist.C
+++ b/src/finiteVolume/fvMesh/wallDist/wallDist/wallDist.C
@@ -56,10 +56,12 @@ void Foam::wallDist::constructn() const
 
     const fvPatchList& patches = mesh().boundary();
 
+    volVectorField::GeometricBoundaryField& nbf = n_.ref().boundaryFieldRef();
+
     forAllConstIter(labelHashSet, patchIDs_, iter)
     {
         label patchi = iter.key();
-        n_.ref().boundaryField()[patchi] == patches[patchi].nf();
+        nbf[patchi] == patches[patchi].nf();
     }
 }
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
index 8b523015e1..1a9b23a807 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
@@ -77,7 +77,7 @@ void Foam::LimitedScheme<Type, Limiter, LimitFunc>::calcLimiter
     }
 
     surfaceScalarField::GeometricBoundaryField& bLim =
-        limiterField.boundaryField();
+        limiterField.boundaryFieldRef();
 
     forAll(bLim, patchi)
     {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/PhiScheme/PhiScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/PhiScheme/PhiScheme.C
index 3101d60893..f3c7bce6d1 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/PhiScheme/PhiScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/PhiScheme/PhiScheme.C
@@ -100,7 +100,7 @@ Foam::PhiScheme<Type, PhiLimiter>::limiter
 
 
     surfaceScalarField::GeometricBoundaryField& bLimiter =
-        Limiter.boundaryField();
+        Limiter.boundaryFieldRef();
 
     forAll(bLimiter, patchI)
     {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationScheme.C
index 138fad0ee7..b3dc8a993b 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationScheme.C
@@ -159,7 +159,7 @@ Foam::limitedSurfaceInterpolationScheme<Type>::weights
     }
 
     surfaceScalarField::GeometricBoundaryField& bWeights =
-        Weights.boundaryField();
+        Weights.boundaryFieldRef();
 
     forAll(bWeights, patchI)
     {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/clippedLinear/clippedLinear.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/clippedLinear/clippedLinear.H
index 59ffd4bac3..d07146ef34 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/clippedLinear/clippedLinear.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/clippedLinear/clippedLinear.H
@@ -156,11 +156,14 @@ public:
             clippedLinearWeights.internalField() =
                 max(min(cdWeights.internalField(), 1 - wfLimit_), wfLimit_);
 
+            surfaceScalarField::GeometricBoundaryField& clwbf =
+                clippedLinearWeights.boundaryFieldRef();
+
             forAll(mesh.boundary(), patchi)
             {
-                if (clippedLinearWeights.boundaryField()[patchi].coupled())
+                if (clwbf[patchi].coupled())
                 {
-                    clippedLinearWeights.boundaryField()[patchi] =
+                    clwbf[patchi] =
                         max
                         (
                             min
@@ -173,8 +176,7 @@ public:
                 }
                 else
                 {
-                    clippedLinearWeights.boundaryField()[patchi] =
-                        cdWeights.boundaryField()[patchi];
+                    clwbf[patchi] = cdWeights.boundaryField()[patchi];
                 }
             }
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubic/cubic.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubic/cubic.H
index 18b767b002..cdca6ca352 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubic/cubic.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubic/cubic.H
@@ -171,11 +171,14 @@ public:
                 );
             }
 
-            forAll(sfCorr.boundaryField(), pi)
+            typename GeometricField<Type, fvsPatchField, surfaceMesh>::
+                GeometricBoundaryField& sfCorrbf = sfCorr.boundaryFieldRef();
+
+            forAll(sfCorrbf, pi)
             {
-                if (!sfCorr.boundaryField()[pi].coupled())
+                if (!sfCorrbf[pi].coupled())
                 {
-                    sfCorr.boundaryField()[pi] = Zero;
+                    sfCorrbf[pi] = Zero;
                 }
             }
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.C
index 6dad60e254..be58785083 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.C
@@ -90,7 +90,7 @@ Foam::linearUpwind<Type>::correction
 
 
     typename GeometricField<Type, fvsPatchField, surfaceMesh>::
-        GeometricBoundaryField& bSfCorr = sfCorr.boundaryField();
+        GeometricBoundaryField& bSfCorr = sfCorr.boundaryFieldRef();
 
     forAll(bSfCorr, patchi)
     {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.C
index f8237e1d35..910569e97a 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.C
@@ -129,7 +129,7 @@ Foam::linearUpwindV<Type>::correction
 
 
     typename GeometricField<Type, fvsPatchField, surfaceMesh>::
-        GeometricBoundaryField& bSfCorr = sfCorr.boundaryField();
+        GeometricBoundaryField& bSfCorr = sfCorr.boundaryFieldRef();
 
     forAll(bSfCorr, patchi)
     {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/localMax/localMax.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/localMax/localMax.H
index e081db91e0..bfa12beb77 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/localMax/localMax.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/localMax/localMax.H
@@ -138,9 +138,12 @@ public:
             );
             GeometricField<Type, fvsPatchField, surfaceMesh>& vff = tvff.ref();
 
-            forAll(vff.boundaryField(), patchi)
+            typename GeometricField<Type, fvsPatchField, surfaceMesh>::
+                GeometricBoundaryField& vffbf = vff.boundaryFieldRef();
+
+            forAll(vffbf, patchi)
             {
-                vff.boundaryField()[patchi] = vf.boundaryField()[patchi];
+                vffbf[patchi] = vf.boundaryField()[patchi];
             }
 
             const labelUList& own = mesh.owner();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/localMin/localMin.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/localMin/localMin.H
index 7099656ddd..56cb815ff8 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/localMin/localMin.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/localMin/localMin.H
@@ -138,9 +138,12 @@ public:
             );
             GeometricField<Type, fvsPatchField, surfaceMesh>& vff = tvff.ref();
 
-            forAll(vff.boundaryField(), patchi)
+            typename GeometricField<Type, fvsPatchField, surfaceMesh>::
+                GeometricBoundaryField& vffbf = vff.boundaryFieldRef();
+
+            forAll(vffbf, patchi)
             {
-                vff.boundaryField()[patchi] = vf.boundaryField()[patchi];
+                vffbf[patchi] = vf.boundaryField()[patchi];
             }
 
             const labelUList& own = mesh.owner();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/midPoint/midPoint.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/midPoint/midPoint.H
index 021abd7685..3b2b8a2710 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/midPoint/midPoint.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/midPoint/midPoint.H
@@ -116,14 +116,14 @@ public:
                 )
             );
 
-            surfaceScalarField::GeometricBoundaryField& awbf =
-                taw.ref().boundaryField();
+            surfaceScalarField::GeometricBoundaryField& awBf =
+                taw.ref().boundaryFieldRef();
 
-            forAll(awbf, patchi)
+            forAll(awBf, patchi)
             {
-                if (!awbf[patchi].coupled())
+                if (!awBf[patchi].coupled())
                 {
-                    awbf[patchi] = 1.0;
+                    awBf[patchi] = 1.0;
                 }
             }
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C
index 303421fcf6..6f665ee419 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C
@@ -104,7 +104,7 @@ correction
 
 
     typename GeometricField<Type, fvsPatchField, surfaceMesh>::
-        GeometricBoundaryField& bSfCorr = tsfCorr.ref().boundaryField();
+        GeometricBoundaryField& bSfCorr = tsfCorr.ref().boundaryFieldRef();
 
     forAll(bSfCorr, patchi)
     {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/reverseLinear/reverseLinear.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/reverseLinear/reverseLinear.H
index e84e66f251..ffc68526b5 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/reverseLinear/reverseLinear.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/reverseLinear/reverseLinear.H
@@ -128,17 +128,19 @@ public:
             reverseLinearWeights.internalField() =
                 1.0 - cdWeights.internalField();
 
+            surfaceScalarField::GeometricBoundaryField& rlwbf =
+                reverseLinearWeights.boundaryFieldRef();
+
+
             forAll(mesh.boundary(), patchI)
             {
-                if (reverseLinearWeights.boundaryField()[patchI].coupled())
+                if (rlwbf[patchI].coupled())
                 {
-                    reverseLinearWeights.boundaryField()[patchI] =
-                        1.0 - cdWeights.boundaryField()[patchI];
+                    rlwbf[patchI] = 1.0 - cdWeights.boundaryField()[patchI];
                 }
                 else
                 {
-                    reverseLinearWeights.boundaryField()[patchI] =
-                        cdWeights.boundaryField()[patchI];
+                    rlwbf[patchI] = cdWeights.boundaryField()[patchI];
                 }
             }
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.C
index 4e7bf60dd4..ad51c93b7c 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.C
@@ -90,11 +90,12 @@ void Foam::skewCorrectionVectors::calcSkewCorrectionVectors()
             Cpf - ((Sf[facei] & Cpf)/(Sf[facei] & d))*d;
     }
 
+    typename surfaceVectorField::GeometricBoundaryField& skewCorrVecsBf =
+        skewCorrectionVectors_.boundaryFieldRef();
 
-    forAll(skewCorrectionVectors_.boundaryField(), patchI)
+    forAll(skewCorrVecsBf, patchi)
     {
-        fvsPatchVectorField& patchSkewCorrVecs =
-            skewCorrectionVectors_.boundaryField()[patchI];
+        fvsPatchVectorField& patchSkewCorrVecs = skewCorrVecsBf[patchi];
 
         if (!patchSkewCorrVecs.coupled())
         {
@@ -104,8 +105,8 @@ void Foam::skewCorrectionVectors::calcSkewCorrectionVectors()
         {
             const fvPatch& p = patchSkewCorrVecs.patch();
             const labelUList& faceCells = p.faceCells();
-            const vectorField& patchFaceCentres = Cf.boundaryField()[patchI];
-            const vectorField& patchSf = Sf.boundaryField()[patchI];
+            const vectorField& patchFaceCentres = Cf.boundaryField()[patchi];
+            const vectorField& patchSf = Sf.boundaryField()[patchi];
             const vectorField patchD(p.delta());
 
             forAll(p, patchFaceI)
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
index 6195180a87..0b7e3edc4b 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
@@ -36,7 +36,7 @@ Description
 
 namespace Foam
 {
-defineTypeNameAndDebug(surfaceInterpolation, 0);
+    defineTypeNameAndDebug(surfaceInterpolation, 0);
 }
 
 
@@ -183,12 +183,12 @@ void Foam::surfaceInterpolation::makeWeights() const
         w[facei] = SfdNei/(SfdOwn + SfdNei);
     }
 
+    typename surfaceScalarField::GeometricBoundaryField& wBf =
+        weights.boundaryFieldRef();
+
     forAll(mesh_.boundary(), patchi)
     {
-        mesh_.boundary()[patchi].makeWeights
-        (
-            weights.boundaryField()[patchi]
-        );
+        mesh_.boundary()[patchi].makeWeights(wBf[patchi]);
     }
 
     if (debug)
@@ -227,7 +227,7 @@ void Foam::surfaceInterpolation::makeDeltaCoeffs() const
         mesh_,
         dimless/dimLength
     );
-    surfaceScalarField& DeltaCoeffs = *deltaCoeffs_;
+    surfaceScalarField& deltaCoeffs = *deltaCoeffs_;
 
 
     // Set local references to mesh data
@@ -237,13 +237,15 @@ void Foam::surfaceInterpolation::makeDeltaCoeffs() const
 
     forAll(owner, facei)
     {
-        DeltaCoeffs[facei] = 1.0/mag(C[neighbour[facei]] - C[owner[facei]]);
+        deltaCoeffs[facei] = 1.0/mag(C[neighbour[facei]] - C[owner[facei]]);
     }
 
-    forAll(DeltaCoeffs.boundaryField(), patchi)
+    typename surfaceScalarField::GeometricBoundaryField& deltaCoeffsBf =
+        deltaCoeffs.boundaryFieldRef();
+
+    forAll(deltaCoeffsBf, patchi)
     {
-        DeltaCoeffs.boundaryField()[patchi] =
-            1.0/mag(mesh_.boundary()[patchi].delta());
+        deltaCoeffsBf[patchi] = 1.0/mag(mesh_.boundary()[patchi].delta());
     }
 }
 
@@ -303,11 +305,14 @@ void Foam::surfaceInterpolation::makeNonOrthDeltaCoeffs() const
         nonOrthDeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
     }
 
-    forAll(nonOrthDeltaCoeffs.boundaryField(), patchi)
+    typename surfaceScalarField::GeometricBoundaryField& nonOrthDeltaCoeffsBf =
+        nonOrthDeltaCoeffs.boundaryFieldRef();
+
+    forAll(nonOrthDeltaCoeffsBf, patchi)
     {
         vectorField delta(mesh_.boundary()[patchi].delta());
 
-        nonOrthDeltaCoeffs.boundaryField()[patchi] =
+        nonOrthDeltaCoeffsBf[patchi] =
             1.0/max(mesh_.boundary()[patchi].nf() & delta, 0.05*mag(delta));
     }
 }
@@ -358,9 +363,12 @@ void Foam::surfaceInterpolation::makeNonOrthCorrectionVectors() const
     // and calculated consistently with internal corrections for
     // coupled patches
 
-    forAll(corrVecs.boundaryField(), patchi)
+    typename surfaceVectorField::GeometricBoundaryField& corrVecsBf =
+        corrVecs.boundaryFieldRef();
+
+    forAll(corrVecsBf, patchi)
     {
-        fvsPatchVectorField& patchCorrVecs = corrVecs.boundaryField()[patchi];
+        fvsPatchVectorField& patchCorrVecs = corrVecsBf[patchi];
 
         if (!patchCorrVecs.coupled())
         {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
index 4f6a79c9a2..8e16fb1f73 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
@@ -190,6 +190,8 @@ Foam::surfaceInterpolationScheme<Type>::interpolate
 
 
     // Interpolate across coupled patches using given lambdas and ys
+    typename GeometricField<Type, fvsPatchField, surfaceMesh>::
+        GeometricBoundaryField& sfbf = sf.boundaryFieldRef();
 
     forAll(lambdas.boundaryField(), pi)
     {
@@ -198,13 +200,13 @@ Foam::surfaceInterpolationScheme<Type>::interpolate
 
         if (vf.boundaryField()[pi].coupled())
         {
-            sf.boundaryField()[pi] =
+            sfbf[pi] =
                 pLambda*vf.boundaryField()[pi].patchInternalField()
               + pY*vf.boundaryField()[pi].patchNeighbourField();
         }
         else
         {
-            sf.boundaryField()[pi] = vf.boundaryField()[pi];
+            sfbf[pi] = vf.boundaryField()[pi];
         }
     }
 
@@ -283,11 +285,14 @@ Foam::surfaceInterpolationScheme<Type>::dotInterpolate
 
     // Interpolate across coupled patches using given lambdas
 
+    typename GeometricField<RetType, fvsPatchField, surfaceMesh>::
+        GeometricBoundaryField& sfbf = sf.boundaryFieldRef();
+
     forAll(lambdas.boundaryField(), pi)
     {
         const fvsPatchScalarField& pLambda = lambdas.boundaryField()[pi];
         const typename SFType::PatchFieldType& pSf = Sf.boundaryField()[pi];
-        fvsPatchField<RetType>& psf = sf.boundaryField()[pi];
+        fvsPatchField<RetType>& psf = sfbf[pi];
 
         if (vf.boundaryField()[pi].coupled())
         {
diff --git a/src/finiteVolume/interpolation/volPointInterpolation/pointConstraintsTemplates.C b/src/finiteVolume/interpolation/volPointInterpolation/pointConstraintsTemplates.C
index 85491c2688..4c4bf2d4e4 100644
--- a/src/finiteVolume/interpolation/volPointInterpolation/pointConstraintsTemplates.C
+++ b/src/finiteVolume/interpolation/volPointInterpolation/pointConstraintsTemplates.C
@@ -91,9 +91,12 @@ void Foam::pointConstraints::setPatchFields
     GeometricField<Type, pointPatchField, pointMesh>& pf
 )
 {
-    forAll(pf.boundaryField(), patchI)
+    typename GeometricField<Type, pointPatchField, pointMesh>::
+        GeometricBoundaryField& pfbf = pf.boundaryFieldRef();
+
+    forAll(pfbf, patchI)
     {
-        pointPatchField<Type>& ppf = pf.boundaryField()[patchI];
+        pointPatchField<Type>& ppf = pfbf[patchI];
 
         if (isA<valuePointPatchField<Type>>(ppf))
         {
diff --git a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
index 0f4750b279..74b2f6b655 100644
--- a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
+++ b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
@@ -86,12 +86,15 @@ void Foam::volPointInterpolation::addSeparated
         Pout<< "volPointInterpolation::addSeparated" << endl;
     }
 
-    forAll(pf.boundaryField(), patchI)
+    typename GeometricField<Type, pointPatchField, pointMesh>::
+        GeometricBoundaryField& pfbf = pf.boundaryFieldRef();
+
+    forAll(pfbf, patchI)
     {
-        if (pf.boundaryField()[patchI].coupled())
+        if (pfbf[patchI].coupled())
         {
             refCast<coupledPointPatchField<Type>>
-                (pf.boundaryField()[patchI]).initSwapAddSeparated
+                (pfbf[patchI]).initSwapAddSeparated
                 (
                     Pstream::nonBlocking,
                     pf.internalField()
@@ -102,12 +105,12 @@ void Foam::volPointInterpolation::addSeparated
     // Block for any outstanding requests
     Pstream::waitRequests();
 
-    forAll(pf.boundaryField(), patchI)
+    forAll(pfbf, patchI)
     {
-        if (pf.boundaryField()[patchI].coupled())
+        if (pfbf[patchI].coupled())
         {
             refCast<coupledPointPatchField<Type>>
-                (pf.boundaryField()[patchI]).swapAddSeparated
+                (pfbf[patchI]).swapAddSeparated
                 (
                     Pstream::nonBlocking,
                     pf.internalField()