diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index cda60d98b2..c778ef40a8 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -41,7 +41,6 @@
         p_rghDDtEqn =
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-          + fvc::div(phiHbyA)
         ==
             fvOptions(psi, p_rgh, rho.name())
         );
@@ -52,6 +51,7 @@
         fvScalarMatrix p_rghEqn
         (
             p_rghDDtEqn()
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rhorAUf, p_rgh)
         );
 
diff --git a/doc/Doxygen/css/doxyMod.css b/doc/Doxygen/css/doxyMod.css
index 46ac2dfc5e..fbb968f594 100644
--- a/doc/Doxygen/css/doxyMod.css
+++ b/doc/Doxygen/css/doxyMod.css
@@ -118,6 +118,8 @@ tr.memlist
 .OFPlainTable tr td {
     height: 20px;
     padding-left: 5px;
+}
+
 div.line,
 span.comment,
 span.keyword,
diff --git a/etc/config.sh/FFTW b/etc/config.sh/FFTW
index 7c0a48861d..f3cb4ea70e 100644
--- a/etc/config.sh/FFTW
+++ b/etc/config.sh/FFTW
@@ -64,7 +64,7 @@ then
     # it is either located within ThirdParty, or a central installation
     # outside of ThirdParty and must be added to the lib-path.
 
-    ending="${FFTW_ARCH_PATH_PATH##*-}"
+    ending="${FFTW_ARCH_PATH##*-}"
     if [ "$ending" != none -a "$ending" != system ]
     then
         _foamAddLib $FFTW_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
index a6d05b7f73..5f22c04431 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
+++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
@@ -129,23 +129,29 @@ Foam::plane::plane(const vector& normalVector)
 }
 
 
-Foam::plane::plane(const point& basePoint, const vector& normalVector)
+Foam::plane::plane
+(
+    const point& basePoint,
+    const vector& normalVector,
+    const bool normalise
+)
 :
     unitVector_(normalVector),
     basePoint_(basePoint)
 {
-    scalar magUnitVector(mag(unitVector_));
+    scalar magSqrUnitVector(magSqr(unitVector_));
 
-    if (magUnitVector > VSMALL)
-    {
-        unitVector_ /= magUnitVector;
-    }
-    else
+    if (magSqrUnitVector < VSMALL)
     {
         FatalErrorInFunction
             << "plane normal has zero length. basePoint:" << basePoint_
             << abort(FatalError);
     }
+
+    if (normalise)
+    {
+        unitVector_ /= Foam::sqrt(magSqrUnitVector);
+    }
 }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
index 44564fa92a..6645fb1b95 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
@@ -131,7 +131,12 @@ public:
         plane(const vector& normalVector);
 
         //- Construct from normal vector and point in plane
-        plane(const point& basePoint, const vector& normalVector);
+        plane
+        (
+            const point& basePoint,
+            const vector& normalVector,
+            const bool normalise = true
+        );
 
         //- Construct from three points in plane
         plane(const point& point1, const point& point2, const point& point3);
diff --git a/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C b/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C
index d43f4c2578..77c2d8f310 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C
+++ b/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C
@@ -47,26 +47,31 @@ void Foam::LESModels::maxDeltaxyz::calcDelta()
     label nD = mesh.nGeometricD();
 
     const cellList& cells = mesh.cells();
+    const vectorField& cellC = mesh.cellCentres();
+    const vectorField& faceC = mesh.faceCentres();
+    const vectorField faceN(mesh.faceAreas()/mag(mesh.faceAreas()));
     scalarField hmax(cells.size());
 
-    forAll(cells,cellI)
+    forAll(cells, celli)
     {
         scalar deltaMaxTmp = 0.0;
-        const labelList& cFaces = mesh.cells()[cellI];
-        const point& centrevector = mesh.cellCentres()[cellI];
+        const labelList& cFaces = cells[celli];
+        const point& cc = cellC[celli];
 
-        forAll(cFaces, cFaceI)
+        forAll(cFaces, cFacei)
         {
-            label faceI = cFaces[cFaceI];
-            const point& facevector = mesh.faceCentres()[faceI];
-            scalar tmp = mag(facevector - centrevector);
+            label facei = cFaces[cFacei];
+            const point& fc = faceC[facei];
+            const vector& n = faceN[facei];
+
+            scalar tmp = magSqr(n*(n & (fc - cc)));
             if (tmp > deltaMaxTmp)
             {
                 deltaMaxTmp = tmp;
             }
         }
 
-        hmax[cellI] = deltaCoeff_*deltaMaxTmp;
+        hmax[celli] = deltaCoeff_*Foam::sqrt(deltaMaxTmp);
     }
 
     if (nD == 3)
@@ -76,7 +81,7 @@ void Foam::LESModels::maxDeltaxyz::calcDelta()
     else if (nD == 2)
     {
         WarningInFunction
-            << "Case is 2D, LES is not strictly applicable\n"
+            << "Case is 2D, LES is not strictly applicable" << nl
             << endl;
 
         delta_.internalField() = hmax;
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C
index 5371eda475..ad6ebca09f 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C
@@ -332,6 +332,20 @@ void Foam::activePressureForceBaffleVelocityFvPatchVectorField::updateCoeffs()
             areaFraction = 1 - openFraction_;
         }
 
+        if (patch().boundaryMesh().mesh().moving())
+        {
+            // All areas have been recalculated already so are consistent
+            // with the new points. Note we already only do this routine
+            // once per timestep. The alternative is to use the upToDate
+            // mechanism for regIOobject + polyMesh::upToDatePoints
+            initWallSf_ = patch().Sf();
+            initCyclicSf_ = patch().boundaryMesh()[cyclicPatchLabel_].Sf();
+            nbrCyclicSf_ =  refCast<const cyclicFvPatch>
+            (
+                patch().boundaryMesh()[cyclicPatchLabel_]
+            ).neighbFvPatch().Sf();
+        }
+
         // Update this wall patch
         vectorField::subField Sfw = patch().patch().faceAreas();
         vectorField newSfw((1 - areaFraction)*initWallSf_);
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
index 3f855d52fa..24714aaf46 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
@@ -48,7 +48,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
 :
     uniformJumpFvPatchField<Type>(p, iF),
     phiName_("phi"),
-    rhoName_("rho")
+    rhoName_("rho"),
+    uniformJump_(false)
 {}
 
 
@@ -62,7 +63,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
 :
     uniformJumpFvPatchField<Type>(p, iF, dict),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    uniformJump_(dict.lookupOrDefault<bool>("uniformJump", "false"))
 {}
 
 
@@ -77,7 +79,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
 :
     uniformJumpFvPatchField<Type>(ptf, p, iF, mapper),
     phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_)
+    rhoName_(ptf.rhoName_),
+    uniformJump_(ptf.uniformJump_)
 {}
 
 
