ENH: cvMesh: Dynamic tree added for surface points.

Now there is a search for the surface points to check whether they are close
to any other surface point. If they are, reject the point unless the angle
between the normals of the surfaces are around 180 degrees, which means that
the points are close but on opposite surfaces (i.e. there is a small gap).

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -658,6 +658,17 @@ private:
             label callCount = 0
         ) const;
 
+        //- Find angle between the normals of two close surface points.
+        scalar angleBetweenSurfacePoints(Foam::point pA, Foam::point pB) const;
+
+        //- Check if a surface point is near another.
+        bool nearSurfacePoint
+        (
+            const pointIndexHit& pHit,
+            DynamicList<Foam::point>& existingSurfacePtLocations,
+            dynamicIndexedOctree<dynamicTreeDataPoint>& surfacePtLocationTree
+        ) const;
+
         //- Check if a location is in the exclusion range of an existing feature
         //- edge conformation location
         bool nearFeatureEdgeLocation
@@ -694,7 +705,9 @@ private:
             DynamicList<label>& featureEdgeFeaturesHit,
             DynamicList<Foam::point>& newEdgeLocations,
             DynamicList<Foam::point>& existingEdgeLocations,
-            dynamicIndexedOctree<dynamicTreeDataPoint>& edgeLocationTree
+            dynamicIndexedOctree<dynamicTreeDataPoint>& edgeLocationTree,
+            DynamicList<Foam::point>& existingSurfacePtLocations,
+            dynamicIndexedOctree<dynamicTreeDataPoint>& surfacePtLocationTree
         ) const;
 
         //- Store the surface conformation with the indices offset to be

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C

@@ -127,23 +127,29 @@ 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)
     );
 
-    overallBb.min() -= Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-    overallBb.max() += Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
     dynamicIndexedOctree<dynamicTreeDataPoint> edgeLocationTree
     (
         dynamicTreeDataPoint(existingEdgeLocations),
         overallBb,  // overall search domain
-        10,         // max levels
+        10,         // max levels; n/a for dynamic tree
-        10.0,       // maximum ratio of cubes v.s. cells
+        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.
 	);
 
@@ -249,7 +255,9 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                         featureEdgeFeaturesHit,
                         newEdgeLocations,
                         existingEdgeLocations,
-                        edgeLocationTree
+                        edgeLocationTree,
+                        existingSurfacePtLocations,
+                        surfacePtLocationTree
                     );
                 }
             }
@@ -466,7 +474,9 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                         featureEdgeFeaturesHit,
                         newEdgeLocations,
                         existingEdgeLocations,
-                        edgeLocationTree
+                        edgeLocationTree,
+                        existingSurfacePtLocations,
+                        surfacePtLocationTree
                     );
                 }
             }
@@ -503,7 +513,9 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                         featureEdgeFeaturesHit,
                         newEdgeLocations,
                         existingEdgeLocations,
-                        edgeLocationTree
+                        edgeLocationTree,
+                        existingSurfacePtLocations,
+                        surfacePtLocationTree
                     );
                 }
             }
@@ -1841,6 +1853,97 @@ void Foam::conformalVoronoiMesh::limitDisplacement
 }
 
 
