From 8667536a2bf56f22465b0b2e7439710ec78155dd Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 2 Dec 2011 10:38:18 +0000
Subject: [PATCH] ENH: various: move boundaryField(), internalField() out of
 loop

---
 .../backgroundMeshDecomposition.C             | 35 +++++++------
 .../GeometricField/MapGeometricFields.H       |  9 ++--
 .../consumptionSpeed/consumptionSpeed.C       | 14 ++++--
 .../dynamicRefineFvMesh/dynamicRefineFvMesh.C | 17 ++++---
 .../fvMeshAdder/fvMeshAdderTemplates.C        | 49 ++++++++++---------
 .../cfdTools/general/adjustPhi/adjustPhi.C    | 10 ++--
 .../CoEulerDdtScheme/CoEulerDdtScheme.C       |  5 +-
 .../finiteVolume/fvc/fvcAverage.C             |  7 ++-
 .../extendedLeastSquaresVectors.C             | 10 ++--
 .../leastSquaresGrad/leastSquaresVectors.C    |  8 +--
 .../fvMatrices/solvers/MULES/MULESTemplates.C |  7 ++-
 .../molecule/mdTools/averageMDFields.H        | 29 ++++++-----
 .../old/mdTools/averageMDFields.H             | 30 +++++++-----
 .../solidMixtureThermo/solidMixtureThermo.C   | 48 ++++++++++++------
 .../ODESolidChemistryModel.C                  |  4 +-
 15 files changed, 169 insertions(+), 113 deletions(-)

diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index ad420d096f5..edafdb50f91 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -146,6 +146,9 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
 
     decompositionMethod& decomposer = decomposerPtr_();
 
+    volScalarField::InternalField& icellWeights = cellWeights.internalField();
+
+
     // For each cell in the mesh has it been determined if it is fully
     // inside, outside, or overlaps the surface
     labelList volumeStatus
@@ -214,10 +217,10 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                         volumeStatus[cellI] = searchableSurface::UNKNOWN;
                     }
 
-                    cellWeights.internalField()[cellI] = max
+                    icellWeights[cellI] = max
                     (
                         1.0,
-                        cellWeights.internalField()[cellI]/8.0
+                        icellWeights[cellI]/8.0
                     );
                 }
 
@@ -405,7 +408,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
             (
                 mesh_,
                 mesh_.cellCentres(),
-                cellWeights
+                icellWeights
             );
 
             fvMeshDistribute distributor(mesh_, mergeDist_);
@@ -629,6 +632,8 @@ Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
     volScalarField& cellWeights
 ) const
 {
+    volScalarField::InternalField& icellWeights = cellWeights.internalField();
+
     labelHashSet cellsToRefine;
 
     // Determine/update the status of each cell
@@ -650,7 +655,7 @@ Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
                 (
                     cellI,
                     volumeStatus[cellI],
-                    cellWeights.internalField()[cellI]
+                    icellWeights[cellI]
                 )
             )
             {
@@ -871,15 +876,17 @@ Foam::backgroundMeshDecomposition::distribute
         mesh_.write();
     }
 
