From 8eef91c5e2b423ff22686a166dcca142421e4746 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Sat, 1 May 2021 16:09:48 +0200
Subject: [PATCH 1/8] ENH: improve edge access for face/triFace

- additional rcEdge(), rcEdges() methods for reverse order walk

- accept generic edge() method as alternative to faceEdge() for
  single edge retrieval.

- edge() method with points -> returns the vector

- reduce the number of operations in edgeDirection methods

DEFEATURE: remove longestEdge global function

- deprecated and replaced by face::longestEdge() method (2017-04)
---
 applications/test/faces/Test-faces.C          | 103 +++++++++++-
 src/OpenFOAM/meshes/meshShapes/face/face.C    | 158 ++++++++----------
 src/OpenFOAM/meshes/meshShapes/face/face.H    |  58 ++++---
 src/OpenFOAM/meshes/meshShapes/face/faceI.H   |  78 +++++----
 .../meshes/meshShapes/triFace/triFace.H       |  54 ++++--
 .../meshes/meshShapes/triFace/triFaceI.H      | 146 +++++++++++-----
 src/conversion/fire/checkFireEdges.C          |   6 +-
 .../polyTopoChange/edgeCollapser.C            |  10 +-
 .../cfdTools/general/levelSet/levelSet.C      |   5 +-
 .../general/levelSet/levelSetTemplates.C      |   5 +-
 src/lagrangian/basic/particle/particle.C      |   2 +-
 .../surfaceIntersectionFuncs.C                |   4 +-
 .../cutting/cuttingSurfaceBaseTemplates.C     |   2 +-
 13 files changed, 412 insertions(+), 219 deletions(-)

diff --git a/applications/test/faces/Test-faces.C b/applications/test/faces/Test-faces.C
index b967f4b5ba..7909cac61b 100644
--- a/applications/test/faces/Test-faces.C
+++ b/applications/test/faces/Test-faces.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,6 +33,8 @@ Description
 
 #include "argList.H"
 #include "labelledTri.H"
+#include "faceList.H"
+#include "triFaceList.H"
 #include "pointList.H"
 #include "ListOps.H"
 
@@ -63,6 +65,35 @@ void testSign
 }
 
 
+template<class Face>
+void testEdges(const Face& f)
+{
+    const label nEdges = f.nEdges();
+    Info<< "face: " << f << nl
+        << "flip: " << f.reverseFace() << nl
+        << "  fc edges:" << flatOutput(f.edges()) << nl
+        << "  rc edges:" << flatOutput(f.rcEdges()) << nl;
+
+    Info<< "  forward edges" << nl;
+    for (label edgei = 0; edgei < nEdges; ++edgei)
+    {
+        Info<< "    " << edgei << " : " << f.edge(edgei) << nl;
+    }
+    Info<< "  reverse edges" << nl;
+    for (label edgei = 0; edgei < nEdges; ++edgei)
+    {
+        Info<< "    " << edgei << " : " << f.rcEdge(edgei) << nl;
+    }
+}
+
+
+void testCompare(const triFace& a, const triFace& b)
+{
+    Info<< "compare: " << a << " with " << b
+        <<  " == " << triFace::compare(a, b) << nl;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 //  Main program:
 
@@ -92,7 +123,7 @@ int main(int argc, char *argv[])
         { 2, 2, 2}
     });
 
-    face f1{ 1, 2, 3, 4 };
+    face f1({1, 2, 3, 4});
     Info<< "face:"; faceInfo(f1, points1); Info << nl;
     testSign(f1, points1, testPoints);
 
@@ -100,7 +131,7 @@ int main(int argc, char *argv[])
     testSign(f1, points2, testPoints);
     Info<< nl;
 
-    triFace t1{ 1, 2, 3 };
+    triFace t1({1, 2, 3});
     Info<< "triFace:"; faceInfo(t1, points1); Info << nl;
     testSign(t1, points1, testPoints);
 
@@ -115,11 +146,12 @@ int main(int argc, char *argv[])
     f1 = t1.triFaceFace();
     Info<< "face:" << f1 << nl;
 
-    // expect these to fail
+    #if 0
+    // Expect failure, but triggers abort which cannot be caught
     const bool throwingError = FatalError.throwExceptions();
     try
     {
-        labelledTri l1{ 1, 2, 3, 10, 24 };
+        labelledTri l1({1, 2, 3, 10, 24});
         Info<< "labelled:" << l1 << nl;
     }
     catch (const Foam::error& err)
@@ -128,11 +160,12 @@ int main(int argc, char *argv[])
             << "Caught FatalError " << err << nl << endl;
     }
     FatalError.throwExceptions(throwingError);
+    #endif
 
-    labelledTri l2{ 1, 2, 3 };
+    labelledTri l2({1, 2, 3});
     Info<< "labelled:" << l2 << nl;
 
-    labelledTri l3{ 1, 2, 3, 10 };
+    labelledTri l3({1, 2, 3, 10});
     Info<< "labelled:" << l3 << nl;
 
     t1.flip();
@@ -141,6 +174,62 @@ int main(int argc, char *argv[])
     Info<< "flip:" << t1 << nl;
     Info<< "flip:" << l3 << nl;
 
+    {
+        triFaceList faceList1
+        ({
+            triFace{1, 2, 3},
+            triFace{4, 2, 100},
+            triFace{1, 3, 2},
+        });
+
+        Info<< nl << "Test edges" << nl;
+
+        for (const auto& f : faceList1)
+        {
+            testEdges(f);
+            Info<< nl;
+        }
+    }
+
+    {
+        faceList faceList1
+        ({
+            face{1, 2, 3, 4},
+            face{1, 4, 3, 2},
+            face{4, 2, 100, 8, 35},
+        });
+
+        Info<< nl << "Test edges" << nl;
+
+        for (const auto& f : faceList1)
+        {
+            testEdges(f);
+            Info<< nl;
+        }
+    }
+
+    {
+        triFaceList faceList1
+        ({
+            triFace{1, 2, 3},
+            triFace{1, 3, 2},
+            triFace{3, 1, 2},
+            triFace{4, 5, 1},
+        });
+
+        Info<< nl << "Test triFace compare" << nl;
+        for (const triFace& a : faceList1)
+        {
+            for (const triFace& b : faceList1)
+            {
+                testCompare(a, b);
+            }
+        }
+        Info<< nl;
+    }
+
+    Info<< "\nEnd\n" << endl;
+
     return 0;
 }
 
diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.C b/src/OpenFOAM/meshes/meshShapes/face/face.C
index d7fd7c3d07..cfd4d85516 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.C
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.C
@@ -31,6 +31,7 @@ License
 #include "triPointRef.H"
 #include "mathematicalConstants.H"
 #include "ConstCirculator.H"
+#include <algorithm>
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -40,25 +41,18 @@ const char* const Foam::face::typeName = "face";
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::tmp<Foam::vectorField>
-Foam::face::calcEdges(const UList<point>& points) const
+Foam::face::calcEdgeVectors(const UList<point>& points) const
 {
-    tmp<vectorField> tedges(new vectorField(size()));
-    vectorField& edges = tedges.ref();
+    auto tedgeVecs = tmp<vectorField>::New(size());
+    auto& edgeVecs = tedgeVecs.ref();
 
     forAll(*this, i)
     {
-        label ni = fcIndex(i);
-
-        point thisPt = points[operator[](i)];
-        point nextPt = points[operator[](ni)];
-
-        vector vec(nextPt - thisPt);
-        vec /= Foam::mag(vec) + VSMALL;
-
-        edges[i] = vec;
+        edgeVecs[i] = vector(points[nextLabel(i)] - points[thisLabel(i)]);
+        edgeVecs[i].normalise();
     }
 
-    return tedges;
+    return tedgeVecs;
 }
 
 
@@ -68,11 +62,11 @@ Foam::scalar Foam::face::edgeCos
     const label index
 ) const
 {
-    label leftEdgeI = left(index);
-    label rightEdgeI = right(index);
+    const vector& leftEdge = edges[rcIndex(index)];
+    const vector& rightEdge = edges[index];
 
     // Note negate on left edge to get correct left-pointing edge.
-    return -(edges[leftEdgeI] & edges[rightEdgeI]);
+    return -(leftEdge & rightEdge);
 }
 
 
