From 7f8af1c4330d2e6cc7fe60952b18360f8e4b2522 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Fri, 1 May 2009 18:47:13 +0100
Subject: [PATCH] Changed findSurfaceNearest to return surface index instead of
 normals directly.

Added searching for largest normal protrusion to surface conformation within an
iterative loop after the initial, closest point on surface conformation.
---
 .../conformalVoronoiMesh.C                    | 339 ++++++++++++++----
 .../conformalVoronoiMesh.H                    |  22 +-
 .../conformalVoronoiMeshI.H                   |   9 +
 .../conformationSurfaces.C                    |  87 +++--
 .../conformationSurfaces.H                    |  20 +-
 .../cvControls/cvControls.C                   |  11 +-
 .../cvControls/cvControls.H                   |   8 +
 .../cvControls/cvControlsI.H                  |   6 +
 .../initialPointsMethod/pointFile/pointFile.C |   2 +
 9 files changed, 381 insertions(+), 123 deletions(-)

diff --git a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index e23ce49130..590647af0b 100644
--- a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -46,7 +46,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
     (
         IOobject
         (
-            "cvSearchableSurfacesDirectory",
+            "cvSearchableSurfaces",
             runTime_.constant(),
             "triSurface",
             runTime_,
@@ -120,40 +120,52 @@ Foam::conformalVoronoiMesh::~conformalVoronoiMesh()
 
 void Foam::conformalVoronoiMesh::insertSurfacePointPairs
 (
-    const List<scalar>& surfacePpDist,
-    const List<point>& surfacePoints,
-    const List<vector>& surfaceNormals,
+    const List<pointIndexHit>& surfaceHits,
+    const List<label>& hitSurfaces,
     const fileName fName
 )
 {
-    if
-    (
-        surfacePpDist.size() != surfacePoints.size()
-     || surfacePpDist.size() != surfaceNormals.size()
-    )
+    if (surfaceHits.size() != hitSurfaces.size())
     {
         FatalErrorIn("Foam::conformalVoronoiMesh::insertPointPairs")
-            << "surfacePpDist, surfacePoints and surfaceNormals are not "
-            << "the same size. Sizes "
-            << surfacePpDist.size() << ' '
-            << surfacePoints.size() << ' '
-            << surfaceNormals.size()
+            << "surfaceHits and hitSurfaces are not the same size. Sizes "
+            << surfaceHits.size() << ' '
+            << hitSurfaces.size()
             << exit(FatalError);
     }
 
-    forAll(surfacePoints, p)
+    forAll(surfaceHits, i)
     {
+        vectorField norm(1);
+
+        allGeometry_[hitSurfaces[i]].getNormal
+        (
+            List<pointIndexHit>(1, surfaceHits[i]),
+            norm
+        );
+
+        const vector& normal = norm[0];
+
+        const point& surfacePt(surfaceHits[i].hitPoint());
+
         insertPointPair
         (
-            surfacePpDist[p],
-            surfacePoints[p],
-            surfaceNormals[p]
+            pointPairDistance(surfacePt),
+            surfacePt,
+            normal
         );
     }
 
     if (fName != fileName::null)
     {
-        writePoints(fName, surfacePoints);
+        List<point> surfacePts(surfaceHits.size());
+
+        forAll(surfacePts, i)
+        {
+            surfacePts[i] = surfaceHits[i].hitPoint();
+        }
+
+        writePoints(fName, surfacePts);
     }
 }
 
@@ -161,7 +173,8 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
 void Foam::conformalVoronoiMesh::insertEdgePointGroups
 (
     const List<pointIndexHit>& edgeHits,
-    const labelList& featuresHit
+    const labelList& featuresHit,
+    const fileName fName
 )
 {
     if (edgeHits.size() != featuresHit.size())
@@ -205,6 +218,18 @@ void Foam::conformalVoronoiMesh::insertEdgePointGroups
             insertMultipleEdgePointGroup(feMesh, edgeHits[i]);
         }
     }
+
+    if (fName != fileName::null)
+    {
+        List<point> edgePts(edgeHits.size());
+
+        forAll(edgePts, i)
+        {
+            edgePts[i] = edgeHits[i].hitPoint();
+        }
+
+        writePoints(fName, edgePts);
+    }
 }
 
 
