diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index f2f2af745e250c84a87b0d6c10f4e28c15e7b3df..720e01f542819e23d79b0d8a22264af4210690cd 100644
--- a/applications/solvers/DNS/dnsFoam/dnsFoam.C
+++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C
@@ -61,7 +61,7 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        force.internalField() = ReImSum
+        force.internalFieldRef() = ReImSum
         (
             fft::reverseTransform
             (
diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.C b/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.C
index e85e41b432c6a109979a1a8689f5007d6587fffe..82bbbede1f3808c39424d1c35ba5bdf4b555c7ea 100644
--- a/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.C
@@ -105,7 +105,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::basicSubGrid::XiEq() const
         mesh,
         dimensionedScalar("zero", Nv.dimensions(), 0.0)
     );
-    N.internalField() = Nv.internalField()*Cw;
+    N.internalFieldRef() = Nv.internalField()*Cw;
 
     volSymmTensorField ns
     (
@@ -125,7 +125,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::basicSubGrid::XiEq() const
             Zero
         )
     );
-    ns.internalField() = nsv.internalField()*Cw;
+    ns.internalFieldRef() = nsv.internalField()*Cw;
 
     volScalarField n(max(N - (Uhat & ns & Uhat), scalar(1e-4)));
     volScalarField b((Uhat & B_ & Uhat)/sqrt(n));
diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C
index 46dee65938d3a66e9583e7962c2665dd47c145e8..20ae2215cdcd9bf21fa0e47c7c9e638e165654b9 100644
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.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
@@ -74,7 +74,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::Gulder::XiEq() const
 
     if (subGridSchelkin_)
     {
-        up.internalField() += calculateSchelkinEffect(uPrimeCoef_);
+        up.internalFieldRef() += calculateSchelkinEffect(uPrimeCoef_);
     }
 
     volScalarField tauEta(sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon))));
diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
index 915567939664d6680efb073b5fc66f5487d5a930..6a585e388182e90b12ec3239c149d4aec4e32a7d 100644
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
@@ -82,7 +82,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
     volScalarField up(sqrt((2.0/3.0)*k));
     if (subGridSchelkin_)
     {
-        up.internalField() += calculateSchelkinEffect(uPrimeCoef_);
+        up.internalFieldRef() += calculateSchelkinEffect(uPrimeCoef_);
     }
 
     volScalarField l(lCoef_*sqrt(3.0/2.0)*up*k/epsilon);
diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C
index 966dae717538dc088b1e947e672e0458f8c70ea3..b436b1b0a8f973be580c395bab748e116865d217 100644
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C
@@ -114,7 +114,7 @@ Foam::XiEqModel::calculateSchelkinEffect(const scalar uPrimeCoef) const
         )
     );
     volScalarField& N = tN.ref();
-    N.internalField() = Nv.internalField()*pow(mesh.V(), 2.0/3.0);
+    N.internalFieldRef() = Nv.internalField()*pow(mesh.V(), 2.0/3.0);
 
     volSymmTensorField ns
     (
@@ -134,7 +134,7 @@ Foam::XiEqModel::calculateSchelkinEffect(const scalar uPrimeCoef) const
             Zero
         )
     );
-    ns.internalField() = nsv.internalField()*pow(mesh.V(), 2.0/3.0);
+    ns.internalFieldRef() = nsv.internalField()*pow(mesh.V(), 2.0/3.0);
 
     const volVectorField Uhat
     (
@@ -150,7 +150,7 @@ Foam::XiEqModel::calculateSchelkinEffect(const scalar uPrimeCoef) const
     const scalarField deltaUp(upLocal*(max(scalar(1.0), pow(nr, 0.5)) - 1.0));
 
     // Re use tN
-    N.internalField() = upLocal*(max(scalar(1.0), pow(nr, 0.5)) - 1.0);
+    N.internalFieldRef() = upLocal*(max(scalar(1.0), pow(nr, 0.5)) - 1.0);
 
     return tN;
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
index f900fba316f967f0344d222ab506b21d0e584cd3..13d95cf01adde926a0087ad7a8a360ad849c20e9 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
@@ -61,7 +61,7 @@
                 )
             );
 
-            aniAlphas[i].internalField() =
+            aniAlphas[i].internalFieldRef() =
                 coordinates[i].R().transformVector(tkappaByCp());
             aniAlphas[i].correctBoundaryConditions();
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
index 9eb71909c16d3d260f1e161b0ae8529042a2587b..d96c34e364b8fa3e17f4b4f6a3b53753ffcaa31d 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
@@ -14,7 +14,7 @@ if (!thermo.isotropic())
     tmp<volVectorField> tkappaByCp = thermo.Kappa()/cp;
     const coordinateSystem& coodSys = coordinates[i];
 
-    aniAlpha.internalField() =
+    aniAlpha.internalFieldRef() =
         coodSys.R().transformVector(tkappaByCp());
     aniAlpha.correctBoundaryConditions();
 
diff --git a/applications/solvers/lagrangian/DPMFoam/DPMFoam.C b/applications/solvers/lagrangian/DPMFoam/DPMFoam.C
index 728e3e45d2a346972f329f510b3d81129f3a09c7..9aa262d9b86971b8a3b66bc4a36b296d8dde4fb9 100644
--- a/applications/solvers/lagrangian/DPMFoam/DPMFoam.C
+++ b/applications/solvers/lagrangian/DPMFoam/DPMFoam.C
@@ -107,7 +107,7 @@ int main(int argc, char *argv[])
             zeroGradientFvPatchVectorField::typeName
         );
 
-        cloudVolSUSu.internalField() = -cloudSU.source()/mesh.V();
+        cloudVolSUSu.internalFieldRef() = -cloudSU.source()/mesh.V();
         cloudVolSUSu.correctBoundaryConditions();
         cloudSU.source() = Zero;
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/segregated/segregated.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/segregated/segregated.C
index ed45643189b03bd570163ff4cd2c0e7c5ad26024..a1e50a16f46c76e047ea7912f8a0f70ff7202a68 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/segregated/segregated.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/segregated/segregated.C
@@ -104,7 +104,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
         dimensionedScalar("L", dimLength, 0),
         zeroGradientFvPatchField<scalar>::typeName
     );
-    L.internalField() = cbrt(mesh.V());
+    L.internalFieldRef() = cbrt(mesh.V());
     L.correctBoundaryConditions();
 
     volScalarField I
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.C
index 079b5d187b61968c246b672d1c8ceb263b8358d9..d7b156da80b7773d9c81751da9aa4dd88f17b205 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.C
@@ -104,7 +104,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
         dimensionedScalar("L", dimLength, 0),
         zeroGradientFvPatchField<scalar>::typeName
     );
-    L.internalField() = cbrt(mesh.V());
+    L.internalFieldRef() = cbrt(mesh.V());
     L.correctBoundaryConditions();
 
     volScalarField I
diff --git a/applications/test/GAMGAgglomeration/Test-GAMGAgglomeration.C b/applications/test/GAMGAgglomeration/Test-GAMGAgglomeration.C
index c388a0b221ea8c06a312baee1e1eec62b6409e47..a3b6e68213c1c54a63857b63e3897155289a77bc 100644
--- a/applications/test/GAMGAgglomeration/Test-GAMGAgglomeration.C
+++ b/applications/test/GAMGAgglomeration/Test-GAMGAgglomeration.C
@@ -87,7 +87,7 @@ int main(int argc, char *argv[])
             mesh,
             dimensionedScalar("aggomeration", dimless, 0.0)
         );
-        scalarField& fld = scalarAgglomeration.internalField();
+        scalarField& fld = scalarAgglomeration.internalFieldRef();
         forAll(fld, celli)
         {
             fld[celli] = cellToCoarse[celli];
@@ -166,7 +166,7 @@ int main(int argc, char *argv[])
                 mesh,
                 dimensionedScalar("aggomeration", dimless, 0.0)
             );
-            scalarField& fld = scalarAgglomeration.internalField();
+            scalarField& fld = scalarAgglomeration.internalFieldRef();
             forAll(fld, celli)
             {
                 fld[celli] = cellToCoarse[celli];
diff --git a/applications/test/hexRef8/Test-hexRef8.C b/applications/test/hexRef8/Test-hexRef8.C
index f79628fd8c01983176a7d9f8be1946ab0c9e5b54..737401a627f477ecad0b24d2e996601b8520b097 100644
--- a/applications/test/hexRef8/Test-hexRef8.C
+++ b/applications/test/hexRef8/Test-hexRef8.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
@@ -164,7 +164,7 @@ int main(int argc, char *argv[])
         dimensionedScalar("one", dimless, 1.0),
         calculatedPointPatchScalarField::typeName
     );
-    pointX.internalField() = mesh.points().component(0);
+    pointX.internalFieldRef() = mesh.points().component(0);
     pointX.correctBoundaryConditions();
     Info<< "Writing x-component field "
         << pointX.name() << " in " << runTime.timeName() << endl;
diff --git a/applications/test/mappedPatch/Test-MappedPatch.C b/applications/test/mappedPatch/Test-MappedPatch.C
index 603c5e20620c145201c1d1726928db49e6c3f660..46577e4c8b76b5ce0be18327545aff1f1f3fc00a 100644
--- a/applications/test/mappedPatch/Test-MappedPatch.C
+++ b/applications/test/mappedPatch/Test-MappedPatch.C
@@ -86,7 +86,7 @@ int main(int argc, char *argv[])
         patchFieldTypes
     );
 
-    cc.internalField() = mesh.C().internalField();
+    cc.internalFieldRef() = mesh.C().internalField();
     cc.boundaryFieldRef().updateCoeffs();
 
     forAll(cc.boundaryField(), patchi)
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 573b48984e719f1a2f0e79116db67fb7a3acb51d..1873af19e268559e5f6c8b8c268fc596f568c71a 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -1137,7 +1137,7 @@ void Foam::conformalVoronoiMesh::writeCellSizes
             zeroGradientFvPatchScalarField::typeName
         );
 
-        scalarField& cellSize = targetCellSize.internalField();
+        scalarField& cellSize = targetCellSize.internalFieldRef();
 
         const vectorField& C = mesh.cellCentres();
 
@@ -1163,7 +1163,7 @@ void Foam::conformalVoronoiMesh::writeCellSizes
         //     zeroGradientFvPatchScalarField::typeName
         // );
 
-        // targetCellVolume.internalField() = pow3(cellSize);
+        // targetCellVolume.internalFieldRef() = pow3(cellSize);
 
         // Info<< nl << "Create actualCellVolume volScalarField" << endl;
 
@@ -1182,7 +1182,7 @@ void Foam::conformalVoronoiMesh::writeCellSizes
         //     zeroGradientFvPatchScalarField::typeName
         // );
 
-        // actualCellVolume.internalField() = mesh.cellVolumes();
+        // actualCellVolume.internalFieldRef() = mesh.cellVolumes();
 
         // Info<< nl << "Create equivalentCellSize volScalarField" << endl;
 
@@ -1201,7 +1201,7 @@ void Foam::conformalVoronoiMesh::writeCellSizes
         //     zeroGradientFvPatchScalarField::typeName
         // );
 
-        // equivalentCellSize.internalField() = pow
+        // equivalentCellSize.internalFieldRef() = pow
         // (
         //     actualCellVolume.internalField(),
         //     1.0/3.0
@@ -1247,7 +1247,7 @@ void Foam::conformalVoronoiMesh::writeCellSizes
     //         pointPatchVectorField::calculatedType()
     //     );
 
-    //     scalarField& cellSize = ptTargetCellSize.internalField();
+    //     scalarField& cellSize = ptTargetCellSize.internalFieldRef();
 
     //     const vectorField& P = tetMesh.points();
 
@@ -1283,7 +1283,7 @@ void Foam::conformalVoronoiMesh::writeCellAlignments
 //        zeroGradientFvPatchTensorField::typeName
 //    );
 //
-//    tensorField& cellAlignment = cellAlignments.internalField();
+//    tensorField& cellAlignment = cellAlignments.internalFieldRef();
 //
 //    const vectorField& C = mesh.cellCentres();
 //
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshTemplates.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshTemplates.C
index 55fd3521d261407072f950ed2378b47ade0df9f7..4167cb61db900cefacfa7883439fe6100e6b1905 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshTemplates.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshTemplates.C
@@ -117,16 +117,14 @@ bool Foam::conformalVoronoiMesh::distributeBackground(const Triangulation& mesh)
             }
         }
 
+        scalarField& cwi = cellWeights.internalFieldRef();
+
         forAll(cellVertices, cI)
         {
             // Give a small but finite weight for empty cells.  Some
             // decomposition methods have difficulty with integer overflows in
             // the sum of the normalised weight field.
-            cellWeights.internalField()[cI] = max
-            (
-                cellVertices[cI],
-                1e-2
-            );
+            cwi[cI] = max(cellVertices[cI], 1e-2);
         }
 
         autoPtr<mapDistributePolyMesh> mapDist = decomposition_().distribute
diff --git a/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C b/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C
index 18b2bf7da4378b9878749e915dff330b66a47353..9887130f4b60b7c54ba1ed8576deeadd74215f6a 100644
--- a/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C
@@ -636,7 +636,7 @@ int main(int argc, char *argv[])
 
         {
             // Internal field
-            cellDistance.internalField() = signedDistance
+            cellDistance.internalFieldRef() = signedDistance
             (
                 distSqr,
                 fvm.C(),
@@ -701,7 +701,7 @@ int main(int argc, char *argv[])
                 -sqr(GREAT)             // null value
             );
 
-            pointDistance.internalField() = signedDistance
+            pointDistance.internalFieldRef() = signedDistance
             (
                 pointDistSqr,
                 fvm.points(),
diff --git a/applications/utilities/postProcessing/dataConversion/foamToTetDualMesh/foamToTetDualMesh.C b/applications/utilities/postProcessing/dataConversion/foamToTetDualMesh/foamToTetDualMesh.C
index 7ee8b33279beedb83eb2c3bebb967d38fb018408..6017e27fe20abc5eb2cfe11ab97d40b863b95bc3 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToTetDualMesh/foamToTetDualMesh.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToTetDualMesh/foamToTetDualMesh.C
@@ -91,7 +91,7 @@ void ReadAndMapFields
             )
         );
 
-        Field<Type>& fld = tetFields[i].internalField();
+        Field<Type>& fld = tetFields[i].internalFieldRef();
 
         // Map from read field. Set unmapped entries to nullValue.
         fld.setSize(map.size(), nullValue);
diff --git a/applications/utilities/preProcessing/boxTurb/boxTurb.C b/applications/utilities/preProcessing/boxTurb/boxTurb.C
index 55a35322abad886d8ba68dea721168ca2027fe8e..3eb161016ef9146282c61fbb40b5778414bdc083 100644
--- a/applications/utilities/preProcessing/boxTurb/boxTurb.C
+++ b/applications/utilities/preProcessing/boxTurb/boxTurb.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,7 +57,7 @@ int main(int argc, char *argv[])
 
     turbGen Ugen(K, Ea, k0);
 