@@ -90,12 +84,12 @@ Foam::label Foam::face::mostConcaveAngle
 
     forAll(edges, i)
     {
-        label leftEdgeI = left(i);
-        label rightEdgeI = right(i);
+        const vector& leftEdge = edges[rcIndex(i)];
+        const vector& rightEdge = edges[i];
 
-        vector edgeNormal = edges[rightEdgeI] ^ edges[leftEdgeI];
+        vector edgeNormal = (rightEdge ^ leftEdge);
 
-        scalar edgeCos = edges[leftEdgeI] & edges[rightEdgeI];
+        scalar edgeCos = (leftEdge & rightEdge);
         scalar edgeAngle = acos(max(-1.0, min(1.0, edgeCos)));
 
         scalar angle;
@@ -133,15 +127,15 @@ Foam::label Foam::face::split
     faceList& quadFaces
 ) const
 {
-    label oldIndices = (triI + quadI);
+    const label oldIndices = (triI + quadI);
 
-    if (size() <= 2)
+    if (size() < 3)
     {
         FatalErrorInFunction
             << "Serious problem: asked to split a face with < 3 vertices"
             << abort(FatalError);
     }
-    if (size() == 3)
+    else if (size() == 3)
     {
         // Triangle. Just copy.
         if (mode == COUNTTRIANGLE || mode == COUNTQUAD)
@@ -166,7 +160,7 @@ Foam::label Foam::face::split
         else if (mode == SPLITTRIANGLE)
         {
             //  Start at point with largest internal angle.
-            const vectorField edges(calcEdges(points));
+            const vectorField edges(calcEdgeVectors(points));
 
             scalar minAngle;
             label startIndex = mostConcaveAngle(points, edges, minAngle);
@@ -197,13 +191,13 @@ Foam::label Foam::face::split
     {
         // General case. Like quad: search for largest internal angle.
 
-        const vectorField edges(calcEdges(points));
+        const vectorField edges(calcEdgeVectors(points));
 
         scalar minAngle = 1;
         label startIndex = mostConcaveAngle(points, edges, minAngle);
 
         scalar bisectAngle = minAngle/2;
-        vector rightEdge = edges[right(startIndex)];
+        const vector& rightEdge = edges[startIndex];
 
         //
         // Look for opposite point which as close as possible bisects angle
@@ -222,7 +216,7 @@ Foam::label Foam::face::split
                 points[operator[](index)]
               - points[operator[](startIndex)]
             );
-            splitEdge /= Foam::mag(splitEdge) + VSMALL;
+            splitEdge.normalise();
 
             const scalar splitCos = splitEdge & rightEdge;
             const scalar splitAngle = acos(max(-1.0, min(1.0, splitCos)));
@@ -786,62 +780,57 @@ Foam::tensor Foam::face::inertia
 
 Foam::edgeList Foam::face::edges() const
 {
-    const labelList& points = *this;
+    const labelList& verts = *this;
+    const label nVerts = verts.size();
 
-    edgeList e(points.size());
+    edgeList theEdges(nVerts);
 
-    for (label pointi = 0; pointi < points.size() - 1; ++pointi)
+    // Last edge closes the polygon
+    theEdges.last().first() = verts.last();
+    theEdges.last().second() = verts[0];
+
+    for (label verti = 0; verti < nVerts - 1; ++verti)
     {
-        e[pointi] = edge(points[pointi], points[pointi + 1]);
+        theEdges[verti].first() = verts[verti];
+        theEdges[verti].second() = verts[verti + 1];
     }
 
-    // Add last edge
-    e.last() = edge(points.last(), points[0]);
-
-    return e;
+    return theEdges;
 }
 
 
-int Foam::face::edgeDirection(const edge& e) const
+Foam::edgeList Foam::face::rcEdges() const
 {
-    forAll(*this, i)
+    const labelList& verts = *this;
+    const label nVerts = verts.size();
+
+    edgeList theEdges(nVerts);
+
+    // First edge closes the polygon
+    theEdges.first().first() = verts[0];
+    theEdges.first().second() = verts.last();
+
+    for (label verti = 1; verti < nVerts; ++verti)
     {
-        if (operator[](i) == e.first())
-        {
-            if (operator[](rcIndex(i)) == e.second())
-            {
-                // Reverse direction
-                return -1;
-            }
-            else if (operator[](fcIndex(i)) == e.second())
-            {
-                // Forward direction
-                return 1;
-            }
-
-            // No match
-            return 0;
-        }
-        else if (operator[](i) == e.second())
-        {
-            if (operator[](rcIndex(i)) == e.first())
-            {
-                // Forward direction
-                return 1;
-            }
-            else if (operator[](fcIndex(i)) == e.first())
-            {
-                // Reverse direction
-                return -1;
-            }
-
-            // No match
-            return 0;
-        }
+        theEdges[verti].first() = verts[nVerts - verti];
+        theEdges[verti].second() = verts[nVerts - verti - 1];
     }
 
-    // Not found
-    return 0;
+    return theEdges;
+}
+
+
+int Foam::face::edgeDirection(const Foam::edge& e) const
+{
+    const label idx = find(e.first());
+
+    if (idx != -1)
+    {
+        if (e.second() == nextLabel(idx)) return 1;  // Forward
+        if (e.second() == prevLabel(idx)) return -1; // Reverse
+    }
+
+    return 0;  // Not found
 }
 
 
@@ -894,29 +883,26 @@ Foam::label Foam::face::trianglesQuads
 
 Foam::label Foam::face::longestEdge(const UList<point>& pts) const
 {
-    const edgeList& eds = this->edges();
+    const labelList& verts = *this;
+    const label nVerts = verts.size();
 
-    label longestEdgeI = -1;
-    scalar longestEdgeLength = -SMALL;
+    // Last edge closes the polygon. Use it to initialize loop
+    label longest = nVerts - 1;
+    scalar longestLen = Foam::edge(verts.first(), verts.last()).mag(pts);
 
-    forAll(eds, edI)
+    // Examine other edges
+    for (label edgei = 0; edgei < nVerts - 1; ++edgei)
     {
-        scalar edgeLength = eds[edI].mag(pts);
+        scalar edgeLen = Foam::edge(verts[edgei], verts[edgei + 1]).mag(pts);
 
-        if (edgeLength > longestEdgeLength)
+        if (longestLen < edgeLen)
         {
-            longestEdgeI = edI;
-            longestEdgeLength = edgeLength;
+            longest = edgei;
+            longestLen = edgeLen;
         }
     }
 
-    return longestEdgeI;
-}
-
-
-Foam::label Foam::longestEdge(const face& f, const UList<point>& pts)
-{
-    return f.longestEdge(pts);
+    return longest;
 }
 
 
diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.H b/src/OpenFOAM/meshes/meshShapes/face/face.H
index 9c5565fcf8..8862402896 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.H
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.H
@@ -76,14 +76,8 @@ class face
 {
     // Private Member Functions
 
-        //- Edge to the right of face vertex i
-        inline label right(const label i) const;
-
-        //- Edge to the left of face vertex i
-        inline label left(const label i) const;
-
         //- Construct list of edge vectors for face
-        tmp<vectorField> calcEdges
+        tmp<vectorField> calcEdgeVectors
         (
             const UList<point>& points
         ) const;
@@ -190,7 +184,7 @@ public:
         void flip();
 
         //- Return the points corresponding to this face
-        inline pointField points(const UList<point>& points) const;
+        inline pointField points(const UList<point>& pts) const;
 
         //- Centre point of face
         point centre(const UList<point>& points) const;
@@ -227,13 +221,14 @@ public:
 
         // Navigation through face vertices
 
-        //- Return true if the point label is found in face.
-        inline bool found(const label pointLabel) const;
-
-        //- Find local index on face for the point label,
+        //- Find local index on face for the point label, same as find()
         //  \return position in face (0,1,2,...) or -1 if not found.
         inline label which(const label pointLabel) const;
 
+        //- The vertex on face - identical to operator[], but with naming
+        //- similar to nextLabel(), prevLabel()
+        inline label thisLabel(const label i) const;
+
         //- Next vertex on face
         inline label nextLabel(const label i) const;
 
@@ -342,25 +337,45 @@ public:
         ) const;
 
         //- Return number of edges
-        inline label nEdges() const;
+        inline label nEdges() const noexcept;
 
-        //- Return edges in face point ordering,
+        //- Return i-th face edge in forward walk order.
+        //- The faceEdge(0) is the edge between [0] and [1]
+        inline Foam::edge faceEdge(const label edgei) const;
+
+        //- Return i-th face edge in forward walk order.
+        //- Identical to faceEdge() but with generic name
+        inline Foam::edge edge(const label edgei) const;
+
+        //- Return vector of i-th face edge in forward walk order.
+        inline vector edge(const label edgei, const UList<point>& pts) const;
+
+        //- Return i-th face edge in reverse walk order.
+        //- The rcEdge(0) is the edge between [0] and [n-1]
+        inline Foam::edge rcEdge(const label edgei) const;
+
+        //- Return vector of i-th face edge in reverse walk order.
+        inline vector rcEdge(const label edgei, const UList<point>& pts) const;
+
+        //- Return list of edges in forward walk order.
         //  i.e. edges()[0] is edge between [0] and [1]
         edgeList edges() const;
 
-        //- Return n-th face edge
-        inline edge faceEdge(const label n) const;
+        //- Return list of edges in reverse walk order.
+        //  i.e. rcEdges()[0] is edge between [0] and [n-1]
+        edgeList rcEdges() const;
 
-        //- The edge direction on the face
+        //- Test the edge direction on the face
         //  \return
         //  -  0: edge not found on the face
         //  - +1: forward (counter-clockwise) on the face
         //  - -1: reverse (clockwise) on the face
-        int edgeDirection(const edge& e) const;
+        int edgeDirection(const Foam::edge& e) const;
 
         //- Find the longest edge on a face.
         label longestEdge(const UList<point>& pts) const;
 
+
         // Face splitting utilities
 
             //- Number of triangles after splitting
@@ -493,13 +508,6 @@ struct offsetOp<face>
 };
 
 
-//- Deprecated(2017-04) find the longest edge on a face.
-//- Face point labels index into pts.
-//  \deprecated(2017-04) use class method instead
-FOAM_DEPRECATED_FOR(2017-04, "use face::longestEdge() method")
-label longestEdge(const face& f, const UList<point>& pts);
-
-
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
 inline bool operator==(const face& a, const face& b);
diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceI.H b/src/OpenFOAM/meshes/meshShapes/face/faceI.H
index 631d8cad45..fe2ca3f754 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/faceI.H
+++ b/src/OpenFOAM/meshes/meshShapes/face/faceI.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,20 +26,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-inline Foam::label Foam::face::right(const label i) const
-{
-    return i;
-}
-
-
-inline Foam::label Foam::face::left(const label i) const
-{
-    return rcIndex(i);
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 inline Foam::face::face(const label sz)
@@ -98,23 +84,19 @@ inline Foam::face::face(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::pointField Foam::face::points
-(
-    const UList<point>& meshPoints
-) const
+inline Foam::pointField Foam::face::points(const UList<point>& pts) const
 {
     // There are as many points as there are labels for them
     pointField p(size());
 
-    // For each point in list, set it to the point in 'pnts' addressed
-    // by 'labs'
-    label i = 0;
+    auto iter = p.begin();
+
     for (const label pointi : *this)
     {
-        p[i++] = meshPoints[pointi];
+        *iter = pts[pointi];
+        ++iter;
     }
 
-    // Return list
     return p;
 }
 
@@ -133,22 +115,54 @@ inline Foam::scalar Foam::face::mag(const UList<point>& p) const
 }
 
 
-inline Foam::label Foam::face::nEdges() const
+inline Foam::label Foam::face::nEdges() const noexcept
 {
     // for a closed polygon a number of edges is the same as number of points
     return size();
 }
 
 
-inline Foam::edge Foam::face::faceEdge(const label n) const
+inline Foam::edge Foam::face::faceEdge(const label edgei) const
 {
-    return edge(operator[](n), operator[](fcIndex(n)));
+    return Foam::edge(thisLabel(edgei), nextLabel(edgei));
 }
 
 
-inline bool Foam::face::found(const label pointLabel) const
+inline Foam::edge Foam::face::edge(const label edgei) const
 {
-    return labelList::found(pointLabel);
+    return faceEdge(edgei);
+}
+
+
+inline Foam::vector Foam::face::edge
+(
+    const label edgei,
+    const UList<point>& pts
+) const
+{
+    return vector(pts[nextLabel(edgei)] - pts[thisLabel(edgei)]);
+}
+
+
+inline Foam::edge Foam::face::rcEdge(const label edgei) const
+{
+    // Edge 0 (forward and reverse) always starts at [0]
+    // for consistency with face flipping
+    const label pointi = edgei ? (nEdges() - edgei) : 0;
+    return Foam::edge(thisLabel(pointi), prevLabel(pointi));
+}
+
+
+inline Foam::vector Foam::face::rcEdge
+(
+    const label edgei,
+    const UList<point>& pts
+) const
+{
+    // Edge 0 (forward and reverse) always starts at [0]
+    // for consistency with face flipping
+    const label pointi = edgei ? (3 - edgei) : 0;
+    return vector(pts[prevLabel(pointi)] - pts[thisLabel(pointi)]);
 }
 
 
@@ -158,6 +172,12 @@ inline Foam::label Foam::face::which(const label pointLabel) const
 }
 
 
+inline Foam::label Foam::face::thisLabel(const label i) const
+{
+    return labelList::operator[](i);
+}
+
+
 inline Foam::label Foam::face::nextLabel(const label i) const
 {
     return labelList::fcValue(i);
diff --git a/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H b/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H
index 50fd55fe87..8083f16e9c 100644
--- a/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H
+++ b/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H
@@ -116,7 +116,7 @@ public:
         inline void flip();
 
         //- Return the points corresponding to this face
-        inline pointField points(const UList<point>& points) const;
+        inline pointField points(const UList<point>& pts) const;
 
         //- Return triangle as a face
         inline face triFaceFace() const;
@@ -148,20 +148,27 @@ public:
         //- Magnitude of face area
         inline scalar mag(const UList<point>& points) const;
 
-        //- Number of triangles after splitting
-        inline label nTriangles() const;
+        //- Number of triangles after splitting == 1
+        inline label nTriangles() const noexcept;
 
         //- Return face with reverse direction
         //  The starting points of the original and reverse face are identical.
         inline triFace reverseFace() const;
 
-        //- Return true if the point label is found in face.
-        inline bool found(const label pointLabel) const;
-
-        //- Find local index on face for the point label.
+        //- Find local index on face for the point label, same as find()
         //  \return position in face (0,1,2) or -1 if not found.
         inline label which(const label pointLabel) const;
 
+        //- Next vertex on face
+        inline label nextLabel(const label i) const;
+
+        //- Previous vertex on face
+        inline label prevLabel(const label i) const;
+
+        //- The vertex on face - identical to operator[], but with naming
+        //- similar to nextLabel(), prevLabel()
+        inline label thisLabel(const label i) const;
+
         //- Return swept-volume from old-points to new-points
         inline scalar sweptVol
         (
@@ -243,22 +250,41 @@ public:
             const scalar tol = SMALL
         ) const;
 
-        //- Return number of edges
-        inline label nEdges() const;
+        //- Return number of edges == 3
+        inline label nEdges() const noexcept;
 
-        //- Return edges in face point ordering,
+        //- Return i-th face edge in forward walk order.
+        //- The faceEdge(0) is the edge between [0] and [1]
+        inline Foam::edge faceEdge(const label edgei) const;
+
+        //- Return i-th face edge in forward walk order.
+        //- Identical to faceEdge() but with generic name
+        inline Foam::edge edge(const label edgei) const;
+
+        //- Return vector of i-th face edge in forward walk order.
+        inline vector edge(const label edgei, const UList<point>& pts) const;
+
+        //- Return i-th face edge in reverse walk order.
+        //- The rcEdge(0) is the edge between [0] to [n-1]
+        inline Foam::edge rcEdge(const label edgei) const;
+
+        //- Return vector of i-th face edge in reverse walk order.
+        inline vector rcEdge(const label edgei, const UList<point>& pts) const;
+
+        //- Return list of edges in forward walk order.
         //  i.e. edges()[0] is edge between [0] and [1]
         inline edgeList edges() const;
 
-        //- Return n-th face edge
-        inline edge faceEdge(const label n) const;
+        //- Return list of edges in reverse walk order.
+        //  i.e. rcEdges()[0] is edge between [0] to [n-1]
+        inline edgeList rcEdges() const;
 
-        //- Return the edge direction on the face
+        //- Test the edge direction on the face
         //  \return
         //  - +1: forward (counter-clockwise) on the face
         //  - -1: reverse (clockwise) on the face
         //  -  0: edge not found on the face
-        inline int edgeDirection(const edge& e) const;
+        inline int edgeDirection(const Foam::edge& e) const;
 
         //- Compare triFaces
         //  \return:
diff --git a/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H b/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H
index e1f2c00869..40d70f246b 100644
--- a/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H
+++ b/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H
@@ -29,7 +29,7 @@ License
 #include "IOstreams.H"
 #include "face.H"
 #include "triPointRef.H"
-#include "Swap.H"
+#include <algorithm>  // For std::swap
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
@@ -144,13 +144,13 @@ inline void Foam::triFace::flip()
 }
 
 
-inline Foam::pointField Foam::triFace::points(const UList<point>& points) const
+inline Foam::pointField Foam::triFace::points(const UList<point>& pts) const
 {
     pointField p(3);
 
-    p[0] = points[operator[](0)];
-    p[1] = points[operator[](1)];
-    p[2] = points[operator[](2)];
+    p[0] = pts[operator[](0)];
+    p[1] = pts[operator[](1)];
+    p[2] = pts[operator[](2)];
 
     return p;
 }
@@ -209,7 +209,7 @@ inline Foam::scalar Foam::triFace::mag(const UList<point>& points) const
 }
 
 
-inline Foam::label Foam::triFace::nTriangles() const
+inline Foam::label Foam::triFace::nTriangles() const noexcept
 {
     return 1;
 }
@@ -222,18 +222,30 @@ inline Foam::triFace Foam::triFace::reverseFace() const
 }
 
 
-inline bool Foam::triFace::found(const label pointLabel) const
-{
-    return FixedList<label, 3>::found(pointLabel);
-}
-
-
 inline Foam::label Foam::triFace::which(const label pointLabel) const
 {
     return FixedList<label, 3>::find(pointLabel);
 }
 
 
+inline Foam::label Foam::triFace::thisLabel(const label i) const
+{
+    return operator[](i);
+}
+
+
+inline Foam::label Foam::triFace::nextLabel(const label i) const
+{
+    return operator[]((i == 2 ? 0 : i+1));
+}
+
+
+inline Foam::label Foam::triFace::prevLabel(const label i) const
+{
+    return operator[]((i ? i-1 : 2));
+}
+
+
 inline Foam::scalar Foam::triFace::sweptVol
 (
     const UList<point>& opts,
@@ -353,61 +365,109 @@ inline int Foam::triFace::sign
 }
 
 
-inline Foam::label Foam::triFace::nEdges() const
+inline Foam::label Foam::triFace::nEdges() const noexcept
 {
     return 3;
 }
 
 
+inline Foam::edge Foam::triFace::faceEdge(const label edgei) const
+{
+    return Foam::edge(thisLabel(edgei), nextLabel(edgei));
+}
+
+
+inline Foam::edge Foam::triFace::edge(const label edgei) const
+{
+    return faceEdge(edgei);
+}
+
+
+inline Foam::vector Foam::triFace::edge
+(
+    const label edgei,
+    const UList<point>& pts
+) const
+{
+    return vector(pts[nextLabel(edgei)] - pts[thisLabel(edgei)]);
+}
+
+
+inline Foam::edge Foam::triFace::rcEdge(const label edgei) const
+{
+    // Edge 0 (forward and reverse) always starts at [0]
+    // for consistency with face flipping
+    const label pointi = edgei ? (3 - edgei) : 0;
+    return Foam::edge(thisLabel(pointi), prevLabel(pointi));
+}
+
+
+inline Foam::vector Foam::triFace::rcEdge
+(
+    const label edgei,
+    const UList<point>& pts
+) const
+{
+    // Edge 0 (forward and reverse) always starts at [0]
+    // for consistency with face flipping
+    const label pointi = edgei ? (3 - edgei) : 0;
+    return vector(pts[prevLabel(pointi)] - pts[thisLabel(pointi)]);
+}
+
+
 inline Foam::edgeList Foam::triFace::edges() const
 {
-    edgeList e(3);
+    edgeList theEdges(3);
 
-    e[0].first()  = operator[](0);
-    e[0].second() = operator[](1);
+    theEdges[0].first()  = operator[](0);
+    theEdges[0].second() = operator[](1);
 
-    e[1].first()  = operator[](1);
-    e[1].second() = operator[](2);
+    theEdges[1].first()  = operator[](1);
+    theEdges[1].second() = operator[](2);
 
-    e[2].first()  = operator[](2);
-    e[2].second() = operator[](0);
+    theEdges[2].first()  = operator[](2);
+    theEdges[2].second() = operator[](0);
 
-    return e;
+    return theEdges;
 }
 
 
-inline Foam::edge Foam::triFace::faceEdge(const label n) const
+inline Foam::edgeList Foam::triFace::rcEdges() const
 {
-    return edge(operator[](n), operator[](fcIndex(n)));
+    edgeList theEdges(3);
+
+    theEdges[0].first()  = operator[](0);
+    theEdges[0].second() = operator[](2);
+
+    theEdges[1].first()  = operator[](2);
+    theEdges[1].second() = operator[](1);
+
+    theEdges[2].first()  = operator[](1);
+    theEdges[2].second() = operator[](0);
+
+    return theEdges;
 }
 
 
-// return
-//  - +1: forward (counter-clockwise) on the face
-//  - -1: reverse (clockwise) on the face
-//  -  0: edge not found on the face
-inline int Foam::triFace::edgeDirection(const edge& e) const
+inline int Foam::triFace::edgeDirection(const Foam::edge& e) const
 {
-    if
-    (
-        (operator[](0) == e.first() && operator[](1) == e.second())
-     || (operator[](1) == e.first() && operator[](2) == e.second())
-     || (operator[](2) == e.first() && operator[](0) == e.second())
-    )
+    if (e.first() == operator[](0))
     {
-        return 1;
+        if (e.second() == operator[](1)) return 1;  // Forward
+        if (e.second() == operator[](2)) return -1; // Reverse
     }
-    else if
-    (
-        (operator[](0) == e.second() && operator[](1) == e.first())
-     || (operator[](1) == e.second() && operator[](2) == e.first())
-     || (operator[](2) == e.second() && operator[](0) == e.first())
-    )
+    if (e.first() == operator[](1))
     {
-        return -1;
+        if (e.second() == operator[](2)) return 1;  // Forward
+        if (e.second() == operator[](0)) return -1; // Reverse
+    }
+    if (e.first() == operator[](2))
+    {
+        if (e.second() == operator[](0)) return 1;  // Forward
+        if (e.second() == operator[](1)) return -1; // Reverse
     }
 
-    return 0;
+    return 0;  // Not found
 }
 
 
