Specialized dotInterpolate for the efficient calculation of flux fields

e.g. (fvc::interpolate(HbyA) & mesh.Sf()) -> fvc::flux(HbyA) This removes the need to create an intermediate face-vector field when computing fluxes which is more efficient, reduces the peak storage and improved cache coherency in addition to providing a simpler and cleaner API.
2016-04-06 20:20:53 +01:00
parent d52546b3ce
commit fda9aadb3a
101 changed files with 523 additions and 451 deletions
--- a/applications/solvers/DNS/dnsFoam/dnsFoam.C
+++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C
@ -92,7 +92,7 @@ int main(int argc, char *argv[])
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
-                (fvc::interpolate(HbyA) & mesh.Sf())
+                fvc::flux(HbyA)
              + rAUf*fvc::ddtCorr(U, phi)
            );

--- a/applications/solvers/basic/potentialFoam/createFields.H
+++ b/applications/solvers/basic/potentialFoam/createFields.H
@ -24,13 +24,13 @@ surfaceScalarField phi
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
-    fvc::interpolate(U) & mesh.Sf()
+    fvc::flux(U)
 );

 if (args.optionFound("initialiseUBCs"))
 {
    U.correctBoundaryConditions();
-    phi = fvc::interpolate(U) & mesh.Sf();
+    phi = fvc::flux(U);
 }


--- a/applications/solvers/basic/potentialFoam/potentialFoam.C
+++ b/applications/solvers/basic/potentialFoam/potentialFoam.C
@ -2,7 +2,7 @@
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -121,8 +121,7 @@ int main(int argc, char *argv[])
    U.correctBoundaryConditions();

    Info<< "Interpolated velocity error = "
-        << (sqrt(sum(sqr((fvc::interpolate(U) & mesh.Sf()) - phi)))
-          /sum(mesh.magSf())).value()
+        << (sqrt(sum(sqr(fvc::flux(U) - phi)))/sum(mesh.magSf())).value()
        << endl;

    // Write U and phi
--- a/applications/solvers/combustion/PDRFoam/pEqn.H
+++ b/applications/solvers/combustion/PDRFoam/pEqn.H
@ -10,7 +10,7 @@ if (pimple.transonic())
        "phid",
        fvc::interpolate(psi)
       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
           /fvc::interpolate(rho)
        )
@ -41,7 +41,7 @@ else
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
        )
    );
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/pEqn.H
@ -11,7 +11,7 @@ if (pimple.transonic())
        "phid",
        fvc::interpolate(psi)
       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
        )
    );
@ -43,7 +43,7 @@ else
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)
        )
    );
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/pEqn.H
@ -12,7 +12,7 @@ if (pimple.transonic())
        fvc::interpolate(psi)
       *(
            (
-                (fvc::interpolate(HbyA) & mesh.Sf())
+                fvc::flux(HbyA)
              + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
            )
        )
@ -46,7 +46,7 @@ else
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
        )
    );
--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/applications/solvers/combustion/fireFoam/pEqn.H
@ -10,7 +10,7 @@ surfaceScalarField phiHbyA
 (
    "phiHbyA",
    (
-        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+        fvc::flux(rho*HbyA)
      + rhorAUf*fvc::ddtCorr(rho, U, phi)
    )
  + phig
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@ -19,7 +19,7 @@ if (pimple.transonic())
        "phid",
        fvc::interpolate(psi)
       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
        )
    );
@ -51,7 +51,7 @@ else
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)
        )
    );
--- a/applications/solvers/combustion/reactingFoam/pcEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pcEqn.H
@ -19,7 +19,7 @@ if (pimple.transonic())
        "phid",
        fvc::interpolate(psi)
       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
           /fvc::interpolate(rho)
        )
@ -63,7 +63,7 @@ else
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
        )
    );
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
@ -15,7 +15,7 @@
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)
        )
      + phig
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
@ -15,7 +15,7 @@
        (
            "phiHbyA",
            (
-                (fvc::interpolate(HbyA) & mesh.Sf())
+                fvc::flux(HbyA)
              + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
            )
        );
