ENH: surfaceFeatures: Add feature points with 2 edges given the angle

src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C

@@ -118,7 +118,11 @@ Foam::List<Foam::surfaceFeatures::edgeStatus> Foam::surfaceFeatures::toStatus()
 
 
 // Set from value for every edge
-void Foam::surfaceFeatures::setFromStatus(const List<edgeStatus>& edgeStat)
+void Foam::surfaceFeatures::setFromStatus
+(
+    const List<edgeStatus>& edgeStat,
+    const scalar includedAngle
+)
 {
     // Count
 
@@ -170,20 +174,24 @@ void Foam::surfaceFeatures::setFromStatus(const List<edgeStatus>& edgeStat)
         }
     }
 
-    // Calculate consistent feature points
+    const scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
-    calcFeatPoints(edgeStat);
+
+    calcFeatPoints(edgeStat, minCos);
 }
 
 
 //construct feature points where more than 2 feature edges meet
 void Foam::surfaceFeatures::calcFeatPoints
 (
-    const List<edgeStatus>& edgeStat
+    const List<edgeStatus>& edgeStat,
+    const scalar minCos
 )
 {
     DynamicList<label> featurePoints(surf_.nPoints()/1000);
 
     const labelListList& pointEdges = surf_.pointEdges();
+    const edgeList& edges = surf_.edges();
+    const pointField& localPoints = surf_.localPoints();
 
     forAll(pointEdges, pointI)
     {
@@ -203,6 +211,27 @@ void Foam::surfaceFeatures::calcFeatPoints
         {
             featurePoints.append(pointI);
         }
+        else if (nFeatEdges == 2)
+        {
+            // Check the angle between the two edges
+            DynamicList<vector> edgeVecs(2);
+
+            forAll(pEdges, i)
+            {
+                const label edgeI = pEdges[i];
+
+                if (edgeStat[edgeI] != NONE)
+                {
+                    edgeVecs.append(edges[edgeI].vec(localPoints));
+                    edgeVecs.last() /= mag(edgeVecs.last());
+                }
+            }
+
+            if (mag(edgeVecs[0] & edgeVecs[1]) < minCos)
+            {
+                featurePoints.append(pointI);
+            }
+        }
     }
 
     featurePoints_.transfer(featurePoints);
@@ -230,7 +259,7 @@ void Foam::surfaceFeatures::classifyFeatureAngles
         if (eFaces.size() != 2)
         {
             // Non-manifold. What to do here? Is region edge? external edge?
-            edgeStat[edgeI] = REGION;
+            edgeStat[edgeI] = NONE;
         }
         else
         {
@@ -464,7 +493,7 @@ Foam::surfaceFeatures::surfaceFeatures
 
     if (minLen > 0 || minElems > 0)
     {
-        trimFeatures(minLen, minElems);
+        trimFeatures(minLen, minElems, includedAngle);
     }
 }
 
@@ -585,7 +614,7 @@ Foam::surfaceFeatures::surfaceFeatures
     edgeStat.clear();
     dynFeatEdges.clear();
 
-    setFromStatus(allEdgeStat);
+    setFromStatus(allEdgeStat, GREAT);
 }
 
 
@@ -664,7 +693,7 @@ void Foam::surfaceFeatures::findFeatures
         geometricTestOnly
     );
 
-    setFromStatus(edgeStat);
+    setFromStatus(edgeStat, includedAngle);
 }
 
 
@@ -673,7 +702,8 @@ void Foam::surfaceFeatures::findFeatures
 Foam::labelList Foam::surfaceFeatures::trimFeatures
 (
     const scalar minLen,
-    const label minElems
+    const label minElems,
+    const scalar includedAngle
 )
 {
     // Convert feature edge list to status per edge.
@@ -797,7 +827,7 @@ Foam::labelList Foam::surfaceFeatures::trimFeatures
     }
 
     // Convert back to edge labels
-    setFromStatus(edgeStat);
+    setFromStatus(edgeStat, includedAngle);
 
     return featLines;
 }

src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H

@@ -133,7 +133,11 @@ private:
 
 
         //- Construct feature points where more than 2 feature edges meet
-        void calcFeatPoints(const List<edgeStatus>&);
+        void calcFeatPoints
+        (
+            const List<edgeStatus>& edgeStat,
+            const scalar minCos
+        );
 
         //- Classify the angles of the feature edges
         void classifyFeatureAngles
@@ -288,13 +292,22 @@ public:
 
             //- Delete small sets of edges. Edges are stringed up and any
             //  string of length < minLen (or nElems < minElems) is deleted.
-            labelList trimFeatures(const scalar minLen, const label minElems);
+            labelList trimFeatures
+            (
+                const scalar minLen,
+                const label minElems,
+                const scalar includedAngle
+            );
 
             //- From member feature edges to status per edge.
             List<edgeStatus> toStatus() const;
 
             //- Set from status per edge
-            void setFromStatus(const List<edgeStatus>&);
+            void setFromStatus
+            (
+                const List<edgeStatus>& edgeStat,
+                const scalar includedAngle
+            );
 
 
         // Find