@@ -225,8 +250,10 @@ void Foam::conformalVoronoiMesh::insertExternalEdgePointGroup
     const vector& nA = feNormals[edNormalIs[0]];
     const vector& nB = feNormals[edNormalIs[1]];
 
+
+
     // Convex. So refPt will be inside domain and hence a master point
-    point refPt = edgePt - ppDist*(nA + nB)/(1 + (nA & nB));
+    point refPt = edgePt - ppDist*(nA + nB)/(1 + (nA & nB) + VSMALL);
 
     // Insert the master point referring the the first slave
     label masterPtIndex = insertPoint(refPt, number_of_vertices() + 1);
@@ -260,7 +287,7 @@ void Foam::conformalVoronoiMesh::insertInternalEdgePointGroup
     const vector& nB = feNormals[edNormalIs[1]];
 
     // Concave. master and reflected points inside the domain.
-    point refPt = edgePt - ppDist*(nA + nB)/(1 + (nA & nB));
+    point refPt = edgePt - ppDist*(nA + nB)/(1 + (nA & nB) + VSMALL );
 
     // Generate reflected master to be outside.
     point reflMasterPt = edgePt + 2*(edgePt - refPt);
@@ -351,7 +378,7 @@ void Foam::conformalVoronoiMesh::conformToFeaturePoints()
 
     insertMixedFeaturePoints();
 
-    Info<< "   Inserted " << number_of_vertices() << " vertices" << endl;
+    Info<< "    Inserted " << number_of_vertices() << " vertices" << endl;
 
     featureVertices_.setSize(number_of_vertices());
 
@@ -516,12 +543,12 @@ void Foam::conformalVoronoiMesh::insertInitialPoints()
 }
 
 