@ -56,7 +56,7 @@
        (
            "phiHbyA",
            (
-                (fvc::interpolate(rho*HbyA) & mesh.Sf())
+                fvc::flux(rho*HbyA)
              + rhorAUf*fvc::ddtCorr(rho, U, phi)
            )
        );
--- a/applications/solvers/compressible/rhoCentralFoam/createFields.H
+++ b/applications/solvers/compressible/rhoCentralFoam/createFields.H
@ -95,7 +95,7 @@ surfaceScalarField neg
    dimensionedScalar("neg", dimless, -1.0)
 );

-surfaceScalarField phi("phi", mesh.Sf() & fvc::interpolate(rhoU));
+surfaceScalarField phi("phi", fvc::flux(rhoU));

 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
@ -208,7 +208,7 @@ int main(int argc, char *argv[])
            "sigmaDotU",
            (
                fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
-              + (mesh.Sf() & fvc::interpolate(tauMC))
+              + fvc::dotInterpolate(mesh.Sf(), tauMC)
            )
            & (a_pos*U_pos + a_neg*U_neg)
        );
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
@ -201,9 +201,9 @@ int main(int argc, char *argv[])
            "sigmaDotU",
            (
                fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
-              + (mesh.Sf() & fvc::interpolate(tauMC))
+              + fvc::dotInterpolate(mesh.Sf(), tauMC)
            )
-            & (a_pos*U_pos + a_neg*U_neg)
+          & (a_pos*U_pos + a_neg*U_neg)
        );

        solve
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@ -19,7 +19,7 @@ if (pimple.transonic())
        "phid",
        fvc::interpolate(psi)
       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
        )
    );
@ -51,7 +51,7 @@ else
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)
        )
    );
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@ -19,7 +19,7 @@ if (pimple.transonic())
        "phid",
        fvc::interpolate(psi)
       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
           /fvc::interpolate(rho)
        )
@ -63,7 +63,7 @@ else
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
        )
    );
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@ -19,7 +19,7 @@ if (pimple.transonic())
        "phid",
        fvc::interpolate(psi)
       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
        )
    );
@ -51,7 +51,7 @@ else
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+        fvc::flux(rho*HbyA)
      + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
    );

--- a/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
@ -38,7 +38,7 @@
        surfaceScalarField phid
        (
            "phid",
-            fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
+            fvc::interpolate(psi)*fvc::flux(U)
        );

        rDeltaT.dimensionedInternalField() = max
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@ -12,7 +12,7 @@
        (
            "phid",
            fvc::interpolate(psi)
-           *(fvc::interpolate(HbyA) & mesh.Sf())
+           *fvc::flux(HbyA)
        );

        MRF.makeRelative(fvc::interpolate(psi), phid);
@ -42,12 +42,7 @@
    }
    else
    {
-        surfaceScalarField phiHbyA
-        (
-            "phiHbyA",
-            fvc::interpolate(rho*HbyA) & mesh.Sf()
-        );
-
+        surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
        MRF.makeRelative(fvc::interpolate(rho), phiHbyA);

        closedVolume = adjustPhi(phiHbyA, U, p);
--- a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
@ -11,7 +11,7 @@ if (simple.transonic())
    (
        "phid",
        fvc::interpolate(psi)
-       *(fvc::interpolate(HbyA) & mesh.Sf())
+       *fvc::flux(HbyA)
    );

    MRF.makeRelative(fvc::interpolate(psi), phid);
@ -52,12 +52,7 @@ if (simple.transonic())
 }
 else
 {
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        fvc::interpolate(rho*HbyA) & mesh.Sf()
-    );
-
+    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);

    closedVolume = adjustPhi(phiHbyA, U, p);
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H
@ -16,12 +16,7 @@

    bool closedVolume = false;

-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        fvc::interpolate(rho*HbyA) & mesh.Sf()
-    );
-
+    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);

    closedVolume = adjustPhi(phiHbyA, U, p);
