ENH: snappyHexMesh: Remove cells inside any limitRegion with level -1. Fixes #852.

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -314,8 +314,9 @@ castellatedMeshControls
     // refinement (from features, refinementSurfaces, refinementRegions)
     // in a given geometric region. The syntax is exactly the same as for the
     // refinementRegions; the cell level now specifies the upper limit
-    // for any cell. Note that it does not override the refinement constraints
-    // given by the nCellsBetweenLevels setting.
+    // for any cell. (a special setting is cell level -1 which will remove
+    // any cells inside the region). Note that it does not override the
+    // refinement constraints given by the nCellsBetweenLevels setting.
     limitRegions
     {
     }

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H

@@ -522,6 +522,16 @@ private:
                 labelList& neiPatch
             ) const;
 
+            autoPtr<mapPolyMesh> splitMesh
+            (
+                const label nBufferLayers,
+                const labelList& globalToMasterPatch,
+                const labelList& globalToSlavePatch,
+                labelList& cellRegion,
+                labelList& ownPatch,
+                labelList& neiPatch
+            );
+
             //- Repatches external face or creates baffle for internal face
             //  with user specified patches (might be different for both sides).
             //  Returns label of added face.
@@ -861,6 +871,12 @@ public:
                 return shells_;
             }
 
+            //- Reference to limit shells (regions)
+            const shellSurfaces& limitShells() const
+            {
+                return limitShells_;
+            }
+
             //- Reference to meshcutting engine
             const hexRef8& meshCutter() const
             {
@@ -1108,6 +1124,17 @@ public:
                 const pointField& locationsOutsideMesh
             );
 
+            //- Remove cells from limitRegions if level -1
+            autoPtr<mapPolyMesh> removeLimitShells
+            (
+                const label nBufferLayers,
+                const label nErodeCellZones,
+                const labelList& globalToMasterPatch,
+                const labelList& globalToSlavePatch,
+                const pointField& locationsInMesh,
+                const wordList& regionsInMesh
+            );
+
             //- Find boundary points that connect to more than one cell
             //  region and split them.
             autoPtr<mapPolyMesh> dupNonManifoldPoints(const localPointRegion&);

src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,7 +48,7 @@ License
 #include "polyMeshAdder.H"
 #include "IOmanip.H"
 #include "refinementParameters.H"
-
+#include "shellSurfaces.H"
 #include "zeroGradientFvPatchFields.H"
 #include "volFields.H"
 
@@ -3874,7 +3874,28 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
         cellRegion
     );
 
+    return splitMesh
+    (
+        nBufferLayers,
+        globalToMasterPatch,
+        globalToSlavePatch,
+        cellRegion,
+        ownPatch,
+        neiPatch
+    );
+}
+
+Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
+(
+    const label nBufferLayers,
+    const labelList& globalToMasterPatch,
+    const labelList& globalToSlavePatch,
 
+    labelList& cellRegion,
+    labelList& ownPatch,
+    labelList& neiPatch
+)
+{
     // Walk out nBufferlayers from region boundary
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // (modifies cellRegion, ownPatch)
@@ -4049,7 +4070,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
     label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
     reduce(nCellsToKeep, sumOp<label>());
 
-    Info<< "Keeping all cells containing points " << locationsInMesh << endl
+    Info<< "Keeping all cells containing inside points" << endl
         << "Selected for keeping : " << nCellsToKeep << " cells." << endl;
 
 
@@ -4090,6 +4111,90 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
 }
 
 
+Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeLimitShells
+(
+    const label nBufferLayers,
+    const label nErodeCellZones,
+    const labelList& globalToMasterPatch,
+    const labelList& globalToSlavePatch,
+    const pointField& locationsInMesh,
+    const wordList& zonesInMesh
+)
+{
+    // Determine patches to put intersections into
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Swap neighbouring cell centres and cell level
+    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+    calcNeighbourData(neiLevel, neiCc);
+
+    // Find intersections with all unnamed(!) surfaces
+    labelList ownPatch, neiPatch;
+    getBafflePatches
+    (
+        nErodeCellZones,
+        globalToMasterPatch,
+
+        locationsInMesh,
+        zonesInMesh,
+
+        neiLevel,
+        neiCc,
+
+        ownPatch,
+        neiPatch
+    );
+
+
+    labelList cellRegion(mesh_.nCells(), 0);
+    // Find any cells inside a limit shell with minLevel -1
+    labelList levelShell;
+    limitShells_.findLevel
+    (
+        mesh_.cellCentres(),
+        labelList(mesh_.nCells(), -1),  // pick up only shells with -1
+        levelShell
+    );
+    forAll(levelShell, celli)
+    {
+        if (levelShell[celli] != -1)
+        {
+            // Mark cell region so it gets deletec
+            cellRegion[celli] = -1;
+        }
+    }
+
+    autoPtr<mapPolyMesh> mapPtr = splitMesh
+    (
+        nBufferLayers,
+        globalToMasterPatch,
+        globalToSlavePatch,
+        cellRegion,
+        ownPatch,
+        neiPatch
+    );
+
+    if (debug&meshRefinement::MESH)
+    {
+        const_cast<Time&>(mesh_.time())++;
+        Pout<< "Writing mesh after removing limitShells "
+            << " to time " << timeName() << endl;
+        write
+        (
+            debugType(debug),
+            writeType
+            (
+                writeLevel()
+              | WRITEMESH
+            ),
+            mesh_.time().path()/timeName()
+        );
+    }
+    return mapPtr;
+}
+
+
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints
 (
     const localPointRegion& regionSide

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C

@@ -1348,6 +1348,12 @@ void Foam::snappyRefineDriver::removeInsideCells
     const label nBufferLayers
 )
 {
+    // Skip if no limitRegion and zero bufferLayers
+    if (meshRefiner_.limitShells().shells().size() == 0 && nBufferLayers == 0)
+    {
+        return;
+    }
+
     Info<< nl
         << "Removing mesh beyond surface intersections" << nl
         << "------------------------------------------" << nl
@@ -1360,6 +1366,19 @@ void Foam::snappyRefineDriver::removeInsideCells
        const_cast<Time&>(mesh.time())++;
     }
 
+    // Remove any cells inside limitShells with level -1
+    meshRefiner_.removeLimitShells
+    (
+        nBufferLayers,
+        1,
+        globalToMasterPatch_,
+        globalToSlavePatch_,
+        refineParams.locationsInMesh(),
+        refineParams.zonesInMesh()
+    );
+
+    // Fix any additional (e.g. locationsOutsideMesh). Note: probably not
+    // nessecary.
     meshRefiner_.splitMesh
     (
         nBufferLayers,                  // nBufferLayers
@@ -2256,6 +2275,14 @@ void Foam::snappyRefineDriver::doRefine
         100    // maxIter
     );
 
+    // Remove any extra cells from limitRegion with level -1, without
+    // adding any buffer layer this time
+    removeInsideCells
+    (
+        refineParams,
+        0       // nBufferLayers
+    );
+
     // Refine any hexes with 5 or 6 faces refined to make smooth edges
     danglingCellRefine
     (