-    U.internalField() = Ugen.U();
+    U.internalFieldRef() = Ugen.U();
     U.correctBoundaryConditions();
 
     Info<< "k("
diff --git a/applications/utilities/preProcessing/setFields/setFields.C b/applications/utilities/preProcessing/setFields/setFields.C
index 422b47249630d00960cb81e363b9a5315b8f57ca..59eebecd48e7a96b6b9344a4e79b15d3fd7e6bf9 100644
--- a/applications/utilities/preProcessing/setFields/setFields.C
+++ b/applications/utilities/preProcessing/setFields/setFields.C
@@ -88,7 +88,7 @@ bool setCellFieldType
 
         if (selectedCells.size() == field.size())
         {
-            field.internalField() = value;
+            field.internalFieldRef() = value;
         }
         else
         {
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
index be12611465c008ffe427591c223ced1d319b6c4f..2d43b495cf060de0e55c6190ee325a57ff904e98 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
@@ -50,7 +50,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::readFields
     const dictionary& dict
 )
 {
-    DimensionedField<Type, GeoMesh>::readField(dict, "internalField");
+    Internal::readField(dict, "internalField");
 
     boundaryField_.readField(*this, dict.subDict("boundaryField"));
 
@@ -181,7 +181,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const word& patchFieldType
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, ds, false),
+    Internal(io, mesh, ds, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -206,7 +206,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const wordList& actualPatchTypes
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, ds, false),
+    Internal(io, mesh, ds, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -230,7 +230,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const word& patchFieldType
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, dt, false),
+    Internal(io, mesh, dt, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -257,7 +257,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const wordList& actualPatchTypes
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, dt, false),
+    Internal(io, mesh, dt, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -274,6 +274,30 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
+(
+    const IOobject& io,
+    const Internal& diField,
+    const PtrList<PatchField<Type>>& ptfl
+)
+:
+    Internal(io, diField),
+    timeIndex_(this->time().timeIndex()),
+    field0Ptr_(NULL),
+    fieldPrevIterPtr_(NULL),
+    boundaryField_(this->mesh().boundary(), *this, ptfl)
+{
+    if (debug)
+    {
+        InfoInFunction
+            << "Constructing from components" << endl << this->info() << endl;
+    }
+
+    readIfPresent();
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 (
@@ -284,7 +308,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const PtrList<PatchField<Type>>& ptfl
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, ds, iField),
+    Internal(io, mesh, ds, iField),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -307,7 +331,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const Mesh& mesh
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, dimless, false),
+    Internal(io, mesh, dimless, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -343,7 +367,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const dictionary& dict
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, dimless, false),
+    Internal(io, mesh, dimless, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -378,7 +402,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const GeometricField<Type, PatchField, GeoMesh>& gf
 )
 :
-    DimensionedField<Type, GeoMesh>(gf),
+    Internal(gf),
     timeIndex_(gf.timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -409,7 +433,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
 )
 :
-    DimensionedField<Type, GeoMesh>
+    Internal
     (
         const_cast<GeometricField<Type, PatchField, GeoMesh>&>(tgf()),
         tgf.isTmp()
@@ -439,7 +463,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const GeometricField<Type, PatchField, GeoMesh>& gf
 )
 :
-    DimensionedField<Type, GeoMesh>(io, gf),
+    Internal(io, gf),
     timeIndex_(gf.timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -471,7 +495,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
 )
 :
-    DimensionedField<Type, GeoMesh>
+    Internal
     (
         io,
         const_cast<GeometricField<Type, PatchField, GeoMesh>&>(tgf()),
@@ -503,7 +527,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const GeometricField<Type, PatchField, GeoMesh>& gf
 )
 :
-    DimensionedField<Type, GeoMesh>(newName, gf),
+    Internal(newName, gf),
     timeIndex_(gf.timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -535,7 +559,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
 )
 :
-    DimensionedField<Type, GeoMesh>
+    Internal
     (
         newName,
         const_cast<GeometricField<Type, PatchField, GeoMesh>&>(tgf()),
@@ -566,7 +590,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const word& patchFieldType
 )
 :
-    DimensionedField<Type, GeoMesh>(io, gf),
+    Internal(io, gf),
     timeIndex_(gf.timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -602,7 +626,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 
 )
 :
-    DimensionedField<Type, GeoMesh>(io, gf),
+    Internal(io, gf),
     timeIndex_(gf.timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -644,7 +668,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const wordList& actualPatchTypes
 )
 :
-    DimensionedField<Type, GeoMesh>
+    Internal
     (
         io,
         const_cast<GeometricField<Type, PatchField, GeoMesh>&>(tgf()),
@@ -701,7 +725,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::dimensionedInternalFieldRef()
 template<class Type, template<class> class PatchField, class GeoMesh>
 typename
 Foam::GeometricField<Type, PatchField, GeoMesh>::Internal::FieldType&
-Foam::GeometricField<Type, PatchField, GeoMesh>::internalField()
+Foam::GeometricField<Type, PatchField, GeoMesh>::internalFieldRef()
 {
     this->setUpToDate();
     storeOldTimes();
@@ -992,7 +1016,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::T() const
         )
     );
 
-    Foam::T(result.ref().internalField(), internalField());
+    Foam::T(result.ref().internalFieldRef(), internalField());
     Foam::T(result.ref().boundaryFieldRef(), boundaryField());
 
     return result;
@@ -1029,7 +1053,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::component
         )
     );
 
-    Foam::component(Component.ref().internalField(), internalField(), d);
+    Foam::component(Component.ref().internalFieldRef(), internalField(), d);
     Foam::component(Component.ref().boundaryFieldRef(), boundaryField(), d);
 
     return Component;
@@ -1048,7 +1072,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::replace
      >& gcf
 )
 {
-    internalField().replace(d, gcf.internalField());
+    internalFieldRef().replace(d, gcf.internalField());
     boundaryFieldRef().replace(d, gcf.boundaryField());
 }
 
@@ -1060,7 +1084,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::replace
     const dimensioned<cmptType>& ds
 )
 {
-    internalField().replace(d, ds.value());
+    internalFieldRef().replace(d, ds.value());
     boundaryFieldRef().replace(d, ds.value());
 }
 
@@ -1071,7 +1095,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::max
     const dimensioned<Type>& dt
 )
 {
-    Foam::max(internalField(), internalField(), dt.value());
+    Foam::max(internalFieldRef(), internalField(), dt.value());
     Foam::max(boundaryFieldRef(), boundaryField(), dt.value());
 }
 
@@ -1082,7 +1106,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::min
     const dimensioned<Type>& dt
 )
 {
-    Foam::min(internalField(), internalField(), dt.value());
+    Foam::min(internalFieldRef(), internalField(), dt.value());
     Foam::min(boundaryFieldRef(), boundaryField(), dt.value());
 }
 
@@ -1090,7 +1114,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::min
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::negate()
 {
-    internalField().negate();
+    internalFieldRef().negate();
     boundaryFieldRef().negate();
 }
 
@@ -1141,7 +1165,7 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=
     this->dimensions() = gf.dimensions();
 
     // Transfer the storage from the tmp
-    internalField().transfer
+    internalFieldRef().transfer
     (
         const_cast<Field<Type>&>(gf.internalField())
     );
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
index 30b095e9b95740bd349e4e450b8bafe51b217995..bf8b97a1fccc3e651948e42f234dfa233f363a83 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
@@ -190,7 +190,7 @@ public:
             //- Read the boundary field
             void readField
             (
-                const DimensionedField<Type, GeoMesh>& field,
+                const Internal& field,
                 const dictionary& dict
             );
 
@@ -330,6 +330,14 @@ public:
             const wordList& actualPatchTypes = wordList()
         );
 
+        //- Constructor from components
+        GeometricField
+        (
+            const IOobject&,
+            const Internal&,
+            const PtrList<PatchField<Type>>&
+        );
+
         //- Constructor from components
         GeometricField
         (
@@ -452,7 +460,7 @@ public:
         //- Return a reference to the internal field
         //  Note: this increments the event counter and checks the
         //  old-time fields; avoid in loops.
-        typename Internal::FieldType& internalField();
+        typename Internal::FieldType& internalFieldRef();
 
         //- Return a const-reference to the  internal field
         inline const typename Internal::FieldType& internalField() const;
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
index dda073816170b507ce81c47ca44b0e3a652f0bbe..70e0e5f29b26b6838e79d386a199c3387f7495ea 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
@@ -49,7 +49,7 @@ void component
     const direction d
 )
 {
-    component(gcf.internalField(), gf.internalField(), d);
+    component(gcf.internalFieldRef(), gf.internalField(), d);
     component(gcf.boundaryFieldRef(), gf.boundaryField(), d);
 }
 
@@ -61,7 +61,7 @@ void T
      const GeometricField<Type, PatchField, GeoMesh>& gf1
 )
 {
-    T(gf.internalField(), gf1.internalField());
+    T(gf.internalFieldRef(), gf1.internalField());
     T(gf.boundaryFieldRef(), gf1.boundaryField());
 }
 
@@ -79,7 +79,7 @@ void pow
     const GeometricField<Type, PatchField, GeoMesh>& gf1
 )
 {
-    pow(gf.internalField(), gf1.internalField(), r);
+    pow(gf.internalFieldRef(), gf1.internalField(), r);
     pow(gf.boundaryFieldRef(), gf1.boundaryField(), r);
 }
 
@@ -173,7 +173,7 @@ void sqr
     const GeometricField<Type, PatchField, GeoMesh>& gf1
 )
 {
-    sqr(gf.internalField(), gf1.internalField());
+    sqr(gf.internalFieldRef(), gf1.internalField());
     sqr(gf.boundaryFieldRef(), gf1.boundaryField());
 }
 
@@ -261,7 +261,7 @@ void magSqr
     const GeometricField<Type, PatchField, GeoMesh>& gf
 )
 {
-    magSqr(gsf.internalField(), gf.internalField());
+    magSqr(gsf.internalFieldRef(), gf.internalField());
     magSqr(gsf.boundaryFieldRef(), gf.boundaryField());
 }
 
@@ -333,7 +333,7 @@ void mag
     const GeometricField<Type, PatchField, GeoMesh>& gf
 )
 {
-    mag(gsf.internalField(), gf.internalField());
+    mag(gsf.internalFieldRef(), gf.internalField());
     mag(gsf.boundaryFieldRef(), gf.boundaryField());
 }
 
@@ -410,7 +410,7 @@ void cmptAv
     const GeometricField<Type, PatchField, GeoMesh>& gf
 )
 {
-    cmptAv(gcf.internalField(), gf.internalField());
+    cmptAv(gcf.internalFieldRef(), gf.internalField());
     cmptAv(gcf.boundaryFieldRef(), gf.boundaryField());
 }
 
@@ -599,13 +599,18 @@ void opFunc                                                                    \
     const GeometricField<Type2, PatchField, GeoMesh>& gf2                      \
 )                                                                              \
 {                                                                              \
-    Foam::opFunc(gf.internalField(), gf1.internalField(), gf2.internalField());\
+    Foam::opFunc                                                               \
+    (                                                                          \
+        gf.internalFieldRef(),                                                 \
+        gf1.internalField(),                                                   \
+        gf2.internalField()                                                    \
+    );                                                                         \
     Foam::opFunc                                                               \
     (                                                                          \
         gf.boundaryFieldRef(),                                                 \
         gf1.boundaryField(),                                                   \
         gf2.boundaryField()                                                    \
-    );\
+    );                                                                         \
 }                                                                              \
                                                                                \
 template                                                                       \
@@ -750,7 +755,7 @@ void opFunc                                                                    \
     const dimensioned<Form>& dvs                                               \
 )                                                                              \
 {                                                                              \
-    Foam::opFunc(gf.internalField(), gf1.internalField(), dvs.value());        \
+    Foam::opFunc(gf.internalFieldRef(), gf1.internalField(), dvs.value());     \
     Foam::opFunc(gf.boundaryFieldRef(), gf1.boundaryField(), dvs.value());     \
 }                                                                              \
                                                                                \
@@ -863,7 +868,7 @@ void opFunc                                                                    \
     const GeometricField<Type, PatchField, GeoMesh>& gf1                       \
 )                                                                              \
 {                                                                              \
-    Foam::opFunc(gf.internalField(), dvs.value(), gf1.internalField());        \
+    Foam::opFunc(gf.internalFieldRef(), dvs.value(), gf1.internalField());     \
     Foam::opFunc(gf.boundaryFieldRef(), dvs.value(), gf1.boundaryField());     \
 }                                                                              \
                                                                                \
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C
index d044d677aa909f81cdd90b0c66354d50e62006b9..774b89c7538cf442f0066acbda142e1447cbf7f9 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C
@@ -41,7 +41,7 @@ void Func                                                                      \
     const GeometricField<Type1, PatchField, GeoMesh>& gf1                      \
 )                                                                              \
 {                                                                              \
-    Foam::Func(res.internalField(), gf1.internalField());                      \
+    Foam::Func(res.internalFieldRef(), gf1.internalField());                   \
     Foam::Func(res.boundaryFieldRef(), gf1.boundaryField());                   \
 }                                                                              \
                                                                                \
@@ -110,7 +110,7 @@ void OpFunc                                                                    \
     const GeometricField<Type1, PatchField, GeoMesh>& gf1                      \
 )                                                                              \
 {                                                                              \
-    Foam::OpFunc(res.internalField(), gf1.internalField());                    \
+    Foam::OpFunc(res.internalFieldRef(), gf1.internalField());                 \
     Foam::OpFunc(res.boundaryFieldRef(), gf1.boundaryField());                 \
 }                                                                              \
                                                                                \
@@ -180,13 +180,18 @@ void Func                                                                      \
     const GeometricField<Type2, PatchField, GeoMesh>& gf2                      \
 )                                                                              \
 {                                                                              \
-    Foam::Func(res.internalField(), gf1.internalField(), gf2.internalField()); \
+    Foam::Func                                                                 \
+    (                                                                          \
+        res.internalFieldRef(),                                                \
+        gf1.internalField(),                                                   \
+        gf2.internalField()                                                    \
+    );                                                                         \
     Foam::Func                                                                 \
     (                                                                          \
         res.boundaryFieldRef(),                                                \
         gf1.boundaryField(),                                                   \
         gf2.boundaryField()                                                    \
-    );\
+    );                                                                         \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -313,7 +318,7 @@ void Func                                                                      \
     const GeometricField<Type2, PatchField, GeoMesh>& gf2                      \
 )                                                                              \
 {                                                                              \
-    Foam::Func(res.internalField(), dt1.value(), gf2.internalField());         \
+    Foam::Func(res.internalFieldRef(), dt1.value(), gf2.internalField());      \
     Foam::Func(res.boundaryFieldRef(), dt1.value(), gf2.boundaryField());      \
 }                                                                              \
                                                                                \