--- a/applications/solvers/compressible/sonicFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/pEqn.H
@ -8,7 +8,7 @@ surfaceScalarField phid
    "phid",
    fvc::interpolate(psi)
   *(
-        (mesh.Sf() & fvc::interpolate(HbyA))
+        fvc::flux(HbyA)
      + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
    )
 );
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
@ -8,7 +8,7 @@ surfaceScalarField phid
    "phid",
    fvc::interpolate(psi)
   *(
-       (mesh.Sf() & fvc::interpolate(HbyA))
+       fvc::flux(HbyA)
     + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
    )
 );
--- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
@ -2,7 +2,7 @@
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -89,7 +89,7 @@ int main(int argc, char *argv[])
                    "phid",
                    psi
                   *(
-                       (fvc::interpolate(U) & mesh.Sf())
+                       fvc::flux(U)
                     + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
                    )
                );
--- a/applications/solvers/electromagnetics/mhdFoam/createPhiB.H
+++ b/applications/solvers/electromagnetics/mhdFoam/createPhiB.H
@ -40,7 +40,7 @@
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
-            (fvc::interpolate(B) & mesh.Sf())
+            fvc::flux(B)
        );
    }

--- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
+++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
@ -101,7 +101,7 @@ int main(int argc, char *argv[])
                surfaceScalarField phiHbyA
                (
                    "phiHbyA",
-                    (fvc::interpolate(HbyA) & mesh.Sf())
+                    fvc::flux(HbyA)
                  + rAUf*fvc::ddtCorr(U, phi)
                );

@ -147,8 +147,7 @@ int main(int argc, char *argv[])
            volScalarField rAB(1.0/BEqn.A());
            surfaceScalarField rABf("rABf", fvc::interpolate(rAB));

-            phiB = (fvc::interpolate(B) & mesh.Sf())
-                + rABf*fvc::ddtCorr(B, phiB);
+            phiB = fvc::flux(B) + rABf*fvc::ddtCorr(B, phiB);

            while (bpiso.correctNonOrthogonal())
            {
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
@ -8,7 +8,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
      + rAUf*fvc::ddtCorr(U, phi)
      + phig
    );
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
@ -10,7 +10,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
    );

    MRF.makeRelative(phiHbyA);
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@ -15,7 +15,7 @@
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)
        )
      + phig
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@ -12,7 +12,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+        fvc::flux(rho*HbyA)
    );

    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@ -14,7 +14,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+        fvc::flux(rho*HbyA)
    );

    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@ -15,7 +15,7 @@
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)
        )
      + phig
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
@ -128,11 +128,7 @@ int main(int argc, char *argv[])
            volScalarField rAU(1.0/UEqn.A());
            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
            tUEqn.clear();
-            surfaceScalarField phiHbyA
-            (
-                "phiHbyA",
-                fvc::interpolate(HbyA) & mesh.Sf()
-            );
+            surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
            adjustPhi(phiHbyA, U, p);

            // Update the pressure BCs to ensure flux consistency
@ -175,7 +171,7 @@ int main(int argc, char *argv[])
            //(
            //    fvc::reconstruct
            //    (
-            //        mesh.magSf()*(fvc::snGrad(Ua) & fvc::interpolate(U))
+            //        mesh.magSf()*fvc::dotInterpolate(fvc::snGrad(Ua), U)
            //    )
            //);

@ -204,11 +200,7 @@ int main(int argc, char *argv[])
            volVectorField HbyAa("HbyAa", Ua);
            HbyAa = rAUa*UaEqn.H();
            tUaEqn.clear();
-            surfaceScalarField phiHbyAa
-            (
-                "phiHbyAa",
-                fvc::interpolate(HbyAa) & mesh.Sf()
-            );
+            surfaceScalarField phiHbyAa("phiHbyAa", fvc::flux(HbyAa));
            adjustPhi(phiHbyAa, Ua, pa);

            // Non-orthogonal pressure corrector loop
