diff --git a/applications/solvers/multiphase/bubbleFoam/pEqn.H b/applications/solvers/multiphase/bubbleFoam/pEqn.H
index 631c9cbba6..af3f6e2822 100644
--- a/applications/solvers/multiphase/bubbleFoam/pEqn.H
+++ b/applications/solvers/multiphase/bubbleFoam/pEqn.H
@@ -14,6 +14,12 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
+    U.boundaryField() =
+        alpha1.boundaryField()*U1.boundaryField()
+      + alpha2.boundaryField()*U2.boundaryField();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
+
     surfaceScalarField phiDrag1
     (
         fvc::interpolate(alpha2/rho1*dragCoef*rAU1)*phi2
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 9e67a47c4f..e0c13b16e0 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -29,6 +29,8 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
index 709d30c398..915e714082 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
@@ -22,6 +22,12 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
+    U.boundaryField() =
+        alpha1.boundaryField()*U1.boundaryField()
+      + alpha2.boundaryField()*U2.boundaryField();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
+
     surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
 
     if (g0.value() > 0.0)
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
index 9c56c4d26f..28de546ce0 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
@@ -5,6 +5,11 @@
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
+    mrfZones.absoluteFlux(phi);
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
+    mrfZones.relativeFlux(phi);
+
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H
index d1bc7bf1ed..a21d5b7bc8 100644
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@@ -4,6 +4,8 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index e3ef0dc9ef..aa3ed1fe85 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@@ -4,6 +4,8 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
index 7bc1eeb12c..5e87fa0319 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
@@ -5,6 +5,11 @@
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
+    mrfZones.absoluteFlux(phi);
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
+    mrfZones.relativeFlux(phi);
+
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
index c1e346ce43..2e26fb538a 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@@ -4,6 +4,8 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H
index 9f27e992a8..99c15656dc 100644
--- a/applications/solvers/multiphase/settlingFoam/pEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/pEqn.H
@@ -5,6 +5,9 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(rho.boundaryField()*U.boundaryField())
+      & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
index ddfca4e3ea..810a753c9d 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
@@ -4,6 +4,8 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 9b543bcd25..d6bde3812b 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -18,6 +18,13 @@
     mrfZones.absoluteFlux(phi1);
     mrfZones.absoluteFlux(phi2.oldTime());
     mrfZones.absoluteFlux(phi2);
+    mrfZones.absoluteFlux(phi);
+
+    U.boundaryField() =
+        alpha1.boundaryField()*U1.boundaryField()
+      + alpha2.boundaryField()*U2.boundaryField();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
 
@@ -56,6 +63,7 @@
     mrfZones.relativeFlux(phi1);
     mrfZones.relativeFlux(phi2.oldTime());
     mrfZones.relativeFlux(phi2);
+    mrfZones.relativeFlux(phi);
 
     surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H.old b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H.old
deleted file mode 100644
index 348cc847d7..0000000000
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H.old
+++ /dev/null
@@ -1,118 +0,0 @@
-{
-    surfaceScalarField alpha1f(fvc::interpolate(alpha1));
-    surfaceScalarField alpha2f(scalar(1) - alpha1f);
-
-    volScalarField rAU1(1.0/U1Eqn.A());
-    volScalarField rAU2(1.0/U2Eqn.A());
-
-    rAU1f = fvc::interpolate(rAU1);
-    surfaceScalarField rAU2f(fvc::interpolate(rAU2));
-
-    volVectorField HbyA1("HbyA1", U1);
-    HbyA1 = rAU1*U1Eqn.H();
-
-    volVectorField HbyA2("HbyA2", U2);
-    HbyA2 = rAU2*U2Eqn.H();
-
-    mrfZones.absoluteFlux(phi1.oldTime());
-    mrfZones.absoluteFlux(phi1);
-
-    mrfZones.absoluteFlux(phi2.oldTime());
-    mrfZones.absoluteFlux(phi2);
-
-    surfaceScalarField phiDrag1
-    (
-        fvc::interpolate(alpha2/rho1*K*rAU1)*phi2 + rAU1f*(g & mesh.Sf())
-    );
-
-    if (g0.value() > 0.0)
-    {
-        phiDrag1 -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
-    }
-
-    if (kineticTheory.on())
-    {
-        phiDrag1 -= rAU1f*fvc::snGrad(kineticTheory.pa()/rho1)*mesh.magSf();
-    }
-
-
-    surfaceScalarField phiDrag2
-    (
-        fvc::interpolate(alpha1/rho2*K*rAU2)*phi1 + rAU2f*(g & mesh.Sf())
-    );
-
-    // Fix for gravity on outlet boundary.
-    forAll(p.boundaryField(), patchi)
-    {
-        if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
-        {
-            phiDrag1.boundaryField()[patchi] = 0.0;
-            phiDrag2.boundaryField()[patchi] = 0.0;
-        }
-    }
-
-    surfaceScalarField phiHbyA1
-    (
-        "phiHbyA1",
-        (fvc::interpolate(HbyA1) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU1, U1, phi1)
-      + phiDrag1
-    );
-    mrfZones.relativeFlux(phiHbyA1);
-
-    surfaceScalarField phiHbyA2
-    (
-        "phiHbyA2",
-        (fvc::interpolate(HbyA2) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU2, U2, phi2)
-      + phiDrag2
-    );
-    mrfZones.relativeFlux(phiHbyA2);
-
-    mrfZones.relativeFlux(phi1.oldTime());
-    mrfZones.relativeFlux(phi1);
-    mrfZones.relativeFlux(phi2.oldTime());
-    mrfZones.relativeFlux(phi2);
-
-    surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
-
-    surfaceScalarField Dp
-    (
-        "Dp",
-        alpha1f*rAU1f/rho1 + alpha2f*rAU2f/rho2
-    );
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
-        );
-
-        pEqn.setReference(pRefCell, pRefValue);
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            surfaceScalarField SfGradp(pEqn.flux()/Dp);
-
-            phi1 = phiHbyA1 - rAU1f*SfGradp/rho1;
-            phi2 = phiHbyA2 - rAU2f*SfGradp/rho2;
-            phi = alpha1f*phi1 + alpha2f*phi2;
-
-            p.relax();
-            SfGradp = pEqn.flux()/Dp;
-
-            U1 = HbyA1 + fvc::reconstruct(phiDrag1 - rAU1f*SfGradp/rho1);
-            U1.correctBoundaryConditions();
-
-            U2 = HbyA2 + fvc::reconstruct(phiDrag2 - rAU2f*SfGradp/rho2);
-            U2.correctBoundaryConditions();
-
-            U = alpha1*U1 + alpha2*U2;
-        }
-    }
-}
-
-#include "continuityErrs.H"
diff --git a/applications/utilities/mesh/manipulation/topoSet/topoSetDict b/applications/utilities/mesh/manipulation/topoSet/topoSetDict
index be90a8fb80..7a2049c686 100644
--- a/applications/utilities/mesh/manipulation/topoSet/topoSetDict
+++ b/applications/utilities/mesh/manipulation/topoSet/topoSetDict
@@ -171,13 +171,13 @@ FoamFile
 //        insidePoint (1 2 3);    // point inside region to select
 //    }
 //
