From 22f4ad32b174e219b57105cbe8b4c4c711a24c3d Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Mon, 25 Apr 2016 16:16:05 +0100 Subject: [PATCH] Completed boundaryField() -> boundaryFieldRef() Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1938 Because C++ does not support overloading based on the return-type there is a problem defining both const and non-const member functions which are resolved based on the const-ness of the object for which they are called rather than the intent of the programmer declared via the const-ness of the returned type. The issue for the "boundaryField()" member function is that the non-const version increments the event-counter and checks the state of the stored old-time fields in case the returned value is altered whereas the const version has no side-effects and simply returns the reference. If the the non-const function is called within the patch-loop the event-counter may overflow. To resolve this it in necessary to avoid calling the non-const form of "boundaryField()" if the results is not altered and cache the reference outside the patch-loop when mutation of the patch fields is needed. The most straight forward way of resolving this problem is to name the const and non-const forms of the member functions differently e.g. the non-const form could be named: mutableBoundaryField() mutBoundaryField() nonConstBoundaryField() boundaryFieldRef() Given that in C++ a reference is non-const unless specified as const: "T&" vs "const T&" the logical convention would be boundaryFieldRef() boundaryFieldConstRef() and given that the const form which is more commonly used is it could simply be named "boundaryField()" then the logical convention is GeometricBoundaryField& boundaryFieldRef(); inline const GeometricBoundaryField& boundaryField() const; This is also consistent with the new "tmp" class for which non-const access to the stored object is obtained using the ".ref()" member function. This new convention for non-const access to the components of GeometricField will be applied to "dimensionedInternalField()" and "internalField()" in the future, i.e. "dimensionedInternalFieldRef()" and "internalFieldRef()". --- .../turbulence/PDRkEpsilon/PDRkEpsilon.C | 4 +-- .../XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C | 4 ++- .../SCOPE/SCOPELaminarFlameSpeed.C | 18 ++++++++---- .../rhoCentralDyMFoam/rhoCentralDyMFoam.C | 6 ++-- .../rhoCentralFoam/rhoCentralFoam.C | 6 ++-- .../multiphaseMixtureThermo.C | 2 +- .../solvers/multiphase/interFoam/alphaEqn.H | 5 +++- .../threePhaseInterfaceProperties.C | 4 +-- .../multiphaseSystem/multiphaseSystem.C | 19 +++++++++---- .../multiphase/multiphaseEulerFoam/pEqn.H | 2 +- .../multiphaseMixture/multiphaseMixture.C | 2 +- .../saturationModels/polynomial/polynomial.C | 4 ++- .../wallLubricationModel.C | 4 ++- .../BlendedInterfacialModel.C | 5 +++- .../ThermalPhaseChangePhaseSystem.C | 4 +-- .../multiphaseSystem/multiphaseSystem.C | 8 ++++-- .../reactingMultiphaseEulerFoam/pU/pEqn.H | 7 +++-- .../reactingTwoPhaseEulerFoam/pU/pEqn.H | 28 ++++++++++++------- .../reactingTwoPhaseEulerFoam/pUf/pEqn.H | 2 +- .../JohnsonJacksonSchaefferFrictionalStress.C | 4 ++- .../Schaeffer/SchaefferFrictionalStress.C | 4 ++- .../kineticTheoryModel/kineticTheoryModel.C | 2 +- .../phasePressureModel/phasePressureModel.C | 4 +-- .../twoPhaseSystem/twoPhaseSystem.C | 8 ++++-- .../multiphase/twoPhaseEulerFoam/pU/pEqn.H | 28 ++++++++++++------- .../multiphase/twoPhaseEulerFoam/pUf/pEqn.H | 2 +- .../kineticTheoryModel/kineticTheoryModel.C | 2 +- .../phasePressureModel/phasePressureModel.C | 4 +-- .../BlendedInterfacialModel.C | 5 +++- .../twoPhaseSystem/twoPhaseSystem.C | 10 ++++--- .../GeometricField/GeometricField.H | 26 ++++++++--------- .../RAS/mixtureKEpsilon/mixtureKEpsilon.C | 2 +- .../constrainPressure/constrainPressure.C | 2 +- .../fvMatrices/solvers/MULES/MULESTemplates.C | 4 +-- 34 files changed, 148 insertions(+), 93 deletions(-) diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C index afe97f7f7f..4a2a35635b 100644 --- a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C +++ b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C @@ -130,7 +130,7 @@ void PDRkEpsilon::correct() tgradU.clear(); // Update espsilon and G at the wall - epsilon_.boundaryField().updateCoeffs(); + epsilon_.boundaryFieldRef().updateCoeffs(); // Add the blockage generation term so that it is included consistently // in both the k and epsilon equations @@ -163,7 +163,7 @@ void PDRkEpsilon::correct() epsEqn.ref().relax(); - epsEqn.ref().boundaryManipulate(epsilon_.boundaryField()); + epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef()); solve(epsEqn); bound(epsilon_, epsilonMin_); diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C index 0208818ee0..6ca24555b4 100644 --- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C +++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C @@ -119,9 +119,11 @@ Foam::tmp Foam::XiEqModels::SCOPEXiEq::XiEq() const } } + volScalarField::GeometricBoundaryField& xieqBf = xieq.boundaryFieldRef(); + forAll(xieq.boundaryField(), patchi) { - scalarField& xieqp = xieq.boundaryField()[patchi]; + scalarField& xieqp = xieqBf[patchi]; const scalarField& Kp = K.boundaryField()[patchi]; const scalarField& Map = Ma.boundaryField()[patchi]; const scalarField& upBySup = upBySu.boundaryField()[patchi]; diff --git a/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C b/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C index fd07b0b3cb..17f92a85f2 100644 --- a/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C +++ b/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C @@ -266,9 +266,11 @@ Foam::tmp Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi); } - forAll(Su0.boundaryField(), patchi) + volScalarField::GeometricBoundaryField& Su0Bf = Su0.boundaryFieldRef(); + + forAll(Su0Bf, patchi) { - scalarField& Su0p = Su0.boundaryField()[patchi]; + scalarField& Su0p = Su0Bf[patchi]; const scalarField& pp = p.boundaryField()[patchi]; const scalarField& Tup = Tu.boundaryField()[patchi]; @@ -313,9 +315,11 @@ Foam::tmp Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi[celli]); } - forAll(Su0.boundaryField(), patchi) + volScalarField::GeometricBoundaryField& Su0Bf = Su0.boundaryFieldRef(); + + forAll(Su0Bf, patchi) { - scalarField& Su0p = Su0.boundaryField()[patchi]; + scalarField& Su0p = Su0Bf[patchi]; const scalarField& pp = p.boundaryField()[patchi]; const scalarField& Tup = Tu.boundaryField()[patchi]; const scalarField& phip = phi.boundaryField()[patchi]; @@ -365,9 +369,11 @@ Foam::tmp Foam::laminarFlameSpeedModels::SCOPE::Ma ma[celli] = Ma(phi[celli]); } - forAll(ma.boundaryField(), patchi) + volScalarField::GeometricBoundaryField& maBf = ma.boundaryFieldRef(); + + forAll(maBf, patchi) { - scalarField& map = ma.boundaryField()[patchi]; + scalarField& map = maBf[patchi]; const scalarField& phip = phi.boundaryField()[patchi]; forAll(map, facei) diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C index 3b43b5dcec..fb0fbc7079 100644 --- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C +++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C @@ -189,7 +189,7 @@ int main(int argc, char *argv[]) rhoU.dimensionedInternalField() /rho.dimensionedInternalField(); U.correctBoundaryConditions(); - rhoU.boundaryField() == rho.boundaryField()*U.boundaryField(); + rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); if (!inviscid) { @@ -223,7 +223,7 @@ int main(int argc, char *argv[]) e = rhoE/rho - 0.5*magSqr(U); e.correctBoundaryConditions(); thermo.correct(); - rhoE.boundaryField() == + rhoE.boundaryFieldRef() == rho.boundaryField()* ( e.boundaryField() + 0.5*magSqr(U.boundaryField()) @@ -244,7 +244,7 @@ int main(int argc, char *argv[]) rho.dimensionedInternalField() /psi.dimensionedInternalField(); p.correctBoundaryConditions(); - rho.boundaryField() == psi.boundaryField()*p.boundaryField(); + rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField(); turbulence->correct(); diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C index 3754632d3d..083a1b1b5c 100644 --- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C +++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C @@ -182,7 +182,7 @@ int main(int argc, char *argv[]) rhoU.dimensionedInternalField() /rho.dimensionedInternalField(); U.correctBoundaryConditions(); - rhoU.boundaryField() == rho.boundaryField()*U.boundaryField(); + rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); if (!inviscid) { @@ -216,7 +216,7 @@ int main(int argc, char *argv[]) e = rhoE/rho - 0.5*magSqr(U); e.correctBoundaryConditions(); thermo.correct(); - rhoE.boundaryField() == + rhoE.boundaryFieldRef() == rho.boundaryField()* ( e.boundaryField() + 0.5*magSqr(U.boundaryField()) @@ -237,7 +237,7 @@ int main(int argc, char *argv[]) rho.dimensionedInternalField() /psi.dimensionedInternalField(); p.correctBoundaryConditions(); - rho.boundaryField() == psi.boundaryField()*p.boundaryField(); + rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField(); turbulence->correct(); diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C index 14bbbe8708..db4b5972e9 100644 --- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C @@ -923,7 +923,7 @@ Foam::tmp Foam::multiphaseMixtureThermo::K { tmp tnHatfv = nHatfv(alpha1, alpha2); - correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryField()); + correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef()); // Simple expression for curvature return -fvc::div(tnHatfv & mesh_.Sf()); diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H index 9d6dd2dc18..f75f4d0f51 100644 --- a/applications/solvers/multiphase/interFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interFoam/alphaEqn.H @@ -54,11 +54,14 @@ phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U)); } + surfaceScalarField::GeometricBoundaryField& phicBf = + phic.boundaryFieldRef(); + // Do not compress interface at non-coupled boundary faces // (inlets, outlets etc.) forAll(phic.boundaryField(), patchi) { - fvsPatchScalarField& phicp = phic.boundaryField()[patchi]; + fvsPatchScalarField& phicp = phicBf[patchi]; if (!phicp.coupled()) { diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C b/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C index 8343ceb790..d5ee2f25f8 100644 --- a/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C +++ b/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.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 @@ -134,7 +134,7 @@ void Foam::threePhaseInterfaceProperties::calculateK() // Face unit interface normal surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_)); - correctContactAngle(nHatfv.boundaryField()); + correctContactAngle(nHatfv.boundaryFieldRef()); // Face unit interface normal flux nHatf_ = nHatfv & Sf; diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index 3f514904ec..62e2a97a30 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -124,11 +124,13 @@ void Foam::multiphaseSystem::solveAlphas() ); } + surfaceScalarField::GeometricBoundaryField& alphaPhiCorrBf = + alphaPhiCorr.boundaryFieldRef(); + // Ensure that the flux at inflow BCs is preserved - forAll(alphaPhiCorr.boundaryField(), patchi) + forAll(alphaPhiCorrBf, patchi) { - fvsPatchScalarField& alphaPhiCorrp = - alphaPhiCorr.boundaryField()[patchi]; + fvsPatchScalarField& alphaPhiCorrp = alphaPhiCorrBf[patchi]; if (!alphaPhiCorrp.coupled()) { @@ -372,7 +374,7 @@ Foam::tmp Foam::multiphaseSystem::K { tmp tnHatfv = nHatfv(phase1, phase2); - correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryField()); + correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef()); // Simple expression for curvature return -fvc::div(tnHatfv & mesh_.Sf()); @@ -666,6 +668,9 @@ Foam::tmp Foam::multiphaseSystem::Svm } } + volVectorField::GeometricBoundaryField& SvmBf = + tSvm.ref().boundaryFieldRef(); + // Remove virtual mass at fixed-flux boundaries forAll(phase.phi().boundaryField(), patchi) { @@ -677,7 +682,7 @@ Foam::tmp Foam::multiphaseSystem::Svm ) ) { - tSvm.ref().boundaryField()[patchi] = Zero; + SvmBf[patchi] = Zero; } } @@ -713,6 +718,8 @@ Foam::multiphaseSystem::dragCoeffs() const ) ).ptr(); + volScalarField::GeometricBoundaryField& Kbf = Kptr->boundaryFieldRef(); + // Remove drag at fixed-flux boundaries forAll(dm.phase1().phi().boundaryField(), patchi) { @@ -724,7 +731,7 @@ Foam::multiphaseSystem::dragCoeffs() const ) ) { - Kptr->boundaryField()[patchi] = 0.0; + Kbf[patchi] = 0.0; } } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H index bbcc36f4e7..f65fdf8de1 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H @@ -203,7 +203,7 @@ setSnGrad ( - p_rgh.boundaryField(), + p_rgh.boundaryFieldRef(), ( phiHbyA.boundaryField() - MRF.relative(phib) )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C index ea08ba2986..f1fed0103d 100644 --- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C +++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C @@ -523,7 +523,7 @@ Foam::tmp Foam::multiphaseMixture::K { tmp tnHatfv = nHatfv(alpha1, alpha2); - correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryField()); + correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef()); // Simple expression for curvature return -fvc::div(tnHatfv & mesh_.Sf()); diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/polynomial/polynomial.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/polynomial/polynomial.C index a2fbcd0bbc..886da7c200 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/polynomial/polynomial.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/polynomial/polynomial.C @@ -118,9 +118,11 @@ Foam::saturationModels::polynomial::Tsat Tsat[celli] = C_.value(p[celli]); } + volScalarField::GeometricBoundaryField& TsatBf = Tsat.boundaryFieldRef(); + forAll(Tsat.boundaryField(), patchi) { - scalarField& Tsatp = Tsat.boundaryField()[patchi]; + scalarField& Tsatp = TsatBf[patchi]; const scalarField& pp = p.boundaryField()[patchi]; forAll(Tsatp, facei) diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C index a2a95cce05..784f323e0c 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C @@ -50,11 +50,13 @@ Foam::tmp Foam::wallLubricationModel::zeroGradWalls volVectorField& Fi = tFi.ref(); const fvPatchList& patches = Fi.mesh().boundary(); + volVectorField::GeometricBoundaryField& FiBf = Fi.boundaryFieldRef(); + forAll(patches, patchi) { if (isA(patches[patchi])) { - fvPatchVectorField& Fiw = Fi.boundaryField()[patchi]; + fvPatchVectorField& Fiw = FiBf[patchi]; Fiw = Fiw.patchInternalField(); } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C index f216859b08..a19b995bc5 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C @@ -36,6 +36,9 @@ void Foam::BlendedInterfacialModel::correctFixedFluxBCs GeometricField& field ) const { + typename GeometricField::GeometricBoundaryField& fieldBf = + field.boundaryFieldRef(); + forAll(phase1_.phi()().boundaryField(), patchi) { if @@ -46,7 +49,7 @@ void Foam::BlendedInterfacialModel::correctFixedFluxBCs ) ) { - field.boundaryField()[patchi] = Zero; + fieldBf[patchi] = Zero; } } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C index ff1e0b0998..8ba5a4a502 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C @@ -353,9 +353,9 @@ void Foam::ThermalPhaseChangePhaseSystem::correctThermo() // Limit the H[12] boundary field to avoid /0 const scalar HLimit = 1e-4; - H1.boundaryField() = + H1.boundaryFieldRef() = max(H1.boundaryField(), phase1.boundaryField()*HLimit); - H2.boundaryField() = + H2.boundaryFieldRef() = max(H2.boundaryField(), phase2.boundaryField()*HLimit); volScalarField mDotL diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index b81e174db8..6504f943a0 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -134,11 +134,13 @@ void Foam::multiphaseSystem::solveAlphas() ); } + surfaceScalarField::GeometricBoundaryField& alphaPhiCorrBf = + alphaPhiCorr.boundaryFieldRef(); + // Ensure that the flux at inflow BCs is preserved forAll(alphaPhiCorr.boundaryField(), patchi) { - fvsPatchScalarField& alphaPhiCorrp = - alphaPhiCorr.boundaryField()[patchi]; + fvsPatchScalarField& alphaPhiCorrp = alphaPhiCorrBf[patchi]; if (!alphaPhiCorrp.coupled()) { @@ -477,7 +479,7 @@ Foam::tmp Foam::multiphaseSystem::K { tmp tnHatfv = nHatfv(phase1, phase2); - correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryField()); + correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef()); // Simple expression for curvature return -fvc::div(tnHatfv & mesh_.Sf()); diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H index 05a049b4c6..e078414a20 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H @@ -182,6 +182,9 @@ while (pimple.correct()) ); surfaceScalarField phiCorrCoeff(pos(alphafBar - 0.99)); + surfaceScalarField::GeometricBoundaryField& phiCorrCoeffBf = + phiCorrCoeff.boundaryFieldRef(); + forAll(mesh.boundary(), patchi) { // Set ddtPhiCorr to 0 on non-coupled boundaries @@ -191,7 +194,7 @@ while (pimple.correct()) || isA(mesh.boundary()[patchi]) ) { - phiCorrCoeff.boundaryField()[patchi] = 0; + phiCorrCoeffBf[patchi] = 0; } } @@ -281,7 +284,7 @@ while (pimple.correct()) setSnGrad ( - p_rgh.boundaryField(), + p_rgh.boundaryFieldRef(), ( phiHbyA.boundaryField() - phib )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H index ba312ab248..008134c21a 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H @@ -149,17 +149,25 @@ while (pimple.correct()) surfaceScalarField phiCorrCoeff1(pos(alphaf1Bar - 0.99)); surfaceScalarField phiCorrCoeff2(pos(0.01 - alphaf1Bar)); - forAll(mesh.boundary(), patchi) { - // Set ddtPhiCorr to 0 on non-coupled boundaries - if - ( - !mesh.boundary()[patchi].coupled() - || isA(mesh.boundary()[patchi]) - ) + surfaceScalarField::GeometricBoundaryField& phiCorrCoeff1Bf = + phiCorrCoeff1.boundaryFieldRef(); + + surfaceScalarField::GeometricBoundaryField& phiCorrCoeff2Bf = + phiCorrCoeff2.boundaryFieldRef(); + + forAll(mesh.boundary(), patchi) { - phiCorrCoeff1.boundaryField()[patchi] = 0; - phiCorrCoeff2.boundaryField()[patchi] = 0; + // Set ddtPhiCorr to 0 on non-coupled boundaries + if + ( + !mesh.boundary()[patchi].coupled() + || isA(mesh.boundary()[patchi]) + ) + { + phiCorrCoeff1Bf[patchi] = 0; + phiCorrCoeff2Bf[patchi] = 0; + } } } @@ -215,7 +223,7 @@ while (pimple.correct()) // Update the fixedFluxPressure BCs to ensure flux consistency setSnGrad ( - p_rgh.boundaryField(), + p_rgh.boundaryFieldRef(), ( phiHbyA.boundaryField() - ( diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H index 444f2de81f..01c120d481 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H @@ -201,7 +201,7 @@ while (pimple.correct()) // Update the fixedFluxPressure BCs to ensure flux consistency setSnGrad ( - p_rgh.boundaryField(), + p_rgh.boundaryFieldRef(), ( phiHbyA.boundaryField() - ( diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C index 91b49b30e0..bb9683c39e 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C @@ -163,11 +163,13 @@ JohnsonJacksonSchaeffer::nu const fvPatchList& patches = phase.mesh().boundary(); const volVectorField& U = phase.U(); + volScalarField::GeometricBoundaryField& nufBf = nuf.boundaryFieldRef(); + forAll(patches, patchi) { if (!patches[patchi].coupled()) { - nuf.boundaryField()[patchi] = + nufBf[patchi] = ( pf.boundaryField()[patchi]*sin(phi_.value()) /( diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C index 8c072802d1..0853ec8450 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C @@ -152,11 +152,13 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::nu const fvPatchList& patches = phase.mesh().boundary(); const volVectorField& U = phase.U(); + volScalarField::GeometricBoundaryField& nufBf = nuf.boundaryFieldRef(); + forAll(patches, patchi) { if (!patches[patchi].coupled()) { - nuf.boundaryField()[patchi] = + nufBf[patchi] = ( pf.boundaryField()[patchi]*sin(phi_.value()) /( diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C index ff3be20ab1..d0a5b023d1 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C @@ -295,7 +295,7 @@ Foam::RASModels::kineticTheoryModel::pPrime() const ); volScalarField::GeometricBoundaryField& bpPrime = - tpPrime.ref().boundaryField(); + tpPrime.ref().boundaryFieldRef(); forAll(bpPrime, patchi) { diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C index 9d3b53c782..0edda3819d 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C @@ -165,7 +165,7 @@ Foam::RASModels::phasePressureModel::pPrime() const ); volScalarField::GeometricBoundaryField& bpPrime = - tpPrime.ref().boundaryField(); + tpPrime.ref().boundaryFieldRef(); forAll(bpPrime, patchi) { @@ -193,7 +193,7 @@ Foam::RASModels::phasePressureModel::pPrimef() const ); surfaceScalarField::GeometricBoundaryField& bpPrime = - tpPrime.ref().boundaryField(); + tpPrime.ref().boundaryFieldRef(); forAll(bpPrime, patchi) { diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C index 597d15e7ff..fc86ede74b 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C @@ -316,11 +316,13 @@ void Foam::twoPhaseSystem::solve() ) ); + surfaceScalarField::GeometricBoundaryField& alphaPhic1Bf = + alphaPhic1.boundaryFieldRef(); + // Ensure that the flux at inflow BCs is preserved - forAll(alphaPhic1.boundaryField(), patchi) + forAll(alphaPhic1Bf, patchi) { - fvsPatchScalarField& alphaPhic1p = - alphaPhic1.boundaryField()[patchi]; + fvsPatchScalarField& alphaPhic1p = alphaPhic1Bf[patchi]; if (!alphaPhic1p.coupled()) { diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H index b418877743..eb3eeb562a 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H @@ -146,17 +146,25 @@ while (pimple.correct()) surfaceScalarField phiCorrCoeff1(pos(alphaf1Bar - 0.99)); surfaceScalarField phiCorrCoeff2(pos(0.01 - alphaf1Bar)); - forAll(mesh.boundary(), patchi) { - // Set ddtPhiCorr to 0 on non-coupled boundaries - if - ( - !mesh.boundary()[patchi].coupled() - || isA(mesh.boundary()[patchi]) - ) + surfaceScalarField::GeometricBoundaryField& phiCorrCoeff1Bf = + phiCorrCoeff1.boundaryFieldRef(); + + surfaceScalarField::GeometricBoundaryField& phiCorrCoeff2Bf = + phiCorrCoeff2.boundaryFieldRef(); + + forAll(mesh.boundary(), patchi) { - phiCorrCoeff1.boundaryField()[patchi] = 0; - phiCorrCoeff2.boundaryField()[patchi] = 0; + // Set ddtPhiCorr to 0 on non-coupled boundaries + if + ( + !mesh.boundary()[patchi].coupled() + || isA(mesh.boundary()[patchi]) + ) + { + phiCorrCoeff1Bf[patchi] = 0; + phiCorrCoeff2Bf[patchi] = 0; + } } } @@ -212,7 +220,7 @@ while (pimple.correct()) // Update the fixedFluxPressure BCs to ensure flux consistency setSnGrad ( - p_rgh.boundaryField(), + p_rgh.boundaryFieldRef(), ( phiHbyA.boundaryField() - ( diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H index 6855b6fb96..bb78b2799d 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H @@ -200,7 +200,7 @@ while (pimple.correct()) // Update the fixedFluxPressure BCs to ensure flux consistency setSnGrad ( - p_rgh.boundaryField(), + p_rgh.boundaryFieldRef(), ( phiHbyA.boundaryField() - ( diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C index 6a5a474b34..8d915bc39e 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C @@ -279,7 +279,7 @@ Foam::RASModels::kineticTheoryModel::pPrime() const ); volScalarField::GeometricBoundaryField& bpPrime = - tpPrime.ref().boundaryField(); + tpPrime.ref().boundaryFieldRef(); forAll(bpPrime, patchi) { diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C index 0c5906d905..3245dec4db 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C @@ -171,7 +171,7 @@ Foam::RASModels::phasePressureModel::pPrime() const ); volScalarField::GeometricBoundaryField& bpPrime = - tpPrime.ref().boundaryField(); + tpPrime.ref().boundaryFieldRef(); forAll(bpPrime, patchi) { @@ -199,7 +199,7 @@ Foam::RASModels::phasePressureModel::pPrimef() const ); surfaceScalarField::GeometricBoundaryField& bpPrime = - tpPrime.ref().boundaryField(); + tpPrime.ref().boundaryFieldRef(); forAll(bpPrime, patchi) { diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C index ecf674e6e0..a3e187e7e3 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C @@ -36,6 +36,9 @@ void Foam::BlendedInterfacialModel::correctFixedFluxBCs GeometricField& field ) const { + typename GeometricField::GeometricBoundaryField& fieldBf = + field.boundaryFieldRef(); + forAll(pair_.phase1().phi().boundaryField(), patchi) { if @@ -46,7 +49,7 @@ void Foam::BlendedInterfacialModel::correctFixedFluxBCs ) ) { - field.boundaryField()[patchi] = Zero; + fieldBf[patchi] = Zero; } } } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C index 76af4a310a..a8e810f73c 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.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 @@ -446,11 +446,13 @@ void Foam::twoPhaseSystem::solve() ) ); + surfaceScalarField::GeometricBoundaryField& alphaPhic1Bf = + alphaPhic1.boundaryFieldRef(); + // Ensure that the flux at inflow BCs is preserved - forAll(alphaPhic1.boundaryField(), patchi) + forAll(alphaPhic1Bf, patchi) { - fvsPatchScalarField& alphaPhic1p = - alphaPhic1.boundaryField()[patchi]; + fvsPatchScalarField& alphaPhic1p = alphaPhic1Bf[patchi]; if (!alphaPhic1p.coupled()) { diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H index 2f1b71f7f1..1d2c3e825a 100644 --- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H +++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H @@ -429,30 +429,28 @@ public: // Member Functions - //- Return dimensioned internal field + //- Return a reference to the dimensioned internal field + // Note: this increments the event counter and checks the + // old-time fields; avoid in loops. DimensionedInternalField& dimensionedInternalField(); - //- Return dimensioned internal field + //- Return a const-reference to the dimensioned internal field inline const DimensionedInternalField& dimensionedInternalField() const; - //- Return internal field + //- Return a reference to the internal field + // Note: this increments the event counter and checks the + // old-time fields; avoid in loops. InternalField& internalField(); - //- Return internal field + //- Return a const-reference to the internal field inline const InternalField& internalField() const; - //- Return reference to GeometricBoundaryField + //- Return a reference to the boundary field + // Note: this increments the event counter and checks the + // old-time fields; avoid in loops. GeometricBoundaryField& boundaryFieldRef(); - //- Return reference to GeometricBoundaryField - #ifndef BOUNDARY_FIELD_REF - GeometricBoundaryField& boundaryField() - { - return boundaryFieldRef(); - } - #endif - - //- Return reference to GeometricBoundaryField for const field + //- Return const-reference to the boundary field inline const GeometricBoundaryField& boundaryField() const; //- Return the time index of the field diff --git a/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C b/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C index 32f4482247..01b55287f1 100644 --- a/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C +++ b/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C @@ -196,7 +196,7 @@ void mixtureKEpsilon::correctInletOutlet const volScalarField& refVsf ) const { - volScalarField::GeometricBoundaryField& bf = vsf.boundaryField(); + volScalarField::GeometricBoundaryField& bf = vsf.boundaryFieldRef(); const volScalarField::GeometricBoundaryField& refBf = refVsf.boundaryField(); diff --git a/src/finiteVolume/cfdTools/general/constrainPressure/constrainPressure.C b/src/finiteVolume/cfdTools/general/constrainPressure/constrainPressure.C index 760a335465..e9f1b40ecf 100644 --- a/src/finiteVolume/cfdTools/general/constrainPressure/constrainPressure.C +++ b/src/finiteVolume/cfdTools/general/constrainPressure/constrainPressure.C @@ -44,7 +44,7 @@ void Foam::constrainPressure { const fvMesh& mesh = p.mesh(); - volScalarField::GeometricBoundaryField& pBf = p.boundaryField(); + volScalarField::GeometricBoundaryField& pBf = p.boundaryFieldRef(); const volVectorField::GeometricBoundaryField& UBf = U.boundaryField(); const surfaceScalarField::GeometricBoundaryField& phiHbyABf = diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C index 6673e5d983..b18852a1ae 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C +++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C @@ -609,7 +609,7 @@ void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs) limitSum(phiPsiCorrsInternal); } - surfaceScalarField::GeometricBoundaryField& bfld = + const surfaceScalarField::GeometricBoundaryField& bfld = phiPsiCorrs[0].boundaryField(); forAll(bfld, patchi) @@ -622,7 +622,7 @@ void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs) phiPsiCorrsPatch.set ( phasei, - &phiPsiCorrs[phasei].boundaryField()[patchi] + &phiPsiCorrs[phasei].boundaryFieldRef()[patchi] ); }