--- a/applications/solvers/incompressible/icoFoam/icoFoam.C
+++ b/applications/solvers/incompressible/icoFoam/icoFoam.C
@ -77,7 +77,7 @@ int main(int argc, char *argv[])
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
-                (fvc::interpolate(HbyA) & mesh.Sf())
+                fvc::flux(HbyA)
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );

--- a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
+++ b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
@ -81,7 +81,7 @@ int main(int argc, char *argv[])
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
-                (fvc::interpolate(HbyA) & mesh.Sf())
+                fvc::flux(HbyA)
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );

--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
@ -5,7 +5,7 @@ HbyA = rAUrel*UrelEqn.H();
 surfaceScalarField phiHbyA
 (
    "phiHbyA",
-    (fvc::interpolate(HbyA) & mesh.Sf())
+    fvc::flux(HbyA)
  + fvc::interpolate(rAUrel)*fvc::ddtCorr(Urel, phi)
 );

--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@ -3,7 +3,7 @@ volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 surfaceScalarField phiHbyA
 (
    "phiHbyA",
-    (fvc::interpolate(HbyA) & mesh.Sf())
+    fvc::flux(HbyA)
  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
 );

--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
@ -3,7 +3,7 @@ volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 surfaceScalarField phiHbyA
 (
    "phiHbyA",
-    (fvc::interpolate(HbyA) & mesh.Sf())
+    fvc::flux(HbyA)
  + fvc::interpolate(rAU)*fvc::ddtCorr(U, Uf)
 );

--- a/applications/solvers/incompressible/pisoFoam/pEqn.H
+++ b/applications/solvers/incompressible/pisoFoam/pEqn.H
@ -3,7 +3,7 @@ volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 surfaceScalarField phiHbyA
 (
    "phiHbyA",
-    (fvc::interpolate(HbyA) & mesh.Sf())
+    fvc::flux(HbyA)
  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
 );

--- a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
+++ b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
@ -2,7 +2,7 @@
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -111,7 +111,7 @@ int main(int argc, char *argv[])
                surfaceScalarField phiHbyA
                (
                    "phiHbyA",
-                    (fvc::interpolate(HbyA) & mesh.Sf())
+                    fvc::flux(HbyA)
                  + fvc::interpolate(rAU)*fvc::ddtCorr(h, hU, phi)
                  - phih0
                );
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
@ -3,7 +3,7 @@
    volVectorField HbyA("HbyA", Urel);
    HbyA = rAUrel*UrelEqn.H();

-    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
+    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
    adjustPhi(phiHbyA, Urel, p);

    tmp<volScalarField> rAtUrel(rAUrel);
--- a/applications/solvers/incompressible/simpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/pEqn.H
@ -1,7 +1,7 @@
 {
    volScalarField rAU(1.0/UEqn.A());
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
+    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
    MRF.makeRelative(phiHbyA);
    adjustPhi(phiHbyA, U, p);

--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
@ -10,7 +10,7 @@ else
 volVectorField& HbyA = tHbyA.ref();

 tUEqn.clear();
-surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
+surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));

 MRF.makeRelative(phiHbyA);

--- a/applications/solvers/lagrangian/DPMFoam/UcEqn.H
+++ b/applications/solvers/lagrangian/DPMFoam/UcEqn.H
@ -14,8 +14,7 @@ surfaceScalarField rAUcf("Dp", fvc::interpolate(rAUc));

 surfaceScalarField phicForces
 (
-   (fvc::interpolate(rAUc*cloudVolSUSu/rhoc) & mesh.Sf())
- + rAUcf*(g & mesh.Sf())
+   fvc::flux(rAUc*cloudVolSUSu/rhoc) + rAUcf*(g & mesh.Sf())
 );

 if (pimple.momentumPredictor())
