ENH: add global topology check in checkMesh and makeFaMesh (#2771)

- detect when boundary patches are multiply connected across edges

STYLE: initialize some faMesh values

applications/utilities/finiteArea/makeFaMesh/checkPatchTopology.H

@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    Check topology of poly patches used for finite-area.
+
+Input
+    mesh (polyMesh)
+    meshDefDict (system/faMeshDefintion)
+
+\*---------------------------------------------------------------------------*/
+
+// Preliminary checks
+{
+    typedef typename uindirectPrimitivePatch::surfaceTopo TopoType;
+
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    wordRes polyPatchNames;
+    meshDefDict.readIfPresent("polyMeshPatches", polyPatchNames);
+
+    const labelList patchIDs
+    (
+        pbm.patchSet
+        (
+            polyPatchNames,
+            false,  // warnNotFound
+            true    // useGroups
+        ).sortedToc()
+    );
+
+    label nFaceLabels = 0;
+    for (const label patchi : patchIDs)
+    {
+        nFaceLabels += pbm[patchi].size();
+    }
+
+    labelList faceLabels(nFaceLabels);
+
+    nFaceLabels = 0;
+    for (const label patchi : patchIDs)
+    {
+        for (const label facei : pbm[patchi].range())
+        {
+            faceLabels[nFaceLabels] = facei;
+            ++nFaceLabels;
+        }
+    }
+
+    uindirectPrimitivePatch pp
+    (
+        UIndirectList<face>(mesh.faces(), faceLabels),
+        mesh.points()
+    );
+
+    // Non-manifold
+    labelHashSet badEdges(pp.nEdges()/20);
+
+    labelHashSet* pointSetPtr = nullptr;
+    labelHashSet* badEdgesPtr = &badEdges;
+
+    bool foundError = false;
+
+    if (returnReduceAnd(pp.empty()))
+    {
+        // Empty
+    }
+    else if (UPstream::parRun())
+    {
+        const labelList meshEdges
+        (
+            pp.meshEdges(mesh.edges(), mesh.pointEdges())
+        );
+
+        // Parallel - use mesh edges
+        // - no check for point-pinch
+        // - no check for consistent orientation (if that is posible to
+        //   check?)
+
+        // Count number of edge/face connections (globally)
+        labelList nEdgeConnections(mesh.nEdges(), Zero);
+
+        const labelListList& edgeFaces = pp.edgeFaces();
+
+        forAll(edgeFaces, edgei)
+        {
+            nEdgeConnections[meshEdges[edgei]] = edgeFaces[edgei].size();
+        }
+
+        // Synchronise across coupled edges
+        syncTools::syncEdgeList
+        (
+            mesh,
+            nEdgeConnections,
+            plusEqOp<label>(),
+            label(0)            // null value
+        );
+
+        label labelTyp = TopoType::MANIFOLD;
+        forAll(meshEdges, edgei)
+        {
+            const label meshEdgei = meshEdges[edgei];
+            const label numNbrs = nEdgeConnections[meshEdgei];
+            if (numNbrs == 1)
+            {
+                //if (pointSetPtr) pointSetPtr->insert(mesh.edges()[meshEdgei]);
+                labelTyp = max(labelTyp, TopoType::OPEN);
+            }
+            else if (numNbrs == 0 || numNbrs > 2)
+            {
+                if (pointSetPtr) pointSetPtr->insert(mesh.edges()[meshEdgei]);
+                if (badEdgesPtr) badEdgesPtr->insert(edgei);
+                labelTyp = max(labelTyp, TopoType::ILLEGAL);
+            }
+        }
+        reduce(labelTyp, maxOp<label>());
+
+        foundError = (labelTyp == TopoType::ILLEGAL);
+    }
+    else
+    {
+        const TopoType pTyp = pp.surfaceType(badEdgesPtr);
+
+        foundError = (pTyp == TopoType::ILLEGAL);
+    }
+
+    if (foundError)
+    {
+        edgeList dumpEdges(pp.edges(), badEdges.sortedToc());
+
+        vtk::lineWriter writer
+        (
+            pp.localPoints(),
+            dumpEdges,
+            fileName
+            (
+                mesh.time().globalPath()
+              / ("faMesh-construct.illegalEdges")
+            )
+        );
+
+        writer.writeGeometry();
+
+        // CellData
+        writer.beginCellData();
+        writer.writeProcIDs();
+
+        Info<< "Wrote "
+            << returnReduce(dumpEdges.size(), sumOp<label>())
+            << " bad edges: " << writer.output().name() << nl;
+        writer.close();
+    }
+}
+
+// ************************************************************************* //

applications/utilities/finiteArea/makeFaMesh/decomposeFaFields.H

@@ -16,10 +16,10 @@ Description
 \*---------------------------------------------------------------------------*/
 
 // Embed do-while to support early termination
-if (doDecompose && Pstream::parRun())
+if (doDecompose && UPstream::parRun())
 do
 {
-    faMeshReconstructor reconstructor(aMesh, IOobjectOption::READ_IF_PRESENT);
+    faMeshReconstructor reconstructor(aMesh, IOobjectOption::LAZY_READ);
 
     if (!reconstructor.good())
     {
@@ -61,16 +61,16 @@ do
         auto oldHandler = fileHandler(fileOperation::NewUncollated());
         fileHandler().distributed(true);
 
-        if (Pstream::master())
+        if (UPstream::master())
         {
-            haveMeshOnProc.set(Pstream::myProcNo());
+            haveMeshOnProc.set(UPstream::myProcNo());
             subsetter.reset(new faMeshSubset(serialMesh));
 
-            const bool oldParRun = Pstream::parRun(false);
+            const bool oldParRun = UPstream::parRun(false);
 
             objects = IOobjectList(serialMesh.time(), runTime.timeName());
 
-            Pstream::parRun(oldParRun);
+            UPstream::parRun(oldParRun);
         }
 
         // Restore old settings

applications/utilities/finiteArea/makeFaMesh/faMeshWriteEdgesOBJ.H

@@ -21,12 +21,12 @@ Description
 
     word outputName("finiteArea-edges.obj");
 
-    if (Pstream::parRun())
+    if (UPstream::parRun())
     {
         outputName = word
         (
             "finiteArea-edges-proc"
-          + Foam::name(Pstream::myProcNo())
+          + Foam::name(UPstream::myProcNo())
           + ".obj"
         );
     }

applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,7 +51,9 @@ Original Authors
 #include "faMeshReconstructor.H"
 #include "faMeshSubset.H"
 #include "PtrListOps.H"
+#include "foamVtkLineWriter.H"
 #include "foamVtkIndPatchWriter.H"
+#include "syncTools.H"
 #include "OBJstream.H"
 
 using namespace Foam;
@@ -128,6 +130,8 @@ int main(int argc, char *argv[])
         meshDefDict.add("emptyPatch", patchName, true);
     }
 
+    // Preliminary checks
+    #include "checkPatchTopology.H"
 
     // Create
     faMesh aMesh(mesh, meshDefDict);

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,6 +36,7 @@ License
 #include "IOmanip.H"
 #include "emptyPolyPatch.H"
 #include "processorPolyPatch.H"
+#include "foamVtkLineWriter.H"
 #include "vtkCoordSetWriter.H"
 #include "vtkSurfaceWriter.H"
 #include "checkTools.H"
@@ -44,24 +45,34 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// NEWER CODE
+// ~~~~~~~~~~
+// Return true on error
 template<class PatchType>
-void Foam::checkPatch
+bool Foam::checkPatch
 (
     const bool allGeometry,
-    const word& name,
+    const std::string& name,
     const polyMesh& mesh,
     const PatchType& pp,
-    const labelList& meshFaces,
+    const labelUList& meshEdges,
-    const labelList& meshEdges,
+    labelHashSet* pointSetPtr,
-    pointSet& points
+    labelHashSet* badEdgesPtr
 )
 {
+    if (badEdgesPtr)
+    {
+        badEdgesPtr->clear();
+    }
+
     typedef typename PatchType::surfaceTopo TopoType;
 
+    bool foundError = false;
+
     const label globalSize = returnReduce(pp.size(), sumOp<label>());
 
     Info<< "    "
-        << setw(20) << name
+        << setw(20) << name.c_str()
         << setw(9) << globalSize
         << setw(9) << returnReduce(pp.nPoints(), sumOp<label>());
 
@@ -69,13 +80,169 @@ void Foam::checkPatch
     {
         Info<< setw(34) << "ok (empty)";
     }
-    else if (Pstream::parRun())
+    else if (UPstream::parRun())
     {
         // Parallel - use mesh edges
         // - no check for point-pinch
         // - no check for consistent orientation (if that is posible to
         //   check?)
 
+        // Count number of edge/face connections (globally)
+        labelList nEdgeConnections(mesh.nEdges(), Zero);
+
+        const labelListList& edgeFaces = pp.edgeFaces();
+
+        forAll(edgeFaces, edgei)
+        {
+            nEdgeConnections[meshEdges[edgei]] = edgeFaces[edgei].size();
+        }
+
+        // Synchronise across coupled edges
+        syncTools::syncEdgeList
+        (
+            mesh,
+            nEdgeConnections,
+            plusEqOp<label>(),
+            label(0)            // null value
+        );
+
+        label labelTyp = TopoType::MANIFOLD;
+        forAll(meshEdges, edgei)
+        {
+            const label meshEdgei = meshEdges[edgei];
+            const label numNbrs = nEdgeConnections[meshEdgei];
+            if (numNbrs == 1)
+            {
+                //if (pointSetPtr) pointSetPtr->insert(mesh.edges()[meshEdgei]);
+                labelTyp = max(labelTyp, TopoType::OPEN);
+            }
+            else if (numNbrs == 0 || numNbrs > 2)
+            {
+                if (pointSetPtr) pointSetPtr->insert(mesh.edges()[meshEdgei]);
+                if (badEdgesPtr) badEdgesPtr->insert(edgei);
+                labelTyp = max(labelTyp, TopoType::ILLEGAL);
+            }
+        }
+        reduce(labelTyp, maxOp<label>());
+
+        if (labelTyp == TopoType::MANIFOLD)
+        {
+            Info<< setw(34) << "ok (closed singly connected)";
+        }
+        else if (labelTyp == TopoType::OPEN)
+        {
+            Info<< setw(34) << "ok (non-closed singly connected)";
+        }
+        else
+        {
+            Info<< setw(34) << "multiply connected (shared edge)";
+        }
+
+        foundError = (labelTyp == TopoType::ILLEGAL);
+    }
+    else
+    {
+        const TopoType pTyp = pp.surfaceType(badEdgesPtr);
+
+        if (pTyp == TopoType::MANIFOLD)
+        {
+            if (pp.checkPointManifold(true, pointSetPtr))
+            {
+                Info<< setw(34) << "multiply connected (shared point)";
+            }
+            else
+            {
+                Info<< setw(34) << "ok (closed singly connected)";
+            }
+
+            if (pointSetPtr)
+            {
+                // Add points on non-manifold edges to make set complete
+                pp.checkTopology(false, pointSetPtr);
+            }
+        }
+        else
+        {
+            if (pointSetPtr)
+            {
+                pp.checkTopology(false, pointSetPtr);
+            }
+
+            if (pTyp == TopoType::OPEN)
+            {
+                Info<< setw(34) << "ok (non-closed singly connected)";
+            }
+            else
+            {
+                Info<< setw(34) << "multiply connected (shared edge)";
+            }
+        }
+
+        foundError = (pTyp == TopoType::ILLEGAL);
+    }
+
+    if (allGeometry)
+    {
+        boundBox bb(pp.box());
+        bb.reduce();
+
+        if (bb.good())
+        {
+            Info<< ' ' << bb;
+        }
+    }
+
+    return foundError;
+}
+
+
+// OLDER CODE
+// ~~~~~~~~~~
+// Return true on error
+template<class PatchType>
+bool Foam::checkPatch
+(
+    const bool allGeometry,
+    const std::string& name,
+    const polyMesh& mesh,
+    const PatchType& pp,
+    const labelUList& meshFaces,
+    const labelUList& meshEdges,
+    labelHashSet* pointSetPtr,
+    labelHashSet* badEdgesPtr
+)
+{
+    if (badEdgesPtr)
+    {
+        badEdgesPtr->clear();
+    }
+
+    typedef typename PatchType::surfaceTopo TopoType;
+
+    bool foundError = false;
+
+    const label globalSize = returnReduce(pp.size(), sumOp<label>());
+
+    Info<< "    "
+        << setw(20) << name.c_str()
+        << setw(9) << globalSize
+        << setw(9) << returnReduce(pp.nPoints(), sumOp<label>());
+
+    if (globalSize == 0)
+    {
+        Info<< setw(34) << "ok (empty)";
+    }
+    else if (UPstream::parRun())
+    {
+        // Parallel - use mesh edges
+        // - no check for point-pinch
+        // - no check for consistent orientation (if that is posible to
+        //   check?)
+
+        // OLDER CODE
+        // ~~~~~~~~~~
+        // Synchronise connected faces using global face numbering
+        //
         // (see addPatchCellLayer::globalEdgeFaces)
         // From mesh edge to global face labels. Non-empty sublists only for
         // pp edges.
@@ -88,12 +255,12 @@ void Foam::checkPatch
 
         forAll(edgeFaces, edgei)
         {
-            label meshEdgei = meshEdges[edgei];
+            const label meshEdgei = meshEdges[edgei];
             const labelList& eFaces = edgeFaces[edgei];
 
             // Store face and processor as unique tag.
             labelList& globalEFaces = globalEdgeFaces[meshEdgei];
-            globalEFaces.setSize(eFaces.size());
+            globalEFaces.resize(eFaces.size());
             forAll(eFaces, i)
             {
                 globalEFaces[i] = globalFaces.toGlobal(meshFaces[eFaces[i]]);
@@ -104,29 +271,30 @@ void Foam::checkPatch
             //    << endl;
         }
 
-        // Synchronise across coupled edges.
+        // Synchronise across coupled edges
         syncTools::syncEdgeList
         (
             mesh,
             globalEdgeFaces,
             ListOps::uniqueEqOp<label>(),
-            labelList()             // null value
+            labelList()         // null value
         );
 
-
         label labelTyp = TopoType::MANIFOLD;
         forAll(meshEdges, edgei)
         {
             const label meshEdgei = meshEdges[edgei];
             const labelList& globalEFaces = globalEdgeFaces[meshEdgei];
-            if (globalEFaces.size() == 1)
+            const label numNbrs = globalEFaces.size();
+            if (numNbrs == 1)
             {
-                //points.insert(mesh.edges()[meshEdgei]);
+                //if (pointSetPtr) pointSetPtr->insert(mesh.edges()[meshEdgei]);
                 labelTyp = max(labelTyp, TopoType::OPEN);
             }
-            else if (globalEFaces.size() == 0 || globalEFaces.size() > 2)
+            else if (numNbrs == 0 || numNbrs > 2)
             {
-                points.insert(mesh.edges()[meshEdgei]);
+                if (pointSetPtr) pointSetPtr->insert(mesh.edges()[meshEdgei]);
+                if (badEdgesPtr) badEdgesPtr->insert(edgei);
                 labelTyp = max(labelTyp, TopoType::ILLEGAL);
             }
         }
@@ -138,61 +306,68 @@ void Foam::checkPatch
         }
         else if (labelTyp == TopoType::OPEN)
         {
-            Info<< setw(34)
+            Info<< setw(34) << "ok (non-closed singly connected)";
-                << "ok (non-closed singly connected)";
         }
         else
         {
-            Info<< setw(34)
+            Info<< setw(34) << "multiply connected (shared edge)";
-                << "multiply connected (shared edge)";
         }
+
+        foundError = (labelTyp == TopoType::ILLEGAL);
     }
     else
     {
-        TopoType pTyp = pp.surfaceType();
+        const TopoType pTyp = pp.surfaceType(badEdgesPtr);
 
         if (pTyp == TopoType::MANIFOLD)
         {
-            if (pp.checkPointManifold(true, &points))
+            if (pp.checkPointManifold(true, pointSetPtr))
             {
-                Info<< setw(34)
+                Info<< setw(34) << "multiply connected (shared point)";
-                    << "multiply connected (shared point)";
             }
             else
             {
                 Info<< setw(34) << "ok (closed singly connected)";
             }
 
-            // Add points on non-manifold edges to make set complete
+            if (pointSetPtr)
-            pp.checkTopology(false, &points);
+            {
+                // Add points on non-manifold edges to make set complete
+                pp.checkTopology(false, pointSetPtr);
+            }
         }
         else
         {
-            pp.checkTopology(false, &points);
+            if (pointSetPtr)
+            {
+                pp.checkTopology(false, pointSetPtr);
+            }
 
             if (pTyp == TopoType::OPEN)
             {
-                Info<< setw(34)
+                Info<< setw(34) << "ok (non-closed singly connected)";
-                    << "ok (non-closed singly connected)";
             }
             else
             {
-                Info<< setw(34)
+                Info<< setw(34) << "multiply connected (shared edge)";
-                    << "multiply connected (shared edge)";
             }
         }
+
+        foundError = (pTyp == TopoType::ILLEGAL);
     }
 
     if (allGeometry)
     {
-        const labelList& mp = pp.meshPoints();
+        boundBox bb(pp.box());
+        bb.reduce();
 
-        if (returnReduceOr(mp.size()))
+        if (bb.good())
         {
-            boundBox bb(pp.points(), mp, true); // reduce
             Info<< ' ' << bb;
         }
     }
+
+    return foundError;
 }
 
 
@@ -524,7 +699,7 @@ Foam::label Foam::checkTopology
                 << mesh.time().timeName()/"cellToRegion"
                 << endl;
 
-            labelIOList ctr
+            IOListRef<label>
             (
                 IOobject
                 (
@@ -532,11 +707,11 @@ Foam::label Foam::checkTopology
                     mesh.time().timeName(),
                     mesh,
                     IOobject::NO_READ,
-                    IOobject::NO_WRITE
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
                 ),
                 rs
-            );
+            ).write();
-            ctr.write();
 
 
             // Points in multiple regions
@@ -651,7 +826,7 @@ Foam::label Foam::checkTopology
     }
 
     // Non-manifold points
