From 5583b78d9d57c5f83bc8a201bf4837ed1d421fb2 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 20 Jun 2016 16:39:34 +0100
Subject: [PATCH] ENH: snappyHexMesh: revert zoneing back to
 e8d73e554686fe8e1efdb3de18dbe3e4bb9feaab

since growing single cells just makes the zoneing inconsistent with the
surface intersections (so you don't have a patch when converting surface
intersections to baffles)
---
 .../meshRefinement/meshRefinement.H           |   3 +
 .../meshRefinement/meshRefinementBaffles.C    | 265 +++++++++++-------
 2 files changed, 166 insertions(+), 102 deletions(-)

diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
index 50523060f9..c0b2724da1 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
@@ -669,6 +669,7 @@ private:
             //- Determines cell zone from cell region information.
             bool calcRegionToZone
             (
+                const label backgroundZoneID,
                 const label surfZoneI,
                 const label ownRegion,
                 const label neiRegion,
@@ -680,6 +681,7 @@ private:
             //  marked in namedSurfaceIndex regarded as blocked.
             void findCellZoneTopo
             (
+                const label backgroundZoneID,
                 const pointField& locationsInMesh,
                 const labelList& allSurfaceIndex,
                 const labelList& namedSurfaceIndex,
@@ -700,6 +702,7 @@ private:
             void zonify
             (
                 const bool allowFreeStandingZoneFaces,
+                const label backgroundZoneID,
                 const pointField& locationsInMesh,
                 const wordList& zonesInMesh,
 
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
index 8f6986d682..ebc759119f 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -311,6 +311,7 @@ void Foam::meshRefinement::getBafflePatches
         zonify
         (
             true,               // allowFreeStandingZoneFaces
+            -2,                 // zone to put unreached cells into
             locationsInMesh,
             zonesInMesh,
 
@@ -1810,8 +1811,75 @@ 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,
@@ -1825,12 +1893,17 @@ 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.
-    // This can occasionally go wrong when surfaces pass through a
-    // cell centre and a cell completely gets blocked off so becomes a
-    // separate region. The problem is to distinguish between a cell which
-    // properly should be deleted (so in the 'background' mesh)
-    // or just a single-cell region from incorrect baffling. Currently the
-    // logic is just to detect cells that are on different regions.
+    // See the logic in calcRegionToZone. The problem is what to do
+    // with unreachable parts of the mesh (cellToZone = -2).
+    // In 23x this routine was only called for actually zoneing an existing
+    // mesh so we had to do something with these cells and they
+    // would get set to -1 (background). However in this version this routine
+    // also gets used much earlier on when after the surface refinement it
+    // 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:
@@ -1920,107 +1993,95 @@ void Foam::meshRefinement::findCellZoneTopo
         }
     }
 
-    // Synchronise any changes due to locationInMesh 
-    Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
-    Pstream::listCombineScatter(regionToCellZone);
 
-    // Get coupled neighbour cellRegion
-    labelList neiCellRegion;
-    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++)
+    // Find correspondence between cell zone and surface
+    // by changing cell zone every time we cross a surface.
+    while (true)
     {
-        label own = mesh_.faceOwner()[faceI];
-        label ownRegion = cellRegion[own];
-        label nei = mesh_.faceNeighbour()[faceI];
-        label neiRegion = cellRegion[nei];
+        // Synchronise regionToCellZone.
+        // Note:
+        // - region numbers are identical on all processors
+        // - keepRegion is identical ,,
+        // - 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
-        if (ownRegion != neiRegion)
+
+        bool changed = false;
+
+        // Internal faces
+
+        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
         {
-            if (regionToCellZone[ownRegion] == -2)
+            label surfI = namedSurfaceIndex[faceI];
+
+            // Connected even if no cellZone defined for surface
+            if (surfI != -1)
             {
-                if (regionToCellZone[neiRegion] != -2)
-                {
-                    if (growCellToZone[own] == -2)
-                    {
-                        // First visit of cell
-                        growCellToZone[own] = regionToCellZone[neiRegion];
-                    }
-                    else if (regionToCellZone[neiRegion] != growCellToZone[own])
-                    {
-                        // Found a cell bordering multiple cellZones. Grow.
-                        cellToZone[own] = growCellToZone[own];
-                    }
-                }
-            }
-            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];
-                }
+                // Calculate region to zone from cellRegions on either side
+                // of internal face.
+                bool changedCell = calcRegionToZone
+                (
+                    backgroundZoneID,
+                    surfaceToCellZone[surfI],
+                    cellRegion[mesh_.faceOwner()[faceI]],
+                    cellRegion[mesh_.faceNeighbour()[faceI]],
+                    regionToCellZone
+                );
+
+                changed = changed | changedCell;
             }
         }
-    }
 
+        // Boundary faces
 
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
+        // Get coupled neighbour cellRegion
+        labelList neiCellRegion;
+        syncTools::swapBoundaryCellList(mesh_, cellRegion, neiCellRegion);
 
-        if (pp.coupled())
+        // Calculate region to zone from cellRegions on either side of coupled
+        // face.
+        forAll(patches, patchI)
         {
-            forAll(pp, i)
+            const polyPatch& pp = patches[patchI];
+
+            if (pp.coupled())
             {
-                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)
+                forAll(pp, i)
                 {
-                    if (regionToCellZone[ownRegion] == -2)
+                    label faceI = pp.start()+i;
+
+                    label surfI = namedSurfaceIndex[faceI];
+
+                    // Connected even if no cellZone defined for surface
+                    if (surfI != -1)
                     {
-                        if (regionToCellZone[neiRegion] != -2)
-                        {
-                            if (growCellToZone[own] == -2)
-                            {
-                                // First visit of cell
-                                growCellToZone[own] =
-                                    regionToCellZone[neiRegion];
-                            }
-                            else if
-                            (
-                                regionToCellZone[neiRegion]
-                             != growCellToZone[own]
-                            )
-                            {
-                                // Found a cell bordering multiple cellZones.
-                                cellToZone[own] = growCellToZone[own];
-                            }
-                        }
+                        bool changedCell = calcRegionToZone
+                        (
+                            backgroundZoneID,
+                            surfaceToCellZone[surfI],
+                            cellRegion[mesh_.faceOwner()[faceI]],
+                            neiCellRegion[faceI-mesh_.nInternalFaces()],
+                            regionToCellZone
+                        );
+
+                        changed = changed | changedCell;
                     }
                 }
             }
         }
+
+        if (!returnReduce(changed, orOp<bool>()))
+        {
+            break;
+        }
     }
 
+
     if (debug)
     {
         forAll(regionToCellZone, regionI)
@@ -2264,6 +2325,7 @@ void Foam::meshRefinement::getIntersections
 void Foam::meshRefinement::zonify
 (
     const bool allowFreeStandingZoneFaces,
+    const label backgroundZoneID,
     const pointField& locationsInMesh,
     const wordList& zonesInMesh,
 
@@ -2486,25 +2548,23 @@ void Foam::meshRefinement::zonify
     }
 
 
-    // 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)
         {
@@ -4269,6 +4329,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     zonify
     (
         allowFreeStandingZoneFaces,
+        -1,                 // Set all cells with cellToZone -2 to -1
         locationsInMesh,
         zonesInMesh,