Solvers: Added support for extrapolated pressure boundary conditions

The boundary conditions of HbyA are now constrained by the new "constrainHbyA"
function which applies the velocity boundary values for patches for which the
velocity cannot be modified by assignment and pressure extrapolation is
not specified via the new
"fixedFluxExtrapolatedPressureFvPatchScalarField".

The new function "constrainPressure" sets the pressure gradient
appropriately for "fixedFluxPressureFvPatchScalarField" and
"fixedFluxExtrapolatedPressureFvPatchScalarField" boundary conditions to
ensure the evaluated flux corresponds to the known velocity values at
the boundary.

The "fixedFluxPressureFvPatchScalarField" boundary condition operates
exactly as before, ensuring the correct flux at fixed-flux boundaries by
compensating for the body forces (gravity in particular) with the
pressure gradient.

The new "fixedFluxExtrapolatedPressureFvPatchScalarField" boundary
condition may be used for cases with or without body-forces to set the
pressure gradient to compensate not only for the body-force but also the
extrapolated "HbyA" which provides a second-order boundary condition for
pressure.  This is useful for a range a problems including impinging
flow, extrapolated inlet conditions with body-forces or for highly
viscous flows, pressure-induced separation etc.  To test this boundary
condition at walls in the motorBike tutorial case set

    lowerWall
    {
        type            fixedFluxExtrapolatedPressure;
    }

    motorBikeGroup
    {
        type            fixedFluxExtrapolatedPressure;
    }

Currently the new extrapolated pressure boundary condition is supported
for all incompressible and sub-sonic compressible solvers except those
providing implicit and tensorial porosity support.  The approach will be
extended to cover these solvers and options in the future.

Note: the extrapolated pressure boundary condition is experimental and
requires further testing to assess the range of applicability,
stability, accuracy etc.

Henry G. Weller
CFD Direct Ltd.

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
@@ -88,9 +88,7 @@ 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",
@@ -98,6 +96,9 @@ int main(int argc, char *argv[])
               + 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/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())
 {

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())
 {
@@ -52,6 +50,9 @@ else
 
     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/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())
 {
@@ -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())
         );
@@ -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/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
@@ -41,7 +41,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());
 
@@ -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/pEqn.H

@@ -5,9 +5,7 @@ rho.relax();
 
 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)
 {
@@ -60,6 +58,9 @@ else
 
     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

@@ -5,8 +5,7 @@ rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
 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)
 {
@@ -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());
 
@@ -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
@@ -36,7 +36,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

@@ -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();
 
     if (pimple.transonic())
     {
@@ -65,6 +63,9 @@
 
         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/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.H

@@ -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

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -5,9 +5,7 @@ rho.relax();
 
 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)
 {
@@ -60,6 +58,9 @@ else
 
     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

@@ -5,8 +5,7 @@ rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
 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)
 {
@@ -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

@@ -5,9 +5,7 @@ rho.relax();
 
 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)
 {
@@ -60,6 +58,9 @@ 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())
     {
         // Pressure corrector

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -1,8 +1,7 @@
 {
     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));
-
     UEqn.clear();
 
     bool closedVolume = false;
@@ -23,7 +22,7 @@
             fvScalarMatrix pEqn
             (
                 fvm::div(phid, p)
-              - fvm::laplacian(rho*rAU, p)
+              - fvm::laplacian(rhorAUf, p)
               ==
                 fvOptions(psi, p, rho.name())
             );
