From c0a8bec3c58ccb04df577f253065940ae436a510 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 2 Feb 2010 15:49:08 +0000
Subject: [PATCH] ENH: balancing before refinement

if locally number of cells > maxLocalCells will switch from
refinement followed by balancing to weighted balancing first. This will
avoid temporary huge imbalance if all cells to be refined are on one proc.
---
 .../mesh/generation/snappyHexMesh/snappyHexMesh.C  |  16 +-
 .../generation/snappyHexMesh/snappyHexMeshDict     |  23 +-
 .../autoHexMeshDriver/autoHexMeshDriver.C          |   1 +
 .../autoHexMeshDriver/autoLayerDriver.C            | 360 ++++++++++++---------
 .../autoHexMeshDriver/autoLayerDriver.H            |   5 +-
 .../autoHexMeshDriver/autoLayerDriverShrink.C      |   9 +-
 .../autoHexMeshDriver/autoRefineDriver.C           | 119 +++++--
 .../layerParameters/layerParameters.H              |  10 +-
 .../autoHexMeshDriver/pointData/pointData.C        |   2 -
 .../autoHexMeshDriver/pointData/pointDataI.H       |  18 +-
 .../refinementParameters/refinementParameters.C    |  22 +-
 .../refinementParameters/refinementParameters.H    |  21 ++
 .../autoHexMesh/meshRefinement/meshRefinement.C    |  45 ++-
 .../autoHexMesh/meshRefinement/meshRefinement.H    |  25 +-
 .../meshRefinement/meshRefinementBaffles.C         |  31 +-
 .../meshRefinement/meshRefinementRefine.C          | 203 +++++++++++-
 .../refinementSurfaces/refinementSurfaces.C        | 252 ++++++++++-----
 .../refinementSurfaces/refinementSurfaces.H        |   6 +-
 .../autoHexMesh/trackedParticle/ExactParticle.C    |   2 +-
 .../autoHexMesh/trackedParticle/ExactParticle.H    |   2 +-
 .../searchableSurfaceCollection.C                  |  15 +-
 .../iglooWithFridges/system/snappyHexMeshDict      |  21 +-
 .../simpleFoam/motorBike/system/snappyHexMeshDict  |  21 +-
 23 files changed, 879 insertions(+), 350 deletions(-)

diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 8fc3aaf3..b2d7acff 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -192,7 +192,7 @@ int main(int argc, char *argv[])
         autoRefineDriver::debug = debug;
         autoSnapDriver::debug = debug;
         autoLayerDriver::debug = debug;
-    }        
+    }
 
 
     // Read geometry
@@ -434,6 +434,19 @@ int main(int argc, char *argv[])
         // Layer addition parameters
         layerParameters layerParams(layerDict, mesh.boundaryMesh());
 
+        //!!! Temporary hack to get access to maxLocalCells
+        bool preBalance;
+        {
+            refinementParameters refineParams(refineDict);
+
+            preBalance = returnReduce
+            (
+                (mesh.nCells() >= refineParams.maxLocalCells()),
+                orOp<bool>()
+            );
+        }
+
+
         if (!overwrite)
         {
             const_cast<Time&>(mesh.time())++;
@@ -444,6 +457,7 @@ int main(int argc, char *argv[])
             layerDict,
             motionDict,
             layerParams,
+            preBalance,
             decomposer,
             distributor
         );
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 29c05a5f..9d2367e5 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    object      autoHexMeshDict;
+    object      snappyHexMeshDict;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -71,10 +71,10 @@ castellatedMeshControls
     // Refinement parameters
     // ~~~~~~~~~~~~~~~~~~~~~
 
-    // While refining maximum number of cells per processor. This is basically
-    // the number of cells that fit on a processor. If you choose this too small
-    // it will do just more refinement iterations to obtain a similar mesh.
-    maxLocalCells 1000000;
+    // If local number of cells is >= maxLocalCells on any processor
+    // switches from from refinement followed by balancing
+    // (current method) to (weighted) balancing before refinement.
+    maxLocalCells 100000;
 
     // Overall cell limit (approximately). Refinement will stop immediately
     // upon reaching this number so a refinement level might not complete.
@@ -89,6 +89,13 @@ castellatedMeshControls
     // (unless the number of cells to refine is 0)
     minRefinementCells 0;
 
+    // Allow a certain level of imbalance during refining
+    // (since balancing is quite expensive)
+    // Expressed as fraction of perfect balance (= overall number of cells /
+    // nProcs). 0=balance always.
+    maxLoadUnbalance 0.10;
+
+
     // Number of buffer layers between different levels.
     // 1 means normal 2:1 refinement restriction, larger means slower
     // refinement.
@@ -184,6 +191,12 @@ castellatedMeshControls
     // NOTE: This point should never be on a face, always inside a cell, even
     // after refinement.
     locationInMesh (5 0.28 0.43);
+
+
+    // Whether any faceZones (as specified in the refinementSurfaces)
+    // are only on the boundary of corresponding cellZones or also allow
+    // free-standing zone faces
+    allowFreeStandingZoneFaces true;
 }
 
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
index cc110c1d..76e589c0 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
@@ -539,6 +539,7 @@ void Foam::autoHexMeshDriver::doMesh()
             shrinkDict,
             motionDict,
             layerParams,
+            true,                   // pre-balance
             decomposerPtr_(),
             distributorPtr_()
         );
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 425adee9..30796066 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -2515,138 +2515,159 @@ void Foam::autoLayerDriver::addLayers
 (
     const layerParameters& layerParams,
     const dictionary& motionDict,
+    const labelList& patchIDs,
     const label nAllowableErrors,
-    motionSmoother& meshMover,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor
 )
 {
     fvMesh& mesh = meshRefiner_.mesh();
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const labelList& meshPoints = pp.meshPoints();
 
-    // Precalculate mesh edge labels for patch edges
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    autoPtr<indirectPrimitivePatch> pp
+    (
+        meshRefinement::makePatch
+        (
+            mesh,
+            patchIDs
+        )
+    );
 
-    labelList meshEdges(pp.nEdges());
-    forAll(pp.edges(), edgeI)
-    {
-        const edge& ppEdge = pp.edges()[edgeI];
-        label v0 = meshPoints[ppEdge[0]];
-        label v1 = meshPoints[ppEdge[1]];
-        meshEdges[edgeI] = meshTools::findEdge
+    // Construct iterative mesh mover.
+    Info<< "Constructing mesh displacer ..." << endl;
+
+    autoPtr<motionSmoother> meshMover
+    (
+        new motionSmoother
         (
-            mesh.edges(),
-            mesh.pointEdges()[v0],
-            v0,
-            v1
-        );
-    }
+            mesh,
+            pp(),
+            patchIDs,
+            meshRefinement::makeDisplacementField
+            (
+                pointMesh::New(mesh),
+                patchIDs
+            ),
+            motionDict
+        )
+    );
+
+
+    // Point-wise extrusion data
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Displacement for all pp.localPoints.
-    vectorField patchDisp(pp.nPoints(), vector::one);
+    vectorField patchDisp(pp().nPoints(), vector::one);
 
     // Number of layers for all pp.localPoints. Note: only valid if
     // extrudeStatus = EXTRUDE.
-    labelList patchNLayers(pp.nPoints(), 0);
+    labelList patchNLayers(pp().nPoints(), 0);
 
     // Whether to add edge for all pp.localPoints.
-    List<extrudeMode> extrudeStatus(pp.nPoints(), EXTRUDE);
-
-
-    // Get number of layer per point from number of layers per patch
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
 