-//    // Cells underneath plane such that volume is reached. Can be used
-//    // in setFields.
+//    // Cells underneath plane such that volume is reached. E.g. for use
+//    // in setFields to set the level given a wanted volume.
 //    source targetVolumeToCell;
 //    sourceInfo
 //    {
 //        volume  2e-05;
-//        normal  (1 1 1);
+//        normal  (0 1 0);          // usually in direction of gravity
 //    }
 //
 //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/checkMeshMoving.H b/applications/utilities/postProcessing/dataConversion/foamToEnsight/checkMeshMoving.H
index 123036f63a..8657a9cc21 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/checkMeshMoving.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/checkMeshMoving.H
@@ -1,21 +1,27 @@
-// check for "points" in all of the result directories
+// check for "points" in any of the result directories
 
-bool meshMoving = true;
-if (Times.size() > 2)
+bool meshMoving = false;
+if (Times.size() > 1)
 {
-    for (label n1=2; n1<Times.size() && meshMoving; ++n1)
+    // We already loaded a mesh (usually from constant). See if any other
+    // points files
+    forAll(Times, timeI)
     {
-        meshMoving = IOobject
-        (
-            "points",
-            Times[n1].name(),
-            polyMesh::meshSubDir,
-            mesh,
-            IOobject::NO_READ
-        ).headerOk();
+        if (Times[timeI].name() != mesh.pointsInstance())
+        {
+            IOobject io
+            (
+                "points",
+                Times[timeI].name(),
+                polyMesh::meshSubDir,
+                mesh,
+                IOobject::NO_READ
+            );
+            if (io.headerOk())
+            {
+                meshMoving = true;
+                break;
+            }
+        }
     }
 }