@@ -53,12 +52,15 @@
 
         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,9 +1,6 @@
 volScalarField rAU(1.0/UEqn().A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
-
+volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
-
 UEqn.clear();
 
 bool closedVolume = false;
@@ -70,6 +67,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/pEqn.H

@@ -1,16 +1,16 @@
 {
     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();
 
     UEqn.clear();
 

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -2,10 +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 phid
 (
     "phid",

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

@@ -2,10 +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 phid
 (
     "phid",

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
@@ -97,10 +97,7 @@ 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",
@@ -108,6 +105,9 @@ int main(int argc, char *argv[])
                   + rAUf*fvc::ddtCorr(U, phi)
                 );
 
+                // Update the pressure BCs to ensure flux consistency
+                constrainPressure(p, U, phiHbyA, rAUf);
+
                 while (piso.correctNonOrthogonal())
                 {
                     fvScalarMatrix pEqn

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

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H

@@ -1,9 +1,7 @@
 {
     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());
 
@@ -17,15 +15,8 @@
 
     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/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
@@ -51,7 +51,6 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "simpleControl.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H

@@ -1,9 +1,8 @@
 {
     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();
     UEqn.clear();
 
     surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
@@ -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
@@ -39,7 +39,6 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "pimpleControl.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -6,35 +6,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())
-          + 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())
-    );
 
     fvScalarMatrix p_rghDDtEqn
     (
@@ -49,7 +39,7 @@
         fvScalarMatrix p_rghEqn
         (
             p_rghDDtEqn
-          - fvm::laplacian(rAUf, p_rgh)
+          - fvm::laplacian(rhorAUf, p_rgh)
         );
 
         p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
@@ -64,7 +54,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/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
@@ -36,7 +36,6 @@ Description
 #include "radiationModel.H"
 #include "simpleControl.H"
 #include "fvOptions.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -4,9 +4,7 @@
 
     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();
     UEqn.clear();
 
     surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
@@ -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
@@ -45,7 +45,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
@@ -38,7 +38,6 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "coordinateSystem.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -6,9 +6,7 @@
 
     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();
     UEqn.clear();
 
     surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
@@ -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/pEqn.H

@@ -7,9 +7,7 @@
 
     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());
 
@@ -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/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
@@ -125,8 +125,7 @@ int main(int argc, char *argv[])
             fvOptions.correct(U);
 
             volScalarField rAU(1.0/UEqn().A());
-            volVectorField HbyA("HbyA", U);
+            volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
-            HbyA = rAU*UEqn().H();
             UEqn.clear();
             surfaceScalarField phiHbyA
             (
@@ -135,6 +134,9 @@ int main(int argc, char *argv[])
             );
             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())
             {

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
@@ -73,9 +73,7 @@ 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",
@@ -85,6 +83,9 @@ int main(int argc, char *argv[])
 
             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
@@ -77,9 +77,7 @@ 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",
@@ -89,6 +87,9 @@ int main(int argc, char *argv[])
 
             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/pEqn.H

@@ -26,6 +26,9 @@ if (pimple.nCorrPISO() <= 1)
     UrelEqn.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/pEqn.H

@@ -1,7 +1,5 @@
 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",
@@ -28,17 +26,8 @@ if (pimple.nCorrPISO() <= 1)
     UEqn.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,7 +1,5 @@
 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",
@@ -33,24 +31,15 @@ if (pimple.nCorrPISO() <= 1)
     UEqn.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
@@ -39,7 +39,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
@@ -39,7 +39,6 @@ Description
 #include "turbulentTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/incompressible/pisoFoam/pEqn.H

@@ -1,6 +1,5 @@
 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",
@@ -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/simpleFoam/SRFSimpleFoam/pEqn.H

@@ -18,6 +18,9 @@
 
     UrelEqn.clear();
 
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, Urel, phiHbyA, rAtUrel());
+
     // Non-orthogonal pressure corrector loop
     while (simple.correctNonOrthogonal())
     {

applications/solvers/incompressible/simpleFoam/pEqn.H

@@ -1,8 +1,6 @@
 {
     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());
     MRF.makeRelative(phiHbyA);
     adjustPhi(phiHbyA, U, p);
@@ -19,6 +17,9 @@
 
     UEqn.clear();
 
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, U, phiHbyA, rAtU(), MRF);
+
     // Non-orthogonal pressure corrector loop
     while (simple.correctNonOrthogonal())
     {

applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H

@@ -1,12 +1,13 @@
-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();
 
 UEqn.clear();
 surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());

applications/solvers/lagrangian/DPMFoam/DPMFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,6 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "PhaseIncompressibleTurbulenceModel.H"
 #include "pimpleControl.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 #ifdef MPPIC
     #include "basicKinematicMPPICCloud.H"

applications/solvers/lagrangian/DPMFoam/pEqn.H

@@ -12,15 +12,8 @@
         )
     );
 
