diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index 238a05c696..f2f2af745e 100644
--- 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)
             );
 
diff --git a/applications/solvers/basic/potentialFoam/createFields.H b/applications/solvers/basic/potentialFoam/createFields.H
index ce912510ec..a4e09d481f 100644
--- 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);
 }
 
 
diff --git a/applications/solvers/basic/potentialFoam/potentialFoam.C b/applications/solvers/basic/potentialFoam/potentialFoam.C
index 2e4a569ae1..2de65d10a9 100644
--- 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
diff --git a/applications/solvers/combustion/PDRFoam/pEqn.H b/applications/solvers/combustion/PDRFoam/pEqn.H
index b8d21b7e38..c7308f8b12 100644
--- 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)
         )
     );
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index 39830d0135..afba129157 100644
--- 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)
         )
     );
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index d84b9bd154..aca8c2c11a 100644
--- 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)
         )
     );
diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index 11f140577f..7b5249d57e 100644
--- 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
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index f5b7141224..ac7107acf0 100644
--- 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)
         )
     );
diff --git a/applications/solvers/combustion/reactingFoam/pcEqn.H b/applications/solvers/combustion/reactingFoam/pcEqn.H
index bf8eb899af..713f443fc5 100644
--- 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)
         )
     );
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
index ad01cd5f83..926e69f7c1 100644
--- 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
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
index 3d02c63381..eff014cf43 100644
--- 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)
             )
         );
diff --git a/applications/solvers/compressible/rhoCentralFoam/createFields.H b/applications/solvers/compressible/rhoCentralFoam/createFields.H
index 632efaf9d1..538e936e39 100644
--- 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
diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
index a9fa5559b8..3b43b5dcec 100644
--- 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)
         );
diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
index 841c9a5162..3754632d3d 100644
--- 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
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index f5b7141224..ac7107acf0 100644
--- 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)
         )
     );
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index bf8eb899af..713f443fc5 100644
--- 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)
         )
     );
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index 4abc5cd4aa..7bd540df40 100644
--- 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)
     );
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H b/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
index ba2db0d04a..0a38e05ee9 100644
--- 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
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index b085a2b9d1..e46f2a6691 100644
--- 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);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
index 51b3865165..67b962f805 100644
--- 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);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H
index b586cd15c0..7b55983161 100644
--- 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);
diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H
index f4f6d5ae2a..bd553928b1 100644
--- 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)
     )
 );
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
index 283e450e00..92e0bb3673 100644
--- 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)
     )
 );
diff --git a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
index c60a331b2d..eb147db6c0 100644
--- 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)
                     )
                 );
diff --git a/applications/solvers/electromagnetics/mhdFoam/createPhiB.H b/applications/solvers/electromagnetics/mhdFoam/createPhiB.H
index 9b8aea9ae2..de0c573259 100644
--- 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)
         );
     }
 
diff --git a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
index 0d21312a30..43e800b48b 100644
--- 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())
             {
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
index a9be0e59f7..51d9461c87 100644
--- 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
     );
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
index 7df0ac46d2..f29075b149 100644
--- 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);
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index 66864139fa..bb371bb9f3 100644
--- 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
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index a744b9d19f..6c1ae382fa 100644
--- 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);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index 3fdadc042a..7fe544d6a5 100644
--- 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);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index 358535efe9..a8a0807f24 100644
--- 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
diff --git a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
index a303ae4fe4..c062af0364 100644
--- 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
diff --git a/applications/solvers/incompressible/icoFoam/icoFoam.C b/applications/solvers/incompressible/icoFoam/icoFoam.C
index e3a7fb67c5..04581bcb20 100644
--- 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)
             );
 
diff --git a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
index 3edbedf39f..bdac3f8336 100644
--- 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)
             );
 
diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
index aaeb12ae06..3cb0580ba4 100644
--- 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)
 );
 
diff --git a/applications/solvers/incompressible/pimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pEqn.H
index ffce137326..c8d02ed281 100644
--- 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)
 );
 
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index 452f13fb1f..a53391a90b 100644
--- 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)
 );
 
diff --git a/applications/solvers/incompressible/pisoFoam/pEqn.H b/applications/solvers/incompressible/pisoFoam/pEqn.H
index 7d3284e0d4..05c34376ec 100644
--- 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)
 );
 
diff --git a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
index 47ec8ac9f6..6575a09f05 100644
--- 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
                 );