+    volScalarField::InternalField& icellWeights = cellWeights.internalField();
+
     while (true)
     {
         // Refine large cells if necessary
 
         label nOccupiedCells = 0;
 
-        forAll(cellWeights.internalField(), cI)
+        forAll(icellWeights, cI)
         {
-            if (cellWeights.internalField()[cI] > 1 - SMALL)
+            if (icellWeights[cI] > 1 - SMALL)
             {
                 nOccupiedCells++;
             }
@@ -910,9 +917,9 @@ Foam::backgroundMeshDecomposition::distribute
 
         labelHashSet cellsToRefine;
 
-        forAll(cellWeights, cWI)
+        forAll(icellWeights, cWI)
         {
-            if (cellWeights.internalField()[cWI] > cellWeightLimit)
+            if (icellWeights[cWI] > cellWeightLimit)
             {
                 cellsToRefine.insert(cWI);
             }
@@ -943,7 +950,7 @@ Foam::backgroundMeshDecomposition::distribute
         {
             label cellI = newCellsToRefine[nCTRI];
 
-            cellWeights.internalField()[cellI] /= 8.0;
+            icellWeights[cellI] /= 8.0;
         }
 
         // Mesh changing engine.
@@ -1075,9 +1082,9 @@ Foam::backgroundMeshDecomposition::distribute
         printMeshData(mesh_);
 
         Pout<< "    Pre distribute sum(cellWeights) "
-            << sum(cellWeights.internalField())
+            << sum(icellWeights)
             << " max(cellWeights) "
-            << max(cellWeights.internalField())
+            << max(icellWeights)
             << endl;
     }
 
@@ -1085,7 +1092,7 @@ Foam::backgroundMeshDecomposition::distribute
     (
         mesh_,
         mesh_.cellCentres(),
-        cellWeights
+        icellWeights
     );
 
     Info<< "    Redistributing background mesh cells" << endl;
@@ -1101,9 +1108,9 @@ Foam::backgroundMeshDecomposition::distribute
         printMeshData(mesh_);
 
         Pout<< "    Post distribute sum(cellWeights) "
-            << sum(cellWeights.internalField())
+            << sum(icellWeights)
             << " max(cellWeights) "
-            << max(cellWeights.internalField())
+            << max(icellWeights)
             << endl;
 
         // const_cast<Time&>(mesh_.time())++;
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H
index c32428287a5..10150f7a64b 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H
@@ -131,17 +131,16 @@ void MapGeometricFields
             );
 
             // Map the patch fields
-            forAll(field.boundaryField(), patchi)
+            typename GeometricField<Type, PatchField, GeoMesh>
+            ::GeometricBoundaryField& bfield = field.boundaryField();
+            forAll(bfield, patchi)
             {
                 // Cannot check sizes for patch fields because of
                 // empty fields in FV and because point fields get their size
                 // from the patch which has already been resized
                 //
 
-                field.boundaryField()[patchi].autoMap
-                (
-                    mapper.boundaryMap()[patchi]
-                );
+                bfield[patchi].autoMap(mapper.boundaryMap()[patchi]);
             }
 
             field.instance() = field.time().timeName();
diff --git a/src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C b/src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C
index 44f1e79ca8c..782f8456c4a 100644
--- a/src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C
+++ b/src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C
@@ -104,16 +104,20 @@ Foam::tmp<Foam::volScalarField> Foam::consumptionSpeed::omega0Sigma
 
     volScalarField& omega0 = tomega0();
 
-    forAll(omega0, celli)
+    volScalarField::InternalField& iomega0 = omega0.internalField();
+
+    forAll(iomega0, celli)
     {
-        omega0[celli] = omega0Sigma(sigma[celli], 1.0);
+        iomega0[celli] = omega0Sigma(sigma[celli], 1.0);
     }
 
-    forAll(omega0.boundaryField(), patchi)
+    volScalarField::GeometricBoundaryField& bomega0 = omega0.boundaryField();
+
+    forAll(bomega0, patchi)
     {
-        forAll(omega0.boundaryField()[patchi], facei)
+        forAll(bomega0[patchi], facei)
         {
-            omega0.boundaryField()[patchi][facei] =
+            bomega0[patchi][facei] =
                 omega0Sigma
                 (
                     sigma.boundaryField()[patchi][facei],
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
index 0af332a336d..3b295f46dd8 100644
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
@@ -360,9 +360,11 @@ Foam::dynamicRefineFvMesh::refine
             }
 
             // Recalculate new boundary faces.
-            forAll(phi.boundaryField(), patchI)
+            surfaceScalarField::GeometricBoundaryField& bphi =
+                phi.boundaryField();
+            forAll(bphi, patchI)
             {
-                fvsPatchScalarField& patchPhi = phi.boundaryField()[patchI];
+                fvsPatchScalarField& patchPhi = bphi[patchI];
                 const fvsPatchScalarField& patchPhiU =
                     phiU.boundaryField()[patchI];
 
@@ -404,8 +406,7 @@ Foam::dynamicRefineFvMesh::refine
                     const fvsPatchScalarField& patchPhiU =
                         phiU.boundaryField()[patchI];
 
-                    fvsPatchScalarField& patchPhi =
-                        phi.boundaryField()[patchI];
+                    fvsPatchScalarField& patchPhi = bphi[patchI];
 
                     patchPhi[i] = patchPhiU[i];
                 }
@@ -549,6 +550,9 @@ Foam::dynamicRefineFvMesh::unrefine
             }
 
             surfaceScalarField& phi = const_cast<surfaceScalarField&>(*iter());
+            surfaceScalarField::GeometricBoundaryField& bphi =
+                phi.boundaryField();
+
             const surfaceScalarField phiU
             (
                 fvc::interpolate
@@ -582,10 +586,7 @@ Foam::dynamicRefineFvMesh::unrefine
 
                             const fvsPatchScalarField& patchPhiU =
                                 phiU.boundaryField()[patchI];
-
-                            fvsPatchScalarField& patchPhi =
-                                phi.boundaryField()[patchI];
-
+                            fvsPatchScalarField& patchPhi = bphi[patchI];
                             patchPhi[i] = patchPhiU[i];
                         }
                     }
diff --git a/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C b/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
index f841e2b74dc..bb01e2fa833 100644
--- a/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
+++ b/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
@@ -81,6 +81,9 @@ void Foam::fvMeshAdder::MapVolField
     // Patch fields from old mesh
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+    typename GeometricField<Type, fvPatchField, volMesh>::
+    GeometricBoundaryField& bfld = fld.boundaryField();
+
     {
         const labelList& oldPatchMap = meshMap.oldPatchMap();
         const labelList& oldPatchStarts = meshMap.oldPatchStarts();
@@ -121,18 +124,18 @@ void Foam::fvMeshAdder::MapVolField
 
 
         // Sort deleted ones last so is now in newPatch ordering
-        fld.boundaryField().reorder(oldToNew);
+        bfld.reorder(oldToNew);
         // Extend to covers all patches
-        fld.boundaryField().setSize(mesh.boundaryMesh().size());
+        bfld.setSize(mesh.boundaryMesh().size());
         // Delete unused patches
         for
         (
             label newPatchI = nUsedPatches;
-            newPatchI < fld.boundaryField().size();
+            newPatchI < bfld.size();
             newPatchI++
         )
         {
-            fld.boundaryField().set(newPatchI, NULL);
+            bfld.set(newPatchI, NULL);
         }
 
 
@@ -167,12 +170,12 @@ void Foam::fvMeshAdder::MapVolField
                 //   value
                 // - hope that field mapping allows aliasing since old and new
                 //   are same memory!
-                fld.boundaryField().set
+                bfld.set
                 (
                     newPatchI,
                     fvPatchField<Type>::New
                     (
-                        fld.boundaryField()[newPatchI], // old field
+                        bfld[newPatchI],                // old field
                         mesh.boundary()[newPatchI],     // new fvPatch
                         fld.dimensionedInternalField(), // new internal field
                         patchMapper                     // mapper (new to old)
@@ -201,7 +204,7 @@ void Foam::fvMeshAdder::MapVolField
                 const polyPatch& oldPatch =
                     fldToAdd.mesh().boundaryMesh()[patchI];
 
-                if (!fld.boundaryField()(newPatchI))
+                if (!bfld(newPatchI))
                 {
                     // First occurrence of newPatchI. Map from existing
                     // patchField
@@ -221,7 +224,7 @@ void Foam::fvMeshAdder::MapVolField
 
                     directFvPatchFieldMapper patchMapper(newToAdded);
 
-                    fld.boundaryField().set
+                    bfld.set
                     (
                         newPatchI,
                         fvPatchField<Type>::New
@@ -255,7 +258,7 @@ void Foam::fvMeshAdder::MapVolField
                     const fvPatchField<Type>& addedFld =
                         fldToAdd.boundaryField()[patchI];
 
-                    fvPatchField<Type>& newFld = fld.boundaryField()[newPatchI];
+                    fvPatchField<Type>& newFld = bfld[newPatchI];
 
                     forAll(newFld, i)
                     {
@@ -357,6 +360,9 @@ void Foam::fvMeshAdder::MapSurfaceField
     const fvMesh& mesh = fld.mesh();
     const labelList& oldPatchStarts = meshMap.oldPatchStarts();
 
+    typename GeometricField<Type, fvsPatchField, surfaceMesh>::
+    GeometricBoundaryField& bfld = fld.boundaryField();
+
     // Internal field
     // ~~~~~~~~~~~~~~
 
@@ -375,9 +381,9 @@ void Foam::fvMeshAdder::MapSurfaceField
         // Faces that were boundary faces but are not anymore.
         // Use owner value (so lowest numbered cell, i.e. from 'old' not 'added'
         // mesh)
-        forAll(fld.boundaryField(), patchI)
+        forAll(bfld, patchI)
         {
-            const fvsPatchField<Type>& pf = fld.boundaryField()[patchI];
+            const fvsPatchField<Type>& pf = bfld[patchI];
 
             label start = oldPatchStarts[patchI];
 
@@ -436,18 +442,18 @@ void Foam::fvMeshAdder::MapSurfaceField
 
 
         // Sort deleted ones last so is now in newPatch ordering
-        fld.boundaryField().reorder(oldToNew);
+        bfld.reorder(oldToNew);
         // Extend to covers all patches
-        fld.boundaryField().setSize(mesh.boundaryMesh().size());
+        bfld.setSize(mesh.boundaryMesh().size());
         // Delete unused patches
         for
         (
             label newPatchI = nUsedPatches;
-            newPatchI < fld.boundaryField().size();
+            newPatchI < bfld.size();
             newPatchI++
         )
         {
-            fld.boundaryField().set(newPatchI, NULL);
+            bfld.set(newPatchI, NULL);
         }
 
 
@@ -478,16 +484,16 @@ void Foam::fvMeshAdder::MapSurfaceField
                 // Create new patchField with same type as existing one.
                 // Note:
                 // - boundaryField already in new order so access with newPatchI
-                // - fld.boundaryField()[newPatchI] both used for type and old
+                // - bfld[newPatchI] both used for type and old
                 //   value
                 // - hope that field mapping allows aliasing since old and new
                 //   are same memory!
-                fld.boundaryField().set
+                bfld.set
                 (
                     newPatchI,
                     fvsPatchField<Type>::New
                     (
-                        fld.boundaryField()[newPatchI], // old field
+                        bfld[newPatchI],                // old field
                         mesh.boundary()[newPatchI],     // new fvPatch
                         fld.dimensionedInternalField(), // new internal field
                         patchMapper                     // mapper (new to old)
@@ -516,7 +522,7 @@ void Foam::fvMeshAdder::MapSurfaceField
                 const polyPatch& oldPatch =
                     fldToAdd.mesh().boundaryMesh()[patchI];
 
-                if (!fld.boundaryField()(newPatchI))
+                if (!bfld(newPatchI))
                 {
                     // First occurrence of newPatchI. Map from existing
                     // patchField
@@ -536,7 +542,7 @@ void Foam::fvMeshAdder::MapSurfaceField
 
                     directFvPatchFieldMapper patchMapper(newToAdded);
 
-                    fld.boundaryField().set
+                    bfld.set
                     (
                         newPatchI,
                         fvsPatchField<Type>::New
@@ -570,8 +576,7 @@ void Foam::fvMeshAdder::MapSurfaceField
                     const fvsPatchField<Type>& addedFld =
                         fldToAdd.boundaryField()[patchI];
 
-                    fvsPatchField<Type>& newFld =
-                        fld.boundaryField()[newPatchI];
+                    fvsPatchField<Type>& newFld = bfld[newPatchI];
 
                     forAll(newFld, i)
                     {
diff --git a/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C b/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
index 7c72beaf99d..206e7c4e95e 100644
--- a/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
+++ b/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
@@ -47,10 +47,12 @@ bool Foam::adjustPhi
         scalar fixedMassOut = 0.0;
         scalar adjustableMassOut = 0.0;
 
-        forAll(phi.boundaryField(), patchi)
+        surfaceScalarField::GeometricBoundaryField& bphi = phi.boundaryField();
+
+        forAll(bphi, patchi)
         {
             const fvPatchVectorField& Up = U.boundaryField()[patchi];
-            const fvsPatchScalarField& phip = phi.boundaryField()[patchi];
+            const fvsPatchScalarField& phip = bphi[patchi];
 
             if (!isA<processorFvsPatchScalarField>(phip))
             {
@@ -123,10 +125,10 @@ bool Foam::adjustPhi
                 << exit(FatalError);
         }
 
-        forAll(phi.boundaryField(), patchi)
+        forAll(bphi, patchi)
         {
             const fvPatchVectorField& Up = U.boundaryField()[patchi];
-            fvsPatchScalarField& phip = phi.boundaryField()[patchi];
+            fvsPatchScalarField& phip = bphi[patchi];
 
             if (!isA<processorFvsPatchScalarField>(phip))
             {
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
index e92de99aec8..d6ef3989cdd 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
@@ -75,7 +75,10 @@ tmp<volScalarField> CoEulerDdtScheme<Type>::CorDeltaT() const
             max(corDeltaT[neighbour[faceI]], cofrDeltaT[faceI]);
     }
 
-    forAll(corDeltaT.boundaryField(), patchi)
+    volScalarField::GeometricBoundaryField& bcorDeltaT =
+        corDeltaT.boundaryField();
+
+    forAll(bcorDeltaT, patchi)
     {
         const fvsPatchScalarField& pcofrDeltaT =
             cofrDeltaT.boundaryField()[patchi];
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcAverage.C b/src/finiteVolume/finiteVolume/fvc/fvcAverage.C
index 91c67d823a7..2756f56de56 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcAverage.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcAverage.C
@@ -73,9 +73,12 @@ average
         surfaceSum(mesh.magSf()*ssf)/surfaceSum(mesh.magSf())
     )().internalField();
 
-    forAll(av.boundaryField(), patchi)
+    typename GeometricField<Type, fvPatchField, volMesh>::
+    GeometricBoundaryField& bav = av.boundaryField();
+
+    forAll(bav, patchi)
     {
-        av.boundaryField()[patchi] = ssf.boundaryField()[patchi];
+        bav[patchi] = ssf.boundaryField()[patchi];
     }
 
     av.correctBoundaryConditions();
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C
index 4bf1a4dd6b5..ebec7bb7d5a 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C
@@ -130,9 +130,11 @@ void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
 
     // Visit the boundaries. Coupled boundaries are taken into account
     // in the construction of d vectors.
-    forAll(lsP.boundaryField(), patchi)
+    surfaceVectorField::GeometricBoundaryField& blsP = lsP.boundaryField();
+
+    forAll(blsP, patchi)
     {
-        const fvPatch& p = lsP.boundaryField()[patchi].patch();
+        const fvPatch& p = blsP[patchi].patch();
         const labelUList& faceCells = p.faceCells();
 
         // Build the d-vectors
@@ -242,9 +244,9 @@ void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
             ((-1.0)/magSqr(d[facei]))*(invDd[neighbour[facei]] & d);
     }
 
-    forAll(lsP.boundaryField(), patchI)
+    forAll(blsP, patchI)
     {
-        fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchI];
+        fvsPatchVectorField& patchLsP = blsP[patchI];
 
         const fvPatch& p = patchLsP.patch();
         const labelUList& faceCells = p.faceCells();
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
index 27f5b72fcfd..86b4bc0ad31 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
@@ -121,7 +121,9 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
     }
 
 
-    forAll(lsP.boundaryField(), patchi)
+    surfaceVectorField::GeometricBoundaryField& blsP = lsP.boundaryField();
+
+    forAll(blsP, patchi)
     {
         const fvsPatchScalarField& pw = w.boundaryField()[patchi];
         const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
@@ -172,9 +174,9 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
         lsN[facei] = -w[facei]*magSfByMagSqrd*(invDd[nei] & d);
     }
 
-    forAll(lsP.boundaryField(), patchi)
+    forAll(blsP, patchi)
     {
-        fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchi];
+        fvsPatchVectorField& patchLsP = blsP[patchi];
 
         const fvsPatchScalarField& pw = w.boundaryField()[patchi];
         const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index 55478e59aa4..98ed5368608 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -668,9 +668,12 @@ void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
         limitSum(phiPsiCorrsInternal);
     }
 
-    forAll(phiPsiCorrs[0].boundaryField(), patchi)
+    surfaceScalarField::GeometricBoundaryField& bfld =
+        phiPsiCorrs[0].boundaryField();
+
+    forAll(bfld, patchi)
     {
-        if (phiPsiCorrs[0].boundaryField()[patchi].coupled())
+        if (bfld[patchi].coupled())
         {
             UPtrList<scalarField> phiPsiCorrsPatch(phiPsiCorrs.size());
             forAll(phiPsiCorrs, phasei)
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
index 215d6bbef8f..a19febe07b6 100644
--- a/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
@@ -76,18 +76,18 @@ if (runTime.outputTime())
         }
     }
 
-    forAll(totalVelocity.internalField(), tV)
+    volVectorField::InternalField& itotalVelocity =
+        totalVelocity.internalField();
+
+    forAll(itotalVelocity, tV)
     {
         if (totalMass_sum[tV] > VSMALL)
         {
-            totalVelocity.internalField()[tV] =
-                totalMomentum_sum[tV]
-                /totalMass_sum[tV];
+            itotalVelocity[tV] = totalMomentum_sum[tV]/totalMass_sum[tV];
         }
         else
         {
-            totalVelocity.internalField()[tV] =
-            vector::zero;
+            itotalVelocity[tV] = vector::zero;
         }
     }
 
@@ -150,17 +150,20 @@ if (runTime.outputTime())
         }
     }
 
-    forAll(totalTemperature.internalField(), tT)
+    volScalarField::InternalField& itotalTemperature =
+        totalTemperature.internalField();
+
+    forAll(itotalTemperature, tT)
     {
         if (totalN_sum[tT] > 0)
         {
-            totalTemperature.internalField()[tT] =
+            itotalTemperature[tT] =
                 totalTemperatureVTerms_sum[tT]
                 /(3.0 * moleculeCloud::kb * totalN_sum[tT]);
         }
         else
         {
-            totalTemperature.internalField()[tT] = 0.0;
+            itotalTemperature[tT] = 0.0;
         }
     }
 
@@ -206,17 +209,19 @@ if (runTime.outputTime())
         }
     }
 
-    forAll(totalMeanKE.internalField(), tMKE)
+    volScalarField::InternalField& itotalMeanKE = totalMeanKE.internalField();
+
+    forAll(itotalMeanKE, tMKE)
     {
         if (totalN_sum[tMKE] > 0)
         {
-            totalMeanKE.internalField()[tMKE] =
+            itotalMeanKE[tMKE] =
                 totalKE_sum[tMKE]
                 /totalN_sum[tMKE];
         }
         else
         {
-            totalMeanKE.internalField()[tMKE] = 0.0;
+            itotalMeanKE[tMKE] = 0.0;
         }
     }
 }
diff --git a/src/lagrangian/molecularDynamics/old/mdTools/averageMDFields.H b/src/lagrangian/molecularDynamics/old/mdTools/averageMDFields.H
index 5a195aa9389..a19febe07b6 100644
--- a/src/lagrangian/molecularDynamics/old/mdTools/averageMDFields.H
+++ b/src/lagrangian/molecularDynamics/old/mdTools/averageMDFields.H
@@ -76,18 +76,18 @@ if (runTime.outputTime())
         }
     }
 
-    forAll(totalVelocity.internalField(), tV)
+    volVectorField::InternalField& itotalVelocity =
+        totalVelocity.internalField();
+
+    forAll(itotalVelocity, tV)
     {
         if (totalMass_sum[tV] > VSMALL)
         {
-            totalVelocity.internalField()[tV] =
-                totalMomentum_sum[tV]
-                /totalMass_sum[tV];
+            itotalVelocity[tV] = totalMomentum_sum[tV]/totalMass_sum[tV];
         }
         else
         {
-            totalVelocity.internalField()[tV] =
-            vector::zero;
+            itotalVelocity[tV] = vector::zero;
         }
     }
 
@@ -150,17 +150,20 @@ if (runTime.outputTime())
         }
     }
 
-    forAll(totalTemperature.internalField(), tT)
+    volScalarField::InternalField& itotalTemperature =
+        totalTemperature.internalField();
+
+    forAll(itotalTemperature, tT)
     {
         if (totalN_sum[tT] > 0)
         {
-            totalTemperature.internalField()[tT] =
+            itotalTemperature[tT] =
                 totalTemperatureVTerms_sum[tT]
                 /(3.0 * moleculeCloud::kb * totalN_sum[tT]);
         }
         else
         {
-            totalTemperature.internalField()[tT] = 0.0;
+            itotalTemperature[tT] = 0.0;
         }
     }
 
@@ -206,18 +209,19 @@ if (runTime.outputTime())
         }
     }
 
-    forAll(totalMeanKE.internalField(), tMKE)
+    volScalarField::InternalField& itotalMeanKE = totalMeanKE.internalField();
+
+    forAll(itotalMeanKE, tMKE)
     {
         if (totalN_sum[tMKE] > 0)
         {
-            totalMeanKE.internalField()[tMKE] =
+            itotalMeanKE[tMKE] =
                 totalKE_sum[tMKE]
                 /totalN_sum[tMKE];
         }
         else
         {
-            totalMeanKE.internalField()[tMKE] = 0.0;
+            itotalMeanKE[tMKE] = 0.0;
         }
     }
 }
-
diff --git a/src/thermophysicalModels/basicSolidThermo/solidMixtureThermo/solidMixtureThermo/solidMixtureThermo.C b/src/thermophysicalModels/basicSolidThermo/solidMixtureThermo/solidMixtureThermo/solidMixtureThermo.C
index 912564714a5..024f6cfcedb 100644
--- a/src/thermophysicalModels/basicSolidThermo/solidMixtureThermo/solidMixtureThermo/solidMixtureThermo.C
+++ b/src/thermophysicalModels/basicSolidThermo/solidMixtureThermo/solidMixtureThermo/solidMixtureThermo.C
@@ -39,16 +39,20 @@ void Foam::solidMixtureThermo<MixtureType>::calculate()
     scalarField& sigmaSCells = sigmaS_.internalField();
     scalarField& emissivityCells = emissivity_.internalField();
 
-    forAll(T_.internalField(), celli)
+    const volScalarField::InternalField& iT = T_.internalField();
+
+    forAll(iT, celli)
     {
-        rhoCells[celli] = MixtureType::rho(T_[celli], celli);
-        kappaCells[celli] = MixtureType::kappa(T_[celli], celli);
-        sigmaSCells[celli] = MixtureType::sigmaS(T_[celli], celli);
-        KCells[celli] = MixtureType::K(T_[celli], celli);
-        emissivityCells[celli] = MixtureType::emissivity(T_[celli], celli);
+        rhoCells[celli] = MixtureType::rho(iT[celli], celli);
+        kappaCells[celli] = MixtureType::kappa(iT[celli], celli);
+        sigmaSCells[celli] = MixtureType::sigmaS(iT[celli], celli);
+        KCells[celli] = MixtureType::K(iT[celli], celli);
+        emissivityCells[celli] = MixtureType::emissivity(iT[celli], celli);
     }
 
-    forAll(T_.boundaryField(), patchI)
+    const volScalarField::GeometricBoundaryField& bT = T_.boundaryField();
+
+    forAll(bT, patchI)
     {
         rho_.boundaryField()[patchI] == this->rho(patchI)();
         K_.boundaryField()[patchI] == this->K(patchI)();
@@ -158,14 +162,18 @@ Foam::solidMixtureThermo<MixtureType>::Cp() const
     );
     volScalarField& Cp = tCp();
 
-    forAll(T_.internalField(), celli)
+    const volScalarField::InternalField& iT = T_.internalField();
+
+    forAll(iT, celli)
     {
         Cp[celli] = MixtureType::Cp(T_[celli], celli);
     }
 
-    forAll(Cp.boundaryField(), patchI)
+    volScalarField::GeometricBoundaryField& bCp = Cp.boundaryField();
+
+    forAll(bCp, patchI)
     {
-        Cp.boundaryField()[patchI] == this->Cp(patchI)();
+        bCp[patchI] == this->Cp(patchI)();
     }
 
     return tCp;
@@ -194,14 +202,18 @@ Foam::solidMixtureThermo<MixtureType>::hs() const
     );
     volScalarField& hs = ths();
 
