BUG: snappyHexMesh: fix faceZone splitting of surfaces

This is a bit complicated. When allocating cells to zones
(meshRefinementBaffles::zonify) in finishes with findCellTopo which
tried to fix the odd cell which wasn't put into the correct region. This
was actually modifying whole regions to be a certain cellZone so if
there was some 'bleeding' it would re-assign a whole region to be e.g.
background and get deleted. Instead it now will only reassign single
cells if these are
- unassigned
- inbetween multiple, differing cellZones

Fixes the simpleFoam/rotorDisk meshing (gitlab #141)

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C

@@ -268,8 +268,6 @@ void Foam::meshRefinement::calcCellCellRays
 // Returns first intersection if there are more than one.
 void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
 {
-    const pointField& cellCentres = mesh_.cellCentres();
-
     // Stats on edges to test. Count proc faces only once.
     PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
 
@@ -309,28 +307,17 @@ void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
     // Collect segments we want to test for
     pointField start(changedFaces.size());
     pointField end(changedFaces.size());
-
-    forAll(changedFaces, i)
     {
-        label faceI = changedFaces[i];
+        labelList minLevel;
-        label own = mesh_.faceOwner()[faceI];
+        calcCellCellRays
-
+        (
-        start[i] = cellCentres[own];
+            neiCc,
-        if (mesh_.isInternalFace(faceI))
+            neiLevel,
-        {
+            changedFaces,
-            end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
+            start,
-        }
+            end,
-        else
+            minLevel
-        {
+        );
-            end[i] = neiCc[faceI-mesh_.nInternalFaces()];
-        }
-    }
-
-    // Extend segments a bit
-    {
-        const vectorField smallVec(ROOTSMALL*(end-start));
-        start -= smallVec;
-        end += smallVec;
     }
 
 
@@ -2964,8 +2951,6 @@ void Foam::meshRefinement::dumpRefinementLevel() const
 void Foam::meshRefinement::dumpIntersections(const fileName& prefix) const
 {
     {
-        const pointField& cellCentres = mesh_.cellCentres();
-
         OFstream str(prefix + "_edges.obj");
         label vertI = 0;
         Pout<< "meshRefinement::dumpIntersections :"
@@ -2986,27 +2971,17 @@ void Foam::meshRefinement::dumpIntersections(const fileName& prefix) const
         // Collect segments we want to test for
         pointField start(intersectionFaces.size());
         pointField end(intersectionFaces.size());
-
-        forAll(intersectionFaces, i)
         {
-            label faceI = intersectionFaces[i];
+            labelList minLevel;
-            start[i] = cellCentres[mesh_.faceOwner()[faceI]];
+            calcCellCellRays
-
+            (
-            if (mesh_.isInternalFace(faceI))
+                neiCc,
-            {
+                labelList(neiCc.size(), -1),
-                end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
+                intersectionFaces,
-            }
+                start,
-            else
+                end,
-            {
+                minLevel
-                end[i] = neiCc[faceI-mesh_.nInternalFaces()];
+            );
-            }
-        }
-
-        // Extend segments a bit
-        {
-            const vectorField smallVec(ROOTSMALL*(end-start));
-            start -= smallVec;
-            end += smallVec;
         }
 
 

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H

@@ -669,7 +669,6 @@ private:
             //- Determines cell zone from cell region information.
             bool calcRegionToZone
             (
-                const label backgroundZoneID,
                 const label surfZoneI,
                 const label ownRegion,
                 const label neiRegion,
@@ -681,7 +680,6 @@ private:
             //  marked in namedSurfaceIndex regarded as blocked.
             void findCellZoneTopo
             (
-                const label backgroundZoneID,
                 const pointField& locationsInMesh,
                 const labelList& allSurfaceIndex,
                 const labelList& namedSurfaceIndex,
@@ -702,7 +700,6 @@ private:
             void zonify
             (
                 const bool allowFreeStandingZoneFaces,
-                const label backgroundZoneID,
                 const pointField& locationsInMesh,
                 const wordList& zonesInMesh,
 

src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -311,7 +311,6 @@ void Foam::meshRefinement::getBafflePatches
         zonify
         (
             true,               // allowFreeStandingZoneFaces
-            -2,                 // zone to put unreached cells into
             locationsInMesh,
             zonesInMesh,
 
@@ -1811,75 +1810,8 @@ void Foam::meshRefinement::findCellZoneInsideWalk
 }
 
 
-bool Foam::meshRefinement::calcRegionToZone
-(
-    const label backgroundZoneID,
-    const label surfZoneI,
-    const label ownRegion,
-    const label neiRegion,
-
-    labelList& regionToCellZone
-) const
-{
-    bool changed = false;
-
-    // Check whether inbetween different regions
-    if (ownRegion != neiRegion)
-    {
-        // Jump. Change one of the sides to my type.
-
-        // 1. Interface between my type and unset region.
-        // Set region to keepRegion
-
-        if (regionToCellZone[ownRegion] == -2)
-        {
-            if (regionToCellZone[neiRegion] == surfZoneI)
-            {
-                // Face between unset and my region. Put unset
-                // region into keepRegion
-                //MEJ: see comment in findCellZoneTopo
-                if (backgroundZoneID != -2)
-                {
-                    regionToCellZone[ownRegion] = backgroundZoneID;
-                    changed = true;
-                }
-            }
-            else if (regionToCellZone[neiRegion] != -2)
-            {
-                // Face between unset and other region.
-                // Put unset region into my region
-                regionToCellZone[ownRegion] = surfZoneI;
-                changed = true;
-            }
-        }
-        else if (regionToCellZone[neiRegion] == -2)
-        {
-            if (regionToCellZone[ownRegion] == surfZoneI)
-            {
-                // Face between unset and my region. Put unset
-                // region into keepRegion
-                if (backgroundZoneID != -2)
-                {
-                    regionToCellZone[neiRegion] = backgroundZoneID;
-                    changed = true;
-                }
-            }
-            else if (regionToCellZone[ownRegion] != -2)
-            {
-                // Face between unset and other region.
-                // Put unset region into my region
-                regionToCellZone[neiRegion] = surfZoneI;
-                changed = true;
-            }
-        }
-    }
-    return changed;
-}
-
-
 void Foam::meshRefinement::findCellZoneTopo
 (
-    const label backgroundZoneID,
     const pointField& locationsInMesh,
     const labelList& allSurfaceIndex,
     const labelList& namedSurfaceIndex,
@@ -1893,17 +1825,12 @@ void Foam::meshRefinement::findCellZoneTopo
     // know which side of the face it relates to. So all we are doing here
     // is get the correspondence between surface/cellZone and regionSplit
     // region.
-    // See the logic in calcRegionToZone. The problem is what to do
+    // This can occasionally go wrong when surfaces pass through a
-    // with unreachable parts of the mesh (cellToZone = -2).
+    // cell centre and a cell completely gets blocked off so becomes a
-    // In 23x this routine was only called for actually zoneing an existing
+    // separate region. The problem is to distinguish between a cell which
-    // mesh so we had to do something with these cells and they
+    // properly should be deleted (so in the 'background' mesh)
-    // would get set to -1 (background). However in this version this routine
+    // or just a single-cell region from incorrect baffling. Currently the
-    // also gets used much earlier on when after the surface refinement it
+    // logic is just to detect cells that are on different regions.
-    // removes unreachable cells ('Removing mesh beyond surface intersections')
-    // and this is when we keep -2 so it gets removed.
-    // So the zone to set these unmarked cells to is provided as argument:
-    // - backgroundZoneID = -2 : do not change so remove cells
-    // - backgroundZoneID = -1 : put into background
 
 
     // Assumes:
@@ -1993,95 +1920,107 @@ void Foam::meshRefinement::findCellZoneTopo
         }
     }
 
+    // Synchronise any changes due to locationInMesh 
+    Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
+    Pstream::listCombineScatter(regionToCellZone);
 
-    // Find correspondence between cell zone and surface
+    // Get coupled neighbour cellRegion
-    // by changing cell zone every time we cross a surface.
+    labelList neiCellRegion;
-    while (true)
+    syncTools::swapBoundaryCellList(mesh_, cellRegion, neiCellRegion);
+
+
+
+
+    // Detect unassigned cell (region) bordering two cellZones. Put these
+    // cells into either neighbouring cellZone.
+    labelList growCellToZone(mesh_.nCells(), -2);
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
     {
-        // Synchronise regionToCellZone.
+        label own = mesh_.faceOwner()[faceI];
-        // Note:
+        label ownRegion = cellRegion[own];
-        // - region numbers are identical on all processors
+        label nei = mesh_.faceNeighbour()[faceI];
-        // - keepRegion is identical ,,
+        label neiRegion = cellRegion[nei];
-        // - cellZones are identical ,,
-        // This done at top of loop to account for geometric matching
-        // not being synchronised.
-        Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
-        Pstream::listCombineScatter(regionToCellZone);
 
-
+        // Check whether inbetween different regions
-        bool changed = false;
+        if (ownRegion != neiRegion)
-
-        // Internal faces
-
-        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
         {
-            label surfI = namedSurfaceIndex[faceI];
+            if (regionToCellZone[ownRegion] == -2)
-
-            // Connected even if no cellZone defined for surface
-            if (surfI != -1)
             {
-                // Calculate region to zone from cellRegions on either side
+                if (regionToCellZone[neiRegion] != -2)
-                // of internal face.
+                {
-                bool changedCell = calcRegionToZone
+                    if (growCellToZone[own] == -2)
-                (
+                    {
-                    backgroundZoneID,
+                        // First visit of cell
-                    surfaceToCellZone[surfI],
+                        growCellToZone[own] = regionToCellZone[neiRegion];
-                    cellRegion[mesh_.faceOwner()[faceI]],
+                    }
-                    cellRegion[mesh_.faceNeighbour()[faceI]],
+                    else if (regionToCellZone[neiRegion] != growCellToZone[own])
-                    regionToCellZone
+                    {
-                );
+                        // Found a cell bordering multiple cellZones. Grow.
-
+                        cellToZone[own] = growCellToZone[own];
-                changed = changed | changedCell;
+                    }
+                }
+            }
+            else if (regionToCellZone[neiRegion] == -2)
+            {
+                if (growCellToZone[nei] == -2)
+                {
+                    // First visit of cell
+                    growCellToZone[nei] = regionToCellZone[ownRegion];
+                }
+                else if (regionToCellZone[ownRegion] != growCellToZone[nei])
+                {
+                    cellToZone[nei] = growCellToZone[nei];
+                }
             }
         }
+    }
 
-        // Boundary faces
 
-        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-        // Get coupled neighbour cellRegion
+    forAll(patches, patchI)
-        labelList neiCellRegion;
+    {
-        syncTools::swapBoundaryCellList(mesh_, cellRegion, neiCellRegion);
+        const polyPatch& pp = patches[patchI];
 
-        // Calculate region to zone from cellRegions on either side of coupled
+        if (pp.coupled())
-        // face.
-        forAll(patches, patchI)
         {
-            const polyPatch& pp = patches[patchI];
+            forAll(pp, i)
-
-            if (pp.coupled())
             {
-                forAll(pp, i)
+                label faceI = pp.start()+i;
+
+                label own = mesh_.faceOwner()[faceI];
+                label ownRegion = cellRegion[own];
+                label neiRegion = neiCellRegion[faceI-mesh_.nInternalFaces()];
+
+                // Check whether inbetween different regions
+                if (ownRegion != neiRegion)
                 {
-                    label faceI = pp.start()+i;
+                    if (regionToCellZone[ownRegion] == -2)
-
-                    label surfI = namedSurfaceIndex[faceI];
-
-                    // Connected even if no cellZone defined for surface
-                    if (surfI != -1)
                     {
-                        bool changedCell = calcRegionToZone
+                        if (regionToCellZone[neiRegion] != -2)
-                        (
+                        {
-                            backgroundZoneID,
+                            if (growCellToZone[own] == -2)
-                            surfaceToCellZone[surfI],
+                            {
-                            cellRegion[mesh_.faceOwner()[faceI]],
+                                // First visit of cell
-                            neiCellRegion[faceI-mesh_.nInternalFaces()],
+                                growCellToZone[own] =
-                            regionToCellZone
+                                    regionToCellZone[neiRegion];
-                        );
+                            }
-
+                            else if
-                        changed = changed | changedCell;
+                            (
+                                regionToCellZone[neiRegion]
+                             != growCellToZone[own]
+                            )
+                            {
+                                // Found a cell bordering multiple cellZones.
+                                cellToZone[own] = growCellToZone[own];
+                            }
+                        }
                     }
                 }
             }
         }
-
-        if (!returnReduce(changed, orOp<bool>()))
-        {
-            break;
-        }
     }
 
-
     if (debug)
     {
         forAll(regionToCellZone, regionI)
@@ -2200,10 +2139,6 @@ void Foam::meshRefinement::getIntersections
     PackedBoolList& posOrientation
 ) const
 {
-    const pointField& cellCentres = mesh_.cellCentres();
-    const labelList& faceOwner = mesh_.faceOwner();
-    const labelList& faceNeighbour = mesh_.faceNeighbour();
-
     namedSurfaceIndex.setSize(mesh_.nFaces());
     namedSurfaceIndex = -1;
 
@@ -2219,30 +2154,19 @@ void Foam::meshRefinement::getIntersections
     pointField start(testFaces.size());
     pointField end(testFaces.size());
 
-    forAll(testFaces, i)
     {
-        label faceI = testFaces[i];
+        labelList minLevel;
-
+        calcCellCellRays
-        if (mesh_.isInternalFace(faceI))
+        (
-        {
+            neiCc,
-            start[i] = cellCentres[faceOwner[faceI]];
+            labelList(neiCc.size(), -1),
-            end[i] = cellCentres[faceNeighbour[faceI]];
+            testFaces,
-        }
+            start,
-        else
+            end,
-        {
+            minLevel
-            start[i] = cellCentres[faceOwner[faceI]];
+        );
-            end[i] = neiCc[faceI-mesh_.nInternalFaces()];
-        }
     }
 
-    // Extend segments a bit
-    {
-        const vectorField smallVec(ROOTSMALL*(end-start));
-        start -= smallVec;
-        end += smallVec;
-    }
-
-
     // Do test for intersections
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
     // Note that we intersect all intersected faces again. Could reuse
@@ -2340,7 +2264,6 @@ void Foam::meshRefinement::getIntersections
 void Foam::meshRefinement::zonify
 (
     const bool allowFreeStandingZoneFaces,
-    const label backgroundZoneID,
     const pointField& locationsInMesh,
     const wordList& zonesInMesh,
 
@@ -2535,24 +2458,53 @@ void Foam::meshRefinement::zonify
         );
     }
 
+    if (debug&MESH)
+    {
+        Pout<< "Writing cell zone allocation on mesh to time "
+            << timeName() << endl;
+        volScalarField volCellToZone
+        (
+            IOobject
+            (
+                "cellToZoneBeforeTopo",
+                timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimless, 0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+
+        forAll(cellToZone, cellI)
+        {
+            volCellToZone[cellI] = cellToZone[cellI];
+        }
+        volCellToZone.write();
+    }
+
+
+    // 5. Find any unassigned regions (from regionSplit). This will
+    //    just walk slightly out to cover those single cells that
+    //    are baffled off into a separate region.
+
+    Info<< "Growing from known cellZones to fix cells inbetween cellZones"
+        << nl << endl;
+
+    findCellZoneTopo
+    (
+        pointField(0),
+        globalRegion1,      // To split up cells
+        namedSurfaceIndex,  // Step across named surfaces to propagate
+        surfaceToCellZone,
+        cellToZone
+    );
 
-    // 5. Find any unassigned regions (from regionSplit)
 
     if (namedSurfaces.size())
     {
-        Info<< "Walking from known cellZones; crossing a faceZone "
-            << "face changes cellZone" << nl << endl;
-
-        findCellZoneTopo
-        (
-            backgroundZoneID,
-            pointField(0),
-            globalRegion1,      // To split up cells
-            namedSurfaceIndex,  // Step across named surfaces to propagate
-            surfaceToCellZone,
-            cellToZone
-        );
-
         // Make sure namedSurfaceIndex is unset inbetween same cell zones.
         if (!allowFreeStandingZoneFaces)
         {
@@ -3621,6 +3573,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
     if (debug&MESH)
     {
+        runTime++;
         Pout<< "Writing subsetted mesh to time " << timeName()
             << endl;
         write
@@ -4316,7 +4269,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     zonify
     (
         allowFreeStandingZoneFaces,
-        -1,                 // Set all cells with cellToZone -2 to -1
         locationsInMesh,
         zonesInMesh,
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C

@@ -1316,6 +1316,8 @@ void Foam::snappyRefineDriver::removeInsideCells
 
     if (debug&meshRefinement::MESH)
     {
+        const_cast<Time&>(mesh.time())++;
+
         Pout<< "Writing subsetted mesh to time "
             << meshRefiner_.timeName() << endl;
         meshRefiner_.write