From 14a575e8a325cdd95ddf91d2be12a39ce449185e Mon Sep 17 00:00:00 2001
From: graham <graham.macpherson@strath.ac.uk>
Date: Wed, 27 Aug 2008 18:58:33 +0100
Subject: [PATCH] Using surfaceFeatures to identify and store feature edges,
 planning to use nearestSurfEdge to identify where to place feature edge
 constraining points.

---
 .../mesh/generation/CV3DMesher/CV3D.H         |  4 +-
 .../mesh/generation/CV3DMesher/CV3DMesher.C   | 32 ++++++++-
 .../mesh/generation/CV3DMesher/controls.C     |  2 +-
 .../CV3DMesher/insertFeaturePoints.C          | 15 ++--
 .../insertSurfaceNearestPointPairs.C          |  2 +-
 .../mesh/generation/CV3DMesher/querySurface.C | 69 ++++++-------------
 .../mesh/generation/CV3DMesher/querySurface.H | 23 ++++---
 7 files changed, 75 insertions(+), 72 deletions(-)

diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3D.H b/applications/utilities/mesh/generation/CV3DMesher/CV3D.H
index 226f3bea5d..92dd5c614a 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/CV3D.H
+++ b/applications/utilities/mesh/generation/CV3DMesher/CV3D.H
@@ -85,8 +85,8 @@ public:
 
         //- The feature angle used to select corners to be
         //  explicitly represented in the mesh.
-        //  0 = all features, 180 = no features
-        scalar featAngle;
+        //  0 = no features, 180 = all features
+        scalar includedAngle;
 
         //- Maximum quadrant angle allowed at a concave corner before
         //  additional "mitering" lines are added
diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C b/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C
index ffa56cfb24..0dd342ebb7 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C
@@ -61,12 +61,42 @@ int main(int argc, char *argv[])
 
     // Read the surface to conform to
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    querySurface surf(args.args()[1]);
+    querySurface surf
+    (
+        args.args()[1],
+        readScalar(controlDict.lookup("includedAngle"))
+    );
 
     Info<< nl
         << "Read surface with " << surf.size() << " triangles from file "
         << args.args()[1] << nl << endl;
 
+    surf.features().writeObj("features");
+
+    pointField samples(3);
+
+    samples[0] = point(0.9, 0.9, 0.9);
+    samples[1] = point(0.85, 0.2, 0.2);
+    samples[2] = point(-0.1, -0.1, 0.25);
+
+    labelList edgeLabel;
+    labelList edgeEndPoint;
+    pointField edgePoint;
+
+    surf.features().nearestSurfEdge
+    (
+        surf.features().featureEdges(),
+        samples,
+        vector(0.25, 0.25, 0.25),
+        edgeLabel,
+        edgeEndPoint,
+        edgePoint
+    );
+
+    Info<< edgeLabel
+        << nl << edgeEndPoint
+        << nl << edgePoint << endl;
+
     // Read and triangulation
     // ~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/applications/utilities/mesh/generation/CV3DMesher/controls.C b/applications/utilities/mesh/generation/CV3DMesher/controls.C
index 3365baf61b..c133415451 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/controls.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/controls.C
@@ -32,7 +32,7 @@ Foam::CV3D::controls::controls(const dictionary& controlDict)
 :
     minCellSize(readScalar(controlDict.lookup("minCellSize"))),
     minCellSize2(Foam::sqr(minCellSize)),
-    featAngle(readScalar(controlDict.lookup("featureAngle"))),
+    includedAngle(readScalar(controlDict.lookup("includedAngle"))),
     maxQuadAngle(readScalar(controlDict.lookup("maxQuadAngle"))),
     insertSurfaceNearestPointPairs
     (
diff --git a/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C b/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
index b09aa20cdb..3146c9c38b 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
@@ -38,17 +38,10 @@ void Foam::CV3D::insertFeaturePoints()
     const edgeList& edges = qSurf_.edges();
     const pointField& localPts = qSurf_.localPoints();
 
-    labelList featPoints(0);
-    labelListList featPointFeatEdges(0);
+    const labelList& featPoints = qSurf_.features().featurePoints();
+    labelListList featPointFeatEdges = qSurf_.featurePointFeatureEdges();
 
-    qSurf_.extractFeatures
-    (
-        controls_.featAngle,
-        featPoints,
-        featPointFeatEdges
-    );
-
-    scalar planeErrorAngle = 0.1*controls_.featAngle;
+    scalar planeErrorAngle = 0.1*(180.0 - controls_.includedAngle);
 
     scalar planeErrorAngleCos = cos(mathematicalConstant::pi*planeErrorAngle/180.0);
 
@@ -57,7 +50,7 @@ void Foam::CV3D::insertFeaturePoints()
         label ptI = featPoints[i];
         const point& featPt = localPts[ptI];
 
-        Info<< nl <<"Feature at " << featPt << endl;;
+        Info<< nl <<"Feature at " << featPt << endl;
 
         const labelList& featEdges = featPointFeatEdges[i];
 
diff --git a/applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C b/applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C
index 940b1d6e0b..2a8771c588 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C
@@ -207,7 +207,7 @@ void Foam::CV3D::insertSurfaceNearestPointPairs()
                     }
                 }
 
