diff --git a/applications/utilities/postProcessing/dataConversion/foamDataToFluent/writeFluentScalarField.C b/applications/utilities/postProcessing/dataConversion/foamDataToFluent/writeFluentScalarField.C
index 5401b662b9e637fb83bdb9cce2ac01e39da36550..017be50b6c2fe1143e070b280ed57cbbcdcab6a0 100644
--- a/applications/utilities/postProcessing/dataConversion/foamDataToFluent/writeFluentScalarField.C
+++ b/applications/utilities/postProcessing/dataConversion/foamDataToFluent/writeFluentScalarField.C
@@ -45,7 +45,7 @@ void writeFluentField
     Ostream& stream
 )
 {
-    const scalarField& phiInternal = phi.internalField();
+    const scalarField& phiInternal = phi;
 
     // Writing cells
     stream
@@ -75,7 +75,7 @@ void writeFluentField
         {
             // Form empty patch field repeat the internal field to
             // allow for the node interpolation in Fluent
-            const scalarField& phiInternal = phi.internalField();
+            const scalarField& phiInternal = phi;
 
             // Get reference to internal cells
             const labelList emptyFaceCells =
diff --git a/applications/utilities/postProcessing/dataConversion/foamDataToFluent/writeFluentVectorField.C b/applications/utilities/postProcessing/dataConversion/foamDataToFluent/writeFluentVectorField.C
index 24fb46b8079ae200a6a59ebae7d1e75b79cabcf1..f5953c2994f6abcea67cf53ee21ce12c6968b288 100644
--- a/applications/utilities/postProcessing/dataConversion/foamDataToFluent/writeFluentVectorField.C
+++ b/applications/utilities/postProcessing/dataConversion/foamDataToFluent/writeFluentVectorField.C
@@ -44,7 +44,7 @@ void writeFluentField
     Ostream& stream
 )
 {
-    const vectorField& phiInternal = phi.internalField();
+    const vectorField& phiInternal = phi;
 
     // Writing cells
     stream
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
index 02348e3e5570a9617602c2e71a6ed86675c0acd4..f42443485676e4fc18d8e4ed5311c536febac472 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
@@ -320,7 +320,7 @@ void Foam::MRFZone::addCoriolis
 
     const labelList& cells = mesh_.cellZones()[cellZoneID_];
     vectorField& ddtUc = ddtU.internalField();
-    const vectorField& Uc = U.internalField();
+    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 165dfe7b8af8d81a5a0f7981322f41bfdd8b8606..8d681da85c5e3acf41fdd43675f282e7f46837ed 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
@@ -43,8 +43,8 @@ void Foam::MRFZone::makeRelativeRhoFlux
 
     const vector Omega = omega_->value(mesh_.time().timeOutputValue())*axis_;
 
-    const vectorField& Cfi = Cf.internalField();
-    const vectorField& Sfi = Sf.internalField();
+    const vectorField& Cfi = Cf;
+    const vectorField& Sfi = Sf;
     scalarField& phii = phi.internalField();
 
     // Internal faces
@@ -143,8 +143,8 @@ void Foam::MRFZone::makeAbsoluteRhoFlux
 
     const vector Omega = omega_->value(mesh_.time().timeOutputValue())*axis_;
 
-    const vectorField& Cfi = Cf.internalField();
-    const vectorField& Sfi = Sf.internalField();
+    const vectorField& Cfi = Cf;
+    const vectorField& Sfi = Sf;
     scalarField& phii = phi.internalField();
 
     // Internal faces
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
index da4c454a506c4afd4ddf2143e6202897eb52fc34..5de6d069e9fb828dfcbfa77ce83ce7a7790f5207 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
@@ -351,7 +351,7 @@ localEulerDdtScheme<Type>::fvmDdt
 
     fvMatrix<Type>& fvm = tfvm.ref();
 
-    const scalarField& rDeltaT = localRDeltaT().internalField();
+    const scalarField& rDeltaT = localRDeltaT();
 
     fvm.diag() = rDeltaT*mesh().Vsc();
 
@@ -386,7 +386,7 @@ localEulerDdtScheme<Type>::fvmDdt
     );
     fvMatrix<Type>& fvm = tfvm.ref();
 
-    const scalarField& rDeltaT = localRDeltaT().internalField();
+    const scalarField& rDeltaT = localRDeltaT();
 
     fvm.diag() = rDeltaT*rho.value()*mesh().Vsc();
 
@@ -423,7 +423,7 @@ localEulerDdtScheme<Type>::fvmDdt
     );
     fvMatrix<Type>& fvm = tfvm.ref();
 