@@ -89,7 +92,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
 :
     uniformJumpFvPatchField<Type>(ptf),
     phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_)
+    rhoName_(ptf.rhoName_),
+    uniformJump_(ptf.uniformJump_)
 {}
 
 
@@ -102,7 +106,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
 :
     uniformJumpFvPatchField<Type>(ptf, iF),
     phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_)
+    rhoName_(ptf.rhoName_),
+    uniformJump_(ptf.uniformJump_)
 {}
 
 
@@ -129,6 +134,10 @@ void Foam::fanFvPatchField<Type>::write(Ostream& os) const
     uniformJumpFvPatchField<Type>::write(os);
     this->template writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     this->template writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    this->template writeEntryIfDifferent<bool>
+    (
+        os, "uniformJump", false, uniformJump_
+    );
 }
 
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
index b1080e264f..ccd25cc8c4 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
@@ -42,6 +42,8 @@ Description
         jumpTable    | jump data, e.g. \c csvFile | yes      |
         phi          | flux field name         | no          | phi
         rho          | density field name      | no          | none
+        uniformJump  | applies a uniform pressure based on the averaged
+                       velocity | no | false
     \endtable
 
     Example of the boundary condition specification:
@@ -53,10 +55,11 @@ Description
         jumpTable       csvFile;
         csvFileCoeffs
         {
-            hasHeaderLine   1;
+            nHeaderLine     1;
             refColumn       0;
             componentColumns 1(1);
             separator       ",";
+            mergeSeparators no;
             fileName        "$FOAM_CASE/constant/pressureVsU";
         }
         value           uniform 0;
@@ -67,7 +70,7 @@ Description
     the jump condition.
 
 Note
-     The underlying \c patchType should be set to \c cyclic
+    The underlying \c patchType should be set to \c cyclic
 
 SeeAlso
     Foam::Function1Types
@@ -109,6 +112,9 @@ class fanFvPatchField
         //  if neccessary
         word rhoName_;
 
+        //- Uniform pressure drop
+        bool uniformJump_;
+
 
     // Private Member Functions
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
index a4f2809df8..fd0b61bdaa 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
@@ -44,6 +44,11 @@ void Foam::fanFvPatchField<Foam::scalar>::calcFanJump()
             patch().patchField<surfaceScalarField, scalar>(phi);
 
         scalarField Un(max(phip/patch().magSf(), scalar(0)));
+        if (uniformJump_)
+        {
+            scalar area = gSum(patch().magSf());
+            Un = gSum(Un*patch().magSf())/area;
+        }
 
         if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
         {
@@ -67,7 +72,8 @@ Foam::fanFvPatchField<Foam::scalar>::fanFvPatchField
 :
     uniformJumpFvPatchField<scalar>(p, iF),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    uniformJump_(dict.lookupOrDefault<bool>("uniformJump", false))
 {
     if (this->cyclicPatch().owner())
     {
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
index b0df476fe0..59646f3ad8 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
@@ -848,7 +848,7 @@ turbulentDFSEMInletFvPatchVectorField
     eddy::debug = debug;
 
     // Set UMean as patch area average value
-    UMean_ = gSum(U_*patch().magSf())/gSum(patch().magSf());
+    UMean_ = gSum(U_*patch().magSf())/(gSum(patch().magSf()) + ROOTVSMALL);
 }
 
 
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C
index 5da9b1486c..904d317eda 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C
@@ -40,6 +40,13 @@ namespace fv
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+template<class Type>
+const convectionScheme<Type>& boundedConvectionScheme<Type>::scheme() const
+{
+    return scheme_();
+}
+
+
 template<class Type>
 tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
 boundedConvectionScheme<Type>::interpolate
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H
index 630e17090a..0573ef4fc2 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H
@@ -106,6 +106,8 @@ public:
 
     // Member Functions
 
+        const convectionScheme<Type>& scheme() const;
+
         tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate
         (
             const surfaceScalarField&,
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
index 209c2f7e37..cc54720dd3 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
@@ -39,12 +39,6 @@ namespace Foam
     defineTypeNameAndDebug(tetOverlapVolume, 0);
 }
 
-// When to consider a tet to be zero volume. We want to avoid doing clipping
-// against negative volume tets. Tet volume can be calculated incorrectly
-// due to truncation errors. The value below works for single and double
-// precision but could probably be tighter for double precision.
-Foam::scalar Foam::tetOverlapVolume::minTetVolume_ = SMALL*SMALL;
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
index b027599218..a77e608c6e 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
@@ -54,12 +54,6 @@ class polyMesh;
 
 class tetOverlapVolume
 {
-    // Private data
-
-        //- Minimum tet volume to skip test
-        static scalar minTetVolume_;
-
-
     // Private classes
 
         //- tetPoints handling : sum resulting volumes
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C b/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
index 3ff950311f..34eded1ff3 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,51 +45,90 @@ void Foam::tetOverlapVolume::tetTetOverlap
     tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
     tetPointRef::dummyOp outside;
 
-    if (tetA.tet().mag() < minTetVolume_ || tetB.tet().mag() < minTetVolume_)
+    // Precompute the tet face areas and exit early if any face area is
+    // too small
+    static FixedList<vector, 4> tetAFaceAreas;
+    static FixedList<scalar, 4> tetAMag2FaceAreas;
+    tetPointRef tetATet = tetA.tet();
+    for (label facei = 0; facei < 4; ++facei)
     {
-        return;
+        tetAFaceAreas[facei] = -tetATet.tri(facei).normal();
+        tetAMag2FaceAreas[facei] = magSqr(tetAFaceAreas[facei]);
+        if (tetAMag2FaceAreas[facei] < ROOTVSMALL)
+        {
+            return;
+        }
     }
 
-    // face0
-    plane pl0(tetB[1], tetB[3], tetB[2]);
-    tetA.tet().sliceWithPlane(pl0, cutInside, outside);
-    if (nCutInside == 0)
+    static FixedList<vector, 4> tetBFaceAreas;
+    static FixedList<scalar, 4> tetBMag2FaceAreas;
+    tetPointRef tetBTet = tetB.tet();
+    for (label facei = 0; facei < 4; ++facei)
     {
-        return;
+        tetBFaceAreas[facei] = -tetBTet.tri(facei).normal();
+        tetBMag2FaceAreas[facei] = magSqr(tetBFaceAreas[facei]);
+        if (tetBMag2FaceAreas[facei] < ROOTVSMALL)
+        {
+            return;
+        }
     }
 
-    // face1
-    plane pl1(tetB[0], tetB[2], tetB[3]);
-    nInside = 0;
-    for (label i = 0; i < nCutInside; i++)
+
+    // Face 0
     {
-        const tetPointRef t = cutInsideTets[i].tet();
-        t.sliceWithPlane(pl1, inside, outside);
-    }
-    if (nInside == 0)
-    {
-        return;
+        vector n = tetBFaceAreas[0]/Foam::sqrt(tetBMag2FaceAreas[0]);
+        plane pl0(tetBTet.tri(0).a(), n, false);
+
+        tetA.tet().sliceWithPlane(pl0, cutInside, outside);
+        if (nCutInside == 0)
+        {
+            return;
+        }
     }
 
-    // face2
-    plane pl2(tetB[0], tetB[3], tetB[1]);
-    nCutInside = 0;
-    for (label i = 0; i < nInside; i++)
+    // Face 1
     {
-        const tetPointRef t = insideTets[i].tet();
-        t.sliceWithPlane(pl2, cutInside, outside);
-    }
-    if (nCutInside == 0)
-    {
-        return;
+        vector n = tetBFaceAreas[1]/Foam::sqrt(tetBMag2FaceAreas[1]);
+        plane pl1(tetBTet.tri(1).a(), n, false);
+
+        nInside = 0;
+        for (label i = 0; i < nCutInside; i++)
+        {
+            const tetPointRef t = cutInsideTets[i].tet();
+            t.sliceWithPlane(pl1, inside, outside);
+        }
+        if (nInside == 0)
+        {
+            return;
+        }
     }
 
-    // face3
-    plane pl3(tetB[0], tetB[1], tetB[2]);
-    for (label i = 0; i < nCutInside; i++)
+    // Face 2
     {
-        const tetPointRef t = cutInsideTets[i].tet();
-        t.sliceWithPlane(pl3, insideOp, outside);
+        vector n = tetBFaceAreas[2]/Foam::sqrt(tetBMag2FaceAreas[2]);
+        plane pl2(tetBTet.tri(2).a(), n, false);
+
+        nCutInside = 0;
+        for (label i = 0; i < nInside; i++)
+        {
+            const tetPointRef t = insideTets[i].tet();
+            t.sliceWithPlane(pl2, cutInside, outside);
+        }
+        if (nCutInside == 0)
+        {
+            return;
+        }
+    }
+
+    // Face 3
+    {
+        vector n = tetBFaceAreas[3]/Foam::sqrt(tetBMag2FaceAreas[3]);
+        plane pl3(tetBTet.tri(3).a(), n, false);
+        for (label i = 0; i < nCutInside; i++)
+        {
+            const tetPointRef t = cutInsideTets[i].tet();
+            t.sliceWithPlane(pl3, insideOp, outside);
+        }
     }
 }
 
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
index e74bfb2f69..1d8855d665 100644
--- a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -335,7 +335,8 @@ Foam::regionSizeDistribution::regionSizeDistribution
     active_(true),
     alphaName_(dict.lookup("field")),
     patchNames_(dict.lookup("patches")),
-    log_(true)
+    log_(true),
+    isoPlanes_(dict.lookupOrDefault<bool>("isoPlanes", false))
 {
     // Check if the available mesh is an fvMesh, otherwise deactivate
     if (isA<fvMesh>(obr_))
@@ -388,6 +389,17 @@ void Foam::regionSizeDistribution::read(const dictionary& dict)
                 << "Transforming all vectorFields with coordinate system "
                 << coordSysPtr_().name() << endl;
         }
+
+        if (isoPlanes_)
+        {
+             dict.lookup("origin") >> origin_;
+             dict.lookup("direction") >> direction_;
+             dict.lookup("maxDiameter") >> maxDiameter_;
+             dict.lookup("nDownstreamBins") >> nDownstreamBins_;
+             dict.lookup("maxDownstream") >> maxDownstream_;
+
+             direction_ /= mag(direction_);
+        }
     }
 }
 
