diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3Ddefs.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3Ddefs.H
index 8efb9ec711..3176e4dec4 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3Ddefs.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3Ddefs.H
@@ -72,6 +72,48 @@ typedef Delaunay::Cell_handle      Cell_handle;
 typedef Delaunay::Point            Point;
 
 
+//- Spatial sort traits to use with a pair of point pointers and an integer.
+//  Taken from a post on the CGAL lists: 2010-01/msg00004.html by
+//  Sebastien Loriot (Geometry Factory).
+template<class Triangulation>
+struct Traits_for_spatial_sort
+:
+    public Triangulation::Geom_traits
+{
+    typedef typename Triangulation::Geom_traits Gt;
+
+    typedef std::pair<const typename Triangulation::Point*, int> Point_3;
+
+    struct Less_x_3
+    {
+        bool operator()(const Point_3& p, const Point_3& q) const
+        {
+            return typename Gt::Less_x_3()(*(p.first), *(q.first));
+        }
+    };
+
+    struct Less_y_3
+    {
+        bool operator()(const Point_3& p, const Point_3& q) const
+        {
+            return typename Gt::Less_y_3()(*(p.first), *(q.first));
+        }
+    };
+
+    struct Less_z_3
+    {
+        bool operator()(const Point_3& p, const Point_3& q) const
+        {
+            return typename Gt::Less_z_3()(*(p.first), *(q.first));
+        }
+    };
+
+    Less_x_3  less_x_3_object () const {return Less_x_3();}
+    Less_y_3  less_y_3_object () const {return Less_y_3();}
+    Less_z_3  less_z_3_object () const {return Less_z_3();}
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 58d061b8f7..f9516c09a7 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -316,16 +316,16 @@ void Foam::conformalVoronoiMesh::insertPoints
     // using the range insert (faster than inserting points one by one)
     insert(points.begin(), points.end());
 
-    // Info<< "USING INDIVIDUAL INSERTION TO DETECT FAILURE" << endl;
-    // for
-    // (
-    //List:list<Point>::iterator pit=points.begin();
-    //     pit != points.end();
-    //     ++pit
-    // )
-    // {
-    //     insertPoint(topoint(*pit), Vb::vtInternal);
-    // }
+//     Info<< "USING INDIVIDUAL INSERTION TO DETECT FAILURE" << endl;
+//     for
+//     (
+//         List<Point>::iterator pit=points.begin();
+//         pit != points.end();
+//         ++pit
+//     )
+//     {
+//         insertPoint(topoint(*pit), Vb::vtInternal);
+//     }
 
     label nInserted(number_of_vertices() - preInsertionSize);
 
@@ -395,31 +395,14 @@ void Foam::conformalVoronoiMesh::insertPoints
         //     << " points in total" << endl;
     }
 
-    // Using index is actually pointless, it is always zero.  Keep for clarity
-    // of code.
-
-    forAll(pts, pI)
-    {
-        // creation of points and indices is done assuming that it will be
-        // relative to the instantaneous number_of_vertices() at insertion.
-
-        label type = types[pI];
-
-        if (type > Vb::vtFar)
-        {
-            // This is a member of a point pair, don't use the type directly
-            // (note that this routine never gets called for referredPoints
-            //  so type will never be -procI)
-            type += number_of_vertices();
-        }
-
-        insertPoint
-        (
-            pts[pI],
-            indices[pI] + number_of_vertices(),
-            type
-        );
-    }
+    rangeInsertWithInfo
+    (
+        pts.begin(),
+        pts.end(),
+        *this,
+        indices,
+        types
+    );
 }
 
 
@@ -1202,6 +1185,8 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
     featureVertices_(),
     featurePointLocations_(),
     featurePointTreePtr_(),
+    edgeLocationTreePtr_(),
+    surfacePtLocationTreePtr_(),
     sizeAndAlignmentLocations_(),
     storedSizes_(),
     storedAlignments_(),
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 2e38a504e1..a51e5a595f 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -48,6 +48,7 @@ SourceFiles
 #define CGAL_INEXACT
 
 #include "CGALTriangulation3Ddefs.H"
