From 3973cf7a839a8d953aa457864990b288b0cafc0e Mon Sep 17 00:00:00 2001
From: Mattijs Janssens <ext-mjanssens@esi-group.com>
Date: Wed, 13 Mar 2024 14:44:13 +0000
Subject: [PATCH] Feature overlapping zones

---
 .../meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C |  94 +-
 .../meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H |  15 +-
 src/dynamicMesh/polyMeshAdder/polyMeshAdder.C | 123 ++-
 src/meshTools/polyTopoChange/polyTopoChange.C | 900 ++++++++++++++++--
 src/meshTools/polyTopoChange/polyTopoChange.H | 126 ++-
 .../polyTopoChange/polyTopoChangeTemplates.C  |  26 +-
 .../mesh/parallel/filter/system/fvOptions     |  23 +
 .../mesh/parallel/filter/system/topoSetDict   |  36 +-
 8 files changed, 1160 insertions(+), 183 deletions(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C b/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C
index 29ead1a45e6..9fa134061a9 100644
--- a/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2023 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,17 +83,77 @@ void Foam::ZoneMesh<ZoneType, MeshType>::calcZoneMap() const
         map.reserve(this->totalSize());
 
         const PtrList<ZoneType>& zones = *this;
+
         label zonei = 0;
 
         for (const ZoneType& zn : zones)
         {
             for (const label id : static_cast<const labelList&>(zn))
             {
-                map.insert(id, zonei);
+                //map.insert(id, zonei);
+                const auto fnd = map.cfind(id);
+                if (!fnd)
+                {
+                    map.insert(id, zonei);
+                }
+                else if (fnd() != zonei)
+                {
+                    // Multiple zones for same id
+
+                    if (!additionalMapPtr_)
+                    {
+                        // First occurrence
+                        label maxIndex = -1;
+                        for (const ZoneType& zn : zones)
+                        {
+                            for
+                            (
+                                const label id
+                              : static_cast<const labelList&>(zn)
+                            )
+                            {
+                                maxIndex = max(maxIndex, id);
+                            }
+                        }
+                        additionalMapPtr_.reset(new labelListList(maxIndex+1));
+                    }
+                    auto& additionalMap = *additionalMapPtr_;
+                    additionalMap[id].push_uniq(zonei);
+                }
             }
 
             ++zonei;
         }
+
+        // Sort such that map contains lowest zoneID
+        if (additionalMapPtr_)
+        {
+            auto& additionalMap = *additionalMapPtr_;
+            forAll(additionalMap, id)
+            {
+                labelList& zones = additionalMap[id];
+
+                if (zones.size())
+                {
+                    stableSort(zones);
+                    const label zonei = map[id];
+                    const label index = findLower(zones, zonei);
+                    if (index == -1)
+                    {
+                        // Already first
+                    }
+                    else
+                    {
+                        map.set(id, zones[0]);
+                        for (label i = 0; i < zones.size() && i <= index; i++)
+                        {
+                            zones[i] = zones[i+1];
+                        }
+                        zones[index+1] = zonei;
+                    }
+                }
+            }
+        }
     }
 }
 
@@ -358,6 +418,35 @@ Foam::label Foam::ZoneMesh<ZoneType, MeshType>::whichZone
 }
 
 
+template<class ZoneType, class MeshType>
+Foam::label Foam::ZoneMesh<ZoneType, MeshType>::whichZones
+(
+    const label objectIndex,
+    DynamicList<label>& zones
+) const
+{
+    zones.clear();
+    const auto fnd = zoneMap().find(objectIndex);
+    if (fnd)
+    {
+        // Add main element
+        zones.push_back(fnd());
+        if (additionalMapPtr_)
+        {
+            const auto& additionalMap = *additionalMapPtr_;
+            if (objectIndex < additionalMap.size())
+            {
+                for (const label zonei : additionalMap[objectIndex])
+                {
+                    zones.push_uniq(zonei);
+                }
+            }
+        }
+    }
+    return zones.size();
+}
+
+
 template<class ZoneType, class MeshType>
 Foam::wordList Foam::ZoneMesh<ZoneType, MeshType>::types() const
 {
@@ -814,6 +903,7 @@ template<class ZoneType, class MeshType>
 void Foam::ZoneMesh<ZoneType, MeshType>::clearLocalAddressing()
 {
     zoneMapPtr_.reset(nullptr);
+    additionalMapPtr_.reset(nullptr);
     groupIDsPtr_.reset(nullptr);
 }
 
diff --git a/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H b/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H
index 90ab8ca1542..fa6d3eb0435 100644
--- a/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2023 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -76,6 +76,10 @@ class ZoneMesh
         //- Demand-driven: map of zone labels for given element
         mutable std::unique_ptr<Map<label>> zoneMapPtr_;
 
+        //- Demand-driven: map of additional zone labels for given element.
+        //  Only populated if overlapping zones
+        mutable std::unique_ptr<labelListList> additionalMapPtr_;
+
         //- Demand-driven: list of zone ids per group
         mutable std::unique_ptr<HashTable<labelList>> groupIDsPtr_;
 
@@ -190,6 +194,15 @@ public:
         //  If object does not belong to any zones, return -1
         label whichZone(const label objectIndex) const;
 
+        //- Given a global object index, return (in argument) the zones it is
+        //  in. Returns number of zones (0 if object does not belong to any
+        //  zones)
+        label whichZones
+        (
+            const label objectIndex,
+            DynamicList<label>& zones
+        ) const;
+
         //- Return a list of zone types
         wordList types() const;
 
diff --git a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
index 85ce1edb056..d5634a7cccb 100644
--- a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
+++ b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
@@ -2698,37 +2698,14 @@ void Foam::polyMeshAdder::add
             const polyMesh& mesh = meshes[meshi];
             const cellZoneMesh& cellZones = mesh.cellZones();
 
-            // Precalc offset zones
-            labelList newZoneID(mesh.nCells(), -1);
-
-            forAll(cellZones, zonei)
-            {
-                const labelList& cellLabels = cellZones[zonei];
-
-                for (const label celli : cellLabels)
-                {
-                    if (newZoneID[celli] != -1)
-                    {
-                        WarningInFunction
-                            << "Cell:" << celli
-                            << " centre:" << mesh.cellCentres()[celli]
-                            << " is in two zones:"
-                            << cellZones[newZoneID[celli]].name()
-                            << " and " << cellZones[zonei].name() << endl
-                            << "    This is not supported."
-                            << " Continuing with first zone only." << endl;
-                    }
-                    else
-                    {
-                        newZoneID[celli] = cellZoneMap[meshi][zonei];
-                    }
-                }
-            }
-
             // Add cells in mesh order
             cellProcAddressing[meshi].setSize(mesh.nCells());
+            DynamicList<label> zones;
             for (label celli = 0; celli < mesh.nCells(); celli++)
             {
+                // Get current zones and renumber
+                cellZones.whichZones(celli, zones);
+                inplaceRenumber(cellZoneMap[meshi], zones);
                 // Add cell from cell
                 cellProcAddressing[meshi][celli] = meshMod.addCell
                 (
@@ -2736,7 +2713,7 @@ void Foam::polyMeshAdder::add
                     -1,
                     -1,
                     celli,
-                    newZoneID[celli]
+                    zones
                 );
             }
         }