diff --git a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
index 9982ddf546..e92838716c 100644
--- 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);
diff --git a/applications/solvers/incompressible/simpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/pEqn.H
index 5d5dde848c..7e394e9c10 100644
--- 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);
 
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
index 68bdcc77e0..c76143714e 100644
--- 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);
 
diff --git a/applications/solvers/lagrangian/DPMFoam/UcEqn.H b/applications/solvers/lagrangian/DPMFoam/UcEqn.H
index 0bd4a7146d..d31bead5e1 100644
--- 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())
diff --git a/applications/solvers/lagrangian/DPMFoam/pEqn.H b/applications/solvers/lagrangian/DPMFoam/pEqn.H
index 6b7084de63..9e46551122 100644
--- 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
         )
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index ce9f819646..4f6eaf3550 100644
--- 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)
         )
     );
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
index 11f140577f..7b5249d57e 100644
--- 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
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index cb5b083b88..f6c8ceec21 100644
--- 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)
         )
     );
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
index 9b48bb939d..055eff6f0b 100644
--- 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);
diff --git a/applications/solvers/lagrangian/sprayFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/pEqn.H
index 1447f5b6be..b15fb81e13 100644
--- 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)
         )
     );
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
index 8522e60918..90033814d3 100644
--- 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)
     );
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index 37892c0e73..51a2da33bd 100644
--- 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);
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index 5b238ccefa..4c337c4269 100644
--- 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));
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index c326215512..859fc9cc4e 100644
--- 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)
     );
 
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index a0a91c3255..a4a3fb2e29 100644
--- 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)
     );
 
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
index 26ffef825a..70ac37c1b8 100644
--- 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)
     );
 
diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
index b586e77fad..e77f66b880 100644
--- 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)
     {
diff --git a/applications/solvers/multiphase/driftFluxFoam/pEqn.H b/applications/solvers/multiphase/driftFluxFoam/pEqn.H
index 6cc93506df..6797c44ccb 100644
--- 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);
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
index e6723f47b2..2c7df14ae7 100644
--- 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);
diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H
index 42cf70fce5..43edce0c5d 100644
--- 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);
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
index f28bb8c2aa..7501f58062 100644
--- 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)
     );
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index fd85d96714..0b0a31c669 100644
--- 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);
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.C
index 9dc9f2bf80..2f5f833d20 100644
--- 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
             )
         );
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index d2616fc556..4cddb12921 100644
--- 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]);
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
index 7d4c5f9039..76a44bd134 100644
--- 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)
 );
 
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H
index a74f6c2664..dac95a1974 100644
--- 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)
     );
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C
index 79e52a2856..7c0bf86c94 100644
--- 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());
 }
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
index c17b649212..a2a95cce05 100644
--- 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());
 }
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
index 74c6b7b46a..de31cefaa3 100644
--- 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
             )
         );
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
index 42e0fe70e4..05a049b4c6 100644
--- 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]
             )
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
index 104bc0047f..ba312ab248 100644
--- 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
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index d618226411..444f2de81f 100644
--- 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()
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
index 870034570a..d801fb8fa5 100644
--- 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);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C
index 79e52a2856..7c0bf86c94 100644
--- 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());
 }
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
index 2f3cb1aaa0..d148d2cefd 100644
--- 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());
 }
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
index c2621ae66a..b418877743 100644
--- 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
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
index 63321bab7d..6855b6fb96 100644
--- 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()
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
index b0b5a591b3..a553a01b5e 100644
--- 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
             )
         );
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 97295969b9..87996b789d 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -390,6 +390,7 @@ laplacianSchemes = finiteVolume/laplacianSchemes
 $(laplacianSchemes)/laplacianScheme/laplacianSchemes.C
 $(laplacianSchemes)/gaussLaplacianScheme/gaussLaplacianSchemes.C
 
+finiteVolume/fvc/fvcFlux.C
 finiteVolume/fvc/fvcMeshPhi.C
 finiteVolume/fvc/fvcSmooth/fvcSmooth.C
 finiteVolume/fvc/fvcReconstructMag.C
diff --git a/src/finiteVolume/cfdTools/incompressible/createPhi.H b/src/finiteVolume/cfdTools/incompressible/createPhi.H
index c0c5842319..f2d8b4f9fb 100644
--- a/src/finiteVolume/cfdTools/incompressible/createPhi.H
+++ b/src/finiteVolume/cfdTools/incompressible/createPhi.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,7 +46,7 @@ surfaceScalarField phi
         IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
