diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index e8cedc49c2cbd5e05256b513c4be1250f8140e4f..958fcf6763ec79332bf98faee94882df8eb0c640 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -240,9 +240,9 @@ void Foam::conformalVoronoiMesh::insertPoints
     {
         label preDistributionSize(points.size());
 
-        DynamicList<Foam::point> transferPoints(points.size()/2);
+        DynamicList<Foam::point> transferPoints;
 
-        DynamicList<Point> pointsOnProcessor(points.size()/2);
+        DynamicList<Point> pointsOnProcessor;
 
         for
         (
@@ -277,12 +277,16 @@ void Foam::conformalVoronoiMesh::insertPoints
             decomposition_().distributePoints(transferPoints)
         );
 
+        const label oldSize = points.size();
+
+        points.setSize(oldSize + transferPoints.size());
+
         forAll(transferPoints, tPI)
         {
-            points.append(toPoint(transferPoints[tPI]));
+            points[tPI + oldSize] = toPoint(transferPoints[tPI]);
         }
 
-        label sizeChange = preDistributionSize - label(points.size());
+        label sizeChange = preDistributionSize - points.size();
 
         // if (mag(sizeChange) > 0)
         // {
@@ -417,7 +421,6 @@ void Foam::conformalVoronoiMesh::insertPoints
 //        );
 //    }
 
-
     rangeInsertWithInfo
     (
         pts.begin(),
@@ -1246,6 +1249,8 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
     // better balance the surface conformation load.
     distributeBackground();
 
+//    conformToSurface();
+
     buildSurfaceConformation(rmCoarse);
 
     // The introduction of the surface conformation may have distorted the
@@ -1313,7 +1318,7 @@ void Foam::conformalVoronoiMesh::move()
         {
             cit->cellIndex() = dualVertI;
 
-            dualVertices[dualVertI] = topoint(dual(cit));
+            dualVertices[dualVertI] = cit->dual();
 
             dualVertI++;
         }
@@ -1511,7 +1516,6 @@ void Foam::conformalVoronoiMesh::move()
                             (
                                 toPoint(0.5*(dVA + dVB))
                             );
-
                         }
                     }
                     else if
@@ -1671,9 +1675,57 @@ void Foam::conformalVoronoiMesh::move()
 
     insertPoints(pointsToInsert);
 
+    // Remove internal points that have been inserted outside the surface.
+    label internalPtIsOutside = 0;
+
+    for
+    (
+        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->internalPoint())
+        {
+            bool inside
+                = geometryToConformTo_.inside(topoint(vit->point()));
+
+            if (!inside)
+            {
+                remove(vit);
+                internalPtIsOutside++;
+            }
+        }
+    }
+
+    Info<< "    " << internalPtIsOutside
+        << " internal points were inserted outside the domain. "
+        << "They have been removed." << endl;
+
+    // Fix points that have not been significantly displaced
+//    for
+//    (
+//        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+//        vit != finite_vertices_end();
+//        ++vit
+//    )
+//    {
+//        if (vit->internalPoint())
+//        {
+//            if
+//            (
+//                mag(displacementAccumulator[vit->index()])
+//              < 0.1*targetCellSize(topoint(vit->point()))
+//            )
+//            {
+//                vit->setVertexFixed();
+//            }
+//        }
+//    }
+
     if (cvMeshControls().objOutput() && runTime_.outputTime())
     {
-        writePoints("points_" + runTime_.timeName() + ".obj", false);
+        writePoints("points_" + runTime_.timeName() + ".obj", true);
     }
 
     timeCheck("Internal points inserted");
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index b590eb31cd67383146216cfbdc4f211a38397137..a514c660a1d96350897478ca826dba7efe5af1d8 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -48,7 +48,6 @@ SourceFiles
 #define CGAL_INEXACT
 
 #include "CGALTriangulation3Ddefs.H"
-#include <CGAL/Spatial_sort_traits_adapter_3.h>
 #include "uint.H"
 #include "ulong.H"
 #include "searchableSurfaces.H"
@@ -595,6 +594,12 @@ private:
             Foam::point& b
         ) const;
 
+        label removeProcessorBoundarySeeds(bool reinsertBoundPts);
+
+        void seedProcessorBoundarySurfaces(bool seedProcessors);
+
+        label numberOfUnusedReferredPoints() const;
+
         //- Build the parallelInterfaces of the mesh
         void buildParallelInterface
         (
@@ -1240,7 +1245,7 @@ public:
                     Traits_for_spatial_sort<Triangulation>()
                 );
 