@@ -691,6 +703,8 @@ void Foam::regionSizeDistribution::write()
         // (patchRegions) and background regions from maps.
         // Note that we have to use mesh volume (allRegionVolume) and not
         // allRegionAlphaVolume since background might not have alpha in it.
+        // Deleting regions where the volume-alpha-weighted is lower than
+        // threshold
         forAllIter(Map<scalar>, allRegionVolume, vIter)
         {
             label regionI = vIter.key();
@@ -698,6 +712,7 @@ void Foam::regionSizeDistribution::write()
             (
                 patchRegions.found(regionI)
              || vIter() >= maxDropletVol
+             || (allRegionAlphaVolume[regionI]/vIter() < threshold_)
             )
             {
                 allRegionVolume.erase(vIter);
@@ -735,6 +750,108 @@ void Foam::regionSizeDistribution::write()
                 )
             );
 
+            vectorField centroids(sortedVols.size(), vector::zero);
+
+            // Check if downstream bins are calculed
+            if (isoPlanes_)
+            {
+                vectorField alphaDistance
+                (
+                    (alpha.internalField()*mesh.V())
+                  * (mesh.C().internalField() - origin_)()
+                );
+
+                Map<vector> allRegionAlphaDistance
+                (
+                    regionSum
+                    (
+                        regions,
+                        alphaDistance
+                    )
+                );
+
+                // 2. centroid
+                vectorField sortedMoment
+                (
+                    extractData
+                    (
+                        sortedRegions,
+                        allRegionAlphaDistance
+                    )
+                );
+
+                centroids = sortedMoment/sortedVols + origin_;
+
+                // Bin according to centroid
+                scalarField distToPlane((centroids - origin_) & direction_);
+
+                vectorField radialDistToOrigin
+                (
+                    (centroids - origin_) - (distToPlane*direction_)
+                );
+
+                const scalar deltaX = maxDownstream_/nDownstreamBins_;
+                labelList downstreamIndices(distToPlane.size(), -1);
+                forAll(distToPlane, i)
+                {
+                    if
+                    (
+                        (mag(radialDistToOrigin[i]) < maxDiameter_)
+                     && (distToPlane[i] < maxDownstream_)
+                    )
+                    {
+                        downstreamIndices[i] = distToPlane[i]/deltaX;
+                    }
+                }
+
+                scalarField binDownCount(nDownstreamBins_, 0.0);
+                forAll(distToPlane, i)
+                {
+                    if (downstreamIndices[i] != -1)
+                    {
+                        binDownCount[downstreamIndices[i]] += 1.0;
+                    }
+                }
+
+                // Write
+                if (Pstream::master())
+                {
+                    // Construct mids of bins for plotting
+                    pointField xBin(nDownstreamBins_);
+
+                    scalar x = 0.5*deltaX;
+                    forAll(xBin, i)
+                    {
+                        xBin[i] = point(x, 0, 0);
+                        x += deltaX;
+                    }
+
+                    const coordSet coords("distance", "x", xBin, mag(xBin));
+                    writeGraph(coords, "isoPlanes", binDownCount);
+                }
+
+                // Write to screen
+                if (log_)
+                {
+                    Info<< "    Iso-planes Bins:" << endl;
+                    Info<< "    " << token::TAB << "Bin"
+                        << token::TAB << "Min distance"
+                        << token::TAB << "Count:"
+                        << endl;
+
+                    scalar delta = 0.0;
+                    forAll(binDownCount, binI)
+                    {
+                        Info<< "    " << token::TAB << binI
+                            << token::TAB << delta
+                            << token::TAB << binDownCount[binI] << endl;
+                        delta += deltaX;
+                    }
+                    Info<< endl;
+
+                }
+            }
+
             // Calculate the diameters
             scalarField sortedDiameters(sortedVols.size());
             forAll(sortedDiameters, i)
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
index 29c0bc5ac1..657077366b 100644
--- a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,6 +35,8 @@ Description
     based on where the field is below the threshold value.  These
     regions ("droplets") can now be analysed.
 
