diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index 238a05c6960767463e9f7a37a7e7c0b6b9e3a6c4..f2f2af745e250c84a87b0d6c10f4e28c15e7b3df 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 ce912510ecbf992542d9ba8d3206ec4bb199d6f0..a4e09d481fb2e4c0c693eb6fb9b423efe233ec96 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 2e4a569ae1a157ec2c2bde2a915ddfd725b2da6b..2de65d10a9e9abcb6999bfc36a3365374188847e 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 b8d21b7e38ac8a75911c4e8820c71fb163821925..c7308f8b122d90a98222f20ea6b7c9403c9cac0a 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 39830d0135193e30937e5c56d0fc4a096c2b81d1..afba129157eabffbd7a24e257dda287d533721a3 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 d84b9bd1542aa8948baebab8e3eb9027846c54f2..aca8c2c11a3660eca549c64d4cafbf7eba14f878 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 11f140577faea67a42ae5658936052c59841b30d..7b5249d57e7bd2ce99db1f702aa863eb56a1adf2 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 f5b71412247bef3a384e7e7f7bc0713268305c0a..ac7107acf066722b286275569b0a8b2b39af3e05 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 bf8eb899aff9b436e5a72bb976d0868d1815f531..713f443fc5dbed1e31266ae52d845518954f36e9 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 ad01cd5f832473a1f8c57478940c9d857c7e25ba..926e69f7c1f7933c402e0289b2e7b0f22cc25bb7 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 3d02c63381e0d423fb4a9c036f4f4885a3b3f344..eff014cf43d847917101add0d0be78cf34b64628 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 632efaf9d15bbae97bde1c02683936a11da7521d..538e936e39f11844ea732f489e41a08fbd188d8d 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 a9fa5559b8b7219404f125644d0e2de63f00e3ba..3b43b5dcec39083ca06b08ecc07871b5c418924a 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 841c9a516202e7f9cf0d747f63a55fb56e9adf9c..3754632d3d4be9724dbc7544c96497c2acf68889 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 f5b71412247bef3a384e7e7f7bc0713268305c0a..ac7107acf066722b286275569b0a8b2b39af3e05 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 bf8eb899aff9b436e5a72bb976d0868d1815f531..713f443fc5dbed1e31266ae52d845518954f36e9 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 4abc5cd4aadc7639da4f6ec347091db20109c2c2..7bd540df40d7e4a832824ddf21e4e6c29f312a50 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 ba2db0d04a42fd01a1275f3b03135a2f7274c179..0a38e05ee9fe6972c1b89aaadfea37e24fc25fd7 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 b085a2b9d1c3077e214600758ca9fbbb8911d0ec..e46f2a6691322c46462cb42961fcbd19ab1c8cc7 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 51b3865165718a89eba453ddb781ee2bf5991bb1..67b962f80509f800763f0676e21f7fecfcf0e15c 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 b586cd15c0d0715449b2effff47c58107483bc04..7b559831614b64ecec5b82ca1fb40b85a1d4e309 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 f4f6d5ae2ae0f47e2adc5d864fdfbcec6311eeed..bd553928b16db9d153091432a2ec9f8c046e0443 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 283e450e00aade4d7fb55d07396f622cf034dbfa..92e0bb3673d63f5b379d2b00f07f7f2aef768b5b 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 c60a331b2d6e117c8873cf74765c6c4bdbb7cf9a..eb147db6c020ff6d96d14f10d7c88919af28a63a 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 9b8aea9ae20361e6fbaa9bf805de05fbeed9e7cd..de0c573259c80f9d301b22dd316b242a05055d19 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 0d21312a300e193c87f3d353ced4c2914d182335..43e800b48b755783175d90fe083bd7ff3bf4fee1 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 a9be0e59f7b3aba8e2eccfa14afa23b046368d21..51d9461c87c7c76797c51198594a4b787de8d3d5 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 7df0ac46d290683e2a53a8d085bc1dad26fad5c0..f29075b14913a955fa84a47f912d072e9d079c05 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 66864139fa866775165687bd36c939502889bbc0..bb371bb9f3f548c01c3c197fca69140a9f46ccc8 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 a744b9d19fbc52fc502dba3ac275c6d27a8ab7e4..6c1ae382fae7ca1ef62022d2ca4eab47880f6da0 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 3fdadc042a2684935f04b5563e8a87a54181d352..7fe544d6a57afe68850f81d39483f79c02f10660 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 358535efe97f9ed21a6dab5ba5cf6675f80c63f6..a8a0807f24da01df06e9e8c390ad93dd7dffd448 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 a303ae4fe4bf4e2b9a512faaab16b79fc8f57496..c062af036458a51fcfbe486ebbf12efb6dcc876d 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 e3a7fb67c57d0d40a73dbb3ae3d95665b184e21d..04581bcb20682f202bdb9c0db4408c1ecfbd892c 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 3edbedf39ff0c9d230e44ce530e161d1fc9d30fb..bdac3f83361b2c353ffd668c3b2681c9bbf60f5f 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 aaeb12ae0671f60aa9f28c1b7347b666252ce873..3cb0580ba465828c513bfb1fb07a4b6af7e9335b 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 ffce1373263fc703f13b32a16e742c5595a22830..c8d02ed281ce706383df8fd39cb914da288daf12 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 452f13fb1f126e7af325f9bad7f17a7d40485a4e..a53391a90bbdf070db5e42dd5cc6a381dd02cc1e 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 7d3284e0d4b356350c7d17bf16535e2bda4c48b6..05c34376ece0c21e97e8558f46daaa64b17c7ff0 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 47ec8ac9f63448dd1844c5db6311776cba23a95a..6575a09f05ed851eefccfabca88fdbb04cbb7209 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 9982ddf546798fbb979a59af5cc7efd99a6477cf..e92838716c32a974284f571527ab4a04b324932d 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 5d5dde848c7d83355d24659f3e21ef9c5e55eba3..7e394e9c10b0633262cbaee38e949c747e785116 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 68bdcc77e0ee43efa5b914150f70b0fd47e6cbfd..c76143714ec228016218207ed8be81e186904c67 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 0bd4a7146d491806f94ddab676670bf331d835ba..d31bead5e10e16ea095e99a2b747bcd7e27c7c9c 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 6b7084de633a1305fff57a30d8776d57eeb41b7c..9e465511221baf25438c93e5b2828c3aba915fd0 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 ce9f819646c185af9af7e50013ce66d3c819f52c..4f6eaf3550cc92ec09f807e6b067b37819b7b6e9 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 11f140577faea67a42ae5658936052c59841b30d..7b5249d57e7bd2ce99db1f702aa863eb56a1adf2 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 cb5b083b885c2a8837e85b1cf2eff787a55a205d..f6c8ceec21d28dc6064fee4c129d1d875d724b02 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 9b48bb939df2b32b8cfc2ab770cade94b13dadd0..055eff6f0b34ec7e30fb3da5d8ac9f6f96386d25 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 1447f5b6be5b71db0523889e4349295c8aac40a7..b15fb81e13b3a26efa12c6b113109134f0f65af2 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 8522e609185a8673b0c96c54e5de1ce103fae54b..90033814d33112aa720d016d1f901077ebd095a2 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 37892c0e73b69292d8009b0c7034a3af370388ba..51a2da33bdd8918dba99277e42b820c194cfa584 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 5b238ccefa33658d074c66dde27b851673ea59dc..4c337c4269a26bb625500529bf5d000cf210c6f0 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 c326215512ee705e7c82819d9879d6da9da3814e..859fc9cc4e8aa240f7d280467dfa42017c67449a 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 a0a91c3255451b09dcfc9b0b8408a144760759cd..a4a3fb2e2914150398b4ea3f6da31ee5e58e8b4e 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 26ffef825a5927cec9358ff1299df282e85e64d9..70ac37c1b80e29174dfa5c7447107873a05df60e 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 b586e77fad5a2928a4aac611f5dae8076c7e2b97..e77f66b880ec5c93cf48dc8027888544e4b45d7b 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 6cc93506dfd2403c37952e246725bb2914f28fa6..6797c44ccb6d6200e7be67f0197cf032d803379e 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 e6723f47b21b7aab408ef3a797c46244f3d5e20d..2c7df14ae7864b32f598331db27d1b3cacefaef9 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 42cf70fce5db074854aa2d777e3aff160d440725..43edce0c5d655d44a31d934370be5062bfe34ff4 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 f28bb8c2aa601433d24071d3954542f108dea88f..7501f5806283e70ca3e0dc70203e64796e7ba0be 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 fd85d96714d3adde630d67ab16280487d91bac3a..0b0a31c669c880fd6b8ab7484c7f8196cadc9471 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 9dc9f2bf8005e55a4b9d93759942abcd93a0ffdc..2f5f833d20d9adff0ad89459320c6a4eace77044 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 d2616fc556788e72e4329ba9fef87a34fb483784..4cddb129210821821d84a7d6873d35dde04627f6 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 7d4c5f9039d24e6053b0a98181d94975f78e7458..76a44bd13489c348aa5720e00f7c2ea67eddc180 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 a74f6c2664b39a9b971f7499dfb621773a1dca2e..dac95a197465195446a47e6148ad796c5a969f14 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 79e52a2856a73b68857a05cd65e69f26a55190cd..7c0bf86c94f61c7461686703c207d6a89a583312 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 c17b649212cb423493c743dfc893c9940b8103f9..a2a95cce05733d902c263dac63511a67309d99d5 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 74c6b7b46a7b96a218e332d29218a4488aef0284..de31cefaa365a8e8f846a9510a0fa7a80d3916a6 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 42e0fe70e49adad0b8fab7715f7bce9a1977b1c7..05a049b4c6e2aa4b4e3993ce812c108f2aad68f0 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 104bc0047ff7094bf8d5c83bbd36ca299ff76ea6..ba312ab248e5c8401ba3311a10caa303a3e74758 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 d618226411a105084d34428bb5d6628373c6249c..444f2de81f81692ceda80ebb1a0b30f8244c0a17 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 870034570ac904548f63198fc3de87877d7d1e18..d801fb8fa5b066b33f9748e2b5521a5dbd4cb076 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 79e52a2856a73b68857a05cd65e69f26a55190cd..7c0bf86c94f61c7461686703c207d6a89a583312 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 2f3cb1aaa032019aba8a4278fe2e973e7f7b8c85..d148d2cefd786c30d19b192d9dde8d819711d9bd 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 c2621ae66af286cf3bff0ba0bcae87a72d650160..b4188777437f2b2243d6ab8786c44adad822d8e4 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 63321bab7d8c8835d3ea958f27986e52e3731489..6855b6fb962539e7c9a22c6565fe13f28a81e2e0 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 b0b5a591b3f751dd30ce2faf896beae468b67e51..a553a01b5e68edd687ff0faf6d32368cc7fe0db6 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 97295969b91d18702d43146787d046f1c8808a67..87996b789d704b038a498ee211ddf25645d35931 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 c0c58423193052d9a75326e3f0697d6f5cd1f825..f2d8b4f9fb8043e051419d223ea70996c7624fdd 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 6a723b29febfb456be16ecd0e3e9e566bc2f80d9..97189c353f5d7e763f2c5bf54ea4b4cfa1cb974e 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 aabc1320f332c77d0c503a9549367137186bf071..4165b621912746066682df913e84d6d08842e039 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 9d4be317ff1e0fe2dc11ecb674f12f38966c9440..67473222f51a0c5679cba8328ae7c6c1b0cc2f54 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 e6c7239f615efbc55e1a01c94de72ceb9d30036c..5a72c9d93201f3147b33ce602fc038ed691be574 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 a3b9b841e1d7e649fee41f2ffce53d3b424b1637..399487792213c29ee6524c6d3d9337fe2dd26eef 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 e3127a453beecbfa83a8da8ae5ed7aacf9c4a2f7..a405a2c3f363641b699f9ab830a36b6fd3471c67 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 62cca7ec3a8693a49f4175ac1a7beb1f94e68487..cd0e1adc7e0ceaf3403767bbde21abb06e7379d5 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 2b36031b9274daf780a5a12155d6bb158489ccb2..f475ed24f47a067ceef7a3fc3f01009a3d243b1f 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 69d707974abbcec8faae131f8c6ae834f1d7a769..da4c454a506c4afd4ddf2143e6202897eb52fc34 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 ab9c346e3fe8f2641cbd91fb70f121171f8b74d1..b1bf586e668ed53500b95cf0277c702f021327a5 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
-(
-    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
+Foam::tmp<Foam::surfaceScalarField> Foam::fvc::flux
 (
-    const tmp<surfaceScalarField>& tphi,
-    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf,
-    const word& name
+    const volVectorField& vvf
 )
 {
-    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
+    return scheme<vector>
     (
-        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()+')'
-    );
+        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 tmp<volVectorField>& tvvf
 )
 {
-    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> Flux
-    (
-        fvc::flux(tphi(), vf)
-    );
-    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
-)
-{
-    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 0e3e3a9fe6b9c6b171adc89c1dc21674743b0f46..e62dcec0cb9c89d30a9fdf4781fb806db956ca50 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 0000000000000000000000000000000000000000..ab9c346e3fe8f2641cbd91fb70f121171f8b74d1
--- /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 f7707e0db7cc3ee4165040aea183b1f763f3d85d..c7c38d79e84f758450fd81c7ab4fbdc4cd991394 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 fd9df516e8f85abd1ca7a0ca14816944b2129332..937a359d85bbb4227234a1099471995111de4645 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 cb23171ea69a4bf21d4bd6820a7df412f67e2dff..6a2af46f669d210aa6affe36fbc4130dea879853 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 389e72018a813b0df00be13b7dc512ef238148c6..185284d21ecc2ad267c59f0327028c6989645702 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.H
@@ -157,20 +157,21 @@ namespace fvc
     );
 
 
-    //- Interpolate tmp field onto faces using 'interpolate(\<name\>)'
+    //- Interpolate field onto faces using 'interpolate(\<name\>)'
     template<class Type>
     static tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate
     (
-        const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf
+        const GeometricField<Type, fvPatchField, volMesh>& tvf
     );
 
-    //- Interpolate field onto faces using 'interpolate(\<name\>)'
+    //- Interpolate tmp field onto faces using 'interpolate(\<name\>)'
     template<class Type>
     static tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate
     (
-        const GeometricField<Type, fvPatchField, volMesh>& tvf
+        const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf
     );
 
+
     //- Interpolate boundary field onto faces (simply a type conversion)
     template<class Type>
     static tmp<FieldField<fvsPatchField, Type>> interpolate
@@ -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 ba28c423c9e9188e9fa61f7253762eaea71c711c..1d4d885be499a08291e2af840ddec034fe9597d4 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 529b95f964f0011a1e6e0d94de995c84ca31ede2..ae5357b1851517aee7a1ec255de693836b2ec3fd 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 4ffc27c8d9ec1f30f377db2cae31967599a4d2f2..e1c65f2ad374bcf987818b54751ab9847ce5fb54 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 309cee03ebbc604c538f4d428e7eb05e26ccf40c..12c112d4e08c5b98258dc62f526d3b9df7739cb5 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 309cee03ebbc604c538f4d428e7eb05e26ccf40c..12c112d4e08c5b98258dc62f526d3b9df7739cb5 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;
 }
 
 // ************************************************************************* //