-    linearInterpolate(U) & mesh.Sf()
+    fvc::flux(U)
 );
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/finiteVolume/cfdTools/incompressible/createPhiv.H b/src/finiteVolume/cfdTools/incompressible/createPhiv.H
index 6a723b29fe..97189c353f 100644
--- a/src/finiteVolume/cfdTools/incompressible/createPhiv.H
+++ b/src/finiteVolume/cfdTools/incompressible/createPhiv.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,7 +46,7 @@ surfaceScalarField phiv
         IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
-    linearInterpolate(U) & mesh.Sf()
+    flux(U)
 );
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
index aabc1320f3..4165b62191 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
@@ -600,9 +600,10 @@ CoEulerDdtScheme<Type>::fvcDdtUfCorr
 {
     const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
 
+    fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
     fluxFieldType phiCorr
     (
-        mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+        phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
     );
 
     return tmp<fluxFieldType>
@@ -615,12 +616,7 @@ CoEulerDdtScheme<Type>::fvcDdtUfCorr
                 mesh().time().timeName(),
                 mesh()
             ),
-            this->fvcDdtPhiCoeff
-            (
-                U.oldTime(),
-                (mesh().Sf() & Uf.oldTime()),
-                phiCorr
-            )
+            this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
            *rDeltaT*phiCorr
         )
     );
@@ -639,7 +635,7 @@ CoEulerDdtScheme<Type>::fvcDdtPhiCorr
 
     fluxFieldType phiCorr
     (
-        phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+        phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
     );
 
     return tmp<fluxFieldType>
@@ -681,10 +677,8 @@ CoEulerDdtScheme<Type>::fvcDdtUfCorr
             rho.oldTime()*U.oldTime()
         );
 
-        fluxFieldType phiCorr
-        (
-            mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
-        );
+        fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
+        fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
 
         return tmp<fluxFieldType>
         (
@@ -697,13 +691,7 @@ CoEulerDdtScheme<Type>::fvcDdtUfCorr
                     mesh().time().timeName(),
                     mesh()
                 ),
-                this->fvcDdtPhiCoeff
-                (
-                    rhoU0,
-                    mesh().Sf() & Uf.oldTime(),
-                    phiCorr
-                )
-               *rDeltaT*phiCorr
+                this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr)*rDeltaT*phiCorr
             )
         );
     }
@@ -750,7 +738,7 @@ CoEulerDdtScheme<Type>::fvcDdtPhiCorr
 
         fluxFieldType phiCorr
         (
-            phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+            phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
         );
 
         return tmp<fluxFieldType>
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
index 9d4be317ff..67473222f5 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
@@ -1227,9 +1227,10 @@ CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
             this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
            *(
                 (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
-              - (
-                    mesh().Sf()
-                  & fvc::interpolate(rDtCoef*U.oldTime() + offCentre_(ddt0()))
+              - fvc::dotInterpolate
+                (
+                    mesh().Sf(),
+                    rDtCoef*U.oldTime() + offCentre_(ddt0())
                 )
             )
         )
@@ -1398,9 +1399,10 @@ CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
                 this->fvcDdtPhiCoeff(rhoU0, phi.oldTime())
                *(
                     (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
-                  - (
-                        mesh().Sf()
-                      & fvc::interpolate(rDtCoef*rhoU0 + offCentre_(ddt0()))
+                  - fvc::dotInterpolate
+                    (
+                        mesh().Sf(),
+                        rDtCoef*rhoU0 + offCentre_(ddt0())
                     )
                 )
             )
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
index e6c7239f61..5a72c9d932 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
@@ -489,9 +489,10 @@ EulerDdtScheme<Type>::fvcDdtUfCorr
 {
     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
 
+    fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
     fluxFieldType phiCorr
     (
-        mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+        phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
     );
 
     return tmp<fluxFieldType>
@@ -504,12 +505,7 @@ EulerDdtScheme<Type>::fvcDdtUfCorr
                 mesh().time().timeName(),
                 mesh()
             ),
-            this->fvcDdtPhiCoeff
-            (
-                U.oldTime(),
-                mesh().Sf() & Uf.oldTime(),
-                phiCorr
-            )
+            this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
            *rDeltaT*phiCorr
         )
     );
@@ -528,7 +524,7 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr
 
     fluxFieldType phiCorr
     (
-        phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+        phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
     );
 
     return tmp<fluxFieldType>
@@ -570,10 +566,8 @@ EulerDdtScheme<Type>::fvcDdtUfCorr
             rho.oldTime()*U.oldTime()
         );
 
