diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
index 1c4aa0a1..d9345cd9 100644
--- a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
@@ -19,16 +19,7 @@
             ).fvmDiv(phi, alpha1)
         );
 
-        solve
-        (
-            alpha1Eqn
-          - fv::gaussLaplacianScheme<scalar, scalar>
-            (
-                mesh,
-                linear<scalar>(mesh),
-                fv::uncorrectedSnGrad<scalar>(mesh)
-            ).fvmLaplacian(fvc::interpolate(turbulence->nut()), alpha1)
-        );
+        solve(alpha1Eqn);
 
         Info<< "Phase-1 volume fraction = "
             << alpha1.weightedAverage(mesh.Vsc()).value()
@@ -42,6 +33,7 @@
         if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
         {
             Info<< "Applying the previous iteration correction flux" << endl;
+
             #ifdef LTSSOLVE
             MULES::LTScorrect
             (
diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
index 8ad3234a..4d889f12 100644
--- a/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
@@ -50,7 +50,6 @@
 
     // Apply the diffusion term separately to allow implicit solution
     // and boundedness of the explicit advection
-    if (!MULESCorr)
     {
         fvScalarMatrix alpha1Eqn
         (
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 86983f2e..82cf0515 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -323,6 +323,9 @@ addLayersControls
     // - overall thickness and firstLayerThickness
     // - overall thickness and finalLayerThickness
     // - overall thickness and expansionRatio
+    //
+    // Note: the mode thus selected is global, i.e. one cannot override the
+    //       mode on a per-patch basis (only the values can be overridden)
 
         // Expansion factor for layer mesh
         expansionRatio 1.0;
diff --git a/applications/utilities/surface/surfaceCheck/surfaceCheck.C b/applications/utilities/surface/surfaceCheck/surfaceCheck.C
index fa3fc28c..8116ab09 100644
--- a/applications/utilities/surface/surfaceCheck/surfaceCheck.C
+++ b/applications/utilities/surface/surfaceCheck/surfaceCheck.C
@@ -636,7 +636,7 @@ int main(int argc, char *argv[])
                     faces,
                     "zone",
                     scalarFaceZone,
-                    true
+                    false               // face based data
                 );
             }
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C
index 8ac95cd8..50f01629 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,9 +45,6 @@ Foam::PatchTools::pointNormals
     const PrimitivePatch<Face, FaceList, PointField, PointType>& p
 )
 {
-    // Assume patch is smaller than the globalData().coupledPatch() (?) so
-    // loop over patch meshPoints.
-
     const globalMeshData& globalData = mesh.globalData();
     const indirectPrimitivePatch& coupledPatch = globalData.coupledPatch();
     const Map<label>& coupledPatchMP = coupledPatch.meshPointMap();
@@ -56,6 +53,108 @@ Foam::PatchTools::pointNormals
         globalData.globalTransforms();
 
 
+
+
+    // Combine normals. Note: do on all master points. Cannot just use
+    // patch points since the master point does not have to be on the
+    // patch!
+
+    pointField coupledPointNormals(map.constructSize(), vector::zero);
+
+    {
+        // Collect local pointFaces (sized on patch points only)
+        List<List<point> > pointFaceNormals(map.constructSize());
+        forAll(p.meshPoints(), patchPointI)
+        {
+            label meshPointI = p.meshPoints()[patchPointI];
+            Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
+            if (fnd != coupledPatchMP.end())
+            {
+                label coupledPointI = fnd();
+
+                List<point>& pNormals = pointFaceNormals[coupledPointI];
+                const labelList& pFaces = p.pointFaces()[patchPointI];
+                pNormals.setSize(pFaces.size());
+                forAll(pFaces, i)
+                {
+                    pNormals[i] = p.faceNormals()[pFaces[i]];
+                }
+            }
+        }
+
+
+        // Pull remote data into local slots
+        map.distribute
+        (
+            transforms,
+            pointFaceNormals,
+            mapDistribute::transform()
+        );
+
+
+        // Combine all face normals (-local, -remote,untransformed,
+        //  -remote,transformed)
+
+        const labelListList& slaves = globalData.globalPointSlaves();
+        const labelListList& transformedSlaves =
+            globalData.globalPointTransformedSlaves();
+
+        forAll(slaves, coupledPointI)
+        {
+            const labelList& slaveSlots = slaves[coupledPointI];
+            const labelList& transformedSlaveSlots =
+                transformedSlaves[coupledPointI];
+
+            point& n = coupledPointNormals[coupledPointI];
+
+            // Local entries
+            const List<point>& local = pointFaceNormals[coupledPointI];
+
+            label nFaces =
+                local.size()
+              + slaveSlots.size()
+              + transformedSlaveSlots.size();
+
+            n = sum(local);
+
+            // Add any remote face normals
+            forAll(slaveSlots, i)
+            {
+                n += sum(pointFaceNormals[slaveSlots[i]]);
+            }
+            forAll(transformedSlaveSlots, i)
+            {
+                n += sum(pointFaceNormals[transformedSlaveSlots[i]]);
+            }
+
+            if (nFaces >= 1)
+            {
+                n /= mag(n)+VSMALL;
+            }
+
+            // Put back into slave slots
+            forAll(slaveSlots, i)
+            {
+                coupledPointNormals[slaveSlots[i]] = n;
+            }
+            forAll(transformedSlaveSlots, i)
+            {
+                coupledPointNormals[transformedSlaveSlots[i]] = n;
+            }
+        }
+
+
+        // Send back
+        map.reverseDistribute
+        (
+            transforms,
+            coupledPointNormals.size(),
+            coupledPointNormals,
+            mapDistribute::transform()
+        );
+    }
+
+
     // 1. Start off with local normals (note:without calculating pointNormals
     //    to avoid them being stored)
 
@@ -78,99 +177,7 @@ Foam::PatchTools::pointNormals
     }
 
 
-    // Collect local pointFaces
-    List<List<point> > pointFaceNormals(map.constructSize());
-    forAll(p.meshPoints(), patchPointI)
-    {
-        label meshPointI = p.meshPoints()[patchPointI];
-        Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
-        if (fnd != coupledPatchMP.end())
-        {
-            label coupledPointI = fnd();
-
-            List<point>& pNormals = pointFaceNormals[coupledPointI];
-            const labelList& pFaces = p.pointFaces()[patchPointI];
-            pNormals.setSize(pFaces.size());
-            forAll(pFaces, i)
-            {
-                pNormals[i] = p.faceNormals()[pFaces[i]];
-            }
-        }
-    }
-
-
-    // Pull remote data into local slots
-    map.distribute
-    (
-        transforms,
-        pointFaceNormals,
-        mapDistribute::transform()
-    );
-
-
-    // Combine normals
-    const labelListList& slaves = globalData.globalPointSlaves();
-    const labelListList& transformedSlaves =
-        globalData.globalPointTransformedSlaves();
-
-
-    pointField coupledPointNormals(map.constructSize(), vector::zero);
-
-    forAll(p.meshPoints(), patchPointI)
-    {
-        label meshPointI = p.meshPoints()[patchPointI];
-        Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
-        if (fnd != coupledPatchMP.end())
-        {
-            label coupledPointI = fnd();
-            const labelList& slaveSlots =
-                slaves[coupledPointI];
-            const labelList& transformedSlaveSlots =
-                transformedSlaves[coupledPointI];
-
-            label nFaces = slaveSlots.size()+transformedSlaveSlots.size();
-            if (nFaces > 0)
-            {
-                // Combine
-                point& n = coupledPointNormals[coupledPointI];
-
-                n += sum(pointFaceNormals[coupledPointI]);
-
-                forAll(slaveSlots, i)
-                {
-                    n += sum(pointFaceNormals[slaveSlots[i]]);
-                }
-                forAll(transformedSlaveSlots, i)
-                {
-                    n += sum(pointFaceNormals[transformedSlaveSlots[i]]);
-                }
-                n /= mag(n)+VSMALL;
-
-                // Put back into slave slots
-                forAll(slaveSlots, i)
-                {
-                    coupledPointNormals[slaveSlots[i]] = n;
-                }
-                forAll(transformedSlaveSlots, i)
-                {
-                    coupledPointNormals[transformedSlaveSlots[i]] = n;
-                }
-            }
-        }
-    }
-
-
-    // Send back
-    map.reverseDistribute
-    (
-        transforms,
-        coupledPointNormals.size(),
-        coupledPointNormals,
-        mapDistribute::transform()
-    );
-
-
-    // Override patch normals
+    // 2. Override patch normals on coupled points
     forAll(p.meshPoints(), patchPointI)
     {
         label meshPointI = p.meshPoints()[patchPointI];
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
index b817777a..d546ff86 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -528,6 +528,16 @@ void Foam::motionSmootherAlgo::setDisplacement
     }
 
 
+    // Combine any coupled points
+    syncTools::syncPointList
+    (
+        mesh,
+        displacement,
+        maxMagEqOp(),           // combine op
+        vector::zero            // null value
+    );
+
+
     // Adapt the fixedValue bc's (i.e. copy internal point data to
     // boundaryField for all affected patches)
     setDisplacementPatchFields(patchIDs, displacement);
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
index e41c9bec..090c096d 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,6 +107,7 @@ void Foam::MULES::correct
         psiMax, psiMin,
         nLimiterIter
     );
