ENH: snappyHexMesh: allow sliding

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -710,7 +710,8 @@ void Foam::autoLayerDriver::setNumLayers
     const indirectPrimitivePatch& pp,
     pointField& patchDisp,
     labelList& patchNLayers,
-    List<extrudeMode>& extrudeStatus
+    List<extrudeMode>& extrudeStatus,
+    label& nAddedCells
 ) const
 {
     const fvMesh& mesh = meshRefiner_.mesh();
@@ -797,6 +798,24 @@ void Foam::autoLayerDriver::setNumLayers
         }
     }
 
+
+    // Calculate number of cells to create
+    nAddedCells = 0;
+    forAll(pp.localFaces(), faceI)
+    {
+        const face& f = pp.localFaces()[faceI];
+
+        // Get max of extrusion per point
+        label nCells = 0;
+        forAll(f, fp)
+        {
+            nCells = max(nCells, patchNLayers[f[fp]]);
+        }
+
+        nAddedCells += nCells;
+    }
+    reduce(nAddedCells, sumOp<label>());
+
     //reduce(nConflicts, sumOp<label>());
     //
     //Info<< "Set displacement to zero for " << nConflicts
@@ -2491,9 +2510,13 @@ void Foam::autoLayerDriver::addLayers
     // extrudeStatus = EXTRUDE.
     labelList patchNLayers(pp().nPoints(), 0);
 
+    // Ideal number of cells added
+    label nIdealAddedCells = 0;
+
     // Whether to add edge for all pp.localPoints.
     List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
 
+
     {
         // Get number of layer per point from number of layers per patch
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -2506,7 +2529,8 @@ void Foam::autoLayerDriver::addLayers
 
             patchDisp,
             patchNLayers,
-            extrudeStatus
+            extrudeStatus,
+            nIdealAddedCells
         );
 
         // Precalculate mesh edge labels for patch edges
@@ -2776,6 +2800,7 @@ void Foam::autoLayerDriver::addLayers
         layerParams.nSmoothNormals(),
         layerParams.nSmoothSurfaceNormals(),
         layerParams.minMedianAxisAngleCos(),
+        layerParams.featureAngle(),
 
         dispVec,
         medialRatio,
@@ -3101,10 +3126,24 @@ void Foam::autoLayerDriver::addLayers
 
         label nExtruded = countExtrusion(pp, extrudeStatus);
         label nTotFaces = returnReduce(pp().size(), sumOp<label>());
+        label nAddedCells = 0;
+        {
+            forAll(flaggedCells, cellI)
+            {
+                if (flaggedCells[cellI])
+                {
+                    nAddedCells++;
+                }
+            }
+            reduce(nAddedCells, sumOp<label>());
+        }
         Info<< "Extruding " << nExtruded
             << " out of " << nTotFaces
             << " faces (" << 100.0*nExtruded/nTotFaces << "%)."
             << " Removed extrusion at " << nTotChanged << " faces."
+            << endl
+            << "Added " << nAddedCells << " out of " << nIdealAddedCells
+            << " cells (" << 100.0*nAddedCells/nIdealAddedCells << "%)."
             << endl;
 
         if (nTotChanged == 0)

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -192,7 +192,8 @@ class autoLayerDriver
                     const indirectPrimitivePatch& pp,
                     pointField& patchDisp,
                     labelList& patchNLayers,
-                    List<extrudeMode>& extrudeStatus
+                    List<extrudeMode>& extrudeStatus,
+                    label& nIdealAddedCells
                 ) const;
 
                 //- Helper function to make a pointVectorField with correct
@@ -469,6 +470,7 @@ class autoLayerDriver
                     const label nSmoothNormals,
                     const label nSmoothSurfaceNormals,
                     const scalar minMedianAxisAngleCos,
+                    const scalar featureAngle,
 
                     pointVectorField& dispVec,
                     pointScalarField& medialRatio,

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C

@@ -33,9 +33,10 @@ Description
 #include "motionSmoother.H"
 #include "pointData.H"
 #include "PointEdgeWave.H"
-#include "OFstream.H"
+#include "OBJstream.H"
 #include "meshTools.H"
 #include "PatchTools.H"
+#include "unitConversion.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -775,6 +776,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
     const label nSmoothNormals,
     const label nSmoothSurfaceNormals,
     const scalar minMedianAxisAngleCos,
+    const scalar featureAngle,
 
     pointVectorField& dispVec,
     pointScalarField& medialRatio,
@@ -904,7 +906,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
                 const edge& e = edges[edgeI];
                 // Approximate medial axis location on edge.
                 //const point medialAxisPt = e.centre(points);