-    setNumLayers
-    (
-        layerParams.numLayers(),    // per patch the num layers
-        meshMover.adaptPatchIDs(),  // patches that are being moved
-        pp,                         // indirectpatch for all faces moving
-
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
-
-    // Disable extrusion on non-manifold points
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    {
+        // Get number of layer per point from number of layers per patch
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    handleNonManifolds
-    (
-        pp,
-        meshEdges,
+        setNumLayers
+        (
+            layerParams.numLayers(),    // per patch the num layers
+            meshMover().adaptPatchIDs(),// patches that are being moved
+            pp,                         // indirectpatch for all faces moving
 
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
+            patchDisp,
+            patchNLayers,
+            extrudeStatus
+        );
 
-    // Disable extrusion on feature angles
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // Precalculate mesh edge labels for patch edges
+        labelList meshEdges(pp().meshEdges(mesh.edges(), mesh.pointEdges()));
 
-    handleFeatureAngle
-    (
-        pp,
-        meshEdges,
-        layerParams.featureAngle()*mathematicalConstant::pi/180.0,
+        // Disable extrusion on non-manifold points
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
+        handleNonManifolds
+        (
+            pp,
+            meshEdges,
 
-    // Disable extrusion on warped faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            patchDisp,
+            patchNLayers,
+            extrudeStatus
+        );
 
-    // Undistorted edge length
-    const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
-    const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
+        // Disable extrusion on feature angles
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    handleWarpedFaces
-    (
-        pp,
-        layerParams.maxFaceThicknessRatio(),
-        edge0Len,
-        cellLevel,
+        handleFeatureAngle
+        (
+            pp,
+            meshEdges,
+            layerParams.featureAngle()*mathematicalConstant::pi/180.0,
 
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
+            patchDisp,
+            patchNLayers,
+            extrudeStatus
+        );
 
-    //// Disable extrusion on cells with multiple patch faces
-    //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    //
-    //handleMultiplePatchFaces
-    //(
-    //    pp,
-    //
-    //    patchDisp,
-    //    patchNLayers,
-    //    extrudeStatus
-    //);
+        // Disable extrusion on warped faces
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+        // Undistorted edge length
+        const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
+        const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
 
-    // Grow out region of non-extrusion
-    for (label i = 0; i < layerParams.nGrow(); i++)
-    {
-        growNoExtrusion
+        handleWarpedFaces
         (
             pp,
+            layerParams.maxFaceThicknessRatio(),
+            edge0Len,
+            cellLevel,
+
             patchDisp,
             patchNLayers,
             extrudeStatus
         );
+
+        //// Disable extrusion on cells with multiple patch faces
+        //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        //
+        //handleMultiplePatchFaces
+        //(
+        //    pp,
+        //
+        //    patchDisp,
+        //    patchNLayers,
+        //    extrudeStatus
+        //);
+
+
+        // Grow out region of non-extrusion
+        for (label i = 0; i < layerParams.nGrow(); i++)
+        {
+            growNoExtrusion
+            (
+                pp,
+                patchDisp,
+                patchNLayers,
+                extrudeStatus
+            );
+        }
     }
 
+
+    // Undistorted edge length
+    const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
+    const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
+
     // Determine (wanted) point-wise layer thickness and expansion ratio
-    scalarField thickness(pp.nPoints());
-    scalarField minThickness(pp.nPoints());
-    scalarField expansionRatio(pp.nPoints());
+    scalarField thickness(pp().nPoints());
+    scalarField minThickness(pp().nPoints());
+    scalarField expansionRatio(pp().nPoints());
     calculateLayerThickness
     (
         pp,
-        meshMover.adaptPatchIDs(),
+        meshMover().adaptPatchIDs(),
         layerParams.expansionRatio(),
 
         layerParams.relativeSizes(),        // thickness relative to cellsize?
@@ -2669,14 +2690,14 @@ void Foam::autoLayerDriver::addLayers
 
         // Find maximum length of a patch name, for a nicer output
         label maxPatchNameLen = 0;
-        forAll(meshMover.adaptPatchIDs(), i)
+        forAll(meshMover().adaptPatchIDs(), i)
         {
-            label patchI = meshMover.adaptPatchIDs()[i];
+            label patchI = meshMover().adaptPatchIDs()[i];
             word patchName = patches[patchI].name();
-            maxPatchNameLen = max(maxPatchNameLen,label(patchName.size()));
+            maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
         }
 
-        Info<< nl 
+        Info<< nl
             << setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
             << setw(0) << " faces    layers avg thickness[m]" << nl
             << setf(ios_base::left) << setw(maxPatchNameLen) << " "
@@ -2684,9 +2705,9 @@ void Foam::autoLayerDriver::addLayers
             << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
             << setw(0) << " -----    ------ --------- -------" << endl;
 
-        forAll(meshMover.adaptPatchIDs(), i)
+        forAll(meshMover().adaptPatchIDs(), i)
         {
-            label patchI = meshMover.adaptPatchIDs()[i];
+            label patchI = meshMover().adaptPatchIDs()[i];
 
             const labelList& meshPoints = patches[patchI].meshPoints();
 
@@ -2697,7 +2718,7 @@ void Foam::autoLayerDriver::addLayers
 
             forAll(meshPoints, patchPointI)
             {
-                label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
+                label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]];
 
                 //maxThickness = max(maxThickness, thickness[ppPointI]);
                 //minThickness = min(minThickness, thickness[ppPointI]);
@@ -2767,7 +2788,7 @@ void Foam::autoLayerDriver::addLayers
             IOobject::NO_WRITE,
             false
         ),
-        meshMover.pMesh(),
+        meshMover().pMesh(),
         dimensionedScalar("pointMedialDist", dimless, 0.0)
     );
 
@@ -2782,7 +2803,7 @@ void Foam::autoLayerDriver::addLayers
             IOobject::NO_WRITE,
             false
         ),
-        meshMover.pMesh(),
+        meshMover().pMesh(),
         dimensionedVector("dispVec", dimless, vector::zero)
     );
 
@@ -2797,7 +2818,7 @@ void Foam::autoLayerDriver::addLayers
             IOobject::NO_WRITE,
             false
         ),
-        meshMover.pMesh(),
+        meshMover().pMesh(),
         dimensionedScalar("medialRatio", dimless, 0.0)
     );
 
@@ -2877,7 +2898,7 @@ void Foam::autoLayerDriver::addLayers
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
         {
-            pointField oldPatchPos(pp.localPoints());
+            pointField oldPatchPos(pp().localPoints());
 
             //// Laplacian displacement shrinking.
             //shrinkMeshDistance
@@ -2892,7 +2913,7 @@ void Foam::autoLayerDriver::addLayers
             // Medial axis based shrinking
             shrinkMeshMedialDistance
             (
-                meshMover,
+                meshMover(),
                 meshQualityDict,
 
                 layerParams.nSmoothThickness(),
@@ -2914,7 +2935,7 @@ void Foam::autoLayerDriver::addLayers
             );
 
             // Update patchDisp (since not all might have been honoured)
-            patchDisp = oldPatchPos - pp.localPoints();
+            patchDisp = oldPatchPos - pp().localPoints();
         }
 
         // Truncate displacements that are too small (this will do internal