-    forAll(T_.internalField(), celli)
+    const volScalarField::InternalField& iT = T_.internalField();
+
+    forAll(iT, celli)
     {
         hs[celli] = MixtureType::hs(T_[celli], celli);
     }
 
-    forAll(hs.boundaryField(), patchI)
+    volScalarField::GeometricBoundaryField& bHs = hs.boundaryField();
+
+    forAll(bHs, patchI)
     {
-        hs.boundaryField()[patchI] == this->hs(patchI)();
+        bHs[patchI] == this->hs(patchI)();
     }
 
     return ths;
@@ -230,14 +242,18 @@ Foam::solidMixtureThermo<MixtureType>::Hf() const
     );
     volScalarField& hf = thF();
 
-    forAll(T_.internalField(), celli)
+    const volScalarField::InternalField& iT = T_.internalField();
+
+    forAll(iT, celli)
     {
         hf[celli] = MixtureType::hf(T_[celli], celli);
     }
 
-    forAll(hf.boundaryField(), patchI)
+    volScalarField::GeometricBoundaryField& bhf = hf.boundaryField();
+
+    forAll(bhf, patchI)
     {
-        hf.boundaryField()[patchI] == this->Hf(patchI)();
+        bhf[patchI] == this->Hf(patchI)();
     }
 
     return thF;
diff --git a/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C
index 1c5cc36d09d..6c68214b63f 100644
--- a/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C
@@ -760,11 +760,11 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::gasHs
         )
     );
 
-    volScalarField& gasHs = tHs();
+    volScalarField::InternalField& gasHs = tHs().internalField();
 
     const GasThermo& mixture = gasThermo_[index];
 
-    forAll(gasHs.internalField(), cellI)
+    forAll(gasHs, cellI)
     {
         gasHs[cellI] = mixture.Hs(T[cellI]);
     }
-- 
GitLab