ENH: refinementSurfaces: get normal on nearest

src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C

@@ -26,21 +26,14 @@ License
 #include "refinementFeatures.H"
 #include "Time.H"
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::refinementFeatures::refinementFeatures
+void Foam::refinementFeatures::read
 (
     const objectRegistry& io,
     const PtrList<dictionary>& featDicts
 )
-:
-    PtrList<featureEdgeMesh>(featDicts.size()),
-    levels_(featDicts.size()),
-    edgeTrees_(featDicts.size()),
-    pointTrees_(featDicts.size())
 {
-    // Read features
-
     forAll(featDicts, i)
     {
         const dictionary& dict = featDicts[i];
@@ -75,13 +68,16 @@ Foam::refinementFeatures::refinementFeatures
             << " (" << eMesh.points().size() << " points, "
             << eMesh.edges().size() << " edges)." << endl;
     }
+}
 
 
-    // Search engines
+void Foam::refinementFeatures::buildTrees
-
+(
-    forAll(*this, i)
+    const label featI,
-    {
+    const labelList& featurePoints
-        const featureEdgeMesh& eMesh = operator[](i);
+)
+{
+    const featureEdgeMesh& eMesh = operator[](featI);
     const pointField& points = eMesh.points();
     const edgeList& edges = eMesh.edges();
 
@@ -99,7 +95,7 @@ Foam::refinementFeatures::refinementFeatures
 
     edgeTrees_.set
     (
-            i,
+        featI,
         new indexedOctree<treeDataEdge>
         (
             treeDataEdge
@@ -116,9 +112,43 @@ Foam::refinementFeatures::refinementFeatures
         )
     );
 
+    pointTrees_.set
+    (
+        featI,
+        new indexedOctree<treeDataPoint>
+        (
+            treeDataPoint(points, featurePoints),
+            bb,     // overall search domain
+            8,      // maxLevel
+            10,     // leafsize
+            3.0     // duplicity
+        )
+    );
+}
 
-        // Detect feature points from edges.
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::refinementFeatures::refinementFeatures
+(
+    const objectRegistry& io,
+    const PtrList<dictionary>& featDicts
+)
+:
+    PtrList<featureEdgeMesh>(featDicts.size()),
+    levels_(featDicts.size()),
+    edgeTrees_(featDicts.size()),
+    pointTrees_(featDicts.size())
+{
+    // Read features
+    read(io, featDicts);
+
+    // Search engines
+    forAll(*this, i)
+    {
+        const featureEdgeMesh& eMesh = operator[](i);
         const labelListList& pointEdges = eMesh.pointEdges();
+
         DynamicList<label> featurePoints;
         forAll(pointEdges, pointI)
         {
@@ -129,21 +159,79 @@ Foam::refinementFeatures::refinementFeatures
         }
 
         Info<< "Detected " << featurePoints.size()
-            << " featurePoints out of " << points.size()
+            << " featurePoints out of " << pointEdges.size()
             << " on feature " << eMesh.name() << endl;
 
-        pointTrees_.set
+        buildTrees(i, featurePoints);
+    }
+}
+
+
+Foam::refinementFeatures::refinementFeatures
+(
+    const objectRegistry& io,
+    const PtrList<dictionary>& featDicts,
+    const scalar minCos
+)
+:
+    PtrList<featureEdgeMesh>(featDicts.size()),
+    levels_(featDicts.size()),
+    edgeTrees_(featDicts.size()),
+    pointTrees_(featDicts.size())
+{
+    // Read features
+    read(io, featDicts);
+
+    // Search engines
+    forAll(*this, i)
+    {
+        const featureEdgeMesh& eMesh = operator[](i);
+        const pointField& points = eMesh.points();
+        const edgeList& edges = eMesh.edges();
+        const labelListList& pointEdges = eMesh.pointEdges();
+
+        DynamicList<label> featurePoints;
+        forAll(pointEdges, pointI)
+        {
+            const labelList& pEdges = pointEdges[pointI];
+            if (pEdges.size() > 2)
+            {
+                featurePoints.append(pointI);
+            }
+            else if (pEdges.size() == 2)
+            {
+                // Check the angle
+                const edge& e0 = edges[pEdges[0]];
+                const edge& e1 = edges[pEdges[1]];
+
+                const point& p = points[pointI];
+                const point& p0 = points[e0.otherVertex(pointI)];
+                const point& p1 = points[e1.otherVertex(pointI)];
+
+                vector v0 = p-p0;
+                scalar v0Mag = mag(v0);
+
+                vector v1 = p1-p;
+                scalar v1Mag = mag(v1);
+
+                if
                 (
-            i,
+                    v0Mag > SMALL
-            new indexedOctree<treeDataPoint>
+                 && v1Mag > SMALL
-            (
+                 && ((v0/v0Mag & v1/v1Mag) < minCos)
-                treeDataPoint(points, featurePoints),
-                bb,     // overall search domain
-                8,      // maxLevel
-                10,     // leafsize
-                3.0     // duplicity
                 )
-        );
+                {
+                    featurePoints.append(pointI);
+                }
+            }
+        }
+
+        Info<< "Detected " << featurePoints.size()
+            << " featurePoints out of " << points.size()
+            << " on feature " << eMesh.name()
+            << " when using feature cos " << minCos << endl;
+
+        buildTrees(i, featurePoints);
     }
 }
 

src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,17 +69,33 @@ private:
 
     // Private Member Functions
 
+        //- Read set of feature edge meshes
+        void read(const objectRegistry&, const PtrList<dictionary>&);
+
+        //- Build edge tree and feature point tree
+        void buildTrees(const label, const labelList&);
+
+
 public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from description
         refinementFeatures
         (
             const objectRegistry& io,
             const PtrList<dictionary>& featDicts
         );
 
+        //- Construct from description and do geometric analysis to determine
+        //  feature points
+        refinementFeatures
+        (
+            const objectRegistry& io,
+            const PtrList<dictionary>& featDicts,
+            const scalar minCos
+        );
+
 
     // Member Functions
 

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -917,6 +917,127 @@ void Foam::refinementSurfaces::findNearestIntersection
 }
 
 