-bool Foam::conformalVoronoiMesh::dualCellSurfaceIntersection
+bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
 (
     const Triangulation::Finite_vertices_iterator& vit
 ) const
 {
-    std::list<Facet>  facets;
+    std::list<Facet> facets;
     incident_facets(vit, std::back_inserter(facets));
 
     for
@@ -534,7 +561,7 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceIntersection
         if
         (
             is_infinite(fit->first)
-            || is_infinite(fit->first->neighbor(fit->second))
+         || is_infinite(fit->first->neighbor(fit->second))
         )
         {
             return true;
@@ -567,6 +594,79 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceIntersection
 }
 
 
+void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
+(
+    const Triangulation::Finite_vertices_iterator& vit,
+    pointIndexHit& surfHitLargest,
+    label& hitSurfaceLargest
+) const
+{
+    std::list<Facet> facets;
+    incident_facets(vit, std::back_inserter(facets));
+
+    for
+    (
+        std::list<Facet>::iterator fit=facets.begin();
+        fit != facets.end();
+        ++fit
+    )
+    {
+        if
+        (
+            !is_infinite(fit->first)
+         && !is_infinite(fit->first->neighbor(fit->second))
+        )
+        {
+            point vert(topoint(vit->point()));
+
+            scalar maxProtrusionDistance = maxSurfaceProtrusion(vert);
+
+            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
+           );
+
+            if (surfHit.hit())
+            {
+                vectorField norm(1);
+
+                allGeometry_[hitSurface].getNormal
+                (
+                    List<pointIndexHit>(1, surfHit),
+                    norm
+                );
+
+                const vector& n = norm[0];
+
+                scalar normalProtrusionDistance =
+                    (edgeMid - surfHit.hitPoint()) & n;
+
+                if (normalProtrusionDistance > maxProtrusionDistance)
+                {
+                    surfHitLargest = surfHit;
+                    hitSurfaceLargest = hitSurface;
+
+                    maxProtrusionDistance = normalProtrusionDistance;
+                }
+            }
+        }
+    }
+}
+
+
 void Foam::conformalVoronoiMesh::calcDualMesh
 (
     pointField& points,
@@ -1203,24 +1303,17 @@ void Foam::conformalVoronoiMesh::conformToSurface()
 
     startOfSurfacePointPairs_ = number_of_vertices();
 
-    Info<< "    EDGE DISTANCE COEFFS HARD-CODED." << endl;
-    scalar edgeSearchDistCoeffSqr = sqr(1.25);
-    scalar surfacePtReplaceDistCoeffSqr = sqr(0.3);
-
-    // Surface protrusion conformation
-
-    label nIterations = 2;
-
-    for(label iterationNo = 0; iterationNo < nIterations; iterationNo++)
+    // Initial surface protrusion conformation - nearest surface point
     {
-        DynamicList<scalar> surfacePpDist;
-        DynamicList<scalar> surfaceSearchDistSqr;
-        DynamicList<point> surfacePoints;
-        DynamicList<vector> surfaceNormals;
+        Info<< "    EDGE DISTANCE COEFFS HARD-CODED." << endl;
+        scalar edgeSearchDistCoeffSqr = sqr(1.1);
+        scalar surfacePtReplaceDistCoeffSqr = sqr(0.5);
+
+        DynamicList<pointIndexHit> surfaceHits;
+        DynamicList<label> hitSurfaces;
 
         DynamicList<pointIndexHit> featureEdgeHits;
         DynamicList<label> featureEdgeFeaturesHit;
-        DynamicList<point> tmpFeatureEdgePoints;
 
         for
         (
@@ -1235,21 +1328,21 @@ void Foam::conformalVoronoiMesh::conformToSurface()
                 point vert(topoint(vit->point()));
                 scalar searchDistanceSqr = surfaceSearchDistanceSqr(vert);
                 pointIndexHit surfHit;
-                vector normal;
+                label hitSurface;
 
-                geometryToConformTo_.findSurfaceNearestAndNormal
+                geometryToConformTo_.findSurfaceNearest
                 (
                     vert,
                     searchDistanceSqr,
                     surfHit,
-                    normal
+                    hitSurface
                 );
 
                 if (surfHit.hit())
                 {
                     vit->setNearBoundary();
 
-                    if (dualCellSurfaceIntersection(vit))
+                    if (dualCellSurfaceAnyIntersection(vit))
                     {
                         // If the point is within a given distance of a feature
                         // edge, shift it to being an edge control point
@@ -1279,8 +1372,6 @@ void Foam::conformalVoronoiMesh::conformToSurface()
                             featureEdgeHits.append(edHit);
                             featureEdgeFeaturesHit.append(featureHit);
 
-                            tmpFeatureEdgePoints.append(edHit.hitPoint());
-
                             if
                             (
                                 magSqr(edHit.hitPoint() - surfHit.hitPoint())
@@ -1293,45 +1384,151 @@ void Foam::conformalVoronoiMesh::conformToSurface()
 
                         if (keepSurfacePoint)
                         {
-                            surfacePpDist.append(pointPairDistance(vert));
-                            surfaceSearchDistSqr.append(searchDistanceSqr);
-                            surfacePoints.append(surfHit.hitPoint());
-                            surfaceNormals.append(normal);
+                            surfaceHits.append(surfHit);
+                            hitSurfaces.append(hitSurface);
                         }
                     }
                 }
             }
         }
 
+        Info<< nl <<"    Initial conformation " << nl
+            << "    number_of_vertices " << number_of_vertices() << nl
+            << "    surfaceHits.size() " << surfaceHits.size() << nl
+            << "    featureEdgeHits.size() " << featureEdgeHits.size()
+            << endl;
+
+        insertSurfacePointPairs
+        (
+            surfaceHits,
+            hitSurfaces,
+            "surfaceConformationLocations_initial.obj"
+        );
+
+        insertEdgePointGroups
+        (
+            featureEdgeHits,
+            featureEdgeFeaturesHit,
+            "edgeGenerationLocations_initial.obj"
+        );
+    }
+
+    label iterationNo = 0;
+    label maxIterations = 4;
+    label totalHits = 0;
+
+    do
+    {
+        Info<< "    EDGE DISTANCE COEFFS HARD-CODED." << endl;
+        scalar edgeSearchDistCoeffSqr = sqr(1.25);
+        scalar surfacePtReplaceDistCoeffSqr = sqr(0.7);
+
+        DynamicList<pointIndexHit> surfaceHits;
+        DynamicList<label> hitSurfaces;
+
+        DynamicList<pointIndexHit> featureEdgeHits;
+        DynamicList<label> featureEdgeFeaturesHit;
+
+        for
+        (
+            Triangulation::Finite_vertices_iterator vit =
+            finite_vertices_begin();
+            vit != finite_vertices_end();
+            vit++
+        )
+        {
+            // The initial surface conformation has already identified the
+            // nearBoundary set of vertices.  Previously inserted boundary
+            // points can also generate protrusions and must be assessed too.
+
+            if (vit->nearBoundary() || vit->ppMaster())
+            {
+                point vert(topoint(vit->point()));
+                pointIndexHit surfHit;
+                label hitSurface;
+
+                dualCellLargestSurfaceProtrusion(vit, surfHit, hitSurface);
+
+                if (surfHit.hit())
+                {
+                    // If the point is within a given distance of a feature
+                    // edge, shift it to being an edge control point
+                    // instead, this will prevent "pits" forming.
+
+                    pointIndexHit edHit;
+                    label featureHit;
+
+                    scalar targetCellSizeSqr = sqr(targetCellSize(vert));
+
+                    geometryToConformTo_.findEdgeNearest
+                    (
+                        surfHit.hitPoint(),
+                        edgeSearchDistCoeffSqr*targetCellSizeSqr,
+                        edHit,
+                        featureHit
+                    );
+
+                    bool keepSurfacePoint = true;
+
+                    if (edHit.hit())
+                    {
+                        // Note that edge type classification can be used
+                        // here if necessary, i.e to avoid over-populating
+                        // internal edges
+
+                        featureEdgeHits.append(edHit);
+                        featureEdgeFeaturesHit.append(featureHit);
+
+                        if
+                        (
+                            magSqr(edHit.hitPoint() - surfHit.hitPoint())
+                          < surfacePtReplaceDistCoeffSqr*targetCellSizeSqr
+                        )
+                        {
+                            keepSurfacePoint = false;
+                        }
+                    }
+
+                    if (keepSurfacePoint)
+                    {
+                        surfaceHits.append(surfHit);
+                        hitSurfaces.append(hitSurface);
+                    }
+                }
+            }
+        }
+
         Info<< nl <<"    iterationNo " << iterationNo << nl
             << "    number_of_vertices " << number_of_vertices() << nl
-            << "    surfacePoints.size() " << surfacePoints.size() << nl
+            << "    surfaceHits.size() " << surfaceHits.size() << nl
             << "    featureEdgeHits.size() " << featureEdgeHits.size()
             << endl;
 
-        insertSurfacePointPairs
-        (
-            surfacePpDist,
-            surfacePoints,
-            surfaceNormals,
-            fileName
+        totalHits = surfaceHits.size() + featureEdgeHits.size();
+
+        if (totalHits > 0)
+        {
+            insertSurfacePointPairs
             (
-                "surfaceConformationLocations_" + name(iterationNo) + ".obj"
-            )
-        );
+                surfaceHits,
+                hitSurfaces,
+                fileName
+                (
+                    "surfaceConformationLocations_" + name(iterationNo) + ".obj"
+                )
+            );
 
-        insertEdgePointGroups
-        (
-            featureEdgeHits,
-            featureEdgeFeaturesHit
-        );
+            insertEdgePointGroups
+            (
+                featureEdgeHits,
+                featureEdgeFeaturesHit,
+                "edgeGenerationLocations_" + name(iterationNo) + ".obj"
+            );
+        }
 
-        writePoints
-        (
-            "edgeGenerationLocations_" + name(iterationNo) + ".obj",
-            tmpFeatureEdgePoints
-        );
-    }
+        iterationNo++;
+
+    } while (totalHits > 0 && iterationNo < maxIterations);
 }
 
 