-        fluxFieldType phiCorr
-        (
-            mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
-        );
+        fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
+        fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
 
         return tmp<fluxFieldType>
         (
@@ -586,13 +580,7 @@ EulerDdtScheme<Type>::fvcDdtUfCorr
                     mesh().time().timeName(),
                     mesh()
                 ),
-                this->fvcDdtPhiCoeff
-                (
-                    rhoU0,
-                    mesh().Sf() & Uf.oldTime(),
-                    phiCorr
-                )
-               *rDeltaT*phiCorr
+                this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr)*rDeltaT*phiCorr
             )
         );
     }
@@ -639,7 +627,7 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr
 
         fluxFieldType phiCorr
         (
-            phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+            phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
         );
 
         return tmp<fluxFieldType>
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
index a3b9b841e1..3994877922 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
@@ -599,9 +599,10 @@ SLTSDdtScheme<Type>::fvcDdtUfCorr
 {
     const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
 
+    fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
     fluxFieldType phiCorr
     (
-        mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+        phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
     );
 
     return tmp<fluxFieldType>
@@ -614,12 +615,7 @@ SLTSDdtScheme<Type>::fvcDdtUfCorr
                 mesh().time().timeName(),
                 mesh()
             ),
-            this->fvcDdtPhiCoeff
-            (
-                U.oldTime(),
-                (mesh().Sf() & Uf.oldTime()),
-                phiCorr
-            )
+            this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
            *rDeltaT*phiCorr
         )
     );
@@ -638,7 +634,7 @@ SLTSDdtScheme<Type>::fvcDdtPhiCorr
 
     fluxFieldType phiCorr
     (
-        phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+        phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
     );
 
     return tmp<fluxFieldType>
@@ -680,10 +676,8 @@ SLTSDdtScheme<Type>::fvcDdtUfCorr
             rho.oldTime()*U.oldTime()
         );
 
-        fluxFieldType phiCorr
-        (
-            mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
-        );
+        fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
+        fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
 
         return tmp<fluxFieldType>
         (
@@ -696,13 +690,7 @@ SLTSDdtScheme<Type>::fvcDdtUfCorr
                     mesh().time().timeName(),
                     mesh()
                 ),
-                this->fvcDdtPhiCoeff
-                (
-                    rhoU0,
-                    mesh().Sf() & Uf.oldTime(),
-                    phiCorr
-                )
-               *rDeltaT*phiCorr
+                this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr)*rDeltaT*phiCorr
             )
         );
     }
@@ -749,7 +737,7 @@ SLTSDdtScheme<Type>::fvcDdtPhiCorr
 
         fluxFieldType phiCorr
         (
-            phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+            phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
         );
 
         return tmp<fluxFieldType>
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
index e3127a453b..a405a2c3f3 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
@@ -748,12 +748,10 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
            *rDeltaT
            *(
                 (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
-              - (
-                    mesh().Sf()
-                  & fvc::interpolate
-                    (
-                        coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
-                    )
+              - fvc::dotInterpolate
+                (
+                    mesh().Sf(),
+                    coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
                 )
             )
         )
@@ -886,9 +884,10 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
                *rDeltaT
                *(
                     (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
-                  - (
-                        mesh().Sf()
-                      & fvc::interpolate(coefft0*rhoU0 - coefft00*rhoU00)
+                  - fvc::dotInterpolate
+                    (
+                        mesh().Sf(),
+                        coefft0*rhoU0 - coefft00*rhoU00
                     )
                 )
             )
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C
index 62cca7ec3a..cd0e1adc7e 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C
@@ -183,7 +183,7 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
     const fluxFieldType& phi
 )
 {
-    return fvcDdtPhiCoeff(U, phi, phi - (mesh().Sf() & fvc::interpolate(U)));
+    return fvcDdtPhiCoeff(U, phi, phi - fvc::dotInterpolate(mesh().Sf(), U));
 }
 
 
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H
index 2b36031b92..f475ed24f4 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H
@@ -288,7 +288,7 @@ tmp<surfaceScalarField> SS<scalar>::fvcDdtUfCorr                               \
     const surfaceScalarField& Uf                                               \
 )                                                                              \
 {                                                                              \
-    NotImplemented;                               \
+    NotImplemented;                                                            \
     return surfaceScalarField::null();                                         \
 }                                                                              \
                                                                                \
