diff --git a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index e9bf2cc5dd..69fdc97015 100644
--- a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -861,7 +861,9 @@ Foam::autoPtr<Foam::mapDistributePolyMesh>
 Foam::backgroundMeshDecomposition::distribute
 (
     volScalarField& cellWeights,
-    List<DynamicList<point> >& cellVertices
+    List<DynamicList<point> >& cellVertices,
+    List<DynamicList<label> >& cellVertexIndices,
+    List<DynamicList<label> >& cellVertexTypes
 )
 {
     if (debug)
@@ -953,19 +955,27 @@ Foam::backgroundMeshDecomposition::distribute
         meshCutter_.updateMesh(map);
 
         {
-            // Map cellVertices
+            // Map cellVertices, cellVertexIndices and cellVertexTypes
 
             meshSearch cellSearch(mesh_);
 
             const labelList& reverseCellMap = map().reverseCellMap();
 
             List<DynamicList<point> > newCellVertices(mesh_.nCells());
+            List<DynamicList<label> > newCellVertexIndices(mesh_.nCells());
+            List<DynamicList<label> > newCellVertexTypes(mesh_.nCells());
 
             forAll(cellVertices, oldCellI)
             {
                 DynamicList<point>& oldCellVertices =
                     cellVertices[oldCellI];
 
+                DynamicList<label>& oldCellVertexIndices =
+                    cellVertexIndices[oldCellI];
+
+                DynamicList<label>& oldCellVertexTypes =
+                    cellVertexTypes[oldCellI];
+
                 if (findIndex(newCellsToRefine, oldCellI) >= 0)
                 {
                     // This old cell was refined so the cell for the vertices
@@ -983,10 +993,24 @@ Foam::backgroundMeshDecomposition::distribute
                                 << v << " "
                                 << oldCellI
                                 << newCellI
-                                << endl;
+                                << " find nearest cellI ";
+
+                            newCellI = cellSearch.findNearestCell(v);
+
+                            Pout<< newCellI << endl;
                         }
 
                         newCellVertices[newCellI].append(v);
+
+                        newCellVertexIndices[newCellI].append
+                        (
+                            oldCellVertexIndices[oPI]
+                        );
+
+                        newCellVertexTypes[newCellI].append
+                        (
+                            oldCellVertexTypes[oPI]
+                        );
                     }
                 }
                 else
@@ -995,14 +1019,27 @@ Foam::backgroundMeshDecomposition::distribute
 
                     forAll(oldCellVertices, oPI)
                     {
-                        const point& v = oldCellVertices[oPI];
+                        newCellVertices[newCellI].append
+                        (
+                            oldCellVertices[oPI]
+                        );
 
-                        newCellVertices[newCellI].append(v);
+                        newCellVertexIndices[newCellI].append
+                        (
+                            oldCellVertexIndices[oPI]
+                        );
+
+                        newCellVertexTypes[newCellI].append
+                        (
+                            oldCellVertexTypes[oPI]
+                        );
                     }
                 }
             }
 
             cellVertices.transfer(newCellVertices);
+            cellVertexIndices.transfer(newCellVertexIndices);
+            cellVertexTypes.transfer(newCellVertexTypes);
         }
 
         if (debug)
@@ -1062,12 +1099,15 @@ Foam::backgroundMeshDecomposition::distribute
     }
 
     mapDist().distributeCellData(cellVertices);
+    mapDist().distributeCellData(cellVertexIndices);
+    mapDist().distributeCellData(cellVertexTypes);
 
     buildPatchAndTree();
 
     return mapDist;
 }
 
+
 Foam::autoPtr<Foam::mapDistribute>
 Foam::backgroundMeshDecomposition::distributePoints
 (
diff --git a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
index a9462bf515..51aef7ae3d 100644
--- a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
+++ b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
@@ -211,7 +211,9 @@ public:
         autoPtr<mapDistributePolyMesh> distribute
         (
             volScalarField& cellWeights,
-            List<DynamicList<point> >& cellVertices
+            List<DynamicList<point> >& cellVertices,
+            List<DynamicList<label> >& cellVertexIndices,
+            List<DynamicList<label> >& cellVertexTypes
         );
 
         //- Distribute supplied the points to the appropriate processor
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 1a687c266b..9993162a93 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -357,11 +357,20 @@ void Foam::conformalVoronoiMesh::insertPoints
         // creation of points and indices is done assuming that it will be
         // relative to the instantaneous number_of_vertices() at insertion.
 
+        label type = types[pI];
+
+        if (type > Vb::ptFarPoint)
+        {
+            // This is a member of a point pair, don't use the type directly
+
+            type += number_of_vertices();
+        }
+
         insertPoint
         (
             pts[pI],
             indices[pI] + number_of_vertices(),
-            types[pI] + number_of_vertices()
+            type
         );
     }
 }