@@ -404,7 +409,7 @@ void Func                                                                      \
     const dimensioned<Type2>& dt2                                              \
 )                                                                              \
 {                                                                              \
-    Foam::Func(res.internalField(), gf1.internalField(), dt2.value());         \
+    Foam::Func(res.internalFieldRef(), gf1.internalField(), dt2.value());      \
     Foam::Func(res.boundaryFieldRef(), gf1.boundaryField(), dt2.value());      \
 }                                                                              \
                                                                                \
@@ -503,9 +508,9 @@ void OpFunc                                                                    \
 )                                                                              \
 {                                                                              \
     Foam::OpFunc                                                               \
-        (res.internalField(), gf1.internalField(), gf2.internalField());       \
+    (res.internalFieldRef(), gf1.internalField(), gf2.internalField());        \
     Foam::OpFunc                                                               \
-        (res.boundaryFieldRef(), gf1.boundaryField(), gf2.boundaryField());    \
+    (res.boundaryFieldRef(), gf1.boundaryField(), gf2.boundaryField());        \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -632,7 +637,7 @@ void OpFunc                                                                    \
     const GeometricField<Type2, PatchField, GeoMesh>& gf2                      \
 )                                                                              \
 {                                                                              \
-    Foam::OpFunc(res.internalField(), dt1.value(), gf2.internalField());       \
+    Foam::OpFunc(res.internalFieldRef(), dt1.value(), gf2.internalField());    \
     Foam::OpFunc(res.boundaryFieldRef(), dt1.value(), gf2.boundaryField());    \
 }                                                                              \
                                                                                \
@@ -723,7 +728,7 @@ void OpFunc                                                                    \
     const dimensioned<Type2>& dt2                                              \
 )                                                                              \
 {                                                                              \
-    Foam::OpFunc(res.internalField(), gf1.internalField(), dt2.value());       \
+    Foam::OpFunc(res.internalFieldRef(), gf1.internalField(), dt2.value());    \
     Foam::OpFunc(res.boundaryFieldRef(), gf1.boundaryField(), dt2.value());    \
 }                                                                              \
                                                                                \
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H
index e465e00bed95236feb7e8639e01f05a8c2c78659..eac89f8142a0de80480f88dadd85a2dac5542cbb 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H
@@ -126,7 +126,7 @@ void MapGeometricFields
             // Map the internal field
             MapInternalField<Type, MeshMapper, GeoMesh>()
             (
-                field.internalField(),
+                field.internalFieldRef(),
                 mapper
             );
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C b/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
index 6b58138dc7c52a919258a0766057f88d15fd4c54..cb3ec9aff7adaede75e370951a5f30665eadb429 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
@@ -43,7 +43,7 @@ void stabilise
     const dimensioned<scalar>& ds
 )
 {
-    stabilise(result.internalField(), gsf.internalField(), ds.value());
+    stabilise(result.internalFieldRef(), gsf.internalField(), ds.value());
     stabilise(result.boundaryFieldRef(), gsf.boundaryField(), ds.value());
 }
 
@@ -126,7 +126,7 @@ void pow
     const GeometricField<scalar, PatchField, GeoMesh>& gsf2
 )
 {
-    pow(Pow.internalField(), gsf1.internalField(), gsf2.internalField());
+    pow(Pow.internalFieldRef(), gsf1.internalField(), gsf2.internalField());
     pow(Pow.boundaryFieldRef(), gsf1.boundaryField(), gsf2.boundaryField());
 }
 
@@ -269,7 +269,7 @@ void pow
     const dimensioned<scalar>& ds
 )
 {
-    pow(tPow.internalField(), gsf.internalField(), ds.value());
+    pow(tPow.internalFieldRef(), gsf.internalField(), ds.value());
     pow(tPow.boundaryFieldRef(), gsf.boundaryField(), ds.value());
 }
 
@@ -358,7 +358,7 @@ void pow
     const GeometricField<scalar, PatchField, GeoMesh>& gsf
 )
 {
-    pow(tPow.internalField(), ds.value(), gsf.internalField());
+    pow(tPow.internalFieldRef(), ds.value(), gsf.internalField());
     pow(tPow.boundaryFieldRef(), ds.value(), gsf.boundaryField());
 }
 
@@ -450,7 +450,7 @@ void atan2
     const GeometricField<scalar, PatchField, GeoMesh>& gsf2
 )
 {
-    atan2(Atan2.internalField(), gsf1.internalField(), gsf2.internalField());
+    atan2(Atan2.internalFieldRef(), gsf1.internalField(), gsf2.internalField());
     atan2(Atan2.boundaryFieldRef(), gsf1.boundaryField(), gsf2.boundaryField());
 }
 
@@ -577,7 +577,7 @@ void atan2
     const dimensioned<scalar>& ds
 )
 {
-    atan2(tAtan2.internalField(), gsf.internalField(), ds.value());
+    atan2(tAtan2.internalFieldRef(), gsf.internalField(), ds.value());
     atan2(tAtan2.boundaryFieldRef(), gsf.boundaryField(), ds.value());
 }
 
@@ -666,7 +666,7 @@ void atan2
     const GeometricField<scalar, PatchField, GeoMesh>& gsf
 )
 {
-    atan2(tAtan2.internalField(), ds.value(), gsf.internalField());
+    atan2(tAtan2.internalFieldRef(), ds.value(), gsf.internalField());
     atan2(tAtan2.boundaryFieldRef(), ds.value(), gsf.boundaryField());
 }
 
@@ -799,7 +799,7 @@ void func                                                                      \
     const GeometricField<scalar, PatchField, GeoMesh>& gsf1                    \
 )                                                                              \
 {                                                                              \
-    func(gsf.internalField(), n, gsf1.internalField());                        \
+    func(gsf.internalFieldRef(), n, gsf1.internalField());                     \
     func(gsf.boundaryFieldRef(), n, gsf1.boundaryField());                     \
 }                                                                              \
                                                                                \
diff --git a/src/OpenFOAM/fields/GeometricFields/transformGeometricField/transformGeometricField.C b/src/OpenFOAM/fields/GeometricFields/transformGeometricField/transformGeometricField.C
index baa2f7041b2c555271e211b914d89f39351e29fa..31c70df1c349b47c08b39f14d1458256a6a237a6 100644
--- a/src/OpenFOAM/fields/GeometricFields/transformGeometricField/transformGeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/transformGeometricField/transformGeometricField.C
@@ -45,7 +45,7 @@ void transform
     const GeometricField<Type, PatchField, GeoMesh>& tf
 )
 {
-    transform(rtf.internalField(), trf.internalField(), tf.internalField());
+    transform(rtf.internalFieldRef(), trf.internalField(), tf.internalField());
     transform(rtf.boundaryFieldRef(), trf.boundaryField(), tf.boundaryField());
 }
 
@@ -131,7 +131,7 @@ void transform
     const GeometricField<Type, PatchField, GeoMesh>& tf
 )
 {
-    transform(rtf.internalField(), t.value(), tf.internalField());
+    transform(rtf.internalFieldRef(), t.value(), tf.internalField());
     transform(rtf.boundaryFieldRef(), t.value(), tf.boundaryField());
 }
 
diff --git a/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/IDDESDelta/IDDESDelta.C b/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/IDDESDelta/IDDESDelta.C
index 9b6f2c4c8b1ea296399e77283b86f719fd764d51..8526b582802e9ace5fdd1dcddc0983f8aa42947e 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/IDDESDelta/IDDESDelta.C
+++ b/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/IDDESDelta/IDDESDelta.C
@@ -66,7 +66,7 @@ void Foam::LESModels::IDDESDelta::calcDelta()
         )
     );
 
-    scalarField& faceToFacenMax = tfaceToFacenMax.ref().internalField();
+    scalarField& faceToFacenMax = tfaceToFacenMax.ref().internalFieldRef();
 
     const cellList& cells = mesh.cells();
     const vectorField& faceCentres = mesh.faceCentres();
@@ -113,7 +113,7 @@ void Foam::LESModels::IDDESDelta::calcDelta()
             << "Case must be either 2D or 3D" << exit(FatalError);
     }
 
-    delta_.internalField() =
+    delta_.internalFieldRef() =
         min
         (
             max
diff --git a/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C b/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C
index 9829fed3c511294a2f4d1ebc9f9fc479d1bb7ffe..c322119d90e0a629f09eb05e919f65f2a176fccc 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C
+++ b/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.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
@@ -48,7 +48,7 @@ void Foam::LESModels::cubeRootVolDelta::calcDelta()
 
     if (nD == 3)
     {
-        delta_.internalField() = deltaCoeff_*pow(mesh.V(), 1.0/3.0);
+        delta_.internalFieldRef() = deltaCoeff_*pow(mesh.V(), 1.0/3.0);
     }
     else if (nD == 2)
     {
@@ -68,7 +68,7 @@ void Foam::LESModels::cubeRootVolDelta::calcDelta()
             }
         }
 
