From 138f149b47eca977215ad7387134d8c7f080398c Mon Sep 17 00:00:00 2001
From: laurence <laurence>
Date: Tue, 20 Mar 2012 17:43:27 +0000
Subject: [PATCH] ENH: surfaceFeatures: Add constructor that takes an edgeMesh

---
 .../algorithms/indexedOctree/indexedOctree.C  |   6 +-
 .../surfaceFeatures/surfaceFeatures.C         | 261 ++++++++++++++----
 .../surfaceFeatures/surfaceFeatures.H         |  35 ++-
 3 files changed, 243 insertions(+), 59 deletions(-)

diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
index 024f7d27ec..f16dd2595c 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
+++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
@@ -2832,7 +2832,7 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
 ) const
 {
     label nearestShapeI = -1;
-    point nearestPoint;
+    point nearestPoint = vector::zero;
 
     if (nodes_.size())
     {
@@ -2847,10 +2847,6 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
             nearestPoint
         );
     }
-    else
-    {
-        nearestPoint = vector::zero;
-    }
 
     return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
 }
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
index e39eae5704..465734ff87 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
@@ -33,11 +33,15 @@ License
 #include "OFstream.H"
 #include "IFstream.H"
 #include "unitConversion.H"
+#include "EdgeMap.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(Foam::surfaceFeatures, 0);
 
+const Foam::scalar Foam::surfaceFeatures::parallelTolerance =
+    sin(degToRad(1.0));
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -170,7 +174,10 @@ void Foam::surfaceFeatures::setFromStatus(const List<edgeStatus>& edgeStat)
 
 
 //construct feature points where more than 2 feature edges meet
-void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
+void Foam::surfaceFeatures::calcFeatPoints
+(
+    const List<edgeStatus>& edgeStat
+)
 {
     DynamicList<label> featurePoints(surf_.nPoints()/1000);
 
@@ -200,6 +207,60 @@ void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
 }
 
 
+void Foam::surfaceFeatures::classifyFeatureAngles
+(
+    const labelListList& edgeFaces,
+    List<edgeStatus>& edgeStat,
+    const scalar minCos
+) const
+{
+    const vectorField& faceNormals = surf_.faceNormals();
+    const pointField& points = surf_.points();
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        if (eFaces.size() != 2)
+        {
+            // Non-manifold. What to do here? Is region edge? external edge?
+            edgeStat[edgeI] = REGION;
+        }
+        else
+        {
+            label face0 = eFaces[0];
+            label face1 = eFaces[1];
+
+            if (surf_[face0].region() != surf_[face1].region())
+            {
+                edgeStat[edgeI] = REGION;
+            }
+            else
+            {
+                if ((faceNormals[face0] & faceNormals[face1]) < minCos)
+                {
+
+                    // Check if convex or concave by looking at angle
+                    // between face centres and normal
+                    vector f0Tof1 =
+                        surf_[face1].centre(points)
+                      - surf_[face0].centre(points);
+
+                    if ((f0Tof1 & faceNormals[face0]) > 0.0)
+                    {
+                        edgeStat[edgeI] = INTERNAL;
+                    }
+                    else
+                    {
+                        edgeStat[edgeI] = EXTERNAL;
+                    }
+                }
+            }
+        }
+    }
+}
+
+
 // Returns next feature edge connected to pointI with correct value.
 Foam::label Foam::surfaceFeatures::nextFeatEdge
 (
@@ -436,6 +497,80 @@ Foam::surfaceFeatures::surfaceFeatures
 }
 
 