-    // Update the fixedFluxPressure BCs to ensure flux consistency
+    // Update the pressure BCs to ensure flux consistency
-    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    constrainPressure(p, Uc, phiHbyA, rAUcf);
-    (
-        p.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-            - (mesh.Sf().boundaryField() & Uc.boundaryField())
-        )/(mesh.magSf().boundaryField()*rAUcf.boundaryField())
-    );
 
     // Non-orthogonal pressure corrector loop
     while (pimple.correctNonOrthogonal())

applications/solvers/lagrangian/coalChemistryFoam/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())
 {
@@ -53,6 +51,9 @@ else
 
     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/lagrangian/reactingParcelFilmFoam/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());
 
@@ -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/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.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 "SLGThermo.H"
 #include "fvOptions.H"
 #include "pimpleControl.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/lagrangian/reactingParcelFoam/pEqn.H

@@ -7,10 +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 phiHbyA
     (
         "phiHbyA",
@@ -22,6 +19,9 @@
 
     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/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H

@@ -4,11 +4,9 @@
     thermo.rho() -= psi*p;
 
     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));
-
     UEqn.clear();
-
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
@@ -17,12 +15,15 @@
 
     MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+    // 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)
          ==
             parcels.Srho()
           + fvOptions(psi, p, rho.name())

applications/solvers/lagrangian/sprayFoam/pEqn.H

