ENH: specification of selection of cellZone in snappyHexMeshDict.

Now allows specification of inside point to help leaky surfaces.

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -143,11 +143,19 @@ castellatedMeshControls
                 }
             }
 
-            // Optional angle to detect small-large cell situation perpendicular
+            //- Optional angle to detect small-large cell situation
-            // to the surface. Is the angle of face w.r.t the local surface
+            //  perpendicular to the surface. Is the angle of face w.r.t.
-            // normal. Use on flat(ish) surfaces only. Otherwise
+            //  the local surface normal. Use on flat(ish) surfaces only.
-            // leave out or set to negative number.
+            //  Otherwise leave out or set to negative number.
             //perpendicularAngle 10;
+
+
+            //- Optional faceZone and (for closed surface) cellZone with
+            //  how to select the cells that are in the cellZone
+            //  (inside / outside / specified insidePoint)
+	    //faceZone sphere;
+            //cellZone sphere;
+            //cellZoneInside inside;  //outside/insidePoint
         }
     }
 

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

@@ -258,6 +258,8 @@ void Foam::meshRefinement::checkData()
     meshCutter_.checkRefinementLevels(1, labelList(0));
 
 
+    label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
+
     Pout<< "meshRefinement::checkData() : Checking synchronization."
         << endl;
 
@@ -267,7 +269,7 @@ void Foam::meshRefinement::checkData()
         pointField::subList boundaryFc
         (
             mesh_.faceCentres(),
-            mesh_.nFaces()-mesh_.nInternalFaces(),
+            nBnd,
             mesh_.nInternalFaces()
         );
 
@@ -292,8 +294,8 @@ void Foam::meshRefinement::checkData()
     // Check meshRefinement
     {
         // Get boundary face centre and level. Coupled aware.
-        labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+        labelList neiLevel(nBnd);
-        pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+        pointField neiCc(nBnd);
         calcNeighbourData(neiLevel, neiCc);
 
         // Collect segments we want to test for
@@ -327,11 +329,22 @@ void Foam::meshRefinement::checkData()
                 surfaceLevel
             );
         }
+        // Get the coupled hit
+        labelList neiHit
+        (
+            SubList<label>
+            (
+                surfaceHit,
+                nBnd,
+                mesh_.nInternalFaces()
+            )
+        );
+        syncTools::swapBoundaryFaceList(mesh_, neiHit, false);
 
         // Check
         forAll(surfaceHit, faceI)
         {
-            if (surfaceHit[faceI] != surfaceIndex_[faceI])
+            if (surfaceIndex_[faceI] != surfaceHit[faceI])
             {
                 if (mesh_.isInternalFace(faceI))
                 {
@@ -346,7 +359,11 @@ void Foam::meshRefinement::checkData()
                         << mesh_.cellCentres()[mesh_.faceNeighbour()[faceI]]
                         << endl;
                 }
-                else
+                else if
+                (
+                    surfaceIndex_[faceI]
+                 != neiHit[faceI-mesh_.nInternalFaces()]
+                )
                 {
                     WarningIn("meshRefinement::checkData()")
                         << "Boundary face:" << faceI
@@ -355,6 +372,7 @@ void Foam::meshRefinement::checkData()
                         << " current:" << surfaceHit[faceI]
                         << " ownCc:"
                         << mesh_.cellCentres()[mesh_.faceOwner()[faceI]]
+                        << " end:" << end[faceI]
                         << endl;
                 }
             }

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -436,6 +436,16 @@ private:
                 labelList& cellToZone
             ) const;
 
+            //- Finds zone per cell for cells inside named surfaces that have
+            //  an inside point specified.
+            void findCellZoneInsideWalk
+            (
+                const labelList& locationSurfaces,
+                const labelList& namedSurfaceIndex,
+                const labelList& surfaceToCellZone,
+                labelList& cellToZone
+            ) const;
+
             //- Determines cell zone from cell region information.
             bool calcRegionToZone
             (

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

@@ -1066,6 +1066,106 @@ void Foam::meshRefinement::findCellZoneGeometric
         false
     );
 }
