ENH: snappyHexMesh: allow medialAxis to shrink 1 or 2 cell thick channel. Fixes #986.

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -665,6 +665,13 @@ addLayersControls
         // Default is 0.5*featureAngle. Set to -180 always attempt extrusion
         //layerTerminationAngle 25;
 
+        // Optional: disable shrinking of edges that have one (or two) points
+        // on an extruded patch.
+        // Default is false to enable single/two cell thick channels to still
+        // have layers. In <=1806 this was true by default. On larger gaps it
+        // should have no effect.
+        //disableWallEdges true;
+
         // Optional: at non-patched sides allow mesh to slip if extrusion
         // direction makes angle larger than slipFeatureAngle. Default is
         // 0.5*featureAngle.

src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C

@@ -58,18 +58,18 @@ bool Foam::medialAxisMeshMover::isMaxEdge
 (
     const List<PointData<vector>>& pointWallDist,
     const label edgeI,
-    const scalar minCos
+    const scalar minCos,
+    const bool disableWallEdges
 ) const
 {
     const pointField& points = mesh().points();
-
-    // Do not mark edges with one side on moving wall.
-
     const edge& e = mesh().edges()[edgeI];
 
+    if (disableWallEdges)
+    {
+        // 1. Do not mark edges with one side on moving wall.
         vector v0(points[e[0]] - pointWallDist[e[0]].origin());
         scalar magV0(mag(v0));
-
         if (magV0 < SMALL)
         {
             return false;
@@ -77,14 +77,23 @@ bool Foam::medialAxisMeshMover::isMaxEdge
 
         vector v1(points[e[1]] - pointWallDist[e[1]].origin());
         scalar magV1(mag(v1));
-
         if (magV1 < SMALL)
         {
             return false;
         }
+    }
 
+    //// 2. Do not mark edges with both sides on a moving wall.
+    //vector v0(points[e[0]] - pointWallDist[e[0]].origin());
+    //scalar magV0(mag(v0));
+    //vector v1(points[e[1]] - pointWallDist[e[1]].origin());
+    //scalar magV1(mag(v1));
+    //if (magV0 < SMALL && magV1 < SMALL)
+    //{
+    //    return false;
+    //}
 
-    //- Detect based on vector to nearest point differing for both endpoints
+    //// 3. Detect based on vector to nearest point differing for both endpoints
     //v0 /= magV0;
     //v1 /= magV1;
     //
@@ -170,6 +179,13 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
         mesh().globalData().nTotalPoints()
     );
 
+    const bool disableWallEdges = coeffDict.lookupOrDefault<bool>
+    (
+        "disableWallEdges",
+        false
+    );
+
+
 
     // Predetermine mesh edges
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -336,7 +352,16 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
                 }
 
             }
-            else if (isMaxEdge(pointWallDist, edgeI, minMedialAxisAngleCos))
+            else if
+            (
+                isMaxEdge
+                (
+                    pointWallDist,
+                    edgeI,
+                    minMedialAxisAngleCos,
+                    disableWallEdges
+                )
+            )
             {
                 // Both end points of edge have very different nearest wall
                 // point. Mark both points as medial axis points.
@@ -351,25 +376,39 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
 
                     // Calculate distance along edge
                     const point& p0 = points[e[0]];
+                    const point& origin0 = pointWallDist[e[0]].origin();
                     const point& p1 = points[e[1]];
-                    scalar dist0 = (p0-pointWallDist[e[0]].origin()) & eVec;
+                    const point& origin1 = pointWallDist[e[1]].origin();
-                    scalar dist1 = (pointWallDist[e[1]].origin()-p1) & eVec;
+                    scalar dist0 = (p0-origin0) & eVec;
+                    scalar dist1 = (origin1-p1) & eVec;
                     scalar s = 0.5*(dist1+eMag+dist0);
 
-                    point medialAxisPt;
+                    point medialAxisPt(vector::max);
                     if (s <= dist0)
+                    {
+                        // Make sure point is not on wall. Note that this
+                        // check used to be inside isMaxEdge.
+                        if (magSqr((p0-origin0)) > Foam::sqr(SMALL))
                         {
                             medialAxisPt = p0;
                         }
+                    }
                     else if (s >= dist0+eMag)
+                    {
+                        // Make sure point is not on wall. Note that this
+                        // check used to be inside isMaxEdge.
+                        if (magSqr((p1-origin1)) > Foam::sqr(SMALL))
                         {
                             medialAxisPt = p1;
                         }
+                    }
                     else
                     {
                         medialAxisPt = p0+(s-dist0)*eVec;
                     }
 
+                    if (medialAxisPt != vector::max)
+                    {
                         forAll(e, ep)
                         {
                             label pointI = e[ep];
@@ -391,6 +430,7 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
                     }
                 }
             }
+        }
 
 
         // 2. Seed non-adapt patches

src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H

@@ -114,7 +114,8 @@ class medialAxisMeshMover
             (
                 const List<PointData<vector>>& pointWallDist,
                 const label edgeI,
-                const scalar minCos
+                const scalar minCos,
+                const bool checkWall
             ) const;
 
             //- Initialise medial axis. Uses pointDisplacement field only