+#include <CGAL/Spatial_sort_traits_adapter_3.h>
 #include "uint.H"
 #include "ulong.H"
 #include "searchableSurfaces.H"
@@ -164,6 +165,14 @@ private:
         //- Search tree for feature point locations
         mutable autoPtr<indexedOctree<treeDataPoint> > featurePointTreePtr_;
 
+        //- Search tree for edge point locations
+        mutable autoPtr<dynamicIndexedOctree<dynamicTreeDataPoint> >
+        edgeLocationTreePtr_;
+
+        //- Search tree for surface point locations
+        mutable autoPtr<dynamicIndexedOctree<dynamicTreeDataPoint> >
+        surfacePtLocationTreePtr_;
+
         //- Store locations where the cell size and alignments will be
         //  pre-calculated and looked up
         pointField sizeAndAlignmentLocations_;
@@ -244,6 +253,12 @@ private:
             const Foam::point& pt
         ) const;
 
+        //- Return the square of the local surface point exclusion distance
+        inline scalar surfacePtExclusionDistanceSqr
+        (
+            const Foam::point& pt
+        ) const;
+
         //- Return the square of the local surface search distance
         inline scalar surfaceSearchDistanceSqr(const Foam::point& pt) const;
 
@@ -664,9 +679,41 @@ private:
         //- Check if a surface point is near another.
         bool nearSurfacePoint
         (
-            const pointIndexHit& pHit,
-            DynamicList<Foam::point>& existingSurfacePtLocations,
-            dynamicIndexedOctree<dynamicTreeDataPoint>& surfacePtLocationTree
+            pointIndexHit& pHit,
+            label& surfaceHit,
+            DynamicList<Foam::point>& existingSurfacePtLocations
+        ) const;
+
+        //- Append a point to the surface point tree and the existing list
+        bool appendToSurfacePtTree
+        (
+            const Foam::point& pt,
+            DynamicList<Foam::point>& existingSurfacePtLocations
+        ) const;
+
+        //- Append a point to the edge location tree and the existing list
+        bool appendToEdgeLocationTree
+        (
+            const Foam::point& pt,
+            DynamicList<Foam::point>& existingEdgeLocations
+        ) const;
+
+        //- Check if a point is near any feature edge points.
+        bool pointIsNearFeatureEdge(const Foam::point& pt) const;
+
+        bool pointIsNearFeatureEdge
+        (
+            const Foam::point& pt,
+            pointIndexHit& info
+        ) const;
+
+        //- Check if a point is near any surface conformation points.
+        bool pointIsNearSurfaceLocation(const Foam::point& pt) const;
+
+        bool pointIsNearSurfaceLocation
+        (
+            const Foam::point& pt,
+            pointIndexHit& info
         ) const;
 
         //- Check if a location is in the exclusion range of an existing feature
@@ -674,17 +721,20 @@ private:
         bool nearFeatureEdgeLocation
         (
             pointIndexHit& pHit,
-            DynamicList<Foam::point>& newEdgeLocations,
-            DynamicList<Foam::point>& existingEdgeLocations,
-            dynamicIndexedOctree<dynamicTreeDataPoint>& edgeLocationTree
+            DynamicList<Foam::point>& existingEdgeLocations
         ) const;
 
-        //- Build or rebuild the edgeLocationTree
-//        void buildEdgeLocationTree
-//        (
-//            autoPtr<dynamicIndexedOctree<dynamicTreeDataPoint> >& edgeLocationTree,
-//            const pointField& existingEdgeLocations
-//        ) const;
+        //- Build or rebuild the edge location tree
+        void buildEdgeLocationTree
+        (
+            const DynamicList<Foam::point>& existingEdgeLocations
+        ) const;
+
+        //- Build or rebuild the surface point location tree
+        void buildSurfacePtLocationTree
+        (
+            const DynamicList<Foam::point>& existingSurfacePtLocations
+        ) const;
 
         //- Build or rebuild the sizeAndAlignmentTree
         void buildSizeAndAlignmentTree() const;
@@ -705,9 +755,7 @@ private:
             DynamicList<label>& featureEdgeFeaturesHit,
             DynamicList<Foam::point>& newEdgeLocations,
             DynamicList<Foam::point>& existingEdgeLocations,