+
+
     Regions:
     - print the regions connected to a user-defined set of patches.
       (in spray calculation these form the liquid core)
@@ -53,6 +55,9 @@ Description
     - write graph of sum, average and deviation of droplet volume per bin
     - write graph of sum, average and deviation of user-defined fields.  For
       volVectorFields these are those of the 3 components and the magnitude.
+    - (optional) write graph of histogram of centroids on iso planes
+      downstream of the injector determined by origin, direction and maxDiameter
+      up to maxDownstream
 
     Example of function object specification:
     \verbatim
@@ -76,6 +81,24 @@ Description
             e3              (0 1 1);
             e1              (1 0 0);
         }
+
+        // Optional downstream iso-plane bins.
+        isoPlanes       true;
+
+            // Plane normal and point definition
+            direction       (1 0 1);
+            origin          (1e-4 0 5e-4);
+
+            // Maximum diameter of the cylinder formed by the origin point
+            // and direction
+            maxDiameter     3e-4;
+
+            // Maximum downstream distance
+            maxDownstream   6e-4;
+
+            // Number of iso-plane bins
+            nDownstreamBins 20;
+
     }
     \endverbatim
 
@@ -177,6 +200,27 @@ class regionSizeDistribution
         //- Optional coordinate system
         autoPtr<coordinateSystem> coordSysPtr_;
 
+        // Optional extra definition of bins on planes downstream to the origin
+        // point and maximum diameter
+
+            //- Switch to enable iso-planes sampling
+            bool isoPlanes_;
+
+            //- Optional origin point
+            vector origin_;
+
+            //- Optional plane normal direction
+            vector direction_;
+
+            //- Optional maximum diameter on plane
+            scalar maxDiameter_;
+
+            //- Optional number of bins for
+            scalar nDownstreamBins_;
+
+            //- Optional maximum downstream coordinate from origin
+            scalar maxDownstream_;
+
 
     // Private Member Functions
 
diff --git a/src/postProcessing/functionObjects/graphics/runTimePostProcessing/fieldVisualisationBase.C b/src/postProcessing/functionObjects/graphics/runTimePostProcessing/fieldVisualisationBase.C
index 40f60a7c86..86fd9b0e73 100644
--- a/src/postProcessing/functionObjects/graphics/runTimePostProcessing/fieldVisualisationBase.C
+++ b/src/postProcessing/functionObjects/graphics/runTimePostProcessing/fieldVisualisationBase.C
@@ -43,6 +43,7 @@ License
 #include "vtkSphereSource.h"
 #include "vtkTextActor.h"
 #include "vtkTextProperty.h"
+#include "vtkCellDataToPointData.h"
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
@@ -151,6 +152,9 @@ void Foam::fieldVisualisationBase::addScalarBar
     const vector textColour = colours_["text"]->value(position);
 
     // Work-around to supply our own scalarbar title
+    // - Default scalar bar title text is scales by the scalar bar box
+    //   dimensions so if the title is a long string, the text is shrunk to fit
+    //   Instead, suppress title and set the title using a vtkTextActor
     vtkSmartPointer<vtkTextActor> titleActor =
         vtkSmartPointer<vtkTextActor>::New();
     sbar->SetTitle(" ");
@@ -170,19 +174,18 @@ void Foam::fieldVisualisationBase::addScalarBar
     titleActor->GetPositionCoordinate()->
         SetCoordinateSystemToNormalizedViewport();
 
-/*
-    sbar->SetTitle(scalarBar_.title_.c_str());
-    sbar->GetTitleTextProperty()->SetColor
-    (
-        textColour[0],
-        textColour[1],
-        textColour[2]
-    );
-    sbar->GetTitleTextProperty()->SetFontSize(scalarBar_.fontSize_);
-    sbar->GetTitleTextProperty()->ShadowOff();
-    sbar->GetTitleTextProperty()->BoldOn();
-    sbar->GetTitleTextProperty()->ItalicOff();
-*/
+    // How to use the standard scalar bar text
+    // sbar->SetTitle(scalarBar_.title_.c_str());
+    // sbar->GetTitleTextProperty()->SetColor
+    // (
+    //     textColour[0],
+    //     textColour[1],
+    //     textColour[2]
+    // );
+    // sbar->GetTitleTextProperty()->SetFontSize(scalarBar_.fontSize_);
+    // sbar->GetTitleTextProperty()->ShadowOff();
+    // sbar->GetTitleTextProperty()->BoldOn();
+    // sbar->GetTitleTextProperty()->ItalicOff();
 
     sbar->GetLabelTextProperty()->SetColor
     (
@@ -217,8 +220,8 @@ void Foam::fieldVisualisationBase::addScalarBar
         sbar->SetWidth(0.75);
         sbar->SetHeight(0.07);
         sbar->SetBarRatio(0.5);
-//        sbar->SetHeight(0.1);
-//        sbar->SetTitleRatio(0.01);
+        // sbar->SetHeight(0.1);
+        // sbar->SetTitleRatio(0.01);
         sbar->SetTextPositionToPrecedeScalarBar();
     }
 
@@ -228,10 +231,10 @@ void Foam::fieldVisualisationBase::addScalarBar
         scalarBar_.position_.second() + sbar->GetHeight()
     );
 
-//    sbar->DrawFrameOn();
-//    sbar->DrawBackgroundOn();
-//    sbar->UseOpacityOff();
-//    sbar->VisibilityOff();
+    // sbar->DrawFrameOn();
+    // sbar->DrawBackgroundOn();
+    // sbar->UseOpacityOff();
+    // sbar->VisibilityOff();
     sbar->VisibilityOn();
 
     renderer->AddActor(sbar);
@@ -266,19 +269,21 @@ void Foam::fieldVisualisationBase::setField
             lut->SetVectorMode(vtkScalarsToColors::MAGNITUDE);
 
             // Configure the mapper
-            mapper->SelectColorArray(colourFieldName.c_str());
-            mapper->SetScalarRange(range_.first(), range_.second());
-
-            // Set to use either cell or point data
             const char* fieldName = colourFieldName.c_str();
