diff --git a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/files b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/files
new file mode 100644
index 0000000000..73871a7532
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/files
@@ -0,0 +1,3 @@
+checkCvMesh.C
+
+EXE = $(FOAM_APPBIN)/checkCvMesh
diff --git a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/options
new file mode 100644
index 0000000000..ba68fd3819
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/options
@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/mesh/autoMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -ldynamicMesh \
+    -ltriSurface \
+    -lautoMesh \
+    -lmeshTools
+
diff --git a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/checkCvMesh.C b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/checkCvMesh.C
new file mode 100644
index 0000000000..df2a3b90f8
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/checkCvMesh.C
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    checkCvMesh
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "autoSnapDriver.H"
+#include "faceSet.H"
+#include "motionSmoother.H"
+#include "timeSelector.H"
+
+
+using namespace Foam;
+
+
+int main(int argc, char *argv[])
+{
+#   include "addOverwriteOption.H"
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+    instantList timeDirs = timeSelector::select0(runTime, args);
+
+#   include "createMesh.H"
+
+    runTime.functionObjects().off();
+
+    forAll(timeDirs, timeI)
+    {
+        runTime.setTime(timeDirs[timeI], timeI);
+
+        Info<< "Create mesh for time = " << runTime.timeName()
+            << nl << endl;
+
+        mesh.readUpdate();
+
+        Info<< "Read mesh in = "
+            << runTime.cpuTimeIncrement() << " s" << endl;
+
+        // Check patches and faceZones are synchronised
+        mesh.boundaryMesh().checkParallelSync(true);
+        meshRefinement::checkCoupledFaceZones(mesh);
+
+        // Read meshing dictionary
+        IOdictionary cvMeshDict
+        (
+           IOobject
+           (
+                "cvMeshDict",
+                runTime.system(),
+                mesh,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE
+           )
+        );
+
+        // mesh motion and mesh quality parameters
+        const dictionary& meshQualityDict
+            = cvMeshDict.subDict("meshQualityControls");
+
+
+        Info<< "Checking initial mesh ..." << endl;
+        faceSet wrongFaces(mesh, "wrongFaces", mesh.nFaces()/100);
+        motionSmoother::checkMesh(false, mesh, meshQualityDict, wrongFaces);
+
+        const label nInitErrors = returnReduce
+        (
+            wrongFaces.size(),
+            sumOp<label>()
+        );
+
+        Info<< "Detected " << nInitErrors << " illegal faces"
+            << " (concave, zero area or negative cell pyramid volume)"
+            << endl;
+
+        if (nInitErrors > 0)
+        {
+            Info<< "Writing " << nInitErrors
+                << " faces in error to set "
+                << wrongFaces.name() << endl;
+
+            wrongFaces.instance() = mesh.pointsInstance();
+            wrongFaces.write();
+        }
+
+        Info<< nl << "End of time " << runTime.timeName() << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+
+}
+
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index edafdb50f9..a2d44f0ee1 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -724,6 +724,8 @@ void Foam::backgroundMeshDecomposition::buildPatchAndTree()
 
     globalBackgroundBounds_ = treeBoundBox(bbMin, bbMax);
 
+    octreeNearestDistances_ = bFTreePtr_().calcNearestDistance();
+
     if (cvMesh_.cvMeshControls().objOutput())
     {
         OFstream fStr
@@ -795,6 +797,7 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
     ),
     boundaryFacesPtr_(),
     bFTreePtr_(),
+    octreeNearestDistances_(),
     allBackgroundMeshBounds_(Pstream::nProcs()),
     globalBackgroundBounds_(),
     decomposeDict_
@@ -1150,6 +1153,8 @@ bool Foam::backgroundMeshDecomposition::positionOnThisProcessor
     const point& pt
 ) const
 {
+//    return bFTreePtr_().findAnyOverlap(pt, 0.0);
+
     return
         bFTreePtr_().getVolumeType(pt)
      == indexedOctree<treeDataBPatch>::INSIDE;
@@ -1176,6 +1181,7 @@ bool Foam::backgroundMeshDecomposition::overlapsThisProcessor
     const treeBoundBox& box
 ) const
 {
+//    return !procBounds().contains(box);
     return !bFTreePtr_().findBox(box).empty();
 }
 