diff --git a/src/conversion/fire/checkFireEdges.C b/src/conversion/fire/checkFireEdges.C
index fe01fa5cc4..d19c2da1af 100644
--- a/src/conversion/fire/checkFireEdges.C
+++ b/src/conversion/fire/checkFireEdges.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2018 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -81,9 +81,9 @@ Foam::label Foam::checkFireEdges
     {
         const face& faceA = faces[faceI];
 
-        forAll(faceA, edgeI)
+        forAll(faceA, edgei)
         {
-            const edge currEdge = faceA.faceEdge(edgeI);
+            const edge currEdge = faceA.edge(edgei);
 
             // all faces attached to the first point
             const labelList& otherFaceIds = pointFaces[currEdge[0]];
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
index 85f5dada83..2e1b557c42 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,7 +37,7 @@ License
 #include "motionSmoother.H"
 #include "OFstream.H"
 
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
@@ -45,6 +45,8 @@ defineTypeNameAndDebug(edgeCollapser, 0);
 }
 
 
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
 Foam::labelHashSet Foam::edgeCollapser::checkBadFaces
 (
     const polyMesh& mesh,
@@ -415,7 +417,7 @@ void Foam::edgeCollapser::faceCollapseAxisAndAspectRatio
     {
         // It is possible that all the points of a face are the same
 
-        collapseAxis = f.edges()[f.longestEdge(pts)].unitVec(pts);
+        collapseAxis = f.edge(f.longestEdge(pts)).unitVec(pts);
 
         // Empirical correlation for high aspect ratio faces
 
@@ -432,7 +434,7 @@ void Foam::edgeCollapser::faceCollapseAxisAndAspectRatio
             // Cannot necessarily determine linearly independent
             // eigenvectors, or any at all, use longest edge direction.
 
-            collapseAxis = f.edges()[f.longestEdge(pts)].unitVec(pts);
+            collapseAxis = f.edge(f.longestEdge(pts)).unitVec(pts);
 
             aspectRatio = 1.0;
         }
diff --git a/src/finiteVolume/cfdTools/general/levelSet/levelSet.C b/src/finiteVolume/cfdTools/general/levelSet/levelSet.C
index 9a100b8726..4e800b71a0 100644
--- a/src/finiteVolume/cfdTools/general/levelSet/levelSet.C
+++ b/src/finiteVolume/cfdTools/general/levelSet/levelSet.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -122,9 +123,9 @@ Foam::tmp<Foam::scalarField> Foam::levelSetFraction
         vector a(Zero);
         vector r(Zero);
 
-        for (label eI = 0; eI < f.size(); ++eI)
+        for (label edgei = 0; edgei < f.nEdges(); ++edgei)
         {
-            const edge e = f.faceEdge(eI);
+            const edge e = f.edge(edgei);
 
             const FixedList<point, 3>
                 tri =
diff --git a/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C b/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C
index 943d7b87fe..f5e154f93b 100644
--- a/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C
+++ b/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -141,9 +142,9 @@ Foam::tmp<Foam::Field<Type>> Foam::levelSetAverage
         vector a(Zero);
         sumType r = Zero;
 
-        for (label eI = 0; eI < f.size(); ++eI)
+        for (label edgei = 0; edgei < f.nEdges(); ++edgei)
         {
-            const edge e = f.faceEdge(eI);
+            const edge e = f.edge(edgei);
 
             const FixedList<point, 3>
                 tri =
diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index a0eaddfc02..a71de211aa 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.C
@@ -297,7 +297,7 @@ void Foam::particle::changeFace(const label tetTriI)
         (
             ;
             edgeI < newFace.size()
-         && edge::compare(sharedEdge, newFace.faceEdge(edgeI)) != edgeComp;
+         && edge::compare(sharedEdge, newFace.edge(edgeI)) != edgeComp;
             ++ edgeI
         );
 
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C
index 7f7e547046..724c8000ee 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -84,7 +84,7 @@ Foam::label Foam::surfaceIntersection::getEdge
     const label fp
 )
 {
-    const edge faceEdge = surf.localFaces()[facei].faceEdge(fp);
+    const edge faceEdge = surf.localFaces()[facei].edge(fp);
 
     const labelList& eLabels = surf.faceEdges()[facei];
 
diff --git a/src/sampling/surface/cutting/cuttingSurfaceBaseTemplates.C b/src/sampling/surface/cutting/cuttingSurfaceBaseTemplates.C
index 5fd36c392b..9438c64eed 100644
--- a/src/sampling/surface/cutting/cuttingSurfaceBaseTemplates.C
+++ b/src/sampling/surface/cutting/cuttingSurfaceBaseTemplates.C
@@ -134,7 +134,7 @@ void Foam::cuttingSurfaceBase::walkCellCuts
 
             forAll(f, fp)
             {
-                edge e(f.faceEdge(fp));
+                edge e(f.edge(fp));
 
                 // Action #1: Orient edge (+ve gradient) and detect intersect
                 if (!edgeOrientIntersect(e))

From c0138ee8b67102e20ef92386de4e0c9fe0bcd053 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 4 May 2021 10:09:55 +0200
Subject: [PATCH 2/8] ENH: add forward/reverse circular index/value accessors
 for indirect lists

- improves interchangeability of List vs IndirectList
---
 .../IndirectListBase/IndirectListBase.H       | 70 ++++++++++--------
 .../IndirectListBase/IndirectListBaseI.H      | 74 ++++++++++++++++++-
 src/OpenFOAM/containers/Lists/UList/UList.H   | 12 +--
 src/OpenFOAM/containers/Lists/UList/UListI.H  | 22 +++---
 4 files changed, 127 insertions(+), 51 deletions(-)

diff --git a/src/OpenFOAM/containers/IndirectLists/IndirectListBase/IndirectListBase.H b/src/OpenFOAM/containers/IndirectLists/IndirectListBase/IndirectListBase.H
index cef22166d9..c9b01667ed 100644
--- a/src/OpenFOAM/containers/IndirectLists/IndirectListBase/IndirectListBase.H
+++ b/src/OpenFOAM/containers/IndirectLists/IndirectListBase/IndirectListBase.H
@@ -123,52 +123,25 @@ public:
     // Access
 
         //- The number of elements in the list
-        inline label size() const
+        inline label size() const noexcept
         {
             return addr_.size();
         }
 
         //- True if the list is empty (ie, size() is zero).
-        inline bool empty() const
+        inline bool empty() const noexcept
         {
             return addr_.empty();
         }
 
-        //- True if all entries have identical values, and list is non-empty
-        inline bool uniform() const;
-
-        //- The first element of the list.
-        inline T& first()
-        {
-            return values_[addr_.first()];
-        }
-
-        //- The first element of the list.
-        inline const T& first() const
-        {
-            return values_[addr_.first()];
-        }
-
-        //- The last element of the list.
-        inline T& last()
-        {
-            return values_[addr_.last()];
-        }
-
-        //- The last element of the list.
-        inline const T& last() const
-        {
-            return values_[addr_.last()];
-        }
-
         //- The list of values (without addressing)
-        inline const UList<T>& values() const
+        inline const UList<T>& values() const noexcept
         {
             return values_;
         }
 
         //- The list of values (without addressing)
-        inline UList<T>& values()
+        inline UList<T>& values() noexcept
         {
             return values_;
         }
@@ -179,6 +152,41 @@ public:
             return addr_;
         }
 
+        //- True if all entries have identical values, and list is non-empty
+        inline bool uniform() const;
+
+        //- The first element of the list.
+        inline const T& first() const;
+
+        //- The first element of the list.
+        inline T& first();
+
+        //- The last element of the list.
+        inline const T& last() const;
+
+        //- The last element of the list.
+        inline T& last();
+
+        //- The forward circular index. The next index in the list
+        //- which returns to the first at the end of the list
+        inline label fcIndex(const label i) const;
+
+        //- The reverse circular index. The previous index in the list
+        //- which returns to the last at the beginning of the list
+        inline label rcIndex(const label i) const;
+
+        //- Return forward circular value (ie, next value in the list)
+        inline const T& fcValue(const label i) const;
+
+        //- Return forward circular value (ie, next value in the list)
+        inline T& fcValue(const label i);
+
+        //- Return reverse circular value (ie, previous value in the list)
+        inline const T& rcValue(const label i) const;
+
+        //- Return reverse circular value (ie, previous value in the list)
+        inline T& rcValue(const label i);
+
 
     // Search
 
diff --git a/src/OpenFOAM/containers/IndirectLists/IndirectListBase/IndirectListBaseI.H b/src/OpenFOAM/containers/IndirectLists/IndirectListBase/IndirectListBaseI.H
index bc1aa07431..a6f87e688a 100644
--- a/src/OpenFOAM/containers/IndirectLists/IndirectListBase/IndirectListBaseI.H
+++ b/src/OpenFOAM/containers/IndirectLists/IndirectListBase/IndirectListBaseI.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -84,9 +84,9 @@ inline bool Foam::IndirectListBase<T, Addr>::uniform() const
         return false;
     }
 
-    const T& val = (*this)[0];
+    const T& val = (*this)[0];  // first
 
-    for (label i=1; i<len; ++i)
+    for (label i = 1; i < len; ++i)
     {
         if (val != (*this)[i])
         {
@@ -109,6 +109,74 @@ inline bool Foam::IndirectListBase<T, Addr>::found
 }
 
 
+template<class T, class Addr>
+inline Foam::label Foam::IndirectListBase<T, Addr>::fcIndex(const label i) const
+{
+    return (i == addr_.size()-1 ? 0 : i+1);
+}
+
+
+template<class T, class Addr>
+inline Foam::label Foam::IndirectListBase<T, Addr>::rcIndex(const label i) const
+{
+    return (i ? i-1 : addr_.size()-1);
+}
+
+
+template<class T, class Addr>
+inline const T& Foam::IndirectListBase<T, Addr>::first() const
+{
+    return values_[addr_.first()];
+}
+
+template<class T, class Addr>
+inline T& Foam::IndirectListBase<T, Addr>::first()
+{
+    return values_[addr_.first()];
+}
+
+
+template<class T, class Addr>
+inline const T& Foam::IndirectListBase<T, Addr>::last() const
+{
+    return values_[addr_.last()];
+}
+
+template<class T, class Addr>
+inline T& Foam::IndirectListBase<T, Addr>::last()
+{
+    return values_[addr_.last()];
+}
+
+
+template<class T, class Addr>
+inline const T& Foam::IndirectListBase<T, Addr>::fcValue(const label i) const
+{
+    return values_[this->fcIndex(i)];
+}
+
+
+template<class T, class Addr>
+inline T& Foam::IndirectListBase<T, Addr>::fcValue(const label i)
+{
+    return values_[this->fcIndex(i)];
+}
+
+
+template<class T, class Addr>
+inline const T& Foam::IndirectListBase<T, Addr>::rcValue(const label i) const
+{
+    return values_[this->rcIndex(i)];
+}
+
+
+template<class T, class Addr>
+inline T& Foam::IndirectListBase<T, Addr>::rcValue(const label i)
+{
+    return values_[this->rcIndex(i)];
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class T, class Addr>
diff --git a/src/OpenFOAM/containers/Lists/UList/UList.H b/src/OpenFOAM/containers/Lists/UList/UList.H
index 9df3e1f912..1493a88637 100644
--- a/src/OpenFOAM/containers/Lists/UList/UList.H
+++ b/src/OpenFOAM/containers/Lists/UList/UList.H
@@ -232,9 +232,13 @@ public:
 
         // Access
 
-            //- Return the forward circular index, i.e. next index
+            //- The forward circular index. The next index in the list
             //- which returns to the first at the end of the list
-            inline label fcIndex(const label i) const;
+            inline label fcIndex(const label i) const noexcept;
+
+            //- The reverse circular index. The previous index in the list
+            //- which returns to the last at the beginning of the list
+            inline label rcIndex(const label i) const noexcept;
 
             //- Return forward circular value (ie, next value in the list)
             inline const T& fcValue(const label i) const;
@@ -242,10 +246,6 @@ public:
             //- Return forward circular value (ie, next value in the list)
             inline T& fcValue(const label i);
 
-            //- Return the reverse circular index, i.e. previous index
-            //- which returns to the last at the beginning of the list
-            inline label rcIndex(const label i) const;
-
             //- Return reverse circular value (ie, previous value in the list)
             inline const T& rcValue(const label i) const;
 
diff --git a/src/OpenFOAM/containers/Lists/UList/UListI.H b/src/OpenFOAM/containers/Lists/UList/UListI.H
index 923e896f23..dc46db9faf 100644
--- a/src/OpenFOAM/containers/Lists/UList/UListI.H
+++ b/src/OpenFOAM/containers/Lists/UList/UListI.H
@@ -57,12 +57,19 @@ inline const Foam::UList<T>& Foam::UList<T>::null()
 
 
 template<class T>
-inline Foam::label Foam::UList<T>::fcIndex(const label i) const
+inline Foam::label Foam::UList<T>::fcIndex(const label i) const noexcept
 {
     return (i == size()-1 ? 0 : i+1);
 }
 
 
+template<class T>
+inline Foam::label Foam::UList<T>::rcIndex(const label i) const noexcept
+{
+    return (i ? i-1 : size()-1);
+}
+
+
 template<class T>
 inline const T& Foam::UList<T>::fcValue(const label i) const
 {
@@ -77,13 +84,6 @@ inline T& Foam::UList<T>::fcValue(const label i)
 }
 
 
-template<class T>
-inline Foam::label Foam::UList<T>::rcIndex(const label i) const
-{
-    return (i ? i-1 : size()-1);
-}
-
-
 template<class T>
 inline const T& Foam::UList<T>::rcValue(const label i) const
 {
@@ -146,14 +146,14 @@ inline bool Foam::UList<T>::uniform() const
 {
     const label len = size();
 
-    if (len == 0)
+    if (!len)
     {
         return false;
     }
 
-    const T& val = first();
+    const T& val = (*this)[0];  // first
 
-    for (label i=1; i<len; ++i)
+    for (label i = 1; i < len; ++i)
     {
         if (val != (*this)[i])
         {

From 40567b844a2037ce931ea8a5b552177af2123b4e Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 4 May 2021 14:25:28 +0200
Subject: [PATCH 3/8] DEFEATURE: remove construct List from two iterators
 (#2083)

- this constructor was added for similarity with std::vector,
  but continues to cause various annoyances.

  The main problem is that the templated parameter tends to grab
  anything that is not a perfect match for other constructors.

  Typically seen with two integers (in 64-bit mode), but various other
  cases as well.

  If required, the ListOps::create() function provides a lengthier
  alternative but one that can also incorporate transformations.
  Eg,

      pointField pts = ....;

      List<scalar> mags
      (
          List<scalar>::create
          (
              pts.begin(),
              pts.end(),
              [](const vector& v){ return magSqr(v); }
      );
---
 applications/test/List/Test-List.C               | 16 ++++++++++++++++
 .../test/zoneDistribute/Test-zoneDistribute.C    |  2 --
 .../containers/Lists/DynamicList/DynamicList.H   |  5 -----
 .../containers/Lists/DynamicList/DynamicListI.H  | 13 -------------
 src/OpenFOAM/containers/Lists/List/List.C        |  8 --------
 src/OpenFOAM/containers/Lists/List/List.H        |  5 -----
 6 files changed, 16 insertions(+), 33 deletions(-)

diff --git a/applications/test/List/Test-List.C b/applications/test/List/Test-List.C
index 16b3210958..4f52872b94 100644
--- a/applications/test/List/Test-List.C
+++ b/applications/test/List/Test-List.C
@@ -55,6 +55,9 @@ See also
 #include <numeric>
 #include <functional>
 
+// see issue #2083
+#undef Foam_constructList_from_iterators
+
 namespace Foam
 {
 
@@ -254,8 +257,12 @@ int main(int argc, char *argv[])
     Info<< "list4: " << list4 << nl
         << "list5: " << list5 << endl;
 
+    #ifdef Foam_constructList_from_iterators
     List<vector> list6(list4.begin(), list4.end());
     Info<< "list6: " << list6 << endl;
+    #else
+    Info<< "NOTE: no construction from two iterators" << endl;
+    #endif
 
     // Subset
     const labelList map{0, 2};
@@ -273,9 +280,13 @@ int main(int argc, char *argv[])
         // scalarList slist = identity(15);
         //
         // More writing, but does work:
+        #ifdef Foam_constructList_from_iterators
         scalarList slist(labelRange().begin(), labelRange(15).end());
 
         Info<<"scalar identity:" << flatOutput(slist) << endl;
+        #else
+        Info<<"No iterator means of creating a scalar identity list" << endl;
+        #endif
 
         printListOutputType<label>("labels") << nl;
 
@@ -384,6 +395,7 @@ int main(int argc, char *argv[])
         }
         Info<< "sub-sorted: " << flatOutput(longLabelList) << nl;
 
+        #ifdef Foam_constructList_from_iterators
         // Construct from a label-range
         labelRange range(25,15);
 
@@ -406,6 +418,10 @@ int main(int argc, char *argv[])
         // Even weird things like this
         List<scalar> sident4(labelRange().begin(), labelRange(8).end());
         Info<<"range-list (scalar)=" << sident4 << nl;
+        #else
+        Info<< "NOTE: no construction of labelList from range pair" << nl
+            << "use identity(...) instead" << endl;
+        #endif
     }
 
     wordReList reLst;
diff --git a/applications/test/zoneDistribute/Test-zoneDistribute.C b/applications/test/zoneDistribute/Test-zoneDistribute.C
index 00505735b6..87493b0eb4 100644
--- a/applications/test/zoneDistribute/Test-zoneDistribute.C
+++ b/applications/test/zoneDistribute/Test-zoneDistribute.C
@@ -38,8 +38,6 @@ Description
 #include "centredCPCCellToCellStencilObject.H"
 #include "zoneDistribute.H"
 
-#include "SortableList.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
diff --git a/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H b/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H
index e5e71f558a..378d54ab71 100644
--- a/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H
+++ b/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H
@@ -122,11 +122,6 @@ public:
         template<unsigned N>
         inline DynamicList(const FixedList<T, N>& lst);
 
-        //- Construct given begin/end iterators.
-        //  Uses std::distance to determine the size.
-        template<class InputIterator>
-        inline DynamicList(InputIterator begIter, InputIterator endIter);
-
         //- Construct from an initializer list. Size set to list size.
         inline explicit DynamicList(std::initializer_list<T> lst);
 
diff --git a/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H b/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
index 79da196a1e..911087b6e1 100644
--- a/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
+++ b/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
@@ -147,19 +147,6 @@ inline Foam::DynamicList<T, SizeMin>::DynamicList
 }
 
 
-template<class T, int SizeMin>
-template<class InputIterator>
-inline Foam::DynamicList<T, SizeMin>::DynamicList
-(
-    InputIterator begIter,
-    InputIterator endIter
-)
-:
-    List<T>(begIter, endIter),
-    capacity_(List<T>::size())
-{}
-
-
 template<class T, int SizeMin>
 inline Foam::DynamicList<T, SizeMin>::DynamicList
 (
diff --git a/src/OpenFOAM/containers/Lists/List/List.C b/src/OpenFOAM/containers/Lists/List/List.C
index db23ad23f2..ef0df3dbca 100644
--- a/src/OpenFOAM/containers/Lists/List/List.C
+++ b/src/OpenFOAM/containers/Lists/List/List.C
@@ -336,14 +336,6 @@ Foam::List<T>::List
 }
 
 
-template<class T>
-template<class InputIterator>
-Foam::List<T>::List(InputIterator begIter, InputIterator endIter)
-:
-    List<T>(begIter, endIter, std::distance(begIter, endIter))
-{}
-
-
 template<class T>
 template<unsigned N>
 Foam::List<T>::List(const FixedList<T, N>& list)
diff --git a/src/OpenFOAM/containers/Lists/List/List.H b/src/OpenFOAM/containers/Lists/List/List.H
index 68da584394..17f33459d4 100644
--- a/src/OpenFOAM/containers/Lists/List/List.H
+++ b/src/OpenFOAM/containers/Lists/List/List.H
@@ -160,11 +160,6 @@ public:
         template<unsigned N>
         List(const UList<T>& list, const FixedList<label, N>& indices);
 
-        //- Construct given begin/end iterators.
-        //  Uses std::distance for the size.
-        template<class InputIterator>
-        List(InputIterator begIter, InputIterator endIter);
-
         //- Construct as copy of FixedList\<T, N\>
         template<unsigned N>
         explicit List(const FixedList<T, N>& list);

From 17e6a147735d7288a3738c41ff8cfcf90eaffa9a Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 10 May 2021 10:02:15 +0200
Subject: [PATCH 4/8] ENH: improve support for assignment from indirect list

- copy assignment from indirect list to SubList and Field
---
 .../test/UIndirectList/Test-UIndirectList.C   | 18 ++++++++++--
 src/OpenFOAM/containers/Lists/List/List.H     |  4 ---
 .../containers/Lists/SubList/SubList.H        | 10 +++++--
 .../containers/Lists/SubList/SubListI.H       | 10 ++++++-
 src/OpenFOAM/containers/Lists/UList/UList.C   | 28 +++++++++++++++++--
 src/OpenFOAM/containers/Lists/UList/UList.H   | 11 +++++++-
 src/OpenFOAM/fields/Fields/Field/Field.H      | 13 ++++++---
 src/OpenFOAM/fields/Fields/Field/FieldI.H     | 13 ++++++++-
 8 files changed, 89 insertions(+), 18 deletions(-)

diff --git a/applications/test/UIndirectList/Test-UIndirectList.C b/applications/test/UIndirectList/Test-UIndirectList.C
index 9dbd0e8b0e..a9d9cabd76 100644
--- a/applications/test/UIndirectList/Test-UIndirectList.C
+++ b/applications/test/UIndirectList/Test-UIndirectList.C
@@ -62,7 +62,9 @@ void testFind(const T& val, const ListType& lst)
 
 int main(int argc, char *argv[])
 {
-    List<label> completeList(20);
+    labelList completeList(20);
+    labelList scratch(20);
+    scratch = -1;
 
     forAll(completeList, i)
     {
@@ -71,12 +73,24 @@ int main(int argc, char *argv[])
 
     Info<< "raw : " << flatOutput(completeList) << nl << endl;
 
-    List<label> addresses{1, 0, 3, 7, 4, 8, 5, 1, 0, 3, 7, 4, 8, 5, };
+    List<label> addresses({ 1, 0, 3, 7, 4, 8, 5, 1, 0, 3, 7, 4, 8, 5 });
 
     labelUIndList idl1(completeList, addresses);
 
     printInfo(idl1);
 
+    labelList::subList slice(scratch, idl1.size());
+    slice = idl1;
+
+    Info<< "sliced: " << flatOutput(slice) << nl;
+    Info<< "scratch: " << flatOutput(scratch) << nl;
+
+    // Again, but offset and using intermediate only
+    scratch = -1;
+    labelList::subList(scratch, idl1.size(), 5) = idl1;
+    Info<< "offset: " << flatOutput(scratch) << nl;
+
+
     for (const label val : { 10, 30, 40, 50, 90, 80, 120 } )
     {
         testFind(val, idl1);
diff --git a/src/OpenFOAM/containers/Lists/List/List.H b/src/OpenFOAM/containers/Lists/List/List.H
index 17f33459d4..ebded60c26 100644
--- a/src/OpenFOAM/containers/Lists/List/List.H
+++ b/src/OpenFOAM/containers/Lists/List/List.H
@@ -53,16 +53,12 @@ namespace Foam
 {
 
 // Forward Declarations
-class Istream;
-class Ostream;
-
 template<class T> class List;
 template<class T, unsigned N> class FixedList;
 template<class T, int SizeMin> class DynamicList;
 
 template<class T> class PtrList;
 template<class T> class SortableList;
-template<class T, class Addr> class IndirectListBase;
 
 template<class T> Istream& operator>>(Istream& is, List<T>& list);
 
diff --git a/src/OpenFOAM/containers/Lists/SubList/SubList.H b/src/OpenFOAM/containers/Lists/SubList/SubList.H
index 061028554c..4eddc1abb3 100644
--- a/src/OpenFOAM/containers/Lists/SubList/SubList.H
+++ b/src/OpenFOAM/containers/Lists/SubList/SubList.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -125,12 +125,16 @@ public:
         //- Allow cast to a const List\<T\>&
         inline operator const Foam::List<T>&() const;
 
-        //- Copy assign entries from given sub-list
+        //- Copy assign entries from given sub-list. Sizes must match!
         inline void operator=(const SubList<T>& list);
 
-        //- Copy assign entries from given list
+        //- Copy assign entries from given list. Sizes must match!
         inline void operator=(const UList<T>& list);
 
+        //- Copy assign entries from given indirect list. Sizes must match!
+        template<class Addr>
+        inline void operator=(const IndirectListBase<T, Addr>& list);
+
         //- Assign all entries to the given value
         inline void operator=(const T& val);
 
diff --git a/src/OpenFOAM/containers/Lists/SubList/SubListI.H b/src/OpenFOAM/containers/Lists/SubList/SubListI.H
index 3ec6aa5fa3..53408e3a76 100644
--- a/src/OpenFOAM/containers/Lists/SubList/SubListI.H
+++ b/src/OpenFOAM/containers/Lists/SubList/SubListI.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -135,6 +135,14 @@ inline void Foam::SubList<T>::operator=(const UList<T>& list)
 }
 
 
+template<class T>
+template<class Addr>
+inline void Foam::SubList<T>::operator=(const IndirectListBase<T, Addr>& list)
+{
+    UList<T>::deepCopy(list);
+}
+
+
 template<class T>
 inline void Foam::SubList<T>::operator=(const T& val)
 {
diff --git a/src/OpenFOAM/containers/Lists/UList/UList.C b/src/OpenFOAM/containers/Lists/UList/UList.C
index 41dbb8f036..3263d6f0c8 100644
--- a/src/OpenFOAM/containers/Lists/UList/UList.C
+++ b/src/OpenFOAM/containers/Lists/UList/UList.C
@@ -110,8 +110,8 @@ void Foam::UList<T>::deepCopy(const UList<T>& list)
     if (len != list.size_)
     {
         FatalErrorInFunction
-            << "ULists have different sizes: "
-            << len << " " << list.size_
+            << "Lists have different sizes: "
+            << len << " != " << list.size() << nl
             << abort(FatalError);
     }
     else if (len)
@@ -138,6 +138,30 @@ void Foam::UList<T>::deepCopy(const UList<T>& list)
 }
 
 
+template<class T>
+template<class Addr>
+void Foam::UList<T>::deepCopy(const IndirectListBase<T, Addr>& list)
+{
+    const label len = this->size_;
+
+    if (len != list.size())
+    {
+        FatalErrorInFunction
+            << "Lists have different sizes: "
+            << len << " != " << list.size() << nl
+            << abort(FatalError);
+    }
+    else if (len)
+    {
+        List_ACCESS(T, (*this), lhs);
+        for (label i = 0; i < len; ++i)
+        {
+            lhs[i] = list[i];
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class T>
diff --git a/src/OpenFOAM/containers/Lists/UList/UList.H b/src/OpenFOAM/containers/Lists/UList/UList.H
index 1493a88637..1fd97cb99e 100644
--- a/src/OpenFOAM/containers/Lists/UList/UList.H
+++ b/src/OpenFOAM/containers/Lists/UList/UList.H
@@ -68,10 +68,15 @@ namespace Foam
 {
 
 // Forward Declarations
+class Istream;
+class Ostream;
 class labelRange;
+
 template<class T> class List;
 template<class T> class SubList;
 template<class T> class UList;
+template<class T, class Addr> class IndirectListBase;
+
 template<class T> Istream& operator>>(Istream&, UList<T>&);
 template<class T> Ostream& operator<<(Ostream&, const UList<T>&);
 
@@ -342,6 +347,10 @@ public:
         //- Copy elements of the given UList
         void deepCopy(const UList<T>& list);
 
+        //- Copy elements of the given indirect list
+        template<class Addr>
+        void deepCopy(const IndirectListBase<T, Addr>& list);
+
 
     // Member Operators
 
@@ -493,7 +502,7 @@ public:
         typename std::enable_if<std::is_same<bool, TypeT>::value, bool>::type
         inline test(const label i) const
         {
-            return (i >= 0 && i < size() && v_[i]);
+            return (i >= 0 && i < size_ && v_[i]);
         }
 
         //- A bitSet::get() method for a list of bool
diff --git a/src/OpenFOAM/fields/Fields/Field/Field.H b/src/OpenFOAM/fields/Fields/Field/Field.H
index 916c0944a1..056265f3b4 100644
--- a/src/OpenFOAM/fields/Fields/Field/Field.H
+++ b/src/OpenFOAM/fields/Fields/Field/Field.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -101,7 +101,7 @@ public:
 
     // Constructors
 
-        //- Construct null
+        //- Default construct
         //  For temporary fields that are initialised after construction
         inline constexpr Field() noexcept;
 
@@ -374,8 +374,13 @@ public:
         //- Copy assignment
         void operator=(const Field<Type>&);
         void operator=(const tmp<Field<Type>>&);
-        inline void operator=(const UList<Type>&);
-        inline void operator=(const SubField<Type>&);
+
+        inline void operator=(const UList<Type>& rhs);
+        inline void operator=(const SubField<Type>& rhs);
+
+        //- Copy assign from IndirectList
+        template<class Addr>
+        inline void operator=(const IndirectListBase<Type, Addr>& rhs);
 
         //- Move assignment
         inline void operator=(Field<Type>&& rhs);
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldI.H b/src/OpenFOAM/fields/Fields/Field/FieldI.H
index 111c87115e..7cfd44d6f9 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldI.H
+++ b/src/OpenFOAM/fields/Fields/Field/FieldI.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -172,6 +172,17 @@ inline void Foam::Field<Type>::operator=(const SubField<Type>& rhs)
 }
 
 
+template<class Type>
+template<class Addr>
+inline void Foam::Field<Type>::operator=
+(
+    const IndirectListBase<Type, Addr>& rhs
+)
+{
+    List<Type>::operator=(rhs);
+}
+
+
 template<class Type>
 inline void Foam::Field<Type>::operator=(Field<Type>&& rhs)
 {

From ea2bf72740f587a10145383df0c18acaf980503c Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 10 May 2021 10:53:46 +0200
Subject: [PATCH 5/8] ENH: support gather of indirect list via globalIndex

---
 .../test/globalIndex/Test-globalIndex.C       |  32 +++-
 .../polyMesh/globalMeshData/globalIndex.H     |  53 ++++--
 .../globalMeshData/globalIndexTemplates.C     | 176 ++++++++++++------
 3 files changed, 187 insertions(+), 74 deletions(-)

diff --git a/applications/test/globalIndex/Test-globalIndex.C b/applications/test/globalIndex/Test-globalIndex.C
index bb436561b0..b7bfe0a3f2 100644
--- a/applications/test/globalIndex/Test-globalIndex.C
+++ b/applications/test/globalIndex/Test-globalIndex.C
@@ -28,14 +28,16 @@ Application
     globalIndexTest
 
 Description
-    Simple demonstration and test application for the globalIndex class.
+    Simple tests for the globalIndex class.
 
 \*---------------------------------------------------------------------------*/
 
 #include "globalIndex.H"
+#include "globalMeshData.H"
 #include "argList.H"
 #include "Time.H"
 #include "polyMesh.H"
+#include "IndirectList.H"
 #include "IOstreams.H"
 #include "Random.H"
 
@@ -218,6 +220,34 @@ int main(int argc, char *argv[])
         }
     }
 
+    Info<< "Gather indirect list of boundary points (patch = 0)\n";
+    {
+        const polyPatch& pp = mesh.boundaryMesh()[0];
+
+        // Map mesh point index to local (compact) point index
+        labelList pointToGlobal;
+        labelList uniqueMeshPointLabels;
+
+        autoPtr<globalIndex> globalPointsPtr =
+            mesh.globalData().mergePoints
+            (
+                pp.meshPoints(),
+                pp.meshPointMap(),
+                pointToGlobal,
+                uniqueMeshPointLabels
+            );
+
+        Info<< "local-sizes: " << globalPointsPtr().sizes() << nl;
+
+        UIndirectList<point> procPoints(mesh.points(), uniqueMeshPointLabels);
+        pointField patchPoints;
+
+        globalPointsPtr().gather(procPoints, patchPoints);
+
+        Info<< "gathered point field = " << patchPoints.size() << " points\n";
+    }
+
+    Info<< "\nEnd\n" << endl;
     return 0;
 }
 
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H
index 9059d2aa17..adb90aa2fd 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H
@@ -237,7 +237,7 @@ public:
         // Other
 
             //- Collect data in processor order on master (== procIDs[0]).
-            //  Needs offsets only on master.
+            //  Offsets needed on master only.
             template<class Container, class Type>
             static void gather
             (
@@ -247,12 +247,25 @@ public:
                 const UList<Type>& fld,
                 List<Type>& allFld,
                 const int tag = UPstream::msgType(),
-                const Pstream::commsTypes commsType =
-                    Pstream::commsTypes::nonBlocking
+                const Pstream::commsTypes = Pstream::commsTypes::nonBlocking
+            );
+
+            //- Collect indirect data in processor order on master
+            //  Offsets needed on master only.
+            template<class Type, class Addr>
+            static void gather
+            (
+                const labelUList& offsets,
+                const label comm,
+                const UList<int>& procIDs,
+                const IndirectListBase<Type, Addr>& fld,
+                List<Type>& allFld,
+                const int tag = UPstream::msgType(),
+                const Pstream::commsTypes = Pstream::commsTypes::scheduled
             );
 
             //- Collect data in processor order on master (== procIDs[0]).
-            //  Needs offsets only on master.
+            //  Offsets needed on master only.
             template<class Container, class Type>
             void gather
             (
@@ -276,8 +289,18 @@ public:
                 const UList<Type>& fld,
                 List<Type>& allFld,
                 const int tag = UPstream::msgType(),
-                const Pstream::commsTypes commsType =
-                    Pstream::commsTypes::nonBlocking
+                const Pstream::commsTypes = Pstream::commsTypes::nonBlocking
+            ) const;
+
+            //- Collect data indirectly in processor order on master.
+            //  Does communication with default communicator and message tag.
+            template<class Type, class Addr>
+            void gather
+            (
+                const IndirectListBase<Type, Addr>& fld,
+                List<Type>& allFld,
+                const int tag = UPstream::msgType(),
+                const Pstream::commsTypes = Pstream::commsTypes::scheduled
             ) const;
 
             //- Collect data in processor order on master.
@@ -288,8 +311,7 @@ public:
                 const UList<Type>& fld,
                 List<Type>& allFld,
                 const int tag = UPstream::msgType(),
-                const Pstream::commsTypes commsType =
-                    Pstream::commsTypes::nonBlocking
+                const Pstream::commsTypes = Pstream::commsTypes::nonBlocking
             );
 
 
@@ -303,8 +325,7 @@ public:
                 const Container& procIDs,
                 List<Type>& fld,
                 const int tag = UPstream::msgType(),
-                const Pstream::commsTypes commsType =
-                    Pstream::commsTypes::nonBlocking
+                const Pstream::commsTypes = Pstream::commsTypes::nonBlocking
             );
 
             //- Inplace collect in processor order on master (== procIDs[0]).
@@ -331,8 +352,7 @@ public:
             (
                 List<Type>& fld,
                 const int tag = UPstream::msgType(),
-                const Pstream::commsTypes commsType =
-                    Pstream::commsTypes::nonBlocking
+                const Pstream::commsTypes = Pstream::commsTypes::nonBlocking
             ) const;
 
             //- Inplace collect data in processor order on master
@@ -343,8 +363,7 @@ public:
             (
                 List<Type>& fld,
                 const int tag = UPstream::msgType(),
-                const Pstream::commsTypes commsType =
-                    Pstream::commsTypes::nonBlocking
+                const Pstream::commsTypes = Pstream::commsTypes::nonBlocking
             );
 
 
@@ -358,8 +377,7 @@ public:
                 const UList<Type>& allFld,
                 UList<Type>& fld,
                 const int tag = UPstream::msgType(),
-                const Pstream::commsTypes commsType =
-                    Pstream::commsTypes::nonBlocking
+                const Pstream::commsTypes = Pstream::commsTypes::nonBlocking
             );
 
             //- Distribute data in processor order. Requires fld to be sized!
@@ -386,8 +404,7 @@ public:
                 const UList<Type>& allFld,
                 UList<Type>& fld,
                 const int tag = UPstream::msgType(),
-                const Pstream::commsTypes commsType =
-                    Pstream::commsTypes::nonBlocking
+                const Pstream::commsTypes = Pstream::commsTypes::nonBlocking
             ) const;
 
             //- Get (potentially remote) data. Elements required given as
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexTemplates.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexTemplates.C
index 127d3688eb..9893418c00 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexTemplates.C
@@ -42,9 +42,21 @@ void Foam::globalIndex::gather
     const Pstream::commsTypes commsType
 )
 {
+    if
+    (
+        !is_contiguous<Type>::value
+     && commsType == Pstream::commsTypes::nonBlocking
+    )
+    {
+        FatalErrorInFunction
+            << "Cannot use nonBlocking with non-contiguous data"
+            << exit(FatalError);
+        // Could also warn and change to scheduled etc...
+    }
+
     if (Pstream::myProcNo(comm) == procIDs[0])
     {
-        allFld.setSize(off.last());
+        allFld.resize(off.last());
 
         // Assign my local data
         SubList<Type>(allFld, fld.size(), 0) = fld;
@@ -73,7 +85,7 @@ void Foam::globalIndex::gather
                 }
                 else
                 {
-                    IPstream fromSlave
+                    IPstream fromProc
                     (
                         commsType,
                         procIDs[i],
@@ -81,22 +93,15 @@ void Foam::globalIndex::gather
                         tag,
                         comm
                     );
-                    fromSlave >> procSlot;
+                    fromProc >> procSlot;
                 }
             }
         }
         else
         {
-            // nonBlocking
+            // nonBlocking && is_contiguous == true (already checked)
 
-            if (!is_contiguous<Type>::value)
-            {
-                FatalErrorInFunction
-                    << "nonBlocking not supported for non-contiguous data"
-                    << exit(FatalError);
-            }
-
-            label startOfRequests = Pstream::nRequests();
+            const label startOfRequests = Pstream::nRequests();
 
             // Set up reads
             for (label i = 1; i < procIDs.size(); ++i)
@@ -108,7 +113,7 @@ void Foam::globalIndex::gather
                     commsType,
                     procIDs[i],
                     reinterpret_cast<char*>(procSlot.data()),
-                    procSlot.byteSize(),
+                    procSlot.size_bytes(),
                     tag,
                     comm
                 );
@@ -153,16 +158,9 @@ void Foam::globalIndex::gather
         }
         else
         {
-            // nonBlocking
+            // nonBlocking && is_contiguous == true (already checked)
 
-            if (!is_contiguous<Type>::value)
-            {
-                FatalErrorInFunction
-                    << "nonBlocking not supported for non-contiguous data"
-                    << exit(FatalError);
-            }
-
-            label startOfRequests = Pstream::nRequests();
+            const label startOfRequests = Pstream::nRequests();
 
             // Set up write
             OPstream::write
@@ -182,6 +180,64 @@ void Foam::globalIndex::gather
 }
 
 
+template<class Type, class Addr>
+void Foam::globalIndex::gather
+(
+    const labelUList& off,
+    const label comm,
+    const UList<int>& procIDs,
+    const IndirectListBase<Type, Addr>& fld,
+    List<Type>& allFld,
+    const int tag,
+    const Pstream::commsTypes commsType
+)
+{
+    if (commsType == Pstream::commsTypes::nonBlocking)
+    {
+        WarningInFunction
+            << "Cannot use nonBlocking with indirect list of data"
+            << exit(FatalError);
+        // Could also warn and change to scheduled etc...
+    }
+
+    if (Pstream::myProcNo(comm) == procIDs[0])
+    {
+        allFld.resize(off.last());
+
+        // Assign my local data
+        SubList<Type>(allFld, fld.size(), 0) = fld;
+
+        // Already verified commsType != nonBlocking
+        for (label i = 1; i < procIDs.size(); ++i)
+        {
+            SubList<Type> procSlot(allFld, off[i+1]-off[i], off[i]);
+
+            IPstream fromProc
+            (
+                commsType,
+                procIDs[i],
+                0,
+                tag,
+                comm
+            );
+            fromProc >> procSlot;
+        }
+    }
+    else
+    {
+        OPstream toMaster
+        (
+            commsType,
+            procIDs[0],
+            0,
+            tag,
+            comm
+        );
+        toMaster << fld;
+    }
+}
+
+
 template<class Type>
 void Foam::globalIndex::gather
 (
@@ -203,6 +259,28 @@ void Foam::globalIndex::gather
 }
 
 
+template<class Type, class Addr>
+void Foam::globalIndex::gather
+(
+    const IndirectListBase<Type, Addr>& fld,
+    List<Type>& allFld,
+    const int tag,
+    const Pstream::commsTypes commsType
+) const
+{
+    gather
+    (
+        offsets_,
+        UPstream::worldComm,
+        UPstream::procID(UPstream::worldComm),
+        fld,
+        allFld,
+        tag,
+        commsType
+    );
+}
+
+
 template<class Type>
 void Foam::globalIndex::gatherOp
 (
@@ -293,6 +371,18 @@ void Foam::globalIndex::scatter
     const Pstream::commsTypes commsType
 )
 {
+    if
+    (
+        !is_contiguous<Type>::value
+     && commsType == Pstream::commsTypes::nonBlocking
+    )
+    {
+        FatalErrorInFunction
+            << "Cannot use nonBlocking with non-contiguous data"
+            << exit(FatalError);
+        // Could also warn and change to scheduled etc...
+    }
+
     if (Pstream::myProcNo(comm) == procIDs[0])
     {
         fld.deepCopy(SubList<Type>(allFld, off[1]-off[0]));
@@ -305,12 +395,7 @@ void Foam::globalIndex::scatter
         {
             for (label i = 1; i < procIDs.size(); ++i)
             {
-                const SubList<Type> procSlot
-                (
-                    allFld,
-                    off[i+1]-off[i],
-                    off[i]
-                );
+                const SubList<Type> procSlot(allFld, off[i+1]-off[i], off[i]);
 
                 if (is_contiguous<Type>::value)
                 {
@@ -326,7 +411,7 @@ void Foam::globalIndex::scatter
                 }
                 else
                 {
-                    OPstream toSlave
+                    OPstream toProc
                     (
                         commsType,
                         procIDs[i],
@@ -334,32 +419,20 @@ void Foam::globalIndex::scatter
                         tag,
                         comm
                     );
-                    toSlave << procSlot;
+                    toProc << procSlot;
                 }
             }
         }
         else
         {
-            // nonBlocking
+            // nonBlocking && is_contiguous == true (already checked)
 
-            if (!is_contiguous<Type>::value)
-            {
-                FatalErrorInFunction
-                    << "nonBlocking not supported for non-contiguous data"
-                    << exit(FatalError);
-            }
-
-            label startOfRequests = Pstream::nRequests();
+            const label startOfRequests = Pstream::nRequests();
 
             // Set up writes
             for (label i = 1; i < procIDs.size(); ++i)
             {
-                const SubList<Type> procSlot
-                (
-                    allFld,
-                    off[i+1]-off[i],
-                    off[i]
-                );
+                const SubList<Type> procSlot(allFld, off[i+1]-off[i], off[i]);
 
                 OPstream::write
                 (
@@ -411,16 +484,9 @@ void Foam::globalIndex::scatter
         }
         else
         {
-            // nonBlocking
+            // nonBlocking && is_contiguous == true (already checked)
 
-            if (!is_contiguous<Type>::value)
-            {
-                FatalErrorInFunction
-                    << "nonBlocking not supported for non-contiguous data"
-                    << exit(FatalError);
-            }
-
-            label startOfRequests = Pstream::nRequests();
+            const label startOfRequests = Pstream::nRequests();
 
             // Set up read
             IPstream::read
@@ -472,7 +538,7 @@ void Foam::globalIndex::get
     const int tag
 ) const
 {
-    allFld.setSize(globalIds.size());
+    allFld.resize(globalIds.size());
     if (globalIds.size())
     {
         // Sort according to processor

From 3bc8ab2839abe06d6ec646903ade1f23c72ad06b Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 10 May 2021 13:38:02 +0200
Subject: [PATCH 6/8] ENH: improve PrimitivePatch access for edges

- return internalEdges() and boundaryEdges() sub lists directly

- calculate and return boundaryFaces() to identify faces attached to
  boundary edges.

- minor code cleanup, and add PrimitivePatchBase class for
  non-templated code.

STYLE: mark unused parameter in globalMeshData mergePoints
---
 .../test/primitivePatch/Test-PrimitivePatch.C |  25 +++-
 src/OpenFOAM/Make/files                       |   2 +-
 .../polyMesh/globalMeshData/globalMeshData.C  |   2 +-
 .../polyMesh/globalMeshData/globalMeshData.H  |  20 +--
 .../PrimitivePatch/PrimitivePatch.C           |  39 ++++--
 .../PrimitivePatch/PrimitivePatch.H           |  51 ++++---
 .../PrimitivePatch/PrimitivePatchAddressing.C |  37 ++---
 ...mitivePatchName.C => PrimitivePatchBase.C} |   5 +-
 .../PrimitivePatch/PrimitivePatchBase.H       |  72 ++++++++++
 .../PrimitivePatch/PrimitivePatchBdryFaces.C  |  53 ++++++++
 .../PrimitivePatch/PrimitivePatchBdryPoints.C |  14 +-
 .../PrimitivePatch/PrimitivePatchEdgeLoops.C  |  16 +--
 .../PrimitivePatchLocalPointOrder.C           |  16 +--
 .../PrimitivePatch/PrimitivePatchMeshData.C   | 128 +++---------------
 .../PrimitivePatch/PrimitivePatchMeshEdges.C  |  77 +++++------
 .../PrimitivePatchPointAddressing.C           |  10 +-
 .../PrimitivePatchProjectPoints.C             |  51 +++----
 17 files changed, 322 insertions(+), 296 deletions(-)
 rename src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/{PrimitivePatchName.C => PrimitivePatchBase.C} (92%)
 create mode 100644 src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBase.H
 create mode 100644 src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryFaces.C

diff --git a/applications/test/primitivePatch/Test-PrimitivePatch.C b/applications/test/primitivePatch/Test-PrimitivePatch.C
index a74deb4bdf..e53f81b409 100644
--- a/applications/test/primitivePatch/Test-PrimitivePatch.C
+++ b/applications/test/primitivePatch/Test-PrimitivePatch.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,6 +34,7 @@ Description
 #include "polyMesh.H"
 #include "primitiveFacePatch.H"
 #include "primitivePatch.H"
+#include "IndirectList.H"
 #include "Fstream.H"
 
 using namespace Foam;
@@ -242,6 +243,28 @@ int main(int argc, char *argv[])
         writeEdgeFaces(localPoints, localFaces, edges, edgeFaces);
 
         writeFaceFaces(localPoints, localFaces, faceFaces);
+
+        const labelList bndFaceIds(pp.boundaryFaces());
+
+        Info<< "Have: " << bndFaceIds.size() << " boundary faces" << nl;
+
+        // Can calculate by hand
+        if (!pp.hasFaceCentres())
+        {
+            pointField boundaryCentres(bndFaceIds.size());
+
+            forAll(bndFaceIds, facei)
+            {
+                boundaryCentres[facei] =
+                    pp[bndFaceIds[facei]].centre(pp.points());
+            }
+            Info << "calc faceCentres:"
+                << boundaryCentres << nl;
+        }
+
+        // But will often use the cached information
+        Info<< "faceCentres:"
+            << UIndirectList<point>(pp.faceCentres(), bndFaceIds) << nl;
     }
 
     // Move construct
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 4ca6562acf..e7b0357b08 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -655,7 +655,7 @@ $(mapPolyMesh)/mapDistribute/IOmapDistribute.C
 $(mapPolyMesh)/mapAddedPolyMesh.C
 
 PrimitivePatch = $(primitiveMesh)/PrimitivePatch
-$(PrimitivePatch)/PrimitivePatchName.C
+$(PrimitivePatch)/PrimitivePatchBase.C
 
 pointMesh = meshes/pointMesh
 $(pointMesh)/pointMesh.C
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
index 3b98b0c996..9ff449fde5 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
@@ -2507,7 +2507,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
 Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
 (
     const labelList& meshPoints,
-    const Map<label>& meshPointMap,
+    const Map<label>& /* unused: meshPointMap */,
     labelList& pointToGlobal,
     labelList& uniqueMeshPoints
 ) const
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
index 2b90e7504e..8a58b914aa 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -342,7 +342,7 @@ public:
         // Access
 
             //- Return the mesh reference
-            const polyMesh& mesh() const
+            const polyMesh& mesh() const noexcept
             {
                 return mesh_;
             }
@@ -351,25 +351,25 @@ public:
             //  not running parallel)
             bool parallel() const
             {
-                return processorPatches_.size() > 0;
+                return !processorPatches_.empty();
             }
 
             //- Return total number of points in decomposed mesh. Not
             //  compensated for duplicate points!
-            label nTotalPoints() const
+            label nTotalPoints() const noexcept
             {
                 return nTotalPoints_;
             }
 
             //- Return total number of faces in decomposed mesh. Not
             //  compensated for duplicate faces!
-            label nTotalFaces() const
+            label nTotalFaces() const noexcept
             {
                 return nTotalFaces_;
             }
 
             //- Return total number of cells in decomposed mesh.
-            label nTotalCells() const
+            label nTotalCells() const noexcept
             {
                 return nTotalCells_;
             }
@@ -379,7 +379,7 @@ public:
 
             //- Return list of processor patch labels
             //  (size of list = number of processor patches)
-            const labelList& processorPatches() const
+            const labelList& processorPatches() const noexcept
             {
                 return processorPatches_;
             }
@@ -387,14 +387,14 @@ public:
             //- Return list of indices into processorPatches_ for each patch.
             //  Index = -1 for non-processor parches.
             //  (size of list = number of patches)
-            const labelList& processorPatchIndices() const
+            const labelList& processorPatchIndices() const noexcept
             {
                 return processorPatchIndices_;
             }
 
             //- Return processorPatchIndices of the neighbours
             //  processor patches. -1 if not running parallel.
-            const labelList& processorPatchNeighbours() const
+            const labelList& processorPatchNeighbours() const noexcept
             {
                 return processorPatchNeighbours_;
             }
@@ -577,7 +577,7 @@ public:
                 autoPtr<globalIndex> mergePoints
                 (
                     const labelList& meshPoints,
-                    const Map<label>& meshPointMap,
+                    const Map<label>& meshPointMap,  //!< currently unused
                     labelList& pointToGlobal,
                     labelList& uniqueMeshPoints
                 ) const;
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C
index d45ce84119..6c401211ad 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -130,7 +130,6 @@ Foam::PrimitivePatch<FaceList, PointField>::PrimitivePatch
     const PrimitivePatch<FaceList, PointField>& pp
 )
 :
-    PrimitivePatchName(),
     FaceList(pp),
     points_(pp.points_),
     edgesPtr_(nullptr),
@@ -173,13 +172,7 @@ Foam::PrimitivePatch<FaceList, PointField>::movePoints
     const Field<point_type>&
 )
 {
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-            << "movePoints() : "
-            << "recalculating PrimitivePatch geometry following mesh motion"
-            << endl;
-    }
+    DebugInFunction << "Recalculating geometry following mesh motion" << endl;
 
     clearGeom();
 }
@@ -198,6 +191,24 @@ Foam::PrimitivePatch<FaceList, PointField>::edges() const
 }
 
 
+template<class FaceList, class PointField>
+const Foam::edgeList::subList
+Foam::PrimitivePatch<FaceList, PointField>::internalEdges() const
+{
+    const edgeList& allEdges = this->edges();  // Force demand-driven
+    return edgeList::subList(allEdges, nInternalEdges());
+}
+
+
+template<class FaceList, class PointField>
+const Foam::edgeList::subList
+Foam::PrimitivePatch<FaceList, PointField>::boundaryEdges() const
+{
+    const edgeList& allEdges = this->edges();  // Force demand-driven
+    return edgeList::subList(allEdges, nBoundaryEdges(), nInternalEdges());
+}
+
+
 template<class FaceList, class PointField>
 Foam::label
 Foam::PrimitivePatch<FaceList, PointField>::nInternalEdges() const
@@ -211,6 +222,15 @@ Foam::PrimitivePatch<FaceList, PointField>::nInternalEdges() const
 }
 
 
+template<class FaceList, class PointField>
+Foam::label
+Foam::PrimitivePatch<FaceList, PointField>::nBoundaryEdges() const
+{
+    const edgeList& allEdges = this->edges();  // Force demand-driven
+    return (allEdges.size() - this->nInternalEdges());
+}
+
+
 template<class FaceList, class PointField>
 const Foam::labelList&
 Foam::PrimitivePatch<FaceList, PointField>::boundaryPoints() const
@@ -496,6 +516,7 @@ Foam::PrimitivePatch<FaceList, PointField>::operator=
 #include "PrimitivePatchAddressing.C"
 #include "PrimitivePatchEdgeLoops.C"
 #include "PrimitivePatchClear.C"
+#include "PrimitivePatchBdryFaces.C"
 #include "PrimitivePatchBdryPoints.C"
 #include "PrimitivePatchLocalPointOrder.C"
 #include "PrimitivePatchMeshData.C"
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
index f2cdc95af5..e05ed53e46 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,16 +36,16 @@ Description
     using List and pointField.
 
 SourceFiles
-    PrimitivePatchAddressing.C
-    PrimitivePatchBdryPoints.C
     PrimitivePatch.C
+    PrimitivePatchAddressing.C
+    PrimitivePatchBdryFaces.C
+    PrimitivePatchBdryPoints.C
     PrimitivePatchCheck.C
     PrimitivePatchClear.C
     PrimitivePatchEdgeLoops.C
     PrimitivePatchLocalPointOrder.C
     PrimitivePatchMeshData.C
     PrimitivePatchMeshEdges.C
-    PrimitivePatchName.C
     PrimitivePatchPointAddressing.C
     PrimitivePatchProjectPoints.C
 
@@ -61,6 +61,7 @@ SourceFiles
 #include "intersection.H"
 #include "HashSet.H"
 #include "objectHit.H"
+#include "PrimitivePatchBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -71,13 +72,6 @@ namespace Foam
 class face;
 template<class T> class Map;
 
-/*---------------------------------------------------------------------------*\
-                        Class PrimitivePatchName Declaration
-\*---------------------------------------------------------------------------*/
-
-TemplateName(PrimitivePatch);
-
-
 /*---------------------------------------------------------------------------*\
                            Class PrimitivePatch Declaration
 \*---------------------------------------------------------------------------*/
@@ -85,7 +79,7 @@ TemplateName(PrimitivePatch);
 template<class FaceList, class PointField>
 class PrimitivePatch
 :
-    public PrimitivePatchName,
+    public PrimitivePatchBase,
     public FaceList
 {
 public:
@@ -129,7 +123,7 @@ private:
         PointField points_;
 
 
-    // Demand driven private data
+    // Demand-driven Private Data
 
         //- Edges of the patch; address into local point list;
         //  sorted with internal edges first in upper-triangular order
@@ -194,10 +188,10 @@ private:
 
     // Private Member Functions
 
-        //- Calculate edges of the patch
-        void calcIntBdryEdges() const;
+        //- Calculate internal points on a patch
+        void calcInternPoints() const;
 
-        //- Calculated boundary points on a patch
+        //- Calculate boundary points on a patch
         void calcBdryPoints() const;
 
         //- Calculate addressing
@@ -308,18 +302,24 @@ public:
             return points_;
         }
 
+        //- Number of faces in the patch
+        label nFaces() const
+        {
+            return FaceList::size();
+        }
+
 
     // Access functions for demand-driven data
 
         // Topological data; no mesh required.
 
-            //- Return number of points supporting patch faces
+            //- Number of points supporting patch faces
             label nPoints() const
             {
                 return meshPoints().size();
             }
 
-            //- Return number of edges in patch
+            //- Number of edges in patch
             label nEdges() const
             {
                 return edges().size();
@@ -328,9 +328,18 @@ public:
             //- Return list of edges, address into LOCAL point list
             const edgeList& edges() const;
 
+            //- Return sub-list of internal edges, address into LOCAL point list
+            const edgeList::subList internalEdges() const;
+
+            //- Return sub-list of boundary edges, address into LOCAL point list
+            const edgeList::subList boundaryEdges() const;
+
             //- Number of internal edges
             label nInternalEdges() const;
 
+            //- Number of boundary edges == (nEdges() - nInternalEdges())
+            label nBoundaryEdges() const;
+
             //- Is internal edge?
             bool isInternalEdge(const label edgei) const
             {
@@ -358,6 +367,10 @@ public:
             //- Return patch faces addressing into local point list
             const List<face_type>& localFaces() const;
 
+            //- Calculate and return list of local boundary faces.
+            //  These are the faces attached to boundary edges.
+            labelList boundaryFaces() const;
+
 
         // Addressing into mesh
 
@@ -382,7 +395,7 @@ public:
 
             //- Given an edge in local point labels, return its
             //- index in the edge list.  If the edge is not found, return -1
-            label whichEdge(const edge&) const;
+            label whichEdge(const edge& e) const;
 
             //- Return labels of patch edges in the global edge list using
             //- cell addressing
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchAddressing.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchAddressing.C
index 7913307c60..1e9292768c 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchAddressing.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchAddressing.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -76,40 +76,29 @@ Foam::PrimitivePatch<FaceList, PointField>::calcAddressing() const
     edgeFacesPtr_.reset(new labelListList(maxEdges));
     auto& edgeFaces = *edgeFacesPtr_;
 
-    // faceFaces created using a dynamic list.  Cannot guess size because
-    // of multiple connections
-    List<DynamicList<label>> ff(locFcs.size());
-
     faceEdgesPtr_.reset(new labelListList(locFcs.size()));
     auto& faceEdges = *faceEdgesPtr_;
 
-    // count the number of face neighbours
-    labelList noFaceFaces(locFcs.size());
-
     // initialise the lists of subshapes for each face to avoid duplication
     edgeListList faceIntoEdges(locFcs.size());
 
     forAll(locFcs, facei)
     {
         faceIntoEdges[facei] = locFcs[facei].edges();
-
-        labelList& curFaceEdges = faceEdges[facei];
-        curFaceEdges.setSize(faceIntoEdges[facei].size());
-
-        forAll(curFaceEdges, faceEdgeI)
-        {
-            curFaceEdges[faceEdgeI] = -1;
-        }
+        faceEdges[facei].resize(faceIntoEdges[facei].size(), -1);
     }
 
+    // faceFaces created using a dynamic list.  Cannot guess size because
+    // of multiple connections
+    List<DynamicList<label>> ff(locFcs.size());
+
+
     // This algorithm will produce a separated list of edges, internal edges
     // starting from 0 and boundary edges starting from the top and
     // growing down.
 
     label nEdges = 0;
 
-    bool found = false;
-
     // Note that faceIntoEdges is sorted acc. to local vertex numbering
     // in face (i.e. curEdges[0] is edge between f[0] and f[1])
 
@@ -132,8 +121,6 @@ Foam::PrimitivePatch<FaceList, PointField>::calcAddressing() const
             // If the edge is already detected, skip
             if (faceEdges[facei][edgeI] >= 0) continue;
 
-            found = false;
-
             // Set reference to the current edge
             const edge& e = curEdges[edgeI];
 
@@ -141,6 +128,8 @@ Foam::PrimitivePatch<FaceList, PointField>::calcAddressing() const
 
             const labelList& nbrFaces = pf[e.start()];
 
+            bool found = false;
+
             forAll(nbrFaces, nbrFacei)
             {
                 // set reference to the current neighbour
@@ -212,7 +201,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcAddressing() const
 
                 // Set edge-face addressing
                 labelList& curEf = edgeFaces[nEdges];
-                curEf.setSize(cnf.size() + 1);
+                curEf.resize(cnf.size() + 1);
                 curEf[0] = facei;
 
                 forAll(cnf, cnfI)
@@ -256,7 +245,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcAddressing() const
 
                 // Add edgeFace
                 labelList& curEf = edgeFaces[nEdges];
-                curEf.setSize(1);
+                curEf.resize(1);
                 curEf[0] = facei;
 
                 nEdges++;
@@ -265,10 +254,10 @@ Foam::PrimitivePatch<FaceList, PointField>::calcAddressing() const
     }
 
     // edges
-    edges.setSize(nEdges);
+    edges.resize(nEdges);
 
     // edgeFaces list
-    edgeFaces.setSize(nEdges);
+    edgeFaces.resize(nEdges);
 
     // faceFaces list
     faceFacesPtr_.reset(new labelListList(locFcs.size()));
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchName.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBase.C
similarity index 92%
rename from src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchName.C
rename to src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBase.C
index e3ffb769a5..822579bb16 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchName.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBase.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,13 +26,13 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "PrimitivePatch.H"
+#include "PrimitivePatchBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-defineTypeNameAndDebug(PrimitivePatchName, 0);
+defineTypeNameAndDebug(PrimitivePatchBase, 0);
 }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBase.H b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBase.H
new file mode 100644
index 0000000000..2a2e0d0c96
--- /dev/null
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBase.H
@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PrimitivePatchBase
+
+Description
+    Non-templated base elements for PrimitivePatch
+
+SourceFiles
+    PrimitivePatchBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PrimitivePatchBase_H
+#define PrimitivePatchBase_H
+
+#include "className.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class PrimitivePatchBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class PrimitivePatchBase
+{
+public:
+
+    //- Runtime type information
+    ClassName("PrimitivePatch");
+
+
+    // Constructors
+
+        //- Default construct
+        PrimitivePatchBase() = default;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryFaces.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryFaces.C
new file mode 100644
index 0000000000..61bba9d782
--- /dev/null
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryFaces.C
@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PrimitivePatch.H"
+#include "HashSet.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class FaceList, class PointField>
+Foam::labelList
+Foam::PrimitivePatch<FaceList, PointField>::boundaryFaces() const
+{
+    // By definition boundary edges have a _single_ face attached,
+    // but a face can easily have multiple boundary edges.
+
+    const labelListList& edgeToFace = edgeFaces();
+
+    labelHashSet bndFaces(2*nBoundaryEdges());
+
+    for (label edgei = nInternalEdges(); edgei < edgeToFace.size(); ++edgei)
+    {
+        bndFaces.insert(edgeToFace[edgei][0]);
+    }
+
+    return bndFaces.sortedToc();
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C
index 3ced498e62..5b58bcf49e 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,16 +45,12 @@ Foam::PrimitivePatch<FaceList, PointField>::calcBdryPoints() const
             << abort(FatalError);
     }
 
-    const edgeList& e = edges();
+    labelHashSet bp(2*nEdges());
 
-    labelHashSet bp(2*e.size());
-
-    for (label edgeI = nInternalEdges_; edgeI < e.size(); edgeI++)
+    for (const edge& e : boundaryEdges())
     {
-        const edge& curEdge = e[edgeI];
-
-        bp.insert(curEdge.start());
-        bp.insert(curEdge.end());
+        bp.insert(e.first());
+        bp.insert(e.second());
     }
 
     boundaryPointsPtr_.reset(new labelList(bp.sortedToc()));
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchEdgeLoops.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchEdgeLoops.C
index 90d9f75af5..8b19fb5320 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchEdgeLoops.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchEdgeLoops.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,12 +49,15 @@ Foam::PrimitivePatch<FaceList, PointField>::calcEdgeLoops() const
     }
 
     const edgeList& patchEdges = edges();
-    label nIntEdges = nInternalEdges();
-    label nBdryEdges = patchEdges.size() - nIntEdges;
+    const label nIntEdges = nInternalEdges();
+    const label nBdryEdges = patchEdges.size() - nIntEdges;
+
+    // Size return list plenty big
+    edgeLoopsPtr_.reset(new labelListList(nBdryEdges));
+    auto& edgeLoops = *edgeLoopsPtr_;
 
     if (nBdryEdges == 0)
     {
-        edgeLoopsPtr_.reset(new labelListList(0));
         return;
     }
 
@@ -68,11 +71,6 @@ Foam::PrimitivePatch<FaceList, PointField>::calcEdgeLoops() const
     // Loop per (boundary) edge.
     labelList loopNumber(nBdryEdges, -1);
 
-    // Size return list plenty big
-    edgeLoopsPtr_.reset(new labelListList(nBdryEdges));
-    auto& edgeLoops = *edgeLoopsPtr_;
-
-
     // Current loop number.
     label loopI = 0;
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchLocalPointOrder.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchLocalPointOrder.C
index abe417f13e..d687718b92 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchLocalPointOrder.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchLocalPointOrder.C
@@ -41,13 +41,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcLocalPointOrder() const
     // Note: Cannot use bandCompressing as point-point addressing does
     // not exist and is not considered generally useful.
 
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-            << "calcLocalPointOrder() : "
-            << "calculating local point order"
-            << endl;
-    }
+    DebugInFunction << "Calculating local point order" << endl;
 
     if (localPointOrderPtr_)
     {
@@ -116,13 +110,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcLocalPointOrder() const
         }
     }
 
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-            << "calcLocalPointOrder() "
-            << "finished calculating local point order"
-            << endl;
-    }
+    DebugInfo << "Calculated local point order" << endl;
 }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C
index df1e23e511..a58efd3385 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C
@@ -35,13 +35,7 @@ template<class FaceList, class PointField>
 void
 Foam::PrimitivePatch<FaceList, PointField>::calcMeshData() const
 {
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcMeshData() : "
-               "calculating mesh data in PrimitivePatch"
-            << endl;
-    }
+    DebugInFunction << "Calculating mesh data" << endl;
 
     if (meshPointsPtr_ || localFacesPtr_)
     {
@@ -117,13 +111,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcMeshData() const
         }
     }
 
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcMeshData() : "
-               "finished calculating mesh data in PrimitivePatch"
-            << endl;
-    }
+    DebugInfo << "Calculated mesh data" << endl;
 }
 
 
