From 133149a1dd221a309084a99bf4b69a517fb901a4 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 14 Sep 2023 16:12:53 +0100
Subject: [PATCH] ENH: snappyHexMesh: mutliple inside points. Fixes #2981

---
 etc/caseDicts/annotated/snappyHexMeshDict     |  2 ++
 .../meshRefinement/meshRefinementBaffles.C    | 23 ++++++++----
 .../refinementSurfaces/surfaceZonesInfo.C     | 36 +++++++++++++++----
 .../refinementSurfaces/surfaceZonesInfo.H     | 20 +++++++----
 4 files changed, 60 insertions(+), 21 deletions(-)

diff --git a/etc/caseDicts/annotated/snappyHexMeshDict b/etc/caseDicts/annotated/snappyHexMeshDict
index c49c1f2941..35a7dd184f 100644
--- a/etc/caseDicts/annotated/snappyHexMeshDict
+++ b/etc/caseDicts/annotated/snappyHexMeshDict
@@ -225,6 +225,8 @@ castellatedMeshControls
             //cellZone sphere;
             //cellZoneInside inside;    // outside/insidePoint
             //insidePoint    (1 1 1);   // if (cellZoneInside == insidePoint)
+            //   or alternative multiple insidePoints:
+            //insidePoints   ((1 1 1));   // if (cellZoneInside == insidePoint)
 
             //- Optional specification of what to do with faceZone faces:
             //      internal : keep them as internal faces (default)
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
index 2c99fecf15..57822b772c 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -3344,24 +3344,33 @@ void Foam::meshRefinement::zonify
             << nl << endl;
 
         // Collect per surface the -insidePoint -cellZone
-        pointField insidePoints(locationSurfaces.size());
-        labelList insidePointCellZoneIDs(locationSurfaces.size(), -1);
+        // Usually only a single inside point per surface so no clever
+        // counting - just use DynamicField
+        DynamicField<point> insidePoints(locationSurfaces.size());
+        DynamicList<label> insidePointCellZoneIDs(locationSurfaces.size());
         forAll(locationSurfaces, i)
         {
-            label surfI = locationSurfaces[i];
-            insidePoints[i] = surfZones[surfI].zoneInsidePoint();
+            const label surfI = locationSurfaces[i];
+            const auto& surfInsidePoints = surfZones[surfI].zoneInsidePoints();
 
             const word& name = surfZones[surfI].cellZoneName();
+            label zoneID = -1;
             if (name != "none")
             {
-                label zoneID = mesh_.cellZones().findZoneID(name);
+                zoneID = mesh_.cellZones().findZoneID(name);
                 if (zoneID == -1)
                 {
                     FatalErrorInFunction
-                        << "problem"
+                        << "Specified non-existing cellZone " << name
+                        << " for surface " << surfaces_.names()[surfI]
                         << abort(FatalError);
                 }
-                insidePointCellZoneIDs[i] = zoneID;
+            }
+
+            for (const auto& pt : surfInsidePoints)
+            {
+                insidePoints.append(pt);
+                insidePointCellZoneIDs.append(zoneID);
             }
         }
 
diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/surfaceZonesInfo.C b/src/mesh/snappyHexMesh/refinementSurfaces/surfaceZonesInfo.C
index 70e267c95f..dae2c0467e 100644
--- a/src/mesh/snappyHexMesh/refinementSurfaces/surfaceZonesInfo.C
+++ b/src/mesh/snappyHexMesh/refinementSurfaces/surfaceZonesInfo.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2015 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,7 +83,7 @@ Foam::surfaceZonesInfo::surfaceZonesInfo
     faceZoneNames_(),
     cellZoneName_(),
     zoneInside_(NONE),