@@ -1183,9 +1189,11 @@ bool Foam::backgroundMeshDecomposition::overlapsThisProcessor
 bool Foam::backgroundMeshDecomposition::overlapsThisProcessor
 (
     const point& centre,
-    scalar radiusSqr
+    const scalar radiusSqr
 ) const
 {
+    //return bFTreePtr_().findAnyOverlap(centre, radiusSqr);
+
     return bFTreePtr_().findNearest(centre, radiusSqr).hit();
 }
 
@@ -1645,6 +1653,7 @@ Foam::labelListList Foam::backgroundMeshDecomposition::overlapsProcessors
 
         // If the sphere finds a nearest element of the patch, then it overlaps
         sphereOverlapsCandidate[sI] = bFTreePtr_().findNearest(c, rSqr).hit();
+        //sphereOverlapsCandidate[sI] = bFTreePtr_().findAnyOverlap(c, rSqr);
     }
 
     map().reverseDistribute
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
index 87ab88af35..d8e6388cfb 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
@@ -111,6 +111,8 @@ class backgroundMeshDecomposition
         //- Search tree for the boundaryFaces_ patch
         autoPtr<indexedOctree<treeDataBPatch> > bFTreePtr_;
 
+        List<scalar> octreeNearestDistances_;
+
         //- The bounds of all background meshes on all processors
         treeBoundBoxList allBackgroundMeshBounds_;
 
@@ -225,16 +227,16 @@ public:
         //- Are the given positions inside the domain of this decomposition
         boolList positionOnThisProcessor(const List<point>& pts) const;
 
-        //- Does the given box overlap the faces of the bounday of this
+        //- Does the given box overlap the faces of the boundary of this
         //  processor
         bool overlapsThisProcessor(const treeBoundBox& box) const;
 
-        //- Does the given sphere overlap the faces of the bounday of this
+        //- Does the given sphere overlap the faces of the boundary of this
         //  processor
         bool overlapsThisProcessor
         (
             const point& centre,
-            scalar radiusSqr
+            const scalar radiusSqr
         ) const;
 
         //- Find nearest intersection of line between start and end, (exposing
@@ -289,6 +291,12 @@ public:
             //- Return access to the underlying mesh
             inline const fvMesh& mesh() const;
 
+            //- Return access to the underlying tree
+            inline const indexedOctree<treeDataBPatch>& tree() const;
+
+            //- Return access to the nearest distance of the octree nodes
+            inline const List<scalar>& octreeNearestDistances() const;
+
             //- Return the boundBox of this processor
             inline const treeBoundBox& procBounds() const;
 };
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecompositionI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecompositionI.H
index df659cfd9c..2897d3fc58 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecompositionI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecompositionI.H
@@ -30,6 +30,18 @@ const Foam::fvMesh& Foam::backgroundMeshDecomposition::mesh() const
     return mesh_;
 }
 