-        delta_.internalField() = deltaCoeff_*sqrt(mesh.V()/thickness);
+        delta_.internalFieldRef() = deltaCoeff_*sqrt(mesh.V()/thickness);
     }
     else
     {
diff --git a/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C b/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C
index a81e954677c503e0a5d111d0bbdf3024c11a1de0..4b0b708a1beeea1ace9482cef26532f5cde40c31 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C
+++ b/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C
@@ -71,7 +71,7 @@ void Foam::LESModels::maxDeltaxyz::calcDelta()
 
     if (nD == 3)
     {
-        delta_.internalField() = hmax;
+        delta_.internalFieldRef() = hmax;
     }
     else if (nD == 2)
     {
@@ -79,7 +79,7 @@ void Foam::LESModels::maxDeltaxyz::calcDelta()
             << "Case is 2D, LES is not strictly applicable\n"
             << endl;
 
-        delta_.internalField() = hmax;
+        delta_.internalFieldRef() = hmax;
     }
     else
     {
diff --git a/src/TurbulenceModels/turbulenceModels/LES/LESfilters/anisotropicFilter/anisotropicFilter.C b/src/TurbulenceModels/turbulenceModels/LES/LESfilters/anisotropicFilter/anisotropicFilter.C
index 0d326caa8cf4bbaf5bb48ac7b5a5739bc1e1f800..75dc5ef9e748a999fd82b7a1660f6305ed35449c 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/LESfilters/anisotropicFilter/anisotropicFilter.C
+++ b/src/TurbulenceModels/turbulenceModels/LES/LESfilters/anisotropicFilter/anisotropicFilter.C
@@ -63,7 +63,7 @@ Foam::anisotropicFilter::anisotropicFilter
 {
     for (direction d=0; d<vector::nComponents; d++)
     {
-        coeff_.internalField().replace
+        coeff_.internalFieldRef().replace
         (
             d,
             (1/widthCoeff_)*
@@ -100,7 +100,7 @@ Foam::anisotropicFilter::anisotropicFilter
 {
     for (direction d=0; d<vector::nComponents; d++)
     {
-        coeff_.internalField().replace
+        coeff_.internalFieldRef().replace
         (
             d,
             (1/widthCoeff_)*
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
index 4e05229233a6dfde5456179fb60224c6cd5ed14f..29b14966eafbad004f246911dccdf6a9482e9aa1 100644
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
@@ -344,7 +344,7 @@ Foam::dynamicRefineFvMesh::refine
 
                 surfaceScalarField& phi = *iter();
 
-                sigFpe::fillNan(phi.internalField());
+                sigFpe::fillNan(phi.internalFieldRef());
 
                 continue;
             }
diff --git a/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C b/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
index 42a9c0648e097b44e64bc5142b77ba533e190022..f5d9c5e821eb2aeb1679cf0478f09b8134461687 100644
--- a/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
+++ b/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
@@ -49,7 +49,7 @@ void Foam::fvMeshAdder::MapVolField
         Field<Type> oldInternalField(fld.internalField());
 
         // Modify internal field
-        Field<Type>& intFld = fld.internalField();
+        Field<Type>& intFld = fld.internalFieldRef();
 
         intFld.setSize(mesh.nCells());
 
@@ -340,7 +340,7 @@ void Foam::fvMeshAdder::MapSurfaceField
         Field<Type> oldField(fld);
 
         // Modify internal field
-        Field<Type>& intFld = fld.internalField();
+        Field<Type>& intFld = fld.internalFieldRef();
 
         intFld.setSize(mesh.nInternalFaces());
 
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
index f42443485676e4fc18d8e4ed5311c536febac472..a92fac8937f447c830815d4f2593a87f130299c0 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
@@ -319,7 +319,7 @@ void Foam::MRFZone::addCoriolis
     }
 
     const labelList& cells = mesh_.cellZones()[cellZoneID_];
-    vectorField& ddtUc = ddtU.internalField();
+    vectorField& ddtUc = ddtU.internalFieldRef();
     const vectorField& Uc = U;
 
     const vector Omega = this->Omega();
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
index 8d681da85c5e3acf41fdd43675f282e7f46837ed..b915f46ca371e155f8d6db2f2f23e249d278e28f 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
@@ -45,7 +45,7 @@ void Foam::MRFZone::makeRelativeRhoFlux
 
     const vectorField& Cfi = Cf;
     const vectorField& Sfi = Sf;
-    scalarField& phii = phi.internalField();
+    scalarField& phii = phi.internalFieldRef();
 
     // Internal faces
     forAll(internalFaces_, i)
@@ -145,7 +145,7 @@ void Foam::MRFZone::makeAbsoluteRhoFlux
 
     const vectorField& Cfi = Cf;
     const vectorField& Sfi = Sf;
-    scalarField& phii = phi.internalField();
+    scalarField& phii = phi.internalFieldRef();
 
     // Internal faces
     forAll(internalFaces_, i)
diff --git a/src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C b/src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C
index e19d96c44d2b35ed1d6f33db68a573fece391378..e8605ac57a0d95e395219bcd30db012652db9bea 100644
--- a/src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C
+++ b/src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C
@@ -228,7 +228,7 @@ Foam::tmp<Foam::volVectorField> Foam::SRF::SRFModel::Uabs() const
     volVectorField& Uabs = tUabs.ref();
 
     // Add SRF contribution to internal field
-    Uabs.internalField() += Urel_.internalField();
+    Uabs.internalFieldRef() += Urel_.internalField();
 
     // Add Urel boundary contributions
     volVectorField::Boundary& Uabsbf = Uabs.boundaryFieldRef();
diff --git a/src/finiteVolume/cfdTools/general/bound/bound.C b/src/finiteVolume/cfdTools/general/bound/bound.C
index 009a14388dc20a6a2a8a1a3c1de81b43e0881bd0..074a66f7a854fb8242a5368187f5bf9a2735d764 100644
--- a/src/finiteVolume/cfdTools/general/bound/bound.C
+++ b/src/finiteVolume/cfdTools/general/bound/bound.C
@@ -42,7 +42,7 @@ Foam::bound(volScalarField& vsf, const dimensionedScalar& lowerBound)
             << " average: " << gAverage(vsf.internalField())
             << endl;
 
-        vsf.internalField() = max
+        vsf.internalFieldRef() = max
         (
             max
             (
diff --git a/src/finiteVolume/cfdTools/general/include/volContinuity.H b/src/finiteVolume/cfdTools/general/include/volContinuity.H
index 68202f13b6d20764bae57a797665c1af5255e0d0..4d53eeebabac3b9ab75b522a9c5342c2e184ba06 100644
--- a/src/finiteVolume/cfdTools/general/include/volContinuity.H
+++ b/src/finiteVolume/cfdTools/general/include/volContinuity.H
@@ -3,7 +3,7 @@
     // The ddt term constructed by hand because it would be wrong for
     // Backward Differencing in time.
 
-    conserve.internalField() +=
+    conserve.internalFieldRef() +=
         (1.0 - mesh.V0()/mesh.V())/runTime.deltaTValue();
 
     scalar sumLocalContErr = runTime.deltaTValue()*
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
index 82bacaceabc17dba8dc3ce755b0145d5771f2518..baedf3de3e1964601d1d56ea9d0f2cc11fba9da6 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
@@ -179,7 +179,7 @@ CoEulerDdtScheme<Type>::fvcDdt
             )
         );
 
-        tdtdt.ref().internalField() =
+        tdtdt.ref().internalFieldRef() =
             rDeltaT.internalField()*dt.value()
            *(1.0 - mesh().Vsc0()/mesh().Vsc());
 
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
index 734e0d903f1e3d521f86e7f03cdace2f92f384cd..fed01ec469b725811e6429684e3fb3a9e686bc54 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
@@ -373,7 +373,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
 
-            ddt0.internalField() =
+            ddt0.internalFieldRef() =
             (
                 rDtCoef0*
                 (
@@ -466,7 +466,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
 
-            ddt0.internalField() =
+            ddt0.internalFieldRef() =
             (
                 rDtCoef0*rho.value()*
                 (
@@ -559,7 +559,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
 
-            ddt0.internalField() =
+            ddt0.internalFieldRef() =
             (
                 rDtCoef0*
                 (
@@ -663,7 +663,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
 
-            ddt0.internalField() =
+            ddt0.internalFieldRef() =
             (
                 rDtCoef0*
                 (
@@ -801,7 +801,7 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
 
-            ddt0.internalField() =
+            ddt0.internalFieldRef() =
             (
                 rDtCoef0*
                 (
@@ -886,7 +886,7 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
 
-            ddt0.internalField() =
+            ddt0.internalFieldRef() =
             (
                 rDtCoef0*rho.value()*
                 (
@@ -972,7 +972,7 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
 
-            ddt0.internalField() =
+            ddt0.internalFieldRef() =
             (
                 rDtCoef0*
                 (
@@ -1067,7 +1067,7 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
         {
             scalar rDtCoef0 = rDtCoef0_(ddt0).value();
 
-            ddt0.internalField() =
+            ddt0.internalFieldRef() =
             (
                 rDtCoef0*
                 (
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
index 5a72c9d93201f3147b33ce602fc038ed691be574..3e0805777290f5028de7e22562e0ae3fb898f276 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
@@ -73,7 +73,7 @@ EulerDdtScheme<Type>::fvcDdt
             )
         );
 
-        tdtdt.ref().internalField() =
+        tdtdt.ref().internalFieldRef() =
             rDeltaT.value()*dt.value()*(1.0 - mesh().Vsc0()/mesh().Vsc());
 
         return tdtdt;
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
index 12512767c7e4a18970219d56771f09f351fec9e6..dbe8ccb9d185dfd6aac11e500c400f0b4d8ba928 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
@@ -117,7 +117,7 @@ tmp<volScalarField> SLTSDdtScheme<Type>::SLrDeltaT() const
 
     if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
     {
-        rDeltaT.internalField() = max
+        rDeltaT.internalFieldRef() = max
         (
             rDeltaT.internalField()/mesh().V(),
             scalar(1)/deltaT.value()
@@ -131,7 +131,7 @@ tmp<volScalarField> SLTSDdtScheme<Type>::SLrDeltaT() const
                 rhoName_
             ).oldTime();
 
-        rDeltaT.internalField() = max
+        rDeltaT.internalFieldRef() = max
         (
             rDeltaT.internalField()/(rho.internalField()*mesh().V()),
             scalar(1)/deltaT.value()
@@ -183,7 +183,7 @@ SLTSDdtScheme<Type>::fvcDdt
             )
         );
 
-        tdtdt.ref().internalField() =
+        tdtdt.ref().internalFieldRef() =
             rDeltaT.internalField()*dt.value()*(1.0 - mesh().V0()/mesh().V());
 
         return tdtdt;
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
index a405a2c3f363641b699f9ab830a36b6fd3471c67..14279eea01a31607fede48eff8e06be88228fd28 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
@@ -111,7 +111,7 @@ backwardDdtScheme<Type>::fvcDdt
             )
         );
 
-        tdtdt.ref().internalField() = rDeltaT.value()*dt.value()*
+        tdtdt.ref().internalFieldRef() = rDeltaT.value()*dt.value()*
         (
             coefft - (coefft0*mesh().V0() - coefft00*mesh().V00())/mesh().V()
         );
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
index 5de6d069e9fb828dfcbfa77ce83ce7a7790f5207..4a4d59827b72a3d0be80f3aac08d6194d15366f2 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
@@ -79,7 +79,7 @@ localEulerDdtScheme<Type>::fvcDdt
             )
         );
 
-        tdtdt.ref().internalField() =
+        tdtdt.ref().internalFieldRef() =
             rDeltaT.internalField()*dt.value()
            *(1.0 - mesh().Vsc0()/mesh().Vsc());
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcAverage.C b/src/finiteVolume/finiteVolume/fvc/fvcAverage.C
index 7304e78573b4fece9d43189d8dcaa2838d24fd10..0fcc4b2fd38119553e5f409a87d9a0ad231a1dda 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcAverage.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcAverage.C
@@ -67,7 +67,7 @@ average
     );
     GeometricField<Type, fvPatchField, volMesh>& av = taverage.ref();
 
-    av.internalField() =
+    av.internalFieldRef() =
     (
         surfaceSum(mesh.magSf()*ssf)().internalField()
        /surfaceSum(mesh.magSf())().internalField()
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSurfaceIntegrate.C b/src/finiteVolume/finiteVolume/fvc/fvcSurfaceIntegrate.C
index 3ce8ab38ec6caa89c798bd5afffbe2bd717e31ac..34e5b37d1ae7ac202f73d7108c03261c38d27b44 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcSurfaceIntegrate.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSurfaceIntegrate.C
@@ -109,7 +109,7 @@ surfaceIntegrate
     );
     GeometricField<Type, fvPatchField, volMesh>& vf = tvf.ref();
 
-    surfaceIntegrate(vf.internalField(), ssf);
+    surfaceIntegrate(vf.internalFieldRef(), ssf);
     vf.correctBoundaryConditions();
 
     return tvf;
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C
index 6fd873c727cda9d082cd82b9f84c1c4cc3cec8aa..e35327363ecd09ff92f17410cd61dee0e89f922d 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C
@@ -183,7 +183,7 @@ Foam::fv::cellLimitedGrad<Foam::scalar>::calcGrad
             << " average: " << gAverage(limiter) << endl;
     }
 
-    g.internalField() *= limiter;
+    g.internalFieldRef() *= limiter;
     g.correctBoundaryConditions();
     gaussGrad<scalar>::correctBoundaryConditions(vsf, g);
 
@@ -336,7 +336,7 @@ Foam::fv::cellLimitedGrad<Foam::vector>::calcGrad
             << " average: " << gAverage(limiter) << endl;
     }
 
-    tensorField& gIf = g.internalField();
+    tensorField& gIf = g.internalFieldRef();
 
     forAll(gIf, celli)
     {
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C
index 40b1dcd1f7a6bdf14154c6c76cae8c150ccf6f8e..070987701cf270e962a8bae32af645a58dc49c61 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C
@@ -165,7 +165,7 @@ Foam::fv::faceLimitedGrad<Foam::scalar>::calcGrad
             << " average: " << gAverage(limiter) << endl;
     }
 
-    g.internalField() *= limiter;
+    g.internalFieldRef() *= limiter;
     g.correctBoundaryConditions();
     gaussGrad<scalar>::correctBoundaryConditions(vsf, g);
 
@@ -326,7 +326,7 @@ Foam::fv::faceLimitedGrad<Foam::vector>::calcGrad
             << " average: " << gAverage(limiter) << endl;
     }
 
-    g.internalField() *= limiter;
+    g.internalFieldRef() *= limiter;
     g.correctBoundaryConditions();
     gaussGrad<vector>::correctBoundaryConditions(vvf, g);
 
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C
index 4e7b8b83f95af8a69d147c7c5e74ba995aea312f..a621a62a4ada7acb8e887066d88a48ecaa15d1de 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C
@@ -68,7 +68,7 @@ Foam::fv::faceCorrectedSnGrad<Type>::fullGradCorrection
         )
     );
 
-    Field<Type>& sfCorr = tsfCorr.ref().internalField();
+    Field<Type>& sfCorr = tsfCorr.ref().internalFieldRef();
 
     const pointField& points = mesh.points();
     const faceList& faces = mesh.faces();
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 9abbb77eaa6b2c29afc3c48a10b3f81116bd1d32..2c54f9bac91bb86aac46de5e51fd980d211005b4 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -191,7 +191,7 @@ void Foam::fvMatrix<Type>::setValuesFromList
         const_cast
         <
             GeometricField<Type, fvPatchField, volMesh>&
-        >(psi_).internalField();
+        >(psi_).internalFieldRef();
 
     forAll(cellLabels, i)
     {
@@ -742,7 +742,7 @@ Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
         )
     );
 
-    tAphi.ref().internalField() = D()/psi_.mesh().V();
+    tAphi.ref().internalFieldRef() = D()/psi_.mesh().V();
     tAphi.ref().correctBoundaryConditions();
 
     return tAphi;
@@ -782,13 +782,13 @@ Foam::fvMatrix<Type>::H() const
         boundaryDiagCmpt.negate();
         addCmptAvBoundaryDiag(boundaryDiagCmpt);
 
-        Hphi.internalField().replace(cmpt, boundaryDiagCmpt*psiCmpt);
+        Hphi.internalFieldRef().replace(cmpt, boundaryDiagCmpt*psiCmpt);
     }
 
-    Hphi.internalField() += lduMatrix::H(psi_.internalField()) + source_;
-    addBoundarySource(Hphi.internalField());
+    Hphi.internalFieldRef() += lduMatrix::H(psi_.internalField()) + source_;
+    addBoundarySource(Hphi.internalFieldRef());
 
-    Hphi.internalField() /= psi_.mesh().V();
+    Hphi.internalFieldRef() /= psi_.mesh().V();
     Hphi.correctBoundaryConditions();
 
     typename Type::labelType validComponents
@@ -834,7 +834,7 @@ Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::H1() const
     );
     volScalarField& H1_ = tH1.ref();
 
-    H1_.internalField() = lduMatrix::H1();
+    H1_.internalFieldRef() = lduMatrix::H1();
 
     forAll(psi_.boundaryField(), patchi)
     {
@@ -851,7 +851,7 @@ Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::H1() const
         }
     }
 
-    H1_.internalField() /= psi_.mesh().V();
+    H1_.internalFieldRef() /= psi_.mesh().V();
     H1_.correctBoundaryConditions();
 
     return tH1;
@@ -894,7 +894,7 @@ flux() const
 
     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
     {
-        fieldFlux.internalField().replace
+        fieldFlux.internalFieldRef().replace
         (
             cmpt,
             lduMatrix::faceH(psi_.internalField().component(cmpt))
@@ -2289,18 +2289,18 @@ Foam::operator&
             scalarField psiCmpt(psi.field().component(cmpt));
             scalarField boundaryDiagCmpt(M.diag());
             M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
-            Mphi.internalField().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
+            Mphi.internalFieldRef().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
         }
     }
     else
     {
-        Mphi.internalField() = Zero;
+        Mphi.internalFieldRef() = Zero;
     }
 
-    Mphi.internalField() += M.lduMatrix::H(psi.field()) + M.source();
-    M.addBoundarySource(Mphi.internalField());
+    Mphi.internalFieldRef() += M.lduMatrix::H(psi.field()) + M.source();
+    M.addBoundarySource(Mphi.internalFieldRef());
 
-    Mphi.internalField() /= -psi.mesh().V();
+    Mphi.internalFieldRef() /= -psi.mesh().V();
     Mphi.correctBoundaryConditions();
 
     return tMphi;
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
index 991d899ad8374b3f5529c58ff624d8d9d1e1c62b..60b29fb78d46feee1bdec54c0b90018cbc9e1e2e 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
@@ -202,7 +202,7 @@ Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::solveSegregated
 
         solverPerfVec.replace(cmpt, solverPerf);
 
-        psi.internalField().replace(cmpt, psiCmpt);
+        psi.internalFieldRef().replace(cmpt, psiCmpt);
         diag() = saveDiag;
     }
 
diff --git a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
index 252ac0bb7937c4c2731681cb34f72cb4f23e39de..c72f271f849234ec6b91603accd3078367b7d302 100644
--- a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
@@ -114,7 +114,7 @@ Foam::solverPerformance Foam::fvMatrix<Foam::scalar>::fvSolver::solve
 
     solverPerformance solverPerf = solver_->solve
     (
-        psi.internalField(),
+        psi.internalFieldRef(),
         totalSource
     );
 
@@ -166,7 +166,7 @@ Foam::solverPerformance Foam::fvMatrix<Foam::scalar>::solveSegregated
         internalCoeffs_,
         psi_.boundaryField().scalarInterfaces(),
         solverControls
-    )->solve(psi.internalField(), totalSource);
+    )->solve(psi.internalFieldRef(), totalSource);
 
     if (solverPerformance::debug)
     {
@@ -229,10 +229,10 @@ Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Foam::scalar>::H() const
     );
     volScalarField& Hphi = tHphi.ref();
 
-    Hphi.internalField() = (lduMatrix::H(psi_.internalField()) + source_);
-    addBoundarySource(Hphi.internalField());
+    Hphi.internalFieldRef() = (lduMatrix::H(psi_.internalField()) + source_);
+    addBoundarySource(Hphi.internalFieldRef());
 
-    Hphi.internalField() /= psi_.mesh().V();
+    Hphi.internalFieldRef() /= psi_.mesh().V();
     Hphi.correctBoundaryConditions();
 
     return tHphi;
@@ -261,10 +261,10 @@ Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Foam::scalar>::H1() const
     );
     volScalarField& H1_ = tH1.ref();
 
-    H1_.internalField() = lduMatrix::H1();
+    H1_.internalFieldRef() = lduMatrix::H1();
     //addBoundarySource(Hphi.internalField());
 
-    H1_.internalField() /= psi_.mesh().V();
+    H1_.internalFieldRef() /= psi_.mesh().V();
     H1_.correctBoundaryConditions();
 
     return tH1;
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
index 6ecd89acf5fe5faaf6cb64191cfc25ad358e7690..38d71faf54f7a2ae523892626964a56331a66870 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
@@ -53,7 +53,7 @@ void Foam::MULES::correct
 
     if (mesh.moving())
     {
-        psi.internalField() =
+        psi.internalFieldRef() =
         (
             rho.field()*psi.internalField()*rDeltaT
           + Su.field()
@@ -62,7 +62,7 @@ void Foam::MULES::correct
     }
     else
     {
-        psi.internalField() =
+        psi.internalFieldRef() =
         (
             rho.field()*psi.internalField()*rDeltaT
           + Su.field()
diff --git a/src/finiteVolume/fvMesh/fvMesh.C b/src/finiteVolume/fvMesh/fvMesh.C
index 026a71a341942a5e79bb64e1345c6e8697caec29..31e5852cda01c0ec03066a031c7bffb90fb4c7e8 100644
--- a/src/finiteVolume/fvMesh/fvMesh.C
+++ b/src/finiteVolume/fvMesh/fvMesh.C
@@ -763,8 +763,8 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
     tmp<scalarField> tsweptVols = polyMesh::movePoints(p);
     scalarField& sweptVols = tsweptVols.ref();
 
-    phi.internalField() = scalarField::subField(sweptVols, nInternalFaces());
-    phi.internalField() *= rDeltaT;
+    phi.internalFieldRef() = scalarField::subField(sweptVols, nInternalFaces());
+    phi.internalFieldRef() *= rDeltaT;
 
     const fvPatchList& patches = boundary();
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
index 4614507065bd8f8517604eb030aa619c8c75e29e..ede6264d6e5cf34d48acbe7271c11faead9d62a6 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
@@ -57,7 +57,7 @@ void Foam::LimitedScheme<Type, Limiter, LimitFunc>::calcLimiter
 
     const vectorField& C = mesh.C();
 
-    scalarField& pLim = limiterField.internalField();
+    scalarField& pLim = limiterField.internalFieldRef();
 
     forAll(pLim, face)
     {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/PhiScheme/PhiScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/PhiScheme/PhiScheme.C
index 497f594abc109cb21bc30b01e2715821155812a6..77cdbcc3580f37f80a7bd4100a0a23c1aaef9902 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/PhiScheme/PhiScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/PhiScheme/PhiScheme.C
@@ -83,7 +83,7 @@ Foam::PhiScheme<Type, PhiLimiter>::limiter
 
     const surfaceScalarField& Uflux = tUflux();
 
-    scalarField& pLimiter = Limiter.internalField();
+    scalarField& pLimiter = Limiter.internalFieldRef();
 
     forAll(pLimiter, face)
     {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationScheme.C
index abce2de4f4ac87afb880f1bef8e44a4934ebdedc..184c8667a695c31403295798ac5cfeef9991a9a8 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationScheme.C
@@ -149,7 +149,7 @@ Foam::limitedSurfaceInterpolationScheme<Type>::weights
     // from which the weight is calculated using the limiter value
     surfaceScalarField& Weights = tLimiter.ref();
 
-    scalarField& pWeights = Weights.internalField();
+    scalarField& pWeights = Weights.internalFieldRef();
 
     forAll(pWeights, face)
     {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cellCoBlended/cellCoBlended.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cellCoBlended/cellCoBlended.H
index fbeb1eef1a7cad895061a1e3281652f3908f9385..3ba12ae60b655619b985e17d5f61a74c697285e0 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cellCoBlended/cellCoBlended.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/cellCoBlended/cellCoBlended.H
@@ -238,7 +238,7 @@ public:
                 fvc::surfaceSum(mag(tUflux))().internalField()
             );
 
-            Co.internalField() =
+            Co.internalFieldRef() =
                 (sumPhi/mesh.V().field())*(0.5*mesh.time().deltaTValue());
             Co.correctBoundaryConditions();
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/clippedLinear/clippedLinear.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/clippedLinear/clippedLinear.H
index 11201303687527cedf742d82cd301f6aacbc64cc..21fb2bdd8c9573cc43b28bf1dbcc86815b8186a2 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/clippedLinear/clippedLinear.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/clippedLinear/clippedLinear.H
@@ -153,7 +153,7 @@ public:
             surfaceScalarField& clippedLinearWeights =
                 tclippedLinearWeights.ref();
 
-            clippedLinearWeights.internalField() =
+            clippedLinearWeights.internalFieldRef() =
                 max(min(cdWeights.internalField(), 1 - wfLimit_), wfLimit_);
 
             surfaceScalarField::Boundary& clwbf =
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C
index 67d865b18b73ca2b430de95609e1a1d35ad14d0e..b1c274571076b5f7306745e8beb8cf89e580ccdb 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C
@@ -48,7 +48,7 @@ correction
     tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tsfCorr =
         linearInterpolate(vf);
 
-    Field<Type>& sfCorr = tsfCorr.ref().internalField();
+    Field<Type>& sfCorr = tsfCorr.ref().internalFieldRef();
 
     const pointField& points = mesh.points();
     const pointField& C = mesh.C();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/reverseLinear/reverseLinear.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/reverseLinear/reverseLinear.H
index 2e83f78c8d6494776a16bcd31cd5dea0dc0f4436..720fd82bb226324d3aa0f448623c2267369e807f 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/reverseLinear/reverseLinear.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/reverseLinear/reverseLinear.H
@@ -125,7 +125,7 @@ public:
             surfaceScalarField& reverseLinearWeights =
                 treverseLinearWeights.ref();
 
-            reverseLinearWeights.internalField() =
+            reverseLinearWeights.internalFieldRef() =
                 1.0 - cdWeights.internalField();
 
             surfaceScalarField::Boundary& rlwbf =
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
index 2d68f05c8805a30989108c7e7c7916d28b099205..c1291ace08d74ff60068414f7d682e6d92f16b89 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
@@ -169,7 +169,7 @@ void Foam::surfaceInterpolation::makeWeights() const
     const vectorField& Sf = mesh_.faceAreas();
 
     // ... and reference to the internal field of the weighting factors
-    scalarField& w = weights.internalField();
+    scalarField& w = weights.internalFieldRef();
 
     forAll(owner, facei)
     {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
index ed7b2576c60a9e28371507303fa8cf5eecdde7db..31e6a122739d4c8b0ecbe2a2ddf65d79e8f513f5 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
@@ -181,7 +181,7 @@ Foam::surfaceInterpolationScheme<Type>::interpolate
     );
     GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf.ref();
 
-    Field<Type>& sfi = sf.internalField();
+    Field<Type>& sfi = sf.internalFieldRef();
 
     for (label fi=0; fi<P.size(); fi++)
     {
@@ -274,7 +274,7 @@ Foam::surfaceInterpolationScheme<Type>::dotInterpolate
     );
     GeometricField<RetType, fvsPatchField, surfaceMesh>& sf = tsf.ref();
 
-    Field<RetType>& sfi = sf.internalField();
+    Field<RetType>& sfi = sf.internalFieldRef();
 
     const typename SFType::Internal& Sfi = Sf();
 
diff --git a/src/finiteVolume/interpolation/volPointInterpolation/pointConstraints.C b/src/finiteVolume/interpolation/volPointInterpolation/pointConstraints.C
index 480d93c97c2bef77e3b5bc2d03d20095cff25f0e..0c0b20d7c0bbe7a3123738cb86bc92125258e565 100644
--- a/src/finiteVolume/interpolation/volPointInterpolation/pointConstraints.C
+++ b/src/finiteVolume/interpolation/volPointInterpolation/pointConstraints.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
@@ -378,7 +378,7 @@ void Foam::pointConstraints::constrainDisplacement
     syncUntransformedData
     (
         pf.mesh()(),
-        pf.internalField(),
+        pf.internalFieldRef(),
         maxMagSqrEqOp<vector>()
     );
 
@@ -390,7 +390,7 @@ void Foam::pointConstraints::constrainDisplacement
     twoDPointCorrector::New(mesh()()).correctDisplacement
     (
         mesh()().points(),
-        pf.internalField()
+        pf.internalFieldRef()
     );
 
     if (overrideFixedValue)
diff --git a/src/finiteVolume/interpolation/volPointInterpolation/pointConstraintsTemplates.C b/src/finiteVolume/interpolation/volPointInterpolation/pointConstraintsTemplates.C
index ff87215ec9bd73c27902f5d450ed7218b5db6e29..4c8eb878bd01bfe2d9a68b277f149462f7e36974 100644
--- a/src/finiteVolume/interpolation/volPointInterpolation/pointConstraintsTemplates.C
+++ b/src/finiteVolume/interpolation/volPointInterpolation/pointConstraintsTemplates.C
@@ -137,7 +137,12 @@ void Foam::pointConstraints::constrain
     pf.correctBoundaryConditions();
 
     // Sync any dangling points
-    syncUntransformedData(mesh()(), pf.internalField(), maxMagSqrEqOp<Type>());
+    syncUntransformedData
+    (
+        mesh()(),
+        pf.internalFieldRef(),
+        maxMagSqrEqOp<Type>()
+    );
 
     // Apply multiple constraints on edge/corner points
     constrainCorners(pf);
diff --git a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
index 5620b989e7974c52f2a28610397445dc6da3205c..94408ef00c5c9e5490761de63af04e774184f13f 100644
--- a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
+++ b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
@@ -86,6 +86,9 @@ void Foam::volPointInterpolation::addSeparated
         Pout<< "volPointInterpolation::addSeparated" << endl;
     }
 
+    typename GeometricField<Type, pointPatchField, pointMesh>::
+        Internal& pfi = pf.dimensionedInternalFieldRef();
+
     typename GeometricField<Type, pointPatchField, pointMesh>::
         Boundary& pfbf = pf.boundaryFieldRef();
 
@@ -97,7 +100,7 @@ void Foam::volPointInterpolation::addSeparated
                 (pfbf[patchi]).initSwapAddSeparated
                 (
                     Pstream::nonBlocking,
-                    pf.internalField()
+                    pfi
                 );
         }
     }
@@ -113,7 +116,7 @@ void Foam::volPointInterpolation::addSeparated
                 (pfbf[patchi]).swapAddSeparated
                 (
                     Pstream::nonBlocking,
-                    pf.internalField()
+                    pfi
                 );
         }
     }
@@ -213,7 +216,7 @@ void Foam::volPointInterpolation::interpolateBoundaryField
 {
     const primitivePatch& boundary = boundaryPtr_();
 
-    Field<Type>& pfi = pf.internalField();
+    Field<Type>& pfi = pf.internalFieldRef();
 
     // Get face data in flat list
     tmp<Field<Type>> tboundaryVals(flatBoundaryField(vf));
diff --git a/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
index 396ef6f9026aa2c8401821c98ae238ffdd7b3639..b3e8befe3c0b00ed3e018c3ad298f7be6f58b2ad 100644
--- a/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
@@ -158,9 +158,9 @@ Foam::displacementComponentLaplacianFvMotionSolver::curPoints() const
 
         // Apply pointLocation_ b.c. to mesh points.
 
-        pointLocation_().internalField() = fvMesh_.points();
+        pointLocation_().internalFieldRef() = fvMesh_.points();
 
-        pointLocation_().internalField().replace
+        pointLocation_().internalFieldRef().replace
         (
             cmpt_,
             points0_ + pointDisplacement_.internalField()
@@ -181,7 +181,7 @@ Foam::displacementComponentLaplacianFvMotionSolver::curPoints() const
             }
         }
 
-        twoDCorrectPoints(pointLocation_().internalField());
+        twoDCorrectPoints(pointLocation_().internalFieldRef());
 
         return tmp<pointField>(pointLocation_().internalField());
     }
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
index 617193e0fad4373aa3defc6385a0f5023510bd6b..853fcd1362dc88260ed6c654477d3b98e66e50c6 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
@@ -173,7 +173,7 @@ Foam::displacementLaplacianFvMotionSolver::curPoints() const
                 << endl;
         }
 
-        pointLocation_().internalField() =
+        pointLocation_().internalFieldRef() =
             points0()
           + pointDisplacement_.internalField();
 
@@ -190,7 +190,7 @@ Foam::displacementLaplacianFvMotionSolver::curPoints() const
             }
         }
 