--- a/applications/solvers/lagrangian/DPMFoam/pEqn.H
+++ b/applications/solvers/lagrangian/DPMFoam/pEqn.H
@ -6,7 +6,7 @@
    (
        "phiHbyA",
        (
-           (fvc::interpolate(HbyA) & mesh.Sf())
+           fvc::flux(HbyA)
         + alphacf*rAUcf*fvc::ddtCorr(Uc, phic)
         + phicForces
        )
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@ -11,7 +11,7 @@ if (pimple.transonic())
        "phid",
        fvc::interpolate(psi)
       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
        )
    );
@ -44,7 +44,7 @@ else
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)
        )
    );
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
@ -10,7 +10,7 @@ surfaceScalarField phiHbyA
 (
    "phiHbyA",
    (
-        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+        fvc::flux(rho*HbyA)
      + rhorAUf*fvc::ddtCorr(rho, U, phi)
    )
  + phig
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@ -12,7 +12,7 @@
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)
        )
    );
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
@ -10,7 +10,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::interpolate(rho)*fvc::flux(HbyA)
    );

    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
--- a/applications/solvers/lagrangian/sprayFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/pEqn.H
@ -19,7 +19,7 @@ if (pimple.transonic())
        "phid",
        fvc::interpolate(psi)
       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
        )
    );
@ -52,7 +52,7 @@ else
    (
        "phiHbyA",
        (
-            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+            fvc::flux(rho*HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)
        )
    );
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
@ -19,7 +19,7 @@ if (pimple.transonic())
        "phid",
        fvc::interpolate(psi)
       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::flux(HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
        )
    );
@ -52,7 +52,7 @@ else
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+        fvc::flux(rho*HbyA)
      + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
    );

--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@ -15,7 +15,7 @@
    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));

-    phi = (fvc::interpolate(HbyA) & mesh.Sf())
+    phi = fvc::flux(HbyA)
         + rhorAUf*fvc::ddtCorr(U, Uf);
    fvc::makeRelative(phi, U);

--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@ -15,7 +15,7 @@
    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));

-    phi = (fvc::interpolate(HbyA) & mesh.Sf())
+    phi = fvc::flux(HbyA)
         + rhorAUf*fvc::ddtCorr(U, phi);

    surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@ -5,7 +5,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
    );

--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@ -5,7 +5,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );

--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
@ -5,7 +5,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );

--- a/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
@ -11,10 +11,7 @@
        dimensionedScalar("0", phi.dimensions(), 0)
    );

-    surfaceScalarField phir
-    (
-        mesh.Sf() & fvc::interpolate(UdmModel.Udm())
-    );
+    surfaceScalarField phir(fvc::flux(UdmModel.Udm()));

    if (nAlphaSubCycles > 1)
    {
--- a/applications/solvers/multiphase/driftFluxFoam/pEqn.H
+++ b/applications/solvers/multiphase/driftFluxFoam/pEqn.H
@ -5,7 +5,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );
    MRF.makeRelative(phiHbyA);
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
@ -5,7 +5,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
    );
    MRF.makeRelative(phiHbyA);
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@ -7,7 +7,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );
    MRF.makeRelative(phiHbyA);
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
@ -5,7 +5,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
    );

--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@ -5,7 +5,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );
    adjustPhi(phiHbyA, U, p_rgh);
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.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
@ -29,6 +29,7 @@ License
 #include "slipFvPatchFields.H"
 #include "partialSlipFvPatchFields.H"
 #include "surfaceInterpolate.H"
+#include "fvcFlux.H"

 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

@ -177,7 +178,7 @@ Foam::phaseModel::phaseModel
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
-                fvc::interpolate(U_) & mesh.Sf(),
+                fvc::flux(U_),
                phiTypes
            )
        );
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@ -97,7 +97,7 @@

        phiHbyAs[phasei] =
        (
-            (fvc::interpolate(HbyAs[phasei]) & mesh.Sf())
+            fvc::flux(HbyAs[phasei])
          + rAlphaAUfs[phasei]*fvc::ddtCorr(phase.U(), phase.phi())
        );
        MRF.makeRelative(phiHbyAs[phasei]);
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
@ -10,7 +10,7 @@ if (pimple.nCorrPISO() <= 1)
 surfaceScalarField phiHbyA
 (
    "phiHbyA",
-    (fvc::interpolate(HbyA) & mesh.Sf())
+    fvc::flux(HbyA)
  + rAUf*fvc::ddtCorr(U, phi)
 );