-            dynamicIndexedOctree<dynamicTreeDataPoint>& edgeLocationTree,
-            DynamicList<Foam::point>& existingSurfacePtLocations,
-            dynamicIndexedOctree<dynamicTreeDataPoint>& surfacePtLocationTree
+            DynamicList<Foam::point>& existingSurfacePtLocations
         ) const;
 
         //- Store the surface conformation with the indices offset to be
@@ -1101,6 +1149,87 @@ public:
             //- Find the cellSet of the boundary cells which have points that
             //  protrude out of the surface beyond a tolerance.
             void findRemainingProtrusionSet(const fvMesh& mesh) const;
+
+
+            //- Function inserting points into a triangulation and setting the
+            //  index and type data of the point in the correct order. This is
+            //  faster than inserting points individually.
+            //
+            //  Adapted from a post on the CGAL lists: 2010-01/msg00004.html by
+            //  Sebastien Loriot (Geometry Factory).
+            //
+            //  @todo problems putting declaration in the .C file. Function
+            //        prototype is not recognised.
+            template<class Triangulation, class Point_iterator>
+            void rangeInsertWithInfo
+            (
+                Point_iterator begin,
+                Point_iterator end,
+                Triangulation& T,
+                DynamicList<label>& indices,
+                DynamicList<label>& types
+            )
+            {
+                typedef std::vector
+                <
+                    std::pair<const typename Triangulation::Point*, label>
+                > vectorPairPointIndex;
+
+                vectorPairPointIndex points;
+                label index = 0;
+
+                for (Point_iterator it = begin; it != end; ++it)
+                {
+                    points.push_back
+                    (
+                        std::make_pair(&(toPoint(*it)), index++)
+                    );
+                }
+
+                std::random_shuffle(points.begin(), points.end());
+
+                spatial_sort
+                (
+                    points.begin(),
+                    points.end(),
+                    Traits_for_spatial_sort<Triangulation>()
+                );
+
+                typename Triangulation::Cell_handle hint;
+
+                for
+                (
+                    typename vectorPairPointIndex::const_iterator
+                    p = points.begin();
+                    p != points.end();
+                    ++p
+                )
+                {
+                    typename Triangulation::Locate_type lt;
+                    typename Triangulation::Cell_handle c;
+                    label li, lj;
+
+                    c = T.locate(*(p->first), lt, li, lj, hint);
+
+                    typename Triangulation::Vertex_handle v
+                        = T.insert(*(p->first), lt, c, li, lj);
+
+                    label oldIndex = p->second;
+
+                    label type = types[oldIndex];
+
+                    if (type > Vb::vtFar)
+                    {
+                        // This is a member of a point pair, don't use the type
+                        // directly (note that this routine never gets called
+                        // for referredPoints so type will never be -procI)
+                        type += T.number_of_vertices() - 1;
+                    }
+
+                    v->index() = indices[oldIndex] + T.number_of_vertices() - 1;
+                    v->type() = type;
+                }
+            }
 };
 
 
@@ -1112,6 +1241,10 @@ public:
 
 #include "conformalVoronoiMeshI.H"
 
+//#ifdef NoRepository
+//#   include "conformalVoronoiMeshTemplates.C"
+//#endif
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index 7d8ec2619c..20ae246d27 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -127,31 +127,13 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
 
     // Initialise containers to store the edge conformation locations
     DynamicList<Foam::point> newEdgeLocations;
+
     DynamicList<Foam::point> existingEdgeLocations;
     DynamicList<Foam::point> existingSurfacePtLocations;
 
-    treeBoundBox overallBb
-    (
-        geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4)
-    );
+    buildEdgeLocationTree(existingEdgeLocations);
+    buildSurfacePtLocationTree(existingSurfacePtLocations);
 