-        twoDCorrectPoints(pointLocation_().internalField());
+        twoDCorrectPoints(pointLocation_().internalFieldRef());
 
         return tmp<pointField>(pointLocation_().internalField());
     }
diff --git a/src/fvMotionSolver/motionDiffusivity/inverseDistance/inverseDistanceDiffusivity.C b/src/fvMotionSolver/motionDiffusivity/inverseDistance/inverseDistanceDiffusivity.C
index 657dac6d988b8f313c26199b8eec6d6c19287484..30e3d68e0d48d6fc07022c402d4fe16ddaa09b7e 100644
--- a/src/fvMotionSolver/motionDiffusivity/inverseDistance/inverseDistanceDiffusivity.C
+++ b/src/fvMotionSolver/motionDiffusivity/inverseDistance/inverseDistanceDiffusivity.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -100,7 +100,7 @@ void Foam::inverseDistanceDiffusivity::correct()
         dimless,
         zeroGradientFvPatchScalarField::typeName
     );
-    y_.internalField() = y();
+    y_.internalFieldRef() = y();
     y_.correctBoundaryConditions();
 
     faceDiffusivity_ = 1.0/fvc::interpolate(y_);
diff --git a/src/fvMotionSolver/motionDiffusivity/inverseVolume/inverseVolumeDiffusivity.C b/src/fvMotionSolver/motionDiffusivity/inverseVolume/inverseVolumeDiffusivity.C
index b768e5726ea3aaf976af9e073cf65f287d223667..aac4df700958a4ea3bd84c4282334a086aacf847 100644
--- a/src/fvMotionSolver/motionDiffusivity/inverseVolume/inverseVolumeDiffusivity.C
+++ b/src/fvMotionSolver/motionDiffusivity/inverseVolume/inverseVolumeDiffusivity.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,7 +85,7 @@ void Foam::inverseVolumeDiffusivity::correct()
         zeroGradientFvPatchScalarField::typeName
     );
 
-    V.internalField() = mesh().V();
+    V.internalFieldRef() = mesh().V();
     V.correctBoundaryConditions();
 
     faceDiffusivity_ = 1.0/fvc::interpolate(V);