-                vector eVec = e.vec(mesh.points());
+                vector eVec = e.vec(points);
                 scalar eMag = mag(eVec);
                 if (eMag > VSMALL)
                 {
@@ -961,42 +963,111 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
 
         labelHashSet adaptPatches(meshMover.adaptPatchIDs());
 
+
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
-            //const pointPatchVectorField& pvf =
+            const pointPatchVectorField& pvf =
-            //    meshMover.displacement().boundaryField()[patchI];
+                meshMover.displacement().boundaryField()[patchI];
 
             if
             (
                 !pp.coupled()
              && !isA<emptyPolyPatch>(pp)
-             //&&  pvf.constraintType() != word::null //exclude fixedValue
              && !adaptPatches.found(patchI)
             )
             {
-                Info<< "Inserting points on patch " << pp.name() << endl;
-
                 const labelList& meshPoints = pp.meshPoints();
 
-                forAll(meshPoints, i)
+                // Check the type of the patchField. The types are
+                //  - fixedValue (0 or more layers) but the >0 layers have
+                //    already been handled in the adaptPatches loop
+                //  - constraint (but not coupled) types, e.g. symmetryPlane,
+                //    slip.
+                if (pvf.fixesValue())
                 {
-                    label pointI = meshPoints[i];
+                    // Disable all movement on fixedValue patchFields
+                    Info<< "Inserting all points on patch " << pp.name()
+                        << endl;
 
-                    if (!pointMedialDist[pointI].valid(dummyTrackData))
+                    forAll(meshPoints, i)
                     {
-                        maxPoints.append(pointI);
+                        label pointI = meshPoints[i];
-                        maxInfo.append
+                        if (!pointMedialDist[pointI].valid(dummyTrackData))
-                        (
+                        {
-                            pointData
+                            maxPoints.append(pointI);
+                            maxInfo.append
                             (
-                                points[pointI],
+                                pointData
-                                0.0,
+                                (
-                                pointI,         // passive data
+                                    points[pointI],
-                                vector::zero    // passive data
+                                    0.0,
-                            )
+                                    pointI,         // passive data
-                        );
+                                    vector::zero    // passive data
-                        pointMedialDist[pointI] = maxInfo.last();
+                                )
+                            );
+                            pointMedialDist[pointI] = maxInfo.last();
+                        }
+                    }
+                }
+                else
+                {
+                    // Based on geometry: analyse angle w.r.t. nearest moving
+                    // point. In the pointWallDist we transported the
+                    // normal as the passive vector. Note that this points
+                    // out of the originating wall so inside of the domain
+                    // on this patch.
+                    Info<< "Inserting points on patch " << pp.name()
+                        << " if angle to nearest layer patch > "
+                        << featureAngle << " degrees." << endl;
+
+                    scalar featureAngleCos = Foam::cos(degToRad(featureAngle));
+                    pointField pointNormals
+                    (
+                        PatchTools::pointNormals
+                        (
+                            mesh,
+                            pp,
+                            identity(pp.size())+pp.start()
+                        )
+                    );
+
+                    forAll(meshPoints, i)
+                    {
+                        label pointI = meshPoints[i];
+
+                        if (!pointMedialDist[pointI].valid(dummyTrackData))
+                        {
+                            // Check if angle not too large.
+                            scalar cosAngle =
+                            (
+                               -pointWallDist[pointI].v()
+                              & pointNormals[i]
+                            );
+                            if (cosAngle > featureAngleCos)
+                            {
+                                // Extrusion direction practically perpendicular
+                                // to the patch. Disable movement at the patch.
+
+                                maxPoints.append(pointI);
+                                maxInfo.append
+                                (
+                                    pointData
+                                    (
+                                        points[pointI],
+                                        0.0,
+                                        pointI,         // passive data
+                                        vector::zero    // passive data
+                                    )
+                                );
+                                pointMedialDist[pointI] = maxInfo.last();
+                            }
+                            else
+                            {
+                                // Extrusion direction makes angle with patch
+                                // so allow sliding - don't insert zero points
+                            }
+                        }
                     }
                 }
             }
@@ -1141,13 +1212,12 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
     // reduce thickness where thickness/medial axis distance large
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    autoPtr<OFstream> str;
+    autoPtr<OBJstream> str;
-    label vertI = 0;
     if (debug)
     {
         str.reset
         (
-            new OFstream
+            new OBJstream
             (
                 mesh.time().path()
               / "thicknessRatioExcludePoints_"
@@ -1159,13 +1229,12 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
             << str().name() << endl;
     }
 
-    autoPtr<OFstream> medialVecStr;
+    autoPtr<OBJstream> medialVecStr;
-    label medialVertI = 0;
     if (debug)
     {
         medialVecStr.reset
         (
-            new OFstream
+            new OBJstream
             (
                 mesh.time().path()
               / "thicknessRatioExcludeMedialVec_"
@@ -1227,21 +1296,19 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
                 if (str.valid())
                 {
                     const point& pt = mesh.points()[pointI];
-                    meshTools::writeOBJ(str(), pt);
+                    str().write(linePointRef(pt, pt+patchDisp[patchPointI]));
-                    vertI++;
-                    meshTools::writeOBJ(str(), pt+patchDisp[patchPointI]);
-                    vertI++;
-                    str()<< "l " << vertI-1 << ' ' << vertI << nl;
                 }
                 if (medialVecStr.valid())
                 {
                     const point& pt = mesh.points()[pointI];
-                    meshTools::writeOBJ(medialVecStr(), pt);
+                    medialVecStr().write
-                    medialVertI++;
+                    (
-                    meshTools::writeOBJ(medialVecStr(), medialVec[pointI]);
+                        linePointRef
-                    medialVertI++;
+                        (
-                    medialVecStr()<< "l " << medialVertI-1
+                            pt,
-                        << ' ' << medialVertI << nl;
+                            medialVec[pointI]
+                        )
+                    );
                 }
             }
         }