@@ -299,7 +299,7 @@ tmp<surfaceScalarField> SS<scalar>::fvcDdtPhiCorr                              \
     const surfaceScalarField& phi                                              \
 )                                                                              \
 {                                                                              \
-    NotImplemented;                              \
+    NotImplemented;                                                            \
     return surfaceScalarField::null();                                         \
 }                                                                              \
                                                                                \
@@ -311,7 +311,7 @@ tmp<surfaceScalarField> SS<scalar>::fvcDdtUfCorr                               \
     const surfaceScalarField& Uf                                               \
 )                                                                              \
 {                                                                              \
-    NotImplemented;                               \
+    NotImplemented;                                                            \
     return surfaceScalarField::null();                                         \
 }                                                                              \
                                                                                \
@@ -323,7 +323,7 @@ tmp<surfaceScalarField> SS<scalar>::fvcDdtPhiCorr                              \
     const surfaceScalarField& phi                                              \
 )                                                                              \
 {                                                                              \
-    NotImplemented;                              \
+    NotImplemented;                                                            \
     return surfaceScalarField::null();                                         \
 }                                                                              \
                                                                                \
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
index 69d707974a..da4c454a50 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
@@ -496,9 +496,10 @@ localEulerDdtScheme<Type>::fvcDdtUfCorr
 {
     const surfaceScalarField rDeltaT(fvc::interpolate(localRDeltaT()));
 
+    fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
     fluxFieldType phiCorr
     (
-        mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+        phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
     );
 
     return tmp<fluxFieldType>
@@ -511,12 +512,7 @@ localEulerDdtScheme<Type>::fvcDdtUfCorr
                 mesh().time().timeName(),
                 mesh()
             ),
-            this->fvcDdtPhiCoeff
-            (
-                U.oldTime(),
-                (mesh().Sf() & Uf.oldTime()),
-                phiCorr
-            )
+            this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
            *rDeltaT*phiCorr
         )
     );
@@ -535,7 +531,7 @@ localEulerDdtScheme<Type>::fvcDdtPhiCorr
 
     fluxFieldType phiCorr
     (
-        phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+        phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
     );
 
     return tmp<fluxFieldType>
@@ -577,10 +573,8 @@ localEulerDdtScheme<Type>::fvcDdtUfCorr
             rho.oldTime()*U.oldTime()
         );
 
-        fluxFieldType phiCorr
-        (
-            mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
-        );
+        fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
+        fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
 
         return tmp<fluxFieldType>
         (
@@ -593,13 +587,7 @@ localEulerDdtScheme<Type>::fvcDdtUfCorr
                     mesh().time().timeName(),
                     mesh()
                 ),
-                this->fvcDdtPhiCoeff
-                (
-                    rhoU0,
-                    mesh().Sf() & Uf.oldTime(),
-                    phiCorr
-                )
-               *rDeltaT*phiCorr
+                this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr)*rDeltaT*phiCorr
             )
         );
     }
@@ -646,7 +634,7 @@ localEulerDdtScheme<Type>::fvcDdtPhiCorr
 
         fluxFieldType phiCorr
         (
-            phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+            phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
         );
 
         return tmp<fluxFieldType>
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcFlux.C b/src/finiteVolume/finiteVolume/fvc/fvcFlux.C
index ab9c346e3f..b1bf586e66 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcFlux.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcFlux.C
@@ -24,167 +24,32 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "fvcFlux.H"
-#include "fvMesh.H"
-#include "convectionScheme.H"
+#include "surfaceInterpolate.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fvc
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
-flux
+Foam::tmp<Foam::surfaceScalarField> Foam::fvc::flux
 (
-    const surfaceScalarField& phi,
-    const GeometricField<Type, fvPatchField, volMesh>& vf,
-    const word& name
+    const volVectorField& vvf
 )
 {
-    return fv::convectionScheme<Type>::New
+    return scheme<vector>
     (
-        vf.mesh(),
-        phi,
-        vf.mesh().divScheme(name)
-    )().flux(phi, vf);
+        vvf.mesh(),
+        "flux(" + vvf.name() + ')'
+    )().dotInterpolate(vvf.mesh().Sf(), vvf);
 }
 
 
-template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
-flux
+Foam::tmp<Foam::surfaceScalarField> Foam::fvc::flux
 (
-    const tmp<surfaceScalarField>& tphi,
-    const GeometricField<Type, fvPatchField, volMesh>& vf,
-    const word& name
+    const tmp<volVectorField>& tvvf
 )
 {
-    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
-    (
-        fvc::flux(tphi(), vf, name)
-    );
-    tphi.clear();
+    tmp<surfaceScalarField> Flux(fvc::flux(tvvf()));
+    tvvf.clear();
     return Flux;
 }
 
 