-    zoneInsidePoint_(point::min),
+    zoneInsidePoints_(),
     faceType_(INTERNAL)
 {
     const label nRegions = surface.regions().size();
@@ -156,9 +156,31 @@ Foam::surfaceZonesInfo::surfaceZonesInfo
             zoneInside_ = areaSelectionAlgoNames[method];
             if (zoneInside_ == INSIDEPOINT)
             {
-                surfacesDict.readEntry("insidePoint", zoneInsidePoint_);
-            }
+                const bool foundPoints = surfacesDict.readIfPresent
+                (
+                    "insidePoints",
+                    zoneInsidePoints_,
+                    keyType::LITERAL
+                );
 
+                if (foundPoints)
+                {
+                    if (surfacesDict.found("insidePoint", keyType::LITERAL))
+                    {
+                        FatalIOErrorInFunction(surfacesDict)
+                            << "Cannot supply both 'insidePoint'"
+                            << " and 'insidePoints'" << exit(FatalIOError);
+                    }
+                }
+                else
+                {
+                    zoneInsidePoints_ = pointField
+                    (
+                        1,
+                        surfacesDict.get<point>("insidePoint", keyType::LITERAL)
+                    );
+                }
+            }
         }
         else
         {
@@ -216,14 +238,14 @@ Foam::surfaceZonesInfo::surfaceZonesInfo
     const wordList& faceZoneNames,
     const word& cellZoneName,
     const areaSelectionAlgo& zoneInside,
-    const point& zoneInsidePoint,
+    const pointField& zoneInsidePoints,
     const faceZoneType& faceType
 )
 :
     faceZoneNames_(faceZoneNames),
     cellZoneName_(cellZoneName),
     zoneInside_(zoneInside),
-    zoneInsidePoint_(zoneInsidePoint),
+    zoneInsidePoints_(zoneInsidePoints),
     faceType_(faceType)
 {}
 
@@ -233,7 +255,7 @@ Foam::surfaceZonesInfo::surfaceZonesInfo(const surfaceZonesInfo& surfZone)
     faceZoneNames_(surfZone.faceZoneNames()),
     cellZoneName_(surfZone.cellZoneName()),
     zoneInside_(surfZone.zoneInside()),
-    zoneInsidePoint_(surfZone.zoneInsidePoint()),
+    zoneInsidePoints_(surfZone.zoneInsidePoints()),
     faceType_(surfZone.faceType())
 {}
 
diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/surfaceZonesInfo.H b/src/mesh/snappyHexMesh/refinementSurfaces/surfaceZonesInfo.H
index f22e635620..fe67e38e8a 100644
--- a/src/mesh/snappyHexMesh/refinementSurfaces/surfaceZonesInfo.H
+++ b/src/mesh/snappyHexMesh/refinementSurfaces/surfaceZonesInfo.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2014 OpenFOAM Foundation
-    Copyright (C) 2015 OpenCFD Ltd.
+    Copyright (C) 2015,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,7 +38,7 @@ SourceFiles
 #define surfaceZonesInfo_H
 
 #include "Enum.H"
-#include "point.H"
+#include "pointField.H"
 #include "word.H"
 #include "PtrList.H"
 #include "labelList.H"
@@ -110,7 +110,7 @@ private:
         areaSelectionAlgo zoneInside_;
 
         //- If zoneInside=location gives the corresponding inside point
-        point zoneInsidePoint_;
+        pointField zoneInsidePoints_;
 
         //- Per 'interface' surface :
         //  What to do with outside
@@ -142,7 +142,7 @@ public:
             const wordList& faceZoneNames,
             const word& cellZoneNames,
             const areaSelectionAlgo& zoneInside,
-            const point& zoneInsidePoints,
+            const pointField& zoneInsidePoints,
             const faceZoneType& faceType
         );
 
@@ -179,10 +179,16 @@ public:
                 return zoneInside_;
             }
 
-            //- Get specified inside locations for surfaces with a cellZone
-            const point& zoneInsidePoint() const
+            //- Get specified inside location for surfaces with a cellZone
+            const point zoneInsidePoint() const
             {
-                return zoneInsidePoint_;
+                return zoneInsidePoints_[0];
+            }
+
+            //- Get specified inside locations for surfaces with a cellZone
+            const pointField& zoneInsidePoints() const
+            {
+                return zoneInsidePoints_;
             }
 
             //- How to handle face of surfaces with a faceZone