+Foam::surfaceFeatures::surfaceFeatures
+(
+    const triSurface& surf,
+    const pointField& points,
+    const edgeList& edges,
+    const scalar mergeTol
+)
+:
+    surf_(surf),
+    featurePoints_(0),
+    featureEdges_(0),
+    externalStart_(0),
+    internalStart_(0)
+{
+    // Match edge mesh edges with the triSurface edges
+
+    const labelListList& surfEdgeFaces = surf_.edgeFaces();
+    const edgeList& surfEdges = surf_.edges();
+
+    scalar mergeTolSqr = sqr(mergeTol);
+
+    EdgeMap<label> dynFeatEdges(2*edges.size());
+    DynamicList<labelList> dynFeatureEdgeFaces(edges.size());
+
+    labelList edgeLabel;
+
+    nearestFeatEdge
+    (
+        edges,
+        points,
+        mergeTolSqr,
+        edgeLabel     // label of surface edge or -1
+    );
+
+    label count = 0;
+    forAll(edgeLabel, sEdgeI)
+    {
+        const label sEdge = edgeLabel[sEdgeI];
+
+        if (sEdge == -1)
+        {
+            continue;
+        }
+
+        dynFeatEdges.insert(surfEdges[sEdge], count++);
+        dynFeatureEdgeFaces.append(surfEdgeFaces[sEdge]);
+    }
+
+    // Find whether an edge is external or internal
+    List<edgeStatus> edgeStat(dynFeatEdges.size(), NONE);
+
+    classifyFeatureAngles(dynFeatureEdgeFaces, edgeStat, GREAT);
+
+    // Transfer the edge status to a list encompassing all edges in the surface
+    // so that calcFeatPoints can be used.
+    List<edgeStatus> allEdgeStat(surf_.nEdges(), NONE);
+
+    forAll(allEdgeStat, eI)
+    {
+        EdgeMap<label>::const_iterator iter = dynFeatEdges.find(surfEdges[eI]);
+
+        if (iter != dynFeatEdges.end())
+        {
+            allEdgeStat[eI] = edgeStat[iter()];
+        }
+    }
+
+    edgeStat.clear();
+    dynFeatEdges.clear();
+
+    setFromStatus(allEdgeStat);
+}
+
+
 //- Construct as copy
 Foam::surfaceFeatures::surfaceFeatures(const surfaceFeatures& sf)
 :
@@ -496,54 +631,10 @@ void Foam::surfaceFeatures::findFeatures(const scalar includedAngle)
 {
     scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
 
-    const labelListList& edgeFaces = surf_.edgeFaces();
-    const vectorField& faceNormals = surf_.faceNormals();
-    const pointField& points = surf_.points();
-
     // Per edge whether is feature edge.
     List<edgeStatus> edgeStat(surf_.nEdges(), NONE);
 
-    forAll(edgeFaces, edgeI)
-    {
-        const labelList& eFaces = edgeFaces[edgeI];
-
-        if (eFaces.size() != 2)
-        {
-            // Non-manifold. What to do here? Is region edge? external edge?
-            edgeStat[edgeI] = REGION;
-        }
-        else
-        {
-            label face0 = eFaces[0];
-            label face1 = eFaces[1];
-
-            if (surf_[face0].region() != surf_[face1].region())
-            {
-                edgeStat[edgeI] = REGION;
-            }
-            else
-            {
-                if ((faceNormals[face0] & faceNormals[face1]) < minCos)
-                {
-
-                    // Check if convex or concave by looking at angle
-                    // between face centres and normal
-                    vector f0Tof1 =
-                        surf_[face1].centre(points)
-                        - surf_[face0].centre(points);
-
-                    if ((f0Tof1 & faceNormals[face0]) > 0.0)
-                    {
-                        edgeStat[edgeI] = INTERNAL;
-                    }
-                    else
-                    {
-                        edgeStat[edgeI] = EXTERNAL;
-                    }
-                }
-            }
-        }
-    }
+    classifyFeatureAngles(surf_.edgeFaces(), edgeStat, minCos);
 
     setFromStatus(edgeStat);
 }
@@ -1147,6 +1238,9 @@ void Foam::surfaceFeatures::nearestSurfEdge
 
     const pointField& localPoints = surf_.localPoints();
 
+    treeBoundBox searchDomain(localPoints);
+    searchDomain.inflate(0.1);
+
     indexedOctree<treeDataEdge> ppTree
     (
         treeDataEdge
@@ -1156,7 +1250,7 @@ void Foam::surfaceFeatures::nearestSurfEdge
             localPoints,
             selectedEdges
         ),          // all information needed to do geometric checks
-        treeBoundBox(localPoints),  // overall search domain
+        searchDomain,  // overall search domain
         8,      // maxLevel
         10,     // leafsize
         3.0     // duplicity
@@ -1216,6 +1310,8 @@ void Foam::surfaceFeatures::nearestSurfEdge
     pointOnEdge.setSize(selectedSampleEdges.size());
     pointOnFeature.setSize(selectedSampleEdges.size());
 
