From 83d81b7441094e0bed754af3e6c37ab8b5d7155f Mon Sep 17 00:00:00 2001
From: laurence <laurence>
Date: Tue, 24 Jan 2012 12:48:38 +0000
Subject: [PATCH] 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).
---
 .../conformalVoronoiMesh.H                    |  15 +-
 .../conformalVoronoiMeshConformToSurface.C    | 190 ++++++++++++++----
 .../conformalVoronoiMeshIO.C                  |   2 +-
 .../conformationSurfaces.C                    |   1 +
 4 files changed, 165 insertions(+), 43 deletions(-)

diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 0ac85c61dd..2e38a504e1 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/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
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index 4a76612dd2..7d8ec2619c 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -127,26 +127,32 @@ 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.0,       // maximum ratio of cubes v.s. cells
+        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);
 
@@ -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 (existingPoint == info.hitPoint())
+//            if (pI == info.index())
 //            {
-//                toRemove[pI] = true;
-//                //Info<< "Replace " << pt << " and " << info.hitPoint() << " with " << newPt << endl;
+//                //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);
+////
+////        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();
 //    }
@@ -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;
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 1b2f46eed0..bc9d04dfb1 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/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()));
         }
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 3d256b2305..b91d40121d 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -729,6 +729,7 @@ void Foam::conformationSurfaces::findEdgeNearestByType
     }
 }
 
+
 void Foam::conformationSurfaces::findAllNearestEdges
 (
     const point& sample,