-                typename Triangulation::Cell_handle hint;
+                typename Triangulation::Vertex_handle hint;
 
                 for
                 (
@@ -1250,16 +1255,9 @@ public:
                     ++p
                 )
                 {
-                    typename Triangulation::Locate_type lt;
-                    typename Triangulation::Cell_handle c;
-                    label li, lj;
-
-                    c = T.locate(*(p->first), lt, li, lj, hint);
-
                     const size_t checkInsertion = T.number_of_vertices();
 
-                    typename Triangulation::Vertex_handle v
-                        = T.insert(*(p->first), lt, c, li, lj);
+                    hint = T.insert(*(p->first), hint);
 
                     if (checkInsertion != T.number_of_vertices() - 1)
                     {
@@ -1278,12 +1276,11 @@ public:
                             // type directly (note that this routine never gets
                             // called for referredPoints so type will never be
                             // -procI
-                            type += T.number_of_vertices() - 1;
+                            type += checkInsertion;
                         }
 
-                        v->index() = indices[oldIndex]
-                                   + T.number_of_vertices() - 1;
-                        v->type() = type;
+                        hint->index() = indices[oldIndex] + checkInsertion;
+                        hint->type() = type;
                     }
                 }
             }
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index 04ef2af045d22226200bce8d1943ec4ed9438696..d103d3c0fcb1ff968452105319a3c304e5a65578 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -1890,7 +1890,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
         {
             cit->cellIndex() = dualVertI;
 
-            pts[dualVertI] = topoint(dual(cit));
+            pts[dualVertI] = cit->dual();
 
             if
             (
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index 847974a2f4ed3c01bdc652fe52baafafcb8bc96c..7f4c7abb6ef2cd34a4513fda2e94a41cf5af0d9b 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -41,6 +41,11 @@ void Foam::conformalVoronoiMesh::conformToSurface()
 {
     reconformationMode reconfMode = reconformationControl();
 
+    if (Pstream::parRun())
+    {
+        seedProcessorBoundarySurfaces(true);
+    }
+
     if (reconfMode == rmNone)
     {
         // Reinsert stored surface conformation
@@ -69,6 +74,16 @@ void Foam::conformalVoronoiMesh::conformToSurface()
         storeSurfaceConformation();
     }
 
+    if (Pstream::parRun())
+    {
+        label nFarPoints = removeProcessorBoundarySeeds(true);
+
+        reduce(nFarPoints, sumOp<label>());
+
+        Info<< "    Removed " << nFarPoints
+            << " far points from the mesh." << endl;
+    }
+
     // reportSurfaceConformationQuality();
 }
 
@@ -140,13 +155,8 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
     buildEdgeLocationTree(existingEdgeLocations);
     buildSurfacePtLocationTree(existingSurfacePtLocations);
 
-
-    // Initialise the edgeLocationTree
-    //buildEdgeLocationTree(edgeLocationTree, existingEdgeLocations);
-
     label initialTotalHits = 0;
 
-
     // Surface protrusion conformation is done in two steps.
     // 1. the dual edges (of all internal vertices) can stretch to
     //    'infinity' so any intersection would be badly behaved. So
@@ -382,7 +392,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
         (
             Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
             vit != finite_vertices_end();
-            vit++
+            ++vit
         )
         {
             // The initial surface conformation has already identified the
@@ -511,14 +521,18 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
 
         timeCheck("Conformation iteration " + name(iterationNo));
 
-        // Update the parallel interface
-        buildParallelInterface
-        (
-            referralVertices,
-            receivedVertices,
-            false,
-            name(iterationNo)
-        );
+        // Only need to update the interface if there are surface/edge hits
+        if (totalHits > 0)
+        {
+            // Update the parallel interface
+            buildParallelInterface
+            (
+                referralVertices,
+                receivedVertices,
+                false,
+                name(iterationNo)
+            );
+        }
 
         iterationNo++;
 
@@ -563,8 +577,8 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
             continue;
         }
 
-        Foam::point dE0 = topoint(dual(fit->first));
-        Foam::point dE1 = topoint(dual(fit->first->neighbor(fit->second)));
+        Foam::point dE0 = fit->first->dual();
+        Foam::point dE1 = fit->first->neighbor(fit->second)->dual();
 
         if (Pstream::parRun())
         {
@@ -631,8 +645,8 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
 
         // Construct the dual edge and search for intersections of the edge
         // with the surface
-        Foam::point dE0 = topoint(dual(fit->first));
-        Foam::point dE1 = topoint(dual(fit->first->neighbor(fit->second)));
+        Foam::point dE0 = fit->first->dual();
+        Foam::point dE1 = fit->first->neighbor(fit->second)->dual();
 
         pointIndexHit infoIntersection;
         label hitSurfaceIntersection = -1;
@@ -801,6 +815,144 @@ void Foam::conformalVoronoiMesh::buildParallelInterface
 }
 
 
+Foam::label Foam::conformalVoronoiMesh::removeProcessorBoundarySeeds
+(
+    bool reinsertBoundPts
+)
+{
+    label nFarPoints = 0;
+
+
+
+    std::list<Vertex_handle> toRemove;
+
+    for
+    (
+        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->farPoint())
+        {
+            //remove(vit);
+            toRemove.push_back(vit);
+            nFarPoints++;
+        }
+    }
+
+    // This function removes the points in the iterator range and then
+    // retriangulates.
+    timeCheck("Start Removing Seeded Points " + name(toRemove.size()));
+
+    remove_cluster(toRemove.begin(), toRemove.end());
+
+
+    // Need to do this to make sure the triangulation is well-behaved
+    if (reinsertBoundPts)
+    {
+        reinsertBoundingPoints();
+    }
+
+
+
+//    DynamicList<Foam::point> toAdd;
+//    DynamicList<label> indices;
+//    DynamicList<label> types;
+//
+//    for
+//    (
+//        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+//        vit != finite_vertices_end();
+//        ++vit
+//    )
+//    {
+//        if (!vit->farPoint())
+//        {
+//            toAdd.append(topoint(vit->point()));
+//            indices.append(vit->index());
+//            types.append(vit->type());
+//
+//            nFarPoints++;
+//        }
+//    }
+//
+//    this->clear();
+//
+//    // Need to do this to make sure the triangulation is well-behaved
+//    if (reinsertBoundPts)
+//    {
+//        reinsertBoundingPoints();
+//    }
+//
+//    forAll(toAdd, pI)
+//    {
+//        insertPoint(toAdd[pI], indices[pI], types[pI]);
+//    }
+
+
+    timeCheck("End Removing Seeded Points");
+
+    return nFarPoints;
+}
+
+
+void Foam::conformalVoronoiMesh::seedProcessorBoundarySurfaces
+(
+    bool seedProcessors
+)
+{
+    //removeProcessorBoundarySeeds(false);
+
+    // Loop over processor patch faces and insert a point in the centre of the
+    // face
+    const fvMesh& mesh = decomposition_().mesh();
+    const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
+
+    DynamicList<Foam::point> pts;
+    DynamicList<label> indices;
+    DynamicList<label> types;
+
+    label nFarPoints = 0;
+
+    const scalar normalDistance = 5.0;
+    const scalar pert = 0.1*(rndGen_.scalar01() - 0.5);
+
+    forAll(bMesh, patchI)
+    {
+        const polyPatch& patch = bMesh[patchI];
+
+        if (!seedProcessors && isA<processorPolyPatch>(patch))
+        {
+            continue;
+        }
+
+        forAll(patch, faceI)
+        {
+            if (faceI % 1 == 0)
+            {
+                const face& f = patch[faceI];
+
+                pts.append
+                (
+                    f.centre(mesh.points())
+                  + pert*normalDistance*f.normal(mesh.points())
+                );
+                indices.append(nFarPoints++);
+                types.append(Vb::vtFar);
+            }
+        }
+    }
+
+    insertPoints(pts, indices, types, false);
+
+    reduce(nFarPoints, sumOp<label>());
+
+    Info<< "    Inserted " << nFarPoints
+        << " far points into the mesh." << endl;
+}
+
+
 void Foam::conformalVoronoiMesh::buildParallelInterface
 (
     List<labelHashSet>& referralVertices,
@@ -835,22 +987,88 @@ void Foam::conformalVoronoiMesh::buildParallelInterface
         timeCheck("After buildParallelInterfaceAll");
     }
 
-    if (initialEdgeReferral)
-    {
-        // Used as an initial pass to localise the vertex referring - find
-        // vertices whose dual edges pierce nearby processor volumes and refer
-        // them to establish a sensible boundary interface region before
-        // running a circumsphere assessment.
 
-        buildParallelInterfaceIntersection
-        (
-            referralVertices,
-            receivedVertices,
-            outputName
-        );
+    // Reject points that are not near the boundary from the subsequent
+    // searches
 
-        timeCheck("After buildParallelInterfaceIntersection");
-    }
+//    label nearProcCount = 0;
+//    label notNearProcCount = 0;
+//    boundBox quickRejectionBox(decomposition_().procBounds());
+//
+////    Pout<< "Processor boundBox: " << quickRejectionBox << endl;
+//
+//    quickRejectionBox.inflate(-0.1);
+//
+//    OFstream str("rejectionBox_" + name(Pstream::myProcNo()) + ".obj");
+//
+//    meshTools::writeOBJ
+//    (
+//        str,
+//        boundBox::faces(),
+//        quickRejectionBox.points()
+//    );
+//
+//    Pout<< "Inflated rejection boundBox: " << quickRejectionBox << endl;
+
+//    for
+//    (
+//        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+//        vit != finite_vertices_end();
+//        vit++
+//    )
+//    {
+//        if (vit->real() && !vit->nearProcBoundary())
+//        {
+//            const Foam::point& pt = topoint(vit->point());
+//            const scalar range = targetCellSize(pt);
+//
+////            if
+////            (
+////                decomposition_().overlapsThisProcessor
+////                (
+////                    pt,
+////                    range
+////                )
+////            )
+//            if (!quickRejectionBox.contains(pt))
+//            {
+//                vit->setNearProcBoundary();
+//                nearProcCount++;
+//            }
+//            else
+//            {
+//                //vit->setNearProcBoundary();
+//                notNearProcCount++;
+//            }
+//        }
+//    }
+//
+//    reduce(nearProcCount, sumOp<label>());
+//    reduce(notNearProcCount, sumOp<label>());
+//
+//    timeCheck
+//    (
+//        "End of potential intersection search. "
+//      + name(nearProcCount) + " are near a processor boundary. "
+//      + name(notNearProcCount) + " are not."
+//    );
+
+//    if (initialEdgeReferral)
+//    {
+//        // Used as an initial pass to localise the vertex referring - find
+//        // vertices whose dual edges pierce nearby processor volumes and refer
+//        // them to establish a sensible boundary interface region before
+//        // running a circumsphere assessment.
+//
+//        buildParallelInterfaceIntersection
+//        (
+//            referralVertices,
+//            receivedVertices,
+//            outputName
+//        );
+//
+//        timeCheck("After buildParallelInterfaceIntersection");
+//    }
 
     buildParallelInterfaceInfluence
     (
@@ -860,6 +1078,56 @@ void Foam::conformalVoronoiMesh::buildParallelInterface
     );
 
     timeCheck("After buildParallelInterface");
+
+    // Check all referred vertices are actually used on the processor.
+    label nUnusedReferred = numberOfUnusedReferredPoints();
+
+    reduce(nUnusedReferred, sumOp<label>());
+
+    Info<< "    Number of referred points that are not used : "
+        << nUnusedReferred << " (approximate)" << endl;
+}
+
+
+Foam::label Foam::conformalVoronoiMesh::numberOfUnusedReferredPoints() const
+{
+    label nUnusedPoints = 0;
+
+    for
+    (
+        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->referred())
+        {
+            std::list<Vertex_handle> adjVertices;
+            finite_adjacent_vertices(vit, std::back_inserter(adjVertices));
+
+            bool isUsed = false;
+
+            for
+            (
+                std::list<Vertex_handle>::iterator adjVit = adjVertices.begin();
+                adjVit != adjVertices.end();
+                ++adjVit
+            )
+            {
+                if ((*adjVit)->real())
+                {
+                    isUsed = true;
+                }
+            }
+
+            if (!isUsed)
+            {
+                nUnusedPoints++;
+            }
+        }
+    }
+
+    return nUnusedPoints;
 }
 
 