+void Foam::refinementSurfaces::findNearestIntersection
+(
+    const labelList& surfacesToTest,
+    const pointField& start,
+    const pointField& end,
+
+    labelList& surface1,
+    List<pointIndexHit>& hit1,
+    labelList& region1,
+    vectorField& normal1,
+
+    labelList& surface2,
+    List<pointIndexHit>& hit2,
+    labelList& region2,
+    vectorField& normal2
+) const
+{
+    // 1. intersection from start to end
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Initialize arguments
+    surface1.setSize(start.size());
+    surface1 = -1;
+    hit1.setSize(start.size());
+    region1.setSize(start.size());
+    normal1.setSize(start.size());
+
+    // Current end of segment to test.
+    pointField nearest(end);
+    // Work array
+    List<pointIndexHit> nearestInfo(start.size());
+    labelList region;
+    vectorField normal;
+
+    forAll(surfacesToTest, testI)
+    {
+        label surfI = surfacesToTest[testI];
+        const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+
+        // See if any intersection between start and current nearest
+        geom.findLine(start, nearest, nearestInfo);
+        geom.getRegion(nearestInfo, region);
+        geom.getNormal(nearestInfo, normal);
+
+        forAll(nearestInfo, pointI)
+        {
+            if (nearestInfo[pointI].hit())
+            {
+                hit1[pointI] = nearestInfo[pointI];
+                surface1[pointI] = surfI;
+                region1[pointI] = region[pointI];
+                normal1[pointI] = normal[pointI];
+                nearest[pointI] = hit1[pointI].hitPoint();
+            }
+        }
+    }
+
+
+    // 2. intersection from end to last intersection
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Find the nearest intersection from end to start. Note that we initialize
+    // to the first intersection (if any).
+    surface2 = surface1;
+    hit2 = hit1;
+    region2 = region1;
+    normal2 = normal1;
+
+    // Set current end of segment to test.
+    forAll(nearest, pointI)
+    {
+        if (hit1[pointI].hit())
+        {
+            nearest[pointI] = hit1[pointI].hitPoint();
+        }
+        else
+        {
+            // Disable testing by setting to end.
+            nearest[pointI] = end[pointI];
+        }
+    }
+
+    forAll(surfacesToTest, testI)
+    {
+        label surfI = surfacesToTest[testI];
+        const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+
+        // See if any intersection between end and current nearest
+        geom.findLine(end, nearest, nearestInfo);
+        geom.getRegion(nearestInfo, region);
+        geom.getNormal(nearestInfo, normal);
+
+        forAll(nearestInfo, pointI)
+        {
+            if (nearestInfo[pointI].hit())
+            {
+                hit2[pointI] = nearestInfo[pointI];
+                surface2[pointI] = surfI;
+                region2[pointI] = region[pointI];
+                normal2[pointI] = normal[pointI];
+                nearest[pointI] = hit2[pointI].hitPoint();
+            }
+        }
+    }
+
+
+    // Make sure that if hit1 has hit something, hit2 will have at least the
+    // same point (due to tolerances it might miss its end point)
+    forAll(hit1, pointI)
+    {
+        if (hit1[pointI].hit() && !hit2[pointI].hit())
+        {
+            hit2[pointI] = hit1[pointI];
+            surface2[pointI] = surface1[pointI];
+            region2[pointI] = region1[pointI];
+            normal2[pointI] = normal1[pointI];
+        }
+    }
+}
+
+
 void Foam::refinementSurfaces::findAnyIntersection
 (
     const pointField& start,

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -317,6 +317,24 @@ public:
                 labelList& region2
             ) const;
 
+            //- findNearestIntersection but also get normals
+            void findNearestIntersection
+            (
+                const labelList& surfacesToTest,
+                const pointField& start,
+                const pointField& end,
+
+                labelList& surface1,
+                List<pointIndexHit>& hit1,
+                labelList& region1,
+                vectorField& normal1,
+
+                labelList& surface2,
+                List<pointIndexHit>& hit2,
+                labelList& region2,
+                vectorField& normal2
+            ) const;
+
             //- Used for debugging only: find intersection of edge.
             void findAnyIntersection
             (