-            if (pData->GetCellData()->HasArray(fieldName) == 1)
-            {
-                mapper->SetScalarModeToUseCellFieldData();
-            }
-            else if (pData->GetPointData()->HasArray(fieldName) == 1)
+            mapper->SelectColorArray(fieldName);
+
+            // Set to use either point or cell data
+            // Note: if both point and cell data exists, preferentially
+            //       choosing point data.  This is often the case when using
+            //       glyphs
+            if (pData->GetPointData()->HasArray(fieldName) == 1)
             {
                 mapper->SetScalarModeToUsePointFieldData();
             }
+            else if (pData->GetCellData()->HasArray(fieldName) == 1)
+            {
+                mapper->SetScalarModeToUseCellFieldData();
+            }
             else
             {
                 WarningInFunction
@@ -286,7 +291,7 @@ void Foam::fieldVisualisationBase::setField
                     << "- assuming point data";
                 mapper->SetScalarModeToUsePointFieldData();
             }
-
+            mapper->SetScalarRange(range_.first(), range_.second());
             mapper->SetColorModeToMapScalars();
             mapper->SetLookupTable(lut);
             mapper->ScalarVisibilityOn();
@@ -322,9 +327,37 @@ void Foam::fieldVisualisationBase::addGlyphs
     glyph->ScalingOn();
     bool ok = true;
 
-    label nComponents =
-        data->GetPointData()->GetArray(scaleFieldName.c_str())
-            ->GetNumberOfComponents();
+    // Determine whether we have scalar or vector data
+    label nComponents = -1;
+    const char* scaleFieldNameChar = scaleFieldName.c_str();
+    if (data->GetPointData()->HasArray(scaleFieldNameChar) == 1)
+    {
+        nComponents =
+            data->GetPointData()->GetArray(scaleFieldNameChar)
+                ->GetNumberOfComponents();
+    }
+    else if (data->GetCellData()->HasArray(scaleFieldNameChar) == 1)
+    {
+        // Need to convert cell data to point data
+        vtkSmartPointer<vtkCellDataToPointData> cellToPoint =
+            vtkSmartPointer<vtkCellDataToPointData>::New();
+        cellToPoint->SetInputData(data);
+        cellToPoint->Update();
+        vtkDataSet* pds = cellToPoint->GetOutput();
+        vtkDataArray* pData = pds->GetPointData()->GetArray(scaleFieldNameChar);
+
+        // Store in main vtkPolyData
+        data->GetPointData()->AddArray(pData);
+
+        nComponents = pData->GetNumberOfComponents();
+    }
+    else
+    {
+        WarningInFunction
+            << "Glyphs can only be added to scalar or vector data. "
+            << "Unable to process field " << scaleFieldName << endl;
+        return;
+    }
 
     if (nComponents == 1)
     {
@@ -332,9 +365,10 @@ void Foam::fieldVisualisationBase::addGlyphs
             vtkSmartPointer<vtkSphereSource>::New();
         sphere->SetCenter(0, 0, 0);
         sphere->SetRadius(0.5);
-// Setting higher resolution slows the rendering significantly
-//        sphere->SetPhiResolution(20);
-//        sphere->SetThetaResolution(20);
+
+        // Setting higher resolution slows the rendering significantly
+        // sphere->SetPhiResolution(20);
+        // sphere->SetThetaResolution(20);
 
         glyph->SetSourceConnection(sphere->GetOutputPort());
 
@@ -342,18 +376,18 @@ void Foam::fieldVisualisationBase::addGlyphs
         {
             double range[2];
 
-// Can use values to find range
-//            vtkDataArray* values =
-//                data->GetPointData()->GetScalars(scaleFieldName.c_str());
-//            values->GetRange(range);
+            // Can use values to find range
+            // vtkDataArray* values =
+            //     data->GetPointData()->GetScalars(scaleFieldNameChar);
+            // values->GetRange(range);
 
-            // set range accoding to user-supplied limits
+            // Set range accoding to user-supplied limits
             range[0] = range_.first();
             range[1] = range_.second();
             glyph->ClampingOn();
             glyph->SetRange(range);
 
-            // if range[0] != min(value), maxGlyphLength behaviour will not
+            // If range[0] != min(value), maxGlyphLength behaviour will not
             // be correct...
             glyph->SetScaleFactor(maxGlyphLength);
         }
@@ -370,7 +404,7 @@ void Foam::fieldVisualisationBase::addGlyphs
             0,
             0,
             vtkDataObject::FIELD_ASSOCIATION_POINTS,
-            scaleFieldName.c_str()
+            scaleFieldNameChar
         );
     }
     else if (nComponents == 3)
@@ -388,21 +422,21 @@ void Foam::fieldVisualisationBase::addGlyphs
         if (maxGlyphLength > 0)
         {
             vtkDataArray* values =
-                data->GetPointData()->GetVectors(scaleFieldName.c_str());
+                data->GetPointData()->GetVectors(scaleFieldNameChar);
+
             double range[6];
             values->GetRange(range);
 
-/*
             // Attempt to set range for vectors...
-            scalar x0 = sqrt(sqr(range_.first())/3.0);
-            scalar x1 = sqrt(sqr(range_.second())/3.0);
-            range[0] = x0;
-            range[1] = x0;
-            range[2] = x0;
-            range[3] = x1;
-            range[4] = x1;
-            range[5] = x1;
-*/
+            // scalar x0 = sqrt(sqr(range_.first())/3.0);
+            // scalar x1 = sqrt(sqr(range_.second())/3.0);
+            // range[0] = x0;
+            // range[1] = x0;
+            // range[2] = x0;
+            // range[3] = x1;
+            // range[4] = x1;
+            // range[5] = x1;
+
             glyph->ClampingOn();
             glyph->SetRange(range);
             glyph->SetScaleFactor(maxGlyphLength);
@@ -421,7 +455,7 @@ void Foam::fieldVisualisationBase::addGlyphs
             0,
             0,
             vtkDataObject::FIELD_ASSOCIATION_POINTS,
-            scaleFieldName.c_str()
+            scaleFieldNameChar
         );
     }
     else
diff --git a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.H b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.H
index f7166470af..be2fe747ad 100644
--- a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.H
+++ b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -88,6 +88,8 @@ SourceFiles
 #include "OFstream.H"
 #include "Switch.H"
 
+#include "convectionScheme.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -140,6 +142,14 @@ class blendingFactor
         //- Disallow default bitwise assignment
         void operator=(const blendingFactor&);
 
+        //- Helper function to calculate the blending factor for the scheme
+        template<class Type>
+        void calcScheme
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& field,
+            const typename fv::convectionScheme<Type>& cs
+        );
+
         //- Calculate the blending factor
         template<class Type>
         void calc();