@@ -131,13 +119,7 @@ template<class FaceList, class PointField>
 void
 Foam::PrimitivePatch<FaceList, PointField>::calcMeshPointMap() const
 {
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcMeshPointMap() : "
-               "calculating mesh point map in PrimitivePatch"
-            << endl;
-    }
+    DebugInFunction << "Calculating mesh point map" << endl;
 
     if (meshPointMapPtr_)
     {
@@ -157,13 +139,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcMeshPointMap() const
         mpMap.insert(mp[i], i);
     }
 
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcMeshPointMap() : "
-               "finished calculating mesh point map in PrimitivePatch"
-            << endl;
-    }
+    DebugInfo << "Calculated mesh point map" << endl;
 }
 
 
@@ -171,13 +147,7 @@ template<class FaceList, class PointField>
 void
 Foam::PrimitivePatch<FaceList, PointField>::calcLocalPoints() const
 {
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcLocalPoints() : "
-               "calculating localPoints in PrimitivePatch"
-            << endl;
-    }
+    DebugInFunction << "Calculating localPoints" << endl;
 
     if (localPointsPtr_)
     {
@@ -197,13 +167,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcLocalPoints() const
         locPts[pointi] = points_[meshPts[pointi]];
     }
 
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-            << "calcLocalPoints() : "
-            << "finished calculating localPoints in PrimitivePatch"
-            << endl;
-    }
+    DebugInfo << "Calculated localPoints" << endl;
 }
 
 
