ENH: snappyHexMesh: add different method for cellZone-faceZone consistency. Fixes #629.

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -378,6 +378,14 @@ castellatedMeshControls
     // Optional: do not refine surface cells with opposite faces of
     //           differing refinement levels
     //interfaceRefine false;
+
+    // Optional: use an erosion instead of region assignment to allocate
+    //           left-over cells to the background region (i.e. make cellZones
+    //           consistent with the intersections of the surface).
+    //           Erosion is specified as a number of erosion iterations.
+    //           Erosion has less chance of bleeding and changing the zone
+    //           for a complete region.
+    //nCellZoneErodeIter 2;
 }
 
 // Settings for the snapping.

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H

@@ -511,6 +511,7 @@ private:
             //- Determine patches for baffles
             void getBafflePatches
             (
+                const label nErodeCellZones,
                 const labelList& globalToMasterPatch,
                 const pointField& locationsInMesh,
                 const wordList& regionsInMesh,
@@ -673,6 +674,17 @@ private:
                 labelList& cellToZone
             ) const;
 
+            //- Opposite of findCellTopo: finds assigned cell connected to
+            //  an unassigned one and puts it in the background zone.
+            void erodeCellZone
+            (
+                const label nErodeCellZones,
+                const label backgroundZoneID,
+                const labelList& unnamedSurfaceRegion,
+                const labelList& namedSurfaceIndex,
+                labelList& cellToZone
+            ) const;
+
             //- Make namedSurfaceIndex consistent with cellToZone
             //  - clear out any blocked faces inbetween same cell zone.
             void makeConsistentFaceIndex
