From 357939a77052e6410c1a14688705ceead6a56486 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 29 Oct 2012 13:20:42 +0000
Subject: [PATCH] ENH: snappyHexMesh: refine by distance to feature

---
 .../snappyHexMesh/snappyHexMeshDict           |  13 +-
 .../autoHexMeshDriver/autoRefineDriver.C      |   3 +
 .../meshRefinement/meshRefinement.H           |   9 +
 .../meshRefinement/meshRefinementRefine.C     | 119 +++++++++++-
 .../refinementFeatures/refinementFeatures.C   | 180 +++++++++++++++++-
 .../refinementFeatures/refinementFeatures.H   |  34 +++-
 6 files changed, 336 insertions(+), 22 deletions(-)

diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index f0cc40a8e4..4bef1736f0 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -96,20 +96,29 @@ castellatedMeshControls
     // refinement.
     nCellsBetweenLevels 1;
 
+
     // Explicit feature edge refinement
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Specifies a level for any cell intersected by explicitly provided
     // edges.
     // This is a featureEdgeMesh, read from constant/triSurface for now.
+    // Specify 'levels' in the same way as the 'distance' mode in the
+    // refinementRegions (see below). The old specification
+    //      level   2;
+    // is equivalent to
+    //      levels  ((0 2));
+
     features
     (
         //{
         //    file "someLine.eMesh";
-        //    level 2;
+        //    //level 2;
+        //    levels ((0.0 2) (1.0 3));
         //}
     );
 
+
     // Surface based refinement
     // ~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -178,7 +187,7 @@ castellatedMeshControls
     // three modes
     // - distance. 'levels' specifies per distance to the surface the
     //   wanted refinement level. The distances need to be specified in
-    //   descending order.
+    //   increasing order.
     // - inside. 'levels' is only one entry and only the level is used. All
     //   cells inside the surface get refined up to the level. The surface
     //   needs to be closed for this to be possible.
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index 985bcfd103..dfd769fd15 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -97,6 +97,7 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
                     refineParams.curvature(),
 
                     true,               // featureRefinement
+                    false,              // featureDistanceRefinement
                     false,              // internalRefinement
                     false,              // surfaceRefinement
                     false,              // curvatureRefinement
@@ -207,6 +208,7 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
                 refineParams.curvature(),
 
                 false,              // featureRefinement
+                false,              // featureDistanceRefinement
                 false,              // internalRefinement
                 true,               // surfaceRefinement
                 true,               // curvatureRefinement
@@ -368,6 +370,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
                 refineParams.curvature(),
 
                 false,              // featureRefinement
+                true,               // featureDistanceRefinement
                 true,               // internalRefinement
                 false,              // surfaceRefinement
                 false,              // curvatureRefinement
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 918aae9000..b98230e412 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -264,6 +264,14 @@ private:
                 label& nRefine
             ) const;
 