@@ -211,13 +175,7 @@ template<class FaceList, class PointField>
 void
 Foam::PrimitivePatch<FaceList, PointField>::calcPointNormals() const
 {
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcPointNormals() : "
-               "calculating pointNormals in PrimitivePatch"
-            << endl;
-    }
+    DebugInFunction << "Calculating pointNormals" << endl;
 
     if (pointNormalsPtr_)
     {
@@ -248,13 +206,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcPointNormals() const
         curNormal.normalise();
     }
 
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcPointNormals() : "
-               "finished calculating pointNormals in PrimitivePatch"
-            << endl;
-    }
+    DebugInfo << "Calculated pointNormals" << endl;
 }
 
 
@@ -262,13 +214,7 @@ template<class FaceList, class PointField>
 void
 Foam::PrimitivePatch<FaceList, PointField>::calcFaceCentres() const
 {
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcFaceCentres() : "
-               "calculating faceCentres in PrimitivePatch"
-            << endl;
-    }
+    DebugInFunction << "Calculating faceCentres" << endl;
 
     if (faceCentresPtr_)
     {
@@ -286,13 +232,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcFaceCentres() const
         c[facei] = this->operator[](facei).centre(points_);
     }
 
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcFaceCentres() : "
-               "finished calculating faceCentres in PrimitivePatch"
-            << endl;
-    }
+    DebugInfo << "Calculated faceCentres" << endl;
 }
 
 