+//XXXXXXXXX
+void Foam::meshRefinement::findCellZoneInsideWalk
+(
+    const labelList& locationSurfaces,  // indices of surfaces with inside point
+    const labelList& namedSurfaceIndex, // per face index of named surface
+    const labelList& surfaceToCellZone, // cell zone index per surface
+
+    labelList& cellToZone
+) const
+{
+    // Analyse regions. Reuse regionsplit
+    boolList blockedFace(mesh_.nFaces());
+
+    forAll(namedSurfaceIndex, faceI)
+    {
+        if (namedSurfaceIndex[faceI] == -1)
+        {
+            blockedFace[faceI] = false;
+        }
+        else
+        {
+            blockedFace[faceI] = true;
+        }
+    }
+    // No need to sync since namedSurfaceIndex already is synced
+
+    // Set region per cell based on walking
+    regionSplit cellRegion(mesh_, blockedFace);
+    blockedFace.clear();
+
+
+    // For all locationSurface find the cell
+    forAll(locationSurfaces, i)
+    {
+        label surfI = locationSurfaces[i];
+        const point& insidePoint = surfaces_.zoneInsidePoints()[surfI];
+
+        Info<< "For surface " << surfaces_.names()[surfI]
+            << " finding inside point " << insidePoint
+            << endl;
+
+        // Find the region containing the insidePoint
+        label keepRegionI = -1;
+
+        label cellI = mesh_.findCell(insidePoint);
+
+        if (cellI != -1)
+        {
+            keepRegionI = cellRegion[cellI];
+        }
+        reduce(keepRegionI, maxOp<label>());
+
+        Info<< "For surface " << surfaces_.names()[surfI]
+            << " found point " << insidePoint << " in cell " << cellI
+            << " in global region " << keepRegionI
+            << " out of " << cellRegion.nRegions() << " regions." << endl;
+
+        if (keepRegionI == -1)
+        {
+            FatalErrorIn
+            (
+                "meshRefinement::findCellZoneInsideWalk"
+                "(const labelList&, const labelList&"
+                ", const labelList&, const labelList&)"
+            )   << "Point " << insidePoint
+                << " is not inside the mesh." << nl
+                << "Bounding box of the mesh:" << mesh_.globalData().bb()
+                << exit(FatalError);
+        }
+
+        // Set all cells with this region
+        forAll(cellRegion, cellI)
+        {
+            if (cellRegion[cellI] == keepRegionI)
+            {
+                if (cellToZone[cellI] == -2)
+                {
+                    cellToZone[cellI] = surfaceToCellZone[surfI];
+                }
+                else if (cellToZone[cellI] != surfaceToCellZone[surfI])
+                {
+                    WarningIn
+                    (
+                        "meshRefinement::findCellZoneInsideWalk"
+                        "(const labelList&, const labelList&"
+                        ", const labelList&, const labelList&)"
+                    )   << "Cell " << cellI
+                        << " at " << mesh_.cellCentres()[cellI]
+                        << " is inside surface " << surfaces_.names()[surfI]
+                        << " but already marked as being in zone "
+                        << cellToZone[cellI] << endl
+                        << "This can happen if your surfaces are not"
+                        << " (sufficiently) closed."
+                        << endl;
+                }
+            }
+        }
+    }
+}
+//XXXXXXXXX
 
 
 bool Foam::meshRefinement::calcRegionToZone
@@ -2297,9 +2397,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     // Put the cells into the correct zone
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Closed surfaces with cellZone specified.
-    labelList closedNamedSurfaces(surfaces_.getClosedNamedSurfaces());
-
     // Zone per cell:
     // -2 : unset
     // -1 : not in any zone
@@ -2310,6 +2407,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     // Set using geometric test
     // ~~~~~~~~~~~~~~~~~~~~~~~~
 
+    // Closed surfaces with cellZone specified.
+    labelList closedNamedSurfaces(surfaces_.getClosedNamedSurfaces());
+
     if (closedNamedSurfaces.size())
     {
         Info<< "Found " << closedNamedSurfaces.size()
@@ -2319,17 +2419,41 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 
         findCellZoneGeometric
         (
-            closedNamedSurfaces,   // indices of closed surfaces
+            closedNamedSurfaces,    // indices of closed surfaces
-            namedSurfaceIndex,     // per face index of named surface
+            namedSurfaceIndex,      // per face index of named surface
-            surfaceToCellZone,     // cell zone index per surface
+            surfaceToCellZone,      // cell zone index per surface
+
             cellToZone
         );
     }
 
