diff --git a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index f4ff9245e8..41a37cbb17 100644
--- a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "backgroundMeshDecomposition.H"
+#include "meshSearch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -58,13 +59,6 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
 
     decompositionMethod& decomposer = decomposerPtr_();
 
-    hexRef8 meshCutter
-    (
-        mesh_,
-        labelList(mesh_.nCells(), 0),
-        labelList(mesh_.nPoints(), 0)
-    );
-
     // For each cell in the mesh has it been determined if it is fully
     // inside, outside, or overlaps the surface
     labelList volumeStatus
@@ -110,7 +104,6 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
             {
                 labelList refCells = selectRefinementCells
                 (
-                    meshCutter,
                     volumeStatus,
                     cellWeights
                 );
@@ -118,7 +111,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                 // Maintain 2:1 ratio
                 labelList newCellsToRefine
                 (
-                    meshCutter.consistentRefinement
+                    meshCutter_.consistentRefinement
                     (
                         refCells,
                         true                  // extend set
@@ -141,12 +134,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                     );
                 }
 
-
-                if
-                (
-                    returnReduce(newCellsToRefine.size(), sumOp<label>())
-                 == 0
-                )
+                if (returnReduce(newCellsToRefine.size(), sumOp<label>()) == 0)
                 {
                     break;
                 }
@@ -155,7 +143,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                 polyTopoChange meshMod(mesh_);
 
                 // Play refinement commands into mesh changer.
-                meshCutter.setRefinement(newCellsToRefine, meshMod);
+                meshCutter_.setRefinement(newCellsToRefine, meshMod);
 
                 // Create mesh, return map from old to new mesh.
                 autoPtr<mapPolyMesh> map = meshMod.changeMesh
@@ -163,7 +151,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                     mesh_,
                     false,  // inflate
                     true,   // syncParallel
-                    true,   // orderCells (to reduce cell motion in scotch)
+                    true,   // orderCells (to reduce cell transfers)
                     false   // orderPoints
                 );
 
@@ -171,7 +159,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                 mesh_.updateMesh(map);
 
                 // Update numbering of cells/vertices.
-                meshCutter.updateMesh(map);
+                meshCutter_.updateMesh(map);
 
                 {
                     // Map volumeStatus
@@ -187,12 +175,11 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                         if (oldCellI == -1)
                         {
                             newVolumeStatus[newCellI] =
-                            searchableSurface::UNKNOWN;;
+                                searchableSurface::UNKNOWN;
                         }
                         else
                         {
-                            newVolumeStatus[newCellI] =
-                            volumeStatus[oldCellI];
+                            newVolumeStatus[newCellI] = volumeStatus[oldCellI];
                         }
                     }
 
@@ -201,7 +188,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
 
                 if (debug)
                 {
-                    Info<< "Refined from "
+                    Info<< "    Background mesh refined from "
                         << returnReduce(map().nOldCells(), sumOp<label>())
                         << " to " << mesh_.globalData().nTotalCells()
                         << " cells." << endl;
@@ -238,6 +225,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                 }
             }
 
+            // Hard code switch for this stage for testing
             bool removeOutsideCells = false;
 
             if (removeOutsideCells)
@@ -277,7 +265,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                     mesh_,
                     false,  // inflate
                     true,   // syncParallel
-                    true,   // orderCells (to reduce cell motion in scotch)
+                    true,   // orderCells (to reduce cell transfers)
                     false   // orderPoints
                 );
 
@@ -285,7 +273,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                 mesh_.updateMesh(map);
 
                 // Update numbering of cells/vertices.
-                meshCutter.updateMesh(map);
+                meshCutter_.updateMesh(map);
                 cellRemover.updateMesh(map);
 
                 {
@@ -323,7 +311,8 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
             if (debug)
             {
                 const_cast<Time&>(mesh_.time())++;
-                meshCutter.write();
+                Info<< "Time " << mesh_.time().timeName() << endl;
+                meshCutter_.write();
                 mesh_.write();
                 cellWeights.write();
             }
@@ -340,7 +329,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
             autoPtr<mapDistributePolyMesh> mapDist =
                 distributor.distribute(newDecomp);
 
-            meshCutter.distribute(mapDist);
+            meshCutter_.distribute(mapDist);
 
             mapDist().distributeCellData(volumeStatus);
 
@@ -349,7 +338,8 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                 printMeshData(mesh_);
 
                 const_cast<Time&>(mesh_.time())++;
-                meshCutter.write();
+                Info<< "Time " << mesh_.time().timeName() << endl;
+                meshCutter_.write();
                 mesh_.write();
                 cellWeights.write();
             }
