GIT: Initial commit after latest foundation merge

.gitignore

@@ -58,8 +58,8 @@ doc/Doxygen/DTAGS
 # Untracked configuration files
 /etc/prefs.csh
 /etc/prefs.sh
-/etc/config/prefs.csh
+/etc/config.csh/prefs.csh
-/etc/config/prefs.sh
+/etc/config.sh/prefs.sh
 
 # Source packages - anywhere
 *.tar.bz2
@@ -83,5 +83,3 @@ doc/Doxygen/DTAGS
 
 # Ignore the test directory
 /tutorialsTest
-
-# end-of-file

Allwmake

@@ -43,4 +43,4 @@ then
 fi
 
 
-# ----------------------------------------------------------------- end-of-file
+#------------------------------------------------------------------------------

applications/Allwmake

@@ -24,4 +24,4 @@ wmake -all utilities $*
 wmake -all solvers $*
 
 
-# ----------------------------------------------------------------- end-of-file
+#------------------------------------------------------------------------------

applications/solvers/DNS/dnsFoam/dnsFoam.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
@@ -91,16 +91,17 @@ int main(int argc, char *argv[])
         {
             volScalarField rAU(1.0/UEqn.A());
             surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-            volVectorField HbyA("HbyA", U);
+            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-            HbyA = rAU*UEqn.H();
-
             surfaceScalarField phiHbyA
             (
                 "phiHbyA",
-                (fvc::interpolate(HbyA) & mesh.Sf())
+                fvc::flux(HbyA)
               + rAUf*fvc::ddtCorr(U, phi)
             );
 
+            // Update the pressure BCs to ensure flux consistency
+            constrainPressure(p, U, phiHbyA, rAUf);
+
             fvScalarMatrix pEqn
             (
                 fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)

applications/solvers/basic/potentialFoam/createFields.H

@@ -13,7 +13,7 @@ volVectorField U
 );
 
 // Initialise the velocity internal field to zero
-U = dimensionedVector("0", U.dimensions(), vector::zero);
+U = dimensionedVector("0", U.dimensions(), Zero);
 
 surfaceScalarField phi
 (
@@ -25,13 +25,13 @@ surfaceScalarField phi
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
-    fvc::interpolate(U) & mesh.Sf()
+    fvc::flux(U)
 );
 
 if (args.optionFound("initialiseUBCs"))
 {
     U.correctBoundaryConditions();
-    phi = fvc::interpolate(U) & mesh.Sf();
+    phi = fvc::flux(U);
 }
 
 

applications/solvers/basic/potentialFoam/potentialFoam.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
@@ -172,8 +172,7 @@ int main(int argc, char *argv[])
     U.correctBoundaryConditions();
 
     Info<< "Interpolated velocity error = "
-        << (sqrt(sum(sqr((fvc::interpolate(U) & mesh.Sf()) - phi)))
+        << (sqrt(sum(sqr(fvc::flux(U) - phi)))/sum(mesh.magSf())).value()
-          /sum(mesh.magSf())).value()
         << endl;
 
     // Write U and phi

applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C

@@ -91,10 +91,10 @@ int main(int argc, char *argv[])
                 fvOptions(T)
             );
 
+            TEqn.relax();
             fvOptions.constrain(TEqn);
-
             TEqn.solve();
-
+            fvOptions.correct(T);
         }
 
         runTime.write();

applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.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
@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "basicXiSubXiEq.H"
-#include "zeroGradientFvPatchFields.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -93,59 +92,43 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::basicSubGrid::XiEq() const
 
     const scalarField Cw = pow(mesh.V(), 2.0/3.0);
 
-    tmp<volScalarField> tN
+    volScalarField N
     (
-        new volScalarField
+        IOobject
         (
-            IOobject
+            "N",
-            (
+            mesh.time().constant(),
-                "tN",
-                mesh.time().constant(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
             mesh,
-            dimensionedScalar("zero", Nv.dimensions(), 0.0),
+            IOobject::NO_READ,
-            zeroGradientFvPatchVectorField::typeName
+            IOobject::NO_WRITE
-        )
+        ),
+        mesh,
+        dimensionedScalar("zero", Nv.dimensions(), 0.0)
     );
-
-    volScalarField& N = tN();
-
     N.internalField() = Nv.internalField()*Cw;
 
-    tmp<volSymmTensorField> tns
+    volSymmTensorField ns
     (
-        new volSymmTensorField
+        IOobject
         (
-            IOobject
+            "ns",
-            (
+            U.mesh().time().timeName(),
-                "tns",
-                U.mesh().time().timeName(),
-                U.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
             U.mesh(),
-            dimensionedSymmTensor
+            IOobject::NO_READ,
-            (
+            IOobject::NO_WRITE
-                "zero",
+        ),
-                nsv.dimensions(),
+        U.mesh(),
-                pTraits<symmTensor>::zero
+        dimensionedSymmTensor
-            ),
+        (
-             zeroGradientFvPatchSymmTensorField::typeName
+            "zero",
+            nsv.dimensions(),
+            Zero
         )
     );
-
-    volSymmTensorField& ns = tns();
-
     ns.internalField() = nsv.internalField()*Cw;
 
     volScalarField n(max(N - (Uhat & ns & Uhat), scalar(1e-4)));
-
     volScalarField b((Uhat & B_ & Uhat)/sqrt(n));
-
     volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
 
     volScalarField XiSubEq

applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.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
@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "basicXiSubG.H"
-#include "zeroGradientFvPatchFields.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -71,32 +70,10 @@ Foam::tmp<Foam::volScalarField> Foam::XiGModels::basicSubGrid::G() const
     const volScalarField& Lobs = db.lookupObject<volScalarField>("Lobs");
 
     tmp<volScalarField> tGtot = XiGModel_->G();
-    volScalarField& Gtot = tGtot();
+    volScalarField& Gtot = tGtot.ref();
 
     const scalarField Cw = pow(Su_.mesh().V(), 2.0/3.0);
-
+    scalarField N(Nv.internalField()*Cw);
-    tmp<volScalarField> tN
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "tN",
-                Su_.mesh().time().timeName(),
-                Su_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            Su_.mesh(),
-            dimensionedScalar("zero", Nv.dimensions(), 0.0),
-            zeroGradientFvPatchVectorField::typeName
-        )
-    );
-
-    volScalarField& N = tN();
-
-    N.internalField() = Nv.internalField()*Cw;
 
     forAll(N, celli)
     {

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.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
@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "basic.H"
-#include "zeroGradientFvPatchFields.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -109,13 +108,12 @@ Foam::tmp<Foam::volSymmTensorField> Foam::PDRDragModels::basic::Dcu() const
             (
                 "zero",
                 dimMass/dimTime/pow(dimLength, 3),
-                pTraits<symmTensor>::zero
+                Zero
-            ),
+            )
-            zeroGradientFvPatchSymmTensorField::typeName
         )
     );
 
-    volSymmTensorField& DragDcu = tDragDcu();
+    volSymmTensorField& DragDcu = tDragDcu.ref();
 
     if (on_)
     {
@@ -145,12 +143,11 @@ Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basic::Gk() const
                 IOobject::NO_WRITE
             ),
             U_.mesh(),
-            dimensionedScalar("zero", dimMass/dimLength/pow(dimTime, 3), 0.0),
+            dimensionedScalar("zero", dimMass/dimLength/pow(dimTime, 3), 0.0)
-            zeroGradientFvPatchVectorField::typeName
         )
     );
 
-    volScalarField& Gk = tGk();
+    volScalarField& Gk = tGk.ref();
 
     if (on_)
     {

applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.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
@@ -55,7 +55,7 @@ PDRkEpsilon::PDRkEpsilon
     const word& modelName
 )
 :