diff --git a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index e524895f51..e5fe7b4ac3 100644
--- a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -125,6 +125,9 @@ class conformalVoronoiMesh
         //- Return the square of the local surface search distance
         inline scalar surfaceSearchDistanceSqr(const point& pt) const;
 
+        //- Return the local maximum surface protrusion distance
+        inline scalar maxSurfaceProtrusion(const point& pt) const;
+
         //- Return the local minimum allowed dual edge length
         inline scalar minimumEdgeLength(const point& pt) const;
 
@@ -154,9 +157,8 @@ class conformalVoronoiMesh
         //  specified spacing
         void insertSurfacePointPairs
         (
-            const List<scalar>& surfacePpDist,
-            const List<point>& surfacePoints,
-            const List<vector>& surfaceNormals,
+            const List<pointIndexHit>& surfaceHits,
+            const List<label>& hitSurfaces,
             const fileName fName = fileName::null
         );
 
@@ -166,7 +168,8 @@ class conformalVoronoiMesh
         void insertEdgePointGroups
         (
             const List<pointIndexHit>& edgeHits,
-            const labelList& featuresHit
+            const labelList& featuresHit,
+            const fileName fName = fileName::null
         );
 
         //- Insert points to conform to an external edge
@@ -225,11 +228,20 @@ class conformalVoronoiMesh
 
         //- Check to see if dual cell specified by given vertex iterator
         //  intersects the boundary and hence reqires a point-pair.
