From 2e558acf65d9978b17c48e6f267a091b89e5df5a Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 22 Sep 2009 15:24:59 +0100
Subject: [PATCH 1/7] TypeName instead of ClassName

---
 applications/test/dictionary/calcEntry/calcEntry.H | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/applications/test/dictionary/calcEntry/calcEntry.H b/applications/test/dictionary/calcEntry/calcEntry.H
index e8a0925f94..119f155643 100644
--- a/applications/test/dictionary/calcEntry/calcEntry.H
+++ b/applications/test/dictionary/calcEntry/calcEntry.H
@@ -64,7 +64,7 @@ class calcEntry
 public:
 
     //- Runtime type information
-    TypeName("calc");
+    ClassName("calc");
 
 
     // Member Functions

From 859d727696e1acf30fcd48a34052756b1757f024 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 22 Sep 2009 15:25:16 +0100
Subject: [PATCH 2/7] added patch non-manifold edge dumping

---
 .../conversion/writeMeshObj/writeMeshObj.C    | 77 ++++++++++++++++---
 1 file changed, 66 insertions(+), 11 deletions(-)

diff --git a/applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C b/applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C
index e1d571e081..d2b23fef16 100644
--- a/applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C
+++ b/applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C
@@ -32,7 +32,8 @@ Description
 
     patch_YYY_XXX.obj : all face centres of patch YYY
 
-    Optional: patch faces (as polygons) : patchFaces_YYY_XXX.obj
+    Optional: - patch faces (as polygons) : patchFaces_YYY_XXX.obj
+              - non-manifold edges : patchEdges_YYY_XXX.obj
 
 \*---------------------------------------------------------------------------*/
 
@@ -51,7 +52,7 @@ using namespace Foam;
 
 void writeOBJ(const point& pt, Ostream& os)
 {
-    os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
+    os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
 }
 
 // All edges of mesh
@@ -75,8 +76,7 @@ void writePoints(const polyMesh& mesh, const fileName& timeName)
     {
         const edge& e = mesh.edges()[edgeI];
 
-        pointStream << "l " << e.start() + 1 << ' ' << e.end() + 1
-            << endl;
+        pointStream << "l " << e.start() + 1 << ' ' << e.end() + 1 << nl;
     }
 }
 
@@ -277,7 +277,47 @@ void writePatchFaces
             {
                 patchFaceStream << ' ' << f[fp]+1;
             }
-            patchFaceStream << endl;
+            patchFaceStream << nl;
+        }
+    }
+}
+
+
+void writePatchBoundaryEdges
+(
+    const polyMesh& mesh,
+    const fileName& timeName
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        fileName edgeFile
+        (
+            mesh.time().path()
+          / "patchEdges_" + pp.name() + '_' + timeName + ".obj"
+        );
+
+        Info << "Writing patch edges to " << edgeFile << endl;
+
+        OFstream patchEdgeStream(edgeFile);
+
+        forAll(pp.localPoints(), pointI)
+        {
+            writeOBJ(pp.localPoints()[pointI], patchEdgeStream);
+        }
+
+        for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
+        {
+            if (pp.edgeFaces()[edgeI].size() == 1)
+            {
+                const edge& e = pp.edges()[edgeI];
+
+                patchEdgeStream<< "l " << e[0]+1 << ' ' << e[1]+1 << nl;
+            }
         }
     }
 }