+            //- Mark cells for distance-to-feature based refinement.
+            label markInternalDistanceToFeatureRefinement
+            (
+                const label nAllowRefine,
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
             //- Mark cells for refinement-shells based refinement.
             label markInternalRefinement
             (
@@ -656,6 +664,7 @@ public:
                 const scalar curvature,
 
                 const bool featureRefinement,
+                const bool featureDistanceRefinement,
                 const bool internalRefinement,
                 const bool surfaceRefinement,
                 const bool curvatureRefinement,
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index 4d9d10ab80..2c15a847ea 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -290,7 +290,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
         forAll(features_, featI)
         {
             const featureEdgeMesh& featureMesh = features_[featI];
-            const label featureLevel = features_.levels()[featI];
+            const label featureLevel = features_.levels()[featI][0];
             const labelListList& pointEdges = featureMesh.pointEdges();
 
             // Find regions on edgeMesh
@@ -511,6 +511,90 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
 }
 
 
+// Mark cells for distance-to-feature based refinement.
+Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
+(
+    const label nAllowRefine,
+
+    labelList& refineCell,
+    label& nRefine
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const pointField& cellCentres = mesh_.cellCentres();
+
+    // Detect if there are any distance shells
+    if (features_.maxDistance() <= 0.0)
+    {
+        return 0;
+    }
+
+    label oldNRefine = nRefine;
+
+    // Collect cells to test
+    pointField testCc(cellLevel.size()-nRefine);
+    labelList testLevels(cellLevel.size()-nRefine);
+    label testI = 0;
+
+    forAll(cellLevel, cellI)
+    {
+        if (refineCell[cellI] == -1)
+        {
+            testCc[testI] = cellCentres[cellI];
+            testLevels[testI] = cellLevel[cellI];
+            testI++;
+        }
+    }
+
+    // Do test to see whether cells is inside/outside shell with higher level
+    labelList maxLevel;
+    features_.findHigherLevel(testCc, testLevels, maxLevel);
+
+    // Mark for refinement. Note that we didn't store the original cellID so
+    // now just reloop in same order.
+    testI = 0;
+    forAll(cellLevel, cellI)
+    {
+        if (refineCell[cellI] == -1)
+        {
+            if (maxLevel[testI] > testLevels[testI])
+            {
+                bool reachedLimit = !markForRefine
+                (
+                    maxLevel[testI],    // mark with any positive value
+                    nAllowRefine,
+                    refineCell[cellI],
+                    nRefine
+                );
+
+                if (reachedLimit)
+                {
+                    if (debug)
+                    {
+                        Pout<< "Stopped refining internal cells"
+                            << " since reaching my cell limit of "
+                            << mesh_.nCells()+7*nRefine << endl;
+                    }
+                    break;
+                }
+            }
+            testI++;
+        }
+    }
+
+    if
+    (
+        returnReduce(nRefine, sumOp<label>())
+      > returnReduce(nAllowRefine, sumOp<label>())
+    )
+    {
+        Info<< "Reached refinement limit." << endl;
+    }
+
+    return returnReduce(nRefine-oldNRefine, sumOp<label>());
+}
+
+
 // Mark cells for non-surface intersection based refinement.
 Foam::label Foam::meshRefinement::markInternalRefinement
 (
@@ -1128,6 +1212,7 @@ Foam::labelList Foam::meshRefinement::refineCandidates
     const scalar curvature,
 
     const bool featureRefinement,
+    const bool featureDistanceRefinement,
     const bool internalRefinement,
     const bool surfaceRefinement,
     const bool curvatureRefinement,
@@ -1196,8 +1281,24 @@ Foam::labelList Foam::meshRefinement::refineCandidates
                 nRefine
             );
 
-            Info<< "Marked for refinement due to explicit features    : "
-                << nFeatures << " cells."  << endl;
+            Info<< "Marked for refinement due to explicit features             "
+                << ": " << nFeatures << " cells."  << endl;
+        }
+
+        // Inside distance-to-feature shells
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        if (featureDistanceRefinement)
+        {
+            label nShell = markInternalDistanceToFeatureRefinement
+            (
+                nAllowRefine,
+
+                refineCell,
+                nRefine
+            );
+            Info<< "Marked for refinement due to distance to explicit features "
+                ": " << nShell << " cells."  << endl;
         }
 
         // Inside refinement shells
@@ -1212,8 +1313,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
                 refineCell,
                 nRefine
             );
-            Info<< "Marked for refinement due to refinement shells    : "
-                << nShell << " cells."  << endl;
+            Info<< "Marked for refinement due to refinement shells             "
+                << ": " << nShell << " cells."  << endl;
         }
 
         // Refinement based on intersection of surface
@@ -1230,8 +1331,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
                 refineCell,
                 nRefine
             );
-            Info<< "Marked for refinement due to surface intersection : "
-                << nSurf << " cells."  << endl;
+            Info<< "Marked for refinement due to surface intersection          "
+                << ": " << nSurf << " cells."  << endl;
         }
 
         // Refinement based on curvature of surface
@@ -1254,8 +1355,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
                 refineCell,
                 nRefine
             );
-            Info<< "Marked for refinement due to curvature/regions    : "
-                << nCurv << " cells."  << endl;
+            Info<< "Marked for refinement due to curvature/regions             "
+                << ": " << nCurv << " cells."  << endl;
         }
 
         // Pack cells-to-refine
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
index 0cdd440d1a..e72f13d726 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
@@ -25,6 +25,7 @@ License
 
 #include "refinementFeatures.H"
 #include "Time.H"
