From 24d09bdc8eea0d106a857fbd73d36aab8f2d864d Mon Sep 17 00:00:00 2001
From: graham <graham.macpherson@strath.ac.uk>
Date: Wed, 10 Sep 2008 19:43:48 +0100
Subject: [PATCH] removed plane.planePlane intersections. processing edges to
 remove groups of close control points.

---
 .../mesh/generation/CV3DMesher/CV3D.H         |   8 +
 .../mesh/generation/CV3DMesher/calcDualMesh.C |  22 +-
 .../CV3DMesher/insertFeaturePoints.C          |   2 -
 .../insertSurfaceNearestPointPairs.C          | 224 +++++++++++++++++-
 4 files changed, 238 insertions(+), 18 deletions(-)

diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3D.H b/applications/utilities/mesh/generation/CV3DMesher/CV3D.H
index c34c6a796f..be62d61f5d 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/CV3D.H
+++ b/applications/utilities/mesh/generation/CV3DMesher/CV3D.H
@@ -236,6 +236,14 @@ private:
             const Triangulation::Finite_vertices_iterator& vit
         ) const;
 
+        //- Smooths the points positioned along a feature edge by collecting
+        //- bunches together and deleting excess points.
+        void smoothEdgePositions
+        (
+            DynamicList<point>& edgePoints,
+            DynamicList<label>& edgeLabels
+        ) const;
+
         //- Insert point-pairs at the nearest points on the surface to the
         //  control vertex of dual-cells which intersect the boundary in order
         //  to provide a boundary-layer mesh.