@@ -880,7 +1148,7 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceAll
     (
         Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
         vit != finite_vertices_end();
-        vit++
+        ++vit
     )
     {
         if (!vit->farPoint())
@@ -966,7 +1234,10 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceIntersection
 
         // If either Delaunay cell at the end of the Dual edge is infinite,
         // skip.
-        if (!is_infinite(c1) && !is_infinite(c2))
+        if
+        (
+            !is_infinite(c1) && !is_infinite(c2)
+        )
         {
             // The Delaunauy cells at either end of the dual edge need to be
             // real, i.e. all vertices form part of the internal or boundary
@@ -977,12 +1248,14 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceIntersection
              && c2->internalOrBoundaryDualVertex()
             )
             {
-                Foam::point a = topoint(dual(c1));
-                Foam::point b = topoint(dual(c2));
+                const Foam::point& a = c1->dual();
+                const Foam::point& b = c2->dual();
 
                 // Only if the dual edge cuts the boundary of this processor is
                 // it going to be counted.
-                if (decomposition_().findLineAny(a, b).hit())
+                pointIndexHit info = decomposition_().findLineAny(a, b);
+
+                if (info.hit())
                 {
                     dE0.append(a);
                     dE1.append(b);
@@ -995,11 +1268,22 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceIntersection
         fIOuter++;
     }
 
+    reduce(fIOuter, sumOp<label>());
+
+    timeCheck
+    (
+        "End of actual intersection search over "
+      + name(fIOuter)
+      + " faces."
+    );
+
     // Preform intersections in both directions, as there is no sense
     // associated with the Dual edge
     List<List<pointIndexHit> > intersectionForward(intersectsProc(dE0, dE1));
     List<List<pointIndexHit> > intersectionReverse(intersectsProc(dE1, dE0));
 
+    timeCheck("End of find processor intersection");
+
     // Reset counter
     fIOuter = 0;
 
@@ -1175,13 +1459,58 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
     DynamicList<Foam::point> circumcentre;
     DynamicList<scalar> circumradiusSqr;
 
-    PackedBoolList testCellInfluence(number_of_cells(), false);
 
     // Index outer (all) Delaunauy cells for whether they are potential
     // overlaps, index (inner) the list of tests an results.
     label cIInner = 0;
     label cIOuter = 0;
 
+
+    label cellIndexCount = 0;
+    for
+    (
+        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
+        cit != finite_cells_end();
+        ++cit
+    )
+    {
+        cit->cellIndex() = cellIndexCount++;
+    }
+
+    timeCheck("End of cell Indexing");
+
+    labelList testCellInfluence(number_of_cells(), 0);
+
+    label nQuickRejections = 0;
+
+    for
+    (
+        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
+        cit != finite_cells_end();
+        ++cit
+    )
+    {
+        const Foam::point& cc = cit->dual();
+
+        const scalar crSqr = magSqr(cc - topoint(cit->vertex(0)->point()));
+
+        if
+        (
+            decomposition_().tree().quickCircumsphereRejection
+            (
+                cc,
+                crSqr,
+                decomposition_().octreeNearestDistances()
+            )
+        )
+        {
+            nQuickRejections++;
+            testCellInfluence[cit->cellIndex()] = -1;
+        }
+    }
+
+    timeCheck("End of octreeNearestDistances calculation");
+
     for
     (
         Delaunay::Finite_cells_iterator cit = finite_cells_begin();
@@ -1196,14 +1525,15 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
 
         // The Delaunay cells to assess have to be real, i.e. all vertices form
         // part of the internal or any part of the boundary definition
-        if (cit->real())
+        if
+        (
+            (testCellInfluence[cit->cellIndex()] == 0)
+         && (cit->real() || cit->hasFarPoint())
+        )
         {
-            Foam::point cc(topoint(dual(cit)));
+            const Foam::point& cc = cit->dual();
 
-            scalar crSqr
-            (
-                magSqr(cc - topoint(cit->vertex(0)->point()))
-            );
+            const scalar crSqr = magSqr(cc - topoint(cit->vertex(0)->point()));
 
             // Only if the circumsphere overlaps the boundary of this processor
             // is there a chance of it overlapping others
@@ -1212,7 +1542,7 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
                 circumcentre.append(cc);
                 circumradiusSqr.append(crSqr);
 
-                testCellInfluence[cIOuter] = true;
+                testCellInfluence[cit->cellIndex()] = 1;
             }
         }
 
@@ -1221,6 +1551,9 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
 
     timeCheck("End of testing cell influence");
 
+    Pout<< "Number of quick rejections = " << nQuickRejections << endl;
+    Pout<< "Number of influences = " << circumcentre.size() << endl;
+
     // Increasing the circumspheres to increase the overlaps and compensate for
     // floating point errors missing some referrals
     labelListList circumsphereOverlaps
@@ -1242,7 +1575,7 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
     )
     {
         // Pre-tested circumsphere potential influence
-        if (testCellInfluence[cIOuter])
+        if (testCellInfluence[cit->cellIndex()] == 1)
         {
             const labelList& citOverlaps = circumsphereOverlaps[cIInner];
 
@@ -1289,6 +1622,70 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
         cIOuter++;
     }
 
+
+//    label nFarPoints = removeProcessorBoundarySeeds(true);
+//
+//    reduce(nFarPoints, sumOp<label>());
+//
+//    Info<< "    Removed " << nFarPoints
+//        << " far points from the mesh." << endl;
+
+
+//    seedProcessorBoundarySurfaces(false);
+
+//    cIInner = 0;
+//    cIOuter = 0;
+
+
+    // Relying on the order of iteration of cells being the same as before
+//    for
+//    (
+//        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
+//        cit != finite_cells_end();
+//        ++cit
+//    )
+//    {
+//        // Pre-tested circumsphere potential influence
+//        if (testCellInfluence[cIOuter])
+//        {
+//            const labelList& citOverlaps = circumsphereOverlaps[cIInner];
+//
+//            forAll(citOverlaps, cOI)
+//            {
+//                label procI = citOverlaps[cOI];
+//
+//                recursiveCircumsphereSearch
+//                (
+//                    cit,
+//                    procI,
+//                    referralVertices,
+//                    checkedCells,
+//                    parallelInfluencePoints,
+//                    parallelInfluenceIndices,
+//                    targetProcessor
+//                );
+//            }
+//
+//            cIInner++;
+//        }
+//
+//        cIOuter++;
+//    }
+
+//    for
+//    (
+//        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+//        vit != finite_vertices_end();
+//        ++vit
+//    )
+//    {
+//        if (vit->referred())
+//        {
+//            //Pout << "REMOVE: " << topoint(vit->point()) << endl;
+//            remove(vit);
+//        }
+//    }
+
     referVertices
     (
         targetProcessor,
@@ -1344,6 +1741,8 @@ void Foam::conformalVoronoiMesh::referVertices
         );
     }
 
+    timeCheck("Start of referVertices " + stageName + " insertion.");
+
     for (label procI = 0; procI < Pstream::nProcs(); procI++)
     {
         const labelList& constructMap = pointMap.constructMap()[procI];
@@ -1397,7 +1796,7 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
     hitSurfaceLargest = -1;
 
     std::list<Facet> facets;
-    incident_facets(vit, std::back_inserter(facets));
+    finite_incident_facets(vit, std::back_inserter(facets));
 
     const Foam::point vert = topoint(vit->point());
 
@@ -1410,52 +1809,45 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
         ++fit
     )
     {
-        if
+        const Foam::point edgeMid =
+            0.5
+           *(
+                fit->first->dual()
+              + fit->first->neighbor(fit->second)->dual()
+            );
+
+        pointIndexHit surfHit;
+        label hitSurface;
+
+        geometryToConformTo_.findSurfaceAnyIntersection
         (
-            !is_infinite(fit->first)
-         && !is_infinite(fit->first->neighbor(fit->second))
-        )
-        {
-            const Foam::point edgeMid =
-                0.5
-               *(
-                    topoint(dual(fit->first))
-                  + topoint(dual(fit->first->neighbor(fit->second)))
-                );
+            vert,
+            edgeMid,
+            surfHit,
+            hitSurface
+        );
 
-            pointIndexHit surfHit;
-            label hitSurface;
+        if (surfHit.hit())
+        {
+            vectorField norm(1);
 
-            geometryToConformTo_.findSurfaceAnyIntersection
+            allGeometry_[hitSurface].getNormal
             (
-                vert,
-                edgeMid,
-                surfHit,
-                hitSurface
+                List<pointIndexHit>(1, surfHit),
+                norm
             );
 
-            if (surfHit.hit())
-            {
-                vectorField norm(1);
-
-                allGeometry_[hitSurface].getNormal
-                (
-                    List<pointIndexHit>(1, surfHit),
-                    norm
-                );
-
-                const vector& n = norm[0];
+            const vector& n = norm[0];
 
-                const scalar normalProtrusionDistance =
-                    (edgeMid - surfHit.hitPoint()) & n;
+            const scalar normalProtrusionDistance =
+                (edgeMid - surfHit.hitPoint()) & n;
 
-                if (normalProtrusionDistance > maxProtrusionDistance)
-                {
-                    surfHitLargest = surfHit;
-                    hitSurfaceLargest = hitSurface;
+            if (normalProtrusionDistance > maxProtrusionDistance)
+            {
+                surfHitLargest = surfHit;
+                hitSurfaceLargest = hitSurface;
 
-                    maxProtrusionDistance = normalProtrusionDistance;
-                }
+                maxProtrusionDistance = normalProtrusionDistance;
             }
         }
     }
@@ -1489,9 +1881,9 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
     hitSurfaceLargest = -1;
 
     std::list<Facet> facets;
-    incident_facets(vit, std::back_inserter(facets));
+    finite_incident_facets(vit, std::back_inserter(facets));
 
-    Foam::point vert(topoint(vit->point()));
+    const Foam::point vert = topoint(vit->point());
 
     scalar minIncursionDistance = -maxSurfaceProtrusion(vert);
 
@@ -1502,57 +1894,50 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
         ++fit
     )
     {
-        if
+        const Foam::point edgeMid =
+            0.5
+           *(
+                fit->first->dual()
+              + fit->first->neighbor(fit->second)->dual()
+            );
+
+        pointIndexHit surfHit;
+        label hitSurface;
+
+        geometryToConformTo_.findSurfaceAnyIntersection
         (
-            !is_infinite(fit->first)
-         && !is_infinite(fit->first->neighbor(fit->second))
-        )
-        {
-            Foam::point edgeMid =
-                0.5
-               *(
-                    topoint(dual(fit->first))
-                  + topoint(dual(fit->first->neighbor(fit->second)))
-                );
+            vert,
+            edgeMid,
+            surfHit,
+            hitSurface
+        );
 
-            pointIndexHit surfHit;
-            label hitSurface;
+        if (surfHit.hit())
+        {
+            vectorField norm(1);
 
-            geometryToConformTo_.findSurfaceAnyIntersection
+            allGeometry_[hitSurface].getNormal
             (
-                vert,
-                edgeMid,
-                surfHit,
-                hitSurface
+                List<pointIndexHit>(1, surfHit),
+                norm
             );
 
-            if (surfHit.hit())
-            {
-                vectorField norm(1);
-
-                allGeometry_[hitSurface].getNormal
-                (
-                    List<pointIndexHit>(1, surfHit),
-                    norm
-                );
-
-                const vector& n = norm[0];
+            const vector& n = norm[0];
 
-                scalar normalIncursionDistance =
-                    (edgeMid - surfHit.hitPoint()) & n;
+            scalar normalIncursionDistance =
+                (edgeMid - surfHit.hitPoint()) & n;
 
-                if (normalIncursionDistance < minIncursionDistance)
-                {
-                    surfHitLargest = surfHit;
-                    hitSurfaceLargest = hitSurface;
+            if (normalIncursionDistance < minIncursionDistance)
+            {
+                surfHitLargest = surfHit;
+                hitSurfaceLargest = hitSurface;
 
-                    minIncursionDistance = normalIncursionDistance;
+                minIncursionDistance = normalIncursionDistance;
 
-                    // Info<< nl << "# Incursion: " << endl;
-                    // meshTools::writeOBJ(Info, vert);
-                    // meshTools::writeOBJ(Info, edgeMid);
-                    // Info<< "l Na Nb" << endl;
-                }
+                // Info<< nl << "# Incursion: " << endl;
+                // meshTools::writeOBJ(Info, vert);
+                // meshTools::writeOBJ(Info, edgeMid);
+                // Info<< "l Na Nb" << endl;
             }
         }
     }
@@ -1948,6 +2333,35 @@ bool Foam::conformalVoronoiMesh::appendToEdgeLocationTree
 }
 
 
+Foam::List<Foam::pointIndexHit>
+Foam::conformalVoronoiMesh::nearestFeatureEdgeLocations
+(
+    const Foam::point& pt
+) const
+{
+    const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
+
+    labelList elems
+        = edgeLocationTreePtr_().findSphere(pt, exclusionRangeSqr);
+
+    DynamicList<pointIndexHit> dynPointHit;
+
+    forAll(elems, elemI)
+    {
+        label index = elems[elemI];
+
+        const Foam::point& pointI
+            = edgeLocationTreePtr_().shapes().shapePoints()[index];
+
+        pointIndexHit nearHit(true, pointI, index);
+
+        dynPointHit.append(nearHit);
+    }
+
+    return dynPointHit;
+}
+
+
 bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation
 (
     const Foam::point& pt
@@ -2009,103 +2423,77 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
     DynamicList<Foam::point>& existingEdgeLocations
 ) const
 {
-    const Foam::point pt = pHit.hitPoint();
+    Foam::point pt = pHit.hitPoint();
 
     const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
 
-    pointIndexHit info;
-
-    bool closeToFeatureEdge = pointIsNearFeatureEdgeLocation(pt, info);
+    bool closeToFeatureEdge = pointIsNearFeatureEdgeLocation(pt);
 
-    if (!closeToFeatureEdge)
-    {
-        appendToEdgeLocationTree(pt, existingEdgeLocations);
-    }
-    else
+    if (closeToFeatureEdge)
     {
-        // Check if the edge location that the new edge location is near to
-        // "might" be on a different edge. If so, add it anyway.
-        pointIndexHit edgeHit;
-        label featureHit = -1;
+        List<pointIndexHit> nearHits = nearestFeatureEdgeLocations(pt);
 
-        geometryToConformTo_.findEdgeNearest
-        (
-            pt,
-            exclusionRangeSqr,
-            edgeHit,
-            featureHit
-        );
+        forAll(nearHits, elemI)
+        {
+            pointIndexHit& info = nearHits[elemI];
 
-        const extendedFeatureEdgeMesh& eMesh
-            = geometryToConformTo_.features()[featureHit];
+            // Check if the edge location that the new edge location is near to
+            // "might" be on a different edge. If so, add it anyway.
+            pointIndexHit edgeHit;
+            label featureHit = -1;
 
-        const vector& edgeDir = eMesh.edgeDirections()[edgeHit.index()];
+            geometryToConformTo_.findEdgeNearest
+            (
+                pt,
+                exclusionRangeSqr,
+                edgeHit,
+                featureHit
+            );
 
-        const vector lineBetweenPoints = pt - info.hitPoint();
+            const extendedFeatureEdgeMesh& eMesh
+                = geometryToConformTo_.features()[featureHit];
 
-        const scalar cosAngle = vectorTools::cosPhi(edgeDir, lineBetweenPoints);
+            const vector& edgeDir = eMesh.edgeDirections()[edgeHit.index()];
 
-        // Allow the point to be added if it is almost at right angles to the
-        // other point. Also check it is not the same point.
-//        Info<< cosAngle<< " "
-//            << radToDeg(acos(cosAngle)) << " "
-//            << searchConeAngle << " "
-//            << radToDeg(acos(searchConeAngle)) << endl;
-        if
-        (
-            mag(cosAngle) < searchConeAngle
-         && mag(lineBetweenPoints) > SMALL
-        )
-        {
-            closeToFeatureEdge = false;
-            appendToEdgeLocationTree(pt, existingEdgeLocations);
-        }
-    }
+            const vector lineBetweenPoints = pt - info.hitPoint();
 
-    return closeToFeatureEdge;
+            const scalar cosAngle
+                = vectorTools::cosPhi(edgeDir, lineBetweenPoints);
 
-    // Searching for the nearest point in existingEdgeLocations using the
-    // indexedOctree
+            // Allow the point to be added if it is almost at right angles to
+            // the other point. Also check it is not the same point.
+    //        Info<< cosAngle<< " "
+    //            << radToDeg(acos(cosAngle)) << " "
+    //            << searchConeAngle << " "
+    //            << radToDeg(acos(searchConeAngle)) << endl;
 
-    // Average the points...
-//    if (info.hit())
-//    {
-//        Foam::point newPt = 0.5*(info.hitPoint() + pt);
-//
-//        pHit.setPoint(newPt);
-//
-//        //boolList toRemove(existingEdgeLocations.size(), false);
-//
-//        forAll(existingEdgeLocations, pI)
-//        {
-//            if (pI == info.index())
-//            {
-//                //toRemove[pI] = true;
-//                edgeLocationTree.remove(pI);
-//            }
-//        }
-////
-////        pointField newExistingEdgeLocations(existingEdgeLocations.size());
-////
-////        label count = 0;
-////        forAll(existingEdgeLocations, pI)
-////        {
-////            if (toRemove[pI] == false)
-////            {
-////                newExistingEdgeLocations[count++] =
-////                    existingEdgeLocations[pI];
-////            }
-////        }
-////
-////        newExistingEdgeLocations.resize(count);
-////
-////        existingEdgeLocations = newExistingEdgeLocations;
-////
-////        existingEdgeLocations.append(newPt);
-//
-//        return !info.hit();
-//    }
+            if
+            (
+                mag(cosAngle) < searchConeAngle
+             && (
+                    mag(lineBetweenPoints)
+                  > cvMeshControls().pointPairDistanceCoeff()*targetCellSize(pt)
+                )
+            )
+            {
+                pt = edgeHit.hitPoint();
+                pHit.setPoint(pt);
+                closeToFeatureEdge = false;
+            }
+            else
+            {
+                closeToFeatureEdge = true;
+                break;
+            }
+        }
+    }
+
+    if (!closeToFeatureEdge)
+    {
+        appendToEdgeLocationTree(pt, existingEdgeLocations);
+    }
 
+    return closeToFeatureEdge;
 }
 
 