@@ -2755,30 +2732,21 @@ void Foam::polyMeshAdder::add
             const pointField& points = mesh.points();
             const pointZoneMesh& pointZones = mesh.pointZones();
 
-            // Precalc offset zones
-            labelList newZoneID(points.size(), -1);
-
-            forAll(pointZones, zonei)
-            {
-                const labelList& pointLabels = pointZones[zonei];
-
-                for (const label pointi : pointLabels)
-                {
-                    newZoneID[pointi] = pointZoneMap[meshi][zonei];
-                }
-            }
-
             // Add points in mesh order
             const labelList& myUniquePoints = uniquePoints[meshi];
+            DynamicList<label> zones;
             for (const label pointi : myUniquePoints)
             {
+                // Get current zones and renumber
+                pointZones.whichZones(pointi, zones);
+                inplaceRenumber(pointZoneMap[meshi], zones);
                 // Rewrite the addressing (should already be compact) just in
                 // case the polyTopoChange does some special ordering
                 pointProcAddressing[meshi][pointi] = meshMod.addPoint
                 (
                     points[pointi],
                     pointi,
-                    newZoneID[pointi],
+                    zones,
                     true
                 );
             }
@@ -2813,21 +2781,8 @@ void Foam::polyMeshAdder::add
             const labelList& faceNeighbour = mesh.faceNeighbour();
             const faceZoneMesh& faceZones = mesh.faceZones();
 
-            // Precalc offset zones
-            labelList newZoneID(faces.size(), -1);
-            boolList zoneFlip(faces.size(), false);
-
-            forAll(faceZones, zonei)
-            {
-                const labelList& faceLabels = faceZones[zonei];
-                const boolList& flipMap = faceZones[zonei].flipMap();
-
-                forAll(faceLabels, facei)
-                {
-                    newZoneID[faceLabels[facei]] = faceZoneMap[meshi][zonei];
-                    zoneFlip[faceLabels[facei]] = flipMap[facei];
-                }
-            }
+            DynamicList<label> zones;
+            DynamicList<bool> flips;
 
             // Add faces in mesh order
 
@@ -2846,14 +2801,28 @@ void Foam::polyMeshAdder::add
                 {
                     newFace[fp] = pointProcAddressing[meshi][f[fp]];
                 }
+
+                // Get current zones and renumber
+                faceZones.whichZones(facei, zones);
+                flips.clear();
+                forAll(zones, i)
+                {
+                    const label zonei = zones[i];
+                    const label index = faceZones[zonei].localID(facei);
+                    flips.append(faceZones[zonei].flipMap()[index]);
+                }
+                inplaceRenumber(faceZoneMap[meshi], zones);
+
                 bool flipFaceFlux = false;
-                bool flip = zoneFlip[facei];
                 if (newNei < newOwn)
                 {
                     std::swap(newOwn, newNei);
                     newFace = newFace.reverseFace();
                     flipFaceFlux = !flipFaceFlux;
-                    flip = !flip;
+                    for (bool& flip : flips)
+                    {
+                        flip = !flip;
+                    }
                 }
 
                 const label newFacei = meshMod.addFace
@@ -2866,8 +2835,8 @@ void Foam::polyMeshAdder::add
                     facei,                      // masterFaceID
                     flipFaceFlux,               // flipFaceFlux
                     -1,                         // patchID
-                    newZoneID[facei],           // zoneID
-                    flip                        // zoneFlip
+                    zones,                      // zoneIDs
+                    flips                       // zoneFlips
                 );
 
                 faceProcAddressing[meshi][facei] = newFacei;
@@ -2892,6 +2861,17 @@ void Foam::polyMeshAdder::add
                 const label newOwn =
                     cellProcAddressing[meshi][faceOwner[facei]];
 
+                faceZones.whichZones(facei, zones);
+                flips.clear();
+                forAll(zones, i)
+                {
+                    const label zonei = zones[i];
+                    const label index = faceZones[zonei].localID(facei);
+                    flips.append(faceZones[zonei].flipMap()[index]);
+                }
+                inplaceRenumber(faceZoneMap[meshi], zones);
+
+
                 // Neighbour mesh face
                 const auto& nbrMesh = meshes[nbrMeshi];
                 const label nbrFacei = nbrMesh.nInternalFaces()+nbrBFacei;
@@ -2912,7 +2892,6 @@ void Foam::polyMeshAdder::add
                         newFace[fp] = pointProcAddressing[meshi][f[fp]];
                     }
                     const bool flipFaceFlux = false;
-                    const bool flip = zoneFlip[facei];
                     const label newFacei = meshMod.addFace
                     (
                         newFace,
@@ -2923,8 +2902,8 @@ void Foam::polyMeshAdder::add
                         facei,                      // masterFaceID
                         flipFaceFlux,               // flipFaceFlux
                         -1,                         // patchID
-                        newZoneID[facei],           // zoneID
-                        flip                        // zoneFlip
+                        zones,                      // zoneIDs
+                        flips                       // zoneFlips
                     );
 
                     faceProcAddressing[meshi][facei] = newFacei;
@@ -2960,6 +2939,18 @@ void Foam::polyMeshAdder::add
                                 newFace[fp] = pointProcAddressing[meshi][f[fp]];
                             }
 
+                            faceZones.whichZones(facei, zones);
+                            flips.clear();
+                            forAll(zones, i)
+                            {
+                                const label zonei = zones[i];
+                                const label index =
+                                    faceZones[zonei].localID(facei);
+                                flips.append(faceZones[zonei].flipMap()[index]);
+                            }
+                            inplaceRenumber(faceZoneMap[meshi], zones);
+
+
                             const label newFacei = meshMod.addFace
                             (
                                 newFace,
@@ -2970,8 +2961,8 @@ void Foam::polyMeshAdder::add
                                 facei,                      // masterFaceID
                                 false,                      // flipFaceFlux
                                 newPatchi,                  // patchID
-                                newZoneID[facei],           // zoneID
-                                zoneFlip[facei]             // zoneFlip
+                                zones,                      // zoneID
+                                flips                       // zoneFlip
                             );
 
                             faceProcAddressing[meshi][facei] = newFacei;
diff --git a/src/meshTools/polyTopoChange/polyTopoChange.C b/src/meshTools/polyTopoChange/polyTopoChange.C
index 9d63159f40c..822ff1efe4a 100644
--- a/src/meshTools/polyTopoChange/polyTopoChange.C
+++ b/src/meshTools/polyTopoChange/polyTopoChange.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2022,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -905,6 +905,11 @@ void Foam::polyTopoChange::reorderCompactFaces
     renumberKey(oldToNew, faceZone_);
     inplaceReorder(oldToNew, faceZoneFlip_);
     faceZoneFlip_.setCapacity(newSize);
+    if (faceAdditionalZones_.size())
+    {
+        inplaceReorder(oldToNew, faceAdditionalZones_);
+        faceAdditionalZones_.setCapacity(newSize);
+    }
 }
 
 
@@ -913,6 +918,7 @@ void Foam::polyTopoChange::shrink()
     points_.shrink();
     pointMap_.shrink();
     reversePointMap_.shrink();
+    pointAdditionalZones_.shrink();
 
     faces_.shrink();
     region_.shrink();