+const Foam::indexedOctree<Foam::treeDataBPatch>&
+Foam::backgroundMeshDecomposition::tree() const
+{
+    return bFTreePtr_();
+}
+
+const Foam::List<Foam::scalar>&
+Foam::backgroundMeshDecomposition::octreeNearestDistances() const
+{
+    return octreeNearestDistances_;
+}
+
 
 const Foam::treeBoundBox& Foam::backgroundMeshDecomposition::procBounds() const
 {
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index e8cedc49c2..958fcf6763 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 b590eb31cd..a514c660a1 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 04ef2af045..d103d3c0fc 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 847974a2f4..7f4c7abb6e 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
-        (
-            !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)))
-                );
-
-            pointIndexHit surfHit;
-            label hitSurface;
-
-            geometryToConformTo_.findSurfaceAnyIntersection
-            (
-                vert,
-                edgeMid,
-                surfHit,
-                hitSurface
+        const Foam::point edgeMid =
+            0.5
+           *(
+                fit->first->dual()
+              + fit->first->neighbor(fit->second)->dual()
             );
 
-            if (surfHit.hit())
+        pointIndexHit surfHit;
+        label hitSurface;
+
+        geometryToConformTo_.findSurfaceAnyIntersection
+        (
+            vert,
+            edgeMid,
+            surfHit,
+            hitSurface
+        );
+
+        if (surfHit.hit())
+        {
+            vectorField norm(1);
+
+            allGeometry_[hitSurface].getNormal
+            (
+                List<pointIndexHit>(1, surfHit),
+                norm
+            );
+
+            const vector& n = norm[0];
+
+            const scalar normalProtrusionDistance =
+                (edgeMid - surfHit.hitPoint()) & n;
+
+            if (normalProtrusionDistance > maxProtrusionDistance)
             {
-                vectorField norm(1);
+                surfHitLargest = surfHit;
+                hitSurfaceLargest = hitSurface;
 
-                allGeometry_[hitSurface].getNormal
-                (
-                    List<pointIndexHit>(1, surfHit),
-                    norm
-                );
-
-                const vector& n = norm[0];
-
-                const scalar normalProtrusionDistance =
-                    (edgeMid - surfHit.hitPoint()) & n;
-
-                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
-        (
-            !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)))
-                );
-
-            pointIndexHit surfHit;
-            label hitSurface;
-
-            geometryToConformTo_.findSurfaceAnyIntersection
-            (
-                vert,
-                edgeMid,
-                surfHit,
-                hitSurface
+        const Foam::point edgeMid =
+            0.5
+           *(
+                fit->first->dual()
+              + fit->first->neighbor(fit->second)->dual()
             );
 
-            if (surfHit.hit())
+        pointIndexHit surfHit;
+        label hitSurface;
+
+        geometryToConformTo_.findSurfaceAnyIntersection
+        (
+            vert,
+            edgeMid,
+            surfHit,
+            hitSurface
+        );
+
+        if (surfHit.hit())
+        {
+            vectorField norm(1);
+
+            allGeometry_[hitSurface].getNormal
+            (
+                List<pointIndexHit>(1, surfHit),
+                norm
+            );
+
+            const vector& n = norm[0];
+
+            scalar normalIncursionDistance =
+                (edgeMid - surfHit.hitPoint()) & n;
+
+            if (normalIncursionDistance < minIncursionDistance)
             {
-                vectorField norm(1);
+                surfHitLargest = surfHit;
+                hitSurfaceLargest = hitSurface;
 
-                allGeometry_[hitSurface].getNormal
-                (
-                    List<pointIndexHit>(1, surfHit),
-                    norm
-                );
+                minIncursionDistance = normalIncursionDistance;
 
-                const vector& n = norm[0];
-
-                scalar normalIncursionDistance =
-                    (edgeMid - surfHit.hitPoint()) & n;
-
-                if (normalIncursionDistance < minIncursionDistance)
-                {
-                    surfHitLargest = surfHit;
-                    hitSurfaceLargest = hitSurface;
-
-                    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);
 
-    bool closeToFeatureEdge = pointIsNearFeatureEdgeLocation(pt, info);
+    if (closeToFeatureEdge)
+    {
+        List<pointIndexHit> nearHits = nearestFeatureEdgeLocations(pt);
+
+        forAll(nearHits, elemI)
+        {
+            pointIndexHit& info = nearHits[elemI];
+
+            // 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;
+
+            geometryToConformTo_.findEdgeNearest
+            (
+                pt,
+                exclusionRangeSqr,
+                edgeHit,
+                featureHit
+            );
+
+            const extendedFeatureEdgeMesh& eMesh
+                = geometryToConformTo_.features()[featureHit];
+
+            const vector& edgeDir = eMesh.edgeDirections()[edgeHit.index()];
+
+            const vector lineBetweenPoints = pt - info.hitPoint();
+
+            const scalar cosAngle
+                = vectorTools::cosPhi(edgeDir, lineBetweenPoints);
+
+            // 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)
+                  > cvMeshControls().pointPairDistanceCoeff()*targetCellSize(pt)
+                )
+            )
+            {
+                pt = edgeHit.hitPoint();
+                pHit.setPoint(pt);
+                closeToFeatureEdge = false;
+            }
+            else
+            {
+                closeToFeatureEdge = true;
+                break;
+            }
+        }
+    }
 
     if (!closeToFeatureEdge)
     {
         appendToEdgeLocationTree(pt, existingEdgeLocations);
     }
-    else
-    {
-        // 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;
-
-        geometryToConformTo_.findEdgeNearest
-        (
-            pt,
-            exclusionRangeSqr,
-            edgeHit,
-            featureHit
-        );
-
-        const extendedFeatureEdgeMesh& eMesh
-            = geometryToConformTo_.features()[featureHit];
-
-        const vector& edgeDir = eMesh.edgeDirections()[edgeHit.index()];
-
-        const vector lineBetweenPoints = pt - info.hitPoint();
-
-        const scalar cosAngle = vectorTools::cosPhi(edgeDir, lineBetweenPoints);
-
-        // 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);
-        }
-    }
 
     return closeToFeatureEdge;