-else
-{
-    meshMoving = false;
-}
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
index b8fe91fcfe..8a44c5411b 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
@@ -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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1014,6 +1014,7 @@ void Foam::ensightMesh::write
     const fileName& postProcPath,
     const word& prepend,
     const label timeIndex,
+    const bool meshMoving,
     Ostream& ensightCaseFile
 ) const
 {
@@ -1027,7 +1028,7 @@ void Foam::ensightMesh::write
     {
         timeFile += "000.";
     }
-    else if (mesh_.moving())
+    else if (meshMoving)
     {
         timeFile += itoa(timeIndex) + '.';
     }
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.H b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.H
index 6d860dfc91..ffcbf61d66 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -380,6 +380,7 @@ public:
             const fileName& postProcPath,
             const word& prepend,
             const label timeIndex,
+            const bool meshMoving,
             Ostream& ensightCaseFile
         ) const;
 
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
index e9a942f38a..05db6c0d05 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
@@ -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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -228,6 +228,13 @@ int main(int argc, char *argv[])
 
 #   include "checkMeshMoving.H"
 
+    if (meshMoving)
+    {
+        Info<< "Detected a moving mesh (multiple polyMesh/points files)."
+            << " Writing meshes for every timestep." << endl;
+    }
+
+
     wordHashSet allCloudNames;
     if (Pstream::master())
     {
@@ -346,13 +353,14 @@ int main(int argc, char *argv[])
             eMesh.correct();
         }
 
-        if (timeIndex == 0 || (meshState != polyMesh::UNCHANGED))
+        if (timeIndex == 0 || meshMoving)
         {
             eMesh.write
             (
                 ensightDir,
                 prepend,
                 timeIndex,
+                meshMoving,
                 ensightCaseFile
             );
         }
diff --git a/src/OpenFOAM/db/IOobjects/IOdictionary/IOdictionaryIO.C b/src/OpenFOAM/db/IOobjects/IOdictionary/IOdictionaryIO.C
index 881250f3eb..e620b3b942 100644
--- a/src/OpenFOAM/db/IOobjects/IOdictionary/IOdictionaryIO.C
+++ b/src/OpenFOAM/db/IOobjects/IOdictionary/IOdictionaryIO.C
@@ -95,14 +95,16 @@ void Foam::IOdictionary::readFile(const bool masterOnly)
             }
 
             // Note: use ASCII for now - binary IO of dictionaries is
-            // not currently supported
+            // not currently supported or rather the primitiveEntries of
+            // the dictionary think they are in binary form whereas they are
+            // not. Could reset all the ITstreams to ascii?
             IPstream fromAbove
             (
                 Pstream::scheduled,
                 myComm.above(),
                 0,
                 Pstream::msgType(),
-                IOstream::BINARY
+                IOstream::ASCII
             );
             IOdictionary::readData(fromAbove);
         }
@@ -121,7 +123,7 @@ void Foam::IOdictionary::readFile(const bool masterOnly)
                 myComm.below()[belowI],
                 0,
                 Pstream::msgType(),
-                IOstream::BINARY
+                IOstream::ASCII
             );
             IOdictionary::writeData(toBelow);
         }
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
index f4f4f54f65..1cb9326e32 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
@@ -333,7 +333,7 @@ tmp<FieldField<fvsPatchField, Type> > interpolate
     const tmp<FieldField<fvPatchField, Type> >& tfvpff
 )
 {
-    tmp<FieldField<fvPatchField, Type> > tfvspff = interpolate(tfvpff());
+    tmp<FieldField<fvsPatchField, Type> > tfvspff = interpolate(tfvpff());
     tfvpff.clear();
     return tfvspff;
 }