@@ -920,10 +926,12 @@ void Foam::polyTopoChange::shrink()
     faceNeighbour_.shrink();
     faceMap_.shrink();
     reverseFaceMap_.shrink();
+    faceAdditionalZones_.shrink();
 
     cellMap_.shrink();
     reverseCellMap_.shrink();
     cellZone_.shrink();
+    cellAdditionalZones_.shrink();
 }
 
 
@@ -1091,6 +1099,10 @@ void Foam::polyTopoChange::compact
 
         renumberKey(localPointMap, pointZone_);
         renumber(localPointMap, retiredPoints_);
+        if (pointAdditionalZones_.size())
+        {
+            reorder(localPointMap, pointAdditionalZones_);
+        }
 
         // Use map to relabel face vertices
         forAll(faces_, facei)
@@ -1191,6 +1203,11 @@ void Foam::polyTopoChange::compact
 
             reorder(localCellMap, cellZone_);
             cellZone_.setCapacity(newCelli);
+            if (cellAdditionalZones_.size())
+            {
+                reorder(localCellMap, cellAdditionalZones_);
+                cellAdditionalZones_.setCapacity(newCelli);
+            }
 
             renumberKey(localCellMap, cellFromPoint_);
             renumberKey(localCellMap, cellFromEdge_);
@@ -1224,6 +1241,14 @@ void Foam::polyTopoChange::compact
                             std::swap(faceOwner_[facei], faceNeighbour_[facei]);
                             flipFaceFlux_.flip(facei);
                             faceZoneFlip_.flip(facei);
+                            if (facei < faceAdditionalZones_.size())
+                            {
+                                for (auto& zas : faceAdditionalZones_[facei])
+                                {
+                                    // Flip sign
+                                    zas = -zas;
+                                }
+                            }
                         }
                     }
                 }
@@ -1582,6 +1607,19 @@ void Foam::polyTopoChange::resetZones
             }
             nPoints[zonei]++;
         }
+        forAll(pointAdditionalZones_, pointi)
+        {
+            for (const label zonei : pointAdditionalZones_[pointi])
+            {
+                if (zonei < 0 || zonei >= pointZones.size())
+                {
+                    FatalErrorInFunction
+                        << "Illegal zoneID " << zonei << " for point "
+                        << pointi << abort(FatalError);
+                }
+                nPoints[zonei]++;
+            }
+        }
 
         // Distribute points per zone
 
@@ -1599,6 +1637,13 @@ void Foam::polyTopoChange::resetZones
 
             addressing[zonei][nPoints[zonei]++] = pointi;
         }
+        forAll(pointAdditionalZones_, pointi)
+        {
+            for (const label zonei : pointAdditionalZones_[pointi])
+            {
+                addressing[zonei][nPoints[zonei]++] = pointi;
+            }
+        }
         // Sort the addressing
         forAll(addressing, zonei)
         {
@@ -1668,6 +1713,20 @@ void Foam::polyTopoChange::resetZones
             }
             nFaces[zonei]++;
         }
+        forAll(faceAdditionalZones_, facei)
+        {
+            for (const label zoneAndSign : faceAdditionalZones_[facei])
+            {
+                const label zonei = mag(zoneAndSign)-1;
+                if (zonei < 0 || zonei >= faceZones.size())
+                {
+                    FatalErrorInFunction
+                        << "Illegal zoneID " << zonei << " for face "
+                        << facei << abort(FatalError);
+                }
+                nFaces[zonei]++;
+            }
+        }
 
         labelListList addressing(faceZones.size());
         boolListList flipMode(faceZones.size());
@@ -1685,9 +1744,21 @@ void Foam::polyTopoChange::resetZones
             const label zonei = iter.val();
 
             const label index = nFaces[zonei]++;
-
             addressing[zonei][index] = facei;
             flipMode[zonei][index] = faceZoneFlip_.test(facei);
+
+            if (facei < faceAdditionalZones_.size())
+            {
+                for (const label zoneAndSign : faceAdditionalZones_[facei])
+                {
+                    const label zonei = mag(zoneAndSign)-1;
+                    const bool flip = (zoneAndSign > 0);
+
+                    const label index = nFaces[zonei]++;
+                    addressing[zonei][index] = facei;
+                    flipMode[zonei][index] = flip;
+                }
+            }
         }
 
         // Sort the addressing
@@ -1784,6 +1855,20 @@ void Foam::polyTopoChange::resetZones
                 nCells[zonei]++;
             }
         }
+        forAll(cellAdditionalZones_, celli)
+        {
+            for (const label zonei : cellAdditionalZones_[celli])
+            {
+                if (zonei < 0 || zonei >= cellZones.size())
+                {
+                    FatalErrorInFunction
+                        << "Illegal zoneID " << zonei << " for cell "
+                        << celli << abort(FatalError);
+                }
+                nCells[zonei]++;
+            }
+        }
+
 
         labelListList addressing(cellZones.size());
         forAll(addressing, zonei)
@@ -1801,6 +1886,13 @@ void Foam::polyTopoChange::resetZones
                 addressing[zonei][nCells[zonei]++] = celli;
             }
         }
+        forAll(cellAdditionalZones_, celli)
+        {
+            for (const label zonei : cellAdditionalZones_[celli])
+            {
+                addressing[zonei][nCells[zonei]++] = celli;
+            }
+        }
         // Sort the addressing
         forAll(addressing, zonei)
         {
@@ -2236,6 +2328,7 @@ void Foam::polyTopoChange::clear()
     pointMap_.clearStorage();
     reversePointMap_.clearStorage();
     pointZone_.clearStorage();
+    pointAdditionalZones_.clearStorage();
     retiredPoints_.clearStorage();
 
     faces_.clearStorage();
@@ -2249,11 +2342,13 @@ void Foam::polyTopoChange::clear()
     flipFaceFlux_.clearStorage();
     faceZone_.clearStorage();
     faceZoneFlip_.clearStorage();
+    faceAdditionalZones_.clearStorage();
     nActiveFaces_ = 0;
 
     cellMap_.clearStorage();
     reverseCellMap_.clearStorage();
     cellZone_.clearStorage();
+    cellAdditionalZones_.clearStorage();
     cellFromPoint_.clearStorage();
     cellFromEdge_.clearStorage();
     cellFromFace_.clearStorage();
@@ -2288,30 +2383,43 @@ void Foam::polyTopoChange::addMesh
         reversePointMap_.setCapacity(reversePointMap_.size() + points.size());
         pointZone_.resize(pointZone_.size() + points.size()/128);
 
-        // Precalc offset zones
-        labelList newZoneID(points.size(), -1);
-
-        forAll(pointZones, zonei)
-        {
-            const labelList& pointLabels = pointZones[zonei];
-
-            for (const label pointi : pointLabels)
-            {
-                newZoneID[pointi] = pointZoneMap[zonei];
-            }
-        }
-
-        // Add points in mesh order
-        for (label pointi = 0; pointi < mesh.nPoints(); pointi++)
+        // Add points without pointZone
+        labelList newPointID(points.size());
+        forAll(newPointID, pointi)
         {
-            addPoint
+            // Add point from point
+            newPointID[pointi] = addPoint
             (
                 points[pointi],
                 pointi,
-                newZoneID[pointi],
+                 -1,
                 true
             );
         }
+
+        // Modify for pointZones
+        DynamicList<label> zones;
+        forAll(pointZones, zonei)
+        {
+            for (const label pointi : pointZones[zonei])
+            {
+                // Get previous zones
+                this->pointZones(newPointID[pointi], zones);
+                for (auto& zonei : zones)
+                {
+                    zonei = pointZoneMap[zonei];
+                }
+                // Add new zone
+                zones.append(pointZoneMap[zonei]);
+                modifyPoint
+                (
+                    newPointID[pointi],
+                    points_[newPointID[pointi]],
+                    zones,
+                    true
+                );
+            }
+        }
     }
 
     // Add cells