-    pointSet points
+    pointSet nonManifoldPoints
     (
         mesh,
         "nonManifoldPoints",
@@ -670,17 +845,15 @@ Foam::label Foam::checkTopology
             << setw(20) << "Patch"
             << setw(9) << "Faces"
             << setw(9) << "Points"
-            << "Surface topology";
+            << "  Surface topology";
         if (allGeometry)
         {
-            Info<< " Bounding box";
+            Info<< "  Bounding box";
         }
         Info<< endl;
 
-        forAll(patches, patchi)
+        for (const polyPatch& pp : patches)
         {
-            const polyPatch& pp = patches[patchi];
-
             if (!isA<processorPolyPatch>(pp))
             {
                 checkPatch
@@ -689,14 +862,78 @@ Foam::label Foam::checkTopology
                     pp.name(),
                     mesh,
                     pp,
-                    identity(pp.size(), pp.start()),
+                    // identity(pp.size(), pp.start()),
                     pp.meshEdges(),
-                    points
+                   &nonManifoldPoints
                 );
                 Info<< endl;
             }
         }
 
+        // All non-processor boundary patches
+        {
+            labelList faceLabels
+            (
+                identity
+                (
+                    (
+                        patches.range(patches.nNonProcessor()-1).end_value()
+                      - patches.start()
+                    ),
+                    patches.start()
+                )
+            );
+
+            uindirectPrimitivePatch pp
+            (
+                UIndirectList<face>(mesh.faces(), faceLabels),
+                mesh.points()
+            );
+
+            // Non-manifold
+            labelHashSet badEdges(pp.nEdges()/20);
+
+            bool hadBadEdges = checkPatch
+            (
+                allGeometry,
+                "## ALL BOUNDARIES ##",
+                 mesh,
+                pp,
+                // faceLabels,
+                pp.meshEdges(mesh.edges(), mesh.pointEdges()),
+                nullptr,  // No point set
+               &badEdges
+            );
+            Info<< nl << endl;
+
+            if (hadBadEdges)
+            {
+                edgeList dumpEdges(pp.edges(), badEdges.sortedToc());
+
+                vtk::lineWriter writer
+                (
+                    pp.localPoints(),
+                    dumpEdges,
+                    fileName
+                    (
+                        mesh.time().globalPath()
+                      / ("checkMesh-illegal-edges")
+                    )
+                );
+
+                writer.writeGeometry();
+
+                // CellData
+                writer.beginCellData();
+                writer.writeProcIDs();
+
+                Info<< "Wrote "
+                    << returnReduce(dumpEdges.size(), sumOp<label>())
+                    << " bad edges: " << writer.output().name() << nl;
+                writer.close();
+            }
+        }
+
         //Info.setf(ios_base::right);
     }
 