-    Foam::RASModels::kEpsilon<EddyDiffusivity<compressible::turbulenceModel> >
+    Foam::RASModels::kEpsilon<EddyDiffusivity<compressible::turbulenceModel>>
     (
         geometricOneField(),
         rho,
@@ -161,9 +161,9 @@ void PDRkEpsilon::correct()
       - fvm::Sp(C2_*betav*rho_*epsilon_/k_, epsilon_)
     );
 
-    epsEqn().relax();
+    epsEqn.ref().relax();
 
-    epsEqn().boundaryManipulate(epsilon_.boundaryField());
+    epsEqn.ref().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
     bound(epsilon_, epsilonMin_);
@@ -182,7 +182,7 @@ void PDRkEpsilon::correct()
       - fvm::Sp(betav*rho_*epsilon_/k_, k_)
     );
 
-    kEqn().relax();
+    kEqn.ref().relax();
     solve(kEqn);
     bound(k_, kMin_);
 

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.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
@@ -108,7 +108,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
             dimensionedScalar("XiEq", dimless, 0.0)
         )
     );
-    volScalarField& xieq = tXiEq();
+    volScalarField& xieq = tXiEq.ref();
 
     forAll(xieq, celli)
     {

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.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
@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "XiEqModel.H"
-#include "zeroGradientFvPatchFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -111,39 +110,30 @@ Foam::XiEqModel::calculateSchelkinEffect(const scalar uPrimeCoef) const
                 false
             ),
             mesh,
-            dimensionedScalar("zero", Nv.dimensions(), 0.0),
+            dimensionedScalar("zero", Nv.dimensions(), 0.0)
-            zeroGradientFvPatchVectorField::typeName
         )
     );
-
+    volScalarField& N = tN.ref();
-    volScalarField& N = tN();
-
     N.internalField() = Nv.internalField()*pow(mesh.V(), 2.0/3.0);
 
-    tmp<volSymmTensorField> tns
+    volSymmTensorField ns
     (
-        new volSymmTensorField
+        IOobject
         (
-            IOobject
+            "tns",
-            (
+            mesh.time().timeName(),
-                "tns",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
             mesh,
-            dimensionedSymmTensor
+            IOobject::NO_READ,
-            (
+            IOobject::NO_WRITE
-                "zero",
+        ),
-                nsv.dimensions(),
+        mesh,
-                pTraits<symmTensor>::zero
+        dimensionedSymmTensor
-            )
+        (
+            "zero",
+            nsv.dimensions(),
+            Zero
         )
     );
-
-    volSymmTensorField& ns = tns();
-
     ns.internalField() = nsv.internalField()*pow(mesh.V(), 2.0/3.0);
 
     const volVectorField Uhat

applications/solvers/combustion/PDRFoam/bEqn.H

@@ -1,4 +1,4 @@
-tmp<fv::convectionScheme<scalar> > mvConvection
+tmp<fv::convectionScheme<scalar>> mvConvection
 (
     fv::convectionScheme<scalar>::New
     (

applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.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
@@ -259,7 +259,7 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
         )
     );
 
-    volScalarField& Su0 = tSu0();
+    volScalarField& Su0 = tSu0.ref();
 
     forAll(Su0, celli)
     {
@@ -306,7 +306,7 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
         )
     );
 
-    volScalarField& Su0 = tSu0();
+    volScalarField& Su0 = tSu0.ref();
 
     forAll(Su0, celli)
     {
@@ -358,7 +358,7 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Ma
         )
     );
 
-    volScalarField& ma = tMa();
+    volScalarField& ma = tMa.ref();
 
     forAll(ma, celli)
     {

applications/solvers/combustion/PDRFoam/pEqn.H

@@ -1,9 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-
+volVectorField HbyA(constrainHbyA(invA & UEqn.H(), U, p));
-volVectorField HbyA("HbyA", U);
-HbyA = invA & UEqn.H();
 
 if (pimple.transonic())
 {
@@ -12,7 +10,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
           + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
            /fvc::interpolate(rho)
         )
@@ -43,7 +41,7 @@ else
     (
         "phiHbyA",
         (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
           + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
         )
     );

applications/solvers/combustion/XiFoam/ftEqn.H

@@ -1,4 +1,4 @@
-tmp<fv::convectionScheme<scalar> > mvConvection
+tmp<fv::convectionScheme<scalar>> mvConvection
 (
     fv::convectionScheme<scalar>::New
     (

applications/solvers/combustion/XiFoam/pEqn.H

@@ -2,9 +2,7 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -13,7 +11,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
           + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
         )
     );
@@ -45,13 +43,16 @@ else
     (
         "phiHbyA",
         (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
           + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
 
     MRF.makeRelative(phiHbyA);
 
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn

applications/solvers/combustion/chemFoam/createFields.H

@@ -80,8 +80,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedVector("zero", dimVelocity, vector::zero),
+        dimensionedVector("zero", dimVelocity, Zero)
-        p.boundaryField().types()
     );
 
     #include "createPhi.H"

applications/solvers/combustion/chemFoam/createSingleCellMesh.H

@@ -25,7 +25,7 @@ fvMesh mesh
         runTime,
         IOobject::READ_IF_PRESENT
     ),
-    xferMove<Field<vector> >(points),
+    xferMove<Field<vector>>(points),
     faces.xfer(),
     owner.xfer(),
     neighbour.xfer()

applications/solvers/combustion/engineFoam/pEqn.H

@@ -2,9 +2,7 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -14,7 +12,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (
-                (fvc::interpolate(HbyA) & mesh.Sf())
+                fvc::flux(HbyA)
               + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
             )
         )
@@ -29,7 +27,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
@@ -48,7 +46,7 @@ else
     (
         "phiHbyA",
         (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
           + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
         )
     );
@@ -56,13 +54,16 @@ else
     fvc::makeRelative(phiHbyA, rho, U);
     MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );

applications/solvers/combustion/fireFoam/YEEqn.H