diff --git a/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.C b/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.C
index c64c29e1e6..a80688a184 100644
--- a/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.C
+++ b/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.C
@@ -27,6 +27,7 @@ License
 #include "polyMesh.H"
 #include "globalMeshData.H"
 #include "plane.H"
+#include "cellSet.H"
 
 #include "addToRunTimeSelectionTable.H"
 
@@ -74,6 +75,7 @@ Foam::scalar Foam::targetVolumeToCell::volumeOfSet
 Foam::label Foam::targetVolumeToCell::selectCells
 (
     const scalar normalComp,
+    const PackedBoolList& maskSet,
     PackedBoolList& selected
 ) const
 {
@@ -86,7 +88,7 @@ Foam::label Foam::targetVolumeToCell::selectCells
     {
         const point& cc = mesh_.cellCentres()[cellI];
 
-        if ((cc&n_) < normalComp)
+        if (maskSet[cellI] && ((cc&n_) < normalComp))
         {
             selected[cellI] = true;
             nSelected++;
@@ -103,14 +105,25 @@ void Foam::targetVolumeToCell::combine(topoSet& set, const bool add) const
         // Select no cells
         return;
     }
-    else if (gSum(mesh_.cellVolumes()) < vol_)
+
+
+    PackedBoolList maskSet(mesh_.nCells(), 1);
+    label nTotCells = mesh_.globalData().nTotalCells();
+    if (maskSetName_.size())
     {
-        // Select all cells
-        forAll(mesh_.cellVolumes(), cellI)
+        // Read cellSet
+        Info<< "    Operating on subset defined by cellSet " << maskSetName_
+            << endl;
+
+        maskSet = 0;
+        cellSet subset(mesh_, maskSetName_);
+
+        forAllConstIter(cellSet, subset, iter)
         {
-            addOrDelete(set, cellI, add);
+            maskSet[iter.key()] = 1;
         }
-        return;
+
+        nTotCells = returnReduce(subset.size(), sumOp<label>());
     }
 
 
@@ -144,16 +157,16 @@ void Foam::targetVolumeToCell::combine(topoSet& set, const bool add) const
         }
 
         PackedBoolList maxSelected(mesh_.nCells());
-        maxCells = selectCells(maxComp, maxSelected);
+        maxCells = selectCells(maxComp, maskSet, maxSelected);
         maxVol = volumeOfSet(maxSelected);
 
         // Check that maxPoint indeed selects all cells
-        if (maxCells != mesh_.globalData().nTotalCells())
+        if (maxCells != nTotCells)
         {
             WarningIn("targetVolumeToCell::combine(topoSet&, const bool) const")
                 << "Plane " << plane(points[maxPointI], n_)
                 << " selects " << maxCells
-                << " cells instead of all " << mesh_.globalData().nTotalCells()
+                << " cells instead of all " << nTotCells
                 << " cells. Results might be wrong." << endl;
         }
     }