-    dynamicIndexedOctree<dynamicTreeDataPoint> edgeLocationTree
-    (
-        dynamicTreeDataPoint(existingEdgeLocations),
-        overallBb,  // overall search domain
-        10,         // max levels; n/a for dynamic tree
-        100.0,      // maximum ratio of cubes v.s. cells
-        100.0       // max. duplicity; n/a since no bounding boxes.
-    );
-
-    dynamicIndexedOctree<dynamicTreeDataPoint> surfacePtLocationTree
-	(
-		dynamicTreeDataPoint(existingSurfacePtLocations),
-		overallBb,  // overall search domain
-		10,         // max levels; n/a for dynamic tree
-		100.0,      // maximum ratio of cubes v.s. cells
-		100.0       // max. duplicity; n/a since no bounding boxes.
-	);
 
     // Initialise the edgeLocationTree
     //buildEdgeLocationTree(edgeLocationTree, existingEdgeLocations);
@@ -255,9 +237,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                         featureEdgeFeaturesHit,
                         newEdgeLocations,
                         existingEdgeLocations,
-                        edgeLocationTree,
-                        existingSurfacePtLocations,
-                        surfacePtLocationTree
+                        existingSurfacePtLocations
                     );
                 }
             }
@@ -474,9 +454,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                         featureEdgeFeaturesHit,
                         newEdgeLocations,
                         existingEdgeLocations,
-                        edgeLocationTree,
-                        existingSurfacePtLocations,
-                        surfacePtLocationTree
+                        existingSurfacePtLocations
                     );
                 }
             }
@@ -513,9 +491,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                         featureEdgeFeaturesHit,
                         newEdgeLocations,
                         existingEdgeLocations,
-                        edgeLocationTree,
-                        existingSurfacePtLocations,
-                        surfacePtLocationTree
+                        existingSurfacePtLocations
                     );
                 }
             }
@@ -1862,7 +1838,7 @@ Foam::scalar Foam::conformalVoronoiMesh::angleBetweenSurfacePoints
     pointIndexHit pAhit;
     label pAsurfaceHit = -1;
 