-
-    // Searching for the nearest point in existingEdgeLocations using the
-    // indexedOctree
-
-    // 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();
-//    }
-
 }
 
 
@@ -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.
-
-                        keepSurfacePoint = false;
-
-                        // NEED TO REMOVE FROM THE SURFACE TREE...
-                    }
-
-                    if
-                    (
-                        !nearFeatureEdgeLocation
+                        if
                         (
-                            edHit,
-                            existingEdgeLocations
+                            magSqr(edHit.hitPoint() - surfHitI.hitPoint())
+                          < surfacePtReplaceDistCoeffSqr*cellSizeSqr
                         )
-                    )
-                    {
-                        // Do not place edge control points too close to a
-                        // feature point or existing edge control points
+                        {
+                            // If the point is within a given distance of a
+                            // feature edge, give control to edge control points
+                            // instead, this will prevent "pits" forming.
 
-                        featureEdgeHits.append(edHit);
-                        featureEdgeFeaturesHit.append(featureHit);
+                            keepSurfacePoint = false;
 
-                        currentEdgeLocations.append(edHit.hitPoint());
+                            // NEED TO REMOVE FROM THE SURFACE TREE...
+                            surfacePtLocationTreePtr_().remove
+                            (
+                                existingSurfacePtLocations.size()
+                            );
+                        }
+
+                        if
+                        (
+                            !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);
+                        }
                     }
                 }
             }
@@ -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 ab7ce33f85..89afb87134 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 7ab4f10853..5fc7e970b1 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 5fbd378f9d..6450b01908 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 9ad21f92fe..2173e7f5fa 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 246020fd17..525c2156e8 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 c02ec4df20..1e225eaa2e 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;
 }
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
index d8e5fd0166..ce56a977aa 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
@@ -285,7 +285,9 @@ public:
             ) const;
 
             //- Find the nearest points on each feature edge that is within
-            //  a given distance from the sample point
+            //  a given distance from the sample point. Will need to check for
+            //  a hit or a miss because near edges may not have a nearest point
+            //  on them which is perpendicular to the sample point.
             void findAllNearestEdges
             (
                 const point& sample,
diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
index 17ef34fb6f..024f7d27ec 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
+++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
@@ -379,6 +379,124 @@ Foam::label Foam::indexedOctree<Type>::compactContents
 }
 
 
+template <class Type>
+bool Foam::indexedOctree<Type>::quickCircumsphereRejection
+(
+    const label nodeI,
+    const point& cc,
+    const scalar crSqr,
+    const List<scalar>& nearestDistances
+) const
+{
+    const node& nod = nodes_[nodeI];
+
+    volumeType nodeType = volumeType(nodeTypes_.get(nodeI<<3));
+
+    //scalar boxDist = nearestDistances[nodeI] + 0.5*nod.bb_.mag();
+    scalar boxDist = crSqr + magSqr(cc - nod.bb_.midpoint());
+
+    if
+    (
+         nodeType == INSIDE
+      //&& (crSqr < sqr(boxDist))
+      && (boxDist < sqr(nearestDistances[nodeI]))
+    )
+    {
+        return true;
+    }
+    else
+    {
+        direction octant = nod.bb_.subOctant(cc);
+
+        labelBits index = nod.subNodes_[octant];
+
+        if (isNode(index))
+        {
+            return quickCircumsphereRejection
+            (
+                getNode(index),
+                cc,
+                crSqr,
+                nearestDistances
+            );
+        }
+        else
+        {
+            return false;
+        }
+    }
+}
+
+
+template <class Type>
+bool Foam::indexedOctree<Type>::quickCircumsphereRejection
+(
+    const point& cc,
+    const scalar crSqr,
+    const List<scalar>& nearestDistances
+) const
+{
+    if (nodes_.size())
+    {
+        return quickCircumsphereRejection
+        (
+            0,
+            cc,
+            crSqr,
+            nearestDistances
+        );
+    }
+
+    return false;
+}
+
+
+template <class Type>
+Foam::scalar
+Foam::indexedOctree<Type>::calcNearestDistance
+(
+    const label nodeI
+) const
+{
+    const node& nod = nodes_[nodeI];
+
+    const point& nodeCentre = nod.bb_.midpoint();
+
+    scalar nearestDistance = 0.0;
+
+    pointIndexHit pHit = findNearest(nodeCentre, sqr(GREAT));
+
+    if (pHit.hit())
+    {
+        nearestDistance = mag(pHit.hitPoint() - nodeCentre);
+    }
+    else
+    {
+        WarningIn("Foam::indexedOctree<Type>::calcNearestDistance(const label)")
+            << "Cannot calculate distance of nearest point on surface from "
+            << "the midpoint of the octree node. Returning distance of zero."
+            << endl;
+    }
+
+    return nearestDistance;
+}
+
+
+template <class Type>
+Foam::List<Foam::scalar>
+Foam::indexedOctree<Type>::calcNearestDistance() const
+{
+    List<scalar> nearestDistances(nodes_.size());
+
+    forAll(nearestDistances, nodeI)
+    {
+        nearestDistances[nodeI] = calcNearestDistance(nodeI);
+    }
+
+    return nearestDistances;
+}
+
+
 // Pre-calculates wherever possible the volume status per node/subnode.
 // Recurses to determine status of lowest level boxes. Level above is
 // combination of octants below.