@@ -729,9 +966,9 @@ Foam::label Foam::checkTopology
                     fz.name(),
                     mesh,
                     fz(),           // patch
-                    fz,             // mesh face labels
+                    // fz,             // mesh face labels
                     fz.meshEdges(), // mesh edge labels
-                    points
+                   &nonManifoldPoints
                 );
                 Info<< endl;
             }
@@ -761,17 +998,20 @@ Foam::label Foam::checkTopology
         }
     }
 
-    const label nPoints = returnReduce(points.size(), sumOp<label>());
+
+    const label nPoints =
+        returnReduce(nonManifoldPoints.size(), sumOp<label>());
 
     if (nPoints)
     {
         Info<< "  <<Writing " << nPoints
-            << " conflicting points to set " << points.name() << endl;
+            << " conflicting points to set " << nonManifoldPoints.name()
-        points.instance() = mesh.pointsInstance();
+            << endl;
-        points.write();
+        nonManifoldPoints.instance() = mesh.pointsInstance();
+        nonManifoldPoints.write();
         if (setWriter && setWriter->enabled())
         {
-            mergeAndWrite(*setWriter, points);
+            mergeAndWrite(*setWriter, nonManifoldPoints);
         }
     }
 

applications/utilities/mesh/manipulation/checkMesh/checkTopology.H