@@ -2922,7 +2943,7 @@ void Foam::autoLayerDriver::addLayers
         faceSet dummySet(mesh, "wrongPatchFaces", 0);
         truncateDisplacement
         (
-            meshMover,
+            meshMover(),
             minThickness,
             dummySet,
             patchDisp,
@@ -2937,7 +2958,7 @@ void Foam::autoLayerDriver::addLayers
             dumpDisplacement
             (
                 mesh.time().path()/"layer",
-                pp,
+                pp(),
                 patchDisp,
                 extrudeStatus
             );
@@ -2961,11 +2982,11 @@ void Foam::autoLayerDriver::addLayers
         // Determine per point/per face number of layers to extrude. Also
         // handles the slow termination of layers when going switching layers
 
-        labelList nPatchPointLayers(pp.nPoints(),-1);
-        labelList nPatchFaceLayers(pp.localFaces().size(),-1);
+        labelList nPatchPointLayers(pp().nPoints(),-1);
+        labelList nPatchFaceLayers(pp().localFaces().size(),-1);
         setupLayerInfoTruncation
         (
-            meshMover,
+            meshMover(),
             patchNLayers,
             extrudeStatus,
             layerParams.nBufferCellsNoExtrude(),
@@ -3004,7 +3025,7 @@ void Foam::autoLayerDriver::addLayers
         addLayer.setRefinement
         (
             invExpansionRatio,
-            pp,
+            pp(),
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point
             firstDisp,          // thickness of layer nearest internal mesh
@@ -3055,8 +3076,8 @@ void Foam::autoLayerDriver::addLayers
         addLayer.updateMesh
         (
             map,
-            identity(pp.size()),
-            identity(pp.nPoints())
+            identity(pp().size()),
+            identity(pp().nPoints())
         );
 
         // Collect layer faces and cells for outside loop.
@@ -3103,7 +3124,7 @@ void Foam::autoLayerDriver::addLayers
         (
             addLayer,
             meshQualityDict,
-            pp,
+            pp(),
             newMesh,
 
             patchDisp,
@@ -3112,7 +3133,7 @@ void Foam::autoLayerDriver::addLayers
         );
 
         Info<< "Extruding " << countExtrusion(pp, extrudeStatus)
-            << " out of " << returnReduce(pp.size(), sumOp<label>())
+            << " out of " << returnReduce(pp().size(), sumOp<label>())
             << " faces. Removed extrusion at " << nTotChanged << " faces."
             << endl;
 
@@ -3122,8 +3143,8 @@ void Foam::autoLayerDriver::addLayers
         }
 
         // Reset mesh points and start again
-        meshMover.movePoints(oldPoints);
-        meshMover.correct();
+        meshMover().movePoints(oldPoints);
+        meshMover().correct();
 
         Info<< endl;
     }
@@ -3178,6 +3199,7 @@ void Foam::autoLayerDriver::addLayers
         (
             false,
             false,
+            scalarField(mesh.nCells(), 1.0),
             decomposer,
             distributor
         );
@@ -3212,11 +3234,12 @@ void Foam::autoLayerDriver::doLayers
     const dictionary& shrinkDict,
     const dictionary& motionDict,
     const layerParameters& layerParams,
+    const bool preBalance,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor
 )
 {
-    fvMesh& mesh = meshRefiner_.mesh();
+    const fvMesh& mesh = meshRefiner_.mesh();
 
     Info<< nl
         << "Shrinking and layer addition phase" << nl
@@ -3250,59 +3273,80 @@ void Foam::autoLayerDriver::doLayers
     }
     else
     {
-        autoPtr<indirectPrimitivePatch> ppPtr
+        // Check that outside of mesh is not multiply connected.
+        checkMeshManifold();
+
+        // Check initial mesh
+        Info<< "Checking initial mesh ..." << endl;
+        labelHashSet wrongFaces(mesh.nFaces()/100);
+        motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
+        const label nInitErrors = returnReduce
         (
-            meshRefinement::makePatch
-            (
-                mesh,
-                patchIDs
-            )
+            wrongFaces.size(),
+            sumOp<label>()
         );
-        indirectPrimitivePatch& pp = ppPtr();
-
-        // Construct iterative mesh mover.
-        Info<< "Constructing mesh displacer ..." << endl;
-
-        {
-            const pointMesh& pMesh = pointMesh::New(mesh);
 
-            motionSmoother meshMover
-            (
-                mesh,
-                pp,
-                patchIDs,
-                meshRefinement::makeDisplacementField(pMesh, patchIDs),
-                motionDict
-            );
-
-            // Check that outside of mesh is not multiply connected.
-            checkMeshManifold();
+        Info<< "Detected " << nInitErrors << " illegal faces"
+            << " (concave, zero area or negative cell pyramid volume)"
+            << endl;
 
-            // Check initial mesh
-            Info<< "Checking initial mesh ..." << endl;
-            labelHashSet wrongFaces(mesh.nFaces()/100);
-            motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
-            const label nInitErrors = returnReduce
-            (
-                wrongFaces.size(),
-                sumOp<label>()
-            );
 
-            Info<< "Detected " << nInitErrors << " illegal faces"
-                << " (concave, zero area or negative cell pyramid volume)"
+        // Balance
+        if (Pstream::parRun() && preBalance)
+        {
+            Info<< nl
+                << "Doing initial balancing" << nl
+                << "-----------------------" << nl
                 << endl;
 
-            // Do all topo changes
-            addLayers
+            scalarField cellWeights(mesh.nCells(), 1);
+            forAll(numLayers, patchI)
+            {
+                if (numLayers[patchI] > 0)
+                {
+                    const polyPatch& pp = mesh.boundaryMesh()[patchI];
+                    forAll(pp.faceCells(), i)
+                    {
+                        cellWeights[pp.faceCells()[i]] += numLayers[patchI];
+                    }
+                }
+            }
+
+            // Balance mesh (and meshRefinement). No restriction on face zones
+            // and baffles.
+            autoPtr<mapDistributePolyMesh> map = meshRefiner_.balance
             (
-                layerParams,
-                motionDict,
-                nInitErrors,
-                meshMover,
+                false,
+                false,
+                cellWeights,
                 decomposer,
                 distributor
             );
+
+            //{
+            //    globalIndex globalCells(mesh.nCells());
+            //
+            //    Info<< "** Distribution after balancing:" << endl;
+            //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+            //    {
+            //        Info<< "    " << procI << '\t'
+            //            << globalCells.localSize(procI) << endl;
+            //    }
+            //    Info<< endl;
+            //}
         }
+
+
+        // Do all topo changes
+        addLayers
+        (
+            layerParams,
+            motionDict,
+            patchIDs,
+            nInitErrors,
+            decomposer,
+            distributor
+        );
     }
 }
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index 12569570..748f26d8 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -396,7 +396,7 @@ class autoLayerDriver
                     const PackedBoolList& isMasterEdge,
                     const labelList& meshEdges,
                     const scalarField& fieldMin,
-                    const label& nSmoothDisp,
+                    const label nSmoothDisp,
                     scalarField& field
                 ) const;
 
@@ -523,8 +523,8 @@ public:
             (
                 const layerParameters& layerParams,
                 const dictionary& motionDict,
+                const labelList& patchIDs,
                 const label nAllowableErrors,
-                motionSmoother& meshMover,
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor
             );
@@ -535,6 +535,7 @@ public:
                 const dictionary& shrinkDict,
                 const dictionary& motionDict,
                 const layerParameters& layerParams,
+                const bool preBalance,              // balance before adding?
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor
             );
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
index cc88231b..c572f2fd 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
@@ -35,9 +35,6 @@ Description
 #include "pointData.H"
 #include "PointEdgeWave.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Calculate inverse sum of edge weights (currently always 1.0)
@@ -94,7 +91,7 @@ void Foam::autoLayerDriver::smoothField
     const PackedBoolList& isMasterEdge,
     const labelList& meshEdges,
     const scalarField& fieldMin,
-    const label& nSmoothDisp,
+    const label nSmoothDisp,
     scalarField& field
 ) const
 {
@@ -240,7 +237,7 @@ void Foam::autoLayerDriver::smoothNormals
     const edgeList& edges = mesh.edges();
 
     // Points that do not change.
-    PackedBoolList isFixedPoint(mesh.nPoints(), 0);
+    PackedBoolList isFixedPoint(mesh.nPoints());
 
     // Internal points that are fixed
     forAll(fixedPoints, i)
@@ -468,7 +465,7 @@ void Foam::autoLayerDriver::findIsolatedRegions
     // Keep count of number of points unextruded
     label nPointCounter = 0;
 
-    while(true)
+    while (true)
     {
         // Stop layer growth where mesh wraps around edge with a
         // large feature angle
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index a9b31b9b..ccc63d26 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -28,14 +28,14 @@ License
 #include "meshRefinement.H"
 #include "fvMesh.H"
 #include "Time.H"
-#include "boundBox.H"
-#include "mapDistributePolyMesh.H"
 #include "cellSet.H"
 #include "syncTools.H"
 #include "refinementParameters.H"
 #include "featureEdgeMesh.H"
 #include "refinementSurfaces.H"
 #include "shellSurfaces.H"
+#include "mapDistributePolyMesh.H"
+#include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -205,13 +205,36 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
                 const_cast<Time&>(mesh.time())++;
             }
 
-            meshRefiner_.refineAndBalance
+
+            if
             (
-                "feature refinement iteration " + name(iter),
-                decomposer_,
-                distributor_,
-                cellsToRefine
-            );
+                returnReduce
+                (
+                    (mesh.nCells() >= refineParams.maxLocalCells()),
+                    orOp<bool>()
+                )
+            )
+            {
+                meshRefiner_.balanceAndRefine
+                (
+                    "feature refinement iteration " + name(iter),
+                    decomposer_,
+                    distributor_,
+                    cellsToRefine,
+                    refineParams.maxLoadUnbalance()
+                );
+            }
+            else
+            {
+                meshRefiner_.refineAndBalance
+                (
+                    "feature refinement iteration " + name(iter),
+                    decomposer_,
+                    distributor_,
+                    cellsToRefine,
+                    refineParams.maxLoadUnbalance()
+                );
+            }
         }
     }
     return iter;