@@ -2332,39 +2440,31 @@ void Foam::polyTopoChange::addMesh
         cellZone_.setCapacity(cellZone_.size() + nAllCells);
 
 
-        // Precalc offset zones
-        labelList newZoneID(nAllCells, -1);
+        // Add cells without cellZone
+        labelList newCellID(nAllCells);
+        for (label celli = 0; celli < nAllCells; celli++)
+        {
+            // Add cell from cell
+            newCellID[celli] = addCell(-1, -1, -1, celli, -1);
+        }
 
+        // Modify for cellZones
+        DynamicList<label> zones;
         forAll(cellZones, zonei)
         {
-            const labelList& cellLabels = cellZones[zonei];
-
-            for (const label celli : cellLabels)
+            for (const label celli : cellZones[zonei])
             {
-                if (newZoneID[celli] != -1)
-                {
-                    WarningInFunction
-                        << "Cell:" << celli
-                        << " centre:" << mesh.cellCentres()[celli]
-                        << " is in two zones:"
-                        << cellZones[newZoneID[celli]].name()
-                        << " and " << cellZones[zonei].name() << endl
-                        << "    This is not supported."
-                        << " Continuing with first zone only." << endl;
-                }
-                else
+                // Get previous zones
+                this->cellZones(newCellID[celli], zones);
+                for (auto& zonei : zones)
                 {
-                    newZoneID[celli] = cellZoneMap[zonei];
+                    zonei = cellZoneMap[zonei];
                 }
+                // Add new zone
+                zones.append(cellZoneMap[zonei]);
+                modifyCell(newCellID[celli], zones);
             }
         }
-
-        // Add cells in mesh order
-        for (label celli = 0; celli < nAllCells; celli++)
-        {
-            // Add cell from cell
-            addCell(-1, -1, -1, celli, newZoneID[celli]);
-        }
     }
 
     // Add faces
@@ -2388,31 +2488,16 @@ void Foam::polyTopoChange::addMesh
         faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/128);
         flipFaceFlux_.setCapacity(faces_.size() + nAllFaces);
         faceZone_.resize(faceZone_.size() + nAllFaces/128);
-        faceZoneFlip_.setCapacity(faces_.size() + nAllFaces);
-
-
-        // Precalc offset zones
-        labelList newZoneID(nAllFaces, -1);
-        boolList zoneFlip(nAllFaces, false);
-
-        forAll(faceZones, zonei)
-        {
-            const labelList& faceLabels = faceZones[zonei];
-            const boolList& flipMap = faceZones[zonei].flipMap();
+        faceZoneFlip_.setCapacity(faceMap_.capacity());
 
-            forAll(faceLabels, facei)
-            {
-                newZoneID[faceLabels[facei]] = faceZoneMap[zonei];
-                zoneFlip[faceLabels[facei]] = flipMap[facei];
-            }
-        }
 
-        // Add faces in mesh order
+        // Add faces (in mesh order) without faceZone
+        labelList newFaceID(nAllFaces);
 
         // 1. Internal faces
         for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
         {
-            addFace
+            newFaceID[facei] = addFace
             (
                 faces[facei],
                 faceOwner[facei],
@@ -2422,8 +2507,8 @@ void Foam::polyTopoChange::addMesh
                 facei,                      // masterFaceID
                 false,                      // flipFaceFlux
                 -1,                         // patchID
-                newZoneID[facei],           // zoneID
-                zoneFlip[facei]             // zoneFlip
+                -1,                         // zoneID
+                false                       // zoneFlip
             );
         }
 
@@ -2446,7 +2531,7 @@ void Foam::polyTopoChange::addMesh
             {
                 const label facei = pp.start() + patchFacei;
 
-                addFace
+                newFaceID[facei] = addFace
                 (
                     faces[facei],
                     faceOwner[facei],
@@ -2456,8 +2541,49 @@ void Foam::polyTopoChange::addMesh
                     facei,                      // masterFaceID
                     false,                      // flipFaceFlux
                     patchMap[patchi],           // patchID
-                    newZoneID[facei],           // zoneID
-                    zoneFlip[facei]             // zoneFlip
+                    -1,                         // zoneID
+                    false                       // zoneFlip
+                );
+            }
+        }
+
+
+        // Modify for faceZones
+
+        DynamicList<label> zones;
+        DynamicList<bool> flips;
+        forAll(faceZones, zonei)
+        {
+            const auto& fz = faceZones[zonei];
+
+            forAll(fz, i)
+            {
+                const label facei = fz[i];
+                const bool flip = fz.flipMap()[i];
+
+                // Modify face index for combined zones
+                const label newFacei = newFaceID[facei];
+
+                // Get previous zones
+                this->faceZones(newFacei, zones, flips);
+                for (auto& zonei : zones)
+                {
+                    zonei = faceZoneMap[zonei];
+                }
+                // Add new zone
+                zones.append(faceZoneMap[zonei]);
+                flips.append(flip);
+
+                modifyFace
+                (
+                    faces_[newFacei],
+                    newFacei,
+                    faceOwner_[newFacei],
+                    faceNeighbour_[newFacei],
+                    flipFaceFlux_(newFacei),
+                    region_[newFacei],
+                    zones,
+                    flips
                 );
             }
         }
@@ -2644,6 +2770,57 @@ Foam::label Foam::polyTopoChange::addPoint
     {
         pointZone_.insert(pointi, zoneID);
     }