+
     correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su);
 }
 
@@ -172,6 +173,13 @@ void Foam::MULES::limiterCorr
 
     const fvMesh& mesh = psi.mesh();
 
+    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
+
+    scalar extremaCoeff
+    (
+        MULEScontrols.lookupOrDefault<scalar>("extremaCoeff", 0.0)
+    );
+
     const labelUList& owner = mesh.owner();
     const labelUList& neighb = mesh.neighbour();
     tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
@@ -283,8 +291,8 @@ void Foam::MULES::limiterCorr
         }
     }
 
-    psiMaxn = min(psiMaxn, psiMax);
-    psiMinn = max(psiMinn, psiMin);
+    psiMaxn = min(psiMaxn + extremaCoeff*(psiMax - psiMin), psiMax);
+    psiMinn = max(psiMinn - extremaCoeff*(psiMax - psiMin), psiMin);
 
     // scalar smooth = 0.5;
     // psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
diff --git a/src/mesh/autoMesh/Make/files b/src/mesh/autoMesh/Make/files
index be8a5fff..2bc9cd80 100644
--- a/src/mesh/autoMesh/Make/files
+++ b/src/mesh/autoMesh/Make/files
@@ -2,7 +2,7 @@ autoHexMesh             = autoHexMesh
 autoHexMeshDriver       = $(autoHexMesh)/autoHexMeshDriver
 
 $(autoHexMeshDriver)/autoLayerDriver.C
-$(autoHexMeshDriver)/autoLayerDriverShrink.C
+/* $(autoHexMeshDriver)/autoLayerDriverShrink.C */
 $(autoHexMeshDriver)/autoSnapDriver.C
 $(autoHexMeshDriver)/autoSnapDriverFeature.C
 $(autoHexMeshDriver)/autoRefineDriver.C
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index ee94b0b0..07bf19ce 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,7 +29,6 @@ Description
 
 SourceFiles
     autoLayerDriver.C
-    autoLayerDriverShrink.C
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index ed05edf4..006ead54 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -1292,7 +1292,15 @@ void Foam::autoSnapDriver::detectNearSurfaces
     }
 
 
-    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
+    const PackedBoolList isPatchMasterPoint
+    (
+        meshRefinement::getMasterPoints
+        (
+            mesh,
+            meshPoints
+        )
+    );
+
     label nOverride = 0;
 
     // 1. All points to non-interface surfaces
@@ -1405,7 +1413,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
                 }
             }
 
-            if (override && isMasterPoint[meshPoints[pointI]])
+            if (override && isPatchMasterPoint[pointI])
             {
                 nOverride++;
             }
@@ -1543,7 +1551,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
                     }
                 }
 
-                if (override && isMasterPoint[meshPoints[pointI]])
+                if (override && isPatchMasterPoint[pointI])
                 {
                     nOverride++;
                 }