+#include "Tuple2.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -34,15 +35,15 @@ void Foam::refinementFeatures::read
     const PtrList<dictionary>& featDicts
 )
 {
-    forAll(featDicts, i)
+    forAll(featDicts, featI)
     {
-        const dictionary& dict = featDicts[i];
+        const dictionary& dict = featDicts[featI];
 
         fileName featFileName(dict.lookup("file"));
 
         set
         (
-            i,
+            featI,
             new featureEdgeMesh
             (
                 IOobject
@@ -58,15 +59,74 @@ void Foam::refinementFeatures::read
             )
         );
 
-        const featureEdgeMesh& eMesh = operator[](i);
+        const featureEdgeMesh& eMesh = operator[](featI);
 
         //eMesh.mergePoints(meshRefiner_.mergeDistance());
-        levels_[i] = readLabel(dict.lookup("level"));
 
-        Info<< "Refinement level " << levels_[i]
-            << " for all cells crossed by feature " << featFileName
-            << " (" << eMesh.points().size() << " points, "
+        if (dict.found("levels"))
+        {
+            List<Tuple2<scalar, label> > distLevels(dict["levels"]);
+
+            if (dict.size() < 1)
+            {
+                FatalErrorIn
+                (
+                    "refinementFeatures::read"
+                    "(const objectRegistry&"
+                    ", const PtrList<dictionary>&)"
+                )   << " : levels should be at least size 1" << endl
+                    << "levels : "  << dict["levels"]
+                    << exit(FatalError);
+            }
+
+            distances_[featI].setSize(distLevels.size());
+            levels_[featI].setSize(distLevels.size());
+
+            forAll(distLevels, j)
+            {
+                distances_[featI][j] = distLevels[j].first();
+                levels_[featI][j] = distLevels[j].second();
+
+                // Check in incremental order
+                if (j > 0)
+                {
+                    if
+                    (
+                        (distances_[featI][j] <= distances_[featI][j-1])
+                     || (levels_[featI][j] > levels_[featI][j-1])
+                    )
+                    {
+                        FatalErrorIn
+                        (
+                            "refinementFeatures::read"
+                            "(const objectRegistry&"
+                            ", const PtrList<dictionary>&)"
+                        )   << " : Refinement should be specified in order"
+                            << " of increasing distance"
+                            << " (and decreasing refinement level)." << endl
+                            << "Distance:" << distances_[featI][j]
+                            << " refinementLevel:" << levels_[featI][j]
+                            << exit(FatalError);
+                    }
+                }
+            }
+        }
+        else
+        {
+            // Look up 'level' for single level
+            levels_[featI] = labelList(1, readLabel(dict.lookup("level")));
+            distances_[featI] = scalarField(1, 0.0);
+        }
+
+        Info<< "Refinement level according to distance to "
+            << featFileName << " (" << eMesh.points().size() << " points, "
             << eMesh.edges().size() << " edges)." << endl;
+        forAll(levels_[featI], j)
+        {
+            Info<< "    level " << levels_[featI][j]
+                << " for all cells within " << distances_[featI][j]
+                << " meter." << endl;
+        }
     }
 }
 
@@ -127,6 +187,80 @@ void Foam::refinementFeatures::buildTrees
 }
 
 