@@ -339,6 +379,7 @@ int main(int argc, char *argv[])
 {
     timeSelector::addOptions();
     argList::validOptions.insert("patchFaces", "");
+    argList::validOptions.insert("patchEdges", "");
     argList::validOptions.insert("cell", "cellI");
     argList::validOptions.insert("face", "faceI");
     argList::validOptions.insert("point", "pointI");
@@ -351,6 +392,7 @@ int main(int argc, char *argv[])
     runTime.functionObjects().off();
 
     bool patchFaces = args.optionFound("patchFaces");
+    bool patchEdges = args.optionFound("patchEdges");
     bool doCell     = args.optionFound("cell");
     bool doPoint    = args.optionFound("point");
     bool doFace     = args.optionFound("face");
@@ -381,19 +423,23 @@ int main(int argc, char *argv[])
             {
                 writePatchFaces(mesh, runTime.timeName());
             }
-            else if (doCell)
+            if (patchEdges)
+            {
+                writePatchBoundaryEdges(mesh, runTime.timeName());
+            }
+            if (doCell)
             {
                 label cellI = args.optionRead<label>("cell");
 
                 writePoints(mesh, cellI, runTime.timeName());
             }
-            else if (doPoint)
+            if (doPoint)
             {
                 label pointI = args.optionRead<label>("point");
 
                 writePointCells(mesh, pointI, runTime.timeName());
             }
-            else if (doFace)
+            if (doFace)
             {
                 label faceI = args.optionRead<label>("face");
 
@@ -415,7 +461,7 @@ int main(int argc, char *argv[])
 
                 meshTools::writeOBJ(str, faceList(1, f), mesh.points());
             }
-            else if (doCellSet)
+            if (doCellSet)
             {
                 word setName(args.option("cellSet"));
 
@@ -427,7 +473,7 @@ int main(int argc, char *argv[])
                 writePoints(mesh, cells.toc(), runTime.timeName());
 
             }
-            else if (doFaceSet)
+            if (doFaceSet)
             {
                 word setName(args.option("faceSet"));
 
@@ -458,7 +504,16 @@ int main(int argc, char *argv[])
                     faces.toc()
                 );
             }