@@ -300,13 +240,7 @@ template<class FaceList, class PointField>
 void
 Foam::PrimitivePatch<FaceList, PointField>::calcMagFaceAreas() const
 {
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcMagFaceAreas() : "
-               "calculating magFaceAreas in PrimitivePatch"
-            << endl;
-    }
+    DebugInFunction << "Calculating magFaceAreas" << endl;
 
     if (magFaceAreasPtr_)
     {
@@ -324,13 +258,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcMagFaceAreas() const
         a[facei] = this->operator[](facei).mag(points_);
     }
 
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcMagFaceAreas() : "
-               "finished calculating magFaceAreas in PrimitivePatch"
-            << endl;
-    }
+    DebugInfo << "Calculated magFaceAreas" << endl;
 }
 
 
@@ -338,13 +266,7 @@ template<class FaceList, class PointField>
 void
 Foam::PrimitivePatch<FaceList, PointField>::calcFaceAreas() const
 {
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcFaceAreas() : "
-               "calculating faceAreas in PrimitivePatch"
-            << endl;
-    }
+    DebugInFunction << "Calculating faceAreas" << endl;
 
     if (faceAreasPtr_)
     {
@@ -362,13 +284,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcFaceAreas() const
         n[facei] = this->operator[](facei).areaNormal(points_);
     }
 
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcFaceAreas() : "
-               "finished calculating faceAreas in PrimitivePatch"
-            << endl;
-    }
+    DebugInfo << "Calculated faceAreas" << endl;
 }
 
 