@@ -726,30 +735,14 @@ void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
 
 void Foam::conformalVoronoiMesh::insertFeaturePoints()
 {
-    Info<< nl << "Creating bounding points" << endl;
-
-    pointField farPts = geometryToConformTo_.globalBounds().points();
-
-    // Shift corners of bounds relative to origin
-    farPts -= geometryToConformTo_.globalBounds().midpoint();
-
-    // Scale the box up
-    farPts *= 2.0;
-
-    // Shift corners of bounds back to be relative to midpoint
-    farPts += geometryToConformTo_.globalBounds().midpoint();
-
-    forAll(farPts, fPI)
-    {
-        insertPoint(farPts[fPI], Vb::ptFarPoint);
-    }
-
     Info<< nl << "Conforming to feature points" << endl;
 
     DynamicList<Foam::point> pts;
     DynamicList<label> indices;
     DynamicList<label> types;
 
+    label preFeaturePointSize = number_of_vertices();
+
     createConvexFeaturePoints(pts, indices, types);
 
     createConcaveFeaturePoints(pts, indices, types);
@@ -763,7 +756,7 @@ void Foam::conformalVoronoiMesh::insertFeaturePoints()
         writePoints("featureVertices.obj", pts);
     }
 
-    label nFeatureVertices = number_of_vertices() - 8;
+    label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
 
     if (Pstream::parRun())
     {
@@ -772,30 +765,31 @@ void Foam::conformalVoronoiMesh::insertFeaturePoints()
 
     if (nFeatureVertices > 0)
     {
-        Info<< "    Inserted " << nFeatureVertices << " vertices" << endl;
+        Info<< "    Inserted " << nFeatureVertices
+            << " feature vertices" << endl;
     }
 
     Info<< "SORT OUT FEATURE POINT DISTRIBUTION AND STORAGE" << endl;
 
-    featureVertices_.setSize(number_of_vertices());
+    // featureVertices_.setSize(number_of_vertices());
 
-    label featPtI = 0;
+    // label featPtI = 0;
 
-    for
-    (
-        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        vit++
-    )
-    {
-        featureVertices_[featPtI] = Vb(vit->point());
+    // for
+    // (
+    //     Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+    //     vit != finite_vertices_end();
+    //     vit++
+    // )
+    // {
+    //     featureVertices_[featPtI] = Vb(vit->point());
 
-        featureVertices_[featPtI].index() = vit->index();
+    //     featureVertices_[featPtI].index() = vit->index();
 
-        featureVertices_[featPtI].type() = vit->type();
+    //     featureVertices_[featPtI].type() = vit->type();
 
-        featPtI++;
-    }
+    //     featPtI++;
+    // }
 
     constructFeaturePointLocations();
 }
@@ -1124,6 +1118,34 @@ bool Foam::conformalVoronoiMesh::nearFeaturePt(const Foam::point& pt) const
 }
 
 
+void Foam::conformalVoronoiMesh::reset()
+{
+    this->clear();
+
+    insertBoundingPoints();
+}
+
+
+void Foam::conformalVoronoiMesh::insertBoundingPoints()
+{
+    pointField farPts = geometryToConformTo_.globalBounds().points();
+
+    // Shift corners of bounds relative to origin
+    farPts -= geometryToConformTo_.globalBounds().midpoint();
+
+    // Scale the box up
+    farPts *= 2.0;
+
+    // Shift corners of bounds back to be relative to midpoint
+    farPts += geometryToConformTo_.globalBounds().midpoint();
+
+    forAll(farPts, fPI)
+    {
+        insertPoint(farPts[fPI], Vb::ptFarPoint);
+    }
+}
+
+
 void Foam::conformalVoronoiMesh::insertInitialPoints()
 {
     startOfInternalPoints_ = number_of_vertices();
@@ -1154,33 +1176,45 @@ void Foam::conformalVoronoiMesh::insertInitialPoints()
 }
 
 