-    const scalar searchDist = 1.0*targetCellSize(pA);
+    const scalar searchDist = 5.0*targetCellSize(pA);
 
     geometryToConformTo_.findSurfaceNearest
     (
@@ -1907,38 +1883,172 @@ Foam::scalar Foam::conformalVoronoiMesh::angleBetweenSurfacePoints
 
 bool Foam::conformalVoronoiMesh::nearSurfacePoint
 (
-    const pointIndexHit& pHit,
-    DynamicList<Foam::point>& existingSurfacePtLocations,
-    dynamicIndexedOctree<dynamicTreeDataPoint>& surfacePtLocationTree
+    pointIndexHit& pHit,
+    label& surfaceHit,
+    DynamicList<Foam::point>& existingSurfacePtLocations
 ) const
 {
-    const Foam::point pt = pHit.hitPoint();
+    const Foam::point& pt = pHit.hitPoint();
 
-    scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
+    pointIndexHit closePoint;
 
-    pointIndexHit info = surfacePtLocationTree.findNearest(pt, exclusionRangeSqr);
+    const bool closeToSurfacePt = pointIsNearSurfaceLocation(pt, closePoint);
 
-
-    // Add the point to the surface point tree if it will be added.
-    label startIndex = existingSurfacePtLocations.size();
-
-    existingSurfacePtLocations.append(pHit.hitPoint());
-
-    label endIndex = existingSurfacePtLocations.size();
-
-    surfacePtLocationTree.insert(startIndex, endIndex);
-
-
-    if (info.hit())
+    if (closeToSurfacePt)
     {
-        const scalar angle = angleBetweenSurfacePoints(pt, info.hitPoint());
+        const scalar angle
+            = angleBetweenSurfacePoints(pt, closePoint.hitPoint());
 
+        // @todo make this tolerance run-time selectable?
         if (angle > 170.0)
         {
-            return false;
+            pointIndexHit pCloseHit;
+            label pCloseSurfaceHit = -1;
+
+            const scalar searchDist = targetCellSize(closePoint.hitPoint());
+
+            geometryToConformTo_.findSurfaceNearest
+            (
+                closePoint.hitPoint(),
+                searchDist,
+                pCloseHit,
+                pCloseSurfaceHit
+            );
+
+            vectorField norm(1);
+
+            allGeometry_[pCloseSurfaceHit].getNormal
+            (
+                List<pointIndexHit>(1, pCloseHit),
+                norm
+            );
+
+            const vector nA = norm[0];
+
+            pointIndexHit oppositeHit;
+            label oppositeSurfaceHit = -1;
+
+            geometryToConformTo_.findSurfaceNearestIntersection
+            (
+                closePoint.hitPoint() + SMALL*nA,
+                closePoint.hitPoint() + mag(pt - closePoint.hitPoint())*nA,
+                oppositeHit,
+                oppositeSurfaceHit
+            );
+
+            if (oppositeHit.hit())
+            {
+                // Replace point
+                pHit = oppositeHit;
+                surfaceHit = oppositeSurfaceHit;
+
+                Foam::point newPt = oppositeHit.hitPoint();
+
+                appendToSurfacePtTree(newPt, existingSurfacePtLocations);
+
+                return !closeToSurfacePt;
+
+                if (debug)
+                {
+                    Info<< "Point " << pt
+                        << " is close to " << closePoint.hitPoint()
+                        << " so will be moved to " << oppositeHit.hitPoint()
+                        << endl;
+                }
+            }
         }
     }
+    else
+    {
+        appendToSurfacePtTree(pt, existingSurfacePtLocations);
+    }
 
+    return closeToSurfacePt;
+}
+
+
+bool Foam::conformalVoronoiMesh::appendToSurfacePtTree
+(
+    const Foam::point& pt,
+    DynamicList<Foam::point>& existingSurfacePtLocations
+) const
+{
+   label startIndex = existingSurfacePtLocations.size();
+
+   existingSurfacePtLocations.append(pt);
+
+   label endIndex = existingSurfacePtLocations.size();
+
+   return surfacePtLocationTreePtr_().insert(startIndex, endIndex);
+}
+
+
+bool Foam::conformalVoronoiMesh::appendToEdgeLocationTree
+(
+    const Foam::point& pt,
+    DynamicList<Foam::point>& existingEdgeLocations
+) const
+{
+   label startIndex = existingEdgeLocations.size();
+
+   existingEdgeLocations.append(pt);
+
+   label endIndex = existingEdgeLocations.size();
+
+   return edgeLocationTreePtr_().insert(startIndex, endIndex);
+}
+
+
+bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdge
+(
+    const Foam::point& pt
+) const
+{
+    const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
+
+    pointIndexHit info
+        = edgeLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
+
+    return info.hit();
+}
+
+
+bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdge
+(
+    const Foam::point& pt,
+    pointIndexHit& info
+) const
+{
+    const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
+
+    info = edgeLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
+
+    return info.hit();
+}
+
+
+bool Foam::conformalVoronoiMesh::pointIsNearSurfaceLocation
+(
+    const Foam::point& pt
+) const
+{
+    pointIndexHit info;
+
+    pointIsNearSurfaceLocation(pt, info);
+
+    return info.hit();
+}
+
+
+bool Foam::conformalVoronoiMesh::pointIsNearSurfaceLocation
+(
+    const Foam::point& pt,
+    pointIndexHit& info
+) const
+{
+    const scalar exclusionRangeSqr = surfacePtExclusionDistanceSqr(pt);
+
+    info = surfacePtLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
 
     return info.hit();
 }
@@ -1947,26 +2057,60 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
 bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
 (
     pointIndexHit& pHit,
-    DynamicList<Foam::point>& newEdgeLocations,
-    DynamicList<Foam::point>& existingEdgeLocations,
-    dynamicIndexedOctree<dynamicTreeDataPoint>& edgeLocationTree
+    DynamicList<Foam::point>& existingEdgeLocations
 ) const
 {
     const Foam::point pt = pHit.hitPoint();
 
-    scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
+    const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
 
-    label startIndex = existingEdgeLocations.size();
+    pointIndexHit info;
 
-    existingEdgeLocations.append(newEdgeLocations);
+    bool closeToFeatureEdge = pointIsNearFeatureEdge(pt, info);
 
-    label endIndex = existingEdgeLocations.size();
+    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;
 
-    edgeLocationTree.insert(startIndex, endIndex);
+        geometryToConformTo_.findEdgeNearest
+        (
+            pt,
+            exclusionRangeSqr,
+            edgeHit,
+            featureHit
+        );
 
-    newEdgeLocations.clear();
+        const extendedFeatureEdgeMesh& eMesh
+            = geometryToConformTo_.features()[featureHit];
 
-    pointIndexHit info = edgeLocationTree.findNearest(pt, exclusionRangeSqr);
+        const vector& edgeDir = eMesh.edgeDirections()[edgeHit.index()];
+
+        const vector lineBetweenPoints = pt - info.hitPoint();
+
+        const scalar angle = degAngleBetween(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.
+        if
+        (
+            angle < 100.0
+         && angle > 80.0
+         && mag(lineBetweenPoints) > SMALL
+        )
+        {
+            closeToFeatureEdge = false;
+            appendToEdgeLocationTree(pt, existingEdgeLocations);
+        }
+    }
+
+    return closeToFeatureEdge;
 
     // Searching for the nearest point in existingEdgeLocations using the
     // indexedOctree
@@ -2009,36 +2153,61 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
 //        return !info.hit();
 //    }
 
-    return info.hit();
 }
 
 
-//void Foam::conformalVoronoiMesh::buildEdgeLocationTree
-//(
-//    autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree,
-//    DynamicList<Foam::point>& existingEdgeLocations
-//) const
-//{
-//    treeBoundBox overallBb
-//    (
-//        geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4)
-//    );
-//
-//    overallBb.min() -= Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-//    overallBb.max() += Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-//
-//    edgeLocationTree.reset
-//    (
-//        new indexedOctree<treeDataPoint>
-//        (
-//            treeDataPoint(existingEdgeLocations),
-//            overallBb,  // overall search domain
-//            10,         // max levels
-//            10.0,       // maximum ratio of cubes v.s. cells
-//            100.0       // max. duplicity; n/a since no bounding boxes.
-//        )
-//    );
-//}
+void Foam::conformalVoronoiMesh::buildEdgeLocationTree
+(
+    const DynamicList<Foam::point>& existingEdgeLocations
+) const
+{
+    treeBoundBox overallBb
+    (
+        geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4)
+    );
+
+    overallBb.min() -= Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    overallBb.max() += Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+    edgeLocationTreePtr_.reset
+    (
+        new dynamicIndexedOctree<dynamicTreeDataPoint>
+        (
+            dynamicTreeDataPoint(existingEdgeLocations),
+            overallBb,  // overall search domain
+            10,         // max levels, n/a
+            20.0,       // maximum ratio of cubes v.s. cells
+            100.0       // max. duplicity; n/a since no bounding boxes.
+        )
+    );
+}
+
+
+void Foam::conformalVoronoiMesh::buildSurfacePtLocationTree
+(
+    const DynamicList<Foam::point>& existingSurfacePtLocations
+) const
+{
+    treeBoundBox overallBb
+    (
+        geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4)
+    );
+
+    overallBb.min() -= Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    overallBb.max() += Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+    surfacePtLocationTreePtr_.reset
+    (
+        new dynamicIndexedOctree<dynamicTreeDataPoint>
+        (
+            dynamicTreeDataPoint(existingSurfacePtLocations),
+            overallBb,  // overall search domain
+            10,         // max levels, n/a
+            20.0,       // maximum ratio of cubes v.s. cells
+            100.0       // max. duplicity; n/a since no bounding boxes.
+        )
+    );
+}
 
 
 void Foam::conformalVoronoiMesh::buildSizeAndAlignmentTree() const