-        bool dualCellSurfaceIntersection
+        bool dualCellSurfaceAnyIntersection
         (
             const Triangulation::Finite_vertices_iterator& vit
         ) const;
 
+        //- Find the "worst" protrusion of a dual cell through the surface,
+        //- subject to the tolerance
+        void dualCellLargestSurfaceProtrusion
+        (
+            const Triangulation::Finite_vertices_iterator& vit,
+            pointIndexHit& surfHit,
+            label& hitSurface
+        ) const;
+
         //- Dual calculation
         void calcDualMesh
         (
diff --git a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index d3e9e8c9ae..bee8e25048 100644
--- a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -62,6 +62,15 @@ inline Foam::scalar Foam::conformalVoronoiMesh::surfaceSearchDistanceSqr
 }
 
 
+inline Foam::scalar Foam::conformalVoronoiMesh::maxSurfaceProtrusion
+(
+    const point& pt
+) const
+{
+    return targetCellSize(pt)*cvMeshControls_.maxSurfaceProtrusionCoeff();
+}
+
+
 inline Foam::scalar Foam::conformalVoronoiMesh::minimumEdgeLength
 (
     const point& pt
diff --git a/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 7b4c52dad7..87e115a543 100644
--- a/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -329,17 +329,51 @@ bool Foam::conformationSurfaces::findSurfaceAnyIntersection
 }
 
 