@@ -1,3 +1,20 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    Global functions for topology checking
+
+\*---------------------------------------------------------------------------*/
+
 #include "labelList.H"
 #include "autoPtr.H"
 #include "ZoneMesh.H"
@@ -5,38 +22,62 @@
 
 namespace Foam
 {
-    // Forward Declarations
-    class polyMesh;
-    class pointSet;
-    class coordSetWriter;
-    class surfaceWriter;
 
-    template<class PatchType>
+// Forward Declarations
-    void checkPatch
+class polyMesh;
-    (
+class pointSet;
-        const bool allGeometry,
+class coordSetWriter;
-        const word& name,
+class surfaceWriter;
-        const polyMesh& mesh,
-        const PatchType& pp,
-        const labelList& meshFaces,
-        const labelList& meshEdges,
-        pointSet& points
-    );
 
-    template<class Zone>
-    label checkZones
-    (
-        const polyMesh& mesh,
-        const ZoneMesh<Zone, polyMesh>& zones,
-        topoSet& set
-    );
 
-    label checkTopology
+// Check patch topology.
-    (
+// In parallel, uses count of edge connections
-        const polyMesh& mesh,
+template<class PatchType>
-        const bool allTopology,
+bool checkPatch
-        const bool allGeometry,
+(
-        autoPtr<surfaceWriter>& surfWriter,
+    const bool allGeometry,
-        autoPtr<coordSetWriter>& setWriter
+    const std::string& name,
-    );
+    const polyMesh& mesh,
-}
+    const PatchType& pp,
+    const labelUList& meshEdges,
+    labelHashSet* pointSetPtr = nullptr,
+    labelHashSet* badEdgesPtr = nullptr
+);
+
+// OLDER CODE
+// ~~~~~~~~~~
+// Check patch topology.
+// In parallel, uses mesh face ids and global face numbering
+template<class PatchType>
+bool checkPatch
+(
+    const bool allGeometry,
+    const std::string& name,
+    const polyMesh& mesh,
+    const PatchType& pp,
+    const labelUList& meshFaces,
+    const labelUList& meshEdges,
+    labelHashSet* pointSetPtr = nullptr,
+    labelHashSet* badEdgesPtr = nullptr
+);
+
+template<class Zone>
+label checkZones
+(
+    const polyMesh& mesh,
+    const ZoneMesh<Zone, polyMesh>& zones,
+    topoSet& set
+);
+
+label checkTopology
+(
+    const polyMesh& mesh,
+    const bool allTopology,
+    const bool allGeometry,
+    autoPtr<surfaceWriter>& surfWriter,
+    autoPtr<coordSetWriter>& setWriter
+);
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteArea/faMesh/faMesh.C