@@ -686,6 +698,7 @@ private:
             void zonify
             (
                 const bool allowFreeStandingZoneFaces,
+                const label nErodeCellZones,
                 const label backgroundZoneID,
                 const pointField& locationsInMesh,
                 const wordList& zonesInMesh,
@@ -1037,6 +1050,7 @@ public:
                 const bool useTopologicalSnapDetection,
                 const bool removeEdgeConnectedCells,
                 const scalarField& perpendicularAngle,
+                const label nErodeCellZones,
                 const dictionary& motionDict,
                 Time& runTime,
                 const labelList& globalToMasterPatch,
@@ -1067,6 +1081,7 @@ public:
             autoPtr<mapPolyMesh> splitMesh
             (
                 const label nBufferLayers,
+                const label nErodeCellZones,
                 const labelList& globalToMasterPatch,
                 const labelList& globalToSlavePatch,
 
@@ -1148,6 +1163,7 @@ public:
             autoPtr<mapPolyMesh> zonify
             (
                 const bool allowFreeStandingZoneFaces,
+                const label nErodeCellZones,
                 const pointField& locationsInMesh,
                 const wordList& regionsInMesh,
                 wordPairHashTable& zonesToFaceZone

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-2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -279,6 +279,7 @@ void Foam::meshRefinement::getIntersections
 
 void Foam::meshRefinement::getBafflePatches
 (
+    const label nErodeCellZones,
     const labelList& globalToMasterPatch,
     const pointField& locationsInMesh,
     const wordList& zonesInMesh,
@@ -309,6 +310,7 @@ void Foam::meshRefinement::getBafflePatches
         zonify
         (
             true,               // allowFreeStandingZoneFaces
+            nErodeCellZones,
             -2,                 // zone to put unreached cells into
             locationsInMesh,
             zonesInMesh,
@@ -331,16 +333,13 @@ void Foam::meshRefinement::getBafflePatches
     labelList neiCellZone;
     syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
 
-    const labelList testFaces(intersectedFaces());
-
     ownPatch.setSize(mesh_.nFaces());
     ownPatch = -1;
     neiPatch.setSize(mesh_.nFaces());
     neiPatch = -1;
-    forAll(testFaces, i)
-    {
-        label faceI = testFaces[i];
 
+    forAll(ownPatch, faceI)
+    {
         if (unnamedRegion1[faceI] != -1 || unnamedRegion2[faceI] != -1)
         {
             label ownMasterPatch = -1;
@@ -2060,6 +2059,111 @@ void Foam::meshRefinement::findCellZoneTopo
 }
 
 
+void Foam::meshRefinement::erodeCellZone
+(
+    const label nErodeCellZones,
+    const label backgroundZoneID,
+    const labelList& unnamedSurfaceRegion,
+    const labelList& namedSurfaceIndex,
+    labelList& cellToZone
+) const
+{
+    // This routine fixes small problems with left over unassigned regions
+    // (after all off the unreachable bits of the mesh have been removed).
+    // The problem is that the cell zone information might be inconsistent
+    // with the face zone information. So what we do here is to erode
+    // any cell zones until we hit a named face.
+    // - backgroundZoneID = -2 : do not change so remove cells
+    // - backgroundZoneID = -1 : put into background
+    // Note that is the opposite of findCellZoneTopo which moves unassigned
+    // regions into a neighbouring region(=cellZone) unless there is an
+    // intersected faces inbetween the two.
+
+    for (label iter = 0; iter < nErodeCellZones; iter++)
+    {
+        label nChanged = 0;
+
+        labelList erodedCellToZone(cellToZone);
+
+        // Do internal faces
+        for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
+        {
+            if
+            (
+                unnamedSurfaceRegion[facei] == -1
+             && namedSurfaceIndex[facei] == -1
+            )
+            {
+                label own = mesh_.faceOwner()[facei];
+                label nei = mesh_.faceNeighbour()[facei];
+                if (cellToZone[own] == -2 && cellToZone[nei] >= -1)
+                {
+                    erodedCellToZone[nei] = backgroundZoneID;
+                    nChanged++;
+                }
+                else if (cellToZone[nei] == -2 && cellToZone[own] >= -1)
+                {
+                    erodedCellToZone[own] = backgroundZoneID;
+                    nChanged++;
+                }
+            }
+        }
+
+        // Do boundary faces
+
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+        // Get coupled neighbour cellRegion
+        labelList neiCellZone;
+        syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
+
+        // Calculate region to zone from cellRegions on either side of coupled
+        // face.
+        forAll(patches, patchi)
+        {
+            const polyPatch& pp = patches[patchi];
+
+            if (pp.coupled())
+            {
+                forAll(pp, i)
+                {
+                    label facei = pp.start()+i;
+                    if
+                    (
+                        unnamedSurfaceRegion[facei] == -1
+                     && namedSurfaceIndex[facei] == -1
+                    )
+                    {
+                        label own = mesh_.faceOwner()[facei];
+                        label bFacei = facei-mesh_.nInternalFaces();
+                        if (neiCellZone[bFacei] == -2 && cellToZone[own] >= -1)
+                        {
+                            erodedCellToZone[own] = backgroundZoneID;
+                            nChanged++;
+                        }
+                    }
+                }
+            }
+        }
+
+        cellToZone.transfer(erodedCellToZone);
+
+        reduce(nChanged, sumOp<label>());
+        if (debug)
+        {
+            Pout<< "erodeCellZone : eroded " << nChanged
+                << " cells (moved from cellZone to background zone "
+                << backgroundZoneID << endl;
+        }
+
+        if (nChanged == 0)
+        {
+            break;
+        }
+    }
+}
+
+
 void Foam::meshRefinement::makeConsistentFaceIndex
 (
     const labelList& surfaceMap,
@@ -2281,6 +2385,7 @@ void Foam::meshRefinement::getIntersections
 void Foam::meshRefinement::zonify
 (
     const bool allowFreeStandingZoneFaces,
+    const label nErodeCellZones,
     const label backgroundZoneID,
     const pointField& locationsInMesh,
     const wordList& zonesInMesh,
@@ -2496,18 +2601,39 @@ void Foam::meshRefinement::zonify
 
     if (namedSurfaces.size())
     {
-        Info<< "Walking from known cellZones; crossing a faceZone "
-            << "face changes cellZone" << nl << endl;
+        if (nErodeCellZones <= 0)
+        {
+            Info<< "Walking from known cellZones; crossing a faceZone "
+                << "face changes cellZone" << nl << endl;
+
+            // Put unassigned regions into any connected cellZone
+            findCellZoneTopo
+            (
+                backgroundZoneID,
+                pointField(0),
+                unnamedRegion1,      // Intersections with unnamed surfaces
+                namedSurfaceIndex,   // Intersections with named surfaces
+                surfaceToCellZone,
+                cellToZone
+            );
+        }
+        else
+        {
+            Info<< "Eroding cellZone cells to make these consistent with"
+                << " faceZone faces" << nl << endl;
+
+            // Erode cell zone cells (connected to an unassigned cell)
+            // and put them into backgroundZone
+            erodeCellZone
+            (
+                nErodeCellZones,
+                backgroundZoneID,
+                unnamedRegion1,
+                namedSurfaceIndex,
+                cellToZone
+            );
+        }
 
-        findCellZoneTopo
-        (
-            backgroundZoneID,
-            pointField(0),
-            unnamedRegion1,      // Intersections with unnamed surfaces
-            namedSurfaceIndex,   // Intersections with named surfaces
-            surfaceToCellZone,
-            cellToZone
-        );
 
         // Make sure namedSurfaceIndex is unset inbetween same cell zones.
         if (!allowFreeStandingZoneFaces)
@@ -3424,6 +3550,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
     const bool useTopologicalSnapDetection,
     const bool removeEdgeConnectedCells,
     const scalarField& perpendicularAngle,
+    const label nErodeCellZones,
     const dictionary& motionDict,
     Time& runTime,
     const labelList& globalToMasterPatch,
@@ -3452,6 +3579,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
     labelList ownPatch, neiPatch;
     getBafflePatches
     (
+        nErodeCellZones,
         globalToMasterPatch,
 
         locationsInMesh,
@@ -3522,6 +3650,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
             labelList ownPatch, neiPatch;
             getBafflePatches
             (
+                nErodeCellZones,
                 globalToMasterPatch,
 
                 locationsInMesh,
@@ -3689,6 +3818,7 @@ void Foam::meshRefinement::mergeFreeStandingBaffles
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
 (
     const label nBufferLayers,
+    const label nErodeCellZones,
     const labelList& globalToMasterPatch,
     const labelList& globalToSlavePatch,
 
@@ -3709,6 +3839,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
     labelList ownPatch, neiPatch;
     getBafflePatches
     (
+        nErodeCellZones,
         globalToMasterPatch,
 
         locationsInMesh,
@@ -4202,6 +4333,7 @@ Foam::meshRefinement::dupNonManifoldBoundaryPoints()
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 (
     const bool allowFreeStandingZoneFaces,
+    const label nErodeCellZones,
     const pointField& locationsInMesh,
     const wordList& zonesInMesh,
     wordPairHashTable& zonesToFaceZone
@@ -4278,6 +4410,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         zonify
         (
             allowFreeStandingZoneFaces,
+            nErodeCellZones,// Use erosion (>0) or regionSplit to clean up
             -1,             // Set all cells with cellToZone -2 to -1
             locationsInMesh,
             zonesInMesh,

src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C

@@ -68,7 +68,8 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
     interfaceRefine_
     (
         dict.lookupOrDefault<Switch>("interfaceRefine", true)
-    )
+    ),
+    nErodeCellZone_(dict.lookupOrDefault<label>("nCellZoneErodeIter", 0))
 {
     point locationInMesh;
     if (dict.readIfPresent("locationInMesh", locationInMesh))

src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.H

@@ -106,6 +106,9 @@ class refinementParameters
 
         Switch interfaceRefine_;
 
+        label nErodeCellZone_;
+
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -212,6 +215,11 @@ public:
                 return interfaceRefine_;
             }
 
+            label nErodeCellZone() const
+            {
+                return nErodeCellZone_;
+            }
+
 
         // Other
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C

@@ -1362,6 +1362,7 @@ void Foam::snappyRefineDriver::removeInsideCells
     meshRefiner_.splitMesh
     (
         nBufferLayers,                  // nBufferLayers
+        refineParams.nErodeCellZone(),
         globalToMasterPatch_,
         globalToSlavePatch_,
         refineParams.locationsInMesh(),
@@ -1605,6 +1606,7 @@ void Foam::snappyRefineDriver::baffleAndSplitMesh
         refineParams.useTopologicalSnapDetection(),
         false,                          // perpendicular edge connected cells
         scalarField(0),                 // per region perpendicular angle
+        refineParams.nErodeCellZone(),
 
         motionDict,
         const_cast<Time&>(mesh.time()),
@@ -1672,6 +1674,7 @@ void Foam::snappyRefineDriver::zonify
         meshRefiner_.zonify
         (
             refineParams.allowFreeStandingZoneFaces(),
+            refineParams.nErodeCellZone(),
             refineParams.locationsInMesh(),
             refineParams.zonesInMesh(),
             zonesToFaceZone
@@ -1731,6 +1734,7 @@ void Foam::snappyRefineDriver::splitAndMergeBaffles
         refineParams.useTopologicalSnapDetection(),
         handleSnapProblems,                 // remove perp edge connected cells
         perpAngle,                          // perp angle
+        refineParams.nErodeCellZone(),
 
         motionDict,
         const_cast<Time&>(mesh.time()),