+    treeBoundBox searchDomain(surf_.localPoints());
+
     indexedOctree<treeDataEdge> ppTree
     (
         treeDataEdge
@@ -1224,11 +1320,11 @@ void Foam::surfaceFeatures::nearestSurfEdge
             surf_.edges(),
             surf_.localPoints(),
             selectedEdges
-        ),          // all information needed to do geometric checks
-        treeBoundBox(surf_.localPoints()),  // overall search domain
-        8,      // maxLevel
-        10,     // leafsize
-        3.0     // duplicity
+        ),              // all information needed to do geometric checks
+        searchDomain,   // overall search domain
+        8,              // maxLevel
+        10,             // leafsize
+        3.0             // duplicity
     );
 
     forAll(selectedSampleEdges, i)
@@ -1262,6 +1358,67 @@ void Foam::surfaceFeatures::nearestSurfEdge
 }
 
 
+void Foam::surfaceFeatures::nearestFeatEdge
+(
+    const edgeList& edges,
+    const pointField& points,
+    scalar searchSpanSqr,   // Search span
+    labelList& edgeLabel
+) const
+{
+    edgeLabel = labelList(surf_.nEdges(), -1);
+
+    treeBoundBox searchDomain(points);
+    searchDomain.inflate(0.1);
+
+    indexedOctree<treeDataEdge> ppTree
+    (
+        treeDataEdge
+        (
+            false,
+            edges,
+            points,
+            identity(edges.size())
+        ),          // all information needed to do geometric checks
+        searchDomain,  // overall search domain
+        8,      // maxLevel
+        10,     // leafsize
+        3.0     // duplicity
+    );
+
+    const edgeList& surfEdges = surf_.edges();
+    const pointField& surfLocalPoints = surf_.localPoints();
+
+    forAll(surfEdges, edgeI)
+    {
+        const edge& sample = surfEdges[edgeI];
+
+        const point& startPoint = surfLocalPoints[sample.start()];
+        const point& midPoint = sample.centre(surfLocalPoints);
+
+        pointIndexHit infoMid = ppTree.findNearest
+        (
+            midPoint,
+            searchSpanSqr
+        );
+
+        if (infoMid.hit())
+        {
+            const vector surfEdgeDir = midPoint - startPoint;
+
+            const edge& featEdge = edges[infoMid.index()];
+            const vector featEdgeDir = featEdge.vec(points);
+
+            // Check that the edges are nearly parallel
+            if (mag(surfEdgeDir ^ featEdgeDir) < parallelTolerance)
+            {
+                edgeLabel[edgeI] = edgeI;
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 void Foam::surfaceFeatures::operator=(const surfaceFeatures& rhs)
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
index 7f7101b662..0fca9798a2 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
@@ -95,6 +95,11 @@ private:
         {}
     };
 
+    // Static data
+
+        //- Tolerance for determining whether two vectors are parallel
+        static const scalar parallelTolerance;
+
 
     // Private data
 
@@ -130,6 +135,14 @@ private:
         //- Construct feature points where more than 2 feature edges meet
         void calcFeatPoints(const List<edgeStatus>&);
 
+        //- Classify the angles of the feature edges
+        void classifyFeatureAngles
+        (
+            const labelListList& edgeFaces,
+            List<edgeStatus>& edgeStat,
+            const scalar minCos
+        ) const;
+
         //- Choose next unset feature edge.
         label nextFeatEdge
         (
@@ -186,6 +199,15 @@ public:
         //- Construct from file
         surfaceFeatures(const triSurface&, const fileName& fName);
 
+        //- Construct from pointField and edgeList (edgeMesh)
+        surfaceFeatures
+        (
+            const triSurface&,
+            const pointField& points,
+            const edgeList& edges,
+            const scalar mergeTol = 1e-6
+        );
+
         //- Construct as copy
         surfaceFeatures(const surfaceFeatures&);
 
@@ -327,9 +349,8 @@ public:
                 pointField& edgePoint
             ) const;
 
-
             //- Find nearest surface edge (out of selectedEdges) for each
-            // sample edge.
+            //  sample edge.
             //  Sets:
             //  - edgeLabel         : label of surface edge.
             //  - pointOnEdge       : exact position of nearest point on edge.
@@ -347,6 +368,16 @@ public:
                 pointField& pointOnFeature  // point on sample edge
             ) const;
 
+            //- Find nearest feature edge to each surface edge. Uses the
+            //  mid-point of the surface edges.
+            void nearestFeatEdge
+            (
+                const edgeList& edges,
+                const pointField& points,
+                scalar searchSpanSqr,
+                labelList& edgeLabel
+            ) const;
+
 
         // Write