@@ -305,13 +328,36 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
             const_cast<Time&>(mesh.time())++;
         }
 
-        meshRefiner_.refineAndBalance
+
+        if
         (
-            "surface refinement iteration " + name(iter),
-            decomposer_,
-            distributor_,
-            cellsToRefine
-        );
+            returnReduce
+            (
+                (mesh.nCells() >= refineParams.maxLocalCells()),
+                orOp<bool>()
+            )
+        )
+        {
+            meshRefiner_.balanceAndRefine
+            (
+                "surface refinement iteration " + name(iter),
+                decomposer_,
+                distributor_,
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
+            );
+        }
+        else
+        {
+            meshRefiner_.refineAndBalance
+            (
+                "surface refinement iteration " + name(iter),
+                decomposer_,
+                distributor_,
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
+            );
+        }
     }
     return iter;
 }
@@ -487,13 +533,35 @@ Foam::label Foam::autoRefineDriver::shellRefine
             const_cast<Time&>(mesh.time())++;
         }
 
-        meshRefiner_.refineAndBalance
+        if
         (
-            "shell refinement iteration " + name(iter),
-            decomposer_,
-            distributor_,
-            cellsToRefine
-        );
+            returnReduce
+            (
+                (mesh.nCells() >= refineParams.maxLocalCells()),
+                orOp<bool>()
+            )
+        )
+        {
+            meshRefiner_.balanceAndRefine
+            (
+                "shell refinement iteration " + name(iter),
+                decomposer_,
+                distributor_,
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
+            );
+        }
+        else
+        {
+            meshRefiner_.refineAndBalance
+            (
+                "shell refinement iteration " + name(iter),
+                decomposer_,
+                distributor_,
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
+            );
+        }
     }
     meshRefiner_.userFaceData().clear();
 
@@ -557,7 +625,11 @@ void Foam::autoRefineDriver::zonify
             const_cast<Time&>(mesh.time())++;
         }
 
-        meshRefiner_.zonify(refineParams.keepPoints()[0]);
+        meshRefiner_.zonify
+        (
+            refineParams.keepPoints()[0],
+            refineParams.allowFreeStandingZoneFaces()
+        );
 
         if (debug)
         {
@@ -776,11 +848,14 @@ void Foam::autoRefineDriver::doRefine
             const_cast<Time&>(mesh.time())++;
         }
 
-        // Do final balancing. Keep zoned faces on one processor.
+        // Do final balancing. Keep zoned faces on one processor since the
+        // snap phase will convert them to baffles and this only works for
+        // internal faces.
         meshRefiner_.balance
         (
             true,
             false,
+            scalarField(mesh.nCells(), 1), // dummy weights
             decomposer_,
             distributor_
         );
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
index 766897f4..46ef4c6f 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
@@ -206,20 +206,20 @@ public:
                 return nGrow_;
             }
 
-            //- Number of smoothing iterations of surface normals 
+            //- Number of smoothing iterations of surface normals
             label nSmoothSurfaceNormals() const
             {
                 return nSmoothSurfaceNormals_;
             }
 
             //- Number of smoothing iterations of interior mesh movement
-            //  direction  
+            //  direction
             label nSmoothNormals() const
             {
                 return nSmoothNormals_;
             }
 
-            //- Stop layer growth on highly warped cells 
+            //- Stop layer growth on highly warped cells
             scalar maxFaceThicknessRatio() const
             {
                 return maxFaceThicknessRatio_;
@@ -236,8 +236,8 @@ public:
                 return nSmoothThickness_;
             }
 
-            //- Reduce layer growth where ratio thickness to medial 
-            //  distance is large 
+            //- Reduce layer growth where ratio thickness to medial
+            //  distance is large
             scalar maxThicknessToMedialRatio() const
             {
                 return maxThicknessToMedialRatio_;
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.C
index 9b4d224a..30b1feef 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.C
@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "pointData.H"
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H
index 9832bbf8..9366afbe 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H
@@ -48,7 +48,7 @@ inline bool Foam::pointData::update
         v_ = w2.v();
 
         return true;
-    }       
+    }
 
 
 //    if (v_ != w2.v())
@@ -82,7 +82,7 @@ inline bool Foam::pointData::update
     }
 }
 
-    
+
 // Update this with w2 (information on same point)
 inline bool Foam::pointData::update
 (
@@ -98,7 +98,7 @@ inline bool Foam::pointData::update
         s_ = w2.s();
         v_ = w2.v();
         return true;
-    }        
+    }
 
 
 //    if (v_ != w2.v())
@@ -148,7 +148,7 @@ inline Foam::pointData::pointData()
 // Construct from origin, distance
 inline Foam::pointData::pointData
 (
-    const point& origin, 
+    const point& origin,
     const scalar distSqr,
     const scalar s,
     const vector& v
@@ -271,14 +271,14 @@ inline bool Foam::pointData::updatePoint
     const scalar tol
 )
 {
-    return 
+    return
         update
         (
             mesh.points()[pointI],
             edgeInfo,
             tol
         );
-    }    
+    }
 
 
 // Update this with new information on same point
@@ -297,7 +297,7 @@ inline bool Foam::pointData::updatePoint
             newPointInfo,
             tol
         );
-}    
+}
 
 
 // Update this with new information on same point. No extra information.
@@ -308,7 +308,7 @@ inline bool Foam::pointData::updatePoint
 )
 {
     return update(newPointInfo, tol);
-}    
+}
 
 
 // Update this with information from connected point
@@ -334,7 +334,7 @@ inline bool Foam::pointData::updateEdge
             pointInfo,
             tol
         );
-}    
+}
 
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
index d63d527f..48ed0281 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
@@ -43,7 +43,16 @@ Foam::refinementParameters::refinementParameters
     minRefineCells_(readLabel(dict.lookup("minimumRefine"))),
     curvature_(readScalar(dict.lookup("curvature"))),
     nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))),
-    keepPoints_(dict.lookup("keepPoints"))
+    keepPoints_(dict.lookup("keepPoints")),
+    allowFreeStandingZoneFaces_
+    (
+        dict.lookupOrDefault<Switch>
+        (
+            "allowFreeStandingZoneFaces",
+            true
+        )
+    ),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
 {}
 
 
@@ -53,7 +62,16 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
     maxLocalCells_(readLabel(dict.lookup("maxLocalCells"))),
     minRefineCells_(readLabel(dict.lookup("minRefinementCells"))),
     nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))),
-    keepPoints_(pointField(1, dict.lookup("locationInMesh")))
+    keepPoints_(pointField(1, dict.lookup("locationInMesh"))),
+    allowFreeStandingZoneFaces_
+    (
+        dict.lookupOrDefault<Switch>
+        (
+            "allowFreeStandingZoneFaces",
+            true
+        )
+    ),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
 {
     scalar featAngle(readScalar(dict.lookup("resolveFeatureAngle")));
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
index 1e12dd30..32d67772 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
@@ -38,6 +38,7 @@ SourceFiles
 
 #include "dictionary.H"
 #include "pointField.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,6 +74,13 @@ class refinementParameters
         //- Areas to keep
         const pointField keepPoints_;
 
+        //- FaceZone faces allowed which have owner and neighbour in same
+        //  cellZone?
+        Switch allowFreeStandingZoneFaces_;
+
+        //- Allowed load unbalance
+        scalar maxLoadUnbalance_;
+
 
     // Private Member Functions
 
@@ -134,6 +142,19 @@ public:
                 return keepPoints_;
             }
 