-            else
+            else if
+            (
+                !patchFaces 
+             && !patchEdges
+             && !doCell
+             && !doPoint
+             && !doFace
+             && !doCellSet
+             && !doFaceSet
+            )
             {
                 // points & edges
                 writePoints(mesh, runTime.timeName());

From 1bdfa97fc73239f7838ad3115b0696cdbdc2e065 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 22 Sep 2009 15:25:53 +0100
Subject: [PATCH 3/7] adapted extrudeProperties for new extrusion model

---
 .../utilities/mesh/generation/extrudeMesh/extrudeProperties  | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties b/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties
index faf875774d..9a047e7822 100644
--- a/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties
+++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties
@@ -71,7 +71,8 @@ linearNormalCoeffs
 
 linearDirectionCoeffs
 {
-    direction       (0 0 1);
+    direction       (0 1 0);
+    thickness       0.05;
 }
 
 linearRadialCoeffs
@@ -89,7 +90,7 @@ sigmaRadialCoeffs
 
 // Do front and back need to be merged? Usually only makes sense for 360
 // degree wedges.
-mergeFaces true;
+mergeFaces false;   //true;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

From da8198a4876506f8cfced3c603033132fc9dd710 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 22 Sep 2009 15:26:19 +0100
Subject: [PATCH 4/7] added region option

---
 .../mesh/manipulation/mergeMeshes/mergeMeshes.C      |  4 ++--
 .../mesh/manipulation/mergeMeshes/setRoots.H         | 12 ++++++++++--
 .../patch/patchIntegrate/patchIntegrate.C            |  3 ++-
 3 files changed, 14 insertions(+), 5 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C b/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C
index 28955c4d6f..b13dfb86ef 100644
--- a/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C
+++ b/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C
@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
     (
         IOobject
         (
-            mergePolyMesh::defaultRegion,
+            masterRegion,
             runTimeMaster.timeName(),
             runTimeMaster
         )
@@ -63,7 +63,7 @@ int main(int argc, char *argv[])
     (
         IOobject
         (
-            mergePolyMesh::defaultRegion,
+            addRegion,
             runTimeToAdd.timeName(),
             runTimeToAdd
         )
diff --git a/applications/utilities/mesh/manipulation/mergeMeshes/setRoots.H b/applications/utilities/mesh/manipulation/mergeMeshes/setRoots.H
index 7ddd408965..cfbf59cf64 100644
--- a/applications/utilities/mesh/manipulation/mergeMeshes/setRoots.H
+++ b/applications/utilities/mesh/manipulation/mergeMeshes/setRoots.H
@@ -2,9 +2,11 @@
 
     argList::validArgs.append("master root");
     argList::validArgs.append("master case");
+    argList::validOptions.insert("masterRegion", "name");
 
     argList::validArgs.append("root to add");
     argList::validArgs.append("case to add");
+    argList::validOptions.insert("addRegion", "name");
 
     argList args(argc, argv);
 
@@ -15,9 +17,15 @@
 
     fileName rootDirMaster(args.additionalArgs()[0]);
     fileName caseDirMaster(args.additionalArgs()[1]);
+    word masterRegion = polyMesh::defaultRegion;
+    args.optionReadIfPresent("masterRegion", masterRegion);
 
     fileName rootDirToAdd(args.additionalArgs()[2]);
     fileName caseDirToAdd(args.additionalArgs()[3]);
+    word addRegion = polyMesh::defaultRegion;
+    args.optionReadIfPresent("addRegion", addRegion);
 
-    Info<< "Master:      " << rootDirMaster << " " << caseDirMaster << nl
-        << "mesh to add: " << rootDirToAdd << " " << caseDirToAdd << endl;
+    Info<< "Master:      " << rootDirMaster << " " << caseDirMaster
+        << "  region " << masterRegion << nl
+        << "mesh to add: " << rootDirToAdd << " " << caseDirToAdd
+        << "  region " << addRegion << endl;
diff --git a/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C b/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C
index 054c09cac5..8d189dfe28 100644
--- a/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C
+++ b/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C
@@ -38,13 +38,14 @@ Description
 
 int main(int argc, char *argv[])
 {
+#   include "addRegionOption.H"
     timeSelector::addOptions();
     argList::validArgs.append("fieldName");
     argList::validArgs.append("patchName");
 #   include "setRootCase.H"
 #   include "createTime.H"
     instantList timeDirs = timeSelector::select0(runTime, args);
-#   include "createMesh.H"
+#   include "createNamedMesh.H"
 
     word fieldName(args.additionalArgs()[0]);
     word patchName(args.additionalArgs()[1]);

From 1c8ad8914599b86a041698a1bea7f9cf57420324 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 22 Sep 2009 15:26:40 +0100
Subject: [PATCH 5/7] added cellZonesOnly option which does not use walk

---
 .../splitMeshRegions/splitMeshRegions.C       | 85 +++++++++++++------
 1 file changed, 60 insertions(+), 25 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index a8d1ff471e..be267f8724 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -43,6 +43,9 @@ Description
     the largest matching region (-sloppyCellZones). This will accept any
     region that covers more than 50% of the zone. It has to be a subset
     so cannot have any cells in any other zone.
+    - useCellZonesOnly does not do a walk and uses the cellZones only. Use
+    this if you don't mind having disconnected domains in a single region.
+    This option requires all cells to be in one (and one only) region.
 
 \*---------------------------------------------------------------------------*/
 
@@ -609,7 +612,7 @@ void getInterfaceSizes
 // Create mesh for region.
 autoPtr<mapPolyMesh> createRegionMesh
 (
-    const regionSplit& cellRegion,
+    const labelList& cellRegion,
     const EdgeMap<label>& interfaceToPatch,
     const fvMesh& mesh,
     const label regionI,
@@ -766,7 +769,7 @@ autoPtr<mapPolyMesh> createRegionMesh
 void createAndWriteRegion
 (
     const fvMesh& mesh,
-    const regionSplit& cellRegion,
+    const labelList& cellRegion,
     const wordList& regionNames,
     const EdgeMap<label>& interfaceToPatch,
     const label regionI,
@@ -1005,7 +1008,8 @@ void createAndWriteRegion
 EdgeMap<label> addRegionPatches
 (
     fvMesh& mesh,
-    const regionSplit& cellRegion,
+    const labelList& cellRegion,
+    const label nCellRegions,
     const edgeList& interfaces,
     const EdgeMap<label>& interfaceSizes,
     const wordList& regionNames
@@ -1016,7 +1020,7 @@ EdgeMap<label> addRegionPatches
 
     Info<< nl << "Adding patches" << nl << endl;
 
-    EdgeMap<label> interfaceToPatch(cellRegion.nRegions());
+    EdgeMap<label> interfaceToPatch(nCellRegions);
 
     forAll(interfaces, interI)
     {
@@ -1108,13 +1112,14 @@ EdgeMap<label> addRegionPatches
 label findCorrespondingRegion
 (
     const labelList& existingZoneID,    // per cell the (unique) zoneID
-    const regionSplit& cellRegion,
+    const labelList& cellRegion,
+    const label nCellRegions,
     const label zoneI,
     const label minOverlapSize
 )
 {
     // Per region the number of cells in zoneI
-    labelList cellsInZone(cellRegion.nRegions(), 0);
+    labelList cellsInZone(nCellRegions, 0);
 
     forAll(cellRegion, cellI)
     {
@@ -1162,7 +1167,8 @@ label findCorrespondingRegion
 //(
 //    const cellZoneMesh& cellZones,
 //    const labelList& existingZoneID,    // per cell the (unique) zoneID
-//    const regionSplit& cellRegion,
+//    const labelList& cellRegion,
+//    const label nCellRegions,
 //    const label zoneI
 //)
 //{
@@ -1227,6 +1233,7 @@ label findCorrespondingRegion
 int main(int argc, char *argv[])
 {
     argList::validOptions.insert("cellZones", "");
+    argList::validOptions.insert("cellZonesOnly", "");
     argList::validOptions.insert("blockedFaces", "faceSet");
     argList::validOptions.insert("makeCellZones", "");
     argList::validOptions.insert("largestOnly", "");
@@ -1249,13 +1256,14 @@ int main(int argc, char *argv[])
             << blockedFacesName << nl << endl;
     }
 
-    bool makeCellZones   = args.optionFound("makeCellZones");
-    bool largestOnly     = args.optionFound("largestOnly");
-    bool insidePoint     = args.optionFound("insidePoint");
-    bool useCellZones    = args.optionFound("cellZones");
-    bool overwrite       = args.optionFound("overwrite");
-    bool detectOnly      = args.optionFound("detectOnly");
-    bool sloppyCellZones = args.optionFound("sloppyCellZones");
+    bool makeCellZones    = args.optionFound("makeCellZones");
+    bool largestOnly      = args.optionFound("largestOnly");
+    bool insidePoint      = args.optionFound("insidePoint");
+    bool useCellZones     = args.optionFound("cellZones");
+    bool useCellZonesOnly = args.optionFound("cellZonesOnly");
+    bool overwrite        = args.optionFound("overwrite");
+    bool detectOnly       = args.optionFound("detectOnly");
+    bool sloppyCellZones  = args.optionFound("sloppyCellZones");
 
     if (insidePoint && largestOnly)
     {
@@ -1370,13 +1378,37 @@ int main(int argc, char *argv[])
     }
 
 
-    // Do the topological walk to determine regions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Determine per cell the region it belongs to
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // regionSplit is the labelList with the region per cell.
-    regionSplit cellRegion(mesh, blockedFace);
+    // cellRegion is the labelList with the region per cell.
+    labelList cellRegion;
+    label nCellRegions = 0;
+    if (useCellZonesOnly)
+    {
+        label unzonedCellI = findIndex(zoneID, -1);
+        if (unzonedCellI != -1)
+        {
+            FatalErrorIn(args.executable())
+                << "For the cellZonesOnly option all cells "
+                << "have to be in a cellZone." << endl
+                << "Cell " << unzonedCellI
+                << " at" << mesh.cellCentres()[unzonedCellI]
+                << " is not in a cellZone. There might be more unzoned cells."
+                << exit(FatalError);
+        }
+        cellRegion = zoneID;
+        nCellRegions = gMax(cellRegion)+1;
+    }
+    else
+    {
+        // Do a topological walk to determine regions
+        regionSplit regions(mesh, blockedFace);
+        nCellRegions = regions.nRegions();
+        cellRegion.transfer(regions);
+    }
 
-    Info<< endl << "Number of regions:" << cellRegion.nRegions() << nl << endl;
+    Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
 
 
     // Write to manual decomposition option
@@ -1429,7 +1461,7 @@ int main(int argc, char *argv[])
     // Sizes per region
     // ~~~~~~~~~~~~~~~~
 
-    labelList regionSizes(cellRegion.nRegions(), 0);
+    labelList regionSizes(nCellRegions, 0);
 
     forAll(cellRegion, cellI)
     {
@@ -1489,9 +1521,9 @@ int main(int argc, char *argv[])
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Region per zone
-    labelList regionToZone(cellRegion.nRegions(), -1);
+    labelList regionToZone(nCellRegions, -1);
     // Name of region
-    wordList regionNames(cellRegion.nRegions());
+    wordList regionNames(nCellRegions);
     // Zone to region
     labelList zoneToRegion(cellZones.size(), -1);
 
@@ -1506,6 +1538,7 @@ int main(int argc, char *argv[])
             (
                 zoneID,
                 cellRegion,
+                nCellRegions,
                 zoneI,
                 label(0.5*zoneSizes[zoneI]) // minimum overlap
             );
@@ -1532,6 +1565,7 @@ int main(int argc, char *argv[])
             (
                 zoneID,
                 cellRegion,
+                nCellRegions,
                 zoneI,
                 1               // minimum overlap
             );
@@ -1660,7 +1694,7 @@ int main(int argc, char *argv[])
     mesh.clearOut();
 
 
-    if (cellRegion.nRegions() == 1)
+    if (nCellRegions == 1)
     {
         Info<< "Only one region. Doing nothing." << endl;
     }
@@ -1671,7 +1705,7 @@ int main(int argc, char *argv[])
 
         // Check if region overlaps with existing zone. If so keep.
 
-        for (label regionI = 0; regionI < cellRegion.nRegions(); regionI++)
+        for (label regionI = 0; regionI < nCellRegions; regionI++)
         {
             label zoneI = regionToZone[regionI];
 
@@ -1759,6 +1793,7 @@ int main(int argc, char *argv[])
             (
                 mesh,
                 cellRegion,
+                nCellRegions,
                 interfaces,
                 interfaceSizes,
                 regionNames
@@ -1837,7 +1872,7 @@ int main(int argc, char *argv[])
         else
         {
             // Split all
-            for (label regionI = 0; regionI < cellRegion.nRegions(); regionI++)
+            for (label regionI = 0; regionI < nCellRegions; regionI++)
             {
                 Info<< nl
                     << "Region " << regionI << nl

From 2f9e3fa6a9f82019a7d6cce1dbe1540dda52f529 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 22 Sep 2009 15:26:57 +0100
Subject: [PATCH 6/7] disable wildcard matching when expanding variables

---
 .../db/dictionary/primitiveEntry/primitiveEntryIO.C        | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C b/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C
index 896dee8d14..64810f54f4 100644
--- a/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C
+++ b/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C
@@ -81,7 +81,12 @@ bool Foam::primitiveEntry::expandVariable
     word varName = w(1, w.size()-1);
 
     // lookup the variable name in the given dictionary....
-    const entry* ePtr = dict.lookupEntryPtr(varName, true, true);
+    // Note: allow wildcards to match? For now disabled since following
+    // would expand internalField to wildcard match and not expected
+    // internalField:
+    //      internalField XXX;
+    //      boundaryField { ".*" {YYY;} movingWall {value $internalField;}
+    const entry* ePtr = dict.lookupEntryPtr(varName, true, false);
 
     // ...if defined insert its tokens into this
     if (ePtr != NULL)

From 909b590652907c77ae7aad0c57ee489c53ffa003 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 22 Sep 2009 15:27:32 +0100
Subject: [PATCH 7/7] balance before refinement

---
 .../autoHexMeshDriver/autoLayerDriver.C       | 228 ++++++++++-------
 .../autoHexMeshDriver/autoLayerDriver.H       |   2 +-
 .../autoHexMeshDriver/autoRefineDriver.C      |  15 +-
 .../refinementParameters.C                    |   6 +-
 .../refinementParameters.H                    |   9 +
 .../meshRefinement/meshRefinement.C           |  37 ++-
 .../meshRefinement/meshRefinement.H           |   9 +-
 .../meshRefinement/meshRefinementRefine.C     | 231 +++++++++++++++---
 8 files changed, 393 insertions(+), 144 deletions(-)

diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 3d0274bc33..74dc49c516 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -2509,43 +2509,55 @@ 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
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    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
+    autoPtr<indirectPrimitivePatch> pp
+    (
+        meshRefinement::makePatch
         (
-            mesh.edges(),
-            mesh.pointEdges()[v0],
-            v0,
-            v1
-        );
-    }
+            mesh,
+            patchIDs
+        )
+    );
+
+    // Construct iterative mesh mover.
+    Info<< "Constructing mesh displacer ..." << endl;
+
+    autoPtr<motionSmoother> meshMover
+    (
+        new motionSmoother
+        (
+            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);
+    List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
 
 
     // Get number of layer per point from number of layers per patch
@@ -2554,7 +2566,7 @@ void Foam::autoLayerDriver::addLayers
     setNumLayers
     (
         layerParams.numLayers(),    // per patch the num layers
-        meshMover.adaptPatchIDs(),  // patches that are being moved
+        meshMover().adaptPatchIDs(),// patches that are being moved
         pp,                         // indirectpatch for all faces moving
 
         patchDisp,
@@ -2562,6 +2574,9 @@ void Foam::autoLayerDriver::addLayers
         extrudeStatus
     );
 
+    // Precalculate mesh edge labels for patch edges
+    labelList meshEdges(pp().meshEdges(mesh.edges(), mesh.pointEdges()));
+
     // Disable extrusion on non-manifold points
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -2582,7 +2597,7 @@ void Foam::autoLayerDriver::addLayers
     (
         pp,
         meshEdges,
-        layerParams.featureAngle()*constant::math::pi/180.0,
+    layerParams.featureAngle()*constant::math::pi/180.0,
 
         patchDisp,
         patchNLayers,
@@ -2633,14 +2648,15 @@ void Foam::autoLayerDriver::addLayers
         );
     }
 
+
     // 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?
@@ -2663,11 +2679,11 @@ 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 
@@ -2678,9 +2694,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();
 
@@ -2691,7 +2707,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]);
@@ -2761,7 +2777,7 @@ void Foam::autoLayerDriver::addLayers
             IOobject::NO_WRITE,
             false
         ),
-        meshMover.pMesh(),
+        meshMover().pMesh(),
         dimensionedScalar("pointMedialDist", dimless, 0.0)
     );
 
@@ -2776,7 +2792,7 @@ void Foam::autoLayerDriver::addLayers
             IOobject::NO_WRITE,
             false
         ),
-        meshMover.pMesh(),
+        meshMover().pMesh(),
         dimensionedVector("dispVec", dimless, vector::zero)
     );
 
@@ -2791,7 +2807,7 @@ void Foam::autoLayerDriver::addLayers
             IOobject::NO_WRITE,
             false
         ),
-        meshMover.pMesh(),
+        meshMover().pMesh(),
         dimensionedScalar("medialRatio", dimless, 0.0)
     );
 
@@ -2871,7 +2887,7 @@ void Foam::autoLayerDriver::addLayers
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
         {
-            pointField oldPatchPos(pp.localPoints());
+            pointField oldPatchPos(pp().localPoints());
 
             //// Laplacian displacement shrinking.
             //shrinkMeshDistance
@@ -2886,7 +2902,7 @@ void Foam::autoLayerDriver::addLayers
             // Medial axis based shrinking
             shrinkMeshMedialDistance
             (
-                meshMover,
+                meshMover(),
                 meshQualityDict,
 
                 layerParams.nSmoothThickness(),
@@ -2908,7 +2924,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
@@ -2916,7 +2932,7 @@ void Foam::autoLayerDriver::addLayers
         faceSet dummySet(mesh, "wrongPatchFaces", 0);
         truncateDisplacement
         (
-            meshMover,
+            meshMover(),
             minThickness,
             dummySet,
             patchDisp,
@@ -2931,7 +2947,7 @@ void Foam::autoLayerDriver::addLayers
             dumpDisplacement
             (
                 mesh.time().path()/"layer",
-                pp,
+                pp(),
                 patchDisp,
                 extrudeStatus
             );
@@ -2955,11 +2971,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(),
@@ -2998,7 +3014,7 @@ void Foam::autoLayerDriver::addLayers
         addLayer.setRefinement
         (
             invExpansionRatio,
-            pp,
+            pp(),
             labelList(0),       // exposed patchIDs, not used for adding layers
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point
@@ -3050,8 +3066,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.
@@ -3098,7 +3114,7 @@ void Foam::autoLayerDriver::addLayers
         (
             addLayer,
             meshQualityDict,
-            pp,
+            pp(),
             newMesh,
 
             patchDisp,
@@ -3107,7 +3123,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;
 
@@ -3117,8 +3133,8 @@ void Foam::autoLayerDriver::addLayers
         }
 
         // Reset mesh points and start again
-        meshMover.movePoints(oldPoints);
-        meshMover.correct();
+        meshMover().movePoints(oldPoints);
+        meshMover().correct();
 
         Info<< endl;
     }
@@ -3173,6 +3189,7 @@ void Foam::autoLayerDriver::addLayers
         (
             false,
             false,
+            scalarField(mesh.nCells(), 1.0),
             decomposer,
             distributor
         );
@@ -3211,7 +3228,7 @@ void Foam::autoLayerDriver::doLayers
     fvMeshDistribute& distributor
 )
 {
-    fvMesh& mesh = meshRefiner_.mesh();
+    const fvMesh& mesh = meshRefiner_.mesh();
 
     Info<< nl
         << "Shrinking and layer addition phase" << nl
@@ -3245,59 +3262,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;
+        Info<< "Detected " << nInitErrors << " illegal faces"
+            << " (concave, zero area or negative cell pyramid volume)"
+            << endl;
 
+
+        // Balance
+        if (Pstream::parRun())
         {
-            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();
-
-            // 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)"
+            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 1256957006..512766c539 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -523,8 +523,8 @@ public:
             (
                 const layerParameters& layerParams,
                 const dictionary& motionDict,
+                const labelList& patchIDs,
                 const label nAllowableErrors,
-                motionSmoother& meshMover,
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor
             );
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index 3ea198e6c4..41fe385314 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -35,7 +35,6 @@ License
 #include "refinementSurfaces.H"
 #include "shellSurfaces.H"
 #include "mapDistributePolyMesh.H"
-#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -210,7 +209,8 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
                 "feature refinement iteration " + name(iter),
                 decomposer_,
                 distributor_,
-                cellsToRefine
+                cellsToRefine,
+                refineParams.maxLoadUnbalance()
             );
         }
     }
@@ -310,7 +310,8 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
             "surface refinement iteration " + name(iter),
             decomposer_,
             distributor_,
-            cellsToRefine
+            cellsToRefine,
+            refineParams.maxLoadUnbalance()
         );
     }
     return iter;
@@ -492,7 +493,8 @@ Foam::label Foam::autoRefineDriver::shellRefine
             "shell refinement iteration " + name(iter),
             decomposer_,
             distributor_,
-            cellsToRefine
+            cellsToRefine,
+            refineParams.maxLoadUnbalance()
         );
     }
     meshRefiner_.userFaceData().clear();
@@ -776,11 +778,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/refinementParameters/refinementParameters.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
index ce3b82bfc6..8c79ab9c6c 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
@@ -43,7 +43,8 @@ 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")),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
 {}
 
 
@@ -53,7 +54,8 @@ 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"))),
+    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 1e12dd30a9..ae09e3fc0e 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
@@ -73,6 +73,9 @@ class refinementParameters
         //- Areas to keep
         const pointField keepPoints_;
 
+        //- Allowed load unbalance
+        scalar maxLoadUnbalance_;
+
 
     // Private Member Functions
 
@@ -134,6 +137,12 @@ public:
                 return keepPoints_;
             }
 
+            //- 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 7e90beae6f..1064aa3682 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -552,13 +552,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 +576,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 +585,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 +930,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 +972,17 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 
     Map<label> globalToLocalRegion;
     pointField localPoints;
+    scalarField localWeights;
     calcLocalRegions
     (
         globalCells,
         globalRegion,
         coupledRegionToMaster,
+        cellWeights,
 
         globalToLocalRegion,
-        localPoints
+        localPoints,
+        localWeights
     );
 
 
@@ -984,7 +994,7 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 
     if (isA<geomDecomp>(decomposer))
     {
-        regionDistribution = decomposer.decompose(localPoints);
+        regionDistribution = decomposer.decompose(localPoints, localWeights);
     }
     else
     {
@@ -998,7 +1008,12 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
             regionRegions
         );
 
-        regionDistribution = decomposer.decompose(regionRegions, localPoints);
+        regionDistribution = decomposer.decompose
+        (
+            regionRegions,
+            localPoints,
+            localWeights
+        );
     }
 
 