diff --git a/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C b/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C
index 018abb495d..8e9e48bc66 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C
@@ -271,7 +271,7 @@ void Foam::CV3D::calcDualMesh
     {
         labelList copyOwner(owner.size());
 
-        forAll (sortingIndices, sI)
+        forAll(sortingIndices, sI)
         {
             copyOwner[sI] = owner[sortingIndices[sI]];
         }
@@ -282,7 +282,7 @@ void Foam::CV3D::calcDualMesh
     {
         labelList copyNeighbour(neighbour.size());
 
-        forAll (sortingIndices, sI)
+        forAll(sortingIndices, sI)
         {
             copyNeighbour[sI] = neighbour[sortingIndices[sI]];
         }
@@ -293,7 +293,7 @@ void Foam::CV3D::calcDualMesh
     {
         faceList copyFaces(faces.size());
 
-        forAll (sortingIndices, sI)
+        forAll(sortingIndices, sI)
         {
             copyFaces[sI] = faces[sortingIndices[sI]];
         }
@@ -320,7 +320,7 @@ void Foam::CV3D::calcDualMesh
 
     ownerCellJumps.shrink();
 
-    forAll (ownerCellJumps, oCJ)
+    forAll(ownerCellJumps, oCJ)
     {
         label start = ownerCellJumps[oCJ];
 
@@ -339,7 +339,7 @@ void Foam::CV3D::calcDualMesh
 
         SortableList<label> sortedNeighbourBlock(neighbourBlock);
 
-        forAll (sortedNeighbourBlock, sNB)
+        forAll(sortedNeighbourBlock, sNB)
         {
             sortingIndices[start + sNB] =
                 sortedNeighbourBlock.indices()[sNB] + start;
@@ -351,7 +351,7 @@ void Foam::CV3D::calcDualMesh
     {
         labelList copyOwner(owner.size());
 
-        forAll (sortingIndices, sI)
+        forAll(sortingIndices, sI)
         {
             copyOwner[sI] = owner[sortingIndices[sI]];
         }
@@ -362,7 +362,7 @@ void Foam::CV3D::calcDualMesh
     {
         labelList copyNeighbour(neighbour.size());
 
-        forAll (sortingIndices, sI)
+        forAll(sortingIndices, sI)
         {
             copyNeighbour[sI] = neighbour[sortingIndices[sI]];
         }
@@ -373,7 +373,7 @@ void Foam::CV3D::calcDualMesh
     {
         faceList copyFaces(faces.size());
 
-        forAll (sortingIndices, sI)
+        forAll(sortingIndices, sI)
         {
             copyFaces[sI] = faces[sortingIndices[sI]];
         }
@@ -385,7 +385,7 @@ void Foam::CV3D::calcDualMesh
 
     label nBoundaryFaces = 0;
 
-    forAll (patchFaces, p)
+    forAll(patchFaces, p)
     {
         patchFaces[p].shrink();
 
@@ -402,9 +402,9 @@ void Foam::CV3D::calcDualMesh
 
     owner.setSize(nInternalFaces + nBoundaryFaces);
 
-    forAll (patchFaces, p)
+    forAll(patchFaces, p)
     {
-        forAll (patchFaces[p], f)
+        forAll(patchFaces[p], f)
         {
             faces[dualFacei] = patchFaces[p][f];
 
diff --git a/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C b/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
index f0a092c938..a9f111e049 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
@@ -377,8 +377,6 @@ void Foam::CV3D::insertFeaturePoints()
                         )
                     )/cos(phi) - 1;
 
-                    Info<< tab << guard << endl;
-
                     point internalPtF = concaveEdgeExternalPt + (2 + guard) *
                         (concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
 
diff --git a/applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C b/applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C
index e0b690f822..d53de95f43 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C
@@ -83,6 +83,208 @@ bool Foam::CV3D::dualCellSurfaceIntersection
 }
 
 
+void Foam::CV3D::smoothEdgePositions
+(
+    DynamicList<point>& edgePoints,
+    DynamicList<label>& edgeLabels
+) const
+{
+    const pointField& localPts = qSurf_.localPoints();
+    const edgeList& edges = qSurf_.edges();
+
+    // Sort by edge label then by distance from the start of the edge
+
+    SortableList<label> sortedEdgeLabels(edgeLabels);
+
+    labelList sortingIndices = sortedEdgeLabels.indices();
+
+    {
+        labelList copyEdgeLabels(edgeLabels.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyEdgeLabels[sI] = edgeLabels[sortingIndices[sI]];
+        }
+
+        edgeLabels.transfer(copyEdgeLabels);
+    }
+
+    {
+        List<point> copyEdgePoints(edgePoints.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyEdgePoints[sI] = edgePoints[sortingIndices[sI]];
+        }
+
+        edgePoints.transfer(copyEdgePoints);
+    }
+
+    List<scalar> edgeDistances(edgePoints.size());
+
+    forAll(edgeDistances, eD)
+    {
+        const point& edgeStart = localPts[edges[edgeLabels[eD]].start()];
+
+        edgeDistances[eD] = mag(edgeStart - edgePoints[eD]);
+    }
+
+    // Sort by edgeDistances in blocks of edgeLabel
+
+    DynamicList<label> edgeLabelJumps;
+
+    // Force first edgeLabel to be a jump
+    edgeLabelJumps.append(0);
+
+    for (label eL = 1; eL < edgeLabels.size(); eL++)
+    {
+        if (edgeLabels[eL] > edgeLabels[eL-1])
+        {
+            edgeLabelJumps.append(eL);
+        }
+    }
+
+    edgeLabelJumps.shrink();
+
+    forAll(edgeLabelJumps, eLJ)
+    {
+        label start = edgeLabelJumps[eLJ];
+
+        label length;
+
+        if (eLJ == edgeLabelJumps.size() - 1)
+        {
+            length = edgeLabels.size() - start;
+        }
+        else
+        {
+            length = edgeLabelJumps[eLJ + 1] - start;
+        }
+
+        SubList<scalar> edgeDistanceBlock(edgeDistances, length, start);
+
+        SortableList<scalar> sortedEdgeDistanceBlock(edgeDistanceBlock);
+
+        forAll(sortedEdgeDistanceBlock, sEDB)
+        {
+            sortingIndices[start + sEDB] =
+                sortedEdgeDistanceBlock.indices()[sEDB] + start;
+        }
+    }
+
+    {
+        List<point> copyEdgePoints(edgePoints.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyEdgePoints[sI] = edgePoints[sortingIndices[sI]];
+        }
+
+        edgePoints.transfer(copyEdgePoints);
+    }
+
+    {
+        List<scalar> copyEdgeDistances(edgeDistances.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyEdgeDistances[sI] = edgeDistances[sortingIndices[sI]];
+        }
+
+        edgeDistances.transfer(copyEdgeDistances);
+    }
+
+    // Process the points along each edge (in blocks of edgeLabel) performing 3
+    // functions:
+    // 1: move points away from feature edges
+    // 2: aggregate tight groups of points into one point
+    // 3: adjust the spacing of remaining points on a pair by pair basis to
+    //    remove excess points
+
+    DynamicList<point> tempEdgePoints(edgePoints.size());
+
+    DynamicList<label> tempEdgeLabels(edgeLabels.size());
+
+    forAll(edgeLabelJumps, eLJ)
+    {
+        label start = edgeLabelJumps[eLJ];
+
+        label edgeI = edgeLabels[start];
+
+        label length;
+
+        if (eLJ == edgeLabelJumps.size() - 1)
+        {
+            length = edgeLabels.size() - start;
+        }
+        else
+        {
+            length = edgeLabelJumps[eLJ + 1] - start;
+        }
+
+        List<point> edgePointBlock(SubList<point>(edgePoints, length, start));
+
+        List<scalar> edgeDistanceBlock(SubList<scalar>(edgeDistances, length, start));
+
+        scalar maxPointGroupSpacing = tols_.ppDist;
+
+        label groupSize = 0;
+
+        point newEdgePoint(vector::zero);
+
+        if (edgeDistanceBlock.size() == 1)
+        {
+            tempEdgePoints.append(edgePointBlock[0]);
+
+            tempEdgeLabels.append(edgeI);
+        }
+        else if ((edgeDistanceBlock[1] - edgeDistanceBlock[0]) > maxPointGroupSpacing)
+        {
+            tempEdgePoints.append(edgePointBlock[0]);
+
+            tempEdgeLabels.append(edgeI);
+        }
+
+        for (label eDB = 1; eDB < edgeDistanceBlock.size(); eDB++)
+        {
+            const scalar& edgeDist = edgeDistanceBlock[eDB];
+            const scalar& previousEdgeDist = edgeDistanceBlock[eDB - 1];
+
+            if ((edgeDist - previousEdgeDist) < maxPointGroupSpacing)
+            {
+                newEdgePoint += edgePointBlock[eDB];
+
+                groupSize++;
+            }
+            else if (groupSize > 0)
+            {
+                // A point group has been formed and has finished
+
+                newEdgePoint /= groupSize;
+
+                tempEdgePoints.append(newEdgePoint);
+
+                tempEdgeLabels.append(edgeI);
+
+                newEdgePoint = vector::zero;
+
+                groupSize = 0;
+            }
+            else
+            {
+                tempEdgePoints.append(edgePointBlock[eDB]);
+
+                tempEdgeLabels.append(edgeI);
+            }
+        }
+    }
+
+    edgePoints.transfer(tempEdgePoints.shrink());
+
+    edgeLabels.transfer(tempEdgeLabels.shrink());
+}
+
+
 void Foam::CV3D::insertPointPairs
 (
     const DynamicList<point>& nearSurfacePoints,
@@ -164,12 +366,20 @@ void Foam::CV3D::insertEdgePointGroups
         plane planeA(edgePt - tols_.ppDist*nA, nA);
         plane planeB(edgePt - tols_.ppDist*nB, nB);
 
-        plane planeF(edgePt, (nA^nB));
+        // Finding the nearest point on the intersecting line to the edge point.
+        // Floating point errors often encountered using planePlaneIntersect
+        // plane planeF(edgePt, (nA^nB));
+        // point refPt = planeF.planePlaneIntersect(planeA,planeB);
 
-        Info<< "CHANGE TO FINDING THE NEAREST POINT ON THE INTERSECTING LINE TO THE EDGE POINT. "
-            << "Floating point errors in planePlane." << endl;
+        plane::ray planeIntersect(planeA.planeIntersect(planeB));
 
-        point refPt = planeF.planePlaneIntersect(planeA,planeB);
+        pointHit refPtHit = linePointRef
+        (
+            planeIntersect.refPoint() + 2*tols_.span*planeIntersect.dir(),
+            planeIntersect.refPoint() - 2*tols_.span*planeIntersect.dir()
+        ).nearestDist(edgePt);
+
+        point refPt = refPtHit.hitPoint();
 
         point faceAVert =
         localPts[triSurfaceTools::oppositeVertex(qSurf_, faceA, edgeI)];
@@ -246,6 +456,8 @@ void Foam::CV3D::insertEdgePointGroups
         }
     }
 
+    Info<< edgePoints.size() << " edge-control locations inserted" << endl;
+
     if (controls_.writeInsertedPointPairs)
     {
         OFstream str(fName);
@@ -323,7 +535,7 @@ void Foam::CV3D::insertSurfaceNearestPointPairs()
     DynamicList<label> edgeLabels(allEdgeLabels.size());
     DynamicList<vector> edgePoints(allEdgePoints.size());
 
-    forAll (allEdgePoints, eP)
+    forAll(allEdgePoints, eP)
     {
         if (allEdgeLabels[eP] >= 0 && allEdgeEndPoints[eP] < 0)
         {
@@ -335,6 +547,8 @@ void Foam::CV3D::insertSurfaceNearestPointPairs()
     edgePoints.shrink();
     edgeLabels.shrink();
 
+    smoothEdgePositions(edgePoints, edgeLabels);
+
     insertEdgePointGroups
     (
         edgePoints,