+            //- Are zone faces allowed only inbetween different cell zones
+            //  or also just free standing?
+            bool allowFreeStandingZoneFaces() const
+            {
+                return allowFreeStandingZoneFaces_;
+            }
+
+            //- Allowed load unbalance
+            scalar maxLoadUnbalance() const
+            {
+                return maxLoadUnbalance_;
+            }
+
 
         // Other
 
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 5f7826cd..c4ea71f4 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -54,6 +54,7 @@ License
 #include "Random.H"
 #include "searchableSurfaces.H"
 #include "treeBoundBox.H"
+#include "zeroGradientFvPatchFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -552,13 +553,16 @@ void Foam::meshRefinement::calcLocalRegions
     const globalIndex& globalCells,
     const labelList& globalRegion,
     const Map<label>& coupledRegionToMaster,
+    const scalarField& cellWeights,
 
     Map<label>& globalToLocalRegion,
-    pointField& localPoints
+    pointField& localPoints,
+    scalarField& localWeights
 ) const
 {
     globalToLocalRegion.resize(globalRegion.size());
     DynamicList<point> localCc(globalRegion.size()/2);
+    DynamicList<scalar> localWts(globalRegion.size()/2);
 
     forAll(globalRegion, cellI)
     {
@@ -573,6 +577,7 @@ void Foam::meshRefinement::calcLocalRegions
                 // I am master. Allocate region for me.
                 globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
                 localCc.append(mesh_.cellCentres()[cellI]);
+                localWts.append(cellWeights[cellI]);
             }
         }
         else
@@ -581,11 +586,13 @@ void Foam::meshRefinement::calcLocalRegions
             if (globalToLocalRegion.insert(globalRegion[cellI], localCc.size()))
             {
                 localCc.append(mesh_.cellCentres()[cellI]);
+                localWts.append(cellWeights[cellI]);
             }
         }
     }
 
     localPoints.transfer(localCc);
+    localWeights.transfer(localWts);
 
     if (localPoints.size() != globalToLocalRegion.size())
     {
@@ -924,6 +931,7 @@ Foam::label Foam::meshRefinement::countHits() const
 // Determine distribution to move connected regions onto one processor.
 Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 (
+    const scalarField& cellWeights,
     const boolList& blockedFace,
     const List<labelPair>& explicitConnections,
     decompositionMethod& decomposer
@@ -965,14 +973,17 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 
     Map<label> globalToLocalRegion;
     pointField localPoints;
+    scalarField localWeights;
     calcLocalRegions
     (
         globalCells,
         globalRegion,
         coupledRegionToMaster,
+        cellWeights,
 
         globalToLocalRegion,
-        localPoints
+        localPoints,
+        localWeights
     );
 
 
@@ -984,7 +995,7 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 
     if (isA<geomDecomp>(decomposer))
     {
-        regionDistribution = decomposer.decompose(localPoints);
+        regionDistribution = decomposer.decompose(localPoints, localWeights);
     }
     else
     {
@@ -998,7 +1009,12 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
             regionRegions
         );
 
-        regionDistribution = decomposer.decompose(regionRegions, localPoints);
+        regionDistribution = decomposer.decompose
+        (
+            regionRegions,
+            localPoints,
+            localWeights
+        );
     }
 
 
@@ -1058,6 +1074,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 (
     const bool keepZoneFaces,
     const bool keepBaffles,
+    const scalarField& cellWeights,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor
 )
@@ -1163,6 +1180,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 
                 distribution = decomposeCombineRegions
                 (
+                    cellWeights,
                     blockedFace,
                     couples,
                     decomposer
@@ -1182,13 +1200,21 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
             else
             {
                 // Normal decomposition
-                distribution = decomposer.decompose(mesh_.cellCentres());
+                distribution = decomposer.decompose
+                (
+                    mesh_.cellCentres(),
+                    cellWeights
+                );
             }
         }
         else
         {
             // Normal decomposition
-            distribution = decomposer.decompose(mesh_.cellCentres());
+            distribution = decomposer.decompose
+            (
+                mesh_.cellCentres(),
+                cellWeights
+            );
         }
 
         if (debug)
@@ -1249,7 +1275,7 @@ Foam::labelList Foam::meshRefinement::intersectedPoints() const
     const faceList& faces = mesh_.faces();
 
     // Mark all points on faces that will become baffles
-    PackedBoolList isBoundaryPoint(mesh_.nPoints(), 0u);
+    PackedBoolList isBoundaryPoint(mesh_.nPoints());
     label nBoundaryPoints = 0;
 
     forAll(surfaceIndex_, faceI)
@@ -1683,7 +1709,7 @@ Foam::label Foam::meshRefinement::addPatch
 
 Foam::label Foam::meshRefinement::addMeshedPatch
 (
-    const word& name,   
+    const word& name,
     const word& type
 )
 {
@@ -2103,7 +2129,8 @@ void Foam::meshRefinement::dumpRefinementLevel() const
             false
         ),
         mesh_,
-        dimensionedScalar("zero", dimless, 0)
+        dimensionedScalar("zero", dimless, 0),
+        zeroGradientFvPatchScalarField::typeName
     );
 
     const labelList& cellLevel = meshCutter_.cellLevel();
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 93ce3cee..7e857c0b 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -214,9 +214,11 @@ private:
                 const globalIndex& globalCells,
                 const labelList& globalRegion,
                 const Map<label>& coupledRegionToMaster,
+                const scalarField& cellWeights,
 
                 Map<label>& globalToLocalRegion,
-                pointField& localPoints
+                pointField& localPoints,
+                scalarField& localWeights
             ) const;
 
             //- Convert region into global index.
@@ -574,6 +576,7 @@ public:
             //  (e.g. to keep baffles together)
             labelList decomposeCombineRegions
             (
+                const scalarField& cellWeights,
                 const boolList& blockedFace,
                 const List<labelPair>& explicitConnections,
                 decompositionMethod&
@@ -587,6 +590,7 @@ public:
             (
                 const bool keepZoneFaces,
                 const bool keepBaffles,
+                const scalarField& cellWeights,
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor
             );
@@ -648,7 +652,18 @@ public:
                 const string& msg,
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor,
-                const labelList& cellsToRefine
+                const labelList& cellsToRefine,
+                const scalar maxLoadUnbalance
+            );
+
+            //- Balance before refining some cells
+            autoPtr<mapDistributePolyMesh> balanceAndRefine
+            (
+                const string& msg,
+                decompositionMethod& decomposer,
+                fvMeshDistribute& distributor,
+                const labelList& cellsToRefine,
+                const scalar maxLoadUnbalance
             );
 
 
@@ -699,7 +714,11 @@ public:
             //- Put faces/cells into zones according to surface specification.
             //  Returns null if no zone surfaces present. Region containing
             //  the keepPoint will not be put into a cellZone.
-            autoPtr<mapPolyMesh> zonify(const point& keepPoint);
+            autoPtr<mapPolyMesh> zonify
+            (
+                const point& keepPoint,
+                const bool allowFreeStandingZoneFaces
+            );
 
 
         // Other topo changes
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 88f806b5..f7adea09 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -2032,7 +2032,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
 // Zoning
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 (
-    const point& keepPoint
+    const point& keepPoint,
+    const bool allowFreeStandingZoneFaces
 )
 {
     const wordList& cellZoneNames = surfaces_.cellZoneNames();
@@ -2344,9 +2345,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     }
 
 
-    //// Make sure namedSurfaceIndex is unset inbetween same cell cell zones.
-    //makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
-
+    // Make sure namedSurfaceIndex is unset inbetween same cell cell zones.
+    if (!allowFreeStandingZoneFaces)
+    {
+        makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
+    }
 
     // Topochange container
     polyTopoChange meshMod(mesh_);
@@ -2411,6 +2414,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     }
     syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
 