-template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
-flux
-(
-    const surfaceScalarField& phi,
-    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf,
-    const word& name
-)
-{
-    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
-    (
-        fvc::flux(phi, tvf(), name)
-    );
-    tvf.clear();
-    return Flux;
-}
-
-
-template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
-flux
-(
-    const tmp<surfaceScalarField>& tphi,
-    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf,
-    const word& name
-)
-{
-    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
-    (
-        fvc::flux(tphi(), tvf(), name)
-    );
-    tphi.clear();
-    tvf.clear();
-    return Flux;
-}
-
-
-template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
-flux
-(
-    const surfaceScalarField& phi,
-    const GeometricField<Type, fvPatchField, volMesh>& vf
-)
-{
-    return fvc::flux
-    (
-        phi, vf, "flux("+phi.name()+','+vf.name()+')'
-    );
-}
-
-
-template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
-flux
-(
-    const tmp<surfaceScalarField>& tphi,
-    const GeometricField<Type, fvPatchField, volMesh>& vf
-)
-{
-    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
-    (
-        fvc::flux(tphi(), vf)
-    );
-    tphi.clear();
-    return Flux;
-}
-
-
-template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
-flux
-(
-    const surfaceScalarField& phi,
-    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf
-)
-{
-    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
-    (
-        fvc::flux(phi, tvf())
-    );
-    tvf.clear();
-    return Flux;
-}
-
-
-template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
-flux
-(
-    const tmp<surfaceScalarField>& tphi,
-    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf
-)
-{
-    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
-    (
-        fvc::flux(tphi(), tvf())
-    );
-    tphi.clear();
-    tvf.clear();
-    return Flux;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fvc
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcFlux.H b/src/finiteVolume/finiteVolume/fvc/fvcFlux.H
index 0e3e3a9fe6..e62dcec0cb 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcFlux.H
+++ b/src/finiteVolume/finiteVolume/fvc/fvcFlux.H
@@ -50,6 +50,19 @@ namespace Foam
 
 namespace fvc
 {
+    //- Return the face-flux field obtained from the given volVectorField
+    tmp<surfaceScalarField> flux(const volVectorField& vvf);
+
+    //- Return the face-flux field obtained from the given tmp volVectorField
+    tmp<surfaceScalarField> flux(const tmp<volVectorField>& tvvf);
+
+    template<class Type>
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> flux
+    (
+        const surfaceScalarField&,
+        const tmp<GeometricField<Type, fvPatchField, volMesh>>&
+    );
+
     template<class Type>
     tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> flux
     (
@@ -120,7 +133,7 @@ namespace fvc
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "fvcFlux.C"
+    #include "fvcFluxTemplates.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcFluxTemplates.C b/src/finiteVolume/finiteVolume/fvc/fvcFluxTemplates.C
new file mode 100644
index 0000000000..ab9c346e3f
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/fvc/fvcFluxTemplates.C
@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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 "fvcFlux.H"
+#include "fvMesh.H"
+#include "convectionScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fvc
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+flux
+(
+    const surfaceScalarField& phi,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
+    const word& name
+)
+{
+    return fv::convectionScheme<Type>::New
+    (
+        vf.mesh(),
+        phi,
+        vf.mesh().divScheme(name)
+    )().flux(phi, vf);
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+flux
+(
+    const tmp<surfaceScalarField>& tphi,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
+    const word& name
+)
+{
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
+    (
+        fvc::flux(tphi(), vf, name)
+    );
+    tphi.clear();
+    return Flux;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+flux
+(
+    const surfaceScalarField& phi,
+    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf,
+    const word& name
+)
+{
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
+    (
+        fvc::flux(phi, tvf(), name)
+    );
+    tvf.clear();
+    return Flux;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+flux
+(
+    const tmp<surfaceScalarField>& tphi,
+    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf,
+    const word& name
+)
+{
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
+    (
+        fvc::flux(tphi(), tvf(), name)
+    );
+    tphi.clear();
+    tvf.clear();
+    return Flux;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+flux
+(
+    const surfaceScalarField& phi,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    return fvc::flux
+    (
+        phi, vf, "flux("+phi.name()+','+vf.name()+')'
+    );
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+flux
+(
+    const tmp<surfaceScalarField>& tphi,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
+    (
+        fvc::flux(tphi(), vf)
+    );
+    tphi.clear();
+    return Flux;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+flux
+(
+    const surfaceScalarField& phi,
+    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf
+)
+{
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
+    (
+        fvc::flux(phi, tvf())
+    );
+    tvf.clear();
+    return Flux;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+flux
+(
+    const tmp<surfaceScalarField>& tphi,
+    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf
+)
+{
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
+    (
+        fvc::flux(tphi(), tvf())
+    );
+    tphi.clear();
+    tvf.clear();
+    return Flux;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fvc
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C b/src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C
index f7707e0db7..c7c38d79e8 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C
@@ -26,6 +26,7 @@ License
 #include "fvcMeshPhi.H"
 #include "fvMesh.H"
 #include "ddtScheme.H"
+#include "surfaceInterpolate.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C
index fd9df516e8..937a359d85 100644
--- a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C
@@ -121,7 +121,7 @@ gaussLaplacianScheme<Type, GType>::gammaSnGradCorr
         tgammaSnGradCorr.ref().replace
         (
             cmpt,
-            SfGammaCorr & fvc::interpolate(fvc::grad(vf.component(cmpt)))
+            fvc::dotInterpolate(SfGammaCorr, fvc::grad(vf.component(cmpt)))
         );
     }
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
index cb23171ea6..6a2af46f66 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
@@ -308,4 +308,66 @@ Foam::fvc::interpolate
 }
 
 
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::innerProduct<Foam::vector, Type>::type,
+        Foam::fvsPatchField,
+        Foam::surfaceMesh
+    >
+>
+Foam::fvc::dotInterpolate
+(
+    const surfaceVectorField& Sf,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    if (surfaceInterpolation::debug)
+    {
+        InfoInFunction
+            << "interpolating GeometricField<Type, fvPatchField, volMesh> "
+            << vf.name() << " using run-time selected scheme"
+            << endl;
+    }
+
+    return scheme<Type>
+    (
+        vf.mesh(),
+        "dotInterpolate(" + Sf.name() + ',' + vf.name() + ')'
+    )().dotInterpolate(Sf, vf);
+}
+
+
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::innerProduct<Foam::vector, Type>::type,
+        Foam::fvsPatchField,
+        Foam::surfaceMesh
+    >
+>
+Foam::fvc::dotInterpolate
+(
+    const surfaceVectorField& Sf,
+    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf
+)
+{
+    tmp
+    <
+        GeometricField
+        <
+            typename Foam::innerProduct<Foam::vector, Type>::type,
+            fvsPatchField,
+            surfaceMesh
+        >
+    > tsf = dotInterpolate(Sf, tvf());
+    tvf.clear();
+    return tsf;
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.H b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.H
index 389e72018a..185284d21e 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.H
@@ -157,6 +157,13 @@ namespace fvc
     );
 
 
+    //- Interpolate field onto faces using 'interpolate(\<name\>)'
+    template<class Type>
+    static tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate
+    (
+        const GeometricField<Type, fvPatchField, volMesh>& tvf
+    );
+
     //- Interpolate tmp field onto faces using 'interpolate(\<name\>)'
     template<class Type>
     static tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate
@@ -164,12 +171,6 @@ namespace fvc
         const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf
     );
 
-    //- Interpolate field onto faces using 'interpolate(\<name\>)'
-    template<class Type>
-    static tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate
-    (
-        const GeometricField<Type, fvPatchField, volMesh>& tvf
-    );
 
     //- Interpolate boundary field onto faces (simply a type conversion)
     template<class Type>
@@ -190,6 +191,43 @@ namespace fvc
     {
         return one();
     }
+
+
+    //- Interpolate field onto faces
+    //  and 'dot' with given surfaceVectorField Sf
+    template<class Type>
+    static
+    tmp
+    <
+        GeometricField
+        <
+            typename innerProduct<vector, Type>::type,
+            fvsPatchField,
+            surfaceMesh
+            >
+    > dotInterpolate
+    (
+        const surfaceVectorField& Sf,
+        const GeometricField<Type, fvPatchField, volMesh>& tvf
+    );
+
+    //- Interpolate tmp field onto faces
+    //  and 'dot' with given surfaceVectorField Sf
+    template<class Type>
+    static
+    tmp
+    <
+        GeometricField
+        <
+            typename innerProduct<vector, Type>::type,
+            fvsPatchField,
+            surfaceMesh
+            >
+    > dotInterpolate
+    (
+        const surfaceVectorField& Sf,
+        const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf
+    );
 }
 
 
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
index ba28c423c9..1d4d885be4 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
@@ -140,9 +140,10 @@ void Foam::displacementSBRStressFvMotionSolver::solve()
         (
             Df
            *(
+               fvc::dotInterpolate
                (
-                   cellDisplacement_.mesh().Sf()
-                 & fvc::interpolate(gradCd.T() - gradCd)
+                   cellDisplacement_.mesh().Sf(),
+                   gradCd.T() - gradCd
                )
 
                // Solid-body rotation "lambda" term
@@ -162,9 +163,10 @@ void Foam::displacementSBRStressFvMotionSolver::solve()
         (
             Df
            *(
+               fvc::dotInterpolate
                (
-                   cellDisplacement_.mesh().Sf()
-                 & fvc::interpolate(gradCd + gradCd.T())
+                   cellDisplacement_.mesh().Sf(),
+                   gradCd + gradCd.T()
                )
 
                // Solid-body rotation "lambda" term
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
index 529b95f964..ae5357b185 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
@@ -30,6 +30,7 @@ License
 #include "fvcSnGrad.H"
 #include "fvcReconstruct.H"
 #include "fvcVolumeIntegrate.H"
+#include "fvcFlux.H"
 #include "addToRunTimeSelectionTable.H"
 #include "mappedWallPolyPatch.H"
 #include "mapDistribute.H"
@@ -325,7 +326,7 @@ tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
                       + fvc::snGrad(pp, "snGrad(p)")*fvc::interpolate(delta_)
                       + fvc::snGrad(delta_)*fvc::interpolate(pp)
                     )
-                  - (fvc::interpolate(rho_*gTan()) & regionMesh().Sf())
+                  - fvc::flux(rho_*gTan())
                 )
             )
         );
@@ -365,15 +366,14 @@ void kinematicSingleLayer::solveThickness
             fvc::snGrad(pu, "snGrad(p)")
           + fvc::snGrad(pp, "snGrad(p)")*fvc::interpolate(delta_)
         )
-      - (fvc::interpolate(rho_*gTan()) & regionMesh().Sf())
+      - fvc::flux(rho_*gTan())
     );
     constrainFilmField(phiAdd, 0.0);
 
     surfaceScalarField phid
     (
         "phid",
-        (fvc::interpolate(U_*rho_) & regionMesh().Sf())
-      - deltarUAf*phiAdd*rhof
+        fvc::flux(U_*rho_) - deltarUAf*phiAdd*rhof
     );
     constrainFilmField(phid, 0.0);
 
@@ -807,7 +807,7 @@ kinematicSingleLayer::kinematicSingleLayer
                 IOobject::AUTO_WRITE,
                 false
             ),
-            fvc::interpolate(deltaRho_*U_) & regionMesh().Sf()
+            fvc::flux(deltaRho_*U_)
         );
 
         phi_ == phi0;
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
index 4ffc27c8d9..e1c65f2ad3 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
@@ -26,6 +26,7 @@ License
 #include "thermoSingleLayer.H"
 #include "fvcDiv.H"
 #include "fvcLaplacian.H"
+#include "fvcFlux.H"
 #include "fvm.H"
 #include "addToRunTimeSelectionTable.H"
 #include "zeroGradientFvPatchFields.H"
@@ -556,7 +557,7 @@ thermoSingleLayer::thermoSingleLayer
                 IOobject::AUTO_WRITE,
                 false
             ),
-            fvc::interpolate(deltaRho_*U_) & regionMesh().Sf()
+            fvc::flux(deltaRho_*U_)
         );
 
         phi_ == phi0;
diff --git a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/system/fvSchemes b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/system/fvSchemes
index 309cee03eb..12c112d4e0 100644
--- a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/system/fvSchemes
+++ b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/system/fvSchemes
@@ -37,7 +37,7 @@ laplacianSchemes
 
 interpolationSchemes
 {
-    default         none;
+    default         linear;
 }
 
 // ************************************************************************* //
diff --git a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/system/fvSchemes b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/system/fvSchemes
index 309cee03eb..12c112d4e0 100644
--- a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/system/fvSchemes
+++ b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/system/fvSchemes
@@ -37,7 +37,7 @@ laplacianSchemes
 
 interpolationSchemes
 {
-    default         none;
+    default         linear;
 }
 
 // ************************************************************************* //