-void Foam::conformationSurfaces::findSurfaceNearestAndNormal
+void Foam::conformationSurfaces::findSurfaceAnyIntersection
 (
-    const point& sample,
-    scalar nearestDistSqr,
-    pointIndexHit& pHit,
-    vector& normal
+    point start,
+    point end,
+    pointIndexHit& surfHit,
+    label& hitSurface
 ) const
 {
     labelList hitSurfaces;
     List<pointIndexHit> hitInfo;
 
+    searchableSurfacesQueries::findAnyIntersection
+    (
+        allGeometry_,
+        surfaces_,
+        pointField(1, start),
+        pointField(1, end),
+        hitSurfaces,
+        hitInfo
+    );
+
+    surfHit = hitInfo[0];
+
+    if (surfHit.hit())
+    {
+        // hitSurfaces has returned the index of the entry in surfaces_ that was
+        // found, not the index of the surface in allGeometry_, translating this
+        // to allGeometry_
+
+        hitSurface = surfaces_[hitSurfaces[0]];
+    }
+}
+
+
+void Foam::conformationSurfaces::findSurfaceNearest
+(
+    const point& sample,
+    scalar nearestDistSqr,
+    pointIndexHit& surfHit,
+    label& hitSurface
+) const
+{
+    labelList hitSurfaces;
+    List<pointIndexHit> surfaceHits;
+
     searchableSurfacesQueries::findNearest
     (
         allGeometry_,
@@ -347,36 +381,30 @@ void Foam::conformationSurfaces::findSurfaceNearestAndNormal
         pointField(1, sample),
         scalarField(1, nearestDistSqr),
         hitSurfaces,
-        hitInfo
+        surfaceHits
     );
 
-    pHit = hitInfo[0];
+    surfHit = surfaceHits[0];
 
-    if (pHit.hit())
+    if (surfHit.hit())
     {
-        vectorField normals;
-
         // hitSurfaces has returned the index of the entry in surfaces_ that was
         // found, not the index of the surface in allGeometry_, translating this
-        // on access to allGeometry_ for the normal lookup.
+        // to allGeometry_
 
-        allGeometry_[surfaces_[hitSurfaces[0]]].getNormal(hitInfo, normals);
-
-        normal = normals[0];
+        hitSurface = surfaces_[hitSurfaces[0]];
     }
 }
 
 
-void Foam::conformationSurfaces::findSurfaceNearestAndNormal
+void Foam::conformationSurfaces::findSurfaceNearest
 (
     const pointField& samples,
     const scalarField& nearestDistSqr,
-    List<pointIndexHit>& hitInfo,
-    vectorField& normals
+    List<pointIndexHit>& surfaceHits,
+    labelList& hitSurfaces
 ) const
 {
-    labelList hitSurfaces;
-
     searchableSurfacesQueries::findNearest
     (
         allGeometry_,
@@ -384,31 +412,18 @@ void Foam::conformationSurfaces::findSurfaceNearestAndNormal
         samples,
         nearestDistSqr,
         hitSurfaces,
-        hitInfo
+        surfaceHits
     );
 
-    forAll(hitInfo, i)
+    forAll(surfaceHits, i)
     {
-        const pointIndexHit& pHit(hitInfo[i]);
-
-        if (pHit.hit())
+        if (surfaceHits[i].hit())
         {
             // hitSurfaces has returned the index of the entry in surfaces_ that
             // was found, not the index of the surface in allGeometry_,
-            // translating this to the surface in allGeometry_ for the normal
-            // lookup.
+            // translating this to the surface in allGeometry_.
 
             hitSurfaces[i] = surfaces_[hitSurfaces[i]];
-
-            vectorField norm(1);
-
-            allGeometry_[hitSurfaces[i]].getNormal
-            (
-                List<pointIndexHit>(1, pHit),
-                norm
-            );
-
-            normals[i] = norm[0];
         }
     }
 }
diff --git a/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H b/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
index 1ab1d91ef0..3f8d847ea8 100644
--- a/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
+++ b/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
@@ -154,22 +154,32 @@ public:
             // Finding if the line joining start and end intersects the surface
             bool findSurfaceAnyIntersection(point start, point end) const;
 
+            //- Finding if the line joining start and end intersects the surface
+            //  and returning the hit and surface information
+            void findSurfaceAnyIntersection
+            (
+                point start,
+                point end,
+                pointIndexHit& surfHit,
+                label& hitSurface
+            ) const;
+
             //- Find the nearest point to the sample and return it to the
             //  pointIndexHit, and the normal at the hit location, if found.
-            void findSurfaceNearestAndNormal
+            void findSurfaceNearest
             (
                 const point& sample,
                 scalar nearestDistSqr,
-                pointIndexHit& pHit,
-                vector& normal
+                pointIndexHit& surfHit,
+                label& hitSurface
             ) const;
 