+Foam::scalar Foam::conformalVoronoiMesh::angleBetweenSurfacePoints
+(
+    Foam::point pA,
+    Foam::point pB
+) const
+{
+    pointIndexHit pAhit;
+    label pAsurfaceHit = -1;
+
+    const scalar searchDist = 1.0*targetCellSize(pA);
+
+    geometryToConformTo_.findSurfaceNearest
+    (
+        pA,
+        searchDist,
+        pAhit,
+        pAsurfaceHit
+    );
+
+    vectorField norm(1);
+
+    allGeometry_[pAsurfaceHit].getNormal
+    (
+        List<pointIndexHit>(1, pAhit),
+        norm
+    );
+
+    const vector nA = norm[0];
+
+    pointIndexHit pBhit;
+    label pBsurfaceHit = -1;
+
+    geometryToConformTo_.findSurfaceNearest
+    (
+        pB,
+        searchDist,
+        pBhit,
+        pBsurfaceHit
+    );
+
+    allGeometry_[pBsurfaceHit].getNormal
+    (
+        List<pointIndexHit>(1, pBhit),
+        norm
+    );
+
+    const vector nB = norm[0];
+
+    return vectorTools::degAngleBetween(nA, nB);
+}
+
+
+bool Foam::conformalVoronoiMesh::nearSurfacePoint
+(
+    const pointIndexHit& pHit,
+    DynamicList<Foam::point>& existingSurfacePtLocations,
+    dynamicIndexedOctree<dynamicTreeDataPoint>& surfacePtLocationTree
+) const
+{
+    const Foam::point pt = pHit.hitPoint();
+
+    scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
+
+    pointIndexHit info = surfacePtLocationTree.findNearest(pt, exclusionRangeSqr);
+
+
+    // 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())
+    {
+        const scalar angle = angleBetweenSurfacePoints(pt, info.hitPoint());
+
+        if (angle > 170.0)
+        {
+            return false;
+        }
+    }
+
+
+    return info.hit();
+}
+
+
 bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
 (
     pointIndexHit& pHit,
@@ -1853,11 +1956,6 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
 
     scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
 
-    // 0.01 and 1000 determined from speed tests, varying the indexedOctree
-    // rebuild frequency and balance of additions to queries.
-
-
-
     label startIndex = existingEdgeLocations.size();
 
     existingEdgeLocations.append(newEdgeLocations);
@@ -1870,12 +1968,9 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
 
     pointIndexHit info = edgeLocationTree.findNearest(pt, exclusionRangeSqr);
 
-
-
     // Searching for the nearest point in existingEdgeLocations using the
     // indexedOctree
 
-
     // Average the points...
 //    if (info.hit())
 //    {
@@ -1883,35 +1978,33 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
 //
 //        pHit.setPoint(newPt);
 //
-//        boolList toRemove(existingEdgeLocations.size(), false);
+//        //boolList toRemove(existingEdgeLocations.size(), false);
 //
 //        forAll(existingEdgeLocations, pI)
 //        {
-//            const Foam::point& existingPoint = existingEdgeLocations[pI];
+//            if (pI == info.index())
-//
-//            if (existingPoint == info.hitPoint())
 //            {
-//                toRemove[pI] = true;
+//                //toRemove[pI] = true;
-//                //Info<< "Replace " << pt << " and " << info.hitPoint() << " with " << newPt << endl;
+//                edgeLocationTree.remove(pI);
 //            }
 //        }
-//
+////
-//        pointField newExistingEdgeLocations(existingEdgeLocations.size());
+////        pointField newExistingEdgeLocations(existingEdgeLocations.size());
-//
+////
-//        label count = 0;
+////        label count = 0;
-//        forAll(existingEdgeLocations, pI)
+////        forAll(existingEdgeLocations, pI)
-//        {
+////        {
-//            if (toRemove[pI] == false)
+////            if (toRemove[pI] == false)
-//            {
+////            {
-//                newExistingEdgeLocations[count++] = existingEdgeLocations[pI];
+////                newExistingEdgeLocations[count++] = existingEdgeLocations[pI];
-//            }
+////            }
-//        }
+////        }
-//
+////
-//        newExistingEdgeLocations.resize(count);
+////        newExistingEdgeLocations.resize(count);
-//
+////
-//        existingEdgeLocations = newExistingEdgeLocations;
+////        existingEdgeLocations = newExistingEdgeLocations;
-//
+////
-//        existingEdgeLocations.append(newPt);
+////        existingEdgeLocations.append(newPt);
 //
 //        return !info.hit();
 //    }
@@ -1994,7 +2087,9 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
     DynamicList<label>& featureEdgeFeaturesHit,
     DynamicList<Foam::point>& newEdgeLocations,
     DynamicList<Foam::point>& existingEdgeLocations,
-    dynamicIndexedOctree<dynamicTreeDataPoint>& edgeLocationTree
+    dynamicIndexedOctree<dynamicTreeDataPoint>& edgeLocationTree,
+    DynamicList<Foam::point>& existingSurfacePtLocations,
+    dynamicIndexedOctree<dynamicTreeDataPoint>& surfacePtLocationTree
 ) const
 {
     bool keepSurfacePoint = true;
@@ -2026,6 +2121,19 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
             }
         }
 
+        if
+        (
+            nearSurfacePoint
+            (
+                surfHitI,
+                existingSurfacePtLocations,
+                surfacePtLocationTree
+            )
+        )
+        {
+            keepSurfacePoint = false;
+        }
+
         List<List<pointIndexHit> > edHitsByFeature;
 
         labelList featuresHit;

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C

@@ -191,7 +191,7 @@ void Foam::conformalVoronoiMesh::writeBoundaryPoints
         ++vit
     )
     {
-        if (vit->pairPoint())
+        if (!vit->internalPoint())
         {
             meshTools::writeOBJ(str, topoint(vit->point()));
         }

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C

@@ -729,6 +729,7 @@ void Foam::conformationSurfaces::findEdgeNearestByType
     }
 }
 
+
 void Foam::conformationSurfaces::findAllNearestEdges
 (
     const point& sample,