+    // Delay extending the pointAdditionalZones
+
+    if (!inCell)
+    {
+        retiredPoints_.insert(pointi);
+    }
+
+    return pointi;
+}
+
+
+Foam::label Foam::polyTopoChange::addPoint
+(
+    const point& pt,
+    const label masterPointID,
+    const labelUList& zoneIDs,
+    const bool inCell
+)
+{
+    const label pointi = points_.size();
+
+    points_.append(pt);
+    pointMap_.append(masterPointID);
+    reversePointMap_.append(pointi);
+
+    if (zoneIDs.size())
+    {
+        const label minIndex = findMin(zoneIDs);
+
+        pointZone_.set(pointi, zoneIDs[minIndex]);
+
+        // Clear any additional storage
+        if (pointi < pointAdditionalZones_.size())
+        {
+            pointAdditionalZones_[pointi].clear();
+        }
+        // (auto-vivify and) store additional zones
+        forAll(zoneIDs, i)
+        {
+            if (i != minIndex)
+            {
+                if (zoneIDs[i] == zoneIDs[minIndex])
+                {
+                    FatalErrorInFunction << "Duplicates in zones "
+                        << flatOutput(zoneIDs) << " for point " << pointi
+                        << exit(FatalError);
+                }
+                pointAdditionalZones_(pointi).push_uniq(zoneIDs[i]);
+            }
+        }
+    }
 
     if (!inCell)
     {
@@ -2659,6 +2836,88 @@ void Foam::polyTopoChange::modifyPoint
     const label pointi,
     const point& pt,
     const label zoneID,
+    const bool inCell,
+    const bool multiZone
+)
+{
+    if (pointi < 0 || pointi >= points_.size())
+    {
+        FatalErrorInFunction
+            << "illegal point label " << pointi << endl
+            << "Valid point labels are 0 .. " << points_.size()-1
+            << abort(FatalError);
+    }
+    if (pointRemoved(pointi) || pointMap_[pointi] == -1)
+    {
+        FatalErrorInFunction
+            << "point " << pointi << " already marked for removal"
+            << abort(FatalError);
+    }
+    points_[pointi] = pt;
+
+    if (!multiZone)
+    {
+        if (zoneID >= 0)
+        {
+            pointZone_.set(pointi, zoneID);
+        }
+        else
+        {
+            pointZone_.erase(pointi);
+        }
+        if (pointi < pointAdditionalZones_.size())
+        {
+            pointAdditionalZones_[pointi].clear();
+        }
+    }
+    else
+    {
+        auto fnd = pointZone_.find(pointi);
+        if (!fnd)
+        {
+            pointZone_.set(pointi, zoneID);
+
+            // Check that no leftover additional zones
+            if (pointAdditionalZones_(pointi).size())
+            {
+                FatalErrorInFunction
+                    << "Additional zones for point:"
+                    << pointAdditionalZones_(pointi)
+                    << abort(FatalError);
+            }
+        }
+        else if (zoneID == -1)
+        {
+            // Clear
+            pointZone_.erase(fnd);
+            if (pointi < pointAdditionalZones_.size())
+            {
+                pointAdditionalZones_[pointi].clear();
+            }
+        }
+        else if (zoneID != fnd())
+        {
+            // New zone. Store in additional storage.
+            pointAdditionalZones_(pointi).push_uniq(zoneID);
+        }
+    }
+
+    if (inCell)
+    {
+        retiredPoints_.erase(pointi);
+    }
+    else
+    {
+        retiredPoints_.insert(pointi);
+    }
+}
+
+
+void Foam::polyTopoChange::modifyPoint
+(
+    const label pointi,
+    const point& pt,
+    const labelUList& zoneIDs,
     const bool inCell
 )
 {
@@ -2677,13 +2936,33 @@ void Foam::polyTopoChange::modifyPoint
     }
     points_[pointi] = pt;
 
-    if (zoneID >= 0)
+    if (zoneIDs.size())
     {
-        pointZone_.set(pointi, zoneID);
+        if (zoneIDs.found(-1))
+        {
+            FatalErrorInFunction
+                << "zones cannot contain -1 : " << flatOutput(zoneIDs)
+                << " for point:" << pointi
+                << abort(FatalError);
+        }
+
+        pointZone_.set(pointi, zoneIDs[0]);
+        if (pointi < pointAdditionalZones_.size())
+        {
+            pointAdditionalZones_[pointi].clear();
+        }
+        for (label i = 1; i < zoneIDs.size(); ++i)
+        {
+            pointAdditionalZones_(pointi).push_uniq(zoneIDs[i]);
+        }
     }
     else
     {
         pointZone_.erase(pointi);
+        if (pointi < pointAdditionalZones_.size())
+        {
+            pointAdditionalZones_[pointi].clear();
+        }
     }
 
     if (inCell)
@@ -2696,7 +2975,6 @@ void Foam::polyTopoChange::modifyPoint
     }
 }
 
-
 void Foam::polyTopoChange::movePoints(const pointField& newPoints)
 {
     if (newPoints.size() != points_.size())
@@ -2759,6 +3037,10 @@ void Foam::polyTopoChange::removePoint
         reversePointMap_[pointi] = -1;
     }
     pointZone_.erase(pointi);
+    if (pointi < pointAdditionalZones_.size())
+    {
+        pointAdditionalZones_[pointi].clear();
+    }
     retiredPoints_.erase(pointi);
 }
 
@@ -2822,6 +3104,100 @@ Foam::label Foam::polyTopoChange::addFace
         faceZone_.insert(facei, zoneID);
     }
     faceZoneFlip_.set(facei, zoneFlip);
+    // Delay extending the faceAdditionalZones
+
+    return facei;
+}
+
+
+Foam::label Foam::polyTopoChange::addFace
+(
+    const face& f,
+    const label own,
+    const label nei,
+    const label masterPointID,
+    const label masterEdgeID,
+    const label masterFaceID,
+    const bool flipFaceFlux,
+    const label patchID,
+    const labelUList& zoneIDs,
+    const UList<bool>& zoneFlips
+)
+{
+    // Check validity
+    if (debug)
+    {
+        // Note: no check on zones
+        checkFace(f, -1, own, nei, patchID, -1);
+    }
+
+    label facei = faces_.size();
+
+    faces_.append(f);
+    region_.append(patchID);
+    faceOwner_.append(own);
+    faceNeighbour_.append(nei);
+
+    if (masterPointID >= 0)
+    {
+        faceMap_.append(-1);
+        faceFromPoint_.insert(facei, masterPointID);
+    }
+    else if (masterEdgeID >= 0)
+    {
+        faceMap_.append(-1);
+        faceFromEdge_.insert(facei, masterEdgeID);
+    }
+    else if (masterFaceID >= 0)
+    {
+        faceMap_.append(masterFaceID);
+    }
+    else
+    {
+        // Allow inflate-from-nothing?
+        //FatalErrorInFunction
+        //    << "Need to specify a master point, edge or face"
+        //    << "face:" << f << " own:" << own << " nei:" << nei
+        //    << abort(FatalError);
+        faceMap_.append(-1);
+    }
+    reverseFaceMap_.append(facei);
+
+    flipFaceFlux_.set(facei, flipFaceFlux);
+
+    if (zoneIDs.size())
+    {
+        const label minIndex = findMin(zoneIDs);
+
+        faceZone_.set(facei, zoneIDs[minIndex]);
+        faceZoneFlip_.set(facei, zoneFlips[minIndex]);
+
+        // Clear any additional storage
+        if (facei < faceAdditionalZones_.size())
+        {
+            faceAdditionalZones_[facei].clear();
+        }
+        // (auto-vivify and) store additional zones
+        forAll(zoneIDs, i)
+        {
+            if (i != minIndex)
+            {
+                if (zoneIDs[i] == zoneIDs[minIndex])
+                {
+                    FatalErrorInFunction << "Duplicates in zones "
+                        << flatOutput(zoneIDs) << " for face " << facei
+                        << exit(FatalError);
+                }
+                const label zoneAndSign
+                (
+                    zoneFlips[i]
+                  ? zoneIDs[i]+1
+                  : -(zoneIDs[i]+1)
+                );
+                faceAdditionalZones_(facei).push_uniq(zoneAndSign);
+            }
+        }
+    }
 
     return facei;
 }
@@ -2836,7 +3212,8 @@ void Foam::polyTopoChange::modifyFace
     const bool flipFaceFlux,
     const label patchID,
     const label zoneID,
-    const bool zoneFlip
+    const bool zoneFlip,
+    const bool multiZone
 )
 {
     // Check validity
@@ -2851,15 +3228,133 @@ void Foam::polyTopoChange::modifyFace
     region_[facei] = patchID;
 
     flipFaceFlux_.set(facei, flipFaceFlux);
-    faceZoneFlip_.set(facei, zoneFlip);
 
-    if (zoneID >= 0)
+    // Note: same logic as modifyCell except storage is now sparse instead
+    //       of marked with -1
+
+    if (!multiZone)
+    {
+        if (zoneID == -1)
+        {
+            faceZone_.erase(facei);
+            faceZoneFlip_.set(facei, zoneFlip);
+        }
+        else
+        {
+            faceZone_.set(facei, zoneID);
+            faceZoneFlip_.set(facei, zoneFlip);
+        }
+        if (facei < faceAdditionalZones_.size())
+        {
+            faceAdditionalZones_[facei].clear();
+        }
+    }
+    else
+    {
+        auto fnd = faceZone_.find(facei);
+        if (!fnd)
+        {
+            faceZone_.set(facei, zoneID);
+            faceZoneFlip_.set(facei, zoneFlip);
+
+            // Check that no leftover additional zones
+            if (faceAdditionalZones_(facei).size())
+            {
+                FatalErrorInFunction
+                    << "Additional zones for face:"
+                    << faceAdditionalZones_(facei)
+                    << abort(FatalError);
+            }
+        }
+        else if (zoneID == -1)
+        {
+            // Clear
+            faceZone_.erase(fnd);
+            faceZoneFlip_.set(facei, false);
+            if (facei < faceAdditionalZones_.size())
+            {
+                faceAdditionalZones_[facei].clear();
+            }
+        }
+        else if (zoneID != fnd())
+        {
+            // New zone. Store in additional storage.
+            const label zoneAndSign
+            (
+                zoneFlip
+              ? zoneID+1
+              : -(zoneID+1)
+            );
+            faceAdditionalZones_(facei).push_uniq(zoneAndSign);
+        }
+    }
+}
+
+
+void Foam::polyTopoChange::modifyFace
+(
+    const face& f,
+    const label facei,
+    const label own,
+    const label nei,
+    const bool flipFaceFlux,
+    const label patchID,
+    const labelUList& zoneIDs,
+    const UList<bool>& zoneFlips
+)
+{
+    // Check validity
+    if (debug)
+    {
+        // Note: no check on zones
+        checkFace(f, facei, own, nei, patchID, -1);
+    }
+
+    faces_[facei] = f;
+    faceOwner_[facei] = own;
+    faceNeighbour_[facei] = nei;
+    region_[facei] = patchID;
+
+    flipFaceFlux_.set(facei, flipFaceFlux);
+
+    // Note: same logic as modifyCell except storage is now sparse instead
+    //       of marked with -1
+
+    if (zoneIDs.size())
     {
-        faceZone_.set(facei, zoneID);
+        if (zoneIDs.found(-1))
+        {
+            FatalErrorInFunction
+                << "zones cannot contain -1 : " << flatOutput(zoneIDs)
+                << " for face:" << facei
+                << abort(FatalError);
+        }
+
+        faceZone_.set(facei, zoneIDs[0]);
+        faceZoneFlip_.set(facei, zoneFlips[0]);
+        if (facei < faceAdditionalZones_.size())
+        {
+            faceAdditionalZones_[facei].clear();
+        }
+        for (label i = 1; i < zoneIDs.size(); ++i)
+        {
+            const label zoneAndSign
+            (
+                zoneFlips[i]
+              ? zoneIDs[i]+1
+              : -(zoneIDs[i]+1)
+            );
+            faceAdditionalZones_(facei).push_uniq(zoneAndSign);
+        }
     }
     else
     {
         faceZone_.erase(facei);
+        faceZoneFlip_.set(facei, false);
+        if (facei < faceAdditionalZones_.size())
+        {
+            faceAdditionalZones_[facei].clear();
+        }
     }
 }
 
@@ -2908,6 +3403,10 @@ void Foam::polyTopoChange::removeFace
     flipFaceFlux_.unset(facei);
     faceZoneFlip_.unset(facei);
     faceZone_.erase(facei);
+    if (facei < faceAdditionalZones_.size())
+    {
+        faceAdditionalZones_[facei].clear();
+    }
 }
 
 
@@ -2942,7 +3441,65 @@ Foam::label Foam::polyTopoChange::addCell
         cellMap_.append(masterCellID);
     }
     reverseCellMap_.append(celli);
+
     cellZone_.append(zoneID);
+    // Delay extending the cellAdditionalZones
+
+    return celli;
+}
+
+
+Foam::label Foam::polyTopoChange::addCell
+(
+    const label masterPointID,
+    const label masterEdgeID,
+    const label masterFaceID,
+    const label masterCellID,
+    const labelUList& zoneIDs
+)
+{
+    label celli = cellMap_.size();
+
+    if (masterPointID >= 0)
+    {
+        cellMap_.append(-1);
+        cellFromPoint_.insert(celli, masterPointID);
+    }
+    else if (masterEdgeID >= 0)
+    {
+        cellMap_.append(-1);
+        cellFromEdge_.insert(celli, masterEdgeID);
+    }
+    else if (masterFaceID >= 0)
+    {
+        cellMap_.append(-1);
+        cellFromFace_.insert(celli, masterFaceID);
+    }
+    else
+    {
+        cellMap_.append(masterCellID);
+    }
+    reverseCellMap_.append(celli);
+
+    if (zoneIDs.size())
+    {
+        cellZone_.append(zoneIDs[0]);
+
+        // Clear any additional storage
+        if (celli < cellAdditionalZones_.size())
+        {
+            cellAdditionalZones_[celli].clear();
+        }
+        // (auto-vivify and) store additional zones
+        for (label i = 1; i < zoneIDs.size(); ++i)
+        {
+            cellAdditionalZones_(celli).push_uniq(zoneIDs[i]);
+        }
+    }
+    else
+    {
+        cellZone_.append(-1);
+    }
 
     return celli;
 }
@@ -2951,10 +3508,96 @@ Foam::label Foam::polyTopoChange::addCell
 void Foam::polyTopoChange::modifyCell
 (
     const label celli,
-    const label zoneID
+    const label zoneID,
+    const bool multiZone
+)
+{
+    if (!multiZone)
+    {
+        cellZone_[celli] = zoneID;
+        if (celli < cellAdditionalZones_.size())
+        {
+            cellAdditionalZones_[celli].clear();
+        }
+    }
+    else
+    {
+        if (cellZone_[celli] == -1)
+        {
+            cellZone_[celli] = zoneID;
+
+            // Check that no leftover additional zones
+            if (cellAdditionalZones_(celli).size())
+            {
+                FatalErrorInFunction
+                    << "Additional zones for cell:"
+                    << cellAdditionalZones_(celli)
+                    << abort(FatalError);
+            }
+        }
+        else
+        {
+            if (zoneID == -1)
+            {
+                // Clear
+                cellZone_[celli] = zoneID;
+                if (celli < cellAdditionalZones_.size())
+                {
+                    cellAdditionalZones_[celli].clear();
+                }
+            }
+            else if (zoneID != cellZone_[celli])
+            {
+                // New zone. Store in additional storage.
+                cellAdditionalZones_(celli).push_uniq(zoneID);
+            }
+        }
+    }
+}
+
+
+void Foam::polyTopoChange::modifyCell
+(
+    const label celli,
+    const labelUList& zoneIDs
 )
 {
-    cellZone_[celli] = zoneID;
+    if (celli < 0 || celli >= cellMap_.size())
+    {
+        FatalErrorInFunction
+            << "illegal cell label " << celli << endl
+            << "Valid cell labels are 0 .. " << cellMap_.size()-1
+            << abort(FatalError);
+    }
+
+    if (zoneIDs.size())
+    {
+        if (zoneIDs.found(-1))
+        {
+            FatalErrorInFunction
+                << "zones cannot contain -1 : " << flatOutput(zoneIDs)
+                << " for cell:" << celli
+                << abort(FatalError);
+        }
+
+        cellZone_[celli] = zoneIDs[0];
+        if (celli < cellAdditionalZones_.size())
+        {
+            cellAdditionalZones_[celli].clear();
+        }
+        for (label i = 1; i < zoneIDs.size(); ++i)
+        {
+            cellAdditionalZones_(celli).push_uniq(zoneIDs[i]);
+        }
+    }
+    else
+    {
+        cellZone_[celli] = -1;
+        if (celli < cellAdditionalZones_.size())
+        {
+            cellAdditionalZones_[celli].clear();
+        }
+    }
 }
 
 
@@ -2993,6 +3636,109 @@ void Foam::polyTopoChange::removeCell
     cellFromEdge_.erase(celli);
     cellFromFace_.erase(celli);
     cellZone_[celli] = -1;
+    if (celli < cellAdditionalZones_.size())
+    {
+        cellAdditionalZones_[celli].clear();
+    }
+}
+
+
+Foam::label Foam::polyTopoChange::pointZones
+(
+    const label pointi,
+    DynamicList<label>& zones
+) const
+{
+    if (pointi < 0 || pointi >= pointMap_.size())
+    {
+        FatalErrorInFunction
+            << "illegal point label " << pointi << endl
+            << "Valid point labels are 0 .. " << pointMap_.size()-1
+            << abort(FatalError);
+    }
+    zones.clear();
+
+    const auto fnd = pointZone_.find(pointi);
+    if (fnd)
+    {
+        zones.push_back(fnd());
+        if (pointi < pointAdditionalZones_.size())
+        {
+            for (const label zonei : pointAdditionalZones_[pointi])
+            {
+                zones.push_uniq(zonei);
+            }
+        }
+    }
+    return zones.size();
+}
+
+
+Foam::label Foam::polyTopoChange::faceZones
+(
+    const label facei,
+    DynamicList<label>& zones,
+    DynamicList<bool>& flips
+) const
+{
+    if (facei < 0 || facei >= faceMap_.size())
+    {
+        FatalErrorInFunction
+            << "illegal face label " << facei << endl
+            << "Valid face labels are 0 .. " << faceMap_.size()-1
+            << abort(FatalError);
+    }
+    zones.clear();
+    flips.clear();
+
+    const auto fnd = faceZone_.find(facei);
+    if (fnd)
+    {
+        zones.push_back(fnd());
+        flips.push_back(flipFaceFlux_[facei]);
+        if (facei < faceAdditionalZones_.size())
+        {
+            for (const label zoneAndSign : faceAdditionalZones_[facei])
+            {
+                const label zonei = mag(zoneAndSign)-1;
+                if (!zones.found(zonei))
+                {
+                    zones.push_back(zonei);
+                    flips.push_back((zoneAndSign > 0));
+                }
+            }
+        }
+    }
+    return zones.size();
+}
+
+
+Foam::label Foam::polyTopoChange::cellZones
+(
+    const label celli,
+    DynamicList<label>& zones
+) const
+{
+    if (celli < 0 || celli >= cellMap_.size())
+    {
+        FatalErrorInFunction
+            << "illegal cell label " << celli << endl
+            << "Valid cell labels are 0 .. " << cellMap_.size()-1
+            << abort(FatalError);
+    }
+    zones.clear();
+    if (cellZone_[celli] != -1)
+    {
+        zones.push_back(cellZone_[celli]);
+    }
+    if (celli < cellAdditionalZones_.size())
+    {
+        for (const label zonei : cellAdditionalZones_[celli])
+        {
+            zones.push_uniq(zonei);
+        }
+    }
+    return zones.size();
 }
 
 
diff --git a/src/meshTools/polyTopoChange/polyTopoChange.H b/src/meshTools/polyTopoChange/polyTopoChange.H
index 9ad146d575b..318142f9da7 100644
--- a/src/meshTools/polyTopoChange/polyTopoChange.H
+++ b/src/meshTools/polyTopoChange/polyTopoChange.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2022,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,6 +55,14 @@ Description
     - coupled patches: the reorderCoupledFaces routine (borrowed from
     the couplePatches utility) reorders coupled patch faces and
     uses the cyclicPolyPatch,processorPolyPatch functionality.
+    - zones are assumed to be non-overlapping by default. If desired to be
+    overlapping either set the multiZone to true when calling
+    modifyCell|Face|Point or use the variants of addCell|Face|Point
+    and modifyCell|Face|Point that take a list of zones.
+    - if overlapping zones:
+        - 'main' zone is the lowest numbered zone. -1 means no zones.
+        - 'additional' zones are stored in incremental ordering (and cannot
+          contain -1)
 
 SourceFiles
     polyTopoChange.C
@@ -127,6 +135,9 @@ class polyTopoChange
             //- Zone of point
             Map<label> pointZone_;
 
+            //- Additional zones of points
+            DynamicList<labelList> pointAdditionalZones_;
+
             //- Retired points
             labelHashSet retiredPoints_;
 
@@ -170,6 +181,9 @@ class polyTopoChange
             //- Orientation of face in zone
             bitSet faceZoneFlip_;
 
+            //- Additional zones of face (stored as signed face)
+            DynamicList<labelList> faceAdditionalZones_;
+
             //- Active faces
             label nActiveFaces_;
 
@@ -193,9 +207,12 @@ class polyTopoChange
             //- Cells added from face
             Map<label> cellFromFace_;
 
-            //- Zone of cell
+            //- First zone of cell (or -1)
             DynamicList<label> cellZone_;
 
+            //- Additional zones of cell
+            DynamicList<labelList> cellAdditionalZones_;
+
 
     // Private Member Functions
 
@@ -207,13 +224,6 @@ class polyTopoChange
             DynamicList<Type>& lst
         );
 
-        template<class Type>
-        static void reorder
-        (
-            const labelUList& oldToNew,
-            List<DynamicList<Type>>& lst
-        );
-
         template<class Type>
         static void renumberKey
         (
@@ -533,6 +543,18 @@ public:
                 const bool inCell
             );
 
+            //- Add point. Return new point label.
+            //  Notes:
+            //  - masterPointID can be < 0 (appended points)
+            //  - inCell = false: add retired point (to end of point list)
+            label addPoint
+            (
+                const point& pt,
+                const label masterPointID,
+                const labelUList& zoneIDs,
+                const bool inCell
+            );
+
             //- Modify coordinate.
             //  Notes:
             //  - zoneID = +ve (add to zoneID), -ve (remove from zones)
@@ -542,12 +564,32 @@ public:
                 const label pointi,
                 const point& pt,
                 const label zoneID,
+                const bool inCell,
+                const bool multiZone = false
+            );
+
+            //- Modify coordinate.
+            //  Notes:
+            //  - zoneIDs = set pointZones
+            //  - inCell = false: add retired point (to end of point list)
+            void modifyPoint
+            (
+                const label pointi,
+                const point& pt,
+                const labelUList& zoneIDs,
                 const bool inCell
             );
 
             //- Remove/merge point.
             void removePoint(const label pointi, const label mergePointi);
 
+            //- Get current cellZone(s). Return number of zones.
+            label pointZones
+            (
+                const label pointi,
+                DynamicList<label>& zones
+            ) const;
+
             //- Add face to cells. Return new face label.
             //  own,nei<0, zoneID>=0 : add inactive face (to end of face list)
             label addFace
@@ -564,6 +606,22 @@ public:
                 const bool zoneFlip
             );
 
+            //- Add face to cells. Return new face label.
+            //  own,nei<0, zoneID>=0 : add inactive face (to end of face list)
+            label addFace
+            (
+                const face& f,
+                const label own,
+                const label nei,
+                const label masterPointID,
+                const label masterEdgeID,
+                const label masterFaceID,
+                const bool flipFaceFlux,
+                const label patchID,
+                const labelUList& zoneIDs,
+                const UList<bool>& zoneFlips
+            );
+
             //- Modify vertices or cell of face.
             void modifyFace
             (
@@ -574,12 +632,34 @@ public:
                 const bool flipFaceFlux,
                 const label patchID,
                 const label zoneID,
-                const bool zoneFlip
+                const bool zoneFlip,
+                const bool multiZone = false
+            );
+
+            //- Modify vertices or cell of face.
+            void modifyFace
+            (
+                const face& f,
+                const label facei,
+                const label own,
+                const label nei,
+                const bool flipFaceFlux,
+                const label patchID,
+                const labelUList& zoneIDs,
+                const UList<bool>& zoneFlips
             );
 
             //- Remove/merge face.
             void removeFace(const label facei, const label mergeFacei);
 
+            //- Get current faceZone(s). Return number of zones.
+            label faceZones
+            (
+                const label facei,
+                DynamicList<label>& zones,
+                DynamicList<bool>& flips
+            ) const;
+
             //- Add cell. Return new cell label.
             label addCell
             (
@@ -590,12 +670,34 @@ public:
                 const label zoneID
             );
 
-            //- Modify zone of cell
-            void modifyCell(const label celli, const label zoneID);
+            //- Add zoned cell (zones cannot be -1). Return new cell label.
+            label addCell
+            (
+                const label masterPointID,
+                const label masterEdgeID,
+                const label masterFaceID,
+                const label masterCellID,
+                const labelUList& zoneIDs
+            );
+
+            //- Modify zone of cell. Optionally add to zone.
+            void modifyCell
+            (
+                const label celli,
+                const label zoneID,
+                const bool multiZone = false
+            );
+
+            //- Set zones of cell
+            void modifyCell(const label celli, const labelUList& zoneIDs);
 
             //- Remove/merge cell.
             void removeCell(const label celli, const label mergeCelli);
 
+            //- Get current cellZone(s). Return number of zones.
+            label cellZones(const label celli, DynamicList<label>& zones) const;
+
+
             //- Explicitly set the number of patches if construct-without-mesh
             //- used.
             inline void setNumPatches(const label nPatches);
diff --git a/src/meshTools/polyTopoChange/polyTopoChangeTemplates.C b/src/meshTools/polyTopoChange/polyTopoChangeTemplates.C
index 9cc1362b439..9016ce6a3a9 100644
--- a/src/meshTools/polyTopoChange/polyTopoChangeTemplates.C
+++ b/src/meshTools/polyTopoChange/polyTopoChangeTemplates.C
@@ -43,35 +43,13 @@ void Foam::polyTopoChange::reorder
     // Create copy
     DynamicList<Type> oldLst(lst);
 
-    forAll(oldToNew, i)
+    forAll(oldLst, i)
     {
         const label newIdx = oldToNew[i];
 
         if (newIdx >= 0)
         {
-            lst[newIdx] = oldLst[i];
-        }
-    }
-}
-
-
-template<class Type>
-void Foam::polyTopoChange::reorder
-(
-    const labelUList& oldToNew,
-    List<DynamicList<Type>>& lst
-)
-{
-    // Create copy
-    List<DynamicList<Type>> oldLst(lst);
-
-    forAll(oldToNew, i)
-    {
-        const label newIdx = oldToNew[i];
-
-        if (newIdx >= 0)
-        {
-            lst[newIdx].transfer(oldLst[i]);
+            lst[newIdx] = std::move(oldLst[i]);
         }
     }
 }
diff --git a/tutorials/mesh/parallel/filter/system/fvOptions b/tutorials/mesh/parallel/filter/system/fvOptions
index b4dec7e9c38..60b12b0abdb 100644
--- a/tutorials/mesh/parallel/filter/system/fvOptions
+++ b/tutorials/mesh/parallel/filter/system/fvOptions
@@ -42,6 +42,29 @@ filter1
 }
 
 
+heater1
+{
+    // Add some additional water inside a cellZone that is nested inside
+    // the filter. This is just to demonstrate overlapping cellZones
+
+    type            scalarSemiImplicitSource;
+    active          yes;
+
+    timeStart       0.0;
+    duration        0.4;
+    volumeMode      specific;
+
+    selectionMode   cellZone;
+    cellZone        heater;
+
+    sources
+    {
+        rho         (1e-2 0); // kg/s/m^3
+        H2O         (1e-2 0); // kg/s/m^3
+    }
+}
+
+
 massSource1
 {
     type            scalarSemiImplicitSource;
diff --git a/tutorials/mesh/parallel/filter/system/topoSetDict b/tutorials/mesh/parallel/filter/system/topoSetDict
index 1229a203daa..722ff8df865 100644
--- a/tutorials/mesh/parallel/filter/system/topoSetDict
+++ b/tutorials/mesh/parallel/filter/system/topoSetDict
@@ -32,7 +32,6 @@ actions
         set     filterCells;
     }
 
-
     {
         name    leftFluidCells;
         type    cellSet;
@@ -114,6 +113,41 @@ actions
         source  setToFaceZone;
         faceSet cycRightFaces;   // name of faceSet
     }
+
+
+    // Additional zones inside filter
+
+    // heater inside filter
+    {
+        name    heaterCells;
+        type    cellSet;
+        action  new;
+        source  boxToCell;
+        box     (1.7 -10 -10) (1.8 10 10);
+    }
+    // Note: could be done directly using boxToCell
+    {
+        name    heater;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        set     heaterCells;
+    }
+    {
+        name    heaterOutsideFaces;
+        type    faceSet;
+        action  new;
+        source  cellToFace;
+        set     heaterCells;
+        option  outside;
+    }
+    {
+        name    heaterOutside;
+        type    faceZoneSet;
+        action  new;
+        source  setToFaceZone;
+        faceSet heaterOutsideFaces;   // name of faceSet
+    }
 );
 
 
-- 
GitLab