--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H
@ -11,7 +11,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
      + rAUf*fvc::ddtCorr(U, Uf)
    );

--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/liftModels/liftModel/liftModel.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
@ -26,6 +26,7 @@ License
 #include "liftModel.H"
 #include "phasePair.H"
 #include "fvcCurl.H"
+#include "fvcFlux.H"
 #include "surfaceInterpolate.H"

 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -78,11 +79,7 @@ Foam::tmp<Foam::volVectorField> Foam::liftModel::F() const

 Foam::tmp<Foam::surfaceScalarField> Foam::liftModel::Ff() const
 {
-    const fvMesh& mesh(this->pair_.phase1().mesh());
-
-    return
-        fvc::interpolate(pair_.dispersed())
-       *(fvc::interpolate(Fi()) & mesh.Sf());
+    return fvc::interpolate(pair_.dispersed())*fvc::flux(Fi());
 }


--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
@ -25,6 +25,7 @@ License

 #include "wallLubricationModel.H"
 #include "phasePair.H"
+#include "fvcFlux.H"
 #include "surfaceInterpolate.H"
 #include "wallFvPatch.H"

@ -91,11 +92,7 @@ Foam::tmp<Foam::volVectorField> Foam::wallLubricationModel::F() const

 Foam::tmp<Foam::surfaceScalarField> Foam::wallLubricationModel::Ff() const
 {
-    const fvMesh& mesh(this->pair_.phase1().mesh());
-
-    return
-        fvc::interpolate(pair_.dispersed())
-       *(fvc::interpolate(Fi()) & mesh.Sf());
+    return fvc::interpolate(pair_.dispersed())*fvc::flux(Fi());
 }


--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
@ -2,7 +2,7 @@
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -35,6 +35,7 @@ License
 #include "fvmSup.H"
 #include "fvcDdt.H"
 #include "fvcDiv.H"
+#include "fvcFlux.H"
 #include "surfaceInterpolate.H"
 #include "fvMatrix.H"

@ -109,7 +110,7 @@ Foam::MovingPhaseModel<BasePhaseModel>::phi(const volVectorField& U) const
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
-                fvc::interpolate(U) & U.mesh().Sf(),
+                fvc::flux(U),
                phiTypes
            )
        );
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@ -51,7 +51,7 @@ PtrList<surfaceScalarField> phiFs(phases.size());
                new surfaceScalarField
                (
                    IOobject::groupName("phiF", phase.name()),
-                    (fvc::interpolate(rAUs[phasei]*Fs()[phasei]) & mesh.Sf())
+                    fvc::flux(rAUs[phasei]*Fs()[phasei])
                )
            );
        }
@ -201,7 +201,7 @@ while (pimple.correct())
            new surfaceScalarField
            (
                IOobject::groupName("phiHbyA", phase.name()),
-                (fvc::interpolate(HbyAs[phasei]) & mesh.Sf())
+                fvc::flux(HbyAs[phasei])
              + phiCorrCoeff
               *fvc::interpolate
                (
@ -209,7 +209,7 @@ while (pimple.correct())
                )
               *(
                   MRF.absolute(phase.phi().oldTime())
-                 - (fvc::interpolate(phase.U().oldTime()) & mesh.Sf())
+                 - fvc::flux(phase.U().oldTime())
                )/runTime.deltaT()
              - phigFs[phasei]
            )
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
@ -65,18 +65,10 @@ tmp<surfaceScalarField> phiF2;
    surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());

    // Phase-1 dispersion, lift and wall-lubrication flux