diff --git a/src/fvOptions/corrections/limitTemperature/limitTemperature.C b/src/fvOptions/corrections/limitTemperature/limitTemperature.C
index 995471019ad6f8e469b7cc8a292c14aff741d06e..af7cbb409c6637875c9979b5d1fa7a0606d48a7a 100644
--- a/src/fvOptions/corrections/limitTemperature/limitTemperature.C
+++ b/src/fvOptions/corrections/limitTemperature/limitTemperature.C
@@ -84,7 +84,7 @@ void Foam::fv::limitTemperature::correct(volScalarField& he)
     scalarField heMin(thermo.he(thermo.p(), Tmin, cells_));
     scalarField heMax(thermo.he(thermo.p(), Tmax, cells_));
 
-    scalarField& hec = he.internalField();
+    scalarField& hec = he.internalFieldRef();
 
     forAll(cells_, i)
     {
diff --git a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
index 299fa07355e2b4268a984801ca4ae9f070a9679d..2dcbbb003a1f6ca3be69eee218de549c60690060 100644
--- a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
+++ b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
@@ -186,7 +186,7 @@ void Foam::fv::rotorDiskSource::writeField
             )
         );
 
-        Field<Type>& field = tfield.ref().internalField();
+        Field<Type>& field = tfield.ref().internalFieldRef();
 
         if (cells_.size() != values.size())
         {
diff --git a/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C b/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
index 8c51b741ecf347159828122f52695d0d4f4b7a62..41f13fcf5882d0259b4918f0189ded8b4484b218 100644
--- a/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
+++ b/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
@@ -162,7 +162,7 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
     (
         U.internalField(),
         plusEqOp<vector>(),
-        UNbr.internalField()
+        UNbr.internalFieldRef()
     );
 
     fvMatrix<vector> nbrEqn(UNbr, eqn.dimensions());
@@ -216,7 +216,7 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
     (
         U.internalField(),
         plusEqOp<vector>(),
-        UNbr.internalField()
+        UNbr.internalFieldRef()
     );
 
     fvMatrix<vector> nbrEqn(UNbr, eqn.dimensions());
@@ -257,7 +257,7 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
     (
         rho.internalField(),
         plusEqOp<scalar>(),
-        rhoNbr.internalField()
+        rhoNbr.internalFieldRef()
     );
 
     // Map local mu onto neighbour region
@@ -265,7 +265,7 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
     (
         mu.internalField(),
         plusEqOp<scalar>(),
-        muNbr.internalField()
+        muNbr.internalFieldRef()
     );
 
     porosityPtr_->addResistance(nbrEqn, rhoNbr, muNbr);
diff --git a/src/fvOptions/sources/interRegion/interRegionHeatTransfer/interRegionHeatTransferModel/interRegionHeatTransferModel.C b/src/fvOptions/sources/interRegion/interRegionHeatTransfer/interRegionHeatTransferModel/interRegionHeatTransferModel.C
index ee213ae246f4718ee07a742e656c43b32ff86388..fe4c425565ca728bd4346e57db4a60099380ce35 100644
--- a/src/fvOptions/sources/interRegion/interRegionHeatTransfer/interRegionHeatTransferModel/interRegionHeatTransferModel.C
+++ b/src/fvOptions/sources/interRegion/interRegionHeatTransfer/interRegionHeatTransferModel/interRegionHeatTransferModel.C
@@ -201,7 +201,7 @@ void Foam::fv::interRegionHeatTransferModel::addSup
     const volScalarField& Tnbr =
         nbrMesh.lookupObject<volScalarField>(TNbrName_);
 
-    interpolate(Tnbr, Tmapped.internalField());
+    interpolate(Tnbr, Tmapped.internalFieldRef());
 
     if (debug)
     {
diff --git a/src/fvOptions/sources/interRegion/interRegionHeatTransfer/tabulatedHeatTransfer/tabulatedHeatTransfer.C b/src/fvOptions/sources/interRegion/interRegionHeatTransfer/tabulatedHeatTransfer/tabulatedHeatTransfer.C
index 42003298365d5f5661dbe60be9db5ad980a5bd81..e97e6bcd22cb23cbe19931b329fcbc4dcd35aacf 100644
--- a/src/fvOptions/sources/interRegion/interRegionHeatTransfer/tabulatedHeatTransfer/tabulatedHeatTransfer.C
+++ b/src/fvOptions/sources/interRegion/interRegionHeatTransfer/tabulatedHeatTransfer/tabulatedHeatTransfer.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
@@ -122,7 +122,7 @@ void Foam::fv::tabulatedHeatTransfer::calculateHtc()
 
     const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
 
-    scalarField& htcc = htc_.internalField();
+    scalarField& htcc = htc_.internalFieldRef();
 
     forAll(htcc, i)
     {
diff --git a/src/fvOptions/sources/interRegion/interRegionHeatTransfer/variableHeatTransfer/variableHeatTransfer.C b/src/fvOptions/sources/interRegion/interRegionHeatTransfer/variableHeatTransfer/variableHeatTransfer.C
index 6b04ed17e702e4783c3230a6adb749efcc78e2b1..4f78964ed49b7c0ff04d1841198d7fa13a1bdabc 100644
--- a/src/fvOptions/sources/interRegion/interRegionHeatTransfer/variableHeatTransfer/variableHeatTransfer.C
+++ b/src/fvOptions/sources/interRegion/interRegionHeatTransfer/variableHeatTransfer/variableHeatTransfer.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
@@ -122,7 +122,7 @@ void Foam::fv::variableHeatTransfer::calculateHtc()
 
     const scalarField htcNbrMapped(interpolate(htcNbr));
 
-    htc_.internalField() = htcNbrMapped*AoV_;
+    htc_.internalFieldRef() = htcNbrMapped*AoV_;
 }
 
 
diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
index 91a6c48af5b656d6cf19939d790534c81eee42db..991e02725679948bd35edb0acd73c5fcfaa94a6a 100644
--- a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
+++ b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
@@ -202,7 +202,7 @@ void Foam::DSMCCloud<ParcelType>::initialise
         mostAbundantType
     );
 
-    sigmaTcRMax_.internalField() = cP.sigmaT()*maxwellianMostProbableSpeed
+    sigmaTcRMax_.internalFieldRef() = cP.sigmaT()*maxwellianMostProbableSpeed
     (
         temperature,
         cP.mass()
@@ -414,13 +414,13 @@ void Foam::DSMCCloud<ParcelType>::resetFields()
 template<class ParcelType>
 void Foam::DSMCCloud<ParcelType>::calculateFields()
 {
-    scalarField& rhoN = rhoN_.internalField();
-    scalarField& rhoM = rhoM_.internalField();
-    scalarField& dsmcRhoN = dsmcRhoN_.internalField();
-    scalarField& linearKE = linearKE_.internalField();
-    scalarField& internalE = internalE_.internalField();
-    scalarField& iDof = iDof_.internalField();
-    vectorField& momentum = momentum_.internalField();
+    scalarField& rhoN = rhoN_.internalFieldRef();
+    scalarField& rhoM = rhoM_.internalFieldRef();
+    scalarField& dsmcRhoN = dsmcRhoN_.internalFieldRef();
+    scalarField& linearKE = linearKE_.internalFieldRef();
+    scalarField& internalE = internalE_.internalFieldRef();
+    scalarField& iDof = iDof_.internalFieldRef();
+    vectorField& momentum = momentum_.internalFieldRef();
 
     forAllConstIter(typename DSMCCloud<ParcelType>, *this, iter)
     {
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index 339f49b1bb8a8e7d8eb704c9749019507299f8e5..b0a38037270b083adc5a0d74e372f615e467e5a3 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -483,7 +483,7 @@ Foam::KinematicCloud<CloudType>::vDotSweep() const
         vDotSweep[celli] += p.nParticle()*p.areaP()*mag(p.U() - U_[celli]);
     }
 
-    vDotSweep.internalField() /= mesh_.V();
+    vDotSweep.internalFieldRef() /= mesh_.V();
     vDotSweep.correctBoundaryConditions();
 
     return tvDotSweep;
@@ -522,7 +522,7 @@ Foam::KinematicCloud<CloudType>::theta() const
         theta[celli] += p.nParticle()*p.volume();
     }
 
-    theta.internalField() /= mesh_.V();
+    theta.internalFieldRef() /= mesh_.V();
     theta.correctBoundaryConditions();
 
     return ttheta;
@@ -551,7 +551,7 @@ Foam::KinematicCloud<CloudType>::alpha() const
         )
     );
 
-    scalarField& alpha = talpha.ref().internalField();
+    scalarField& alpha = talpha.ref().internalFieldRef();
     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
     {
         const parcelType& p = iter();
@@ -588,7 +588,7 @@ Foam::KinematicCloud<CloudType>::rhoEff() const
         )
     );
 