-    const scalarField& rDeltaT = localRDeltaT().internalField();
+    const scalarField& rDeltaT = localRDeltaT();
 
     fvm.diag() = rDeltaT*rho.internalField()*mesh().Vsc();
 
@@ -463,7 +463,7 @@ localEulerDdtScheme<Type>::fvmDdt
     );
     fvMatrix<Type>& fvm = tfvm.ref();
 
-    const scalarField& rDeltaT = localRDeltaT().internalField();
+    const scalarField& rDeltaT = localRDeltaT();
 
     fvm.diag() = rDeltaT*alpha.internalField()*rho.internalField()*mesh().Vsc();
 
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C
index b0eb52a1080a74996ba2fef2cd6d8fbb5461dc18..4e7b8b83f95af8a69d147c7c5e74ba995aea312f 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/faceCorrectedSnGrad/faceCorrectedSnGrad.C
@@ -72,9 +72,9 @@ Foam::fv::faceCorrectedSnGrad<Type>::fullGradCorrection
 
     const pointField& points = mesh.points();
     const faceList& faces = mesh.faces();
-    const vectorField& Sf = mesh.Sf().internalField();
-    const vectorField& C = mesh.C().internalField();
-    const scalarField& magSf = mesh.magSf().internalField();
+    const vectorField& Sf = mesh.Sf();
+    const vectorField& C = mesh.C();
+    const scalarField& magSf = mesh.magSf();
     const labelList& owner = mesh.owner();
     const labelList& neighbour = mesh.neighbour();
 
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C b/src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C
index 5994fbdcf931a231f41c708fe31bc634e58cae57..2ad75944f449153e66d68ad9e3c115d8d6368b8e 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C
@@ -126,7 +126,7 @@ snGradScheme<Type>::snGrad
     GeometricField<Type, fvsPatchField, surfaceMesh>& ssf = tsf.ref();
 
     // set reference to difference factors array
-    const scalarField& deltaCoeffs = tdeltaCoeffs().internalField();
+    const scalarField& deltaCoeffs = tdeltaCoeffs();
 
     // owner/neighbour addressing
     const labelUList& owner = mesh.owner();
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C
index 429b36ba1b591cf0c70ce7c0088ad38aab73da69..67d865b18b73ca2b430de95609e1a1d35ad14d0e 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C
@@ -51,9 +51,9 @@ correction
     Field<Type>& sfCorr = tsfCorr.ref().internalField();
 
     const pointField& points = mesh.points();
-    const pointField& C = mesh.C().internalField();
+    const pointField& C = mesh.C();
     const faceList& faces = mesh.faces();
-    const scalarField& w = mesh.weights().internalField();
+    const scalarField& w = mesh.weights();
     const labelList& owner = mesh.owner();
     const labelList& neighbour = mesh.neighbour();
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
index 718949afbd8cb30e85a54c0949b950d185746a90..ed7b2576c60a9e28371507303fa8cf5eecdde7db 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C
@@ -157,9 +157,9 @@ Foam::surfaceInterpolationScheme<Type>::interpolate
     const surfaceScalarField& lambdas = tlambdas();
     const surfaceScalarField& ys = tys();
 
-    const Field<Type>& vfi = vf.internalField();
-    const scalarField& lambda = lambdas.internalField();
-    const scalarField& y = ys.internalField();
+    const Field<Type>& vfi = vf;
+    const scalarField& lambda = lambdas;
+    const scalarField& y = ys;
 
     const fvMesh& mesh = vf.mesh();
     const labelUList& P = mesh.owner();
@@ -251,8 +251,8 @@ Foam::surfaceInterpolationScheme<Type>::dotInterpolate
 
     const surfaceScalarField& lambdas = tlambdas();
 
-    const Field<Type>& vfi = vf.internalField();
-    const scalarField& lambda = lambdas.internalField();
+    const Field<Type>& vfi = vf;
+    const scalarField& lambda = lambdas;
 
     const fvMesh& mesh = vf.mesh();
     const labelUList& P = mesh.owner();