@@ -376,13 +292,7 @@ template<class FaceList, class PointField>
 void
 Foam::PrimitivePatch<FaceList, PointField>::calcFaceNormals() const
 {
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcFaceNormals() : "
-               "calculating faceNormals in PrimitivePatch"
-            << endl;
-    }
+    DebugInFunction << "Calculating faceNormals" << endl;
 
     if (faceNormalsPtr_)
     {
@@ -400,13 +310,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcFaceNormals() const
         n[facei] = this->operator[](facei).unitNormal(points_);
     }
 
-    if (debug)
-    {
-        Pout<< "PrimitivePatch<FaceList, PointField>::"
-               "calcFaceNormals() : "
-               "finished calculating faceNormals in PrimitivePatch"
-            << endl;
-    }
+    DebugInfo << "Calculated faceNormals" << endl;
 }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshEdges.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshEdges.C
index 9af7dcea2a..cf1a5af2fd 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshEdges.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshEdges.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,47 +43,38 @@ meshEdges
     DebugInFunction
         << "Calculating labels of patch edges in mesh edge list" << nl;
 
-    // get reference to the list of edges on the patch
+    // The list of edges on the patch
     const edgeList& PatchEdges = edges();
 
+    // The output storage
+    labelList meshEdges(PatchEdges.size());
+
     const labelListList& EdgeFaces = edgeFaces();
 
-    // create the storage
-    labelList meshEdges(PatchEdges.size());
-
-    bool found = false;
-
-    // get reference to the points on the patch
+    // The mesh points associated with the patch
     const labelList& pp = meshPoints();
 
     // WARNING: Remember that local edges address into local point list;
     // local-to-global point label translation is necessary
-    forAll(PatchEdges, edgeI)
+    forAll(PatchEdges, edgei)
     {
-        const edge curEdge
-            (pp[PatchEdges[edgeI].start()], pp[PatchEdges[edgeI].end()]);
+        bool found = false;
 
-        found = false;
+        const edge globalEdge(pp, PatchEdges[edgei]);
 
-        // get the patch faces sharing the edge
-        const labelList& curFaces = EdgeFaces[edgeI];
-
-        forAll(curFaces, facei)
+        // For each patch face sharing the edge
+        for (const label patchFacei : EdgeFaces[edgei])
         {
-            // get the cell next to the face
-            label curCell = faceCells[curFaces[facei]];
+            // The cell next to the face
+            const label curCelli = faceCells[patchFacei];
 
-            // get reference to edges on the cell
-            const labelList& ce = cellEdges[curCell];
-
-            forAll(ce, cellEdgeI)
+            // Check the cell edges
+            for (const label cellEdgei : cellEdges[curCelli])
             {
-                if (allEdges[ce[cellEdgeI]] == curEdge)
+                if (allEdges[cellEdgei] == globalEdge)
                 {
                     found = true;
-
-                    meshEdges[edgeI] = ce[cellEdgeI];
-
+                    meshEdges[edgei] = cellEdgei;
                     break;
                 }
             }
@@ -107,29 +98,27 @@ Foam::PrimitivePatch<FaceList, PointField>::meshEdges
     DebugInFunction
         << "Calculating labels of patch edges in mesh edge list" << nl;
 
-    // get reference to the list of edges on the patch
+    // The list of edges on the patch
     const edgeList& PatchEdges = edges();
 
-    // create the storage
+    // The output storage
     labelList meshEdges(PatchEdges.size());
 
-    // get reference to the points on the patch
+    // The mesh points associated with the patch
     const labelList& pp = meshPoints();
 
     // WARNING: Remember that local edges address into local point list;
     // local-to-global point label translation is necessary
-    forAll(PatchEdges, edgeI)
+    forAll(PatchEdges, edgei)
     {
-        const label globalPointi = pp[PatchEdges[edgeI].start()];
-        const edge curEdge(globalPointi, pp[PatchEdges[edgeI].end()]);
+        const edge globalEdge(pp, PatchEdges[edgei]);
 
-        const labelList& pe = pointEdges[globalPointi];
-
-        forAll(pe, i)
+        // Check the attached edges
+        for (const label patchEdgei : pointEdges[globalEdge.start()])
         {
-            if (allEdges[pe[i]] == curEdge)
+            if (allEdges[patchEdgei] == globalEdge)
             {
-                meshEdges[edgeI] = pe[i];
+                meshEdges[edgei] = patchEdgei;
                 break;
             }
         }
@@ -148,18 +137,16 @@ Foam::PrimitivePatch<FaceList, PointField>::whichEdge
     const edge& e
 ) const
 {
-    // Get pointEdges from the starting point and search all the candidates
-    const edgeList& Edges = edges();
-
-    if (e.start() > -1 && e.start() < nPoints())
+    if (e.start() >= 0 && e.start() < nPoints())
     {
-        const labelList& pe = pointEdges()[e.start()];
+        // Get pointEdges from the starting point and search all the candidates
+        const edgeList& Edges = edges();
 
-        forAll(pe, peI)
+        for (const label patchEdgei : pointEdges()[e.start()])
         {
-            if (e == Edges[pe[peI]])
+            if (e == Edges[patchEdgei])
             {
-                return pe[peI];
+                return patchEdgei;
             }
         }
     }
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchPointAddressing.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchPointAddressing.C
index d564e3e939..20b72272ea 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchPointAddressing.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchPointAddressing.C
@@ -54,8 +54,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcPointEdges() const
 
     invertManyToMany(pe.size(), edges(), pe);
 
-    DebugInfo
-        << "    Finished." << endl;
+    DebugInfo << "    Finished." << endl;
 }
 
 
@@ -81,9 +80,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcPointFaces() const
 
     forAll(locFcs, facei)
     {
-        const face_type& curPoints = locFcs[facei];
-
-        for (const label pointi : curPoints)
+        for (const label pointi : locFcs[facei])
         {
             pointFcs[pointi].append(facei);
         }
@@ -98,8 +95,7 @@ Foam::PrimitivePatch<FaceList, PointField>::calcPointFaces() const
         pf[pointi] = pointFcs[pointi];
     }
 
-    DebugInfo
-        << "    Finished." << endl;
+    DebugInfo << "    Finished." << endl;
 }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchProjectPoints.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchProjectPoints.C
index 886a606200..7ca7f5d079 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchProjectPoints.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchProjectPoints.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -200,7 +200,7 @@ Foam::PrimitivePatch<FaceList, PointField>::projectPoints
                     }
                 }
 
-                if (debug) Info<< ".";
+                DebugInfo << '.';
             } while (closer);
         }
 
@@ -211,10 +211,7 @@ Foam::PrimitivePatch<FaceList, PointField>::projectPoints
         {
             nNSquaredSearches++;
 
-            if (debug)
-            {
-                Info<< "p " << curLocalPointLabel << ": ";
-            }
+            DebugInfo << "p " << curLocalPointLabel << ": ";
 
             result[curLocalPointLabel] = objectHit(false, -1);
             scalar minDistance = GREAT;
@@ -254,23 +251,18 @@ Foam::PrimitivePatch<FaceList, PointField>::projectPoints
                 }
             }
 
-            if (debug)
-            {
-                Info<< result[curLocalPointLabel] << nl;
-            }
+            DebugInfo << result[curLocalPointLabel] << nl;
         }
         else
         {
-            if (debug) Info<< "x";
+            DebugInfo << 'x';
         }
     }
 
-    if (debug)
-    {
-        Info<< nl << "Executed " << nNSquaredSearches
-            << " n-squared searches out of total of "
-            << nPoints() << endl;
-    }
+    DebugInfo
+        << nl << "Executed " << nNSquaredSearches
+        << " n-squared searches out of total of "
+        << nPoints() << endl;
 
     return result;
 }
@@ -439,7 +431,7 @@ Foam::PrimitivePatch<FaceList, PointField>::projectFaceCentres
                     }
                 }
 
-                if (debug) Info<< ".";
+                DebugInfo << '.';
             } while (closer);
         }
 
@@ -447,10 +439,7 @@ Foam::PrimitivePatch<FaceList, PointField>::projectFaceCentres
         {
             nNSquaredSearches++;
 
-            if (debug)
-            {
-                Info<< "p " << curLocalFaceLabel << ": ";
-            }
+            DebugInfo << "p " << curLocalFaceLabel << ": ";
 
             result[curLocalFaceLabel] = objectHit(false, -1);
             scalar minDistance = GREAT;
@@ -490,23 +479,19 @@ Foam::PrimitivePatch<FaceList, PointField>::projectFaceCentres
                 }
             }
 
-            if (debug)
-            {
-                Info<< result[curLocalFaceLabel] << nl;
-            }
+            DebugInfo << result[curLocalFaceLabel] << nl;
         }
         else
         {
-            if (debug) Info<< "x";
+            DebugInfo << 'x';
         }
     }
 
-    if (debug)
-    {
-        Info<< nl << "Executed " << nNSquaredSearches
-            << " n-squared searches out of total of "
-            << this->size() << endl;
-    }
+    DebugInfo
+        << nl
+        << "Executed " << nNSquaredSearches
+        << " n-squared searches out of total of "
+        << this->size() << endl;
 
     return result;
 }

From 2a438385a313ffb12ad7bb9d0ceb0b78b579f490 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 11 May 2021 08:48:18 +0200
Subject: [PATCH 7/8] STYLE: make polyMesh constructor explicit

- use polyMesh::defaultRegion instead of fvMesh::defaultRegion and
  others via static inheritance. (92 vs ~25 occurrences)
---
 .../overRhoSimpleFoam/createUpdatedDynamicFvMesh.H            | 2 +-
 .../simpleFoam/overSimpleFoam/createUpdatedDynamicFvMesh.H    | 2 +-
 applications/test/mesh/Test-mesh.C                            | 2 +-
 .../utilities/mesh/advanced/refinementLevel/refinementLevel.C | 2 +-
 .../utilities/mesh/conversion/foamMeshToFluent/fluentFvMesh.H | 2 +-
 .../mesh/conversion/foamMeshToFluent/foamMeshToFluent.C       | 2 +-
 .../mesh/generation/extrude/extrudeMesh/extrudeMesh.C         | 4 ++--
 .../utilities/mesh/manipulation/mirrorMesh/mirrorMesh.C       | 2 +-
 .../parallelProcessing/redistributePar/redistributePar.C      | 2 +-
 .../preProcessing/changeDictionary/changeDictionary.C         | 2 +-
 .../createExternalCoupledPatchGeometry.C                      | 2 +-
 applications/utilities/preProcessing/mapFields/mapFields.C    | 4 ++--
 .../utilities/preProcessing/mapFieldsPar/mapFieldsPar.C       | 4 ++--
 src/OpenFOAM/include/createMeshNoClear.H                      | 2 +-
 src/OpenFOAM/meshes/polyMesh/polyMesh.H                       | 4 ++--
 src/dynamicFvMesh/dynamicFvMesh/dynamicFvMeshNew.C            | 4 ++--
 src/dynamicFvMesh/include/createNamedDynamicFvMesh.H          | 2 +-
 src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C           | 2 +-
 src/engine/include/createEngineMesh.H                         | 2 +-
 .../fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.C       | 4 ++--
 .../fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.H       | 2 +-
 src/functionObjects/field/streamLine/streamLineBase.C         | 2 +-
 src/sampling/sampledSet/sampledSets/sampledSets.C             | 4 ++--
 23 files changed, 30 insertions(+), 30 deletions(-)