@@ -178,7 +191,7 @@ void Foam::targetVolumeToCell::combine(topoSet& set, const bool add) const
     {
         scalar mid = 0.5*(low + high);
 
-        nSelected = selectCells(mid, selected);
+        nSelected = selectCells(mid, maskSet, selected);
         selectedVol = volumeOfSet(selected);
 
         //Pout<< "High:" << high << " low:" << low << " mid:" << mid << nl
@@ -191,7 +204,7 @@ void Foam::targetVolumeToCell::combine(topoSet& set, const bool add) const
             low = mid;
 
             PackedBoolList highSelected(mesh_.nCells());
-            label nHigh = selectCells(high, selected);
+            label nHigh = selectCells(high, maskSet, selected);
             if (nSelected == nHigh)
             {
                 break;
@@ -202,7 +215,7 @@ void Foam::targetVolumeToCell::combine(topoSet& set, const bool add) const
             high = mid;
 
             PackedBoolList lowSelected(mesh_.nCells());
-            label nLow = selectCells(low, selected);
+            label nLow = selectCells(low, maskSet, selected);
             if (nSelected == nLow)
             {
                 break;
@@ -210,7 +223,7 @@ void Foam::targetVolumeToCell::combine(topoSet& set, const bool add) const
         }
     }
 
-    nSelected = selectCells(high, selected);
+    nSelected = selectCells(high, maskSet, selected);
     selectedVol = volumeOfSet(selected);
 
     if (selectedVol < vol_)
@@ -219,7 +232,7 @@ void Foam::targetVolumeToCell::combine(topoSet& set, const bool add) const
     }
     else
     {
-        nSelected = selectCells(low, selected);
+        nSelected = selectCells(low, maskSet, selected);
         selectedVol = volumeOfSet(selected);
 
         if (selectedVol < vol_)
@@ -278,7 +291,8 @@ Foam::targetVolumeToCell::targetVolumeToCell
 :
     topoSetSource(mesh),
     vol_(readScalar(dict.lookup("volume"))),
-    n_(dict.lookup("normal"))
+    n_(dict.lookup("normal")),
+    maskSetName_(dict.lookupOrDefault<word>("set", ""))
 {}
 
 
diff --git a/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.H b/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.H
index 398e49c9f9..076f15e8aa 100644
--- a/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.H
+++ b/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.H
@@ -37,7 +37,6 @@ SourceFiles
 #define targetVolumeToCell_H
 
 #include "topoSetSource.H"
-//#include "plane.H"
 #include "PackedBoolList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -63,8 +62,10 @@ class targetVolumeToCell
         const scalar vol_;
 
         //- Normal of plane to sweep
-        vector n_;
-        //plane pl_;
+        const vector n_;
+
+        //- Optional name of cellSet to calculate volume in
+        const word maskSetName_;
 
 
     // Private Member Functions
@@ -74,6 +75,7 @@ class targetVolumeToCell
         label selectCells
         (
             const scalar normalComp,
+            const PackedBoolList&,
             PackedBoolList& selected
         ) const;
 
@@ -92,7 +94,6 @@ public:
         (
             const polyMesh& mesh,
             const scalar vol,
-            //const plane&
             const vector&
         );
 
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
index 8dbf0ec5a8..4aa1596b5a 100644
--- a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
@@ -547,7 +547,7 @@ void Foam::regionSizeDistribution::write()
 
 
         // Sum all regions
-        const scalarField alphaVol = alpha.internalField()*mesh.V();
+        const scalarField alphaVol(alpha.internalField()*mesh.V());
         Map<scalar> allRegionVolume(regionSum(regions, mesh.V()));
         Map<scalar> allRegionAlphaVolume(regionSum(regions, alphaVol));
         Map<label> allRegionNumCells
diff --git a/tutorials/compressible/rhoPimpleFoam/les/pitzDaily/0/B b/tutorials/compressible/rhoPimpleFoam/les/pitzDaily/0/B
deleted file mode 100644
index 39354d72f9..0000000000
--- a/tutorials/compressible/rhoPimpleFoam/les/pitzDaily/0/B
+++ /dev/null
@@ -1,52 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volTensorField;
-    object      B;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 2 -2 0 0 0 0];
-
-internalField   uniform (0 0 0 0 0 0 0 0 0);
-
-boundaryField
-{
-    inlet
-    {
-        type            fixedValue;
-        value           uniform (0 0 0 0 0 0 0 0 0);
-    }
-
-    outlet
-    {
-        type            inletOutlet;
-        inletValue      uniform (0 0 0 0 0 0 0 0 0);
-        value           uniform (0 0 0 0 0 0 0 0 0);
-    }
-
-    upperWall
-    {
-        type            zeroGradient;
-    }
-
-    lowerWall
-    {
-        type            zeroGradient;
-    }
-
-    frontAndBack
-    {
-        type            empty;
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/channel395/0.org/B b/tutorials/incompressible/pimpleFoam/channel395/0.org/B
deleted file mode 100644
index 78d2540b2e..0000000000
--- a/tutorials/incompressible/pimpleFoam/channel395/0.org/B
+++ /dev/null
@@ -1,67 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volTensorField;
-    location    "1";
-    object      B;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [ 0 2 -2 0 0 0 0 ];
-
-internalField   uniform ( 0 0 0 0 0 0 0 0 0 );
-
-boundaryField
-{
-    bottomWall
-    {
-        type            zeroGradient;
-    }
-    topWall
-    {
-        type            zeroGradient;
-    }
-    sides1_half0
-    {
-        type            cyclic;
-    }
-    sides2_half0
-    {
-        type            cyclic;
-    }
-    inout1_half0
-    {
-        type            cyclic;
-    }
-    inout2_half0
-    {
-        type            cyclic;
-    }
-    sides2_half1
-    {
-        type            cyclic;
-    }
-    sides1_half1
-    {
-        type            cyclic;
-    }
-    inout1_half1
-    {
-        type            cyclic;
-    }
-    inout2_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/pisoFoam/les/pitzDaily/0/B b/tutorials/incompressible/pisoFoam/les/pitzDaily/0/B
deleted file mode 100644
index 39354d72f9..0000000000
--- a/tutorials/incompressible/pisoFoam/les/pitzDaily/0/B
+++ /dev/null
@@ -1,52 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volTensorField;
-    object      B;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 2 -2 0 0 0 0];
-
-internalField   uniform (0 0 0 0 0 0 0 0 0);
-
-boundaryField
-{
-    inlet
-    {
-        type            fixedValue;
-        value           uniform (0 0 0 0 0 0 0 0 0);
-    }
-
-    outlet
-    {
-        type            inletOutlet;
-        inletValue      uniform (0 0 0 0 0 0 0 0 0);
-        value           uniform (0 0 0 0 0 0 0 0 0);
-    }
-
-    upperWall
-    {
-        type            zeroGradient;
-    }
-
-    lowerWall
-    {
-        type            zeroGradient;
-    }
-
-    frontAndBack
-    {
-        type            empty;
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/pisoFoam/les/pitzDailyMapped/0/B b/tutorials/incompressible/pisoFoam/les/pitzDailyMapped/0/B
deleted file mode 100644
index 3cca215a87..0000000000
--- a/tutorials/incompressible/pisoFoam/les/pitzDailyMapped/0/B
+++ /dev/null
@@ -1,55 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volTensorField;
-    object      B;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 2 -2 0 0 0 0];
-
-internalField   uniform (0 0 0 0 0 0 0 0 0);
-
-boundaryField
-{
-    inlet
-    {
-        type                mapped;
-        value               uniform (0 0 0 0 0 0 0 0 0);
-        interpolationScheme cell;
-        setAverage          false;
-        average             (0 0 0 0 0 0 0 0 0);
-    }
-
-    outlet
-    {
-        type                inletOutlet;
-        inletValue          uniform (0 0 0 0 0 0 0 0 0);
-        value               uniform (0 0 0 0 0 0 0 0 0);
-    }
-
-    upperWall
-    {
-        type                zeroGradient;
-    }
-
-    lowerWall
-    {
-        type                zeroGradient;
-    }
-
-    frontAndBack
-    {
-        type                empty;
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/multiphase/interFoam/les/nozzleFlow2D/0/B b/tutorials/multiphase/interFoam/les/nozzleFlow2D/0/B
deleted file mode 100644
index d84b782054..0000000000
--- a/tutorials/multiphase/interFoam/les/nozzleFlow2D/0/B
+++ /dev/null
@@ -1,57 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volTensorField;
-    object      B;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 2 -2 0 0 0 0];
-
-internalField   uniform (0 0 0 0 0 0 0 0 0);
-
-boundaryField
-{
-    axis
-    {
-        type            empty;
-    }
-
-    inlet
-    {
-        type            fixedValue;
-        value           uniform (0 0 0 0 0 0 0 0 0);
-    }
-
-    wall
-    {
-        type            zeroGradient;
-    }
-
-    atmosphere
-    {
-        type            inletOutlet;
-        inletValue      uniform (0 0 0 0 0 0 0 0 0);
-        value           uniform (0 0 0 0 0 0 0 0 0);
-    }
-
-    front
-    {
-        type            wedge;
-    }
-
-    back
-    {
-        type            wedge;
-    }
-}
-
-// ************************************************************************* //