+
+    // Set using provided locations
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList locationSurfaces(surfaces_.getInsidePointNamedSurfaces());
+    if (locationSurfaces.size())
+    {
+        Info<< "Found " << locationSurfaces.size()
+            << " named surfaces with a provided inside point."
+            << " Assigning cells inside these surfaces"
+            << " to the corresponding cellZone."
+            << nl << endl;
+
+        findCellZoneInsideWalk
+        (
+            locationSurfaces,       // indices of closed surfaces
+            namedSurfaceIndex,      // per face index of named surface
+            surfaceToCellZone,      // cell zone index per surface
+
+            cellToZone
+        );
+    }
+
+
     // Set using walking
     // ~~~~~~~~~~~~~~~~~
 
-    //if (!allowFreeStandingZoneFaces)
     {
         Info<< "Walking from location-in-mesh " << keepPoint
             << " to assign cellZones "
@@ -2341,6 +2465,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
             keepPoint,
             namedSurfaceIndex,
             surfaceToCellZone,
+
             cellToZone
         );
     }

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -32,6 +32,21 @@ License
 #include "searchableSurfacesQueries.H"
 #include "UPtrList.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<>
+const char*
+Foam::NamedEnum<Foam::refinementSurfaces::areaSelectionAlgo, 4>::names[] =
+{
+    "inside",
+    "outside",
+    "insidePoint",
+    "none"
+};
+
+const Foam::NamedEnum<Foam::refinementSurfaces::areaSelectionAlgo, 4>
+    Foam::refinementSurfaces::areaSelectionAlgoNames;
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -47,7 +62,8 @@ Foam::refinementSurfaces::refinementSurfaces
     names_(surfaceDicts.size()),
     faceZoneNames_(surfaceDicts.size()),
     cellZoneNames_(surfaceDicts.size()),