-               if (!internalFeatureEdge && dualCellSurfaceIntersection(vit))
+                if (!internalFeatureEdge && dualCellSurfaceIntersection(vit))
                 {
                     nearSurfacePoints.append(vert);
                     surfacePoints.append(pHit.hitPoint());
diff --git a/applications/utilities/mesh/generation/CV3DMesher/querySurface.C b/applications/utilities/mesh/generation/CV3DMesher/querySurface.C
index 4575abc1cb..66a2096792 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/querySurface.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/querySurface.C
@@ -29,7 +29,11 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::querySurface::querySurface(const fileName& surfaceFileName)
+Foam::querySurface::querySurface
+(
+    const fileName& surfaceFileName,
+    const scalar& includedAngle
+)
 :
     triSurface(surfaceFileName),
     rndGen_(12345),
@@ -41,7 +45,8 @@ Foam::querySurface::querySurface(const fileName& surfaceFileName)
         8,                         // maxLevel
         4, //10,                   // leafsize
         10.0 //3.0                 // duplicity
-    )
+    ),
+    sFeat_(*this, includedAngle)
 {}
 
 
@@ -53,71 +58,41 @@ Foam::querySurface::~querySurface()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::querySurface::extractFeatures
-(
-    const scalar featAngle,
-    labelList& featPoints,
-    labelListList& featPointFeatEdges
-) const
+Foam::labelListList Foam::querySurface::featurePointFeatureEdges() const
 {
-    scalar featCos = cos(mathematicalConstant::pi*featAngle/180.0);
+    const labelList& featPoints = features().featurePoints();
 
-    const labelListList& edgeFaces = this->edgeFaces();
-    const labelListList& pointEdges = this->pointEdges();
-    const pointField& localPoints = this->localPoints();
-    // const edgeList& edges = this->edges();
-    const vectorField& faceNormals = this->faceNormals();
+    const labelList& featEdges = features().featureEdges();
 
-    DynamicList<label> tempFeatPoints(localPoints.size());
+    const edgeList& edges = this->edges();
 
-    DynamicList<DynamicList<label> > tempFeatPointFeatEdges(localPoints.size());
+    List<DynamicList<label> > tempFeatPointFeatEdges(featPoints.size());
 
-    forAll(pointEdges, ptI)
+    forAll(featPoints, pI)
     {
-        // const point& p = localPoints[ptI];
+        label fP = featPoints[pI];
 
-        // const vector& pNormal = pointNormals[ptI];
-
-        const labelList& pEdges = pointEdges[ptI];
-
-        DynamicList<label> tempFeatEdges(0);
-
-        forAll(pEdges, pE)
+        forAll(featEdges, eI)
         {
-            label edgeI = pEdges[pE];
+            label fE = featEdges[eI];
 
-            const labelList& eFaces = edgeFaces[edgeI];
+            const edge& e(edges[fE]);
 
-            if
-            (
-                eFaces.size() == 2
-                && (faceNormals[eFaces[0]] & faceNormals[eFaces[1]]) < featCos
-            )
+            if (e.start() == fP || e.end() == fP)
             {
-                tempFeatEdges.append(edgeI);
+                tempFeatPointFeatEdges[pI].append(fE);
             }
         }
-
-        tempFeatEdges.shrink();
-
-        if(tempFeatEdges.size() > 2)
-        {
-            tempFeatPoints.append(ptI);
-
-            tempFeatPointFeatEdges.append(tempFeatEdges);
-        }
     }
 
-    featPoints.transfer(tempFeatPoints.shrink());
-
-    tempFeatPointFeatEdges.shrink();
-
-    featPointFeatEdges.setSize(tempFeatPointFeatEdges.size());
+    labelListList featPointFeatEdges(tempFeatPointFeatEdges.size());
 
     forAll(featPointFeatEdges, fPFE)
     {
         featPointFeatEdges[fPFE].transfer(tempFeatPointFeatEdges[fPFE].shrink());
     }
+
+    return featPointFeatEdges;
 }
 
 Foam::indexedOctree<Foam::treeDataTriSurface>::volumeType
diff --git a/applications/utilities/mesh/generation/CV3DMesher/querySurface.H b/applications/utilities/mesh/generation/CV3DMesher/querySurface.H
index 44e89380bf..71dfd2354c 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/querySurface.H
+++ b/applications/utilities/mesh/generation/CV3DMesher/querySurface.H
@@ -39,6 +39,7 @@ SourceFiles
 #include "triSurface.H"
 #include "treeDataTriSurface.H"
 #include "indexedOctree.H"
+#include "surfaceFeatures.H"
 #include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -64,6 +65,8 @@ class querySurface
         // Search engine on surface
         indexedOctree<treeDataTriSurface> tree_;
 
+        surfaceFeatures sFeat_;
+
 
     // Private Member Functions
 
@@ -79,7 +82,11 @@ public:
     // Constructors
 
         //- Construct given file name of the surface
-        querySurface(const fileName& surfaceFileName);
+        querySurface
+        (
+            const fileName& surfaceFileName,
+            const scalar& includedAngle
+        );
 
 
     // Destructor
@@ -101,17 +108,15 @@ public:
                 return tree_;
             }
 
+            const surfaceFeatures& features() const
+            {
+                return sFeat_;
+            }
+
 
         // Query
 
-            //- Extract feature edges/points
-            //  using the given feature angle in deg
-            void extractFeatures
-            (
-                const scalar featAngle,
-                labelList& featPoints,
-                labelListList& featPointFeatEdges
-            ) const;
+            labelListList featurePointFeatureEdges() const;
 
             //- Returns inside, outside or mixed (= close to surface)
             indexedOctree<Foam::treeDataTriSurface>::volumeType insideOutside