@@ -359,6 +349,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
     if (debug)
     {
         const_cast<Time&>(mesh_.time())++;
+        Info<< "Time " << mesh_.time().timeName() << endl;
         cellWeights.write();
         mesh_.write();
     }
@@ -424,7 +415,6 @@ void Foam::backgroundMeshDecomposition::printMeshData
 
 bool Foam::backgroundMeshDecomposition::refineCell
 (
-    const polyMesh& mesh,
     label cellI,
     label volType,
     scalar& weightEstimate
@@ -437,10 +427,10 @@ bool Foam::backgroundMeshDecomposition::refineCell
 
     treeBoundBox cellBb
     (
-        mesh.cells()[cellI].points
+        mesh_.cells()[cellI].points
         (
-            mesh.faces(),
-            mesh.points()
+            mesh_.faces(),
+            mesh_.points()
         )
     );
 
@@ -549,21 +539,18 @@ bool Foam::backgroundMeshDecomposition::refineCell
 
 Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
 (
-    const hexRef8& meshCutter,
     labelList& volumeStatus,
     volScalarField& cellWeights
 ) const
 {
     labelHashSet cellsToRefine;
 
-    const polyMesh& mesh = meshCutter.mesh();
-
     // Determine/update the status of each cell
     forAll(volumeStatus, cellI)
     {
         if (volumeStatus[cellI] == searchableSurface::MIXED)
         {
-            if (meshCutter.cellLevel()[cellI] < minLevels_)
+            if (meshCutter_.cellLevel()[cellI] < minLevels_)
             {
                 cellsToRefine.insert(cellI);
             }
@@ -575,7 +562,6 @@ Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
             (
                 refineCell
                 (
-                    mesh,
                     cellI,
                     volumeStatus[cellI],
                     cellWeights.internalField()[cellI]
@@ -693,6 +679,12 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
             IOobject::MUST_READ
         )
     ),
+    meshCutter_
+    (
+        mesh_,
+        labelList(mesh_.nCells(), 0),
+        labelList(mesh_.nPoints(), 0)
+    ),
     boundaryFacesPtr_(),
     bFTreePtr_(),
     decomposeDict_
@@ -714,7 +706,8 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
         coeffsDict_.lookupOrDefault<scalar>("minCellSizeLimit", 0.0)
     ),
     minLevels_(readLabel(coeffsDict_.lookup("minLevels"))),
-    volRes_(readLabel(coeffsDict_.lookup("sampleResolution")))
+    volRes_(readLabel(coeffsDict_.lookup("sampleResolution"))),
+    maxCellWeightCoeff_(readScalar(coeffsDict_.lookup("maxCellWeightCoeff")))
 {
     if (!Pstream::parRun())
     {
@@ -761,16 +754,176 @@ Foam::backgroundMeshDecomposition::~backgroundMeshDecomposition()
 Foam::autoPtr<Foam::mapDistributePolyMesh>
 Foam::backgroundMeshDecomposition::distribute
 (
-    volScalarField& cellWeights
+    volScalarField& cellWeights,
+    List<DynamicList<point> >& cellVertices
 )
 {
     if (debug)
     {
         const_cast<Time&>(mesh_.time())++;
+        Info<< "Time " << mesh_.time().timeName() << endl;
         cellWeights.write();
         mesh_.write();
     }
 
+    while (true)
+    {
+        // Refine large cells if necessary
+
+        scalar cellWeightLimit =
+            maxCellWeightCoeff_
+           *sum(cellWeights).value()
+           /mesh_.globalData().nTotalCells();
+
+        if (debug)
+        {
+            Info<< "    cellWeightLimit " << cellWeightLimit << endl;
+
+            Pout<< "    sum(cellWeights) "
+                << sum(cellWeights.internalField())
+                << " max(cellWeights) "
+                << max(cellWeights.internalField())
+                << endl;
+        }
+
+        labelHashSet cellsToRefine;
+
+        forAll(cellWeights, cWI)
+        {
+            if (cellWeights.internalField()[cWI] > cellWeightLimit)
+            {
+                cellsToRefine.insert(cWI);
+            }
+        }
+
+        if (returnReduce(cellsToRefine.size(), sumOp<label>()) == 0)
+        {
+            break;
+        }
+
+        // Maintain 2:1 ratio
+        labelList newCellsToRefine
+        (
+            meshCutter_.consistentRefinement
+            (
+                cellsToRefine.toc(),
+                true                  // extend set
+            )
+        );
+
+        if (debug && !cellsToRefine.empty())
+        {
+            Pout<< "    cellWeights too large in " << cellsToRefine.size()
+                << " cells" << endl;
+        }
+
+        forAll(newCellsToRefine, nCTRI)
+        {
+            label cellI = newCellsToRefine[nCTRI];
+
+            cellWeights.internalField()[cellI] /= 8.0;
+        }
+
+        // Mesh changing engine.
+        polyTopoChange meshMod(mesh_);
+
+        // Play refinement commands into mesh changer.
+        meshCutter_.setRefinement(newCellsToRefine, meshMod);
+
+        // Create mesh, return map from old to new mesh.
+        autoPtr<mapPolyMesh> map = meshMod.changeMesh
+        (
+            mesh_,
+            false,  // inflate
+            true,   // syncParallel
+            true,   // orderCells (to reduce cell motion)
+            false   // orderPoints
+        );
+
+        // Update fields
+        mesh_.updateMesh(map);
+
+        // Update numbering of cells/vertices.
+        meshCutter_.updateMesh(map);
+
+        {
+            // Map cellVertices
+
+            meshSearch cellSearch(mesh_);
+
+            const labelList& reverseCellMap = map().reverseCellMap();
+
+            List<DynamicList<point> > newCellVertices(mesh_.nCells());
+
+            forAll(cellVertices, oldCellI)
+            {
+                DynamicList<point>& oldCellVertices =
+                    cellVertices[oldCellI];
+
+                if (findIndex(newCellsToRefine, oldCellI) >= 0)
+                {
+                    // This old cell was refined so the cell for the vertices
+                    // in the new mesh needs to be searched for.
+
+                    forAll (oldCellVertices, oPI)
+                    {
+                        const point& v = oldCellVertices[oPI];
+
+                        label newCellI = cellSearch.findCell(v);
+
+                        if (newCellI == -1)
+                        {
+                            Pout<< "findCell backgroundMeshDecomposition "
+                                << v << " "
+                                << oldCellI
+                                << newCellI
+                                << endl;
+                        }
+
+                        newCellVertices[newCellI].append(v);
+                    }
+                }
+                else
+                {
+                    label newCellI = reverseCellMap[oldCellI];
+
+                    forAll(oldCellVertices, oPI)
+                    {
+                        const point& v = oldCellVertices[oPI];
+
+                        newCellVertices[newCellI].append(v);
+                    }
+                }
+            }
+
+            cellVertices.transfer(newCellVertices);
+        }
+
+        if (debug)
+        {
+            Info<< "    Background mesh refined from "
+                << returnReduce(map().nOldCells(), sumOp<label>())
+                << " to " << mesh_.globalData().nTotalCells()
+                << " cells." << endl;
+
+            const_cast<Time&>(mesh_.time())++;
+            Info<< "Time " << mesh_.time().timeName() << endl;
+            cellWeights.write();
+            mesh_.write();
+        }
+    }
+
+    if (debug)
+    {
+        printMeshData(mesh_);
+
+        Pout<< "    Pre distribute sum(cellWeights) "
+            << sum(cellWeights.internalField())
+            << " max(cellWeights) "
+            << max(cellWeights.internalField())
+            << endl;
+    }
+
     labelList newDecomp = decomposerPtr_().decompose
     (
         mesh_,
@@ -778,20 +931,32 @@ Foam::backgroundMeshDecomposition::distribute
         cellWeights
     );
 
+    Info<< "    Redistributing background mesh cells" << endl;
+
     fvMeshDistribute distributor(mesh_, mergeDist_);
 
-    autoPtr<mapDistributePolyMesh> mapDist =
-        distributor.distribute(newDecomp);
+    autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute(newDecomp);
+
+    meshCutter_.distribute(mapDist);
 
     if (debug)
     {
         printMeshData(mesh_);
 
+        Pout<< "    Post distribute sum(cellWeights) "
+            << sum(cellWeights.internalField())
+            << " max(cellWeights) "
+            << max(cellWeights.internalField())
+            << endl;
+
         const_cast<Time&>(mesh_.time())++;
+        Info<< "Time " << mesh_.time().timeName() << endl;
         mesh_.write();
         cellWeights.write();
     }
 
+    mapDist().distributeCellData(cellVertices);
+
     buildPatchAndTree();
 
     return mapDist;
diff --git a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
index c6d77b89bd..f6a0b448a6 100644
--- a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
+++ b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
@@ -101,6 +101,9 @@ class backgroundMeshDecomposition
         //  is responsible for.
         fvMesh mesh_;
 
+        //- Refinement object
+        hexRef8 meshCutter_;
+
         //- Patch containing an independent representation of the surface of the
         //  mesh of this processor
         autoPtr<bPatch> boundaryFacesPtr_;
@@ -131,6 +134,10 @@ class backgroundMeshDecomposition
         //  investigate the local cell size
         label volRes_;
 
+        //- Allowed factor above the average cell weight before a background
+        //  cell needs to be split
+        scalar maxCellWeightCoeff_;
+
 
     // Private Member Functions
 
@@ -144,7 +151,6 @@ class backgroundMeshDecomposition
         //  it
         bool refineCell
         (
-            const polyMesh& mesh,
             label cellI,
             label volType,
             scalar& weightEstimate
@@ -154,7 +160,6 @@ class backgroundMeshDecomposition
         //  meshed
         labelList selectRefinementCells
         (
-            const hexRef8& meshCutter,
             labelList& volumeStatus,
             volScalarField& cellWeights
         ) const;
@@ -195,7 +200,8 @@ public:
         //  returning a map to use to redistribute vertices.
         autoPtr<mapDistributePolyMesh> distribute
         (
-            volScalarField& cellWeights
+            volScalarField& cellWeights,
+            List<DynamicList<point> >& cellVertices
         );
 
         //- Is the given position inside the domain of this decomposition
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 08f2f0d596..a1dcbbec6d 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -859,111 +859,150 @@ void Foam::conformalVoronoiMesh::insertInitialPoints()
 
 void Foam::conformalVoronoiMesh::distribute()
 {
+    if (!Pstream::parRun())
+    {
+        return;
+    }
+
     Info<< nl << "Redistributing points" << endl;
 
     timeCheck("Before distribute");
 
-    const fvMesh& bMesh = decomposition_().mesh();
-
-    volScalarField cellWeights
-    (
-        IOobject
-        (
-            "cellWeights",
-            bMesh.time().timeName(),
-            bMesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        bMesh,
-        dimensionedScalar("weight", dimless, 0),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    meshSearch cellSearch(bMesh);
-
-    List<DynamicList<Foam::point> > cellVertices(bMesh.nCells());
-
-    for
-    (
-        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        vit++
-    )
+    while (true)
     {
-        if (vit->real())
-        {
-            Foam::point v = topoint(vit->point());
-
-            label cellI = cellSearch.findCell(v);
-
-            if (cellI == -1)
-            {
-                Pout<< "findCell "
-                    << vit->type() << " "
-                    << vit->index() << " "
-                    << v << " "
-                    << cellI << endl;
-            }
-            else
-            {
-                cellVertices[cellI].append(topoint(vit->point()));
-            }
-        }
-    }
-
-    forAll(cellVertices, cI)
-    {
-        cellWeights.internalField()[cI] = min
+        scalar globalNumberOfVertices = returnReduce
         (
-            max(cellVertices[cI].size(), 1e-3),
-            100
+            label(number_of_vertices()),
+            sumOp<label>()
         );
-    }
 
-    autoPtr<mapDistributePolyMesh> mapDist = decomposition_().distribute
-    (
-        cellWeights
-    );
+        scalar unbalance = returnReduce
+        (
+            mag
+            (
+                1.0
+              - label(number_of_vertices())
+               /(globalNumberOfVertices/Pstream::nProcs())
+            ),
+            maxOp<scalar>()
+        );
 
-    Pout<< "number_of_vertices() before distribution "
-        << label(number_of_vertices()) << endl;
+        Info<< "    Processor unbalance " << unbalance << endl;
 
-    mapDist().distributeCellData(cellVertices);
-
-    // Remove the entire tessellation
-    this->clear();
-
-    // Info<< "NEED TO MAP FEATURE POINTS" << endl;
-    // reinsertFeaturePoints();
-
-    startOfInternalPoints_ = number_of_vertices();
-
-    timeCheck("Distribution performed");
-
-    Info<< nl << "Inserting distributed tessellation" << endl;
-
-    std::vector<Point> pointsToInsert;
-
-    forAll(cellVertices, cI)
-    {
-        forAll(cellVertices[cI], cVPI)
+        if (unbalance <= cvMeshControls().maxLoadUnbalance())
         {
-            pointsToInsert.push_back(toPoint(cellVertices[cI][cVPI]));
+            break;
         }
+
+        Info<< "    Total number of vertices before redistribution "
+            << globalNumberOfVertices
+            << endl;
+
+        // Pout<< "    number_of_vertices() before distribution "
+        //     << label(number_of_vertices()) << endl;
+
+        const fvMesh& bMesh = decomposition_().mesh();
+
+        volScalarField cellWeights
+        (
+            IOobject
+            (
+                "cellWeights",
+                bMesh.time().timeName(),
+                bMesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            bMesh,
+            dimensionedScalar("weight", dimless, 1e-3),
+            zeroGradientFvPatchScalarField::typeName
+        );
+
+        meshSearch cellSearch(bMesh);
+
+        List<DynamicList<Foam::point> > cellVertices(bMesh.nCells());
+
+        for
+        (
+            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+            vit != finite_vertices_end();
+            vit++
+        )
+        {
+            if (vit->real())
+            {
+                Foam::point v = topoint(vit->point());
+
+                label cellI = cellSearch.findCell(v);
+
+                if (cellI == -1)
+                {
+                    Pout<< "findCell conformalVoronoiMesh::distribute findCell "
+                        << vit->type() << " "
+                        << vit->index() << " "
+                        << v << " "
+                        << cellI << endl;
+                }
+                else
+                {
+                    cellVertices[cellI].append(topoint(vit->point()));
+                }
+            }
+        }
+
+        forAll(cellVertices, cI)
+        {
+            cellWeights.internalField()[cI] = max
+            (
+                cellVertices[cI].size(),
+                1e-2       // Small but finite weight for empty cells
+            );
+        }
+
+        autoPtr<mapDistributePolyMesh> mapDist = decomposition_().distribute
+        (
+            cellWeights,
+            cellVertices
+        );
+
+        // Remove the entire tessellation
+        this->clear();
+
+        Info<< "NEED TO MAP FEATURE POINTS" << endl;
+        // reinsertFeaturePoints();
+
+        startOfInternalPoints_ = number_of_vertices();
+
+        timeCheck("Distribution performed");
+
+        Info<< nl << "    Inserting distributed tessellation" << endl;
+
+        std::vector<Point> pointsToInsert;
+
+        forAll(cellVertices, cI)
+        {
+            forAll(cellVertices[cI], cVPI)
+            {
+                pointsToInsert.push_back(toPoint(cellVertices[cI][cVPI]));
+            }
+        }
+
+        insertPoints(pointsToInsert);
+
+        Info<< "    Total number of vertices after redistribution "
+            << returnReduce(label(number_of_vertices()), sumOp<label>())
+            << endl;
+
+        // Pout<< "    number_of_vertices() after distribution "
+        //     << label(number_of_vertices()) << endl;
+
+        if(cvMeshControls().objOutput())
+        {
+            writePoints("distributedPoints.obj", true);
+        }
+
+        timeCheck("After distribute");
     }
-
-    insertPoints(pointsToInsert);
-
-    Pout<< "number_of_vertices() after distribution "
-        << label(number_of_vertices()) << endl;
-
-    if(cvMeshControls().objOutput())
-    {
-        writePoints("distributedPoints.obj", true);
-    }
-
-    timeCheck("After distribute");
 }
 
 
@@ -1377,6 +1416,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         );
     }
 
+    Info<< "CREATEFEATUREPOINTS MOVED" << endl;
     // createFeaturePoints();
 
     if (cvMeshControls().objOutput())