-    scalarField& rhoEff = trhoEff.ref().internalField();
+    scalarField& rhoEff = trhoEff.ref().internalFieldRef();
     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
     {
         const parcelType& p = iter();
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
index 1c9907ecdb3702c48ea0982143789a96d802eba8..141a8d4515226919c4346e20b7a9848448801e03 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
@@ -129,7 +129,7 @@ inline Foam::tmp<Foam::fvScalarMatrix> Foam::ReactingCloud<CloudType>::SYi
 
             volScalarField& sourceField = trhoTrans.ref();
 
-            sourceField.internalField() =
+            sourceField.internalFieldRef() =
                 rhoTrans_[i]/(this->db().time().deltaTValue()*this->mesh().V());
 
             const dimensionedScalar YiSMALL("YiSMALL", dimless, SMALL);
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
index d23371cb8213b3b9d03e3bed796e9ad560b24920..a3859eeb2a054e129517a63d48ac84f0de299f93 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
@@ -282,7 +282,7 @@ inline Foam::tmp<Foam::volScalarField> Foam::ThermoCloud<CloudType>::Ep() const
 
     if (radiation_)
     {
-        scalarField& Ep = tEp.ref().internalField();
+        scalarField& Ep = tEp.ref().internalFieldRef();
         const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
@@ -318,7 +318,7 @@ inline Foam::tmp<Foam::volScalarField> Foam::ThermoCloud<CloudType>::ap() const
 
     if (radiation_)
     {
-        scalarField& ap = tap.ref().internalField();
+        scalarField& ap = tap.ref().internalFieldRef();
         const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
@@ -355,7 +355,7 @@ Foam::ThermoCloud<CloudType>::sigmap() const
 
     if (radiation_)
     {
-        scalarField& sigmap = tsigmap.ref().internalField();
+        scalarField& sigmap = tsigmap.ref().internalFieldRef();
         const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
index 9310d99ff2358398a3daac054b9672451eaeeaff..25f526ae82eba2896dd966ad553547c3d90bf488 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
@@ -131,7 +131,7 @@ void Foam::ParticleErosion<CloudType>::preEvolve()
 {
     if (QPtr_.valid())
     {
-        QPtr_->internalField() = 0.0;
+        QPtr_->internalFieldRef() = 0.0;
     }
     else
     {
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C
index 01ecb354808837bb1874b9e28df0c108d1555dcb..9108677236f76110177845963609a755dfe0672a 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C
@@ -82,7 +82,7 @@ void Foam::VoidFraction<CloudType>::preEvolve()
 {
     if (thetaPtr_.valid())
     {
-        thetaPtr_->internalField() = 0.0;
+        thetaPtr_->internalFieldRef() = 0.0;
     }
     else
     {
@@ -115,7 +115,7 @@ void Foam::VoidFraction<CloudType>::postEvolve()
 
     const fvMesh& mesh = this->owner().mesh();
 
-    theta.internalField() /= mesh.time().deltaTValue()*mesh.V();
+    theta.internalFieldRef() /= mesh.time().deltaTValue()*mesh.V();
 
     CloudFunctionObject<CloudType>::postEvolve();
 }
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C
index b1df842fbe2735a6b64927df6c211aa37fafbd23..a1b8ae028e78bc46ab1d115fe51c45e22fbb284a 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C
@@ -231,10 +231,10 @@ bool Foam::AveragingMethod<Type>::write() const
     }
 
     // average
-    cellValue.internalField() /= mesh_.V();
-    cellGrad.internalField() /= mesh_.V();
-    pointValue.internalField() /= pointVolume;
-    pointGrad.internalField() /= pointVolume;
+    cellValue.internalFieldRef() /= mesh_.V();
+    cellGrad.internalFieldRef() /= mesh_.V();
+    pointValue.internalFieldRef() /= pointVolume;
+    pointGrad.internalFieldRef() /= pointVolume;
 
     // write
     if (!cellValue.write()) return false;
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C
index b921b0042d46c14ee830acdbee0a9aceb23c1593..126181accfd4fa3dcfb64711872c076439e68f4a 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C
@@ -80,7 +80,7 @@ void Foam::AveragingMethods::Basic<Type>::updateGrad()
         dimensioned<Type>("zero", dimless, Zero),
         zeroGradientFvPatchField<Type>::typeName
     );
-    tempData.internalField() = data_;
+    tempData.internalFieldRef() = data_;
     tempData.correctBoundaryConditions();
     dataGrad_ = fvc::grad(tempData)->internalField();
 }
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
index 6aaf8024dfde5eaac1b5d168342dc8ceb406d038..f95fc10f0f926db6742ccf5893b6d5e3da797fd4 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
@@ -150,7 +150,7 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
             dimensionedScalar("zero", dimDensity, 0),
             zeroGradientFvPatchField<scalar>::typeName
         );
-        rho.internalField() = max(rhoAverage.internalField(), rhoMin_);
+        rho.internalFieldRef() = max(rhoAverage.internalField(), rhoMin_);
         rho.correctBoundaryConditions();
 
         // Stress field
@@ -172,7 +172,7 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
             zeroGradientFvPatchField<scalar>::typeName
         );
 
-        tauPrime.internalField() =
+        tauPrime.internalFieldRef() =
             this->particleStressModel_->dTaudTheta
             (
                 alpha_.internalField(),
@@ -256,7 +256,7 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
                 dimensionedVector("zero", dimVelocity, Zero),
                 fixedValueFvPatchField<vector>::typeName
             );
-            U.internalField() = uAverage.internalField();
+            U.internalFieldRef() = uAverage.internalField();
             U.correctBoundaryConditions();
 
             surfaceScalarField phi
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
index 5c8fde3a124613c75762868e45191daa987cfa7e..457a6ed864df8d002e7993fa97d4d1bfa122e1bf 100644
--- a/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
@@ -8,7 +8,7 @@ if (runTime.outputTime())
 
     forAll(allSpeciesRhoN, rN)
     {
-        allSpeciesRhoN[rN].internalField() =
+        allSpeciesRhoN[rN].internalFieldRef() =
             allSpeciesN_RU[rN]
             /mesh.cellVolumes()
             /nAveragingSteps;
@@ -16,7 +16,7 @@ if (runTime.outputTime())
         totalRhoN_sum += allSpeciesRhoN[rN].internalField();
     }
 
-    totalRhoN.internalField() = totalRhoN_sum;
+    totalRhoN.internalFieldRef() = totalRhoN_sum;
 
 
     /*-----------------------------------------------------------------------*\
@@ -27,7 +27,7 @@ if (runTime.outputTime())
 
     forAll(allSpeciesRhoM, rM)
     {
-        allSpeciesRhoM[rM].internalField() =
+        allSpeciesRhoM[rM].internalFieldRef() =
             allSpeciesM_RU[rM]
             /mesh.cellVolumes()
             /nAveragingSteps;
@@ -35,7 +35,7 @@ if (runTime.outputTime())
         totalRhoM_sum += allSpeciesRhoM[rM].internalField();
     }
 
-    totalRhoM.internalField() = totalRhoM_sum;
+    totalRhoM.internalFieldRef() = totalRhoM_sum;
 
     /*-----------------------------------------------------------------------*\
         Bulk velocity
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H
index 83bcaa8ec754d6f0f16cbaa027b621e42a8ced43..d999d28cda849958d14dc5e4068f4103856ae60f 100644
--- a/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H
@@ -3,25 +3,25 @@
 List<scalarField> allSpeciesN_RU
 (
     molecules.potential().nIds(),
-    scalarField (mesh.nCells(), 0.0)
+    scalarField(mesh.nCells(), 0.0)
 );
 
 List<scalarField> allSpeciesM_RU
 (
     molecules.potential().nIds(),
-    scalarField (mesh.nCells(), 0.0)
+    scalarField(mesh.nCells(), 0.0)
 );
 
 List<vectorField> allSpeciesVelocitySum_RU
 (
     molecules.potential().nIds(),
-    vectorField (mesh.nCells(), Zero)
+    vectorField(mesh.nCells(), Zero)
 );
 
 List<scalarField> allSpeciesVelocityMagSquaredSum_RU
 (
     molecules.potential().nIds(),
-    scalarField (mesh.nCells(), 0.0)
+    scalarField(mesh.nCells(), 0.0)
 );
 
 // Geometric Fields for IO
@@ -60,7 +60,7 @@ forAll(allSpeciesRhoN, rN)
             "zeroGradient"
         )
     );
-    allSpeciesRhoN[rN].internalField() = scalarField (mesh.nCells(), 0.0);
+    allSpeciesRhoN[rN].internalFieldRef() = 0;
     allSpeciesRhoN[rN].correctBoundaryConditions();
 }
 
@@ -80,7 +80,7 @@ volScalarField totalRhoN
     dimless/dimVolume,
     "zeroGradient"
 );
-totalRhoN.internalField() = scalarField (mesh.nCells(), 0.0);
+totalRhoN.internalFieldRef() = 0;
 totalRhoN.correctBoundaryConditions();
 
 /*---------------------------------------------------------------------------*\
@@ -115,7 +115,7 @@ forAll(allSpeciesRhoM, rM)
             "zeroGradient"
         )
     );
-    allSpeciesRhoM[rM].internalField() = scalarField (mesh.nCells(), 0.0);
+    allSpeciesRhoM[rM].internalFieldRef() = 0;
     allSpeciesRhoM[rM].correctBoundaryConditions();
 }
 
@@ -135,7 +135,7 @@ volScalarField totalRhoM
     dimMass/dimVolume,
     "zeroGradient"
 );
-totalRhoM.internalField() = scalarField (mesh.nCells(), 0.0);
+totalRhoM.internalFieldRef() = 0;
 totalRhoM.correctBoundaryConditions();
 
 /*---------------------------------------------------------------------------*\
@@ -170,8 +170,7 @@ forAll(allSpeciesVelocity, v)
             "zeroGradient"
         )
     );
-    allSpeciesVelocity[v].internalField() =
-        vectorField (mesh.nCells(), Zero);
+    allSpeciesVelocity[v].internalFieldRef() = Zero;
     allSpeciesVelocity[v].correctBoundaryConditions();
 }
 
@@ -191,7 +190,7 @@ Info << "    Creating total velocity field" << endl;
 //     dimVelocity,
 //     "zeroGradient"
 // );
-// totalVelocity.internalField() = vectorField (mesh.nCells(), Zero);
+// totalVelocity.internalFieldRef() = Zero;
 // totalVelocity.correctBoundaryConditions();
 
 
@@ -244,7 +243,7 @@ forAll(allSpeciesTemperature, t)
             "zeroGradient"
         )
     );
-    allSpeciesTemperature[t].internalField() = scalarField (mesh.nCells(), 0.0);
+    allSpeciesTemperature[t].internalFieldRef() = 0;
     allSpeciesTemperature[t].correctBoundaryConditions();
 }
 
@@ -264,7 +263,7 @@ volScalarField totalTemperature
     dimTemperature,
     "zeroGradient"
 );
-totalTemperature.internalField() = scalarField (mesh.nCells(), 0.0);
+totalTemperature.internalFieldRef() = 0;
 totalTemperature.correctBoundaryConditions();
 
 /*---------------------------------------------------------------------------*\
@@ -300,7 +299,7 @@ forAll(allSpeciesMeanKE, mKE)
             "zeroGradient"
         )
     );
-    allSpeciesMeanKE[mKE].internalField() = scalarField (mesh.nCells(), 0.0);
+    allSpeciesMeanKE[mKE].internalFieldRef() = 0;
     allSpeciesMeanKE[mKE].correctBoundaryConditions();
 }
 
@@ -320,5 +319,5 @@ volScalarField totalMeanKE
     dimensionSet(1, 2, -2, 0, 0, 0, 0),
     "zeroGradient"
 );
-totalMeanKE.internalField() = scalarField (mesh.nCells(), 0.0);
+totalMeanKE.internalFieldRef() = 0;
 totalMeanKE.correctBoundaryConditions();
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/constantFilmThermo/constantFilmThermo.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/constantFilmThermo/constantFilmThermo.C
index 3e7602868b52223b195129caa8016fb1396786c9..d8e4e78f7dabbd5af58d4d4975f5babc267a95d5 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/constantFilmThermo/constantFilmThermo.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/constantFilmThermo/constantFilmThermo.C
@@ -275,7 +275,7 @@ tmp<volScalarField> constantFilmThermo::rho() const
         )
     );
 
-    trho.ref().internalField() = this->rho(0, 0);
+    trho.ref().internalFieldRef() = this->rho(0, 0);
     trho.ref().correctBoundaryConditions();
 
     return trho;
@@ -302,7 +302,7 @@ tmp<volScalarField> constantFilmThermo::mu() const
         )
     );
 
-    tmu.ref().internalField() = this->mu(0, 0);
+    tmu.ref().internalFieldRef() = this->mu(0, 0);
     tmu.ref().correctBoundaryConditions();
 
     return tmu;
@@ -329,7 +329,7 @@ tmp<volScalarField> constantFilmThermo::sigma() const
         )
     );
 
-    tsigma.ref().internalField() = this->sigma(0, 0);
+    tsigma.ref().internalFieldRef() = this->sigma(0, 0);
     tsigma.ref().correctBoundaryConditions();
 
     return tsigma;
@@ -356,7 +356,7 @@ tmp<volScalarField> constantFilmThermo::Cp() const
         )
     );
 
-    tCp.ref().internalField() = this->Cp(0, 0);
+    tCp.ref().internalFieldRef() = this->Cp(0, 0);
     tCp.ref().correctBoundaryConditions();
 
     return tCp;
@@ -383,7 +383,7 @@ tmp<volScalarField> constantFilmThermo::kappa() const
         )
     );
 
-    tkappa.ref().internalField() = this->kappa(0, 0);
+    tkappa.ref().internalFieldRef() = this->kappa(0, 0);
     tkappa.ref().correctBoundaryConditions();
 
     return tkappa;
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/liquidFilmThermo/liquidFilmThermo.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/liquidFilmThermo/liquidFilmThermo.C
index f4b1fbad9416cfc2491c9302263d4d7b38ed1c4b..55b21144b1a817cd4528bb2255fb75e8ff081a4d 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/liquidFilmThermo/liquidFilmThermo.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/liquidFilmThermo/liquidFilmThermo.C
@@ -255,7 +255,7 @@ tmp<volScalarField> liquidFilmThermo::rho() const
         )
     );
 
-    scalarField& rho = trho.ref().internalField();
+    scalarField& rho = trho.ref().internalFieldRef();
 
     if (useReferenceValues_)
     {
@@ -303,7 +303,7 @@ tmp<volScalarField> liquidFilmThermo::mu() const
         )
     );
 
-    scalarField& mu = tmu.ref().internalField();
+    scalarField& mu = tmu.ref().internalFieldRef();
 
     if (useReferenceValues_)
     {
@@ -351,7 +351,7 @@ tmp<volScalarField> liquidFilmThermo::sigma() const
         )
     );
 
-    scalarField& sigma = tsigma.ref().internalField();
+    scalarField& sigma = tsigma.ref().internalFieldRef();
 
     if (useReferenceValues_)
     {
@@ -399,7 +399,7 @@ tmp<volScalarField> liquidFilmThermo::Cp() const
         )
     );
 
-    scalarField& Cp = tCp.ref().internalField();
+    scalarField& Cp = tCp.ref().internalFieldRef();
 
     if (useReferenceValues_)
     {
@@ -447,7 +447,7 @@ tmp<volScalarField> liquidFilmThermo::kappa() const
         )
     );
 
-    scalarField& kappa = tkappa.ref().internalField();
+    scalarField& kappa = tkappa.ref().internalFieldRef();
 
     if (useReferenceValues_)
     {
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C
index 45c007087a553c5a1b53510cdd60a082efd533d2..57352ab7fb29d7974bda1a12583a68c9b6ad9053 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C
@@ -153,7 +153,7 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
         )
     );
 
-    vectorField& force = tForce.ref().internalField();
+    vectorField& force = tForce.ref().internalFieldRef();
 
     const labelUList& own = owner_.regionMesh().owner();
     const labelUList& nbr = owner_.regionMesh().neighbour();
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
index 5c6d904231f1fa21223b4ae77d72699ac25a59d2..81ed5ec25802beb33442792f3c80fd7d03996b8f 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
@@ -78,7 +78,7 @@ tmp<volScalarField> curvatureSeparation::calcInvR1
     );
 
 
-    scalarField& invR1 = tinvR1.ref().internalField();
+    scalarField& invR1 = tinvR1.ref().internalFieldRef();
 
     // apply defined patch radii
     const scalar rMin = 1e-6;
@@ -207,7 +207,7 @@ tmp<scalarField> curvatureSeparation::calcCosAngle
             dimensionedScalar("zero", dimless, 0.0),
             zeroGradientFvPatchScalarField::typeName
         );
-        volCosAngle.internalField() = cosAngle;
+        volCosAngle.internalFieldRef() = cosAngle;
         volCosAngle.correctBoundaryConditions();
         volCosAngle.write();
     }
@@ -349,7 +349,7 @@ void curvatureSeparation::correct
             dimensionedScalar("zero", dimForce, 0.0),
             zeroGradientFvPatchScalarField::typeName
         );
-        volFnet.internalField() = Fnet;
+        volFnet.internalFieldRef() = Fnet;
         volFnet.correctBoundaryConditions();
         volFnet.write();
     }
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
index ab51ea57e8b24376d8e322606b3079cb2595b12a..f96aa0ce9307f8304edc2e998f685cf6e09e0add 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
@@ -132,7 +132,7 @@ tmp<volScalarField> standardRadiation::Shs()
     Shs = beta_*QinP*alpha*(1.0 - exp(-kappaBar_*delta));
 
     // Update net Qr on local region
-    QrNet_.internalField() = QinP - Shs;
+    QrNet_.internalFieldRef() = QinP - Shs;
     QrNet_.correctBoundaryConditions();
 
     return tShs;
diff --git a/src/rigidBodyMeshMotion/rigidBodyMeshMotion.C b/src/rigidBodyMeshMotion/rigidBodyMeshMotion.C
index a207cfa0c00664d7bfdfd7994f72c289c7dfa6db..27f3333291fa6bb7c0f9dff34027dd0a3ec62f3a 100644
--- a/src/rigidBodyMeshMotion/rigidBodyMeshMotion.C
+++ b/src/rigidBodyMeshMotion/rigidBodyMeshMotion.C
@@ -153,7 +153,7 @@ Foam::rigidBodyMeshMotion::rigidBodyMeshMotion
         pointScalarField& scale = bodyMeshes_[bi].weight_;
 
         // Scaling: 1 up to di then linear down to 0 at do away from patches
-        scale.internalField() =
+        scale.internalFieldRef() =
             min
             (
                 max
@@ -166,7 +166,7 @@ Foam::rigidBodyMeshMotion::rigidBodyMeshMotion
             );
 
         // Convert the scale function to a cosine
-        scale.internalField() =
+        scale.internalFieldRef() =
             min
             (
                 max
@@ -281,7 +281,7 @@ void Foam::rigidBodyMeshMotion::solve()
     // Update the displacements
     if (bodyMeshes_.size() == 1)
     {
-        pointDisplacement_.internalField() = model_.transformPoints
+        pointDisplacement_.internalFieldRef() = model_.transformPoints
         (
             bodyMeshes_[0].bodyID_,
             bodyMeshes_[0].weight_,
@@ -298,7 +298,7 @@ void Foam::rigidBodyMeshMotion::solve()
             weights[bi] = &bodyMeshes_[bi].weight_;
         }
 
-        pointDisplacement_.internalField() =
+        pointDisplacement_.internalFieldRef() =
             model_.transformPoints(bodyIDs, weights, points0()) - points0();
     }
 
diff --git a/src/sampling/meshToMesh/meshToMeshTemplates.C b/src/sampling/meshToMesh/meshToMeshTemplates.C
index 11ca92fe7176a73f6f63b66de86499e40cd9ab56..000aa35ea58ecb45ac5a582d768ccf4667318a23 100644
--- a/src/sampling/meshToMesh/meshToMeshTemplates.C
+++ b/src/sampling/meshToMesh/meshToMeshTemplates.C
@@ -325,7 +325,7 @@ void Foam::meshToMesh::mapSrcToTgt
     GeometricField<Type, fvPatchField, volMesh>& result
 ) const
 {
-    mapSrcToTgt(field, cop, result.internalField());
+    mapSrcToTgt(field, cop, result.internalFieldRef());
 
     const PtrList<AMIPatchToPatchInterpolation>& AMIList = patchAMIs();
 
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
index 9b30a4827bd928804cfa4df869f51acd0c6de80d..91ccabd825f79a2d266e3cde49c27a867fc2430f 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
@@ -83,7 +83,7 @@ void Foam::distanceSurface::createGeometry()
     // Internal field
     {
         const pointField& cc = fvm.C();
-        scalarField& fld = cellDistance.internalField();
+        scalarField& fld = cellDistance.internalFieldRef();
 
         List<pointIndexHit> nearest;
         surfPtr_().findNearest
@@ -257,7 +257,7 @@ void Foam::distanceSurface::createGeometry()
             pointMesh::New(fvm),
             dimensionedScalar("zero", dimLength, 0)
         );
-        pDist.internalField() = pointDistance_;
+        pDist.internalFieldRef() = pointDistance_;
 
         Pout<< "Writing point distance:" << pDist.objectPath() << endl;
         pDist.write();
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
index c89c04fa2bb6da66ae1129f1b5ae84eee0497a81..ec8f04ce135a32e76cff717f309800bae5258dc3 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -125,7 +125,7 @@ void Foam::sampledCuttingPlane::createGeometry()
     // Internal field
     {
         const pointField& cc = fvm.cellCentres();
-        scalarField& fld = cellDistance.internalField();
+        scalarField& fld = cellDistance.internalFieldRef();
 
         forAll(cc, i)
         {
@@ -217,7 +217,7 @@ void Foam::sampledCuttingPlane::createGeometry()
             pointMesh::New(fvm),
             dimensionedScalar("zero", dimLength, 0)
         );
-        pDist.internalField() = pointDistance_;
+        pDist.internalFieldRef() = pointDistance_;
 
         Pout<< "Writing point distance:" << pDist.objectPath() << endl;
         pDist.write();
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
index 4bc236df58b29c23e5985cb644a798b19ccbe7b2..3946c76e4ed2de05d0da9b57306fe79f18f671f8 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
@@ -118,7 +118,7 @@ Foam::sixDoFRigidBodyMotionSolver::sixDoFRigidBodyMotionSolver
         pointPatchDist pDist(pMesh, patchSet_, points0());
 
         // Scaling: 1 up to di then linear down to 0 at do away from patches
-        scale_.internalField() =
+        scale_.internalFieldRef() =
             min
             (
                 max
@@ -130,7 +130,7 @@ Foam::sixDoFRigidBodyMotionSolver::sixDoFRigidBodyMotionSolver
             );
 
         // Convert the scale function to a cosine
-        scale_.internalField() =
+        scale_.internalFieldRef() =
             min
             (
                 max
@@ -241,7 +241,7 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
     }
 
     // Update the displacements
-    pointDisplacement_.internalField() =
+    pointDisplacement_.internalFieldRef() =
         motion_.transform(points0(), scale_) - points0();
 
     // Displacement has changed. Update boundary conditions
diff --git a/src/thermophysicalModels/basic/heThermo/heThermo.C b/src/thermophysicalModels/basic/heThermo/heThermo.C
index 969f6ca1fa418951fe59dd78ba21a82bfdbe9380..f613fb77a67d6dadaf153b8d07a9b45ecc765086 100644
--- a/src/thermophysicalModels/basic/heThermo/heThermo.C
+++ b/src/thermophysicalModels/basic/heThermo/heThermo.C
@@ -54,7 +54,7 @@ heBoundaryCorrection(volScalarField& h)
 template<class BasicThermo, class MixtureType>
 void Foam::heThermo<BasicThermo, MixtureType>::init()
 {
-    scalarField& heCells = he_.internalField();
+    scalarField& heCells = he_.internalFieldRef();
     const scalarField& pCells = this->p_;
     const scalarField& TCells = this->T_;
 
@@ -187,7 +187,7 @@ Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::he
     );
 
     volScalarField& he = the.ref();
-    scalarField& heCells = he.internalField();
+    scalarField& heCells = he.internalFieldRef();
     const scalarField& pCells = p;
     const scalarField& TCells = T;
 
@@ -282,7 +282,7 @@ Foam::heThermo<BasicThermo, MixtureType>::hc() const
     );
 
     volScalarField& hcf = thc.ref();
-    scalarField& hcCells = hcf.internalField();
+    scalarField& hcCells = hcf.internalFieldRef();
 
     forAll(hcCells, celli)
     {
diff --git a/src/thermophysicalModels/basic/psiThermo/hePsiThermo.C b/src/thermophysicalModels/basic/psiThermo/hePsiThermo.C
index 99f228a1596840910cf7d8aa2e0ebc46852a9403..481fe17ec7b52f03ca7ca2847eea30e758a36b04 100644
--- a/src/thermophysicalModels/basic/psiThermo/hePsiThermo.C
+++ b/src/thermophysicalModels/basic/psiThermo/hePsiThermo.C
@@ -33,10 +33,10 @@ void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()
     const scalarField& hCells = this->he_;
     const scalarField& pCells = this->p_;
 
-    scalarField& TCells = this->T_.internalField();
-    scalarField& psiCells = this->psi_.internalField();
-    scalarField& muCells = this->mu_.internalField();
-    scalarField& alphaCells = this->alpha_.internalField();
+    scalarField& TCells = this->T_.internalFieldRef();
+    scalarField& psiCells = this->psi_.internalFieldRef();
+    scalarField& muCells = this->mu_.internalFieldRef();
+    scalarField& alphaCells = this->alpha_.internalFieldRef();
 
     forAll(TCells, celli)
     {
diff --git a/src/thermophysicalModels/basic/rhoThermo/heRhoThermo.C b/src/thermophysicalModels/basic/rhoThermo/heRhoThermo.C
index febff08db1cbdbaefe1f1fee1561a51921d2a570..a6d26271149e7e556e2f122234b11e3399ff13ad 100644
--- a/src/thermophysicalModels/basic/rhoThermo/heRhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/heRhoThermo.C
@@ -33,11 +33,11 @@ void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::calculate()
     const scalarField& hCells = this->he();
     const scalarField& pCells = this->p_;
 
-    scalarField& TCells = this->T_.internalField();
-    scalarField& psiCells = this->psi_.internalField();
-    scalarField& rhoCells = this->rho_.internalField();
-    scalarField& muCells = this->mu_.internalField();
-    scalarField& alphaCells = this->alpha_.internalField();
+    scalarField& TCells = this->T_.internalFieldRef();
+    scalarField& psiCells = this->psi_.internalFieldRef();
+    scalarField& rhoCells = this->rho_.internalFieldRef();
+    scalarField& muCells = this->mu_.internalFieldRef();
+    scalarField& alphaCells = this->alpha_.internalFieldRef();
 
     forAll(TCells, celli)
     {
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
index a2776a561cdb93de9b4f00833e6912f3f5eb8463..045dde224fc251acda9124f12258100ccf56de9b 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
@@ -213,7 +213,7 @@ Foam::radiation::greyMeanAbsorptionEmission::aCont(const label bandI) const
         )
     );
 
-    scalarField& a = ta.ref().internalField();
+    scalarField& a = ta.ref().internalFieldRef();
 
     forAll(a, celli)
     {
@@ -300,11 +300,11 @@ Foam::radiation::greyMeanAbsorptionEmission::ECont(const label bandI) const
 
         if (dQ.dimensions() == dimEnergy/dimTime)
         {
-            E.ref().internalField() = EhrrCoeff_*dQ/mesh_.V();
+            E.ref().internalFieldRef() = EhrrCoeff_*dQ/mesh_.V();
         }
         else if (dQ.dimensions() == dimEnergy/dimTime/dimVolume)
         {
-            E.ref().internalField() = EhrrCoeff_*dQ;
+            E.ref().internalFieldRef() = EhrrCoeff_*dQ;
         }
         else
         {
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C
index 0e2238eeecf3a8f5328bfb1c153dd73bade7c861..bca25380bc78b187ec24947d282b311ab4b3a515 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C
@@ -162,7 +162,7 @@ calc(const label propertyId) const
         )
     );
 
-    scalarField& a = ta.ref().internalField();
+    scalarField& a = ta.ref().internalFieldRef();
 
     forAllConstIter(HashTable<label>, speciesNames_, iter)
     {
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
index 673043c39bebba79e6343a789259a08af98b6929..cd1a916f45c2a09e7d44f2031ffa88afce454d21 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
@@ -177,7 +177,7 @@ Foam::radiation::wideBandAbsorptionEmission::aCont(const label bandI) const
         )
     );
 
-    scalarField& a = ta.ref().internalField();
+    scalarField& a = ta.ref().internalFieldRef();
 
     forAll(a, i)
     {
@@ -255,7 +255,7 @@ Foam::radiation::wideBandAbsorptionEmission::ECont(const label bandI) const
 
         if (dQ.dimensions() == dimEnergy/dimTime)
         {
-            E.ref().internalField() =
+            E.ref().internalFieldRef() =
                 iEhrrCoeffs_[bandI]
                *dQ.internalField()
                *(iBands_[bandI][1] - iBands_[bandI][0])
@@ -264,7 +264,7 @@ Foam::radiation::wideBandAbsorptionEmission::ECont(const label bandI) const
         }
         else if (dQ.dimensions() == dimEnergy/dimTime/dimVolume)
         {
-            E.ref().internalField() =
+            E.ref().internalFieldRef() =
                 iEhrrCoeffs_[bandI]
                *dQ.internalField()
                *(iBands_[bandI][1] - iBands_[bandI][0])
@@ -291,9 +291,9 @@ void Foam::radiation::wideBandAbsorptionEmission::correct
 
     for (label j=0; j<nBands_; j++)
     {
-        aLambda[j].internalField() = this->a(j);
+        aLambda[j].internalFieldRef() = this->a(j);
 
-        a.internalField() +=
+        a.internalFieldRef() +=
             aLambda[j].internalField()
            *(iBands_[j][1] - iBands_[j][0])
            /totalWaveLength_;
diff --git a/src/thermophysicalModels/reactionThermo/psiuReactionThermo/heheuPsiThermo.C b/src/thermophysicalModels/reactionThermo/psiuReactionThermo/heheuPsiThermo.C
index 703e088005c882a071b4f62c7fd1c15f3921b026..b36603cedbde90b94ff6b9621ef231bcd3898170 100644
--- a/src/thermophysicalModels/reactionThermo/psiuReactionThermo/heheuPsiThermo.C
+++ b/src/thermophysicalModels/reactionThermo/psiuReactionThermo/heheuPsiThermo.C
@@ -36,11 +36,11 @@ void Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::calculate()
     const scalarField& heuCells = this->heu_;
     const scalarField& pCells = this->p_;
 
-    scalarField& TCells = this->T_.internalField();
-    scalarField& TuCells = this->Tu_.internalField();
-    scalarField& psiCells = this->psi_.internalField();
-    scalarField& muCells = this->mu_.internalField();
-    scalarField& alphaCells = this->alpha_.internalField();
+    scalarField& TCells = this->T_.internalFieldRef();
+    scalarField& TuCells = this->Tu_.internalFieldRef();
+    scalarField& psiCells = this->psi_.internalFieldRef();
+    scalarField& muCells = this->mu_.internalFieldRef();
+    scalarField& alphaCells = this->alpha_.internalFieldRef();
 
     forAll(TCells, celli)
     {
@@ -176,7 +176,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::heheuPsiThermo
         this->heuBoundaryTypes()
     )
 {
-    scalarField& heuCells = this->heu_.internalField();
+    scalarField& heuCells = this->heu_.internalFieldRef();
     const scalarField& pCells = this->p_;
     const scalarField& TuCells = this->Tu_;
 
@@ -308,7 +308,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::Tb() const
     );
 
     volScalarField& Tb_ = tTb.ref();
-    scalarField& TbCells = Tb_.internalField();
+    scalarField& TbCells = Tb_.internalFieldRef();
     const scalarField& pCells = this->p_;
     const scalarField& TCells = this->T_;
     const scalarField& hCells = this->he_;
@@ -368,7 +368,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::psiu() const
     );
 
     volScalarField& psiu = tpsiu.ref();
-    scalarField& psiuCells = psiu.internalField();
+    scalarField& psiuCells = psiu.internalFieldRef();
     const scalarField& TuCells = this->Tu_;
     const scalarField& pCells = this->p_;
 
@@ -422,7 +422,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::psib() const
     );
 
     volScalarField& psib = tpsib.ref();
-    scalarField& psibCells = psib.internalField();
+    scalarField& psibCells = psib.internalFieldRef();
     const volScalarField Tb_(Tb());
     const scalarField& TbCells = Tb_;
     const scalarField& pCells = this->p_;
@@ -477,7 +477,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::muu() const
     );
 
     volScalarField& muu_ = tmuu.ref();
-    scalarField& muuCells = muu_.internalField();
+    scalarField& muuCells = muu_.internalFieldRef();
     const scalarField& pCells = this->p_;
     const scalarField& TuCells = this->Tu_;
 
@@ -535,7 +535,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::mub() const
     );
 
     volScalarField& mub_ = tmub.ref();
-    scalarField& mubCells = mub_.internalField();
+    scalarField& mubCells = mub_.internalFieldRef();
     const volScalarField Tb_(Tb());
     const scalarField& pCells = this->p_;
     const scalarField& TbCells = Tb_;
diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
index 383b3a3b0a0788f4275fc1b36081145447001c49..0f70a25a166cfacf44911c5d8047e3beb92c7849 100644
--- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
@@ -110,7 +110,7 @@ pyrolysisChemistryModel
             );
 
             // Calculate inital values of Ysi0 = rho*delta*Yi
-            Ys0_[fieldI].internalField() =
+            Ys0_[fieldI].internalFieldRef() =
                 this->solidThermo().rho()
                *max(this->Ys_[fieldI], scalar(0.001))*mesh.V();
         }
diff --git a/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C b/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C
index b6b4bea3be2339dadccefb091612822e736e8166..e0812d4a78427ba7e8d762087755d6331217ae94 100644
--- a/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C
+++ b/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C
@@ -31,12 +31,12 @@ License
 template<class BasicSolidThermo, class MixtureType>
 void Foam::heSolidThermo<BasicSolidThermo, MixtureType>::calculate()
 {
-    scalarField& TCells = this->T_.internalField();
+    scalarField& TCells = this->T_.internalFieldRef();
 
     const scalarField& hCells = this->he_;
     const scalarField& pCells = this->p_;
-    scalarField& rhoCells = this->rho_.internalField();
-    scalarField& alphaCells = this->alpha_.internalField();
+    scalarField& rhoCells = this->rho_.internalFieldRef();
+    scalarField& alphaCells = this->alpha_.internalFieldRef();
 
     forAll(TCells, celli)
     {
@@ -218,7 +218,7 @@ Foam::heSolidThermo<BasicSolidThermo, MixtureType>::Kappa() const
     );
 
     volVectorField& Kappa = tKappa.ref();
-    vectorField& KappaCells = Kappa.internalField();
+    vectorField& KappaCells = Kappa.internalFieldRef();
     const scalarField& TCells = this->T_;
     const scalarField& pCells = this->p_;