@@ -1,4 +1,4 @@
-tmp<fv::convectionScheme<scalar> > mvConvection
+tmp<fv::convectionScheme<scalar>> mvConvection
 (
     fv::convectionScheme<scalar>::New
     (

applications/solvers/combustion/fireFoam/fireFoam.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
@@ -44,7 +44,6 @@ Description
 #include "psiCombustionModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/combustion/fireFoam/pEqn.H

@@ -2,9 +2,7 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn.H();
 
 surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
@@ -12,7 +10,7 @@ surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (
-        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+        fvc::flux(rho*HbyA)
       + rhorAUf*fvc::ddtCorr(rho, U, phi)
     )
   + phig
@@ -20,16 +18,8 @@ surfaceScalarField phiHbyA
 
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-// Update the fixedFluxPressure BCs to ensure flux consistency
+// Update the pressure BCs to ensure flux consistency
-setSnGrad<fixedFluxPressureFvPatchScalarField>
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-(
-    p_rgh.boundaryField(),
-    (
-        phiHbyA.boundaryField()
-      - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
-       *rho.boundaryField()
-    )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
-);
 
 while (pimple.correctNonOrthogonal())
 {

applications/solvers/combustion/reactingFoam/UEqn.H

@@ -2,7 +2,7 @@
 
 MRF.correctBoundaryVelocity(U);
 
-tmp<fvVectorMatrix> UEqn
+tmp<fvVectorMatrix> tUEqn
 (
     fvm::ddt(rho, U) + fvm::div(phi, U)
   + MRF.DDt(rho, U)
@@ -10,14 +10,15 @@ tmp<fvVectorMatrix> UEqn
  ==
     fvOptions(rho, U)
 );
+fvVectorMatrix& UEqn = tUEqn.ref();
 
-UEqn().relax();
+UEqn.relax();
 
-fvOptions.constrain(UEqn());
+fvOptions.constrain(UEqn);
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);

applications/solvers/combustion/reactingFoam/YEqn.H

@@ -1,4 +1,4 @@
-tmp<fv::convectionScheme<scalar> > mvConvection
+tmp<fv::convectionScheme<scalar>> mvConvection
 (
     fv::convectionScheme<scalar>::New
     (

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -3,15 +3,13 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 if (pimple.transonic())
@@ -21,7 +19,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
           + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
         )
     );
@@ -53,13 +51,16 @@ else
     (
         "phiHbyA",
         (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
           + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
 
     MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn

applications/solvers/combustion/reactingFoam/pcEqn.H

@@ -3,14 +3,13 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
-volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
+volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
-volVectorField HbyA("HbyA", U);
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 if (pimple.transonic())
@@ -20,7 +19,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
           + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
            /fvc::interpolate(rho)
         )
@@ -64,7 +63,7 @@ else
     (
         "phiHbyA",
         (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
           + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
         )
     );
@@ -76,6 +75,9 @@ else
 
     volScalarField rhorAtU("rhorAtU", rho*rAtU);
 
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H

@@ -7,9 +7,7 @@
 
     volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-    volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn.H();
 
     surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
@@ -17,7 +15,7 @@
     (
         "phiHbyA",
         (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
           + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
       + phig
@@ -25,16 +23,8 @@
 
     MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-    // Update the fixedFluxPressure BCs to ensure flux consistency
+    // Update the pressure BCs to ensure flux consistency
-    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
-           *rho.boundaryField()
-        )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
-    );
 
     fvScalarMatrix p_rghDDtEqn
     (

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,7 +39,6 @@ Description
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 

applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H

@@ -5,11 +5,9 @@
     // pressure solution - done in 2 parts. Part 1:
     thermo.rho() -= psi*p;
 
-    volScalarField rAU(1.0/UEqn().A());
+    volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-    volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn().H();
 
     if (pimple.transonic())
     {
@@ -17,7 +15,7 @@
         (
             "phiHbyA",
             (
-                (fvc::interpolate(HbyA) & mesh.Sf())
+                fvc::flux(HbyA)
               + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
             )
         );
@@ -58,13 +56,16 @@
         (
             "phiHbyA",
             (
-                (fvc::interpolate(rho*HbyA) & mesh.Sf())
+                fvc::flux(rho*HbyA)
               + rhorAUf*fvc::ddtCorr(rho, U, phi)
             )
         );
 
         MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+        // Update the pressure BCs to ensure flux consistency
+        constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
         fvScalarMatrix pDDtEqn
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p))

applications/solvers/combustion/reactingFoam/setRDeltaT.H

@@ -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
@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 {
-    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rDeltaT = trDeltaT.ref();
 
     const dictionary& pimpleDict = pimple.dict();
 

applications/solvers/compressible/rhoCentralFoam/Allwclean

@@ -6,4 +6,4 @@ wclean libso BCs
 wclean
 wclean rhoCentralDyMFoam
 
-# ----------------------------------------------------------------- end-of-file
+#------------------------------------------------------------------------------

applications/solvers/compressible/rhoCentralFoam/Allwmake

@@ -1,7 +1,11 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
+
+# Parse arguments for library compilation
+targetType=libso
+. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
 set -x
 
-(wmake libso BCs && wmake && wmake rhoCentralDyMFoam)
+(wmake $targetType BCs && wmake && wmake rhoCentralDyMFoam)
 
-# ----------------------------------------------------------------- end-of-file
+#------------------------------------------------------------------------------

applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,15 +26,10 @@ License
 #include "mixedFixedValueSlipFvPatchField.H"
 #include "symmTransformField.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>
-mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
 (
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF
@@ -47,7 +42,21 @@ mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
 
 
 template<class Type>
-mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    transformFvPatchField<Type>(p, iF),
+    refValue_("refValue", dict, p.size()),
+    valueFraction_("valueFraction", dict, p.size())
+{}
+
+
+template<class Type>
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
 (
     const mixedFixedValueSlipFvPatchField<Type>& ptf,
     const fvPatch& p,
@@ -62,21 +71,7 @@ mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
 
 
 template<class Type>
-mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
-(
-    const fvPatch& p,
-    const DimensionedField<Type, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    transformFvPatchField<Type>(p, iF),
-    refValue_("refValue", dict, p.size()),
-    valueFraction_("valueFraction", dict, p.size())
-{}
-
-
-template<class Type>
-mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
 (
     const mixedFixedValueSlipFvPatchField<Type>& ptf
 )
@@ -86,8 +81,9 @@ mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
     valueFraction_(ptf.valueFraction_)
 {}
 
+
 template<class Type>
-mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
 (
     const mixedFixedValueSlipFvPatchField<Type>& ptf,
     const DimensionedField<Type, volMesh>& iF
@@ -101,9 +97,8 @@ mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// Map from self
 template<class Type>
-void mixedFixedValueSlipFvPatchField<Type>::autoMap
+void Foam::mixedFixedValueSlipFvPatchField<Type>::autoMap
 (
     const fvPatchFieldMapper& m
 )
@@ -114,9 +109,8 @@ void mixedFixedValueSlipFvPatchField<Type>::autoMap
 }
 
 
-// Reverse-map the given fvPatchField onto this fvPatchField
 template<class Type>
-void mixedFixedValueSlipFvPatchField<Type>::rmap
+void Foam::mixedFixedValueSlipFvPatchField<Type>::rmap
 (
     const fvPatchField<Type>& ptf,
     const labelList& addr
@@ -125,16 +119,16 @@ void mixedFixedValueSlipFvPatchField<Type>::rmap
     transformFvPatchField<Type>::rmap(ptf, addr);
 
     const mixedFixedValueSlipFvPatchField<Type>& dmptf =
-        refCast<const mixedFixedValueSlipFvPatchField<Type> >(ptf);
+        refCast<const mixedFixedValueSlipFvPatchField<Type>>(ptf);
 
     refValue_.rmap(dmptf.refValue_, addr);
     valueFraction_.rmap(dmptf.valueFraction_, addr);
 }
 
 
-// Return gradient at boundary
 template<class Type>
-tmp<Field<Type> > mixedFixedValueSlipFvPatchField<Type>::snGrad() const
+Foam::tmp<Foam::Field<Type>>
+Foam::mixedFixedValueSlipFvPatchField<Type>::snGrad() const
 {
     tmp<vectorField> nHat = this->patch().nf();
     Field<Type> pif(this->patchInternalField());
@@ -147,9 +141,11 @@ tmp<Field<Type> > mixedFixedValueSlipFvPatchField<Type>::snGrad() const
 }
 
 
-// Evaluate the field on the patch
 template<class Type>
-void mixedFixedValueSlipFvPatchField<Type>::evaluate(const Pstream::commsTypes)
+void Foam::mixedFixedValueSlipFvPatchField<Type>::evaluate
+(
+    const Pstream::commsTypes
+)
 {
     if (!this->updated())
     {
@@ -170,10 +166,9 @@ void mixedFixedValueSlipFvPatchField<Type>::evaluate(const Pstream::commsTypes)
 }
 
 
-// Return defining fields
 template<class Type>
-tmp<Field<Type> >
+Foam::tmp<Foam::Field<Type>>
-mixedFixedValueSlipFvPatchField<Type>::snGradTransformDiag() const
+Foam::mixedFixedValueSlipFvPatchField<Type>::snGradTransformDiag() const
 {
     vectorField nHat(this->patch().nf());
     vectorField diag(nHat.size());
@@ -189,9 +184,8 @@ mixedFixedValueSlipFvPatchField<Type>::snGradTransformDiag() const
 }
 
 
-// Write
 template<class Type>
-void mixedFixedValueSlipFvPatchField<Type>::write(Ostream& os) const
+void Foam::mixedFixedValueSlipFvPatchField<Type>::write(Ostream& os) const
 {
     transformFvPatchField<Type>::write(os);
     refValue_.writeEntry("refValue", os);
@@ -199,8 +193,4 @@ void mixedFixedValueSlipFvPatchField<Type>::write(Ostream& os) const
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.H

@@ -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
@@ -101,9 +101,9 @@ public:
         );
 
         //- Construct and return a clone
-        virtual tmp<fvPatchField<Type> > clone() const
+        virtual tmp<fvPatchField<Type>> clone() const
         {
-            return tmp<fvPatchField<Type> >
+            return tmp<fvPatchField<Type>>
             (
                 new mixedFixedValueSlipFvPatchField<Type>(*this)
             );
@@ -117,12 +117,12 @@ public:
         );
 
         //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchField<Type> > clone
+        virtual tmp<fvPatchField<Type>> clone
         (
             const DimensionedField<Type, volMesh>& iF
         ) const
         {
-            return tmp<fvPatchField<Type> >
+            return tmp<fvPatchField<Type>>
             (
                 new mixedFixedValueSlipFvPatchField<Type>(*this, iF)
             );
@@ -130,6 +130,15 @@ public:
 
     // Member functions
 
+        // Access
+
+            //- Return false: this patch field is not altered by assignment
+            virtual bool assignable() const
+            {
+                return false;
+            }
+
+
         // Mapping functions
 
             //- Map (and resize as needed) from self given a mapping object
@@ -170,7 +179,7 @@ public:
         // Evaluation functions
 
             //- Return gradient at boundary
-            virtual tmp<Field<Type> > snGrad() const;
+            virtual tmp<Field<Type>> snGrad() const;
 
             //- Evaluate the patch field
             virtual void evaluate
@@ -179,7 +188,7 @@ public:
             );
 
             //- Return face-gradient transform diagonal
-            virtual tmp<Field<Type> > snGradTransformDiag() const;
+            virtual tmp<Field<Type>> snGradTransformDiag() const;
 
 
         //- Write

applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,7 +33,7 @@ Description
     This boundary condition provides a fixed density inlet condition for
     compressible solvers, where the density of calculated using:
 
-    \f{
+    \f[
         \rho = \psi p
     \f]
 

applications/solvers/compressible/rhoCentralFoam/createFields.H

@@ -95,7 +95,7 @@ surfaceScalarField neg
     dimensionedScalar("neg", dimless, -1.0)
 );
 
-surfaceScalarField phi("phi", mesh.Sf() & fvc::interpolate(rhoU));
+surfaceScalarField phi("phi", fvc::flux(rhoU));
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H

@@ -3,14 +3,14 @@ namespace Foam
 
 //- Interpolate field vf according to direction dir
 template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf,
     const surfaceScalarField& dir,
     const word& reconFieldName = word::null
 )
 {
-    tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsf
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tsf
     (
         fvc::interpolate
         (
@@ -22,7 +22,7 @@ tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate
         )
     );
 
-    GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf();
+    GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf.ref();
 
     sf.rename(vf.name() + '_' + dir.name());
 

applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.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
@@ -37,7 +37,6 @@ Description
 #include "dynamicFvMesh.H"
 #include "psiThermo.H"
 #include "turbulentFluidThermoModel.H"
-#include "zeroGradientFvPatchFields.H"
 #include "fixedRhoFvPatchScalarField.H"
 #include "directionInterpolate.H"
 #include "motionSolver.H"
@@ -212,7 +211,7 @@ int main(int argc, char *argv[])
             "sigmaDotU",
             (
                 fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
-              + (mesh.Sf() & fvc::interpolate(tauMC))
+              + fvc::dotInterpolate(mesh.Sf(), tauMC)
             )
             & (a_pos*U_pos + a_neg*U_neg)
         );

applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.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
@@ -36,7 +36,6 @@ Description
 #include "fvCFD.H"
 #include "psiThermo.H"
 #include "turbulentFluidThermoModel.H"
-#include "zeroGradientFvPatchFields.H"
 #include "fixedRhoFvPatchScalarField.H"
 #include "directionInterpolate.H"
 #include "localEulerDdtScheme.H"
@@ -205,9 +204,9 @@ int main(int argc, char *argv[])
             "sigmaDotU",
             (
                 fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
-              + (mesh.Sf() & fvc::interpolate(tauMC))
+              + fvc::dotInterpolate(mesh.Sf(), tauMC)
             )
-            & (a_pos*U_pos + a_neg*U_neg)
+          & (a_pos*U_pos + a_neg*U_neg)
         );
 
         solve

applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H

@@ -1,5 +1,5 @@
 {
-    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rDeltaT = trDeltaT.ref();
 
     scalar rDeltaTSmoothingCoeff
     (

applications/solvers/compressible/rhoPimpleFoam/UEqn.H

@@ -2,7 +2,7 @@
 
 MRF.correctBoundaryVelocity(U);
 
-tmp<fvVectorMatrix> UEqn
+tmp<fvVectorMatrix> tUEqn
 (
     fvm::ddt(rho, U) + fvm::div(phi, U)
   + MRF.DDt(rho, U)
@@ -10,14 +10,15 @@ tmp<fvVectorMatrix> UEqn
  ==
     fvOptions(rho, U)
 );
+fvVectorMatrix& UEqn = tUEqn.ref();
 
-UEqn().relax();
+UEqn.relax();
 
-fvOptions.constrain(UEqn());
+fvOptions.constrain(UEqn);
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -3,15 +3,13 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 if (pimple.transonic())
@@ -21,7 +19,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
           + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
         )
     );
@@ -53,13 +51,16 @@ else
     (
         "phiHbyA",
         (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
           + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
 
     MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn

applications/solvers/compressible/rhoPimpleFoam/pcEqn.H

@@ -3,14 +3,13 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
-volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
+volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
-volVectorField HbyA("HbyA", U);
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 if (pimple.transonic())
@@ -20,7 +19,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
           + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
            /fvc::interpolate(rho)
         )
@@ -64,7 +63,7 @@ else
     (
         "phiHbyA",
         (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
           + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
         )
     );
@@ -76,6 +75,9 @@ else
 
     volScalarField rhorAtU("rhorAtU", rho*rAtU);
 
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H

@@ -3,15 +3,13 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 if (pimple.transonic())
@@ -21,7 +19,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
           + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
         )
     );
@@ -53,13 +51,16 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+        fvc::flux(rho*HbyA)
       + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
     );
 
     fvc::makeRelative(phiHbyA, rho, U);
     MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
     while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector

applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H

@@ -1,5 +1,5 @@
 {
-    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rDeltaT = trDeltaT.ref();
 
     const dictionary& pimpleDict = pimple.dict();
 
@@ -38,7 +38,7 @@
         surfaceScalarField phid
         (
             "phid",
-            fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
+            fvc::interpolate(psi)*fvc::flux(U)
         );
 
         rDeltaT.dimensionedInternalField() = max

applications/solvers/compressible/rhoSimpleFoam/UEqn.H

@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(rho, U)
@@ -10,11 +10,12 @@
      ==
         fvOptions(rho, U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
-    fvOptions.constrain(UEqn());
+    fvOptions.constrain(UEqn);
 
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -1,9 +1,8 @@
 {
-    volScalarField rAU(1.0/UEqn().A());
+    volScalarField rAU(1.0/UEqn.A());
-    volVectorField HbyA("HbyA", U);
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    HbyA = rAU*UEqn().H();
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-
+    tUEqn.clear();
-    UEqn.clear();
 
     bool closedVolume = false;
 
@@ -13,7 +12,7 @@
         (
             "phid",
             fvc::interpolate(psi)
-           *(fvc::interpolate(HbyA) & mesh.Sf())
+           *fvc::flux(HbyA)
         );
 
         MRF.makeRelative(fvc::interpolate(psi), phid);
@@ -23,7 +22,7 @@
             fvScalarMatrix pEqn
             (
                 fvm::div(phid, p)
-              - fvm::laplacian(rho*rAU, p)
+              - fvm::laplacian(rhorAUf, p)
               ==
                 fvOptions(psi, p, rho.name())
             );
@@ -43,22 +42,20 @@
     }
     else
     {
-        surfaceScalarField phiHbyA
+        surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-        (
-            "phiHbyA",
-            fvc::interpolate(rho*HbyA) & mesh.Sf()
-        );
-
         MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
         closedVolume = adjustPhi(phiHbyA, U, p);
 
+        // Update the pressure BCs to ensure flux consistency
+        constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
         while (simple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
                 fvc::div(phiHbyA)
-              - fvm::laplacian(rho*rAU, p)
+              - fvm::laplacian(rhorAUf, p)
               ==
                 fvOptions(psi, p, rho.name())
             );

applications/solvers/compressible/rhoSimpleFoam/pcEqn.H

@@ -1,10 +1,7 @@
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
-volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
+volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
-
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-volVectorField HbyA("HbyA", U);
+tUEqn.clear();
-HbyA = rAU*UEqn().H();
-
-UEqn.clear();
 
 bool closedVolume = false;
 
@@ -14,7 +11,7 @@ if (simple.transonic())
     (
         "phid",
         fvc::interpolate(psi)
-       *(fvc::interpolate(HbyA) & mesh.Sf())
+       *fvc::flux(HbyA)
     );
 
     MRF.makeRelative(fvc::interpolate(psi), phid);
@@ -55,12 +52,7 @@ if (simple.transonic())
 }
 else
 {
-    surfaceScalarField phiHbyA
+    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-    (
-        "phiHbyA",
-        fvc::interpolate(rho*HbyA) & mesh.Sf()
-    );
-
     MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     closedVolume = adjustPhi(phiHbyA, U, p);
@@ -70,6 +62,9 @@ else
 
     volScalarField rhorAtU("rhorAtU", rho*rAtU);
 
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
+
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/UEqn.H

@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(rho, U)
@@ -10,8 +10,9 @@
      ==
         fvOptions(rho, U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
     // Include the porous media resistance and solve the momentum equation
     // either implicit in the tensorial resistance or transport using by
@@ -22,18 +23,18 @@
 
     if (pressureImplicitPorosity)
     {
-        tmp<volTensorField> tTU = tensor(I)*UEqn().A();
+        tmp<volTensorField> tTU = tensor(I)*UEqn.A();
-        pZones.addResistance(UEqn(), tTU());
+        pZones.addResistance(UEqn, tTU.ref());
         trTU = inv(tTU());
-        trTU().rename("rAU");
+        trTU.ref().rename("rAU");
 
-        fvOptions.constrain(UEqn());
+        fvOptions.constrain(UEqn);
 
         volVectorField gradp(fvc::grad(p));
 
         for (int UCorr=0; UCorr<nUCorr; UCorr++)
         {
-            U = trTU() & (UEqn().H() - gradp);
+            U = trTU() & (UEqn.H() - gradp);
         }
         U.correctBoundaryConditions();
 
@@ -41,14 +42,14 @@
     }
     else
     {
-        pZones.addResistance(UEqn());
+        pZones.addResistance(UEqn);
 
-        fvOptions.constrain(UEqn());
+        fvOptions.constrain(UEqn);
 
-        solve(UEqn() == -fvc::grad(p));
+        solve(UEqn == -fvc::grad(p));
 
         fvOptions.correct(U);
 
-        trAU = 1.0/UEqn().A();
+        trAU = 1.0/UEqn.A();
-        trAU().rename("rAU");
+        trAU.ref().rename("rAU");
     }

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H

@@ -68,9 +68,9 @@ dimensionedScalar rhoMin
 );
 
 Info<< "Creating turbulence model\n" << endl;
-autoPtr<compressible::RASModel> turbulence
+autoPtr<compressible::turbulenceModel> turbulence
 (
-    compressible::New<compressible::RASModel>
+    compressible::turbulenceModel::New
     (
         rho,
         U,

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H

@@ -1,27 +1,22 @@
 {
     const volScalarField& psi = thermo.psi();
 
-    volVectorField HbyA("HbyA", U);
+    tmp<volVectorField> tHbyA;
-
     if (pressureImplicitPorosity)
     {
-        HbyA = trTU() & UEqn().H();
+        tHbyA = constrainHbyA(trTU()&UEqn.H(), U, p);
     }
     else
     {
-        HbyA = trAU()*UEqn().H();
+        tHbyA = constrainHbyA(trAU()*UEqn.H(), U, p);
     }
+    volVectorField& HbyA = tHbyA.ref();
 
-    UEqn.clear();
+    tUEqn.clear();
 
     bool closedVolume = false;
 
-    surfaceScalarField phiHbyA
+    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-    (
-        "phiHbyA",
-        fvc::interpolate(rho*HbyA) & mesh.Sf()
-    );
-
     MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     closedVolume = adjustPhi(phiHbyA, U, p);
@@ -51,13 +46,15 @@
             );
         }
 
-        tpEqn().setReference(pRefCell, pRefValue);
+        fvScalarMatrix& pEqn = tpEqn.ref();
 
-        tpEqn().solve();
+        pEqn.setReference(pRefCell, pRefValue);
+
+        pEqn.solve();
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi = phiHbyA - tpEqn().flux();
+            phi = phiHbyA - pEqn.flux();
         }
     }
 

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -2,16 +2,13 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn.H();
-
 surfaceScalarField phid
 (
     "phid",
     fvc::interpolate(psi)
    *(
-        (mesh.Sf() & fvc::interpolate(HbyA))
+        fvc::flux(HbyA)
       + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
     )
 );

applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H

@@ -2,16 +2,13 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn.H();
-
 surfaceScalarField phid
 (
     "phid",
     fvc::interpolate(psi)
    *(
-       (mesh.Sf() & fvc::interpolate(HbyA))
+       fvc::flux(HbyA)
      + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
     )
 );

applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.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
@@ -92,7 +92,7 @@ int main(int argc, char *argv[])
                     "phid",
                     psi
                    *(
-                       (fvc::interpolate(U) & mesh.Sf())
+                       fvc::flux(U)
                      + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
                     )
                 );

applications/solvers/electromagnetics/magneticFoam/magnet.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,7 +72,7 @@ public:
         inline magnet()
         :
             remanence_("Mr", dimensionSet(0, -1, 0, 0, 0, 1, 0), 0),
-            orientation_(vector::zero)
+            orientation_(Zero)
         {}
 
         //- Construct from components
@@ -94,7 +94,7 @@ public:
         inline magnet(Istream& is)
         :
             remanence_("Mr", dimensionSet(0, -1, 0, 0, 0, 1, 0), 0),
-            orientation_(vector::zero)
+            orientation_(Zero)
         {
             is >> *this;
         }

applications/solvers/electromagnetics/mhdFoam/createPhiB.H

@@ -40,7 +40,7 @@
                 IOobject::NO_READ,
                 IOobject::AUTO_WRITE
             ),
-            (fvc::interpolate(B) & mesh.Sf())
+            fvc::flux(B)
         );
     }
 

applications/solvers/electromagnetics/mhdFoam/mhdFoam.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
@@ -100,17 +100,17 @@ int main(int argc, char *argv[])
             {
                 volScalarField rAU(1.0/UEqn.A());
                 surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-
+                volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-                volVectorField HbyA("HbyA", U);
-                HbyA = rAU*UEqn.H();
-
                 surfaceScalarField phiHbyA
                 (
                     "phiHbyA",
-                    (fvc::interpolate(HbyA) & mesh.Sf())
+                    fvc::flux(HbyA)
                   + rAUf*fvc::ddtCorr(U, phi)
                 );
 
+                // Update the pressure BCs to ensure flux consistency
+                constrainPressure(p, U, phiHbyA, rAUf);
+
                 while (piso.correctNonOrthogonal())
                 {
                     fvScalarMatrix pEqn
@@ -150,8 +150,7 @@ int main(int argc, char *argv[])
             volScalarField rAB(1.0/BEqn.A());
             surfaceScalarField rABf("rABf", fvc::interpolate(rAB));
 
-            phiB = (fvc::interpolate(B) & mesh.Sf())
+            phiB = fvc::flux(B) + rABf*fvc::ddtCorr(B, phiB);
-                + rABf*fvc::ddtCorr(B, phiB);
 
             while (bpiso.correctNonOrthogonal())
             {

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.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
@@ -54,7 +54,6 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "pimpleControl.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H

@@ -1,31 +1,22 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-    volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn.H();
 
     surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
       + rAUf*fvc::ddtCorr(U, phi)
       + phig
     );
 
     MRF.makeRelative(phiHbyA);
 
-    // Update the fixedFluxPressure BCs to ensure flux consistency
+    // Update the pressure BCs to ensure flux consistency
-    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
-        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-    );
 
     while (pimple.correctNonOrthogonal())
     {

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H

@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(U)
@@ -10,16 +10,17 @@
      ==
         fvOptions(U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
-    fvOptions.constrain(UEqn());
+    fvOptions.constrain(UEqn);
 
     if (simple.momentumPredictor())
     {
         solve
         (
-            UEqn()
+            UEqn
           ==
             fvc::reconstruct
             (

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.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
@@ -54,7 +54,6 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "simpleControl.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H

@@ -1,17 +1,16 @@
 {
-    volScalarField rAU("rAU", 1.0/UEqn().A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
 
-    volVectorField HbyA("HbyA", U);
+    tUEqn.clear();
-    HbyA = rAU*UEqn().H();
-    UEqn.clear();
 
     surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
     );
 
     MRF.makeRelative(phiHbyA);
@@ -20,15 +19,8 @@
 
     phiHbyA += phig;
 
-    // Update the fixedFluxPressure BCs to ensure flux consistency
+    // Update the pressure BCs to ensure flux consistency
-    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
-        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-    );
 
     while (simple.correctNonOrthogonal())
     {

applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.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
@@ -42,7 +42,6 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "pimpleControl.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -11,35 +11,25 @@
     thermo.rho() -= psi*p_rgh;
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf("rAUf", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
 
-    volVectorField HbyA("HbyA", U);
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-    HbyA = rAU*UEqn.H();
-
-    surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
-          + rAUf*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
       + phig
     );
 
     MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-    // Update the fixedFluxPressure BCs to ensure flux consistency
+    // Update the pressure BCs to ensure flux consistency
-    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
-           *rho.boundaryField()
-        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-    );
 
     tmp<fvScalarMatrix> p_rghDDtEqn
     (
@@ -51,6 +41,7 @@
         p_rghDDtEqn =
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+          + fvc::div(phiHbyA)
         ==
             fvOptions(psi, p_rgh, rho.name())
         );
@@ -61,8 +52,7 @@
         fvScalarMatrix p_rghEqn
         (
             p_rghDDtEqn()
-          + fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAUf, p_rgh)
-          - fvm::laplacian(rAUf, p_rgh)
         );
 
         p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
@@ -77,7 +67,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
             K = 0.5*magSqr(U);

applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H

@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(rho, U)
@@ -10,16 +10,17 @@
      ==
         fvOptions(rho, U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
-    fvOptions.constrain(UEqn());
+    fvOptions.constrain(UEqn);
 
     if (simple.momentumPredictor())
     {
         solve
         (
-            UEqn()
+            UEqn
          ==
             fvc::reconstruct
             (

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.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
@@ -39,7 +39,6 @@ Description
 #include "radiationModel.H"
 #include "simpleControl.H"
 #include "fvOptions.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -2,19 +2,17 @@
     rho = thermo.rho();
     rho.relax();
 
-    volScalarField rAU("rAU", 1.0/UEqn().A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-    volVectorField HbyA("HbyA", U);
+    tUEqn.clear();
-    HbyA = rAU*UEqn().H();
-    UEqn.clear();
 
     surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+        fvc::flux(rho*HbyA)
     );
 
     MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
@@ -23,16 +21,8 @@
 
     phiHbyA += phig;
 
-    // Update the fixedFluxPressure BCs to ensure flux consistency
+    // Update the pressure BCs to ensure flux consistency
-    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
-           *rho.boundaryField()
-        )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
-    );
 
 
     while (simple.correctNonOrthogonal())

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.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
@@ -49,7 +49,6 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "coordinateSystem.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.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
@@ -41,7 +41,6 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "coordinateSystem.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H

@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(rho, U)
@@ -10,14 +10,15 @@
      ==
         fvOptions(rho, U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
-    fvOptions.constrain(UEqn());
+    fvOptions.constrain(UEqn);
 
     solve
     (
-        UEqn()
+        UEqn
      ==
         fvc::reconstruct
         (

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H

@@ -4,19 +4,17 @@
     rho = min(rho, rhoMax[i]);
     rho.relax();
 
-    volScalarField rAU("rAU", 1.0/UEqn().A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-    volVectorField HbyA("HbyA", U);
+    tUEqn.clear();
-    HbyA = rAU*UEqn().H();
-    UEqn.clear();
 
     surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+        fvc::flux(rho*HbyA)
     );
 
     MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
@@ -25,16 +23,8 @@
 
     phiHbyA += phig;
 
-    // Update the fixedFluxPressure BCs to ensure flux consistency
+    // Update the pressure BCs to ensure flux consistency
-    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
-           *rho.boundaryField()
-        )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
-    );
 
     dimensionedScalar compressibility = fvc::domainIntegrate(psi);
     bool compressible = (compressibility.value() > SMALL);

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H

@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::ddt(rho, U) + fvm::div(phi, U)
       + MRF.DDt(rho, U)
@@ -10,16 +10,17 @@
      ==
         fvOptions(rho, U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
-    fvOptions.constrain(UEqn());
+    fvOptions.constrain(UEqn);
 
     if (momentumPredictor)
     {
         solve
         (
-            UEqn()
+            UEqn
           ==
             fvc::reconstruct
             (

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H

@@ -15,7 +15,7 @@
         );
     }
 /*
-    forAll (porousFluidRegions, porousI)
+    forAll(porousFluidRegions, porousI)
     {
         CoNum = max
         (

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

@@ -5,11 +5,9 @@
 
     rho = thermo.rho();
 
-    volScalarField rAU("rAU", 1.0/UEqn().A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-    volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn().H();
 
     surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
@@ -17,7 +15,7 @@
     (
         "phiHbyA",
         (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
           + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
       + phig
@@ -25,16 +23,8 @@
 
     MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-    // Update the fixedFluxPressure BCs to ensure flux consistency
+    // Update the pressure BCs to ensure flux consistency
-    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
-           *rho.boundaryField()
-        )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
-    );
 
     tmp<fvScalarMatrix> p_rghDDtEqn
     (

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H

@@ -55,7 +55,7 @@
                     (
                         "zero",
                         tkappaByCp().dimensions(),
-                        symmTensor::zero
+                        Zero
                     ),
                     zeroGradientFvPatchSymmTensorField::typeName
                 )

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H

@@ -6,7 +6,7 @@ if (finalIter)
 {
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        tmp<fvScalarMatrix> hEqn
+        fvScalarMatrix hEqn
         (
             fvm::ddt(betav*rho, h)
           - (
@@ -18,11 +18,11 @@ if (finalIter)
             fvOptions(rho, h)
         );
 
-        hEqn().relax();
+        hEqn.relax();
 
-        fvOptions.constrain(hEqn());
+        fvOptions.constrain(hEqn);
 
-        hEqn().solve(mesh.solver(h.select(finalIter)));
+        hEqn.solve(mesh.solver(h.select(finalIter)));
 
         fvOptions.correct(h);
     }

applications/solvers/heatTransfer/thermoFoam/setAlphaEff.H

@@ -1,12 +1,14 @@
     Info<< "Creating turbulence model\n" << endl;
     tmp<volScalarField> talphaEff;
 
-    IOobject turbulenceHeader
+    IOobject turbulencePropertiesHeader
     (
         "turbulenceProperties",
         runTime.constant(),
         mesh,
-        IOobject::MUST_READ
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE,
+        false
     );
 
     IOobject RASHeader

applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.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
@@ -63,7 +63,7 @@ void zeroCells
 {
     forAll(cells, i)
     {
-        vf[cells[i]] = pTraits<Type>::zero;
+        vf[cells[i]] = Zero;
     }
 }
 
@@ -110,7 +110,7 @@ int main(int argc, char *argv[])
         {
             // Momentum predictor
 
-            tmp<fvVectorMatrix> UEqn
+            tmp<fvVectorMatrix> tUEqn
             (
                 fvm::div(phi, U)
               + turbulence->divDevReff(U)
@@ -118,26 +118,25 @@ int main(int argc, char *argv[])
              ==
                 fvOptions(U)
             );
+            fvVectorMatrix& UEqn = tUEqn.ref();
 
-            UEqn().relax();
+            UEqn.relax();
 
-            fvOptions.constrain(UEqn());
+            fvOptions.constrain(UEqn);
 
-            solve(UEqn() == -fvc::grad(p));
+            solve(UEqn == -fvc::grad(p));
 
             fvOptions.correct(U);
 
-            volScalarField rAU(1.0/UEqn().A());
+            volScalarField rAU(1.0/UEqn.A());
-            volVectorField HbyA("HbyA", U);
+            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-            HbyA = rAU*UEqn().H();
+            tUEqn.clear();
-            UEqn.clear();
+            surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
-            surfaceScalarField phiHbyA
-            (
-                "phiHbyA",
-                fvc::interpolate(HbyA) & mesh.Sf()
-            );
             adjustPhi(phiHbyA, U, p);
 
+            // Update the pressure BCs to ensure flux consistency
+            constrainPressure(p, U, phiHbyA, rAU);
+
             // Non-orthogonal pressure corrector loop
             while (simple.correctNonOrthogonal())
             {
@@ -175,13 +174,13 @@ int main(int argc, char *argv[])
             //(
             //    fvc::reconstruct
             //    (
-            //        mesh.magSf()*(fvc::snGrad(Ua) & fvc::interpolate(U))
+            //        mesh.magSf()*fvc::dotInterpolate(fvc::snGrad(Ua), U)
             //    )
             //);
 
             zeroCells(adjointTransposeConvection, inletCells);
 
-            tmp<fvVectorMatrix> UaEqn
+            tmp<fvVectorMatrix> tUaEqn
             (
                 fvm::div(-phi, Ua)
               - adjointTransposeConvection
@@ -190,24 +189,21 @@ int main(int argc, char *argv[])
              ==
                 fvOptions(Ua)
             );
+            fvVectorMatrix& UaEqn = tUaEqn.ref();
 
-            UaEqn().relax();
+            UaEqn.relax();
 
-            fvOptions.constrain(UaEqn());
+            fvOptions.constrain(UaEqn);
 
-            solve(UaEqn() == -fvc::grad(pa));
+            solve(UaEqn == -fvc::grad(pa));
 
             fvOptions.correct(Ua);
 
-            volScalarField rAUa(1.0/UaEqn().A());
+            volScalarField rAUa(1.0/UaEqn.A());
             volVectorField HbyAa("HbyAa", Ua);
-            HbyAa = rAUa*UaEqn().H();
+            HbyAa = rAUa*UaEqn.H();
-            UaEqn.clear();
+            tUaEqn.clear();
-            surfaceScalarField phiHbyAa
+            surfaceScalarField phiHbyAa("phiHbyAa", fvc::flux(HbyAa));
-            (
-                "phiHbyAa",
-                fvc::interpolate(HbyAa) & mesh.Sf()
-            );
             adjustPhi(phiHbyAa, Ua, pa);
 
             // Non-orthogonal pressure corrector loop

applications/solvers/incompressible/adjointShapeOptimizationFoam/createFields.H

@@ -80,9 +80,9 @@ mesh.setFluxRequired(pa.name());
 
 singlePhaseTransportModel laminarTransport(U, phi);
 
-autoPtr<incompressible::RASModel> turbulence
+autoPtr<incompressible::turbulenceModel> turbulence
 (
-    incompressible::New<incompressible::RASModel>(U, phi, laminarTransport)
+    incompressible::turbulenceModel::New(U, phi, laminarTransport)
 );
 
 

applications/solvers/incompressible/boundaryFoam/createFields.H

@@ -31,9 +31,9 @@ surfaceScalarField phi
 
 singlePhaseTransportModel laminarTransport(U, phi);
 
-autoPtr<incompressible::RASModel> turbulence
+autoPtr<incompressible::turbulenceModel> turbulence
 (
-    incompressible::New<incompressible::RASModel>(U, phi, laminarTransport)
+    incompressible::turbulenceModel::New(U, phi, laminarTransport)
 );
 
 dimensionedVector Ubar("Ubar", dimVelocity, laminarTransport);
@@ -45,5 +45,5 @@ dimensionedVector gradP
 (
     "gradP",
     dimensionSet(0, 1, -2, 0, 0),
-    vector::zero
+    Zero
 );

applications/solvers/incompressible/boundaryFoam/interrogateWallPatches.H

@@ -3,7 +3,7 @@
 label faceId(-1);
 label patchId(-1);
 label nWallFaces(0);
-vector wallNormal(vector::zero);
+vector wallNormal(Zero);
 
 const fvPatchList& patches = mesh.boundary();
 

applications/solvers/incompressible/icoFoam/icoFoam.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
@@ -104,18 +104,19 @@ int main(int argc, char *argv[])
         while (piso.correct())
         {
             volScalarField rAU(1.0/UEqn.A());
-
+            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-            volVectorField HbyA("HbyA", U);
-            HbyA = rAU*UEqn.H();
             surfaceScalarField phiHbyA
             (
                 "phiHbyA",
-                (fvc::interpolate(HbyA) & mesh.Sf())
+                fvc::flux(HbyA)
               + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
             );
 
             adjustPhi(phiHbyA, U, p);
 
+            // Update the pressure BCs to ensure flux consistency
+            constrainPressure(p, U, phiHbyA, rAU);
+
             // Non-orthogonal pressure corrector loop
             while (piso.correctNonOrthogonal())
             {

applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.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
@@ -80,18 +80,19 @@ int main(int argc, char *argv[])
         while (piso.correct())
         {
             volScalarField rAU(1.0/UEqn.A());
-
+            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-            volVectorField HbyA("HbyA", U);
-            HbyA = rAU*UEqn.H();
             surfaceScalarField phiHbyA
             (
                 "phiHbyA",
-                (fvc::interpolate(HbyA) & mesh.Sf())
+                fvc::flux(HbyA)
               + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
             );
 
             adjustPhi(phiHbyA, U, p);
 
+            // Update the pressure BCs to ensure flux consistency
+            constrainPressure(p, U, phiHbyA, rAU);
+
             // Non-orthogonal pressure corrector loop
             while (piso.correctNonOrthogonal())
             {

applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/UrelEqn.H

@@ -1,5 +1,5 @@
     // Relative momentum predictor
-    tmp<fvVectorMatrix> UrelEqn
+    tmp<fvVectorMatrix> tUrelEqn
     (
         fvm::ddt(Urel)
       + fvm::div(phi, Urel)
@@ -8,11 +8,12 @@
      ==
         fvOptions(Urel)
     );
+    fvVectorMatrix& UrelEqn = tUrelEqn.ref();
 
-    UrelEqn().relax();
+    UrelEqn.relax();
 
-    fvOptions.constrain(UrelEqn());
+    fvOptions.constrain(UrelEqn);
 
-    solve(UrelEqn() == -fvc::grad(p));
+    solve(UrelEqn == -fvc::grad(p));
 
     fvOptions.correct(Urel);

applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H

@@ -1,11 +1,11 @@
-volScalarField rAUrel(1.0/UrelEqn().A());
+volScalarField rAUrel(1.0/UrelEqn.A());
 volVectorField HbyA("HbyA", Urel);
-HbyA = rAUrel*UrelEqn().H();
+HbyA = rAUrel*UrelEqn.H();
 
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
-    (fvc::interpolate(HbyA) & mesh.Sf())
+    fvc::flux(HbyA)
   + fvc::interpolate(rAUrel)*fvc::ddtCorr(Urel, phi)
 );
 
@@ -15,7 +15,7 @@ tmp<volScalarField> rAtUrel(rAUrel);
 
 if (pimple.consistent())
 {
-    rAtUrel = 1.0/max(1.0/rAUrel - UrelEqn().H1(), 0.1/rAUrel);
+    rAtUrel = 1.0/max(1.0/rAUrel - UrelEqn.H1(), 0.1/rAUrel);
     phiHbyA +=
         fvc::interpolate(rAtUrel() - rAUrel)*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAUrel - rAtUrel())*fvc::grad(p);
@@ -23,9 +23,12 @@ if (pimple.consistent())
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UrelEqn.clear();
+    tUrelEqn.clear();
 }
 
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, Urel, phiHbyA, rAtUrel());
+
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
 {

applications/solvers/incompressible/pimpleFoam/UEqn.H

@@ -2,7 +2,7 @@
 
 MRF.correctBoundaryVelocity(U);
 
-tmp<fvVectorMatrix> UEqn
+tmp<fvVectorMatrix> tUEqn
 (
     fvm::ddt(U) + fvm::div(phi, U)
   + MRF.DDt(U)
@@ -10,14 +10,15 @@ tmp<fvVectorMatrix> UEqn
  ==
     fvOptions(U)
 );
+fvVectorMatrix& UEqn = tUEqn.ref();
 
-UEqn().relax();
+UEqn.relax();
 
-fvOptions.constrain(UEqn());
+fvOptions.constrain(UEqn);
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);
 }

applications/solvers/incompressible/pimpleFoam/pEqn.H

@@ -1,11 +1,9 @@
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
-volVectorField HbyA("HbyA", U);
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-HbyA = rAU*UEqn().H();
-
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
-    (fvc::interpolate(HbyA) & mesh.Sf())
+    fvc::flux(HbyA)
   + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
 );
 
@@ -17,7 +15,7 @@ tmp<volScalarField> rAtU(rAU);
 
 if (pimple.consistent())
 {
-    rAtU = 1.0/max(1.0/rAU - UEqn().H1(), 0.1/rAU);
+    rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
     phiHbyA +=
         fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU())*fvc::grad(p);
@@ -25,20 +23,11 @@ if (pimple.consistent())
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
-surfaceScalarField rAUf("rAUf", fvc::interpolate(rAtU()));
+// Update the pressure BCs to ensure flux consistency
-
+constrainPressure(p, U, phiHbyA, rAtU(), MRF);
-// Update the fixedFluxPressure BCs to ensure flux consistency
-setSnGrad<fixedFluxPressureFvPatchScalarField>
-(
-    p.boundaryField(),
-    (
-        phiHbyA.boundaryField()
-      - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
-    )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-);
 
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
@@ -46,7 +35,7 @@ while (pimple.correctNonOrthogonal())
     // Pressure corrector
     fvScalarMatrix pEqn
     (
-        fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
+        fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
     );
 
     pEqn.setReference(pRefCell, pRefValue);

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H

@@ -1,11 +1,9 @@
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
-volVectorField HbyA("HbyA", U);
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-HbyA = rAU*UEqn().H();
-
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
-    (fvc::interpolate(HbyA) & mesh.Sf())
+    fvc::flux(HbyA)
   + fvc::interpolate(rAU)*fvc::ddtCorr(U, Uf)
 );
 
@@ -22,7 +20,7 @@ tmp<volScalarField> rAtU(rAU);
 
 if (pimple.consistent())
 {
-    rAtU = 1.0/max(1.0/rAU - UEqn().H1(), 0.1/rAU);
+    rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
     phiHbyA +=
         fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU())*fvc::grad(p);
@@ -30,27 +28,18 @@ if (pimple.consistent())
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
-surfaceScalarField rAUf("rAUf", fvc::interpolate(rAtU()));
+// Update the pressure BCs to ensure flux consistency
-
+constrainPressure(p, U, phiHbyA, rAtU(), MRF);
-// Update the fixedFluxPressure BCs to ensure flux consistency
-setSnGrad<fixedFluxPressureFvPatchScalarField>
-(
-    p.boundaryField(),
-    (
-        phiHbyA.boundaryField()
-      - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
-    )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-);
 
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix pEqn
     (
-        fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
+        fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
     );
 
     pEqn.setReference(pRefCell, pRefValue);

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.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
@@ -42,7 +42,6 @@ Description
 #include "pimpleControl.H"
 #include "CorrectPhi.H"
 #include "fvOptions.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/incompressible/pimpleFoam/pimpleFoam.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
@@ -72,7 +72,6 @@ Description
 #include "turbulentTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/incompressible/pisoFoam/pEqn.H

@@ -1,10 +1,9 @@
 volScalarField rAU(1.0/UEqn.A());
-volVectorField HbyA("HbyA", U);
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-HbyA = rAU*UEqn.H();
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
-    (fvc::interpolate(HbyA) & mesh.Sf())
+    fvc::flux(HbyA)
   + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
 );
 
@@ -12,6 +11,9 @@ MRF.makeRelative(phiHbyA);
 
 adjustPhi(phiHbyA, U, p);
 
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, U, phiHbyA, rAU, MRF);
+
 // Non-orthogonal pressure corrector loop
 while (piso.correctNonOrthogonal())
 {

applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.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
@@ -114,7 +114,7 @@ int main(int argc, char *argv[])
                 surfaceScalarField phiHbyA
                 (
                     "phiHbyA",
-                    (fvc::interpolate(HbyA) & mesh.Sf())
+                    fvc::flux(HbyA)
                   + fvc::interpolate(rAU)*fvc::ddtCorr(h, hU, phi)
                   - phih0
                 );

applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H

@@ -1,6 +1,6 @@
     // Relative momentum predictor
 
-    tmp<fvVectorMatrix> UrelEqn
+    tmp<fvVectorMatrix> tUrelEqn
     (
         fvm::div(phi, Urel)
       + turbulence->divDevReff(Urel)
@@ -8,11 +8,15 @@
      ==
         fvOptions(Urel)
     );
+    fvVectorMatrix& UrelEqn = tUrelEqn.ref();
 
-    UrelEqn().relax();
+    UrelEqn.relax();
 
-    fvOptions.constrain(UrelEqn());
+    fvOptions.constrain(UrelEqn);
 
-    solve(UrelEqn() == -fvc::grad(p));
+    if (simple.momentumPredictor())
+    {
+        solve(UrelEqn == -fvc::grad(p));
 
-    fvOptions.correct(Urel);
+        fvOptions.correct(Urel);
+    }

applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/createFields.H

@@ -64,7 +64,7 @@ volVectorField U
 
 singlePhaseTransportModel laminarTransport(U, phi);
 
-autoPtr<incompressible::RASModel> turbulence
+autoPtr<incompressible::turbulenceModel> turbulence
 (
-    incompressible::New<incompressible::RASModel>(U, phi, laminarTransport)
+    incompressible::turbulenceModel::New(U, phi, laminarTransport)
 );