@@ -163,6 +163,11 @@ void Foam::faMesh::initPatch() const
             mesh().points()
         )
     );
+    // Could set some basic primitive data here...
+    // nEdges_ = patchPtr_->nEdges();
+    // nInternalEdges_ = patchPtr_->nInternalEdges();
+    // nFaces_ = patchPtr_->size();
+    // nPoints_ = patchPtr_->nPoints();
     bndConnectPtr_.reset(nullptr);
     haloMapPtr_.reset(nullptr);
     haloFaceCentresPtr_.reset(nullptr);
@@ -512,7 +517,13 @@ Foam::faMesh::faMesh
     haloMapPtr_(nullptr),
     haloFaceCentresPtr_(nullptr),
     haloFaceNormalsPtr_(nullptr)
-{}
+{
+    // Not yet much for primitive mesh data possible...
+    nPoints_ = 0;
+    nEdges_ = 0;
+    nInternalEdges_ = 0;
+    nFaces_ = faceLabels_.size();
+}
 
 
 Foam::faMesh::faMesh
@@ -593,7 +604,13 @@ Foam::faMesh::faMesh
     haloMapPtr_(nullptr),
     haloFaceCentresPtr_(nullptr),
     haloFaceNormalsPtr_(nullptr)
-{}
+{
+    // Not yet much for primitive mesh data possible...
+    nPoints_ = 0;
+    nEdges_ = 0;
+    nInternalEdges_ = 0;
+    nFaces_ = faceLabels_.size();
+}
 
 
 Foam::faMesh::faMesh(const polyPatch& pp, const bool doInit)