@@ -540,6 +658,67 @@ Foam::indexedOctree<Type>::getSide
 // ~~~~~~~~~~~~~~
 //
 
+
+//template <class Type>
+//bool Foam::indexedOctree<Type>::findAnyOverlap
+//(
+//    const label nodeI,
+//    const point& sample,
+//    const scalar nearestDistSqr
+//) const
+//{
+//    const node& nod = nodes_[nodeI];
+//
+//    // Determine order to walk through octants
+//    FixedList<direction, 8> octantOrder;
+//    nod.bb_.searchOrder(sample, octantOrder);
+//
+//    // Go into all suboctants (one containing sample first) and update
+//    // nearest.
+//    for (direction i = 0; i < 8; i++)
+//    {
+//        direction octant = octantOrder[i];
+//
+//        labelBits index = nod.subNodes_[octant];
+//
+//        if (isNode(index))
+//        {
+//            label subNodeI = getNode(index);
+//
+//            const treeBoundBox& subBb = nodes_[subNodeI].bb_;
+//
+//            if (overlaps(subBb.min(), subBb.max(), nearestDistSqr, sample))
+//            {
+//                return findAnyOverlap
+//                (
+//                    subNodeI,
+//                    sample,
+//                    nearestDistSqr
+//                );
+//            }
+//        }
+//        else if (isContent(index))
+//        {
+//            if
+//            (
+//                overlaps
+//                (
+//                    nod.bb_,
+//                    octant,
+//                    nearestDistSqr,
+//                    sample
+//                )
+//            )
+//            {
+//                return true;
+//            }
+//        }
+//    }
+//
+//    return false;
+//}
+
+
 // Find nearest point starting from nodeI
 template <class Type>
 void Foam::indexedOctree<Type>::findNearest
@@ -1614,7 +1793,6 @@ void Foam::indexedOctree<Type>::traverseNode
         }
     }
 
-
     const node& nod = nodes_[nodeI];
 
     labelBits index = nod.subNodes_[octant];
@@ -1781,6 +1959,11 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
     label i = 0;
     for (; i < 100000; i++)
     {
+//        if (isLineInsideOrOutside(nodeI, treeStart, treeEnd))
+//        {
+//            return hitInfo;
+//        }
+
         // Ray-trace to end of current node. Updates point (either on triangle
         // in case of hit or on node bounding box in case of miss)
 
@@ -1935,6 +2118,38 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
 }
 
 
+//template <class Type>
+//bool Foam::indexedOctree<Type>::isLineInsideOrOutside
+//(
+//    const label nodeI,
+//    const point& start,
+//    const point& end
+//) const
+//{
+//    const node& nod = nodes_[nodeI];
+//
+//    direction startOctant = nod.bb_.subOctant(start);
+//    direction endOctant = nod.bb_.subOctant(end);
+//
+//    if (startOctant == endOctant)
+//    {
+//        volumeType startOctantType
+//            = volumeType(nodeTypes_.get((nodeI<<3) + startOctant));
+//
+//        if
+//        (
+//            startOctantType == INSIDE || startOctantType == OUTSIDE
+//        )
+//        {
+//            //Info<< nodeI << " | " << start << " " << end << endl;
+//            return true;
+//        }
+//    }
+//
+//    return false;
+//}
+
+
 // Find first intersection
 template <class Type>
 Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