diff --git a/src/fvOptions/sources/derived/meanVelocityForce/meanVelocityForce.C b/src/fvOptions/sources/derived/meanVelocityForce/meanVelocityForce.C
index 24ca385704752b6a78b22be66e936249929cfa5c..64f73074382d3a1c68dbe0bb00f91288673b8193 100644
--- a/src/fvOptions/sources/derived/meanVelocityForce/meanVelocityForce.C
+++ b/src/fvOptions/sources/derived/meanVelocityForce/meanVelocityForce.C
@@ -147,7 +147,7 @@ Foam::scalar Foam::fv::meanVelocityForce::magUbarAve
 
 void Foam::fv::meanVelocityForce::correct(volVectorField& U)
 {
-    const scalarField& rAU = rAPtr_().internalField();
+    const scalarField& rAU = rAPtr_();
 
     // Integrate flow variables over cell set
     scalar rAUave = 0.0;
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
index e4f28920a40a2aa4584f1a108ee5926643e4bd1d..fa478e3f30440b4ff00595e6a34899d5d62eb25d 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
@@ -1081,8 +1081,8 @@ void kinematicSingleLayer::info()
 {
     Info<< "\nSurface film: " << type() << endl;
 
-    const scalarField& deltaInternal = delta_.internalField();
-    const vectorField& Uinternal = U_.internalField();
+    const scalarField& deltaInternal = delta_;
+    const vectorField& Uinternal = U_;
     scalar addedMassTotal = 0.0;
     outputProperties().readIfPresent("addedMassTotal", addedMassTotal);
     addedMassTotal += returnReduce(addedMassTotal_, sumOp<scalar>());
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.C
index 489b453cb90df3022efb1fb398b940638d5774a4..807ca3ad9111b04d59dee6e5522d6d4b537e3745 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.C
@@ -125,8 +125,8 @@ tmp<volScalarField> constantRadiation::Shs()
     if ((time >= timeStart_) && (time <= timeStart_ + duration_))
     {
         scalarField& Shs = tShs.ref();
-        const scalarField& Qr = QrConst_.internalField();
-        const scalarField& alpha = owner_.alpha().internalField();
+        const scalarField& Qr = QrConst_;
+        const scalarField& alpha = owner_.alpha();
 
         Shs = mask_*Qr*alpha*absorptivity_;
     }
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C
index a6350c3ac18131e368eafd92be0a39caa5d3fb68..55ea60d3224cf3455e7020ce0e0b4459e093edb0 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C
@@ -108,8 +108,8 @@ tmp<volScalarField> primaryRadiation::Shs()
     );
 
     scalarField& Shs = tShs.ref();
-    const scalarField& QinP = QinPrimary_.internalField();
-    const scalarField& alpha = owner_.alpha().internalField();
+    const scalarField& QinP = QinPrimary_;
+    const scalarField& alpha = owner_.alpha();
 
     Shs = QinP*alpha;
 
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
index 2d6e8b77b8f484ac4d3884f9a0fb6a909ca27624..ab51ea57e8b24376d8e322606b3079cb2595b12a 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
@@ -125,9 +125,9 @@ tmp<volScalarField> standardRadiation::Shs()
     );
 
     scalarField& Shs = tShs.ref();