src/finiteArea/faMesh/faMeshTopology.C

@@ -445,9 +445,9 @@ Foam::faMesh::getBoundaryEdgeConnections() const
                 // A boundary finiteEdge edge (known from this side)
 
                 auto& gathered = gatheredConnections[cppEdgei];
-                gathered.setCapacity(2);
+                gathered.setCapacity_nocopy(2);
-                gathered.resize(1);
+                gathered.resize_nocopy(1);
-                auto& tuple = gathered.back();
+                auto& tuple = gathered.front();
 
                 tuple = bndEdgeConnections[bndEdgei].first();
             }
@@ -502,9 +502,9 @@ Foam::faMesh::getBoundaryEdgeConnections() const
                 const label meshFacei = patchFacei + pp.start();
 
                 auto& gathered = gatheredConnections[cppEdgei];
-                gathered.setCapacity(2);
+                gathered.setCapacity_nocopy(2);
-                gathered.resize(1);
+                gathered.resize_nocopy(1);
-                auto& tuple = gathered.back();
+                auto& tuple = gathered.front();
 
                 tuple.procNo(UPstream::myProcNo());
                 tuple.patchi(patchi);
@@ -545,6 +545,7 @@ Foam::faMesh::getBoundaryEdgeConnections() const
         << "Collating sync information" << endl;
 
     // Pick out gathered connections and add into primary bookkeeping