@@ -2559,6 +2774,27 @@ Foam::scalar& Foam::indexedOctree<Type>::perturbTol()
 }
 
 
+//template <class Type>
+//bool Foam::indexedOctree<Type>::findAnyOverlap
+//(
+//    const point& sample,
+//    const scalar startDistSqr
+//) const
+//{
+//    if (nodes_.size())
+//    {
+//        return findAnyOverlap
+//        (
+//            0,
+//            sample,
+//            startDistSqr
+//        );
+//    }
+//
+//    return false;
+//}
+
+
 template <class Type>
 Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
 (
diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H
index aa493f34e9..40ee7a8a32 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H
+++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H
@@ -201,6 +201,16 @@ private:
                 label& compactI
             );
 
+            scalar calcNearestDistance(const label nodeI) const;
+
+            bool quickCircumsphereRejection
+            (
+                const label nodeI,
+                const point& cc,
+                const scalar crSqr,
+                const List<scalar>& nearestDistances
+            ) const;
+
             //- Determine inside/outside per node (mixed if cannot be
             //  determined). Only valid for closed shapes.
             volumeType calcVolumeType(const label nodeI) const;
@@ -320,6 +330,13 @@ private:
                 const bool verbose = false
             ) const;
 
+//            bool isLineInsideOrOutside
+//            (
+//                const label nodeI,
+//                const point& start,
+//                const point& end
+//            ) const;
+
             //- Find any or nearest intersection of line between start and end.
             pointIndexHit findLine
             (
@@ -532,6 +549,19 @@ public:
                 const scalar nearestDistSqr
             ) const;
 
+//            bool findAnyOverlap
+//            (
+//                const point& sample,
+//                const scalar nearestDistSqr
+//            ) const;
+//
+//            bool findAnyOverlap
+//            (
+//                const label nodeI,
+//                const point& sample,
+//                const scalar nearestDistSqr
+//            ) const;
+
             //- Low level: calculate nearest starting from subnode.
             void findNearest
             (
@@ -628,6 +658,16 @@ public:
                 CompareOp& cop
             ) const;
 
+            //- Return a list containing the nearest distance of nodes to any
+            //  shapes
+            List<scalar> calcNearestDistance() const;
+
+            bool quickCircumsphereRejection
+            (
+                const point& cc,
+                const scalar crSqr,
+                const List<scalar>& nearestDistances
+            ) const;
 
         // Write
 
diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
index 0ebb7b58ab..81f853aa41 100644
--- a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
@@ -546,6 +546,49 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
     }
 
     featurePointNormals_ = featurePointNormals;
+
+
+    // Create featurePointEdges_
+    // For each feature point, stores a list of edges which are arranged in
+    // order according to their connectivity
+
+//    featurePointEdges_.setSize(nonFeatureStart_);
+//
+//    Info<< sFeatEds.size() << " " << surf.pointEdges().size() << " "
+//        << edges().size() << endl;
+//
+//    forAll(sFeatEds, eI)
+//    {
+//        Info<< "Edge " << eI << " " << sFeatEds[eI] << " "
+//            << edges()[eI] << endl;
+//    }
+//
+//    const edgeList& edges = eds;
+//
+//    for (label i = 0; i < nonFeatureStart_; i++)
+//    {
+//        const labelList& ptEds = surf.pointEdges()[i];
+//
+//        DynamicList<label> tmpFtEdges;
+//
+//        forAll(ptEds, eI)
+//        {
+//            const label edgeI = ptEds[eI];
+//            const edge& e = sFeatEds[edgeI];
+//
+//            Info<< "Edges: " << edgeI << " " << e << endl;
+//
+//            forAll(sFeatEds, fEdgeI)
+//            {
+//                if (edges[fEdgeI] == e)
+//                {
+//                    tmpFtEdges.append(fEdgeI);
+//                }
+//            }
+//        }
+//
+//        featurePointEdges_[i] = tmpFtEdges;
+//    }
 }
 
 
@@ -821,7 +864,7 @@ void Foam::extendedFeatureEdgeMesh::allNearestFeatureEdges
 
             if (!hitPoint.hit())
             {
-                nearHit = pointIndexHit(true, hitPoint.missPoint(), hitIndex);
+                nearHit = pointIndexHit(false, hitPoint.missPoint(), hitIndex);
             }
             else
             {