diff --git a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorTemplates.C b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorTemplates.C
index 26832a892d..7031607bd4 100644
--- a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorTemplates.C
+++ b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "gaussConvectionScheme.H"
+#include "boundedConvectionScheme.H"
 #include "blendedSchemeBase.H"
 #include "fvcCellReduce.H"
 #include "zeroGradientFvPatchFields.H"
@@ -31,39 +32,35 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::blendingFactor::calc()
+void Foam::blendingFactor::calcScheme
+(
+    const GeometricField<Type, fvPatchField, volMesh>& field,
+    const typename fv::convectionScheme<Type>& cs
+)
 {
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-
-    if (!obr_.foundObject<fieldType>(fieldName_))
+    if (!isA<fv::gaussConvectionScheme<Type>>(cs))
     {
+        WarningInFunction
+            << "Scheme for field " << field.name() << " is not a "
+            << fv::gaussConvectionScheme<Type>::typeName
+            << " scheme. Not calculating " << resultName_ << endl;
+
         return;
     }
 
-    const fvMesh& mesh = refCast<const fvMesh>(obr_);
-
-    const fieldType& field = mesh.lookupObject<fieldType>(fieldName_);
-
-    const word divScheme("div(" + phiName_ + ',' + fieldName_ + ')');
-    ITstream& its = mesh.divScheme(divScheme);
-
-    const surfaceScalarField& phi =
-        mesh.lookupObject<surfaceScalarField>(phiName_);
-
-    tmp<fv::convectionScheme<Type>> cs =
-        fv::convectionScheme<Type>::New(mesh, phi, its);
-
     const fv::gaussConvectionScheme<Type>& gcs =
-        refCast<const fv::gaussConvectionScheme<Type>>(cs());
+        refCast<const fv::gaussConvectionScheme<Type>>(cs);
 
-    const surfaceInterpolationScheme<Type>& interpScheme =
-        gcs.interpScheme();
+
+    const surfaceInterpolationScheme<Type>& interpScheme = gcs.interpScheme();
 
     if (!isA<blendedSchemeBase<Type>>(interpScheme))
     {
-        FatalErrorInFunction
-            << interpScheme.typeName << " is not a blended scheme"
-            << exit(FatalError);
+        WarningInFunction
+            << interpScheme.type() << " is not a blended scheme"
+            << ". Not calculating " << resultName_ << endl;
+
+        return;
     }
 
     // Retrieve the face-based blending factor
@@ -123,4 +120,46 @@ void Foam::blendingFactor::calc()
 }
 
 
+template<class Type>
+void Foam::blendingFactor::calc()
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (!obr_.foundObject<fieldType>(fieldName_))
+    {
+        return;
+    }
+
+    const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
+    const fieldType& field = mesh.lookupObject<fieldType>(fieldName_);
+
+    const word divScheme("div(" + phiName_ + ',' + fieldName_ + ')');
+    ITstream& its = mesh.divScheme(divScheme);
+
+    const surfaceScalarField& phi =
+        mesh.lookupObject<surfaceScalarField>(phiName_);
+
+    tmp<fv::convectionScheme<Type>> tcs =
+        fv::convectionScheme<Type>::New(mesh, phi, its);
+
+    const fv::convectionScheme<Type>& cs = tcs();
+
+    if (isA<fv::boundedConvectionScheme<Type>>(cs))
+    {
+        const fv::boundedConvectionScheme<Type>& bcs =
+            refCast<const fv::boundedConvectionScheme<Type>>(cs);
+
+        calcScheme(field, bcs.scheme());
+    }
+    else
+    {
+        const fv::gaussConvectionScheme<Type>& gcs =
+            refCast<const fv::gaussConvectionScheme<Type>>(cs);
+
+        calcScheme(field, gcs);
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/utilities/mapFields/mapFields.H b/src/postProcessing/functionObjects/utilities/mapFields/mapFields.H
index 80c2ebe119..9343616cfe 100644
--- a/src/postProcessing/functionObjects/utilities/mapFields/mapFields.H
+++ b/src/postProcessing/functionObjects/utilities/mapFields/mapFields.H
@@ -122,6 +122,13 @@ class mapFieldsFO
         //- Helper function to create the mesh-to-mesh interpolation
         void createInterpolation(const dictionary& dict);
 
+        //- Helper function to evaluate constraint patches after mapping
+        template<class Type>
+        void evaluateConstraintTypes
+        (
+            GeometricField<Type, fvPatchField, volMesh>& fld
+        ) const;
+
         //- Helper function to interpolate and write the fied
         template<class Type>
         bool writeFieldType() const;
diff --git a/src/postProcessing/functionObjects/utilities/mapFields/mapFieldsTemplates.C b/src/postProcessing/functionObjects/utilities/mapFields/mapFieldsTemplates.C
index de269ea867..dd1086049b 100644
--- a/src/postProcessing/functionObjects/utilities/mapFields/mapFieldsTemplates.C
+++ b/src/postProcessing/functionObjects/utilities/mapFields/mapFieldsTemplates.C
@@ -25,6 +25,91 @@ License
 
 #include "meshToMesh.H"
 
+template<class Type>
+void Foam::mapFieldsFO::evaluateConstraintTypes
+(
+    GeometricField<Type, fvPatchField, volMesh>& fld
+) const
+{
+    typename GeometricField<Type, fvPatchField, volMesh>::
+        GeometricBoundaryField& fldBf = fld.boundaryField();
+
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
+    {
+        label nReq = Pstream::nRequests();
+
+        forAll(fldBf, patchi)
+        {
+            fvPatchField<Type>& tgtField = fldBf[patchi];
+
+            if
+            (
+                tgtField.type() == tgtField.patch().patch().type()
+             && polyPatch::constraintType(tgtField.patch().patch().type())
+            )
+            {
+                tgtField.initEvaluate(Pstream::defaultCommsType);
+            }
+        }
+
+        // Block for any outstanding requests
+        if
+        (
+            Pstream::parRun()
+         && Pstream::defaultCommsType == Pstream::nonBlocking
+        )
+        {
+            Pstream::waitRequests(nReq);
+        }
+
+        forAll(fldBf, patchi)
+        {
+            fvPatchField<Type>& tgtField = fldBf[patchi];
+
+            if
+            (
+                tgtField.type() == tgtField.patch().patch().type()
+             && polyPatch::constraintType(tgtField.patch().patch().type())
+            )
+            {
+                tgtField.evaluate(Pstream::defaultCommsType);
+            }
+        }
+    }
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
+    {
+        const lduSchedule& patchSchedule =
+            fld.mesh().globalData().patchSchedule();
+
+        forAll(patchSchedule, patchEvali)
+        {
+            label patchi = patchSchedule[patchEvali].patch;
+            fvPatchField<Type>& tgtField = fldBf[patchi];
+
+            if
+            (
+                tgtField.type() == tgtField.patch().patch().type()
+             && polyPatch::constraintType(tgtField.patch().patch().type())
+            )
+            {
+                if (patchSchedule[patchEvali].init)
+                {
+                    tgtField.initEvaluate(Pstream::scheduled);
+                }
+                else
+                {
+                    tgtField.evaluate(Pstream::scheduled);
+                }
+            }
+        }
+    }
+}
+
+
 template<class Type>
 bool Foam::mapFieldsFO::writeFieldType() const
 {
@@ -53,6 +138,9 @@ bool Foam::mapFieldsFO::writeFieldType() const
         if (log_) Info<< ": interpolated";
 
         FieldType fieldMapRegion(mapRegionIO, tfieldMapRegion);
+
+        evaluateConstraintTypes(fieldMapRegion);
+
         fieldMapRegion.write();
 
         if (log_) Info<< " and written" << nl;
diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.C b/src/randomProcesses/noise/noiseFFT/noiseFFT.C
index 3a61c9b0dd..e82d25b90f 100644
--- a/src/randomProcesses/noise/noiseFFT/noiseFFT.C
+++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.C
@@ -127,10 +127,13 @@ void Foam::noiseFFT::octaveBandInfo
 
     fBandIDs = bandIDs.toc();
 
-    // Remove the last centre frequency (beyond upper frequency limit)
-    fc.remove();
+    if (fc.size())
+    {
+        // Remove the last centre frequency (beyond upper frequency limit)
+        fc.remove();
 
-    fCentre.transfer(fc);
+        fCentre.transfer(fc);
+    }
 }
 
 
diff --git a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C
index b8e312c0ce..92a3df57fd 100644
--- a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C
+++ b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C
@@ -494,8 +494,21 @@ void surfaceNoise::calculate()
         octave13FreqCentre
     );
 