-            void findSurfaceNearestAndNormal
+            void findSurfaceNearest
             (
                 const pointField& samples,
                 const scalarField& nearestDistSqr,
                 List<pointIndexHit>& hitInfo,
-                vectorField& normals
+                labelList& hitSurfaces
             ) const;
 
 
diff --git a/src/conformalVoronoiMesh/cvControls/cvControls.C b/src/conformalVoronoiMesh/cvControls/cvControls.C
index 7af66078f0..597913c7fc 100644
--- a/src/conformalVoronoiMesh/cvControls/cvControls.C
+++ b/src/conformalVoronoiMesh/cvControls/cvControls.C
@@ -51,18 +51,17 @@ Foam::cvControls::cvControls
         surfDict.lookup("surfaceSearchDistanceCoeff")
     );
 
-    maxQuadAngle_= readScalar
+    maxSurfaceProtrusionCoeff_ = readScalar
     (
-        surfDict.lookup("maxQuadAngle")
+        surfDict.lookup("maxSurfaceProtrusionCoeff")
     );
 
+    maxQuadAngle_= readScalar(surfDict.lookup("maxQuadAngle"));
+
     // Motion control controls
     const dictionary& motionDict(cvMeshDict_.subDict("motionControl"));
 
-    defaultCellSize_ = readScalar
-    (
-        motionDict.lookup("defaultCellSize")
-    );
+    defaultCellSize_ = readScalar(motionDict.lookup("defaultCellSize"));
 
     const dictionary& insertionDict
     (
diff --git a/src/conformalVoronoiMesh/cvControls/cvControls.H b/src/conformalVoronoiMesh/cvControls/cvControls.H
index 09aa1c6f0f..5d34e5d99e 100644
--- a/src/conformalVoronoiMesh/cvControls/cvControls.H
+++ b/src/conformalVoronoiMesh/cvControls/cvControls.H
@@ -70,6 +70,11 @@ class cvControls
             //  target cell size
             scalar surfaceSearchDistanceCoeff_;
 
+            //- Maximum allowable protrusion through the surface before
+            //  conformation points are added - fraction of the local target
+            //  cell size
+            scalar maxSurfaceProtrusionCoeff_;
+
             //- Maximum quadrant angle allowed at a concave edge before
             //  additional "mitering" lines are added
             scalar maxQuadAngle_;
@@ -130,6 +135,9 @@ public:
             //- Return the surfaceSearchDistanceCoeff
             inline scalar surfaceSearchDistanceCoeff() const;
 
+            //- Return the maxSurfaceProtrusionCoeff
+            inline scalar maxSurfaceProtrusionCoeff() const;
+
             //- Return the maxQuadAngle
             inline scalar maxQuadAngle() const;
 
diff --git a/src/conformalVoronoiMesh/cvControls/cvControlsI.H b/src/conformalVoronoiMesh/cvControls/cvControlsI.H
index f5ff56a459..6f735d1815 100644
--- a/src/conformalVoronoiMesh/cvControls/cvControlsI.H
+++ b/src/conformalVoronoiMesh/cvControls/cvControlsI.H
@@ -52,6 +52,12 @@ inline Foam::scalar Foam::cvControls::surfaceSearchDistanceCoeff() const
 }
 
 
+inline Foam::scalar Foam::cvControls::maxSurfaceProtrusionCoeff() const
+{
+    return maxSurfaceProtrusionCoeff_;
+}
+
+
 inline Foam::scalar Foam::cvControls::maxQuadAngle() const
 {
     return maxQuadAngle_;
diff --git a/src/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C b/src/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
index 0d72fda33a..55c7f58fdb 100644
--- a/src/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
+++ b/src/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
@@ -66,6 +66,8 @@ std::vector<Vb::Point> pointFile::initialPoints() const
         )
     );
 
+    Info<< "    Inserting points from file " << pointFileName_ << endl;
+
     std::vector<Vb::Point> initialPoints;
 
     Field<bool> insidePoints = cvMesh_.geometryToConformTo().wellInside