-void Foam::conformalVoronoiMesh::distribute()
+bool Foam::conformalVoronoiMesh::distributeBackground()
 {
     if (!Pstream::parRun())
     {
-        return;
+        return false;
     }
 
     Info<< nl << "Redistributing points" << endl;
 
     timeCheck("Before distribute");
 
+    label iteration = 0;
+
     while (true)
     {
-        scalar globalNumberOfVertices = returnReduce
+        label nRealVertices = 0;
+
+        for
         (
-            label(number_of_vertices()),
+            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+            vit != finite_vertices_end();
+            vit++
+        )
+        {
+            if (vit->real())
+            {
+                nRealVertices++;
+            }
+        }
+
+        scalar globalNRealVertices = returnReduce
+        (
+            nRealVertices,
             sumOp<label>()
         );
 
         scalar unbalance = returnReduce
         (
-            mag
-            (
-                1.0
-              - label(number_of_vertices())
-               /(globalNumberOfVertices/Pstream::nProcs())
-            ),
+            mag(1.0 - nRealVertices/(globalNRealVertices/Pstream::nProcs())),
             maxOp<scalar>()
         );
 
@@ -1188,15 +1222,17 @@ void Foam::conformalVoronoiMesh::distribute()
 
         if (unbalance <= cvMeshControls().maxLoadUnbalance())
         {
-            break;
+            // If this is the first iteration, return false, if it was a
+            // subsequent one, return true;
+            return iteration != 0;
         }
 
         Info<< "    Total number of vertices before redistribution "
-            << globalNumberOfVertices
+            << globalNRealVertices
             << endl;
 
-        // Pout<< "    number_of_vertices() before distribution "
-        //     << label(number_of_vertices()) << endl;
+        // Pout<< "    Real vertices before distribution "
+        //     << nRealVertices << endl;
 
         const fvMesh& bMesh = decomposition_().mesh();
 
@@ -1218,6 +1254,8 @@ void Foam::conformalVoronoiMesh::distribute()
         meshSearch cellSearch(bMesh);
 
         List<DynamicList<Foam::point> > cellVertices(bMesh.nCells());
+        List<DynamicList<label> > cellVertexIndices(bMesh.nCells());
+        List<DynamicList<label> > cellVertexTypes(bMesh.nCells());
 
         for
         (
@@ -1238,12 +1276,17 @@ void Foam::conformalVoronoiMesh::distribute()
                         << vit->type() << " "
                         << vit->index() << " "
                         << v << " "
-                        << cellI << endl;
-                }
-                else
-                {
-                    cellVertices[cellI].append(topoint(vit->point()));
+                        << cellI
+                        << " find nearest cellI ";
+
+                    cellI = cellSearch.findNearestCell(v);
+
+                    Pout<< cellI << endl;
                 }
+
+                cellVertices[cellI].append(topoint(vit->point()));
+                cellVertexIndices[cellI].append(vit->index());
+                cellVertexTypes[cellI].append(vit->type());
             }
         }
 
@@ -1262,11 +1305,13 @@ void Foam::conformalVoronoiMesh::distribute()
         autoPtr<mapDistributePolyMesh> mapDist = decomposition_().distribute
         (
             cellWeights,
-            cellVertices
+            cellVertices,
+            cellVertexIndices,
+            cellVertexTypes
         );
 
         // Remove the entire tessellation
-        this->clear();
+        reset();
 
         Info<< "NEED TO MAP FEATURE POINTS" << endl;
         // reinsertFeaturePoints();
@@ -1277,26 +1322,43 @@ void Foam::conformalVoronoiMesh::distribute()
 
         Info<< nl << "    Inserting distributed tessellation" << endl;
 
-        std::list<Point> pointsToInsert;
+        DynamicList<Foam::point> pointsToInsert;
+        DynamicList<label> indices;
+        DynamicList<label> types;
 
         forAll(cellVertices, cI)
         {
             forAll(cellVertices[cI], cVPI)
             {
-                pointsToInsert.push_back(toPoint(cellVertices[cI][cVPI]));
+                pointsToInsert.append(cellVertices[cI][cVPI]);
+                indices.append(0);
+
+                label type = cellVertexTypes[cI][cVPI];
+
+                if (type > Vb::ptFarPoint)
+                {
+                    // This is a member of a point pair, don't use the type
+                    // directly, make type relative to the index in preparation
+                    // for insertion.
+
+                    type -= cellVertexIndices[cI][cVPI];
+                }
+
+                types.append(type);
             }
         }
 
         // Assume that the distribution made the correct decision for which
         // processor each point should be on, so give distribute = false