@@ -5,9 +5,7 @@ rho.relax();
 
 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)
 {
@@ -61,6 +59,9 @@ else
 
     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/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H

@@ -5,9 +5,7 @@ rho.relax();
 
 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)
 {
@@ -61,6 +59,9 @@ 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

applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H

@@ -13,9 +13,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();
 
     phi = (fvc::interpolate(HbyA) & mesh.Sf())
          + rhorAUf*fvc::ddtCorr(U, Uf);

applications/solvers/multiphase/cavitatingFoam/pEqn.H

@@ -13,9 +13,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();
 
     phi = (fvc::interpolate(HbyA) & mesh.Sf())
          + rhorAUf*fvc::ddtCorr(U, phi);

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.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
@@ -47,7 +47,6 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "pimpleControl.H"
 #include "CorrectPhi.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H

@@ -1,10 +1,7 @@
 {
     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 phiHbyA
     (
         "phiHbyA",
@@ -22,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);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - (mesh.Sf().boundaryField() & U.boundaryField())
-        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-    );
 
     // Make the fluxes relative to the mesh motion
     fvc::makeRelative(phiHbyA, U);

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.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 "twoPhaseMixtureThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "pimpleControl.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -1,10 +1,7 @@
 {
     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 phiHbyA
     (
         "phiHbyA",
@@ -22,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);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - (mesh.Sf().boundaryField() & U.boundaryField())
-        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-    );
 
     tmp<fvScalarMatrix> p_rghEqnComp1;
     tmp<fvScalarMatrix> p_rghEqnComp2;

applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C

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

applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H

@@ -1,10 +1,7 @@
 {
     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 phiHbyA
     (
         "phiHbyA",
@@ -22,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);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - (mesh.Sf().boundaryField() & U.boundaryField())
-        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-    );
 
     PtrList<fvScalarMatrix> p_rghEqnComps(mixture.phases().size());
 

applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.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 "CompressibleTurbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 #include "gaussLaplacianScheme.H"
 #include "uncorrectedSnGrad.H"
 

applications/solvers/multiphase/driftFluxFoam/pEqn.H

@@ -1,10 +1,7 @@
 {
     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 phiHbyA
     (
         "phiHbyA",
@@ -23,15 +20,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 (pimple.correctNonOrthogonal())
     {

applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.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 "pimpleControl.H"
 #include "fvOptions.H"
 #include "CorrectPhi.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 

applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H

@@ -1,10 +1,7 @@
 {
     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 phiHbyA
     (
         "phiHbyA",
@@ -30,15 +27,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);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - (mesh.Sf().boundaryField() & U.boundaryField())
-        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-    );
 
     while (pimple.correctNonOrthogonal())
     {

applications/solvers/multiphase/interFoam/interFoam.C

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

applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.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
@@ -38,7 +38,6 @@ Description
 #include "pimpleControl.H"
 #include "fvOptions.H"
 #include "CorrectPhi.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 

applications/solvers/multiphase/interFoam/pEqn.H

@@ -2,8 +2,7 @@
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
-    volVectorField HbyA("HbyA", U);
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-    HbyA = rAU*UEqn.H();
 
     surfaceScalarField phiHbyA
     (
@@ -24,15 +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, 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/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.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
@@ -51,7 +51,6 @@ Description
 #include "pimpleControl.H"
 #include "fvOptions.H"
 #include "CorrectPhi.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H

@@ -1,10 +1,7 @@
 {
     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 phiHbyA
     (
         "phiHbyA",
@@ -29,15 +26,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);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - (mesh.Sf().boundaryField() & U.boundaryField())
-        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-    );
 
     Pair<tmp<volScalarField>> vDotP = mixture->vDotP();
     const volScalarField& vDotcP = vDotP[0]();

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C

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

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -1,10 +1,7 @@
 {
     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 phiHbyA
     (
         "phiHbyA",
@@ -23,15 +20,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);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - (mesh.Sf().boundaryField() & U.boundaryField())
-        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-    );
 
     Pair<tmp<volScalarField>> vDotP = mixture->vDotP();
     const volScalarField& vDotcP = vDotP[0]();

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.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
@@ -40,7 +40,6 @@ Description
 #include "pimpleControl.H"
 #include "IOMRFZoneList.H"
 #include "CorrectPhi.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterDyMFoam/multiphaseInterDyMFoam.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 "pimpleControl.H"
 #include "fvOptions.H"
 #include "CorrectPhi.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.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
@@ -38,7 +38,6 @@ Description
 #include "pimpleControl.H"
 #include "fvOptions.H"
 #include "CorrectPhi.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H

@@ -1,8 +1,6 @@
 volScalarField rAU(1.0/UEqn().A());
 surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-
+volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p_gh));
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -19,15 +17,8 @@ surfaceScalarField phiHbyA
 MRF.makeRelative(phiHbyA);
 adjustPhi(phiHbyA, U, p_gh);
 
-// Update the fixedFluxPressure BCs to ensure flux consistency
+// Update the pressure BCs to ensure flux consistency
-setSnGrad<fixedFluxPressureFvPatchScalarField>
+constrainPressure(p_gh, U, phiHbyA, rAUf);
-(
-    p_gh.boundaryField(),
-    (
-        phiHbyA.boundaryField()
-      - (mesh.Sf().boundaryField() & U.boundaryField())
-    )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-);
 
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())

applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H

@@ -1,9 +1,7 @@
 {
     rAU = 1.0/UEqn().A();
     surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-
+    volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p_gh));
-    volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn().H();
 
     if (pimple.nCorrPISO() <= 1)
     {
@@ -26,15 +24,8 @@
         fvc::makeAbsolute(phiHbyA, U);
     }
 
-    // Update the fixedFluxPressure BCs to ensure flux consistency
+    // Update the pressure BCs to ensure flux consistency
-    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    constrainPressure(p_gh, U, phiHbyA, rAUf);
-    (
-        p_gh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - (mesh.Sf().boundaryField() & U.boundaryField())
-        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-    );
 
     // Non-orthogonal pressure corrector loop
     while (pimple.correctNonOrthogonal())

applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/potentialFreeSurfaceDyMFoam.C

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

applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.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 "turbulentTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/reactingMultiphaseEulerFoam.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
@@ -35,7 +35,6 @@ Description
 
 #include "fvCFD.H"
 #include "multiphaseSystem.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 #include "pimpleControl.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.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 "twoPhaseSystem.H"
 #include "phaseCompressibleTurbulenceModel.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 #include "pimpleControl.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"

applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H

@@ -1,10 +1,7 @@
 {
     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 phiHbyA
     (
         "phiHbyA",
@@ -20,15 +17,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);
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - (mesh.Sf().boundaryField() & U.boundaryField())
-        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
-    );
 
     while (pimple.correctNonOrthogonal())
     {

applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.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 "incompressibleTwoPhaseMixture.H"
 #include "turbulentTransportModel.H"
 #include "pimpleControl.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.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
@@ -35,7 +35,6 @@ Description
 #include "PhaseCompressibleTurbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
 #include "fixedValueFvsPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteVolume/Make/files

@@ -138,6 +138,7 @@ $(derivedFvPatchFields)/externalCoupledMixed/externalCoupledMixedFvPatchFields.C
 $(derivedFvPatchFields)/fan/fanFvPatchFields.C
 $(derivedFvPatchFields)/fanPressure/fanPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
+$(derivedFvPatchFields)/fixedFluxExtrapolatedPressure/fixedFluxExtrapolatedPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/fixedInternalValueFvPatchField/fixedInternalValueFvPatchFields.C
 $(derivedFvPatchFields)/fixedJump/fixedJumpFvPatchFields.C
 $(derivedFvPatchFields)/fixedJumpAMI/fixedJumpAMIFvPatchFields.C
@@ -394,6 +395,7 @@ finiteVolume/fvc/fvcReconstructMag.C
 
 general = cfdTools/general
 $(general)/findRefCell/findRefCell.C
+$(general)/constrainHbyA/constrainHbyA.C
 $(general)/adjustPhi/adjustPhi.C
 $(general)/bound/bound.C
 $(general)/CorrectPhi/correctUphiBCs.C

src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.C

@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "constrainHbyA.H"
+#include "volFields.H"
+#include "fixedFluxExtrapolatedPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volVectorField> Foam::constrainHbyA
+(
+    const tmp<volVectorField>& tHbyA,
+    const volVectorField& U,
+    const volScalarField& p
+)
+{
+    tmp<volVectorField> tHbyANew;
+
+    if (tHbyA.isTmp())
+    {
+        volVectorField* HbyAPtr = tHbyA.ptr();
+        HbyAPtr->rename("HbyA");
+        tHbyANew = HbyAPtr;
+    }
+    else
+    {
+        tHbyANew = new volVectorField("HbyA", tHbyA);
+    }
+
+    volVectorField& HbyA = tHbyANew();
+
+    forAll(U.boundaryField(), patchi)
+    {
+        if
+        (
+           !U.boundaryField()[patchi].assignable()
+        && !isA<fixedFluxExtrapolatedPressureFvPatchScalarField>
+            (
+                p.boundaryField()[patchi]
+            )
+        )
+        {
+            HbyA.boundaryField()[patchi] = U.boundaryField()[patchi];
+        }
+    }
+
+    return tHbyANew;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.H

@@ -0,0 +1,61 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::constrainHbyA
+
+Description
+
+SourceFiles
+    constrainHbyA.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constrainHbyA_H
+#define constrainHbyA_H
+
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+tmp<volVectorField> constrainHbyA
+(
+    const tmp<volVectorField>& tHbyA,
+    const volVectorField& U,
+    const volScalarField& p
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/constrainPressure/constrainPressure.C

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "constrainPressure.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "geometricOneField.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class RhoType, class RAUType, class MRFType>
+void Foam::constrainPressure
+(
+    volScalarField& p,
+    const RhoType& rho,
+    const volVectorField& U,
+    const surfaceScalarField& phiHbyA,
+    const RAUType& rhorAU,
+    const MRFType& MRF
+)
+{
+    const fvMesh& mesh = p.mesh();
+    volScalarField::GeometricBoundaryField& pBf = p.boundaryField();
+
+    bool hasFixedFluxBCs = false;
+    forAll(pBf, patchi)
+    {
+        if (isA<fixedFluxPressureFvPatchScalarField>(pBf[patchi]))
+        {
+            hasFixedFluxBCs = true;
+            break;
+        }
+    }
+
+    if (hasFixedFluxBCs)
+    {
+        const surfaceScalarField::GeometricBoundaryField& phiHbyABf =
+            phiHbyA.boundaryField();
+        const typename RAUType::GeometricBoundaryField& rhorAUBf =
+            rhorAU.boundaryField();
+        const surfaceScalarField::GeometricBoundaryField& magSfBf =
+            mesh.magSf().boundaryField();
+
+        // Pre-compute tho relative flux for all patches: currently MRFZone does
+        // not support computing the relative flux for individual patches
+        FieldField<fvsPatchField, scalar> phiRelBf
+        (
+            MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
+        );
+
+        forAll(pBf, patchi)
+        {
+            if (isA<fixedFluxPressureFvPatchScalarField>(pBf[patchi]))
+            {
+                refCast<fixedFluxPressureFvPatchScalarField>(pBf[patchi])
+               .updateCoeffs
+                (
+                    (
+                        phiHbyABf[patchi]
+                      - rho.boundaryField()[patchi]*phiRelBf[patchi]
+                    )
+                   /(magSfBf[patchi]*rhorAUBf[patchi])
+                );
+            }
+        }
+    }
+}
+
+
+template<class RAUType>
+void Foam::constrainPressure
+(
+    volScalarField& p,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& phiHbyA,
+    const RAUType& rAU
+)
+{
+    constrainPressure(p, rho, U, phiHbyA, rAU, NullMRF());
+}
+
+
+template<class RAUType, class MRFType>
+void Foam::constrainPressure
+(
+    volScalarField& p,
+    const volVectorField& U,
+    const surfaceScalarField& phiHbyA,
+    const RAUType& rAU,
+    const MRFType& MRF
+)
+{
+    constrainPressure(p, geometricOneField(), U, phiHbyA, rAU, MRF);
+}
+
+
+template<class RAUType>
+void Foam::constrainPressure
+(
+    volScalarField& p,
+    const volVectorField& U,
+    const surfaceScalarField& phiHbyA,
+    const RAUType& rAU
+)
+{
+    constrainPressure(p, U, phiHbyA, rAU, NullMRF());
+}
+
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/constrainPressure/constrainPressure.H

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::constrainPressure
+
+Description
+
+SourceFiles
+    constrainPressure.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constrainPressure_H
+#define constrainPressure_H
+
+#include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
+#include "IOMRFZoneList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+/*---------------------------------------------------------------------------*\
+                           Class NullMRF Declaration
+\*---------------------------------------------------------------------------*/
+
+class NullMRF
+{
+public:
+
+    // Constructors
+
+        //- Construct null
+        NullMRF()
+        {}
+
+
+    // Member Functions
+
+        //- Return the argument unchanged
+        template<class Type>
+        inline const Type& relative(const Type& U) const
+        {
+            return U;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class RhoType, class RAUType, class MRFType>
+void constrainPressure
+(
+    volScalarField& p,
+    const RhoType& rho,
+    const volVectorField& U,
+    const surfaceScalarField& phiHbyA,
+    const RAUType& rhorAU,
+    const MRFType& MRF
+);
+
+template<class RAUType>
+void constrainPressure
+(
+    volScalarField& p,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& phiHbyA,
+    const RAUType& rhorAU
+);
+
+template<class RAUType, class MRFType>
+void constrainPressure
+(
+    volScalarField& p,
+    const volVectorField& U,
+    const surfaceScalarField& phiHbyA,
+    const RAUType& rAU,
+    const MRFType& MRF
+);
+
+template<class RAUType>
+void constrainPressure
+(
+    volScalarField& p,
+    const volVectorField& U,
+    const surfaceScalarField& phiHbyA,
+    const RAUType& rAU
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "constrainPressure.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/include/fvCFD.H

@@ -12,6 +12,9 @@
 #include "uniformDimensionedFields.H"
 #include "calculatedFvPatchFields.H"
 #include "fixedValueFvPatchFields.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+#include "constrainHbyA.H"
+#include "constrainPressure.H"
 #include "adjustPhi.H"
 #include "findRefCell.H"
 #include "IOMRFZoneList.H"

src/finiteVolume/fields/fvPatchFields/derived/fixedFluxExtrapolatedPressure/fixedFluxExtrapolatedPressureFvPatchScalarField.C

@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fixedFluxExtrapolatedPressureFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fixedFluxExtrapolatedPressureFvPatchScalarField::
+fixedFluxExtrapolatedPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedFluxPressureFvPatchScalarField(p, iF)
+{}
+
+
+Foam::fixedFluxExtrapolatedPressureFvPatchScalarField::
+fixedFluxExtrapolatedPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedFluxPressureFvPatchScalarField(p, iF, dict)
+{}
+
+
+Foam::fixedFluxExtrapolatedPressureFvPatchScalarField::
+fixedFluxExtrapolatedPressureFvPatchScalarField
+(
+    const fixedFluxExtrapolatedPressureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedFluxPressureFvPatchScalarField(ptf, p, iF, mapper)
+{}
+
+
+Foam::fixedFluxExtrapolatedPressureFvPatchScalarField::
+fixedFluxExtrapolatedPressureFvPatchScalarField
+(
+    const fixedFluxExtrapolatedPressureFvPatchScalarField& wbppsf
+)
+:
+    fixedFluxPressureFvPatchScalarField(wbppsf)
+{}
+
+
+Foam::fixedFluxExtrapolatedPressureFvPatchScalarField::
+fixedFluxExtrapolatedPressureFvPatchScalarField
+(
+    const fixedFluxExtrapolatedPressureFvPatchScalarField& wbppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedFluxPressureFvPatchScalarField(wbppsf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        fixedFluxExtrapolatedPressureFvPatchScalarField
+    );
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/fixedFluxExtrapolatedPressure/fixedFluxExtrapolatedPressureFvPatchScalarField.H

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fixedFluxExtrapolatedPressureFvPatchScalarField
+
+Group
+    grpInletBoundaryConditions grpWallBoundaryConditions
+
+Description
+    This boundary condition sets the pressure gradient to the provided value
+    such that the flux on the boundary is that specified by the velocity
+    boundary condition.
+
+    Example of the boundary condition specification:
+    \verbatim
+    myPatch
+    {
+        type            fixedFluxExtrapolatedPressure;
+    }
+    \endverbatim
+
+SeeAlso
+    Foam::fixedGradientFvPatchField
+
+SourceFiles
+    fixedFluxExtrapolatedPressureFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedFluxExtrapolatedPressureFvPatchScalarFields_H
+#define fixedFluxExtrapolatedPressureFvPatchScalarFields_H
+
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+             Class fixedFluxExtrapolatedPressureFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class fixedFluxExtrapolatedPressureFvPatchScalarField
+:
+    public fixedFluxPressureFvPatchScalarField
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("fixedFluxExtrapolatedPressure");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        fixedFluxExtrapolatedPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        fixedFluxExtrapolatedPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  fixedFluxExtrapolatedPressureFvPatchScalarField onto a new patch
+        fixedFluxExtrapolatedPressureFvPatchScalarField
+        (
+            const fixedFluxExtrapolatedPressureFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        fixedFluxExtrapolatedPressureFvPatchScalarField
+        (
+            const fixedFluxExtrapolatedPressureFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new fixedFluxExtrapolatedPressureFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        fixedFluxExtrapolatedPressureFvPatchScalarField
+        (
+            const fixedFluxExtrapolatedPressureFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new fixedFluxExtrapolatedPressureFvPatchScalarField(*this, iF)
+            );
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/fvSchemes

@@ -38,13 +38,7 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
+    default         Gauss linear uncorrected;
-    laplacian((rho*nuEff),U) Gauss linear uncorrected;
-    laplacian(rhorAUf,p_rgh) Gauss linear uncorrected;
-    laplacian(alphaEff,h) Gauss linear uncorrected;
-    laplacian((rho*DkEff),k) Gauss linear uncorrected;
-    laplacian((rho*DepsilonEff),epsilon) Gauss linear uncorrected;
-    laplacian((rho*DomegaEff),omega) Gauss linear uncorrected;
 }
 
 interpolationSchemes

tutorials/incompressible/pimpleFoam/pitzDaily/system/fvSchemes

@@ -23,8 +23,6 @@ ddtSchemes
 gradSchemes
 {
     default         Gauss linear;
-    grad(p)         Gauss linear;
-    grad(U)         Gauss linear;
 }
 
 divSchemes
@@ -41,19 +39,12 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
+    default         Gauss linear corrected;
-    laplacian(nuEff,U) Gauss linear corrected;
-    laplacian(rAUf,p)  Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian(DREff,R) Gauss linear corrected;
-    laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
 }
 
 interpolationSchemes
 {
     default         linear;
-    interpolate(U)  linear;
 }
 
 snGradSchemes