-    List<scalarField> surfPSD13f(octave13BandIDs.size() - 1);
-    List<scalarField> surfPrms13f2(octave13BandIDs.size() - 1);
+    label bandSize = 0;
+    if (octave13BandIDs.empty())
+    {
+        WarningInFunction
+            << "Ocatve band calculation failed (zero sized). "
+            << "please check your input data"
+            << endl;
+    }
+    else
+    {
+        bandSize = octave13BandIDs.size() - 1;
+    }
+
+    List<scalarField> surfPSD13f(bandSize);
+    List<scalarField> surfPrms13f2(bandSize);
     forAll(surfPSD13f, freqI)
     {
         surfPSD13f[freqI].setSize(nLocalFace);
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
index af8bf5563a..2d3535d873 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
@@ -335,6 +335,10 @@ void reactingOneDim::solveEnergy()
         hEqn += fvc::div(phiQr);
     }
 
+/* NOTE: The moving mesh option is only correct for reaction such as
+   Solid -> Gas, thus the ddt term is compesated exaclty by chemistrySh and
+   the mesh flux is not necessary.
+   
     if (regionMesh().moving())
     {
         surfaceScalarField phihMesh
@@ -344,7 +348,7 @@ void reactingOneDim::solveEnergy()
 
         hEqn -= fvc::div(phihMesh);
     }
-
+*/
     hEqn.relax();
     hEqn.solve();
 }
diff --git a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.C b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.C
index 3e1b422bde..f90d05b85d 100644
--- a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.C
+++ b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.C
@@ -34,7 +34,7 @@ namespace Foam
 }
 
 
-void Foam::ensightSurfaceReader::skip(const label n, IFstream& is) const
+void Foam::ensightSurfaceReader::skip(const label n, Istream& is) const
 {
     label i = 0;
     token t;
@@ -59,6 +59,40 @@ void Foam::ensightSurfaceReader::skip(const label n, IFstream& is) const
 }
 
 
+void Foam::ensightSurfaceReader::readLine(IFstream& is, string& buffer) const
+{
+    buffer = "";
+    while (is.good() && buffer == "")
+    {
+        is.getLine(buffer);
+    }
+}
+
+
+void Foam::ensightSurfaceReader::debugSection
+(
+    const word& expected,
+    IFstream& is
+) const
+{
+    string actual = "";
+    readLine(is, actual);
+
+    if (expected != actual)
+    {
+        FatalIOErrorInFunction(is)
+            << "Expected section header '" << expected
+            << "' but read the word '" << actual << "'"
+            << exit(FatalIOError);
+    }
+
+    if (debug)
+    {
+        Info<< "Read section header: " << expected << endl;
+    }
+}
+
+
 void Foam::ensightSurfaceReader::readGeometryHeader(ensightReadFile& is) const
 {
     // Binary flag string if applicable
@@ -101,29 +135,6 @@ void Foam::ensightSurfaceReader::readGeometryHeader(ensightReadFile& is) const
 }
 
 
-void Foam::ensightSurfaceReader::debugSection
-(
-    const word& expected,
-    IFstream& is
-) const
-{
-    word actual(is);
-
-    if (expected != actual)
-    {
-        FatalIOErrorInFunction(is)
-            << "Expected section header '" << expected
-            << "' but read the word '" << actual << "'"
-            << exit(FatalIOError);
-    }
-
-    if (debug)
-    {
-        Info<< "Read section header: " << expected << endl;
-    }
-}
-
-
 void Foam::ensightSurfaceReader::readCase(IFstream& is)
 {
     if (debug)
@@ -141,31 +152,54 @@ void Foam::ensightSurfaceReader::readCase(IFstream& is)
     string buffer;
 
     // Read the file
+
     debugSection("FORMAT", is);
-    skip(3, is); // type: ensight gold
+    readLine(is, buffer);  // type: ensight gold
 
     debugSection("GEOMETRY", is);
-    readSkip(is, 2, meshFileName_);
+    readLine(is, buffer);
+    readFromLine(2, buffer, meshFileName_); // model: 1 xxx.0000.mesh
 
     debugSection("VARIABLE", is);
 
+    // Read the field description
     DynamicList<word> fieldNames(10);
-    DynamicList<word> fieldFileNames(10);
+    DynamicList<string> fieldFileNames(10);
     word fieldName;
-    word fieldFileName;
+    string fieldFileName;
     while (is.good())
     {
-        word primitiveType(is); // scalar, vector
+        readLine(is, buffer);
 
-        if (primitiveType == "TIME")
+        if (buffer == "TIME")
         {
             break;
         }
 
-        readSkip(is, 3, fieldName); // p, U etc
+        IStringStream iss(buffer);
+
+        // Read the field name, e.g. p U etc
+        readFromLine(4, iss, fieldName);
         fieldNames.append(fieldName);
 
-        is >> fieldFileName; // surfaceName.****.fieldName
+        // Field file name may contain /'s e.g.
+        //   surfaceName.****.fieldName
+        // This is not parser friendly - simply take remainder of buffer
+        label iPos = iss.stdStream().tellg();
+        fieldFileName = buffer(iPos, buffer.size() - iPos);
+        size_t p0  = fieldFileName.find_first_not_of(' ');
+        if (p0 == string::npos)
+        {
+            WarningInFunction
+                << "Error reading field file name. "
+                << "Current buffer: " << buffer
+                << endl;
+        }
+        else
+        {
+            size_t p1  = fieldFileName.find_last_not_of(' ');
+            fieldFileName = fieldFileName.substr(p0, p1 - p0 + 1);
+        }
         fieldFileNames.append(fieldFileName);
     }
     fieldNames_.transfer(fieldNames);
@@ -178,10 +212,14 @@ void Foam::ensightSurfaceReader::readCase(IFstream& is)
     }
 
     // Start reading time information