+    // Get per face whether is it master (of a coupled set of faces)
+    PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
+
     // Set owner as no-flip
     forAll(patches, patchI)
     {
@@ -2428,14 +2434,27 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
                 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
 
                 bool flip;
-                if (ownZone == max(ownZone, neiZone))
+
+                label maxZone = max(ownZone, neiZone);
+
+                if (maxZone == -1)
                 {
                     flip = false;
                 }
-                else
+                else if (ownZone == neiZone)
+                {
+                    // Free-standing zone face or coupled boundary. Keep master
+                    // face unflipped.
+                    flip = !isMasterFace[faceI];
+                }
+                else if (neiZone == maxZone)
                 {
                     flip = true;
                 }
+                else
+                {
+                    flip = false;
+                }
 
                 meshMod.setAction
                 (
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index e35607a5..fde478a9 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -37,6 +37,7 @@ License
 #include "mapDistributePolyMesh.H"
 #include "featureEdgeMesh.H"
 #include "Cloud.H"
+//#include "globalIndex.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -1252,7 +1253,8 @@ Foam::meshRefinement::refineAndBalance
     const string& msg,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor,
-    const labelList& cellsToRefine
+    const labelList& cellsToRefine,
+    const scalar maxLoadUnbalance
 )
 {
     // Do all refinement
@@ -1287,19 +1289,162 @@ Foam::meshRefinement::refineAndBalance
 
     if (Pstream::nProcs() > 1)
     {
-        distMap = balance
+        scalar nIdealCells =
+            mesh_.globalData().nTotalCells()
+          / Pstream::nProcs();
+
+        scalar unbalance = returnReduce
         (
-            false,  //keepZoneFaces
-            false,  //keepBaffles
-            decomposer,
-            distributor
+            mag(1.0-mesh_.nCells()/nIdealCells),
+            maxOp<scalar>()
         );
 
-        Info<< "Balanced mesh in = "
-            << mesh_.time().cpuTimeIncrement() << " s" << endl;
+        if (unbalance <= maxLoadUnbalance)
+        {
+            Info<< "Skipping balancing since max unbalance " << unbalance
+                << " is less than allowable " << maxLoadUnbalance
+                << endl;
+        }
+        else
+        {
+            scalarField cellWeights(mesh_.nCells(), 1);
 
-        printMeshInfo(debug, "After balancing " + msg);
+            distMap = balance
+            (
+                false,  //keepZoneFaces
+                false,  //keepBaffles
+                cellWeights,
+                decomposer,
+                distributor
+            );
+
+            Info<< "Balanced mesh in = "
+                << mesh_.time().cpuTimeIncrement() << " s" << endl;
+
+            printMeshInfo(debug, "After balancing " + msg);
+
+
+            if (debug)
+            {
+                Pout<< "Writing balanced " << msg
+                    << " mesh to time " << timeName() << endl;
+                write
+                (
+                    debug,
+                    mesh_.time().path()/timeName()
+                );
+                Pout<< "Dumped debug data in = "
+                    << mesh_.time().cpuTimeIncrement() << " s" << endl;
+
+                // test all is still synced across proc patches
+                checkData();
+            }
+        }
+    }
+
+    return distMap;
+}
+
+
+// Do load balancing followed by refinement of consistent set of cells.
+Foam::autoPtr<Foam::mapDistributePolyMesh>
+Foam::meshRefinement::balanceAndRefine
+(
+    const string& msg,
+    decompositionMethod& decomposer,
+    fvMeshDistribute& distributor,
+    const labelList& initCellsToRefine,
+    const scalar maxLoadUnbalance
+)
+{
+    labelList cellsToRefine(initCellsToRefine);
+
+    //{
+    //    globalIndex globalCells(mesh_.nCells());
+    //
+    //    Info<< "** Distribution before balancing/refining:" << endl;
+    //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+    //    {
+    //        Info<< "    " << procI << '\t'
+    //            << globalCells.localSize(procI) << endl;
+    //    }
+    //    Info<< endl;
+    //}
+    //{
+    //    globalIndex globalCells(cellsToRefine.size());
+    //
+    //    Info<< "** Cells to be refined:" << endl;
+    //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+    //    {
+    //        Info<< "    " << procI << '\t'
+    //            << globalCells.localSize(procI) << endl;
+    //    }
+    //    Info<< endl;
+    //}
+
+
+    // Load balancing
+    // ~~~~~~~~~~~~~~
+
+    autoPtr<mapDistributePolyMesh> distMap;
+
+    if (Pstream::nProcs() > 1)
+    {
+        // First check if we need to balance at all. Precalculate number of
+        // cells after refinement and see what maximum difference is.
+        scalar nNewCells = scalar(mesh_.nCells() + 7*cellsToRefine.size());
+        scalar nIdealNewCells =
+            returnReduce(nNewCells, sumOp<scalar>())
+          / Pstream::nProcs();
+        scalar unbalance = returnReduce
+        (
+            mag(1.0-nNewCells/nIdealNewCells),
+            maxOp<scalar>()
+        );
+
+        if (unbalance <= maxLoadUnbalance)
+        {
+            Info<< "Skipping balancing since max unbalance " << unbalance
+                << " is less than allowable " << maxLoadUnbalance
+                << endl;
+        }
+        else
+        {
+            scalarField cellWeights(mesh_.nCells(), 1);
+            forAll(cellsToRefine, i)
+            {
+                cellWeights[cellsToRefine[i]] += 7;
+            }
+
+            distMap = balance
+            (
+                false,  //keepZoneFaces
+                false,  //keepBaffles
+                cellWeights,
+                decomposer,
+                distributor
+            );
 
+            // Update cells to refine
+            distMap().distributeCellIndices(cellsToRefine);
+
+            Info<< "Balanced mesh in = "
+                << mesh_.time().cpuTimeIncrement() << " s" << endl;
+        }
+
+        //{
+        //    globalIndex globalCells(mesh_.nCells());
+        //
+        //    Info<< "** Distribution after balancing:" << endl;
+        //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+        //    {
+        //        Info<< "    " << procI << '\t'
+        //            << globalCells.localSize(procI) << endl;
+        //    }
+        //    Info<< endl;
+        //}
+
+        printMeshInfo(debug, "After balancing " + msg);
 
         if (debug)
         {
@@ -1318,6 +1463,46 @@ Foam::meshRefinement::refineAndBalance
         }
     }
 
+
+    // Refinement
+    // ~~~~~~~~~~
+
+    refine(cellsToRefine);
+
+    if (debug)
+    {
+        Pout<< "Writing refined " << msg
+            << " mesh to time " << timeName() << endl;
+        write
+        (
+            debug,
+            mesh_.time().path()
+           /timeName()
+        );
+        Pout<< "Dumped debug data in = "
+            << mesh_.time().cpuTimeIncrement() << " s" << endl;
+
+        // test all is still synced across proc patches
+        checkData();
+    }
+
+    Info<< "Refined mesh in = "
+        << mesh_.time().cpuTimeIncrement() << " s" << endl;
+
+    //{
+    //    globalIndex globalCells(mesh_.nCells());
+    //
+    //    Info<< "** After refinement distribution:" << endl;
+    //    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+    //    {
+    //        Info<< "    " << procI << '\t'
+    //            << globalCells.localSize(procI) << endl;
+    //    }
+    //    Info<< endl;
+    //}
+
+    printMeshInfo(debug, "After refinement " + msg);
+
     return distMap;
 }
 
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index ea2adab0..a23c1ef5 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -338,7 +338,7 @@ Foam::refinementSurfaces::refinementSurfaces
                 {
                     if (regionsDict.found(regionNames[regionI]))
                     {
-                        // Get the dictionary for region 
+                        // Get the dictionary for region
                         const dictionary& regionDict = regionsDict.subDict
                         (
                             regionNames[regionI]
@@ -496,23 +496,63 @@ Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const
 }
 
 
-// Count number of triangles per surface region
-Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s)
-{
-    const geometricSurfacePatchList& regions = s.patches();
-
-    labelList nTris(regions.size(), 0);
-
-    forAll(s, triI)
-    {
-        nTris[s[triI].region()]++;
-    }
-    return nTris;
-}
+// // Count number of triangles per surface region
+// Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s)
+// {
+//     const geometricSurfacePatchList& regions = s.patches();
+// 
+//     labelList nTris(regions.size(), 0);
+// 
+//     forAll(s, triI)
+//     {
+//         nTris[s[triI].region()]++;
+//     }
+//     return nTris;
+// }
+
+
+// // Pre-calculate the refinement level for every element
+// void Foam::refinementSurfaces::wantedRefinementLevel
+// (
+//     const shellSurfaces& shells,
+//     const label surfI,
+//     const List<pointIndexHit>& info,    // Indices
+//     const pointField& ctrs,             // Representative coordinate
+//     labelList& minLevelField
+// ) const
+// {
+//     const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+// 
+//     // Get per element the region
+//     labelList region;
+//     geom.getRegion(info, region);
+// 
+//     // Initialise fields to region wise minLevel
+//     minLevelField.setSize(ctrs.size());
+//     minLevelField = -1;
+// 
+//     forAll(minLevelField, i)
+//     {
+//         if (info[i].hit())
+//         {
+//             minLevelField[i] = minLevel(surfI, region[i]);
+//         }
+//     }
+// 
+//     // Find out if triangle inside shell with higher level
+//     // What level does shell want to refine fc to?
+//     labelList shellLevel;
+//     shells.findHigherLevel(ctrs, minLevelField, shellLevel);
+// 
+//     forAll(minLevelField, i)
+//     {
+//         minLevelField[i] = max(minLevelField[i], shellLevel[i]);
+//     }
+// }
 
 
 // Precalculate the refinement level for every element of the searchable
-// surface. This is currently hardcoded for triSurfaceMesh only.
+// surface.
 void Foam::refinementSurfaces::setMinLevelFields
 (
     const shellSurfaces& shells
@@ -522,58 +562,55 @@ void Foam::refinementSurfaces::setMinLevelFields
     {
         const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
 
-        if (isA<triSurfaceMesh>(geom))
+        // Precalculation only makes sense if there are different regions
+        // (so different refinement levels possible) and there are some
+        // elements. Possibly should have 'enough' elements to have fine
+        // enough resolution but for now just make sure we don't catch e.g.
+        // searchableBox (size=6)
+        if (geom.regions().size() > 1 && geom.globalSize() > 10)
         {
-            const triSurfaceMesh& triMesh = refCast<const triSurfaceMesh>(geom);
+            // Representative local coordinates
+            const pointField ctrs = geom.coordinates();
 
-            autoPtr<triSurfaceLabelField> minLevelFieldPtr
-            (
-                new triSurfaceLabelField
+            labelList minLevelField(ctrs.size(), -1);
+            {
+                // Get the element index in a roundabout way. Problem is e.g.
+                // distributed surface where local indices differ from global
+                // ones (needed for getRegion call)
+                List<pointIndexHit> info;
+                geom.findNearest
                 (
-                    IOobject
-                    (
-                        "values",
-                        triMesh.objectRegistry::time().timeName(),  // instance
-                        "triSurface",                               // local
-                        triMesh,
-                        IOobject::NO_READ,
-                        IOobject::AUTO_WRITE
-                    ),
-                    triMesh,
-                    dimless
-                )
-            );
-            triSurfaceLabelField& minLevelField = minLevelFieldPtr();
+                    ctrs,
+                    scalarField(ctrs.size(), sqr(GREAT)),
+                    info
+                );
 
-            const triSurface& s = static_cast<const triSurface&>(triMesh);
+                // Get per element the region
+                labelList region;
+                geom.getRegion(info, region);
 
-            // Initialise fields to region wise minLevel
-            forAll(s, triI)
-            {
-                minLevelField[triI] = minLevel(surfI, s[triI].region());
+                // From the region get the surface-wise refinement level
+                forAll(minLevelField, i)
+                {
+                    if (info[i].hit())
+                    {
+                        minLevelField[i] = minLevel(surfI, region[i]);
+                    }
+                }
             }
 
             // Find out if triangle inside shell with higher level
-            pointField fc(s.size());
-            forAll(s, triI)
-            {
-                fc[triI] = s[triI].centre(s.points());
-            }
             // What level does shell want to refine fc to?
             labelList shellLevel;
-            shells.findHigherLevel(fc, minLevelField, shellLevel);
+            shells.findHigherLevel(ctrs, minLevelField, shellLevel);
 
-            forAll(minLevelField, triI)
+            forAll(minLevelField, i)
             {
-                minLevelField[triI] = max
-                (
-                    minLevelField[triI],
-                    shellLevel[triI]
-                );
+                minLevelField[i] = max(minLevelField[i], shellLevel[i]);
             }
 
-            // Store field on triMesh
-            minLevelFieldPtr.ptr()->store();
+            // Store minLevelField on surface
+            const_cast<searchableSurface&>(geom).setField(minLevelField);
         }
     }
 }
@@ -601,44 +638,88 @@ void Foam::refinementSurfaces::findHigherIntersection
         return;
     }
 
-    // Work arrays
-    labelList hitMap(identity(start.size()));
-    pointField p0(start);
-    pointField p1(end);
-    List<pointIndexHit> hitInfo(start.size());
-
-    forAll(surfaces_, surfI)
+    if (surfaces_.size() == 1)
     {
+        // Optimisation: single segmented surface. No need to duplicate
+        // point storage.
+
+        label surfI = 0;
+
         const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
 
-        geom.findLineAny(p0, p1, hitInfo);
+        // Do intersection test
+        List<pointIndexHit> intersectionInfo(start.size());
+        geom.findLineAny(start, end, intersectionInfo);
 
+        // See if a cached level field available
         labelList minLevelField;
-        if (isA<triSurfaceMesh>(geom))
+        geom.getField(intersectionInfo, minLevelField);
+        bool haveLevelField =
+        (
+            returnReduce(minLevelField.size(), sumOp<label>())
+          > 0
+        );
+
+        if (!haveLevelField && geom.regions().size() == 1)
         {
-            const triSurfaceMesh& triMesh = refCast<const triSurfaceMesh>(geom);
-            triMesh.getField
+            minLevelField = labelList
             (
-                //"minLevel",
-                hitInfo,
-                minLevelField
+                intersectionInfo.size(),
+                minLevel(surfI, 0)
             );
+            haveLevelField = true;
+        }
+
+        if (haveLevelField)
+        {
+            forAll(intersectionInfo, i)
+            {
+                if
+                (
+                    intersectionInfo[i].hit()
+                 && minLevelField[i] > currentLevel[i]
+                )
+                {
+                    surfaces[i] = surfI;    // index of surface
+                    surfaceLevel[i] = minLevelField[i];
+                }
+            }
+            return;
         }
+    }
+
 
 
-        // Copy all hits into arguments, continue with misses
-        label newI = 0;
-        forAll(hitInfo, hitI)
+    // Work arrays
+    pointField p0(start);
+    pointField p1(end);
+    labelList intersectionToPoint(identity(start.size()));
+    List<pointIndexHit> intersectionInfo(start.size());
+
+    forAll(surfaces_, surfI)
+    {
+        const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+
+        // Do intersection test
+        geom.findLineAny(p0, p1, intersectionInfo);
+
+        // See if a cached level field available
+        labelList minLevelField;
+        geom.getField(intersectionInfo, minLevelField);
+
+        // Copy all hits into arguments, In-place compact misses.
+        label missI = 0;
+        forAll(intersectionInfo, i)
         {
             // Get the minLevel for the point
             label minLocalLevel = -1;
 
-            if (hitInfo[hitI].hit())
+            if (intersectionInfo[i].hit())
             {
                 // Check if minLevelField for this surface.
                 if (minLevelField.size())
                 {
-                    minLocalLevel = minLevelField[hitI];
+                    minLocalLevel = minLevelField[i];
                 }
                 else
                 {
@@ -648,36 +729,35 @@ void Foam::refinementSurfaces::findHigherIntersection
                 }
             }
 
-            label pointI = hitMap[hitI];
+
+            label pointI = intersectionToPoint[i];
 
             if (minLocalLevel > currentLevel[pointI])
             {
+                // Mark point for refinement
                 surfaces[pointI] = surfI;
                 surfaceLevel[pointI] = minLocalLevel;
             }
             else
             {
-                if (hitI != newI)
-                {
-                    hitMap[newI] = hitMap[hitI];
-                    p0[newI] = p0[hitI];
-                    p1[newI] = p1[hitI];
-                }
-                newI++;
+                p0[missI] = start[pointI];
+                p1[missI] = end[pointI];
+                intersectionToPoint[missI] = pointI;
+                missI++;
             }
         }
 
         // All done? Note that this decision should be synchronised
-        if (returnReduce(newI, sumOp<label>()) == 0)
+        if (returnReduce(missI, sumOp<label>()) == 0)
         {
             break;
         }
 
-        // Trim and continue
-        hitMap.setSize(newI);
-        p0.setSize(newI);
-        p1.setSize(newI);
-        hitInfo.setSize(newI);
+        // Trim misses
+        p0.setSize(missI);
+        p1.setSize(missI);
+        intersectionToPoint.setSize(missI);
+        intersectionInfo.setSize(missI);
     }
 }
 
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index cee48bbe..b1188f1d 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -218,8 +218,8 @@ public:
                 const shellSurfaces& shells
             );
 