diff --git a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createUpdatedDynamicFvMesh.H b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createUpdatedDynamicFvMesh.H
index db58f90828..ea5823cc1c 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createUpdatedDynamicFvMesh.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createUpdatedDynamicFvMesh.H
@@ -7,7 +7,7 @@
         (
             IOobject
             (
-                dynamicFvMesh::defaultRegion,
+                polyMesh::defaultRegion,
                 runTime.timeName(),
                 runTime,
                 IOobject::MUST_READ
diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createUpdatedDynamicFvMesh.H b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createUpdatedDynamicFvMesh.H
index db58f90828..ea5823cc1c 100644
--- a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createUpdatedDynamicFvMesh.H
+++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createUpdatedDynamicFvMesh.H
@@ -7,7 +7,7 @@
         (
             IOobject
             (
-                dynamicFvMesh::defaultRegion,
+                polyMesh::defaultRegion,
                 runTime.timeName(),
                 runTime,
                 IOobject::MUST_READ
diff --git a/applications/test/mesh/Test-mesh.C b/applications/test/mesh/Test-mesh.C
index 386e5d8ede..49c64c4b51 100644
--- a/applications/test/mesh/Test-mesh.C
+++ b/applications/test/mesh/Test-mesh.C
@@ -47,7 +47,7 @@ int main(int argc, char *argv[])
     (
         IOobject
         (
-            fvMesh::defaultRegion,
+            polyMesh::defaultRegion,
             runTime.timeName(),
             runTime,
             IOobject::MUST_READ
diff --git a/applications/utilities/mesh/advanced/refinementLevel/refinementLevel.C b/applications/utilities/mesh/advanced/refinementLevel/refinementLevel.C
index a2332508f4..8d6a3c3870 100644
--- a/applications/utilities/mesh/advanced/refinementLevel/refinementLevel.C
+++ b/applications/utilities/mesh/advanced/refinementLevel/refinementLevel.C
@@ -212,7 +212,7 @@ int main(int argc, char *argv[])
     (
         IOobject
         (
-            fvMesh::defaultRegion,
+            polyMesh::defaultRegion,
             runTime.timeName(),
             runTime,
             IOobject::NO_READ,
diff --git a/applications/utilities/mesh/conversion/foamMeshToFluent/fluentFvMesh.H b/applications/utilities/mesh/conversion/foamMeshToFluent/fluentFvMesh.H
index f3cc64d45d..cfaf760e2b 100644
--- a/applications/utilities/mesh/conversion/foamMeshToFluent/fluentFvMesh.H
+++ b/applications/utilities/mesh/conversion/foamMeshToFluent/fluentFvMesh.H
@@ -56,7 +56,7 @@ public:
     // Constructors
 
         //- Construct from IOobject
-        fluentFvMesh(const IOobject& io);
+        explicit fluentFvMesh(const IOobject& io);
 
 
     // Member functions
diff --git a/applications/utilities/mesh/conversion/foamMeshToFluent/foamMeshToFluent.C b/applications/utilities/mesh/conversion/foamMeshToFluent/foamMeshToFluent.C
index de485f99e7..72a854d20f 100644
--- a/applications/utilities/mesh/conversion/foamMeshToFluent/foamMeshToFluent.C
+++ b/applications/utilities/mesh/conversion/foamMeshToFluent/foamMeshToFluent.C
@@ -59,7 +59,7 @@ int main(int argc, char *argv[])
     (
         IOobject
         (
-            fluentFvMesh::defaultRegion,
+            polyMesh::defaultRegion,
             runTime.constant(),
             runTime
         )
diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
index 8f9c0b9d7b..bdb9ce9270 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
@@ -233,7 +233,7 @@ int main(int argc, char *argv[])
     }
     else
     {
-        regionName = fvMesh::defaultRegion;
+        regionName = polyMesh::defaultRegion;
         Info<< "Create mesh for time = "
             << runTimeExtruded.timeName() << nl << endl;
     }
@@ -787,7 +787,7 @@ int main(int argc, char *argv[])
             (
                 IOobject
                 (
-                    extrudedMesh::defaultRegion,
+                    polyMesh::defaultRegion,
                     runTimeExtruded.constant(),
                     runTimeExtruded
                 ),
diff --git a/applications/utilities/mesh/manipulation/mirrorMesh/mirrorMesh.C b/applications/utilities/mesh/manipulation/mirrorMesh/mirrorMesh.C
index 5806a316f9..978cb86662 100644
--- a/applications/utilities/mesh/manipulation/mirrorMesh/mirrorMesh.C
+++ b/applications/utilities/mesh/manipulation/mirrorMesh/mirrorMesh.C
@@ -73,7 +73,7 @@ int main(int argc, char *argv[])
     (
         IOobject
         (
-            mirrorFvMesh::defaultRegion,
+            polyMesh::defaultRegion,
             runTime.constant(),
             runTime
         ),
diff --git a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
index 5e6fd78416..af9a933237 100644
--- a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
+++ b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
@@ -2905,7 +2905,7 @@ int main(int argc, char *argv[])
             const word& regionName = regionNames[regioni];
             const fileName meshSubDir
             (
-                regionName == fvMesh::defaultRegion
+                regionName == polyMesh::defaultRegion
               ? fileName(polyMesh::meshSubDir)
               : regionNames[regioni]/polyMesh::meshSubDir
             );
diff --git a/applications/utilities/preProcessing/changeDictionary/changeDictionary.C b/applications/utilities/preProcessing/changeDictionary/changeDictionary.C
index 8951d73ba7..458dc61e8a 100644
--- a/applications/utilities/preProcessing/changeDictionary/changeDictionary.C
+++ b/applications/utilities/preProcessing/changeDictionary/changeDictionary.C
@@ -499,7 +499,7 @@ int main(int argc, char *argv[])
 
 
         fileName regionPrefix;
-        if (regionName != fvMesh::defaultRegion)
+        if (regionName != polyMesh::defaultRegion)
         {
             regionPrefix = regionName;
         }
diff --git a/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C b/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C
index eb340016bc..3ebce53977 100644
--- a/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C
+++ b/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C
@@ -89,7 +89,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
 
-    wordList regionNames(1, fvMesh::defaultRegion);
+    wordList regionNames(1, polyMesh::defaultRegion);
     if (!args.readIfPresent("region", regionNames.first()))
     {
         args.readIfPresent("regions", regionNames);
diff --git a/applications/utilities/preProcessing/mapFields/mapFields.C b/applications/utilities/preProcessing/mapFields/mapFields.C
index 78645b20b4..047e613209 100644
--- a/applications/utilities/preProcessing/mapFields/mapFields.C
+++ b/applications/utilities/preProcessing/mapFields/mapFields.C
@@ -291,7 +291,7 @@ int main(int argc, char *argv[])
     const fileName caseDirSource = casePath.name();
 
     Info<< "Source: " << rootDirSource << " " << caseDirSource << endl;
-    word sourceRegion = fvMesh::defaultRegion;
+    word sourceRegion = polyMesh::defaultRegion;
     if (args.found("sourceRegion"))
     {
         sourceRegion = args["sourceRegion"];
@@ -299,7 +299,7 @@ int main(int argc, char *argv[])
     }
 
     Info<< "Target: " << rootDirTarget << " " << caseDirTarget << endl;
-    word targetRegion = fvMesh::defaultRegion;
+    word targetRegion = polyMesh::defaultRegion;
     if (args.found("targetRegion"))
     {
         targetRegion = args["targetRegion"];
diff --git a/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C b/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C
index 025c400b26..4bb265894b 100644
--- a/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C
+++ b/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C
@@ -219,14 +219,14 @@ int main(int argc, char *argv[])
     const fileName caseDirSource = casePath.name();
 
     Info<< "Source: " << rootDirSource << " " << caseDirSource << endl;
-    word sourceRegion = fvMesh::defaultRegion;
+    word sourceRegion = polyMesh::defaultRegion;
     if (args.readIfPresent("sourceRegion", sourceRegion))
     {
         Info<< "Source region: " << sourceRegion << endl;
     }
 
     Info<< "Target: " << rootDirTarget << " " << caseDirTarget << endl;
-    word targetRegion = fvMesh::defaultRegion;
+    word targetRegion = polyMesh::defaultRegion;
     if (args.readIfPresent("targetRegion", targetRegion))
     {
         Info<< "Target region: " << targetRegion << endl;
diff --git a/src/OpenFOAM/include/createMeshNoClear.H b/src/OpenFOAM/include/createMeshNoClear.H
index 321120aae8..5ccde9ca78 100644
--- a/src/OpenFOAM/include/createMeshNoClear.H
+++ b/src/OpenFOAM/include/createMeshNoClear.H
@@ -33,7 +33,7 @@ Foam::fvMesh mesh
 (
     Foam::IOobject
     (
-        Foam::fvMesh::defaultRegion,
+        Foam::polyMesh::defaultRegion,
         runTime.timeName(),
         runTime,
         Foam::IOobject::MUST_READ
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
index cda2419cbd..940de59331 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017, 2020 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -325,7 +325,7 @@ public:
     // Constructors
 
         //- Read construct from IOobject
-        polyMesh(const IOobject& io, const bool doInit = true);
+        explicit polyMesh(const IOobject& io, const bool doInit = true);
 
         //- Construct from IOobject or as zero-sized mesh
         //  Boundary is added using addPatches() member function
diff --git a/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMeshNew.C b/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMeshNew.C
index 482ff63e4c..487534b1d2 100644
--- a/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMeshNew.C
+++ b/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMeshNew.C
@@ -154,7 +154,7 @@ Foam::autoPtr<Foam::dynamicFvMesh> Foam::dynamicFvMesh::New
             (
                 IOobject
                 (
-                    dynamicFvMesh::defaultRegion,
+                    polyMesh::defaultRegion,
                     runTime.timeName(),
                     runTime,
                     IOobject::MUST_READ
@@ -168,7 +168,7 @@ Foam::autoPtr<Foam::dynamicFvMesh> Foam::dynamicFvMesh::New
             (
                 IOobject
                 (
-                    dynamicFvMesh::defaultRegion,
+                    polyMesh::defaultRegion,
                     runTime.timeName(),
                     runTime,
                     IOobject::MUST_READ
diff --git a/src/dynamicFvMesh/include/createNamedDynamicFvMesh.H b/src/dynamicFvMesh/include/createNamedDynamicFvMesh.H
index 21b6c092d8..6cb1444055 100644
--- a/src/dynamicFvMesh/include/createNamedDynamicFvMesh.H
+++ b/src/dynamicFvMesh/include/createNamedDynamicFvMesh.H
@@ -1,4 +1,4 @@
-    Foam::word regionName = Foam::fvMesh::defaultRegion;
+    Foam::word regionName = Foam::polyMesh::defaultRegion;
 
     if (args.readIfPresent("region", regionName))
     {
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
index f8b89e90b3..9745cb61db 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
@@ -1756,7 +1756,7 @@ Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
     (
         IOobject
         (
-            fvMesh::defaultRegion,
+            polyMesh::defaultRegion,
             runTime.timeName(),
             runTime,
             IOobject::NO_READ
diff --git a/src/engine/include/createEngineMesh.H b/src/engine/include/createEngineMesh.H
index 263ebed5ff..7a9daed5f4 100644
--- a/src/engine/include/createEngineMesh.H
+++ b/src/engine/include/createEngineMesh.H
@@ -7,7 +7,7 @@
         (
             IOobject
             (
-                engineMesh::defaultRegion,
+                polyMesh::defaultRegion,
                 runTime.timeName(),
                 runTime,
                 Foam::IOobject::MUST_READ
diff --git a/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.C b/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.C
index 39a0da94ed..c373ab314b 100644
--- a/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.C
+++ b/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.C
@@ -103,7 +103,7 @@ bool Foam::simplifiedMeshes::columnFvMeshInfo::setPatchEntries
         (
             runTime,
             runTime.timeName(),
-            (regionName_ == fvMesh::defaultRegion ? "" : regionName_)
+            (regionName_ == polyMesh::defaultRegion ? "" : regionName_)
         );
 
         if (objects.empty())
@@ -411,7 +411,7 @@ Foam::simplifiedMeshes::columnFvMeshInfo::columnFvMeshInfo
     regionName_(regionName),
     regionPrefix_
     (
-        regionName_ == fvMesh::defaultRegion
+        regionName_ == polyMesh::defaultRegion
       ? ""
       : regionName_ + '/'
     ),
diff --git a/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.H b/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.H
index 3bfd009bf7..4e345ac3ac 100644
--- a/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.H
+++ b/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.H
@@ -149,7 +149,7 @@ public:
     columnFvMesh
     (
         const Time& runTime,
-        const word& regionName = fvMesh::defaultRegion
+        const word& regionName = polyMesh::defaultRegion
     );
 };
 
diff --git a/src/functionObjects/field/streamLine/streamLineBase.C b/src/functionObjects/field/streamLine/streamLineBase.C
index 93b9ff90ac..184173bdca 100644
--- a/src/functionObjects/field/streamLine/streamLineBase.C
+++ b/src/functionObjects/field/streamLine/streamLineBase.C
@@ -663,7 +663,7 @@ bool Foam::functionObjects::streamLineBase::writeToFile()
         (
             time_.globalPath()/functionObject::outputPrefix/"sets"/name()
         );
-        if (mesh_.name() != fvMesh::defaultRegion)
+        if (mesh_.name() != polyMesh::defaultRegion)
         {
             vtkPath = vtkPath/mesh_.name();
         }
diff --git a/src/sampling/sampledSet/sampledSets/sampledSets.C b/src/sampling/sampledSet/sampledSets/sampledSets.C
index 965cad242f..2420bc7ac5 100644
--- a/src/sampling/sampledSet/sampledSets/sampledSets.C
+++ b/src/sampling/sampledSet/sampledSets/sampledSets.C
@@ -104,7 +104,7 @@ Foam::sampledSets::sampledSets
         mesh_.time().globalPath()/functionObject::outputPrefix/name
     );
 
-    if (mesh_.name() != fvMesh::defaultRegion)
+    if (mesh_.name() != polyMesh::defaultRegion)
     {
         outputPath_ = outputPath_/mesh_.name();
     }
@@ -137,7 +137,7 @@ Foam::sampledSets::sampledSets
         mesh_.time().globalPath()/functionObject::outputPrefix/name
     );
 
-    if (mesh_.name() != fvMesh::defaultRegion)
+    if (mesh_.name() != polyMesh::defaultRegion)
     {
         outputPath_ = outputPath_/mesh_.name();
     }

From 5022591c368a4b3746983fde71604b03a8d3f41d Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 11 May 2021 13:17:43 +0200
Subject: [PATCH 8/8] ENH: return old file handler on change

- intended for the following type of use:

      auto oldHandler = fileHandler(fileOperation::NewUncollated());

      ... do something that only works with uncollated

      // Restore previous (if any)
      if (oldHandler)
      {
          fileHandler(std::move(oldHandler));
      }

ENH: make fileOperation distributed(bool) mutable

- use is "static-like" and akin to Pstream::parRun(bool),
  thus allow toggling of the switch without a const_cast
---
 .../redistributePar/redistributePar.C         |  2 +-
 .../utilities/postProcessing/noise/noise.C    |  4 +--
 .../fileOperation/fileOperation.C             | 28 +++++++++++++++++--
 .../fileOperation/fileOperation.H             | 25 +++++++++--------
 4 files changed, 42 insertions(+), 17 deletions(-)

diff --git a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
index af9a933237..7180ede66b 100644
--- a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
+++ b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
@@ -2303,7 +2303,7 @@ int main(int argc, char *argv[])
     //                than it writes to
     // - reconstruct - reads parallel, write on master only and to parent
     //                 directory
-    const_cast<fileOperation&>(fileHandler()).distributed(true);
+    fileHandler().distributed(true);
 
 
     #include "foamDlOpenLibs.H"
diff --git a/applications/utilities/postProcessing/noise/noise.C b/applications/utilities/postProcessing/noise/noise.C
index 505e67b77e..19f2eb74fc 100644
--- a/applications/utilities/postProcessing/noise/noise.C
+++ b/applications/utilities/postProcessing/noise/noise.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -116,7 +116,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
 
     // As much as possible avoid synchronised operation
-    const_cast<fileOperation&>(fileHandler()).distributed(true);
+    fileHandler().distributed(true);
 
     #include "createTime.H"
 
diff --git a/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.C b/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.C
index 60c396b559..3919c92f96 100644
--- a/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.C
+++ b/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.C
@@ -729,6 +729,14 @@ Foam::fileOperation::New
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::fileOperation::distributed(bool on) const noexcept
+{
+    bool old(distributed_);
+    distributed_ = on;
+    return old;
+}
+
+
 Foam::fileName Foam::fileOperation::objectPath
 (
     const IOobject& io,
@@ -1461,6 +1469,17 @@ Foam::label Foam::fileOperation::detectProcessorPath(const fileName& fName)
 }
 
 
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::fileOperation> Foam::fileOperation::NewUncollated()
+{
+    return autoPtr<fileOperation>
+    (
+        new fileOperations::uncollatedFileOperation(false)
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
 const Foam::fileOperation& Foam::fileHandler()
@@ -1481,7 +1500,8 @@ const Foam::fileOperation& Foam::fileHandler()
 }
 
 
-void Foam::fileHandler(autoPtr<fileOperation>&& newHandler)
+Foam::autoPtr<Foam::fileOperation>
+Foam::fileHandler(autoPtr<fileOperation>&& newHandler)
 {
     if
     (
@@ -1490,10 +1510,14 @@ void Foam::fileHandler(autoPtr<fileOperation>&& newHandler)
      && newHandler->type() == fileOperation::fileHandlerPtr_->type()
     )
     {
-        return;
+        return nullptr;  // No change
     }
 
+    autoPtr<fileOperation> old(std::move(fileOperation::fileHandlerPtr_));
+
     fileOperation::fileHandlerPtr_ = std::move(newHandler);
+
+    return old;
 }
 
 
diff --git a/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.H b/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.H
index 0b659cdc71..2cc6ad7f00 100644
--- a/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.H
+++ b/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.H
@@ -112,7 +112,7 @@ protected:
         const label comm_;
 
         //- Distributed roots (parallel run)
-        bool distributed_;
+        mutable bool distributed_;
 
         //- Detected processors directories
         mutable HashTable<dirIndexList> procsDirs_;
@@ -181,6 +181,9 @@ public:
     //- Static fileOperation
     static autoPtr<fileOperation> fileHandlerPtr_;
 
+    //- Static construct the commonly used uncollatedFileOperation
+    static autoPtr<fileOperation> NewUncollated();
+
 
     // Constructors
 
@@ -239,13 +242,9 @@ public:
             return distributed_;
         }
 
-        //- Set distributed roots on/off, return old value
-        bool distributed(bool on) noexcept
-        {
-            bool old(distributed_);
-            distributed_ = on;
-            return old;
-        }
+        //- Set distributed roots on/off (mutable)
+        //  \return old value
+        bool distributed(bool on) const noexcept;
 
 
         // OSSpecific equivalents
@@ -595,14 +594,16 @@ public:
 };
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Global declarations: defined in fileOperation.C
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+// Note: defined in fileOperation.C
 
 //- Get current file handler
 const fileOperation& fileHandler();
 
-//- Replace, reset file handler
-void fileHandler(autoPtr<fileOperation>&& newHandler);
+//- Replace, reset file handler.
+//  \return old handler on change, null otherwise
+autoPtr<fileOperation> fileHandler(autoPtr<fileOperation>&& newHandler);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //