From 5677e10d907a98c103eb0cb0ca2334c20586d565 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 3 Oct 2022 11:41:42 +0100
Subject: [PATCH] ENH: checkMesh: report overlapping zones. See #2521

---
 .../manipulation/checkMesh/checkTopology.C    | 131 +++++++++++++++++-
 .../manipulation/checkMesh/checkTopology.H    |  10 ++
 2 files changed, 137 insertions(+), 4 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index 4bdb0cf1b1..82af0aa152 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -104,9 +104,6 @@ void Foam::checkPatch
             //    << endl;
         }
 
-        //DebugVar(globalEdgeFaces);
-
-
         // Synchronise across coupled edges.
         syncTools::syncEdgeList
         (
@@ -116,7 +113,6 @@ void Foam::checkPatch
             labelList()             // null value
         );
 
-        //DebugVar(globalEdgeFaces);
         
         label labelTyp = TopoType::MANIFOLD;
         forAll(meshEdges, edgei)
@@ -200,6 +196,35 @@ void Foam::checkPatch
 }
 
 
+template<class Zone>
+Foam::label Foam::checkZones
+(
+    const polyMesh& mesh,
+    const ZoneMesh<Zone, polyMesh>& zones,
+    topoSet& set
+)
+{
+    labelList zoneID(set.maxSize(mesh), -1);
+    for (const auto& zone : zones)
+    {
+        for (const label elem : zone)
+        {
+            if
+            (
+                zoneID[elem] != -1
+             && zoneID[elem] != zone.index()
+            )
+            {
+                set.insert(elem);
+            }
+            zoneID[elem] = zone.index();
+        }
+    }
+
+    return returnReduce(set.size(), sumOp<label>());
+}
+
+
 Foam::label Foam::checkTopology
 (
     const polyMesh& mesh,
@@ -714,6 +739,25 @@ Foam::label Foam::checkTopology
                 );
                 Info<< endl;
             }
+
+            // Check for duplicates
+            if (allTopology)
+            {
+                faceSet mzFaces(mesh, "multiZoneFaces", mesh.nFaces()/100);
+                const label nMulti = checkZones(mesh, faceZones, mzFaces);
+                if (nMulti)
+                {
+                    Info<< "  <<Writing " << nMulti
+                        << " faces that are in multiple zones"
+                        << " to set " << mzFaces.name() << endl;
+                    mzFaces.instance() = mesh.pointsInstance();
+                    mzFaces.write();
+                    if (surfWriter && surfWriter->enabled())
+                    {
+                        mergeAndWrite(*surfWriter, mzFaces);
+                    }
+                }
+            }
         }
         else
         {
@@ -791,6 +835,26 @@ Foam::label Foam::checkTopology
                     << returnReduce(v, sumOp<scalar>())
                     << ' ' << bb << endl;
             }
+
+
+            // Check for duplicates
+            if (allTopology)
+            {
+                cellSet mzCells(mesh, "multiZoneCells", mesh.nCells()/100);
+                const label nMulti = checkZones(mesh, cellZones, mzCells);
+                if (nMulti)
+                {
+                    Info<< "  <<Writing " << nMulti
+                        << " cells that are in multiple zones"
+                        << " to set " << mzCells.name() << endl;
+                    mzCells.instance() = mesh.pointsInstance();
+                    mzCells.write();
+                    if (surfWriter && surfWriter->enabled())
+                    {
+                        mergeAndWrite(*surfWriter, mzCells);
+                    }
+                }
+            }
         }
         else
         {
@@ -798,6 +862,65 @@ Foam::label Foam::checkTopology
         }
     }
 
+
+    {
+        Info<< "\nChecking basic pointZone addressing..." << endl;
+
+        Pout.setf(ios_base::left);
+
+        const pointZoneMesh& pointZones = mesh.pointZones();
+
+        if (pointZones.size())
+        {
+            Info<< "    "
+                << setw(20) << "PointZone"
+                << setw(8) << "Points"
+                << "BoundingBox" << endl;
+
+            for (const auto& zone : pointZones)
+            {
+                boundBox bb;
+                for (const label pointi : zone)
+                {
+                    bb.add(mesh.points()[pointi]);
+                }
+
+                bb.reduce();  // Global min/max
+
+                Info<< "    "
+                    << setw(20) << zone.name()
+                    << setw(8)
+                    << returnReduce(zone.size(), sumOp<label>())
+                    << bb << endl;
+            }
+
+
+            // Check for duplicates
+            if (allTopology)
+            {
+                pointSet mzPoints(mesh, "multiZonePoints", mesh.nPoints()/100);
+                const label nMulti = checkZones(mesh, pointZones, mzPoints);
+                if (nMulti)
+                {
+                    Info<< "  <<Writing " << nMulti
+                        << " points that are in multiple zones"
+                        << " to set " << mzPoints.name() << endl;
+                    mzPoints.instance() = mesh.pointsInstance();
+                    mzPoints.write();
+                    if (setWriter && setWriter->enabled())
+                    {
+                        mergeAndWrite(*setWriter, mzPoints);
+                    }
+                }
+            }
+        }
+        else
+        {
+            Info<< "    No pointZones found."<<endl;
+        }
+    }
+
+
     // Force creation of all addressing if requested.
     // Errors will be reported as required
     if (allTopology)
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.H b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.H
index e179399257..0ebf1e7497 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.H
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.H
@@ -1,5 +1,7 @@
 #include "labelList.H"
 #include "autoPtr.H"
+#include "ZoneMesh.H"
+#include "topoSet.H"
 
 namespace Foam
 {
@@ -21,6 +23,14 @@ namespace Foam
         pointSet& points
     );
 
+    template<class Zone>
+    label checkZones
+    (
+        const polyMesh& mesh,
+        const ZoneMesh<Zone, polyMesh>& zones,
+        topoSet& set
+    );
+
     label checkTopology
     (
         const polyMesh& mesh,