+    badEdges.clear();
     for (label cppEdgei = 0; cppEdgei < nCoupledEdges; ++cppEdgei)
     {
         const auto& gathered = gatheredConnections[cppEdgei];
@@ -552,25 +553,61 @@ Foam::faMesh::getBoundaryEdgeConnections() const
         const label bndEdgei =
             edgeToBoundaryIndex.lookup(cpp.meshEdge(cppEdgei), -1);
 
-        if
+        if (bndEdgei != -1)
-        (
-            // A boundary finiteEdge edge (known from this side)
-            bndEdgei != -1
-
-            // Gathered a one-to-one connection
-         && gathered.size() == 2
-        )
         {
-            const auto& a = gathered[0];
+            // A boundary finiteEdge edge (known from this side)
-            const auto& b = gathered[1];
-
-            // Copy second side of connection
             auto& connection = bndEdgeConnections[bndEdgei];
 
-            connection.second() = (connection.first() == b) ? a : b;
+            if (gathered.size() == 2)
+            {
+                // Copy second side of connection
+                const auto& a = gathered[0];
+                const auto& b = gathered[1];
+
+                connection.second() = (connection.first() == b) ? a : b;
+            }
+            else if (gathered.size() > 2)
+            {
+                // Multiply connected!! - this needs to be addressed
+                badEdges.insert(cppEdgei);
+            }
         }
     }
 
+    if (returnReduceOr(badEdges.size()))
+    {
+        WarningInFunction
+            << nl << "Multiply connected edges detected" << endl;
+
+        // Print out edges as point pairs
+        // These are globally synchronised - so only output on master
+        constexpr label maxOutput = 10;
+
+        label nOutput = 0;
+
+        for (const label cppEdgei : badEdges.sortedToc())
+        {
+            const edge e(cpp.meshEdge(cppEdgei));
+
+            const auto& gathered = gatheredConnections[cppEdgei];
+
+            Info<< "connection: ";
+            gathered.writeList(Info) << nl;
+
+            Info<<"    edge  : "
+                << cpp.points()[e.first()] << ' '
+                << cpp.points()[e.second()] << nl;
+
+            ++nOutput;
+            if (maxOutput > 0 && nOutput >= maxOutput)
+            {
+                Info<< " ... suppressing further output" << nl;
+                break;
+            }
+        }
+    }
+
+
     // Check missing/invalid
     badEdges.clear();
     for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei)