ENH: autoHexMesh: allow location-in-mesh to be made into a cellZone

Normally the location-in-mesh is equivalent to a cellZone -1. This can
now be overridden by a surface-specified cellZone as before.

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -58,6 +58,7 @@ Description
 #include "globalIndex.H"
 #include "IOmanip.H"
 #include "decompositionModel.H"
+#include "fvMeshTools.H"
 
 using namespace Foam;
 
@@ -813,6 +814,7 @@ int main(int argc, char *argv[])
         readScalar(meshDict.lookup("mergeTolerance"))
     );
 
+    const Switch keepPatches(meshDict.lookupOrDefault("keepPatches", false));
 
 
     // Read decomposePar dictionary
@@ -1517,6 +1519,12 @@ int main(int argc, char *argv[])
             motionDict
         );
 
+        // Remove zero sized patches originating from faceZones
+        if (!keepPatches && !wantSnap && !wantLayers)
+        {
+            fvMeshTools::removeEmptyPatches(mesh, true);
+        }
+
         writeMesh
         (
             "Refined mesh",
@@ -1559,6 +1567,12 @@ int main(int argc, char *argv[])
             snapParams
         );
 
+        // Remove zero sized patches originating from faceZones
+        if (!keepPatches && !wantLayers)
+        {
+            fvMeshTools::removeEmptyPatches(mesh, true);
+        }
+
         writeMesh
         (
             "Snapped mesh",
@@ -1609,6 +1623,12 @@ int main(int argc, char *argv[])
             distributor
         );
 
+        // Remove zero sized patches originating from faceZones
+        if (!keepPatches)
+        {
+            fvMeshTools::removeEmptyPatches(mesh, true);
+        }
+
         writeMesh
         (
             "Layer mesh",
@@ -1622,6 +1642,7 @@ int main(int argc, char *argv[])
     }
 
 
+
     {
         // Check final mesh
         Info<< "Checking final mesh ..." << endl;

src/dynamicMesh/fvMeshTools/fvMeshTools.C

@@ -355,6 +355,69 @@ void Foam::fvMeshTools::reorderPatches
 }
 
 
+Foam::labelList Foam::fvMeshTools::removeEmptyPatches
+(
+    fvMesh& mesh,
+    const bool validBoundary
+)
+{
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    labelList newToOld(pbm.size());
+    labelList oldToNew(pbm.size(), -1);
+    label newI = 0;
+
+
+    // Assumes all non-coupled boundaries are on all processors!
+    forAll(pbm, patchI)
+    {
+        const polyPatch& pp = pbm[patchI];
+
+        if (!isA<processorPolyPatch>(pp))
+        {
+            label nFaces = pp.size();
+            if (validBoundary)
+            {
+                reduce(nFaces, sumOp<label>());
+            }
+
+            if (nFaces > 0)
+            {
+                newToOld[newI] = patchI;
+                oldToNew[patchI] = newI++;
+            }
+        }
+    }
+
+    // Same for processor patches (but need no reduction)
+    forAll(pbm, patchI)
+    {
+        const polyPatch& pp = pbm[patchI];
+
+        if (isA<processorPolyPatch>(pp) && pp.size())
+        {
+            newToOld[newI] = patchI;
+            oldToNew[patchI] = newI++;
+        }
+    }
+
+    newToOld.setSize(newI);
+
+    // Move all deleteable patches to the end
+    forAll(oldToNew, patchI)
+    {
+        if (oldToNew[patchI] == -1)
+        {
+            oldToNew[patchI] = newI++;
+        }
+    }
+
+    reorderPatches(mesh, oldToNew, newToOld.size(), validBoundary);
+
+    return newToOld;
+}
+
+
 Foam::autoPtr<Foam::fvMesh> Foam::fvMeshTools::newMesh
 (
     const IOobject& io,

src/dynamicMesh/fvMeshTools/fvMeshTools.H

@@ -113,7 +113,7 @@ public:
     static void zeroPatchFields(fvMesh& mesh, const label patchI);
 
     //- Reorder and remove trailing patches. If validBoundary call is parallel
-    //  synced and all add the same patch with same settings
+    //  synced
     static void reorderPatches
     (
         fvMesh&,
@@ -122,6 +122,11 @@ public:
         const bool validBoundary
     );
 
+    //- Remove zero sized patches. All but processor patches are
+    //  assumed to be present on all processors (so size will be reduced
+    //  if validBoundary). Return map from new
+    //  to old patches
+    static labelList removeEmptyPatches(fvMesh&, const bool validBoundary);
 
     //- Read mesh or create dummy mesh (0 cells, >0 patches). Works in two
     //  modes according to masterOnlyReading:
@@ -133,7 +138,6 @@ public:
         const IOobject& io,
         const bool masterOnlyReading
     );
-
 };
 
 

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -318,6 +318,47 @@ void Foam::meshRefinement::getBafflePatches
             namedSurfaceIndex,
             posOrientation
         );
+
+
+        // Some stats
+        if (debug)
+        {
+            label nZones = gMax(cellToZone)+1;
+
+            label nUnvisited = 0;
+            label nBackgroundCells = 0;
+            labelList nZoneCells(nZones, 0);
+            forAll(cellToZone, cellI)
+            {
+                label zoneI = cellToZone[cellI];
+                if (zoneI >= 0)
+                {
+                    nZoneCells[zoneI]++;
+                }
+                else if (zoneI == -1)
+                {
+                    nBackgroundCells++;
+                }
+                else if (zoneI == -2)
+                {
+                    nUnvisited++;
+                }
+                else
+                {
+                    FatalErrorIn("meshRefinement::getBafflePatches()")
+                        << "problem" << exit(FatalError);
+                }
+            }
+            reduce(nUnvisited, sumOp<label>());
+            reduce(nBackgroundCells, sumOp<label>());
+            forAll(nZoneCells, zoneI)
+            {
+                reduce(nZoneCells[zoneI], sumOp<label>());
+            }
+            Info<< "nUnvisited      :" << nUnvisited << endl;
+            Info<< "nBackgroundCells:" << nBackgroundCells << endl;
+            Info<< "nZoneCells      :" << nZoneCells << endl;
+        }
     }
 
 
@@ -342,9 +383,16 @@ void Foam::meshRefinement::getBafflePatches
     syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
 
 
-    // 3. Baffle all faces on outside of unvisited cells (cellToZone = -2)
+    // 3. Baffle all boundary faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Baffle all boundary faces except those on outside of unvisited cells
+    // (cellToZone = -2)
     // Use patch according to globalRegion1,2
+    // Note: the behaviour is slightly different from version3.0 and earlier
+    //       in that it will not create baffles which are outside the meshed
+    //       domain. On a medium testcase (motorBike tutorial geometry) this
+    //       selects about 20 cells less (out of 120000). These cells are where
+    //       there might e.g. be a single cell which is fully unreachable.
 
     ownPatch.setSize(mesh_.nFaces());
     ownPatch = -1;
@@ -365,12 +413,7 @@ void Foam::meshRefinement::getBafflePatches
             neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
         }
 
-
+        if (ownZone != -2 || neiZone != -2)
-        if
-        (
-            (ownZone == -2 && neiZone != -2)
-         || (ownZone != -2 && neiZone == -2)
-        )
         {
             if (globalRegion1[faceI] != -1)
             {
@@ -383,7 +426,6 @@ void Foam::meshRefinement::getBafflePatches
         }
     }
 
-
     // No need to parallel sync since intersection data (surfaceIndex_ etc.)
     // already guaranteed to be synced...
     // However:
@@ -1390,14 +1432,21 @@ void Foam::meshRefinement::findCellZoneGeometric
 
     forAll(insideSurfaces, cellI)
     {
-        if (cellToZone[cellI] == -2)
+        label surfI = insideSurfaces[cellI];
-        {
-            label surfI = insideSurfaces[cellI];
 
-            if (surfI != -1)
+        if (surfI != -1)
+        {
+            if (cellToZone[cellI] == -2)
             {
                 cellToZone[cellI] = surfaceToCellZone[surfI];
             }
+            else if (cellToZone[cellI] == -1)
+            {
+                // ? Allow named surface to override background zone (-1)
+                // This is used in the multiRegionHeater tutorial where the
+                // locationInMesh is inside a named surface.
+                cellToZone[cellI] = surfaceToCellZone[surfI];
+            }
         }
     }