-    const scalarField& QinP = QinPrimary_.internalField();
-    const scalarField& delta = owner_.delta().internalField();
-    const scalarField& alpha = owner_.alpha().internalField();
+    const scalarField& QinP = QinPrimary_;
+    const scalarField& delta = owner_.delta();
+    const scalarField& alpha = owner_.alpha();
 
     Shs = beta_*QinP*alpha*(1.0 - exp(-kappaBar_*delta));
 
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
index 59cfd243eac6ead41f8af4770579bbdf47faf64c..ed59cb375792432a20e05dc774ccc93cf85a855f 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
@@ -727,7 +727,7 @@ void thermoSingleLayer::info()
 {
     kinematicSingleLayer::info();
 
-    const scalarField& Tinternal = T_.internalField();
+    const scalarField& Tinternal = T_;
 
     Info<< indent << "min/mean/max(T)    = "
         << gMin(Tinternal) << ", "
diff --git a/src/sampling/meshToMesh0/calculateMeshToMesh0Weights.C b/src/sampling/meshToMesh0/calculateMeshToMesh0Weights.C
index 277e5e6d91d605b345499d9c16e5ba809adc502f..651897ba416edbcd9affede521a9f142d76cff9b 100644
--- a/src/sampling/meshToMesh0/calculateMeshToMesh0Weights.C
+++ b/src/sampling/meshToMesh0/calculateMeshToMesh0Weights.C
@@ -51,8 +51,8 @@ void Foam::meshToMesh0::calculateInverseDistanceWeights() const
 
     // get reference to source mesh data
     const labelListList& cc = fromMesh_.cellCells();
-    const vectorField& centreFrom = fromMesh_.C().internalField();
-    const vectorField& centreTo = toMesh_.C().internalField();
+    const vectorField& centreFrom = fromMesh_.C();
+    const vectorField& centreTo = toMesh_.C();
 
     forAll(cellAddressing_, celli)
     {
diff --git a/src/thermophysicalModels/basic/heThermo/heThermo.C b/src/thermophysicalModels/basic/heThermo/heThermo.C
index 9ac59bfedc250bc45f92eda2f4795f78c960261a..969f6ca1fa418951fe59dd78ba21a82bfdbe9380 100644
--- a/src/thermophysicalModels/basic/heThermo/heThermo.C
+++ b/src/thermophysicalModels/basic/heThermo/heThermo.C
@@ -55,8 +55,8 @@ template<class BasicThermo, class MixtureType>
 void Foam::heThermo<BasicThermo, MixtureType>::init()
 {
     scalarField& heCells = he_.internalField();
-    const scalarField& pCells = this->p_.internalField();
-    const scalarField& TCells = this->T_.internalField();
+    const scalarField& pCells = this->p_;
+    const scalarField& TCells = this->T_;
 
     forAll(heCells, celli)
     {
@@ -188,8 +188,8 @@ Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::he
 
     volScalarField& he = the.ref();
     scalarField& heCells = he.internalField();
-    const scalarField& pCells = p.internalField();
-    const scalarField& TCells = T.internalField();
+    const scalarField& pCells = p;
+    const scalarField& TCells = T;
 
     forAll(heCells, celli)
     {
diff --git a/src/thermophysicalModels/basic/psiThermo/hePsiThermo.C b/src/thermophysicalModels/basic/psiThermo/hePsiThermo.C
index b73f3207d2e6ffd33c491cfbc171dfbb7030c9c1..99f228a1596840910cf7d8aa2e0ebc46852a9403 100644
--- a/src/thermophysicalModels/basic/psiThermo/hePsiThermo.C
+++ b/src/thermophysicalModels/basic/psiThermo/hePsiThermo.C
@@ -30,8 +30,8 @@ License
 template<class BasicPsiThermo, class MixtureType>
 void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()
 {
-    const scalarField& hCells = this->he_.internalField();
-    const scalarField& pCells = this->p_.internalField();
+    const scalarField& hCells = this->he_;
+    const scalarField& pCells = this->p_;
 
     scalarField& TCells = this->T_.internalField();
     scalarField& psiCells = this->psi_.internalField();
diff --git a/src/thermophysicalModels/basic/rhoThermo/heRhoThermo.C b/src/thermophysicalModels/basic/rhoThermo/heRhoThermo.C
index ec0a220d6da1f2c83b1da1743f5edfbd38a35f13..febff08db1cbdbaefe1f1fee1561a51921d2a570 100644
--- a/src/thermophysicalModels/basic/rhoThermo/heRhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/heRhoThermo.C
@@ -30,8 +30,8 @@ License
 template<class BasicPsiThermo, class MixtureType>
 void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::calculate()
 {
-    const scalarField& hCells = this->he().internalField();
-    const scalarField& pCells = this->p_.internalField();
+    const scalarField& hCells = this->he();
+    const scalarField& pCells = this->p_;
 
     scalarField& TCells = this->T_.internalField();
     scalarField& psiCells = this->psi_.internalField();
diff --git a/src/thermophysicalModels/reactionThermo/psiuReactionThermo/heheuPsiThermo.C b/src/thermophysicalModels/reactionThermo/psiuReactionThermo/heheuPsiThermo.C
index 6c3d0f1ff6d10d9eb68268fa2b36ecbedb3b4074..703e088005c882a071b4f62c7fd1c15f3921b026 100644
--- a/src/thermophysicalModels/reactionThermo/psiuReactionThermo/heheuPsiThermo.C
+++ b/src/thermophysicalModels/reactionThermo/psiuReactionThermo/heheuPsiThermo.C
@@ -32,9 +32,9 @@ License
 template<class BasicPsiThermo, class MixtureType>
 void Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::calculate()
 {
-    const scalarField& hCells = this->he_.internalField();
-    const scalarField& heuCells = this->heu_.internalField();
-    const scalarField& pCells = this->p_.internalField();
+    const scalarField& hCells = this->he_;
+    const scalarField& heuCells = this->heu_;
+    const scalarField& pCells = this->p_;
 
     scalarField& TCells = this->T_.internalField();
     scalarField& TuCells = this->Tu_.internalField();
@@ -177,8 +177,8 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::heheuPsiThermo
     )
 {
     scalarField& heuCells = this->heu_.internalField();
-    const scalarField& pCells = this->p_.internalField();
-    const scalarField& TuCells = this->Tu_.internalField();
+    const scalarField& pCells = this->p_;
+    const scalarField& TuCells = this->Tu_;
 
     forAll(heuCells, celli)
     {
@@ -309,9 +309,9 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::Tb() const
 
     volScalarField& Tb_ = tTb.ref();
     scalarField& TbCells = Tb_.internalField();
-    const scalarField& pCells = this->p_.internalField();
-    const scalarField& TCells = this->T_.internalField();
-    const scalarField& hCells = this->he_.internalField();
+    const scalarField& pCells = this->p_;
+    const scalarField& TCells = this->T_;
+    const scalarField& hCells = this->he_;
 
     forAll(TbCells, celli)
     {
@@ -369,8 +369,8 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::psiu() const
 
     volScalarField& psiu = tpsiu.ref();
     scalarField& psiuCells = psiu.internalField();
-    const scalarField& TuCells = this->Tu_.internalField();
-    const scalarField& pCells = this->p_.internalField();
+    const scalarField& TuCells = this->Tu_;
+    const scalarField& pCells = this->p_;
 
     forAll(psiuCells, celli)
     {
@@ -424,8 +424,8 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::psib() const
     volScalarField& psib = tpsib.ref();
     scalarField& psibCells = psib.internalField();
     const volScalarField Tb_(Tb());
-    const scalarField& TbCells = Tb_.internalField();
-    const scalarField& pCells = this->p_.internalField();
+    const scalarField& TbCells = Tb_;
+    const scalarField& pCells = this->p_;
 
     forAll(psibCells, celli)
     {
@@ -478,8 +478,8 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::muu() const
 
     volScalarField& muu_ = tmuu.ref();
     scalarField& muuCells = muu_.internalField();
-    const scalarField& pCells = this->p_.internalField();
-    const scalarField& TuCells = this->Tu_.internalField();
+    const scalarField& pCells = this->p_;
+    const scalarField& TuCells = this->Tu_;
 
     forAll(muuCells, celli)
     {
@@ -537,8 +537,8 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::mub() const
     volScalarField& mub_ = tmub.ref();
     scalarField& mubCells = mub_.internalField();
     const volScalarField Tb_(Tb());
-    const scalarField& pCells = this->p_.internalField();
-    const scalarField& TbCells = Tb_.internalField();
+    const scalarField& pCells = this->p_;
+    const scalarField& TbCells = Tb_;
 
     forAll(mubCells, celli)
     {
diff --git a/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C b/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C
index 7dc39c3adcb11045b901405eb1a7ae92eff20ee6..b6b4bea3be2339dadccefb091612822e736e8166 100644
--- a/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C
+++ b/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C
@@ -33,8 +33,8 @@ void Foam::heSolidThermo<BasicSolidThermo, MixtureType>::calculate()
 {
     scalarField& TCells = this->T_.internalField();
 
-    const scalarField& hCells = this->he_.internalField();
-    const scalarField& pCells = this->p_.internalField();
+    const scalarField& hCells = this->he_;
+    const scalarField& pCells = this->p_;
     scalarField& rhoCells = this->rho_.internalField();
     scalarField& alphaCells = this->alpha_.internalField();
 
@@ -219,8 +219,8 @@ Foam::heSolidThermo<BasicSolidThermo, MixtureType>::Kappa() const
 
     volVectorField& Kappa = tKappa.ref();
     vectorField& KappaCells = Kappa.internalField();
-    const scalarField& TCells = this->T_.internalField();
-    const scalarField& pCells = this->p_.internalField();
+    const scalarField& TCells = this->T_;
+    const scalarField& pCells = this->p_;
 
     forAll(KappaCells, celli)
     {