+
+// Find maximum level of a shell.
+void Foam::refinementFeatures::findHigherLevel
+(
+    const pointField& pt,
+    const label featI,
+    labelList& maxLevel
+) const
+{
+    const labelList& levels = levels_[featI];
+
+    const scalarField& distances = distances_[featI];
+
+    // Collect all those points that have a current maxLevel less than
+    // (any of) the shell. Also collect the furthest distance allowable
+    // to any shell with a higher level.
+
+    pointField candidates(pt.size());
+    labelList candidateMap(pt.size());
+    scalarField candidateDistSqr(pt.size());
+    label candidateI = 0;
+
+    forAll(maxLevel, pointI)
+    {
+        forAllReverse(levels, levelI)
+        {
+            if (levels[levelI] > maxLevel[pointI])
+            {
+                candidates[candidateI] = pt[pointI];
+                candidateMap[candidateI] = pointI;
+                candidateDistSqr[candidateI] = sqr(distances[levelI]);
+                candidateI++;
+                break;
+            }
+        }
+    }
+    candidates.setSize(candidateI);
+    candidateMap.setSize(candidateI);
+    candidateDistSqr.setSize(candidateI);
+
+    // Do the expensive nearest test only for the candidate points.
+    const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
+
+    List<pointIndexHit> nearInfo(candidates.size());
+    forAll(candidates, candidateI)
+    {
+        nearInfo[candidateI] = tree.findNearest
+        (
+            candidates[candidateI],
+            candidateDistSqr[candidateI]
+        );
+    }
+
+    // Update maxLevel
+    forAll(nearInfo, candidateI)
+    {
+        if (nearInfo[candidateI].hit())
+        {
+            // Check which level it actually is in.
+            label minDistI = findLower
+            (
+                distances,
+                mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
+            );
+
+            label pointI = candidateMap[candidateI];
+
+            // pt is inbetween shell[minDistI] and shell[minDistI+1]
+            maxLevel[pointI] = levels[minDistI+1];
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::refinementFeatures::refinementFeatures
@@ -136,6 +270,7 @@ Foam::refinementFeatures::refinementFeatures
 )
 :
     PtrList<featureEdgeMesh>(featDicts.size()),
+    distances_(featDicts.size()),
     levels_(featDicts.size()),
     edgeTrees_(featDicts.size()),
     pointTrees_(featDicts.size())
@@ -175,6 +310,7 @@ Foam::refinementFeatures::refinementFeatures
 )
 :
     PtrList<featureEdgeMesh>(featDicts.size()),
+    distances_(featDicts.size()),
     levels_(featDicts.size()),
     edgeTrees_(featDicts.size()),
     pointTrees_(featDicts.size())
@@ -336,4 +472,32 @@ void Foam::refinementFeatures::findNearestPoint
 }
 
 
+void Foam::refinementFeatures::findHigherLevel
+(
+    const pointField& pt,
+    const labelList& ptLevel,
+    labelList& maxLevel
+) const
+{
+    // Maximum level of any shell. Start off with level of point.
+    maxLevel = ptLevel;
+
+    forAll(*this, featI)
+    {
+        findHigherLevel(pt, featI, maxLevel);
+    }
+}
+
+
+Foam::scalar Foam::refinementFeatures::maxDistance() const
+{
+    scalar overallMax = -GREAT;
+    forAll(distances_, featI)
+    {
+        overallMax = max(overallMax, max(distances_[featI]));
+    }
+    return overallMax;
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
index 5442e140dd..aa2807bac1 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
@@ -57,8 +57,11 @@ private:
 
     // Private data
 
-        //- Refinement levels
-        labelList levels_;
+        //- Per shell the list of ranges
+        List<scalarField> distances_;
+
+        //- Per shell per distance the refinement level
+        labelListList levels_;
 
         //- Edge
         PtrList<indexedOctree<treeDataEdge> > edgeTrees_;
@@ -75,6 +78,13 @@ private:
         //- Build edge tree and feature point tree
         void buildTrees(const label, const labelList&);
 
+        //- Find shell level higher than ptLevel
+        void findHigherLevel
+        (
+            const pointField& pt,
+            const label featI,
+            labelList& maxLevel
+        ) const;
 
 public:
 
@@ -101,11 +111,18 @@ public:
 
         // Access
 
-            const labelList& levels() const
+            //- Per featureEdgeMesh the list of level
+            const labelListList& levels() const
             {
                 return levels_;
             }
 
+            //- Per featureEdgeMesh the list of ranges
+            const List<scalarField>& distances() const
+            {
+                return distances_;
+            }
+
             const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
             {
                 return edgeTrees_;
@@ -119,6 +136,9 @@ public:
 
         // Query
 
+            //- Highest distance of all features
+            scalar maxDistance() const;
+
             //- Find nearest point on nearest feature edge
             void findNearestEdge
             (
@@ -141,6 +161,14 @@ public:
                 labelList& nearIndex
             ) const;
 
+            //- Find shell level higher than ptLevel
+            void findHigherLevel
+            (
+                const pointField& pt,
+                const labelList& ptLevel,
+                labelList& maxLevel
+            ) const;
+
 };