-        insertPoints(pointsToInsert, false);
+        insertPoints(pointsToInsert, indices, types, false);
 
+        // Adjust by 8 because of insertion of bounding points by reset
         Info<< "    Total number of vertices after redistribution "
-            << returnReduce(label(number_of_vertices()), sumOp<label>())
+            << returnReduce(label(number_of_vertices() - 8), sumOp<label>())
             << endl;
 
-        // Pout<< "    number_of_vertices() after distribution "
-        //     << label(number_of_vertices()) << endl;
+        // Pout<< "    Real vertices after distribution "
+        //     << label(number_of_vertices() - 8) << endl;
 
         if(cvMeshControls().objOutput())
         {
@@ -1304,7 +1366,11 @@ void Foam::conformalVoronoiMesh::distribute()
         }
 
         timeCheck("After distribute");
+
+        iteration++;
     }
+
+    return true;
 }
 
 
@@ -1726,22 +1792,34 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         );
     }
 
-    Info<< "INSERTFEATUREPOINTS MOVED" << endl;
-    // insertFeaturePoints();
-
     if (cvMeshControls().objOutput())
     {
         geometryToConformTo_.writeFeatureObj("cvMesh");
     }
 
+    insertBoundingPoints();
+
     insertInitialPoints();
 
-    distribute();
+    // Improve the guess that the backgroundMeshDecomposition makes with the
+    // initial positions.  Use before building the surface conformation to
+    // better balance the surface conformation load.
+    distributeBackground();
 
     insertFeaturePoints();
 
     buildSurfaceConformation(rmCoarse);
 
+    // The introduction of the surface conformation may have distorted the
+    // balance of vertices, distribute if necessary.
+    if (distributeBackground())
+    {
+        // distributeBackground has destroyed all referred vertices, so the
+        // parallel interface needs to be rebuilt.
+
+        buildParallelInterface("rebuild");
+    }
+
     if(cvMeshControls().objOutput())
     {
         writePoints("allInitialPoints.obj", false);
@@ -2102,7 +2180,7 @@ void Foam::conformalVoronoiMesh::move()
     }
 
     // Remove the entire tessellation
-    this->clear();
+    reset();
 
     reinsertFeaturePoints();
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 0e4578bfcc..ffc87006d6 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -428,13 +428,22 @@ private:
         //- Check if a location is in exclusion range around a feature point
         bool nearFeaturePt(const Foam::point& pt) const;
 
+        //- Clear the entire tesselation and insert bounding points
+        void reset();
+
+        //- Insert far points in a large bounding box to avoid dual edges
+        //  spanning huge distances
+        void insertBoundingPoints();
+
         //- Insert the initial points into the triangulation, based on the
         //  initialPointsMethod
         void insertInitialPoints();
 
         //- In parallel redistribute the backgroundMeshDecomposition and
-        //- vertices to balance the number of vertices on each processor
-        void distribute();
+        //  vertices to balance the number of vertices on each processor.
+        //  Returns true if the background mesh changes as this removes all
+        //  referred vertices, so the parallel interface may need rebuilt.
+        bool distributeBackground();
 
         //- Store data for sizeAndAlignmentLocations_, storedSizes_ and
         //  storedAlignments_ and initialise the sizeAndAlignmentTreePtr_
@@ -509,6 +518,13 @@ private:
 
         //- Build the parallelInterfaces of the mesh
         void buildParallelInterface
+        (
+            const word& outputName
+        );
+
+        //- Build the parallelInterfaces of the mesh, supply hash sets
+        //  externally to allow updates
+        void buildParallelInterface
         (
             List<labelHashSet>& referralVertices,
             List<labelHashSet>& receivedVertices,
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index af791a589b..1b88564fa6 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -551,6 +551,24 @@ bool Foam::conformalVoronoiMesh::clipLineToProc
 }
 
 
+void Foam::conformalVoronoiMesh::buildParallelInterface
+(
+    const word& outputName
+)
+{
+    List<labelHashSet> referralVertices(Pstream::nProcs());
+    List<labelHashSet> receivedVertices(Pstream::nProcs());
+
+    buildParallelInterface
+    (
+        referralVertices,
+        receivedVertices,
+        true,
+        outputName
+    );
+}
+
+
 void Foam::conformalVoronoiMesh::buildParallelInterface
 (
     List<labelHashSet>& referralVertices,