-            //- Helper: count number of triangles per region
-            static labelList countRegions(const triSurface&);
+            ////- Helper: count number of triangles per region
+            //static labelList countRegions(const triSurface&);
 
 
         // Searching
@@ -291,7 +291,7 @@ public:
             (
                 const labelList& surfacesToTest,
                 const pointField& samples,
-                const  scalarField& nearestDistSqr,
+                const scalarField& nearestDistSqr,
                 labelList& hitSurface,
                 labelList& hitRegion
             ) const;
diff --git a/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C b/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C
index 237e3527..0bda7cd1 100644
--- a/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C
+++ b/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C
@@ -39,7 +39,7 @@ Foam::label Foam::ExactParticle<ParticleType>::track
     this->facei_ = -1;
 
     // Tracks to endPosition or stop on boundary
-    while(!this->onBoundary() && this->stepFraction_ < 1.0 - SMALL)
+    while (!this->onBoundary() && this->stepFraction_ < 1.0 - SMALL)
     {
         this->stepFraction_ +=
             trackToFace(endPosition, td)
diff --git a/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.H b/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.H
index d62c7e08..64bcbd3b 100644
--- a/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.H
+++ b/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.H
@@ -84,7 +84,7 @@ public:
         :
             Particle<ParticleType>(cloud, position, celli)
         {}
-            
+
 
         //- Construct from Istream
         ExactParticle
diff --git a/src/meshTools/searchableSurface/searchableSurfaceCollection.C b/src/meshTools/searchableSurface/searchableSurfaceCollection.C
index 9279087c..bef7f79e 100644
--- a/src/meshTools/searchableSurface/searchableSurfaceCollection.C
+++ b/src/meshTools/searchableSurface/searchableSurfaceCollection.C
@@ -681,19 +681,22 @@ void Foam::searchableSurfaceCollection::getField
         List<List<label> > infoMap;
         sortHits(info, surfInfo, infoMap);
 
-        values.setSize(info.size());
-        //?Misses do not get set? values = 0;
-
         // Do surface tests
         forAll(surfInfo, surfI)
         {
             labelList surfValues;
             subGeom_[surfI].getField(surfInfo[surfI], surfValues);
 
-            const labelList& map = infoMap[surfI];
-            forAll(map, i)
+            if (surfValues.size())
             {
-                values[map[i]] = surfValues[i];
+                // Size values only when we have a surface that supports it.
+                values.setSize(info.size());
+
+                const labelList& map = infoMap[surfI];
+                forAll(map, i)
+                {
+                    values[map[i]] = surfValues[i];
+                }
             }
         }
     }
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
index 714a95eb..71714217 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
@@ -116,9 +116,9 @@ castellatedMeshControls
     // Refinement parameters
     // ~~~~~~~~~~~~~~~~~~~~~
 