-    phiF1 =
-    (
-        DbyA1()*snGradAlpha1
-      + (fvc::interpolate(rAU1*F) & mesh.Sf())
-    );
+    phiF1 = DbyA1()*snGradAlpha1 + fvc::flux(rAU1*F);

    // Phase-2 dispersion, lift and wall-lubrication flux
-    phiF2 =
-    (
-      - DbyA2()*snGradAlpha1
-      - (fvc::interpolate(rAU2*F) & mesh.Sf())
-    );
+    phiF2 = - DbyA2()*snGradAlpha1 - fvc::flux(rAU2*F);

    // Cache the phase diffusivities for implicit treatment in the
    // phase-fraction equation
@ -175,11 +167,11 @@ while (pimple.correct())
    surfaceScalarField phiHbyA1
    (
        IOobject::groupName("phiHbyA", phase1.name()),
-        (fvc::interpolate(HbyA1) & mesh.Sf())
+        fvc::flux(HbyA1)
      + phiCorrCoeff1*fvc::interpolate(alpha1.oldTime()*rho1.oldTime()*rAU1)
       *(
            MRF.absolute(phi1.oldTime())
-          - (fvc::interpolate(U1.oldTime()) & mesh.Sf())
+          - fvc::flux(U1.oldTime())
        )/runTime.deltaT()
      - phiF1()
      - phig1
@ -189,11 +181,11 @@ while (pimple.correct())
    surfaceScalarField phiHbyA2
    (
        IOobject::groupName("phiHbyA", phase2.name()),
-        (fvc::interpolate(HbyA2) & mesh.Sf())
+        fvc::flux(HbyA2)
      + phiCorrCoeff2*fvc::interpolate(alpha2.oldTime()*rho2.oldTime()*rAU2)
       *(
            MRF.absolute(phi2.oldTime())
-          - (fvc::interpolate(U2.oldTime()) & mesh.Sf())
+          - fvc::flux(U2.oldTime())
        )/runTime.deltaT()
      - phiF2()
      - phig2
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@ -157,7 +157,7 @@ while (pimple.correct())
       *(
            (alphaRhof10 + Vmf)
           *MRF.absolute(phi1.oldTime())/runTime.deltaT()
-          + (fvc::interpolate(U1Eqn.H()) & mesh.Sf())
+          + fvc::flux(U1Eqn.H())
          + Vmf*ddtPhi2
          + Kdf*MRF.absolute(phi2)
          - Ff1()
@ -175,7 +175,7 @@ while (pimple.correct())
       *(
            (alphaRhof20 + Vmf)
           *MRF.absolute(phi2.oldTime())/runTime.deltaT()
-          + (fvc::interpolate(U2Eqn.H()) & mesh.Sf())
+          + fvc::flux(U2Eqn.H())
          + Vmf*ddtPhi1
          + Kdf*MRF.absolute(phi1)
          - Ff2()
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
@ -5,7 +5,7 @@
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::flux(HbyA)
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );
    adjustPhi(phiHbyA, U, p_rgh);
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.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
@ -26,6 +26,7 @@ License
 #include "liftModel.H"
 #include "phasePair.H"
 #include "fvcCurl.H"
+#include "fvcFlux.H"
 #include "surfaceInterpolate.H"

 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -78,11 +79,7 @@ Foam::tmp<Foam::volVectorField> Foam::liftModel::F() const

 Foam::tmp<Foam::surfaceScalarField> Foam::liftModel::Ff() const
 {
-    const fvMesh& mesh(this->pair_.phase1().mesh());
-
-    return
-        fvc::interpolate(pair_.dispersed())
-       *(fvc::interpolate(Fi()) & mesh.Sf());
+    return fvc::interpolate(pair_.dispersed())*fvc::flux(Fi());
 }


--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.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
@ -25,6 +25,7 @@ License

 #include "wallLubricationModel.H"
 #include "phasePair.H"
+#include "fvcFlux.H"
 #include "surfaceInterpolate.H"

 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -67,11 +68,7 @@ Foam::tmp<Foam::volVectorField> Foam::wallLubricationModel::F() const

 Foam::tmp<Foam::surfaceScalarField> Foam::wallLubricationModel::Ff() const
 {
-    const fvMesh& mesh(this->pair_.phase1().mesh());
-
-    return
-        fvc::interpolate(pair_.dispersed())
-       *(fvc::interpolate(Fi()) & mesh.Sf());
+    return fvc::interpolate(pair_.dispersed())*fvc::flux(Fi());
 }


--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
@ -70,18 +70,10 @@ tmp<surfaceScalarField> phiF2;
    surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());

    // Phase-1 dispersion, lift and wall-lubrication flux
-    phiF1 =
-    (
-        Df1*snGradAlpha1
-      + (fvc::interpolate(rAU1*F) & mesh.Sf())
-    );
+    phiF1 = Df1*snGradAlpha1 + fvc::flux(rAU1*F);

    // Phase-1 dispersion, lift and wall-lubrication flux
-    phiF2 =
-    (
-      - Df2*snGradAlpha1
-      - (fvc::interpolate(rAU2*F) & mesh.Sf())
-    );
+    phiF2 = - Df2*snGradAlpha1 - fvc::flux(rAU2*F);
 }


@ -172,11 +164,11 @@ while (pimple.correct())
    surfaceScalarField phiHbyA1
    (
        IOobject::groupName("phiHbyA", phase1.name()),
-        (fvc::interpolate(HbyA1) & mesh.Sf())
+        fvc::flux(HbyA1)
      + phiCorrCoeff1*fvc::interpolate(alpha1.oldTime()*rho1.oldTime()*rAU1)
       *(
            MRF.absolute(phi1.oldTime())
-          - (fvc::interpolate(U1.oldTime()) & mesh.Sf())
+          - fvc::flux(U1.oldTime())
        )/runTime.deltaT()
      - phiF1()
      - phig1
@ -186,11 +178,11 @@ while (pimple.correct())
    surfaceScalarField phiHbyA2
    (
        IOobject::groupName("phiHbyA", phase2.name()),
-        (fvc::interpolate(HbyA2) & mesh.Sf())
+        fvc::flux(HbyA2)
      + phiCorrCoeff2*fvc::interpolate(alpha2.oldTime()*rho2.oldTime()*rAU2)
       *(
            MRF.absolute(phi2.oldTime())
-          - (fvc::interpolate(U2.oldTime()) & mesh.Sf())
+          - fvc::flux(U2.oldTime())
        )/runTime.deltaT()
      - phiF2()
      - phig2
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
@ -156,7 +156,7 @@ while (pimple.correct())
       *(
            (alphaRhof10 + Vmf)
           *MRF.absolute(phi1.oldTime())/runTime.deltaT()
-          + (fvc::interpolate(U1Eqn.H()) & mesh.Sf())
+          + fvc::flux(U1Eqn.H())
          + Vmf*ddtPhi2
          + Kdf*MRF.absolute(phi2)
          - Ff1()
@ -174,7 +174,7 @@ while (pimple.correct())
       *(
            (alphaRhof20 + Vmf)
           *MRF.absolute(phi2.oldTime())/runTime.deltaT()
-          + (fvc::interpolate(U2Eqn.H()) & mesh.Sf())
+          + fvc::flux(U2Eqn.H())
          + Vmf*ddtPhi1
          + Kdf*MRF.absolute(phi1)
          - Ff2()
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.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
@ -34,6 +34,7 @@ License
 #include "fixedValueFvPatchFields.H"
 #include "slipFvPatchFields.H"
 #include "partialSlipFvPatchFields.H"
+#include "fvcFlux.H"
 #include "surfaceInterpolate.H"


@ -175,7 +176,7 @@ Foam::phaseModel::phaseModel
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
-                fvc::interpolate(U_) & fluid_.mesh().Sf(),
+                fvc::flux(U_),
                phiTypes
            )
        );