@@ -2265,12 +2653,6 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
             featuresHit
         );
 
-        // Gather edge locations but do not add them to newEdgeLocations inside
-        // the loop as they will prevent nearby edge locations of different
-        // types being conformed to.
-
-        DynamicList<Foam::point> currentEdgeLocations;
-
         forAll(edHitsByFeature, i)
         {
             const label featureHit = featuresHit[i];
@@ -2281,39 +2663,43 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
             {
                 pointIndexHit& edHit = edHits[eHitI];
 
-                if (!nearFeaturePt(edHit.hitPoint()) && keepSurfacePoint)
+                if (edHit.hit())
                 {
-                    if
-                    (
-                        magSqr(edHit.hitPoint() - surfHitI.hitPoint())
-                      < surfacePtReplaceDistCoeffSqr*cellSizeSqr
-                    )
+                    if (!nearFeaturePt(edHit.hitPoint()))
                     {
-                        // If the point is within a given distance of a feature
-                        // edge, give control to edge control points instead,
-                        // this will prevent "pits" forming.
+                        if
+                        (
+                            magSqr(edHit.hitPoint() - surfHitI.hitPoint())
+                          < surfacePtReplaceDistCoeffSqr*cellSizeSqr
+                        )
+                        {
+                            // If the point is within a given distance of a
+                            // feature edge, give control to edge control points
+                            // instead, this will prevent "pits" forming.
 
-                        keepSurfacePoint = false;
+                            keepSurfacePoint = false;
 
-                        // NEED TO REMOVE FROM THE SURFACE TREE...
-                    }
+                            // NEED TO REMOVE FROM THE SURFACE TREE...
+                            surfacePtLocationTreePtr_().remove
+                            (
+                                existingSurfacePtLocations.size()
+                            );
+                        }
 
-                    if
-                    (
-                        !nearFeatureEdgeLocation
+                        if
                         (
-                            edHit,
-                            existingEdgeLocations
+                            !nearFeatureEdgeLocation
+                            (
+                                edHit,
+                                existingEdgeLocations
+                            )
                         )
-                    )
-                    {
-                        // Do not place edge control points too close to a
-                        // feature point or existing edge control points
-
-                        featureEdgeHits.append(edHit);
-                        featureEdgeFeaturesHit.append(featureHit);
-
-                        currentEdgeLocations.append(edHit.hitPoint());
+                        {
+                            // Do not place edge control points too close to a
+                            // feature point or existing edge control points
+                            featureEdgeHits.append(edHit);
+                            featureEdgeFeaturesHit.append(featureHit);
+                        }
                     }
                 }
             }
@@ -2322,7 +2708,6 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
         if (keepSurfacePoint)
         {
             surfaceHits.append(surfHitI);
-
             hitSurfaces.append(hitSurfaceI);
         }
     }
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index ab7ce33f8563760f144c6aff33f6893a1a7e0078..89afb87134a98ce7b7c9c5ec37afe4764ed4f4b2 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -125,7 +125,7 @@ void Foam::conformalVoronoiMesh::drawDelaunayCell
     os  << "# cell index: " << label(c->cellIndex()) << endl;
 
     os  << "# circumradius "
-        << mag(topoint(dual(c)) - topoint(c->vertex(0)->point()))
+        << mag(c->dual() - topoint(c->vertex(0)->point()))
         << endl;
 
     for (int i = 0; i < 4; i++)
@@ -144,7 +144,7 @@ void Foam::conformalVoronoiMesh::drawDelaunayCell
 
     os  << "# cicumcentre " << endl;
 
-    meshTools::writeOBJ(os, topoint(dual(c)));
+    meshTools::writeOBJ(os, c->dual());
 
     os  << "l " << 1 + offset << " " << 5 + offset << endl;
 }
@@ -167,7 +167,7 @@ void Foam::conformalVoronoiMesh::writePoints
         ++vit
     )
     {
-        if (!internalOnly || vit->internalOrBoundaryPoint())
+        if (!internalOnly || vit->internalPoint())
         {
             meshTools::writeOBJ(str, topoint(vit->point()));
         }
@@ -241,7 +241,7 @@ void Foam::conformalVoronoiMesh::writeProcessorInterface
     {
         if (!cit->farCell())
         {
-            points[cit->cellIndex()] = topoint(dual(cit));
+            points[cit->cellIndex()] = cit->dual();
         }
     }
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
index 7ab4f10853427619285ce87430f1ccf14cd31e42..5fc7e970b1131c78607d13f71c91ff3b90aeb276 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
@@ -51,6 +51,7 @@ SourceFiles
 #include "Swap.H"
 #include "InfoProxy.H"
 #include "tetCell.H"
+#include "typeInfo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -154,6 +155,8 @@ public:
 
         inline int cellIndex() const;
 
+        inline const Foam::point& dual();
+
         inline bool farCell() const;
 
         inline int& filterCount();
@@ -163,6 +166,11 @@ public:
         //- Is the Delaunay cell real, i.e. any real vertex
         inline bool real() const;
 
+        //- Does the Delaunay cell have a far point
+        inline bool hasFarPoint() const;
+
+        inline bool hasInternalPoint() const;
+
         //- Does the Dual vertex form part of a processor patch
         inline bool parallelDualVertex() const;
 
@@ -190,6 +198,8 @@ public:
         //  least one Delaunay vertex outside and at least one inside
         inline bool boundaryDualVertex() const;
 
+        inline bool nearProcBoundary() const;
+
 
     // Info
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H
index 5fbd378f9d9a220845245127070d9624f2ac9bba..6450b01908a44149dbde6d8abc6d8e927f6083b1 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H
@@ -86,6 +86,24 @@ int CGAL::indexedCell<Gt, Cb>::cellIndex() const
 }
 
 