-    zoneInside_(surfaceDicts.size()),
+    zoneInside_(surfaceDicts.size(), NONE),
+    zoneInsidePoints_(surfaceDicts.size()),
     regionOffset_(surfaceDicts.size())
 {
     labelList globalMinLevel(surfaceDicts.size(), 0);
@@ -74,19 +90,51 @@ Foam::refinementSurfaces::refinementSurfaces
         globalMaxLevel[surfI] = readLabel(dict.lookup("maxRefinementLevel"));
 
         // Global zone names per surface
-        if (dict.found("faceZone"))
+        if (dict.readIfPresent("faceZone", faceZoneNames_[surfI]))
         {
-            dict.lookup("faceZone") >> faceZoneNames_[surfI];
+            // Read optional entry to determine inside of faceZone
-            bool hasSide = dict.readIfPresent("zoneInside", zoneInside_[surfI]);
+
+            word method;
+            bool hasSide = dict.readIfPresent("cellZoneInside", method);
+            if (hasSide)
+            {
+                zoneInside_[surfI] = areaSelectionAlgoNames[method];
+                if (zoneInside_[surfI] == INSIDEPOINT)
+                {
+                    dict.lookup("insidePoint") >> zoneInsidePoints_[surfI];
+                }
+            }
+            else
+            {
+                // Check old syntax
+                bool inside;
+                if (dict.readIfPresent("zoneInside", inside))
+                {
+                    hasSide = true;
+                    zoneInside_[surfI] = (inside ? INSIDE : OUTSIDE);
+                }
+            }
+
+            // Read optional cellZone name
+
             if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
             {
-                if (hasSide && !allGeometry_[surfaces_[surfI]].hasVolumeType())
+                if
+                (
+                    (
+                        zoneInside_[surfI] == INSIDE
+                     || zoneInside_[surfI] == OUTSIDE
+                    )
+                 && !allGeometry_[surfaces_[surfI]].hasVolumeType()
+                )
                 {
                     IOWarningIn
                     (
                         "refinementSurfaces::refinementSurfaces(..)",
                         dict
-                    )   << "Unused entry zoneInside for faceZone "
+                    )   << "Illegal entry zoneInside "
+                        << areaSelectionAlgoNames[zoneInside_[surfI]]
+                        << " for faceZone "
                         << faceZoneNames_[surfI]
                         << " since surface " << names_[surfI]
                         << " is not closed." << endl;
@@ -282,7 +330,8 @@ Foam::refinementSurfaces::refinementSurfaces
     names_(surfacesDict.size()),
     faceZoneNames_(surfacesDict.size()),
     cellZoneNames_(surfacesDict.size()),
-    zoneInside_(surfacesDict.size()),
+    zoneInside_(surfacesDict.size(), NONE),
+    zoneInsidePoints_(surfacesDict.size()),
     regionOffset_(surfacesDict.size())
 {
     // Wilcard specification : loop over all surface, all regions
@@ -305,7 +354,7 @@ Foam::refinementSurfaces::refinementSurfaces
     names_.setSize(surfI);
     faceZoneNames_.setSize(surfI);
     cellZoneNames_.setSize(surfI);
-    zoneInside_.setSize(surfI);
+    zoneInside_.setSize(surfI, NONE);
     regionOffset_.setSize(surfI);
 
     labelList globalMinLevel(surfI, 0);
@@ -332,19 +381,42 @@ Foam::refinementSurfaces::refinementSurfaces
             globalMaxLevel[surfI] = refLevel[1];
 
             // Global zone names per surface
-            if (dict.found("faceZone"))
+            if (dict.readIfPresent("faceZone", faceZoneNames_[surfI]))
             {
-                dict.lookup("faceZone") >> faceZoneNames_[surfI];
+                // Read optional entry to determine inside of faceZone
-                bool hasSide = dict.readIfPresent
+
-                (
+                word method;
-                    "zoneInside",
+                bool hasSide = dict.readIfPresent("cellZoneInside", method);
-                    zoneInside_[surfI]
+                if (hasSide)
-                );
+                {
+                    zoneInside_[surfI] = areaSelectionAlgoNames[method];
+                    if (zoneInside_[surfI] == INSIDEPOINT)
+                    {
+                        dict.lookup("insidePoint") >> zoneInsidePoints_[surfI];
+                    }
+
+                }
+                else
+                {
+                    // Check old syntax
+                    bool inside;
+                    if (dict.readIfPresent("zoneInside", inside))
+                    {
+                        hasSide = true;
+                        zoneInside_[surfI] = (inside ? INSIDE : OUTSIDE);
+                    }
+                }
+
+                // Read optional cellZone name
+
                 if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
                 {
                     if
                     (
-                        hasSide
+                        (
+                            zoneInside_[surfI] == INSIDE
+                         || zoneInside_[surfI] == OUTSIDE
+                        )
                     && !allGeometry_[surfaces_[surfI]].hasVolumeType()
                     )
                     {
@@ -352,7 +424,9 @@ Foam::refinementSurfaces::refinementSurfaces
                         (
                             "refinementSurfaces::refinementSurfaces(..)",
                             dict
-                        )   << "Unused entry zoneInside for faceZone "
+                        )   << "Illegal entry zoneInside "
+                            << areaSelectionAlgoNames[zoneInside_[surfI]]
+                            << " for faceZone "
                             << faceZoneNames_[surfI]
                             << " since surface " << names_[surfI]
                             << " is not closed." << endl;
@@ -533,12 +607,36 @@ Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const
     label closedI = 0;
     forAll(cellZoneNames_, surfI)
     {
-        if (cellZoneNames_[surfI].size())
+        if
+        (
+            cellZoneNames_[surfI].size()
+         && (
+                zoneInside_[surfI] == INSIDE
+             || zoneInside_[surfI] == OUTSIDE
+            )
+         && allGeometry_[surfaces_[surfI]].hasVolumeType()
+        )
         {
-            if (allGeometry_[surfaces_[surfI]].hasVolumeType())
+            closed[closedI++] = surfI;
-            {
+        }
-                closed[closedI++] = surfI;
+    }
-            }
+    closed.setSize(closedI);
+
+    return closed;
+}
+
+
+// Get indices of named surfaces with a 
+Foam::labelList Foam::refinementSurfaces::getInsidePointNamedSurfaces() const
+{
+    labelList closed(cellZoneNames_.size());
+
+    label closedI = 0;
+    forAll(cellZoneNames_, surfI)
+    {
+        if (cellZoneNames_[surfI].size() && zoneInside_[surfI] == INSIDEPOINT)
+        {
+            closed[closedI++] = surfI;
         }
     }
     closed.setSize(closedI);
@@ -1199,6 +1297,16 @@ void Foam::refinementSurfaces::findInside
     {
         label surfI = testSurfaces[i];
 
+        if (zoneInside_[surfI] != INSIDE && zoneInside_[surfI] != OUTSIDE)
+        {
+            FatalErrorIn("refinementSurfaces::findInside(..)")
+                << "Trying to use surface "
+                << allGeometry_[surfaces_[surfI]].name()
+                << " which has non-geometric inside selection method "
+                << areaSelectionAlgoNames[zoneInside_[surfI]]
+                << exit(FatalError);
+        }
+
         if (allGeometry_[surfaces_[surfI]].hasVolumeType())
         {
             List<searchableSurface::volumeType> volType;
@@ -1212,11 +1320,11 @@ void Foam::refinementSurfaces::findInside
                     (
                         (
                             volType[pointI] == triSurfaceMesh::INSIDE
-                         && zoneInside_[surfI]
+                         && zoneInside_[surfI] == INSIDE
                         )
                      || (
                             volType[pointI] == triSurfaceMesh::OUTSIDE
-                         && !zoneInside_[surfI]
+                         && zoneInside_[surfI] == OUTSIDE
                         )
                     )
                     {

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -57,6 +57,21 @@ class triSurfaceMesh;
 
 class refinementSurfaces
 {
+public:
+
+    //- Types of selection of area
+    enum areaSelectionAlgo
+    {
+        INSIDE,
+        OUTSIDE,
+        INSIDEPOINT,
+        NONE
+    };
+
+    static const NamedEnum<areaSelectionAlgo, 4> areaSelectionAlgoNames;
+
+private:
+
     // Private data
 
         //- Reference to all geometry.
@@ -75,9 +90,12 @@ class refinementSurfaces
         wordList cellZoneNames_;
 
         //- Per 'interface' surface : (only used if surface is closed)
-        //  whether to zone cells inside or outside surface.
+        //  How to select zone cells : surface inside or outside or given
-        boolList zoneInside_;
+        //  inside location.
+        List<areaSelectionAlgo> zoneInside_;
 
+        //- If zoneInside=location gives the corresponding inside point
+        pointField zoneInsidePoints_;
 
         //- From local region number to global region number
         labelList regionOffset_;
@@ -159,9 +177,20 @@ public:
             //- Get indices of named surfaces (surfaces with faceZoneName)
             labelList getNamedSurfaces() const;
 
-            //- Get indices of closed surfaces with a cellZone
+            //- Get indices of surfaces with a cellZone that are closed and
+            //  have 'inside' or 'outside' selection.
             labelList getClosedNamedSurfaces() const;
 
+            //- Get indices of surfaces with a cellZone that have 'insidePoint'
+            //  section.
+            labelList getInsidePointNamedSurfaces() const;
+
+            //- Get specified inside locations for surfaces with a cellZone
+            const pointField& zoneInsidePoints() const
+            {
+                return zoneInsidePoints_;
+            }
+
             //- From local region number to global region number
             const labelList& regionOffset() const
             {

tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/snappyHexMeshDict

@@ -130,7 +130,7 @@ castellatedMeshControls
 
 	    faceZone bottomAir;
             cellZone bottomAir;
-            zoneInside true;
+            cellZoneInside inside;
         }
 
         topAir
@@ -140,7 +140,7 @@ castellatedMeshControls
 
 	    faceZone topAir;
             cellZone topAir;
-            zoneInside true;
+            cellZoneInside inside;
         }
 
         heater
@@ -150,7 +150,7 @@ castellatedMeshControls
 
 	    faceZone heater;
             cellZone heater;
-            zoneInside true;
+            cellZoneInside inside;
         }
 
         leftSolid
@@ -160,7 +160,7 @@ castellatedMeshControls
 
 	    faceZone leftSolid;
             cellZone leftSolid;
-            zoneInside true;
+            cellZoneInside inside;
         }
 
         rightSolid
@@ -170,7 +170,7 @@ castellatedMeshControls
 
 	    faceZone rightSolid;
             cellZone rightSolid;
-            zoneInside true;
+            cellZoneInside inside;
         }
     }