@@ -2087,9 +2256,7 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
     DynamicList<label>& featureEdgeFeaturesHit,
     DynamicList<Foam::point>& newEdgeLocations,
     DynamicList<Foam::point>& existingEdgeLocations,
-    dynamicIndexedOctree<dynamicTreeDataPoint>& edgeLocationTree,
-    DynamicList<Foam::point>& existingSurfacePtLocations,
-    dynamicIndexedOctree<dynamicTreeDataPoint>& surfacePtLocationTree
+    DynamicList<Foam::point>& existingSurfacePtLocations
 ) const
 {
     bool keepSurfacePoint = true;
@@ -2099,8 +2266,8 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
 
     forAll(surfHit, sI)
     {
-        const pointIndexHit& surfHitI = surfHit[sI];
-        const label hitSurfaceI = hitSurface[sI];
+        pointIndexHit surfHitI = surfHit[sI];
+        label hitSurfaceI = hitSurface[sI];
 
         if (!surfHitI.hit())
         {
@@ -2112,7 +2279,8 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
             keepSurfacePoint = false;
 
             // If the triggering Vertex is part of a feature point, allow it to
-            // conform to the surface
+            // conform to the surface.
+            // @todo Is this needed?! Shouldn't be any feature points here...
             if (vit->index() < startOfInternalPoints_)
             {
                 surfaceHits.append(surfHitI);
@@ -2126,8 +2294,8 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
             nearSurfacePoint
             (
                 surfHitI,
-                existingSurfacePtLocations,
-                surfacePtLocationTree
+                hitSurfaceI,
+                existingSurfacePtLocations
             )
         )
         {
@@ -2184,9 +2352,7 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
                         !nearFeatureEdgeLocation
                         (
                             edHit,
-                            newEdgeLocations,
-                            existingEdgeLocations,
-                            edgeLocationTree
+                            existingEdgeLocations
                         )
                     )
                     {
@@ -2202,8 +2368,6 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
             }
         }
 
-        newEdgeLocations.append(currentEdgeLocations);
-
         if (keepSurfacePoint)
         {
             surfaceHits.append(surfHitI);
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index 0cad8731f7..a18b173ecd 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -186,6 +186,20 @@ inline Foam::scalar Foam::conformalVoronoiMesh::featureEdgeExclusionDistanceSqr
 }
 
 
+inline Foam::scalar Foam::conformalVoronoiMesh::surfacePtExclusionDistanceSqr
+(
+    const Foam::point& pt
+) const
+{
+    return
+        sqr
+        (
+            targetCellSize(pt)
+           *cvMeshControls().surfacePtExclusionDistanceCoeff()
+        );
+}
+
+
 inline Foam::scalar Foam::conformalVoronoiMesh::surfaceSearchDistanceSqr
 (
     const Foam::point& pt
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
index 1058b6f133..9b59ce9b9b 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
@@ -99,6 +99,11 @@ Foam::cvControls::cvControls
         coarseDict.subDict("iteration")
     );
 
+    surfacePtExclusionDistanceCoeff_ = readScalar
+    (
+        coarseInitialDict.lookup("surfacePtExclusionDistanceCoeff")
+    );
+
     edgeSearchDistCoeffSqr_coarse_initial_ = sqr
     (
         readScalar
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
index 28e92c05bd..f6bdb75559 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
@@ -99,6 +99,11 @@ class cvControls
 
             // Controls for coarse surface conformation
 
+                //- Distance to an existing surface conformation point location
+                //  within which other surface point locations are excluded
+                //  - fraction of the local target cell size
+                scalar surfacePtExclusionDistanceCoeff_;
+
                 //- Distance to search for feature edges near to
                 //  surface protrusions - fraction of the local target
                 //  cell size. Coarse conformation, initial protrusion
@@ -318,6 +323,9 @@ public:
             //- Return the featureEdgeExclusionDistanceCoeff
             inline scalar featureEdgeExclusionDistanceCoeff() const;
 
+            //- Return the surfacePtExclusionDistanceCoeff
+            inline scalar surfacePtExclusionDistanceCoeff() const;
+
             //- Return whether to use specialised feature points
             inline Switch specialiseFeaturePoints() const;
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
index 53e79ca1cd..dc3f391100 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
@@ -54,6 +54,11 @@ inline Foam::scalar Foam::cvControls::featureEdgeExclusionDistanceCoeff() const
     return featureEdgeExclusionDistanceCoeff_;
 }
 
+inline Foam::scalar Foam::cvControls::surfacePtExclusionDistanceCoeff() const
+{
+    return surfacePtExclusionDistanceCoeff_;
+}
+
 inline Foam::Switch Foam::cvControls::specialiseFeaturePoints() const
 {
     return specialiseFeaturePoints_;