+template<class Gt, class Cb>
+const Foam::point& CGAL::indexedCell<Gt, Cb>::dual()
+{
+#ifdef CGAL_INEXACT
+    return reinterpret_cast<const Foam::point&>(this->circumcenter());
+#else
+    const typename Gt::Point_3& P = this->circumcenter();
+
+    return
+    (
+        CGAL::to_double(P.x()),
+        CGAL::to_double(P.y()),
+        CGAL::to_double(P.z())
+    );
+#endif
+}
+
+
 template<class Gt, class Cb>
 inline bool CGAL::indexedCell<Gt, Cb>::farCell() const
 {
@@ -112,10 +130,45 @@ inline bool CGAL::indexedCell<Gt, Cb>::real() const
 {
     return
     (
-        this->vertex(0)->real()
-     || this->vertex(1)->real()
-     || this->vertex(2)->real()
-     || this->vertex(3)->real()
+        (
+            this->vertex(0)->real()
+         || this->vertex(1)->real()
+         || this->vertex(2)->real()
+         || this->vertex(3)->real()
+        )
+        &&
+        !(
+            this->vertex(0)->farPoint()
+         || this->vertex(1)->farPoint()
+         || this->vertex(2)->farPoint()
+         || this->vertex(3)->farPoint()
+        )
+    );
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::hasFarPoint() const
+{
+    return
+    (
+        this->vertex(0)->farPoint()
+     || this->vertex(1)->farPoint()
+     || this->vertex(2)->farPoint()
+     || this->vertex(3)->farPoint()
+    );
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::hasInternalPoint() const
+{
+    return
+    (
+        this->vertex(0)->internalPoint()
+     || this->vertex(1)->internalPoint()
+     || this->vertex(2)->internalPoint()
+     || this->vertex(3)->internalPoint()
     );
 }
 
@@ -285,4 +338,17 @@ inline bool CGAL::indexedCell<Gt, Cb>::boundaryDualVertex() const
 }
 
 
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::nearProcBoundary() const
+{
+    return
+    (
+        this->vertex(0)->nearProcBoundary()
+     || this->vertex(1)->nearProcBoundary()
+     || this->vertex(2)->nearProcBoundary()
+     || this->vertex(3)->nearProcBoundary()
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.C
index 9ad21f92feade712950e9aed5d5cb6d1a21a3eaf..2173e7f5fa7b6be6a0469b004a46e9148f4195f7 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.C
@@ -29,7 +29,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "indexedVertex.H"
-//#include "conformalVoronoiMesh.H"
 #include "point.H"
 
 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
index 246020fd172bbeb127089dc9bcc085c0aa7908e2..525c2156e82b60682ca93eb20d44fc4d34d051d7 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
@@ -45,6 +45,7 @@ SourceFiles
 #include <CGAL/Triangulation_3.h>
 #include "tensor.H"
 #include "InfoProxy.H"
+#include "point.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -101,6 +102,8 @@ class indexedVertex
         //- Specify whether the vertex is fixed or movable.
         bool vertexFixed_;
 
+        bool nearProcBoundary_;
+
 
 public:
 
@@ -198,6 +201,12 @@ public:
         //- Set the point to be near the boundary
         inline void setNearBoundary();
 
+        //- Is point internal and near a proc boundary
+        inline bool nearProcBoundary() const;
+
+        //- Set the point to be near a proc boundary
+        inline void setNearProcBoundary();
+
         //- Either master or slave of pointPair.
         inline bool pairPoint() const;
 
@@ -227,7 +236,7 @@ public:
         inline bool isVertexFixed() const;
 
         //- Fix the vertex so that it can't be moved
-        inline void setVertexFixed() const;
+        inline void setVertexFixed();
 
     // inline void operator=(const Delaunay::Finite_vertices_iterator vit)
     // {
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H
index c02ec4df2009489efa9fe91221a666b9f9ae0050..1e225eaa2ed02a2674b57a713fbc142f05d451b4 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H
@@ -38,7 +38,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex()
     type_(vtInternal),
     alignment_(),
     targetCellSize_(0.0),
-    vertexFixed_(false)
+    vertexFixed_(false),
+    nearProcBoundary_(false)
 {}
 
 
@@ -50,7 +51,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p)
     type_(vtInternal),
     alignment_(),
     targetCellSize_(0.0),
-    vertexFixed_(false)
+    vertexFixed_(false),
+    nearProcBoundary_(false)
 {}
 
 
@@ -67,7 +69,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex
     type_(type),
     alignment_(),
     targetCellSize_(0.0),
-    vertexFixed_(false)
+    vertexFixed_(false),
+    nearProcBoundary_(false)
 {}
 
 
@@ -79,7 +82,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p, Cell_handle f)
     type_(vtInternal),
     alignment_(),
     targetCellSize_(0.0),
-    vertexFixed_(false)
+    vertexFixed_(false),
+    nearProcBoundary_(false)
 {}
 
 
@@ -91,7 +95,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(Cell_handle f)
     type_(vtInternal),
     alignment_(),
     targetCellSize_(0.0),
-    vertexFixed_(false)
+    vertexFixed_(false),
+    nearProcBoundary_(false)
 {}
 
 
@@ -251,6 +256,20 @@ inline void CGAL::indexedVertex<Gt, Vb>::setNearBoundary()
 }
 
 
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::nearProcBoundary() const
+{
+    return nearProcBoundary_;
+}
+
+
+template<class Gt, class Vb>
+inline void CGAL::indexedVertex<Gt, Vb>::setNearProcBoundary()
+{
+    nearProcBoundary_ = true;
+}
+
+
 template<class Gt, class Vb>
 inline bool CGAL::indexedVertex<Gt, Vb>::pairPoint() const
 {
@@ -331,7 +350,7 @@ inline bool CGAL::indexedVertex<Gt, Vb>::isVertexFixed() const
 
 
 template<class Gt, class Vb>
-inline void CGAL::indexedVertex<Gt, Vb>::setVertexFixed() const
+inline void CGAL::indexedVertex<Gt, Vb>::setVertexFixed()
 {
     vertexFixed_ = true;
 }