@@ -1058,6 +1073,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 (
     const bool keepZoneFaces,
     const bool keepBaffles,
+    const scalarField& cellWeights,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor
 )
@@ -1151,6 +1167,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 
                 distribution = decomposeCombineRegions
                 (
+                    cellWeights,
                     blockedFace,
                     couples,
                     decomposer
@@ -1170,13 +1187,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)
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 93ce3ceeb8..5b603ddbac 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,8 @@ public:
                 const string& msg,
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor,
-                const labelList& cellsToRefine
+                const labelList& cellsToRefine,
+                const scalar maxLoadUnbalance
             );
 
 
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index e35607a598..e8c1d43448 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  * * * * * * * * * * * //
 
@@ -1244,40 +1245,122 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
 }
 
 
-// Do refinement of consistent set of cells followed by truncation and
-// load balancing.
+//// Do refinement of consistent set of cells followed by truncation and
+//// load balancing.
+//Foam::autoPtr<Foam::mapDistributePolyMesh>
+//Foam::meshRefinement::refineAndBalance
+//(
+//    const string& msg,
+//    decompositionMethod& decomposer,
+//    fvMeshDistribute& distributor,
+//    const labelList& cellsToRefine
+//)
+//{
+//    // Do all refinement
+//    refine(cellsToRefine);
+//
+//    if (debug)
+//    {
+//        Pout<< "Writing refined but unbalanced " << 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;
+//    printMeshInfo(debug, "After refinement " + msg);
+//
+//
+//    // Load balancing
+//    // ~~~~~~~~~~~~~~
+//
+//    autoPtr<mapDistributePolyMesh> distMap;
+//
+//    if (Pstream::nProcs() > 1)
+//    {
+//        scalarField cellWeights(mesh_.nCells(), 1);
+//
+//        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::refineAndBalance
 (
     const string& msg,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor,
-    const labelList& cellsToRefine
+    const labelList& initCellsToRefine,
+    const scalar maxLoadUnbalance
 )
 {
-    // Do all refinement
-    refine(cellsToRefine);
+    labelList cellsToRefine(initCellsToRefine);
 
-    if (debug)
-    {
-        Pout<< "Writing refined but unbalanced " << 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;
-    printMeshInfo(debug, "After refinement " + msg);
+    //{
+    //    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
@@ -1287,20 +1370,62 @@ Foam::meshRefinement::refineAndBalance
 
     if (Pstream::nProcs() > 1)
     {
-        distMap = balance
+        // 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
         (
-            false,  //keepZoneFaces
-            false,  //keepBaffles
-            decomposer,
-            distributor
+            mag(1.0-nNewCells/nIdealNewCells),
+            maxOp<scalar>()
         );
 
-        Info<< "Balanced mesh in = "
-            << mesh_.time().cpuTimeIncrement() << " s" << endl;
+        if (unbalance <= maxLoadUnbalance)
+        {
+            Info<< "Skipping balancing since max unbalance " << unbalance
+                << " in = "
+                << mesh_.time().cpuTimeIncrement() << " s" << 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)
         {
             Pout<< "Writing balanced " << msg
@@ -1318,6 +1443,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;
 }