@@ -1758,16 +1766,19 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
         }
 
         {
+            const PackedBoolList isPatchMasterPoint
+            (
+                meshRefinement::getMasterPoints
+                (
+                    mesh,
+                    pp.meshPoints()
+                )
+            );
+
             scalarField magDisp(mag(patchDisp));
 
             Info<< "Wanted displacement : average:"
-                <<  meshRefinement::gAverage
-                    (
-                        mesh,
-                        syncTools::getMasterPoints(mesh),
-                        pp.meshPoints(),
-                        magDisp
-                    )
+                <<  meshRefinement::gAverage(isPatchMasterPoint, magDisp)
                 << " min:" << gMin(magDisp)
                 << " max:" << gMax(magDisp) << endl;
         }
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index c7a00df0..b1cc58a0 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -126,7 +126,7 @@ bool Foam::autoSnapDriver::isFeaturePoint
 
 void Foam::autoSnapDriver::smoothAndConstrain
 (
-    const PackedBoolList& isMasterEdge,
+    const PackedBoolList& isPatchMasterEdge,
     const indirectPrimitivePatch& pp,
     const labelList& meshEdges,
     const List<pointConstraint>& constraints,
@@ -168,7 +168,7 @@ void Foam::autoSnapDriver::smoothAndConstrain
                 {
                     label edgeI = pEdges[i];
 
-                    if (isMasterEdge[meshEdges[edgeI]])
+                    if (isPatchMasterEdge[edgeI])
                     {
                         label nbrPointI = edges[edgeI].otherVertex(pointI);
                         if (constraints[nbrPointI].first() >= nConstraints)
@@ -3231,20 +3231,24 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         patchConstraints
     );
 
-    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
+    //const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
+    const PackedBoolList isPatchMasterPoint
+    (
+        meshRefinement::getMasterPoints
+        (
+            mesh,
+            pp.meshPoints()
+        )
+    );
     {
         vector avgPatchDisp = meshRefinement::gAverage
         (
-            mesh,
-            isMasterPoint,
-            pp.meshPoints(),
+            isPatchMasterPoint,
             patchDisp
         );
         vector avgPatchAttr = meshRefinement::gAverage
         (
-            mesh,
-            isMasterPoint,
-            pp.meshPoints(),
+            isPatchMasterPoint,
             patchAttraction
         );
 
@@ -3280,8 +3284,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
     // Count
     {
-        const labelList& meshPoints = pp.meshPoints();
-
         label nMasterPoints = 0;
         label nPlanar = 0;
         label nEdge = 0;
@@ -3289,7 +3291,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
         forAll(patchConstraints, pointI)
         {
-            if (isMasterPoint[meshPoints[pointI]])
+            if (isPatchMasterPoint[pointI])
             {
                 nMasterPoints++;
 
@@ -3337,7 +3339,19 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
     if (featureAttract < 1-0.001)
     {
-        const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
+        //const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
+        const labelList meshEdges
+        (
+            pp.meshEdges(mesh.edges(), mesh.pointEdges())
+        );
+        const PackedBoolList isPatchMasterEdge
+        (
+            meshRefinement::getMasterEdges
+            (
+                mesh,
+                meshEdges
+            )
+        );
 
         const vectorField pointNormals
         (
@@ -3347,17 +3361,13 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
                 pp
             )
         );
-        const labelList meshEdges
-        (
-            pp.meshEdges(mesh.edges(), mesh.pointEdges())
-        );
 
 
         // 1. Smoothed all displacement
         vectorField smoothedPatchDisp = patchDisp;
         smoothAndConstrain
         (
-            isMasterEdge,
+            isPatchMasterEdge,
             pp,
             meshEdges,
             patchConstraints,
@@ -3370,7 +3380,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
         smoothAndConstrain
         (
-            isMasterEdge,
+            isPatchMasterEdge,
             pp,
             meshEdges,
             patchConstraints,
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C
index 55e60904..2067b24c 100644
--- a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -49,7 +49,6 @@ Foam::externalDisplacementMeshMover::externalDisplacementMeshMover
 {}
 
 
-
 // * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
 
 Foam::autoPtr<Foam::externalDisplacementMeshMover>
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
index b605ae67..df17b0ec 100644
--- a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
@@ -129,9 +129,8 @@ Foam::medialAxisMeshMover::getPatch
 void Foam::medialAxisMeshMover::smoothPatchNormals
 (
     const label nSmoothDisp,
-    const PackedBoolList& isMasterPoint,
-    const PackedBoolList& isMasterEdge,
-    const labelList& meshEdges,
+    const PackedBoolList& isPatchMasterPoint,
+    const PackedBoolList& isPatchMasterEdge,
     pointField& normals
 ) const
 {
@@ -142,13 +141,12 @@ void Foam::medialAxisMeshMover::smoothPatchNormals
     // Get smoothly varying internal normals field.
     Info<< typeName << " : Smoothing normals ..." << endl;
 
-    scalarField edgeWeights(meshEdges.size());
+    scalarField edgeWeights(edges.size());
     scalarField invSumWeight(meshPoints.size());
     meshRefinement::calculateEdgeWeights
     (
         mesh(),
-        isMasterEdge,
-        meshEdges,
+        isPatchMasterEdge,
         meshPoints,
         edges,
         edgeWeights,
@@ -162,8 +160,7 @@ void Foam::medialAxisMeshMover::smoothPatchNormals
         meshRefinement::weightedSum
         (
             mesh(),
-            isMasterEdge,
-            meshEdges,
+            isPatchMasterEdge,
             meshPoints,
             edges,
             edgeWeights,
@@ -177,9 +174,7 @@ void Foam::medialAxisMeshMover::smoothPatchNormals
         {
             scalar resid = meshRefinement::gAverage
             (
-                mesh(),
-                isMasterPoint,
-                meshPoints,
+                isPatchMasterPoint,
                 mag(normals-average)()
             );
             Info<< "    Iteration " << iter << "   residual " << resid << endl;
@@ -201,8 +196,8 @@ void Foam::medialAxisMeshMover::smoothPatchNormals
 void Foam::medialAxisMeshMover::smoothNormals
 (
     const label nSmoothDisp,
-    const PackedBoolList& isMasterPoint,
-    const PackedBoolList& isMasterEdge,
+    const PackedBoolList& isMeshMasterPoint,
+    const PackedBoolList& isMeshMasterEdge,
     const labelList& fixedPoints,
     pointVectorField& normals
 ) const
@@ -229,18 +224,16 @@ void Foam::medialAxisMeshMover::smoothNormals
 
 
     // Correspondence between local edges/points and mesh edges/points
-    const labelList meshEdges(identity(mesh().nEdges()));
     const labelList meshPoints(identity(mesh().nPoints()));
 
     // Calculate inverse sum of weights
 
-    scalarField edgeWeights(meshEdges.size());
+    scalarField edgeWeights(mesh().nEdges());
     scalarField invSumWeight(meshPoints.size());
     meshRefinement::calculateEdgeWeights
     (
         mesh(),
-        isMasterEdge,
-        meshEdges,
+        isMeshMasterEdge,
         meshPoints,
         edges,
         edgeWeights,
@@ -253,8 +246,7 @@ void Foam::medialAxisMeshMover::smoothNormals
         meshRefinement::weightedSum
         (
             mesh(),
-            isMasterEdge,
-            meshEdges,
+            isMeshMasterEdge,
             meshPoints,
             edges,
             edgeWeights,
@@ -268,8 +260,7 @@ void Foam::medialAxisMeshMover::smoothNormals
         {
             scalar resid = meshRefinement::gAverage
             (
-                mesh(),
-                isMasterPoint,
+                isMeshMasterPoint,
                 mag(normals-average)()
             );
             Info<< "    Iteration " << iter << "   residual " << resid << endl;
@@ -411,9 +402,9 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
     // Predetermine mesh edges
     // ~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Precalulate master point/edge (only relevant for shared points/edges)
-    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh()));
-    const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh()));
+    // Precalulate (mesh) master point/edge (only relevant for shared pts/edges)
+    const PackedBoolList isMeshMasterPoint(syncTools::getMasterPoints(mesh()));
+    const PackedBoolList isMeshMasterEdge(syncTools::getMasterEdges(mesh()));
     // Precalculate meshEdge per pp edge
     const labelList meshEdges
     (
@@ -424,6 +415,23 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
         )
     );
 
+    // Precalulate (patch) master point/edge
+    const PackedBoolList isPatchMasterPoint
+    (
+        meshRefinement::getMasterPoints
+        (
+            mesh(),
+            meshPoints
+        )
+    );
+    const PackedBoolList isPatchMasterEdge
+    (
+        meshRefinement::getMasterEdges
+        (
+            mesh(),
+            meshEdges
+        )
+    );
 
     // Determine pointNormal
     // ~~~~~~~~~~~~~~~~~~~~~
@@ -434,9 +442,8 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
     smoothPatchNormals
     (
         nSmoothSurfaceNormals,
-        isMasterPoint,
-        isMasterEdge,
-        meshEdges,
+        isPatchMasterPoint,
+        isPatchMasterEdge,
         pointNormals
     );
 
@@ -795,8 +802,8 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
     smoothNormals
     (
         nSmoothNormals,
-        isMasterPoint,
-        isMasterEdge,
+        isMeshMasterPoint,
+        isMeshMasterEdge,
         meshPoints,
         dispVec_
     );
@@ -999,9 +1006,8 @@ void Foam::medialAxisMeshMover::syncPatchDisplacement
 void Foam::medialAxisMeshMover::minSmoothField
 (
     const label nSmoothDisp,
-    const PackedBoolList& isMasterPoint,
-    const PackedBoolList& isMasterEdge,
-    const labelList& meshEdges,
+    const PackedBoolList& isPatchMasterPoint,
+    const PackedBoolList& isPatchMasterEdge,
     const scalarField& fieldMin,
     scalarField& field
 ) const
@@ -1010,13 +1016,12 @@ void Foam::medialAxisMeshMover::minSmoothField
     const edgeList& edges = pp.edges();
     const labelList& meshPoints = pp.meshPoints();
 
-    scalarField edgeWeights(meshEdges.size());
+    scalarField edgeWeights(edges.size());
     scalarField invSumWeight(meshPoints.size());
     meshRefinement::calculateEdgeWeights
     (
         mesh(),
-        isMasterEdge,
-        meshEdges,
+        isPatchMasterEdge,
         meshPoints,
         edges,
         edgeWeights,
@@ -1032,8 +1037,7 @@ void Foam::medialAxisMeshMover::minSmoothField
         meshRefinement::weightedSum
         (
             mesh(),
-            isMasterEdge,
-            meshEdges,
+            isPatchMasterEdge,
             meshPoints,
             edges,
             edgeWeights,
@@ -1064,9 +1068,7 @@ void Foam::medialAxisMeshMover::minSmoothField
         {
             scalar resid = meshRefinement::gAverage
             (
-                mesh(),
-                isMasterPoint,
-                meshPoints,
+                isPatchMasterPoint,
                 mag(field-average)()
             );
             Info<< "    Iteration " << iter << "   residual " << resid << endl;
@@ -1080,7 +1082,7 @@ void Foam::medialAxisMeshMover::
 handleFeatureAngleLayerTerminations
 (
     const scalar minCos,
-    const PackedBoolList& isMasterPoint,
+    const PackedBoolList& isPatchMasterPoint,
     const labelList& meshEdges,
     List<autoLayerDriver::extrudeMode>& extrudeStatus,
     pointField& patchDisp,
@@ -1122,7 +1124,6 @@ handleFeatureAngleLayerTerminations
 
     const labelListList& edgeFaces = pp.edgeFaces();
     const vectorField& faceNormals = pp.faceNormals();
-    const labelList& meshPoints = pp.meshPoints();
 
     forAll(edgeFaces, edgeI)
     {
@@ -1183,14 +1184,14 @@ handleFeatureAngleLayerTerminations
                     {
                         if (unmarkExtrusion(v0, patchDisp, extrudeStatus))
                         {
-                            if (isMasterPoint[meshPoints[v0]])
+                            if (isPatchMasterPoint[v0])
                             {
                                 nPointCounter++;
                             }
                         }
                         if (unmarkExtrusion(v1, patchDisp, extrudeStatus))
                         {
-                            if (isMasterPoint[meshPoints[v1]])
+                            if (isPatchMasterPoint[v1])
                             {
                                 nPointCounter++;
                             }
@@ -1213,8 +1214,8 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
 (
     const scalar minCosLayerTermination,
     const bool detectExtrusionIsland,
-    const PackedBoolList& isMasterPoint,
-    const PackedBoolList& isMasterEdge,
+    const PackedBoolList& isPatchMasterPoint,
+    const PackedBoolList& isPatchMasterEdge,
     const labelList& meshEdges,
     const scalarField& minThickness,
     List<autoLayerDriver::extrudeMode>& extrudeStatus,
@@ -1223,13 +1224,32 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
 {
     const indirectPrimitivePatch& pp = adaptPatchPtr_();
     const labelListList& pointFaces = pp.pointFaces();
-
+    const labelList& meshPoints = pp.meshPoints();
 
     Info<< typeName << " : Removing isolated regions ..." << endl;
 
     // Keep count of number of points unextruded
     label nPointCounter = 0;
 
+
+    autoPtr<OBJstream> str;
+    if (debug)
+    {
+        str.reset
+        (
+            new OBJstream
+            (
+                mesh().time().path()
+              / "islandExcludePoints_"
+              + mesh().time().timeName()
+              + ".obj"
+            )
+        );
+        Info<< typeName
+            << " : Writing points surrounded by non-extruded points to "
+            << str().name() << endl;
+    }
+
     while (true)
     {
         // Stop layer growth where mesh wraps around edge with a
@@ -1237,7 +1257,7 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
         handleFeatureAngleLayerTerminations
         (
             minCosLayerTermination,
-            isMasterPoint,
+            isPatchMasterPoint,
             meshEdges,
 
             extrudeStatus,
@@ -1353,7 +1373,7 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
         syncTools::syncPointList
         (
             mesh(),
-            pp.meshPoints(),
+            meshPoints,
             keptPoints,
             orEqOp<bool>(),
             false               // null value
@@ -1369,6 +1389,11 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
                 {
                     nPointCounter++;
                     nChanged++;
+
+                    if (str.valid())
+                    {
+                        str().write(pp.points()[meshPoints[patchPointI]]);
+                    }
                 }
             }
         }
@@ -1390,7 +1415,7 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
 
     forAll(edges, edgeI)
     {
-        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
+        if (isPatchMasterEdge[edgeI])
         {
             const edge& e = edges[edgeI];
 
@@ -1411,7 +1436,7 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
     syncTools::syncPointList
     (
         mesh(),
-        pp.meshPoints(),
+        meshPoints,
         isolatedPoint,
         plusEqOp<label>(),
         label(0)        // null value
@@ -1458,6 +1483,11 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
                     )
                     {
                         nPointCounter++;
+
+                        if (str.valid())
+                        {
+                            str().write(pp.points()[meshPoints[f[fp]]]);
+                        }
                     }
                 }
             }
@@ -1474,25 +1504,23 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
 void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
 (
     const label nSmoothDisp,
-    const PackedBoolList& isMasterPoint,
-    const PackedBoolList& isMasterEdge,
+    const PackedBoolList& isMeshMasterPoint,
+    const PackedBoolList& isMeshMasterEdge,
     vectorField& displacement
 ) const
 {
     const edgeList& edges = mesh().edges();
 
     // Correspondence between local edges/points and mesh edges/points
-    const labelList meshEdges(identity(mesh().nEdges()));
     const labelList meshPoints(identity(mesh().nPoints()));
 
     // Calculate inverse sum of weights
-    scalarField edgeWeights(meshEdges.size());
+    scalarField edgeWeights(mesh().nEdges());
     scalarField invSumWeight(meshPoints.size());
     meshRefinement::calculateEdgeWeights
     (
         mesh(),
-        isMasterEdge,
-        meshEdges,
+        isMeshMasterEdge,
         meshPoints,
         edges,
         edgeWeights,
@@ -1512,8 +1540,7 @@ void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
         meshRefinement::weightedSum
         (
             mesh(),
-            isMasterEdge,
-            meshEdges,
+            isMeshMasterEdge,
             meshPoints,
             edges,
             edgeWeights,
@@ -1533,8 +1560,7 @@ void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
         meshRefinement::weightedSum
         (
             mesh(),
-            isMasterEdge,
-            meshEdges,
+            isMeshMasterEdge,
             meshPoints,
             edges,
             edgeWeights,
@@ -1558,8 +1584,7 @@ void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
         {
             scalar resid = meshRefinement::gAverage
             (
-                mesh(),
-                isMasterPoint,
+                isMeshMasterPoint,
                 mag(displacement-average)()
             );
             Info<< "    Iteration " << iter << "   residual " << resid << endl;
@@ -1735,8 +1760,8 @@ void Foam::medialAxisMeshMover::calculateDisplacement
 
 
     // Precalulate master points/edge (only relevant for shared points/edges)
-    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh()));
-    const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh()));
+    const PackedBoolList isMeshMasterPoint(syncTools::getMasterPoints(mesh()));
+    const PackedBoolList isMeshMasterEdge(syncTools::getMasterEdges(mesh()));
     // Precalculate meshEdge per pp edge
     const labelList meshEdges
     (
@@ -1747,6 +1772,24 @@ void Foam::medialAxisMeshMover::calculateDisplacement
         )
     );
 
+    // Precalulate (patch) master point/edge
+    const PackedBoolList isPatchMasterPoint
+    (
+        meshRefinement::getMasterPoints
+        (
+            mesh(),
+            meshPoints
+        )
+    );
+    const PackedBoolList isPatchMasterEdge
+    (
+        meshRefinement::getMasterEdges
+        (
+            mesh(),
+            meshEdges
+        )
+    );
+
 
     scalarField thickness(patchDisp.size());
 
@@ -1846,7 +1889,7 @@ void Foam::medialAxisMeshMover::calculateDisplacement
 
                 patchDisp[patchPointI] = thickness[patchPointI]*n;
 
-                if (isMasterPoint[pointI])
+                if (isPatchMasterPoint[patchPointI])
                 {
                     numThicknessRatioExclude++;
                 }
@@ -1885,8 +1928,8 @@ void Foam::medialAxisMeshMover::calculateDisplacement
         minCosLayerTermination,
         detectExtrusionIsland,
 
-        isMasterPoint,
-        isMasterEdge,
+        isPatchMasterPoint,
+        isPatchMasterEdge,
         meshEdges,
         minThickness,
 
@@ -1908,9 +1951,8 @@ void Foam::medialAxisMeshMover::calculateDisplacement
     minSmoothField
     (
         nSmoothPatchThickness,
-        isMasterPoint,
-        isMasterEdge,
-        meshEdges,
+        isPatchMasterPoint,
+        isPatchMasterEdge,
         minThickness,
 
         thickness
@@ -1991,8 +2033,8 @@ void Foam::medialAxisMeshMover::calculateDisplacement
         smoothLambdaMuDisplacement
         (
             nSmoothDisplacement,
-            isMasterPoint,
-            isMasterEdge,
+            isMeshMasterPoint,
+            isMeshMasterEdge,
             displacement
         );
     }
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
index 12c2d814..c060e4d2 100644
--- a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
@@ -117,7 +117,6 @@ class medialAxisMeshMover
                 const label nSmoothDisp,
                 const PackedBoolList& isMasterPoint,
                 const PackedBoolList& isMasterEdge,
-                const labelList& meshEdges,
                 pointField& normals
             ) const;
 
@@ -175,7 +174,6 @@ class medialAxisMeshMover
                 const label nSmoothDisp,
                 const PackedBoolList& isMasterPoint,
                 const PackedBoolList& isMasterEdge,
-                const labelList& meshEdges,
                 const scalarField& fieldMin,
                 scalarField& field
             ) const;
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index e90ccce3..5d96aafc 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -1954,7 +1954,6 @@ void Foam::meshRefinement::calculateEdgeWeights
 (
     const polyMesh& mesh,
     const PackedBoolList& isMasterEdge,
-    const labelList& meshEdges,
     const labelList& meshPoints,
     const edgeList& edges,
     scalarField& edgeWeights,
@@ -1964,7 +1963,7 @@ void Foam::meshRefinement::calculateEdgeWeights
     const pointField& pts = mesh.points();
 
     // Calculate edgeWeights and inverse sum of edge weights
-    edgeWeights.setSize(meshEdges.size());
+    edgeWeights.setSize(isMasterEdge.size());
     invSumWeight.setSize(meshPoints.size());
 
     forAll(edges, edgeI)
@@ -1987,7 +1986,6 @@ void Foam::meshRefinement::calculateEdgeWeights
     (
         mesh,
         isMasterEdge,
-        meshEdges,
         meshPoints,
         edges,
         edgeWeights,
@@ -2639,8 +2637,82 @@ bool Foam::meshRefinement::write() const
 }
 
 
+Foam::PackedBoolList Foam::meshRefinement::getMasterPoints
+(
+    const polyMesh& mesh,
+    const labelList& meshPoints
+)
+{
+    const globalIndex globalPoints(meshPoints.size());
+
+    labelList myPoints(meshPoints.size());
+    forAll(meshPoints, pointI)
+    {
+        myPoints[pointI] = globalPoints.toGlobal(pointI);
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        meshPoints,
+        myPoints,
+        minEqOp<label>(),
+        labelMax
+    );
+
+
+    PackedBoolList isPatchMasterPoint(meshPoints.size());
+    forAll(meshPoints, pointI)
+    {
+        if (myPoints[pointI] == globalPoints.toGlobal(pointI))
+        {
+            isPatchMasterPoint[pointI] = true;
+        }
+    }
+
+    return isPatchMasterPoint;
+}
+
+
+Foam::PackedBoolList Foam::meshRefinement::getMasterEdges
+(
+    const polyMesh& mesh,
+    const labelList& meshEdges
+)
+{
+    const globalIndex globalEdges(meshEdges.size());
+
+    labelList myEdges(meshEdges.size());
+    forAll(meshEdges, edgeI)
+    {
+        myEdges[edgeI] = globalEdges.toGlobal(edgeI);
+    }
+
+    syncTools::syncEdgeList
+    (
+        mesh,
+        meshEdges,
+        myEdges,
+        minEqOp<label>(),
+        labelMax
+    );
+
+
+    PackedBoolList isMasterEdge(meshEdges.size());
+    forAll(meshEdges, edgeI)
+    {
+        if (myEdges[edgeI] == globalEdges.toGlobal(edgeI))
+        {
+            isMasterEdge[edgeI] = true;
+        }
+    }
+
+    return isMasterEdge;
+}
+
+
 void Foam::meshRefinement::printMeshInfo(const bool debug, const string& msg)
- const
+const
 {
     const globalMeshData& pData = mesh_.globalData();
 
@@ -2664,11 +2736,11 @@ void Foam::meshRefinement::printMeshInfo(const bool debug, const string& msg)
             }
         }
 
-        PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh_));
+        PackedBoolList isMeshMasterPoint(syncTools::getMasterPoints(mesh_));
         label nMasterPoints = 0;
-        forAll(isMasterPoint, i)
+        forAll(isMeshMasterPoint, i)
         {
-            if (isMasterPoint[i])
+            if (isMeshMasterPoint[i])
             {
                 nMasterPoints++;
             }
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 770b5025..1b1ff73b 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -733,7 +733,6 @@ public:
             (
                 const polyMesh& mesh,
                 const PackedBoolList& isMasterEdge,
-                const labelList& meshEdges,
                 const labelList& meshPoints,
                 const edgeList& edges,
                 scalarField& edgeWeights,
@@ -747,7 +746,6 @@ public:
             (
                 const polyMesh& mesh,
                 const PackedBoolList& isMasterEdge,
-                const labelList& meshEdges,
                 const labelList& meshPoints,
                 const edgeList& edges,
                 const scalarField& edgeWeights,
@@ -1069,6 +1067,22 @@ public:
                 const UList<T>& data
             );
 
+            //- Determine master point for subset of points. If coupled
+            //  chooses only one
+            static PackedBoolList getMasterPoints
+            (
+                const polyMesh& mesh,
+                const labelList& meshPoints
+            );
+
+            //- Determine master edge for subset of edges. If coupled
+            //  chooses only one
+            static PackedBoolList getMasterEdges
+            (
+                const polyMesh& mesh,
+                const labelList& meshEdges
+            );
+
             //- Print some mesh stats.
             void printMeshInfo(const bool, const string&) const;
 
@@ -1100,22 +1114,10 @@ public:
             template<class T>
             static T gAverage
             (
-                const polyMesh& mesh,
                 const PackedBoolList& isMasterElem,
                 const UList<T>& values
             );
 
-            //- Helper: calculate average over selected elements
-            template<class T>
-            static T gAverage
-            (
-                const polyMesh& mesh,
-                const PackedBoolList& isMasterElem,
-                const labelList& meshPoints,
-                const UList<T>& values
-            );
-
-
             //- Get/set write level
             static writeType writeLevel();
             static void writeLevel(const writeType);
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
index 528e8121..2cbc1f69 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,7 +57,6 @@ template<class T> void Foam::meshRefinement::updateList
 template<class T>
 T Foam::meshRefinement::gAverage
 (
-    const polyMesh& mesh,
     const PackedBoolList& isMasterElem,
     const UList<T>& values
 )
@@ -68,7 +67,6 @@ T Foam::meshRefinement::gAverage
         (
             "meshRefinement::gAverage\n"
             "(\n"
-            "    const polyMesh&,\n"
             "    const PackedBoolList& isMasterElem,\n"
             "    const UList<T>& values\n"
             ")\n"
@@ -104,58 +102,6 @@ T Foam::meshRefinement::gAverage
 }
 
 
-template<class T>
-T Foam::meshRefinement::gAverage
-(
-    const polyMesh& mesh,
-    const PackedBoolList& isMasterElem,
-    const labelList& meshElems,
-    const UList<T>& values
-)
-{
-    if (values.size() != meshElems.size())
-    {
-        FatalErrorIn
-        (
-            "meshRefinement::gAverage\n"
-            "(\n"
-            "    const polyMesh&,\n"
-            "    const labelList&,\n"
-            "    const PackedBoolList& isMasterElem,\n"
-            "    const UList<T>& values\n"
-            ")\n"
-        )   << "Number of elements in list " << values.size()
-            << " does not correspond to number of elements in meshElems "
-            << meshElems.size()
-            << exit(FatalError);
-    }
-
-    T sum = pTraits<T>::zero;
-    label n = 0;
-
-    forAll(values, i)
-    {
-        if (isMasterElem[meshElems[i]])
-        {
-            sum += values[i];
-            n++;
-        }
-    }
-
-    reduce(sum, sumOp<T>());
-    reduce(n, sumOp<label>());
-
-    if (n > 0)
-    {
-        return sum/n;
-    }
-    else
-    {
-        return pTraits<T>::max;
-    }
-}
-
-
 // Compare two lists over all boundary faces
 template<class T>
 void Foam::meshRefinement::testSyncBoundaryFaceList
@@ -343,7 +289,6 @@ void Foam::meshRefinement::weightedSum
 (
     const polyMesh& mesh,
     const PackedBoolList& isMasterEdge,
-    const labelList& meshEdges,
     const labelList& meshPoints,
     const edgeList& edges,
     const scalarField& edgeWeights,
@@ -353,15 +298,13 @@ void Foam::meshRefinement::weightedSum
 {
     if
     (
-        mesh.nEdges() != isMasterEdge.size()
-     || edges.size() != meshEdges.size()
+        edges.size() != isMasterEdge.size()
      || edges.size() != edgeWeights.size()
      || meshPoints.size() != pointData.size()
     )
     {
         FatalErrorIn("medialAxisMeshMover::weightedSum(..)")
             << "Inconsistent sizes for edge or point data:"
-            << " meshEdges:" << meshEdges.size()
             << " isMasterEdge:" << isMasterEdge.size()
             << " edgeWeights:" << edgeWeights.size()
             << " edges:" << edges.size()
@@ -375,7 +318,7 @@ void Foam::meshRefinement::weightedSum
 
     forAll(edges, edgeI)
     {
-        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
+        if (isMasterEdge[edgeI])
         {
             const edge& e = edges[edgeI];
 
diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
index 6ab8930a..ea24a1c8 100644
--- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
+++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
@@ -41,6 +41,7 @@ License
 #include "triPointRef.H"
 #include "syncTools.H"
 #include "treeDataCell.H"
+#include "DynamicField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,14 +53,15 @@ namespace Foam
     const char* Foam::NamedEnum
     <
         Foam::mappedPatchBase::sampleMode,
-        5
+        6
     >::names[] =
     {
         "nearestCell",
         "nearestPatchFace",
         "nearestPatchFaceAMI",
         "nearestPatchPoint",
-        "nearestFace"
+        "nearestFace",
+        "nearestOnlyCell"
     };
 
     template<>
@@ -76,7 +78,7 @@ namespace Foam
 }
 
 
-const Foam::NamedEnum<Foam::mappedPatchBase::sampleMode, 5>
+const Foam::NamedEnum<Foam::mappedPatchBase::sampleMode, 6>
     Foam::mappedPatchBase::sampleModeNames_;
 
 const Foam::NamedEnum<Foam::mappedPatchBase::offsetMode, 3>
@@ -187,6 +189,7 @@ void Foam::mappedPatchBase::collectSamples
 // for samples being found in two processors.
 void Foam::mappedPatchBase::findSamples
 (
+    const sampleMode mode,
     const pointField& samples,
     labelList& sampleProcs,
     labelList& sampleIndices,
@@ -199,7 +202,7 @@ void Foam::mappedPatchBase::findSamples
     // All the info for nearest. Construct to miss
     List<nearInfo> nearest(samples.size());
 
-    switch (mode_)
+    switch (mode)
     {
         case NEARESTCELL:
         {
@@ -210,7 +213,7 @@ void Foam::mappedPatchBase::findSamples
                     "mappedPatchBase::findSamples(const pointField&,"
                     " labelList&, labelList&, pointField&) const"
                 )   << "No need to supply a patch name when in "
-                    << sampleModeNames_[mode_] << " mode." << exit(FatalError);
+                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
             }
 
             //- Note: face-diagonal decomposition
@@ -244,6 +247,36 @@ void Foam::mappedPatchBase::findSamples
             break;
         }
 
+        case NEARESTONLYCELL:
+        {
+            if (samplePatch_.size() && samplePatch_ != "none")
+            {
+                FatalErrorIn
+                (
+                    "mappedPatchBase::findSamples(const pointField&,"
+                    " labelList&, labelList&, pointField&) const"
+                )   << "No need to supply a patch name when in "
+                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
+            }
+
+            //- Note: face-diagonal decomposition
+            const indexedOctree<Foam::treeDataCell>& tree = mesh.cellTree();
+
+            forAll(samples, sampleI)
+            {
+                const point& sample = samples[sampleI];
+
+                nearest[sampleI].first() = tree.findNearest(sample, sqr(GREAT));
+                nearest[sampleI].second().first() = magSqr
+                (
+                    nearest[sampleI].first().hitPoint()
+                   -sample
+                );
+                nearest[sampleI].second().second() = Pstream::myProcNo();
+            }
+            break;
+        }
+
         case NEARESTPATCHFACE:
         {
             Random rndGen(123456);
@@ -398,7 +431,7 @@ void Foam::mappedPatchBase::findSamples
                     "mappedPatchBase::findSamples(const pointField&,"
                     " labelList&, labelList&, pointField&) const"
                 )   << "No need to supply a patch name when in "
-                    << sampleModeNames_[mode_] << " mode." << exit(FatalError);
+                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
             }
 
             //- Note: face-diagonal decomposition
@@ -565,7 +598,7 @@ void Foam::mappedPatchBase::calcMapping() const
     labelList sampleProcs;
     labelList sampleIndices;
     pointField sampleLocations;
-    findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
+    findSamples(mode_, samples, sampleProcs, sampleIndices, sampleLocations);
 
     // Check for samples that were not found. This will only happen for
     // NEARESTCELL since finds cell containing a location
@@ -607,32 +640,36 @@ void Foam::mappedPatchBase::calcMapping() const
                 hasWarned = true;
             }
 
-            // Reset the samples that cannot be found to the cell centres.
-            pointField patchCc;
-            {
-                List<pointField> globalCc(Pstream::nProcs());
-                globalCc[Pstream::myProcNo()] = patch_.faceCellCentres();
-                Pstream::gatherList(globalCc);
-                Pstream::scatterList(globalCc);
-                patchCc = ListListOps::combine<pointField>
-                (
-                    globalCc,
-                    accessOp<pointField>()
-                );
-            }
+            // Collect the samples that cannot be found
+            DynamicList<label> subMap;
+            DynamicField<point> subSamples;
 
             forAll(sampleProcs, sampleI)
             {
                 if (sampleProcs[sampleI] == -1)
                 {
-                    // Reset to cell centres
-                    samples[sampleI] = patchCc[sampleI];
+                    subMap.append(sampleI);
+                    subSamples.append(samples[sampleI]);
                 }
             }
 
-            // And re-search. Note: could be optimised to only search missing
-            // points.
-            findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
+            // And re-search for pure nearest (should not fail)
+            labelList subSampleProcs;
+            labelList subSampleIndices;
+            pointField subSampleLocations;
+            findSamples
+            (
+                NEARESTONLYCELL,
+                subSamples,
+                subSampleProcs,
+                subSampleIndices,
+                subSampleLocations
+            );
+
+            // Insert
+            UIndirectList<label>(sampleProcs, subMap) = subSampleProcs;
+            UIndirectList<label>(sampleIndices, subMap) = subSampleIndices;
+            UIndirectList<point>(sampleLocations, subMap) = subSampleLocations;
         }
     }
 
diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H
index 7dd98c18..e6fcbb54 100644
--- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H
+++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,8 @@ Description
 
         // What to sample:
         // - nearestCell         : sample cell containing point
+        // - nearestOnlyCell     : nearest sample cell (even if not containing
+        //                         point)
         // - nearestPatchFace    : nearest face on selected patch
         // - nearestPatchFaceAMI : nearest face on selected patch
                                    - patches need not conform
@@ -107,11 +109,12 @@ public:
         //- Mesh items to sample
         enum sampleMode
         {
-            NEARESTCELL,         // nearest cell
+            NEARESTCELL,         // nearest cell containing sample
             NEARESTPATCHFACE,    // nearest face on selected patch
             NEARESTPATCHFACEAMI, // nearest patch face + AMI interpolation
             NEARESTPATCHPOINT,   // nearest point on selected patch
-            NEARESTFACE          // nearest face
+            NEARESTFACE,         // nearest face
+            NEARESTONLYCELL      // nearest cell (even if not containing cell)
         };
 
         //- How to project face centres
@@ -122,7 +125,7 @@ public:
             NORMAL              // use face normal + distance
         };
 
-        static const NamedEnum<sampleMode, 5> sampleModeNames_;
+        static const NamedEnum<sampleMode, 6> sampleModeNames_;
 
         static const NamedEnum<offsetMode, 3> offsetModeNames_;
 
@@ -255,6 +258,7 @@ protected:
         //- Find cells/faces containing samples
         void findSamples
         (
+            const sampleMode mode,      // search mode
             const pointField&,
             labelList& sampleProcs,     // processor containing sample
             labelList& sampleIndices,   // local index of cell/face
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
index 441aa3be..90b37e8d 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -700,4 +700,38 @@ Foam::label Foam::edgeIntersections::removeDegenerates
 }
 
 
+void Foam::edgeIntersections::replace
+(
+    const edgeIntersections& subInfo,
+    const labelList& edgeMap,
+    const labelList& faceMap
+)
+{
+    forAll(subInfo, subI)
+    {
+        const List<pointIndexHit>& subHits = subInfo[subI];
+        const labelList& subClass = subInfo.classification()[subI];
+
+        label edgeI = edgeMap[subI];
+        List<pointIndexHit>& intersections = operator[](edgeI);
+        labelList& intersectionTypes = classification_[edgeI];
+
+        intersections.setSize(subHits.size());
+        intersectionTypes.setSize(subHits.size());
+
+        forAll(subHits, i)
+        {
+            const pointIndexHit& subHit = subHits[i];
+            intersections[i] = pointIndexHit
+            (
+                subHit.hit(),
+                subHit.rawPoint(),
+                faceMap[subHit.index()]
+            );
+            intersectionTypes[i] = subClass[i];
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.H b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.H
index dc25f2e8..a2e088fd 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.H
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -197,6 +197,15 @@ public:
                 const scalarField& surf1PointTol,
                 pointField& points1
             );
+
+            //- Replace edge intersection for a subset (given as edge map and
+            //  face map - for face indices stored in pointIndexHit.index())
+            void replace
+            (
+                const edgeIntersections&,
+                const labelList& edgeMap,
+                const labelList& faceMap
+            );
 };
 
 
diff --git a/src/postProcessing/functionObjects/field/nearWallFields/nearWallFields.C b/src/postProcessing/functionObjects/field/nearWallFields/nearWallFields.C
index 1737e1de..ef4be716 100644
--- a/src/postProcessing/functionObjects/field/nearWallFields/nearWallFields.C
+++ b/src/postProcessing/functionObjects/field/nearWallFields/nearWallFields.C
@@ -52,7 +52,6 @@ void Foam::nearWallFields::calcAddressing()
     }
 
     // Global indexing
-    globalIndex globalCells(mesh.nCells());
     globalIndex globalWalls(nPatchFaces);
 
     if (debug)
diff --git a/src/sampling/sampledSet/face/faceOnlySet.C b/src/sampling/sampledSet/face/faceOnlySet.C
index 1583261d..85fb7bc7 100644
--- a/src/sampling/sampledSet/face/faceOnlySet.C
+++ b/src/sampling/sampledSet/face/faceOnlySet.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,7 @@ namespace Foam
 
 bool Foam::faceOnlySet::trackToBoundary
 (
+    passiveParticleCloud& particleCloud,
     passiveParticle& singleParticle,
     DynamicList<point>& samplingPts,
     DynamicList<label>& samplingCells,
@@ -55,7 +56,6 @@ bool Foam::faceOnlySet::trackToBoundary
     const vector smallVec = tol*offset;
     const scalar smallDist = mag(smallVec);
 
-    passiveParticleCloud particleCloud(mesh());
     particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
 
     // Alias
@@ -115,8 +115,9 @@ void Foam::faceOnlySet::calcSamples
     const vector smallVec = tol*offset;
     const scalar smallDist = mag(smallVec);
 
-    // Force calculation of minimum-tet decomposition.
-    (void) mesh().tetBasePtIs();
+    // Force calculation of cloud addressing on all processors
+    const bool oldMoving = const_cast<polyMesh&>(mesh()).moving(false);
+    passiveParticleCloud particleCloud(mesh());
 
     // Get all boundary intersections
     List<pointIndexHit> bHits = searchEngine().intersections
@@ -218,6 +219,7 @@ void Foam::faceOnlySet::calcSamples
 
         bool reachedBoundary = trackToBoundary
         (
+            particleCloud,
             singleParticle,
             samplingPts,
             samplingCells,
@@ -286,6 +288,8 @@ void Foam::faceOnlySet::calcSamples
 
         startSegmentI = samplingPts.size();
     }
+
+    const_cast<polyMesh&>(mesh()).moving(oldMoving);
 }
 
 
diff --git a/src/sampling/sampledSet/face/faceOnlySet.H b/src/sampling/sampledSet/face/faceOnlySet.H
index f50900e6..fdb47d4d 100644
--- a/src/sampling/sampledSet/face/faceOnlySet.H
+++ b/src/sampling/sampledSet/face/faceOnlySet.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,7 @@ class faceOnlySet
         //  reached
         bool trackToBoundary
         (
+            passiveParticleCloud& particleCloud,
             passiveParticle& singleParticle,
             DynamicList<point>& samplingPts,
             DynamicList<label>& samplingCells,
diff --git a/src/sampling/sampledSet/polyLine/polyLineSet.C b/src/sampling/sampledSet/polyLine/polyLineSet.C
index 5ba23c52..023acc1d 100644
--- a/src/sampling/sampledSet/polyLine/polyLineSet.C
+++ b/src/sampling/sampledSet/polyLine/polyLineSet.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,7 @@ namespace Foam
 
 bool Foam::polyLineSet::trackToBoundary
 (
+    passiveParticleCloud& particleCloud,
     passiveParticle& singleParticle,
     label& sampleI,
     DynamicList<point>& samplingPts,
@@ -51,7 +52,6 @@ bool Foam::polyLineSet::trackToBoundary
     DynamicList<scalar>& samplingCurveDist
 ) const
 {
-    passiveParticleCloud particleCloud(mesh());
     particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
 
     // Alias
@@ -157,8 +157,9 @@ void Foam::polyLineSet::calcSamples
         oldPoint = sampleCoords_[sampleI];
     }
 
-    // Force calculation of minimum-tet decomposition.
-    (void) mesh().tetBasePtIs();
+    // Force calculation of cloud addressing on all processors
+    const bool oldMoving = const_cast<polyMesh&>(mesh()).moving(false);
+    passiveParticleCloud particleCloud(mesh());
 
     // current segment number
     label segmentI = 0;
@@ -267,6 +268,7 @@ void Foam::polyLineSet::calcSamples
 
         bool bReached = trackToBoundary
         (
+            particleCloud,
             singleParticle,
             sampleI,
             samplingPts,
@@ -306,6 +308,8 @@ void Foam::polyLineSet::calcSamples
 
         startSegmentI = samplingPts.size();
     }
+
+    const_cast<polyMesh&>(mesh()).moving(oldMoving);
 }
 
 
diff --git a/src/sampling/sampledSet/polyLine/polyLineSet.H b/src/sampling/sampledSet/polyLine/polyLineSet.H
index b95e2e0c..0c1339ee 100644
--- a/src/sampling/sampledSet/polyLine/polyLineSet.H
+++ b/src/sampling/sampledSet/polyLine/polyLineSet.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,6 +65,7 @@ class polyLineSet
         //  Returns false if end of samples reached.
         bool trackToBoundary
         (
+            passiveParticleCloud& particleCloud,
             passiveParticle& singleParticle,
             label& sampleI,
             DynamicList<point>& samplingPts,
diff --git a/src/sampling/sampledSet/uniform/uniformSet.C b/src/sampling/sampledSet/uniform/uniformSet.C
index 480d373f..9689218a 100644
--- a/src/sampling/sampledSet/uniform/uniformSet.C
+++ b/src/sampling/sampledSet/uniform/uniformSet.C
@@ -77,6 +77,7 @@ bool Foam::uniformSet::nextSample
 
 bool Foam::uniformSet::trackToBoundary
 (
+    passiveParticleCloud& particleCloud,
     passiveParticle& singleParticle,
     point& samplePt,
     label& sampleI,
@@ -94,7 +95,6 @@ bool Foam::uniformSet::trackToBoundary
     // Alias
     const point& trackPt = singleParticle.position();
 
-    passiveParticleCloud particleCloud(mesh());
     particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
 
     while(true)
@@ -105,7 +105,7 @@ bool Foam::uniformSet::trackToBoundary
             // no more samples.
             if (debug)
             {
-                Info<< "trackToBoundary : Reached end : samplePt now:"
+                Pout<< "trackToBoundary : Reached end : samplePt now:"
                     << samplePt << "  sampleI now:" << sampleI << endl;
             }
             return false;
@@ -116,7 +116,7 @@ bool Foam::uniformSet::trackToBoundary
             // trackPt corresponds with samplePt. Store and use next samplePt
             if (debug)
             {
-                Info<< "trackToBoundary : samplePt corresponds to trackPt : "
+                Pout<< "trackToBoundary : samplePt corresponds to trackPt : "
                     << "  trackPt:" << trackPt << "  samplePt:" << samplePt
                     << endl;
             }
@@ -132,7 +132,7 @@ bool Foam::uniformSet::trackToBoundary
                 // no more samples.
                 if (debug)
                 {
-                    Info<< "trackToBoundary : Reached end : "
+                    Pout<< "trackToBoundary : Reached end : "
                         << "  samplePt now:" << samplePt
                         << "  sampleI now:" << sampleI
                         << endl;
@@ -145,7 +145,7 @@ bool Foam::uniformSet::trackToBoundary
 
         if (debug)
         {
-            Info<< "Searching along trajectory from "
+            Pout<< "Searching along trajectory from "
                 << "  trackPt:" << trackPt
                 << "  trackCellI:" << singleParticle.cell()
                 << "  to:" << samplePt << endl;
@@ -160,7 +160,7 @@ bool Foam::uniformSet::trackToBoundary
 
             if (debug)
             {
-                Info<< "Result of tracking "
+                Pout<< "Result of tracking "
                     << "  trackPt:" << trackPt
                     << "  trackCellI:" << singleParticle.cell()
                     << "  trackFaceI:" << singleParticle.face()
@@ -178,10 +178,10 @@ bool Foam::uniformSet::trackToBoundary
 
         if (singleParticle.onBoundary())
         {
-            //Info<< "trackToBoundary : reached boundary" << endl;
+            //Pout<< "trackToBoundary : reached boundary" << endl;
             if (mag(trackPt - samplePt) < smallDist)
             {
-                //Info<< "trackToBoundary : boundary is also sampling point"
+                //Pout<< "trackToBoundary : boundary is also sampling point"
                 //    << endl;
                 // Reached samplePt on boundary
                 samplingPts.append(trackPt);
@@ -193,7 +193,7 @@ bool Foam::uniformSet::trackToBoundary
             return true;
         }
 
-        //Info<< "trackToBoundary : reached internal sampling point" << endl;
+        //Pout<< "trackToBoundary : reached internal sampling point" << endl;
         // Reached samplePt in cell or on internal face
         samplingPts.append(trackPt);
         samplingCells.append(singleParticle.cell());
@@ -226,16 +226,14 @@ void Foam::uniformSet::calcSamples
             << exit(FatalError);
     }
 
-    // Dummy cloud to force cyclicAMI and min-tet decomposition construction
-    // even on those processors where there are no points
-    passiveParticleCloud dummyCloud(mesh());
-
-
     const vector offset = (end_ - start_)/(nPoints_ - 1);
     const vector normOffset = offset/mag(offset);
     const vector smallVec = tol*offset;
     const scalar smallDist = mag(smallVec);
 
+    // Force calculation of cloud addressing on all processors
+    const bool oldMoving = const_cast<polyMesh&>(mesh()).moving(false);
+    passiveParticleCloud particleCloud(mesh());
 
     // Get all boundary intersections
     List<pointIndexHit> bHits = searchEngine().intersections
@@ -312,6 +310,7 @@ void Foam::uniformSet::calcSamples
 
         bool reachedBoundary = trackToBoundary
         (
+            particleCloud,
             singleParticle,
             samplePt,
             sampleI,
@@ -332,7 +331,7 @@ void Foam::uniformSet::calcSamples
         {
             if (debug)
             {
-                Info<< "calcSamples : Reached end of samples: "
+                Pout<< "calcSamples : Reached end of samples: "
                     << "  samplePt now:" << samplePt
                     << "  sampleI now:" << sampleI
                     << endl;
@@ -351,7 +350,7 @@ void Foam::uniformSet::calcSamples
 
             if (debug)
             {
-                Info<< "Finding next boundary : "
+                Pout<< "Finding next boundary : "
                     << "bPoint:" << bHits[bHitI].hitPoint()
                     << "  tracking:" << singleParticle.position()
                     << "  dist:" << dist
@@ -385,6 +384,8 @@ void Foam::uniformSet::calcSamples
 
         startSegmentI = samplingPts.size();
     }
+
+    const_cast<polyMesh&>(mesh()).moving(oldMoving);
 }
 
 
@@ -422,6 +423,7 @@ void Foam::uniformSet::genSamples()
     );
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::uniformSet::uniformSet
@@ -444,7 +446,7 @@ Foam::uniformSet::uniformSet
 
     if (debug)
     {
-        write(Info);
+        write(Pout);
     }
 }
 
@@ -466,7 +468,7 @@ Foam::uniformSet::uniformSet
 
     if (debug)
     {
-        write(Info);
+        write(Pout);
     }
 }
 
diff --git a/src/sampling/sampledSet/uniform/uniformSet.H b/src/sampling/sampledSet/uniform/uniformSet.H
index 215b347f..edd5212e 100644
--- a/src/sampling/sampledSet/uniform/uniformSet.H
+++ b/src/sampling/sampledSet/uniform/uniformSet.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -82,6 +82,7 @@ class uniformSet
         //  reached
         bool trackToBoundary
         (
+            passiveParticleCloud& particleCloud,
             passiveParticle& singleParticle,
             point& samplePt,
             label& sampleI,
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C
index d1ff2be0..60ad8271 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C
@@ -55,6 +55,14 @@ void fWallFunctionFvPatchScalarField::checkType()
 }
 
 
+void fWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
+{
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+}
+
+
 scalar fWallFunctionFvPatchScalarField::yPlusLam
 (
     const scalar kappa,
@@ -233,11 +241,8 @@ void fWallFunctionFvPatchScalarField::evaluate
 
 void fWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
+    writeLocalEntries(os);
     fixedValueFvPatchField<scalar>::write(os);
-    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
-    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
-    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
-    writeEntry("value", os);
 }
 
 
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H
index b78613ee..6163d05a 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,6 +103,9 @@ protected:
         //- Check the type of the patch
         virtual void checkType();
 
+        //- Write local wall function variables
+        virtual void writeLocalEntries(Ostream&) const;
+
         //- Calculate the Y+ at the edge of the laminar sublayer
         scalar yPlusLam(const scalar kappa, const scalar E);
 
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.C
index 639dbafd..6759772f 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.C
@@ -55,6 +55,14 @@ void v2WallFunctionFvPatchScalarField::checkType()
 }
 
 
+void v2WallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
+{
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+}
+
+
 scalar v2WallFunctionFvPatchScalarField::yPlusLam
 (
     const scalar kappa,
@@ -224,11 +232,8 @@ void v2WallFunctionFvPatchScalarField::evaluate
 
 void v2WallFunctionFvPatchScalarField::write(Ostream& os) const
 {
+    writeLocalEntries(os);
     fixedValueFvPatchField<scalar>::write(os);
-    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
-    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
-    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
-    writeEntry("value", os);
 }
 
 
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.H
index bfedbbdd..dc64df93 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,6 +105,9 @@ protected:
         //- Check the type of the patch
         virtual void checkType();
 
+        //- Write local wall function variables
+        virtual void writeLocalEntries(Ostream&) const;
+
         //- Calculate the Y+ at the edge of the laminar sublayer
         scalar yPlusLam(const scalar kappa, const scalar E);
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C
index cb905c80..c24937c5 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,6 +56,14 @@ void fWallFunctionFvPatchScalarField::checkType()
 }
 
 
+void fWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
+{
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+}
+
+
 scalar fWallFunctionFvPatchScalarField::yPlusLam
 (
     const scalar kappa,
@@ -232,11 +240,8 @@ void fWallFunctionFvPatchScalarField::evaluate
 
 void fWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
+    writeLocalEntries(os);
     fixedValueFvPatchField<scalar>::write(os);
-    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
-    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
-    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
-    writeEntry("value", os);
 }
 
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H
index 771da70d..b5d92d94 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,6 +103,9 @@ protected:
         //- Check the type of the patch
         virtual void checkType();
 
+        //- Write local wall function variables
+        virtual void writeLocalEntries(Ostream&) const;
+
         //- Calculate the Y+ at the edge of the laminar sublayer
         scalar yPlusLam(const scalar kappa, const scalar E);
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.C
index b2a3e38d..19ae103d 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,6 +55,14 @@ void v2WallFunctionFvPatchScalarField::checkType()
 }
 
 
+void v2WallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
+{
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+}
+
+
 scalar v2WallFunctionFvPatchScalarField::yPlusLam
 (
     const scalar kappa,
@@ -222,11 +230,8 @@ void v2WallFunctionFvPatchScalarField::evaluate
 
 void v2WallFunctionFvPatchScalarField::write(Ostream& os) const
 {
+    writeLocalEntries(os);
     fixedValueFvPatchField<scalar>::write(os);
-    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
-    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
-    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
-    writeEntry("value", os);
 }
 
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.H
index 68bb2daa..a0a8cde8 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,6 +105,9 @@ protected:
         //- Check the type of the patch
         virtual void checkType();
 
+        //- Write local wall function variables
+        virtual void writeLocalEntries(Ostream&) const;
+
         //- Calculate the Y+ at the edge of the laminar sublayer
         scalar yPlusLam(const scalar kappa, const scalar E);