-    skip(3, is); // time set: 1
-    readSkip(is, 3, nTimeSteps_);
-    readSkip(is, 3, timeStartIndex_);
-    readSkip(is, 2, timeIncrement_);
+    readLine(is, buffer); // time set: 1
+
+    readLine(is, buffer);
+    readFromLine(3, buffer, nTimeSteps_);
+    readLine(is, buffer);
+    readFromLine(3, buffer, timeStartIndex_);
+    readLine(is, buffer);
+    readFromLine(2, buffer, timeIncrement_);
 
     if (debug)
     {
@@ -191,7 +229,7 @@ void Foam::ensightSurfaceReader::readCase(IFstream& is)
     }
 
     // Read the time values
-    skip(2, is);
+    readLine(is, buffer); // time values:
     timeValues_.setSize(nTimeSteps_);
     for (label i = 0; i < nTimeSteps_; i++)
     {
diff --git a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.H b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.H
index 8c99cfd466..0ba7e1c49b 100644
--- a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.H
+++ b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.H
@@ -69,7 +69,7 @@ protected:
         List<word> fieldNames_;
 
         //- Field file names
-        List<word> fieldFileNames_;
+        List<string> fieldFileNames_;
 
         //- Number of time steps
         label nTimeSteps_;
@@ -91,21 +91,38 @@ protected:
 
     // Protected Member Functions
 
+        //- Helper function to skip forward n steps in stream
+        void skip(const label n, Istream& is) const;
+
+        //- Helper function to read an ascii line from file
+        void readLine(IFstream& is, string& buffer) const;
+
         //- Read and check a section header
         void debugSection(const word& expected, IFstream& is) const;
 
-        //- Read the case file
-        void readCase(IFstream& is);
-
-        //- Helper function to skip forward n steps in stream
-        void skip(const label n, IFstream& is) const;
-
         //- Read (and throw away) geometry file header
         void readGeometryHeader(ensightReadFile& is) const;
 
+        //- Read the case file
+        void readCase(IFstream& is);
+
         //- Helper function to return Type after skipping n tokens
         template<class Type>
-        void readSkip(IFstream& is, const label nSkip, Type& value) const;
+        void readFromLine
+        (
+            const label nSkip,
+            IStringStream& is,
+            Type& value
+        ) const;
+
+        //- Helper function to return Type after skipping n tokens
+        template<class Type>
+        void readFromLine
+        (
+            const label nSkip,
+            const string& buffer,
+            Type& value
+        ) const;
 
         //- Helper function to return a field
         template<class Type>
diff --git a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReaderTemplates.C b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReaderTemplates.C
index 04a593575a..dace070320 100644
--- a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReaderTemplates.C
+++ b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReaderTemplates.C
@@ -29,10 +29,10 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::ensightSurfaceReader::readSkip
+void Foam::ensightSurfaceReader::readFromLine
 (
-    IFstream& is,
     const label nSkip,
+    IStringStream& is,
     Type& value
 ) const
 {
@@ -42,6 +42,20 @@ void Foam::ensightSurfaceReader::readSkip
 }
 
 
+template<class Type>
+void Foam::ensightSurfaceReader::readFromLine
+(
+    const label nSkip,
+    const string& buffer,
+    Type& value
+) const
+{
+    IStringStream is(buffer);
+
+    readFromLine(nSkip, is, value);
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::ensightSurfaceReader::readField
 (
@@ -60,9 +74,19 @@ Foam::tmp<Foam::Field<Type> > Foam::ensightSurfaceReader::readField
     fileName fieldFileName(fieldFileNames_[fieldIndex]);
 
     std::ostringstream oss;
-    oss << std::setfill('0') << std::setw(4) << fileIndex;
+    label nMask = 0;
+    for (size_t chari = 0; chari < fieldFileName.size(); chari++)
+    {
+        if (fieldFileName[chari] == '*')
+        {
+            nMask++;
+        }
+    }
+
+    const std::string maskStr(nMask, '*');
+    oss << std::setfill('0') << std::setw(nMask) << fileIndex;
     const word indexStr = oss.str();
-    fieldFileName.replace("****", indexStr);
+    fieldFileName.replace(maskStr, indexStr);
 
     
     ensightReadFile is(baseDir_/fieldFileName, streamFormat_);
diff --git a/src/thermophysicalModels/doc/thermophysicalModels.dox b/src/thermophysicalModels/doc/thermophysicalModels.dox
index a8f7f48ae9..9df1640a0b 100644
--- a/src/thermophysicalModels/doc/thermophysicalModels.dox
+++ b/src/thermophysicalModels/doc/thermophysicalModels.dox
@@ -23,7 +23,7 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-\page pageThermophsyicalModels Thermophsyical Models
+\page pageThermophsyicalModels Thermophysical Models
 
 \section secSchemes Overview
 The available thermophysical models are grouped into the following categories:
diff --git a/tutorials/incompressible/icoFoam/cavityMappingTest/Allrun-serial b/tutorials/incompressible/icoFoam/cavityMappingTest/Allrun
similarity index 100%
rename from tutorials/incompressible/icoFoam/cavityMappingTest/Allrun-serial
rename to tutorials/incompressible/icoFoam/cavityMappingTest/Allrun
diff --git a/tutorials/mesh/parallel/filter/Allrun b/tutorials/mesh/parallel/filter/Allrun
index 8a0978f43f..4be640030b 100755
--- a/tutorials/mesh/parallel/filter/Allrun
+++ b/tutorials/mesh/parallel/filter/Allrun
@@ -4,6 +4,8 @@ cd ${0%/*} || exit 1    # Run from this directory
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
+application=$(getApplication)
+
 # Create mesh
 runApplication blockMesh
 
@@ -16,7 +18,7 @@ runApplication topoSet
 # Create baffles and fields
 runApplication createBaffles -overwrite
 
-runApplication $(getApplication)
+runApplication $application
 
 #- RedistributePar to do decomposition
 runParallel redistributePar -decompose -cellDist