-    // While refining maximum number of cells per processor. This is basically
-    // the number of cells that fit on a processor. If you choose this too small
-    // it will do just more refinement iterations to obtain a similar mesh.
+    // If local number of cells is >= maxLocalCells on any processor
+    // switches from from refinement followed by balancing
+    // (current method) to (weighted) balancing before refinement.
     maxLocalCells 1000000;
 
     // Overall cell limit (approximately). Refinement will stop immediately
@@ -249,6 +249,8 @@ snapControls
 // Settings for the layer addition.
 addLayersControls
 {
+    // Are the thickness parameters below relative to the undistorted
+    // size of the refined cell outside layer (true) or absolute sizes (false).
     relativeSizes true;
 
     // Per final patch (so not geometry!) the layer information
@@ -268,14 +270,13 @@ addLayersControls
     expansionRatio 1.0;
 
     //- Wanted thickness of final added cell layer. If multiple layers
-    //  is the
-    //  thickness of the layer furthest away from the wall.
-    //  Relative to undistorted size of cell outside layer.
+    //  is the thickness of the layer furthest away from the wall.
+    //  See relativeSizes parameter.
     finalLayerThickness 0.5;
 
     //- Minimum thickness of cell layer. If for any reason layer
     //  cannot be above minThickness do not add layer.
-    //  Relative to undistorted size of cell outside layer.
+    //  See relativeSizes parameter.
     minThickness 0.25;
 
     //- If points get not extruded do nGrow layers of connected faces that are
@@ -316,7 +317,10 @@ addLayersControls
     // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
 
-    // Overall max number of layer addition iterations
+
+    // Overall max number of layer addition iterations. The mesher will exit
+    // if it reaches this number of iterations; possibly with an illegal
+    // mesh.
     nLayerIter 50;
 }
 
@@ -342,6 +346,7 @@ meshQualityControls
     minFlatness 0.5;
 
     //- Minimum pyramid volume. Is absolute volume of cell pyramid.
+    //  Set to a sensible fraction of the smallest cell volume expected.
     //  Set to very negative number (e.g. -1E30) to disable.
     minVol 1e-13;
 
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict b/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
index 70593293..fefd88a4 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
@@ -51,9 +51,9 @@ castellatedMeshControls
     // Refinement parameters
     // ~~~~~~~~~~~~~~~~~~~~~
 
-    // While refining maximum number of cells per processor. This is basically
-    // the number of cells that fit on a processor. If you choose this too small
-    // it will do just more refinement iterations to obtain a similar mesh.
+    // If local number of cells is >= maxLocalCells on any processor
+    // switches from from refinement followed by balancing
+    // (current method) to (weighted) balancing before refinement.
     maxLocalCells 1000000;
 
     // Overall cell limit (approximately). Refinement will stop immediately
@@ -69,6 +69,13 @@ castellatedMeshControls
     // (unless the number of cells to refine is 0)
     minRefinementCells 10;
 
+    // Allow a certain level of imbalance during refining
+    // (since balancing is quite expensive)
+    // Expressed as fraction of perfect balance (= overall number of cells /
+    // nProcs). 0=balance always.
+    maxLoadUnbalance 0.10;
+
+
     // Number of buffer layers between different levels.
     // 1 means normal 2:1 refinement restriction, larger means slower
     // refinement.
@@ -174,6 +181,8 @@ snapControls
 // Settings for the layer addition.
 addLayersControls
 {
+    // Are the thickness parameters below relative to the undistorted
+    // size of the refined cell outside layer (true) or absolute sizes (false).
     relativeSizes true;
 
     // Per final patch (so not geometry!) the layer information
@@ -457,9 +466,8 @@ addLayersControls
     expansionRatio 1.0;
 
     //- Wanted thickness of final added cell layer. If multiple layers
-    //  is the
-    //  thickness of the layer furthest away from the wall.
-    //  Relative to undistorted size of cell outside layer.
+    //  is the thickness of the layer furthest away from the wall.
+    //  See relativeSizes parameter.
     finalLayerThickness 0.3;
 
     //- Minimum thickness of cell layer. If for any reason layer
@@ -531,6 +539,7 @@ meshQualityControls
     minFlatness 0.5;
 
     //- Minimum pyramid volume. Is absolute volume of cell pyramid.
+    //  Set to a sensible fraction of the smallest cell volume expected.
     //  Set to very negative number (e.g. -1E30) to disable.
     minVol 1e-13;
 
-- 
2.11.4.GIT