diff --git a/applications/utilities/mesh/generation/blockMesh/blockMesh.C b/applications/utilities/mesh/generation/blockMesh/blockMesh.C
index 8cdb4faf0ed0de114b6dc4109f41b1957ee50245..520f4b2105ee13f03244455030cd73e6448d5da9 100644
--- a/applications/utilities/mesh/generation/blockMesh/blockMesh.C
+++ b/applications/utilities/mesh/generation/blockMesh/blockMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -96,23 +96,24 @@ int main(int argc, char *argv[])
     argList::addNote
     (
         "Block mesh generator.\n"
+        " \n"
         "  The ordering of vertex and face labels within a block as shown "
         "below.\n"
         "  For the local vertex numbering in the sequence 0 to 7:\n"
         "    Faces 0, 1 (x-direction) are left, right.\n"
         "    Faces 2, 3 (y-direction) are front, back.\n"
         "    Faces 4, 5 (z-direction) are bottom, top.\n"
-        "\n"
+        " \n"
         "                        7 ---- 6\n"
-        "                 f5     |\\     |\\     f3\n"
+        "                 f5     |\\     :\\     f3\n"
         "                 |      | 4 ---- 5     \\\n"
-        "                 |      3 |--- 2 |      \\\n"
+        "                 |      3.|....2 |      \\\n"
         "                 |       \\|     \\|      f2\n"
         "                 f4       0 ---- 1\n"
         "    Y  Z\n"
         "     \\ |                f0 ------ f1\n"
         "      \\|\n"
-        "       O--- X\n"
+        "       o--- X\n"
     );
 
     argList::noParallel();
@@ -129,7 +130,11 @@ int main(int argc, char *argv[])
     argList::addBoolOption
     (
         "write-vtk",
-        "Write topology as VTK file and exit"
+        "Write topology as VTU file and exit"
+    );
+    argList::addVerboseOption
+    (
+        "Force verbose output. (Can be used multiple times)"
     );
 
     argList::addBoolOption
@@ -170,12 +175,12 @@ int main(int argc, char *argv[])
     const bool writeCellSets = args.found("sets");
 
     // Default merge (topology), unless otherwise specified
-    blockMesh::mergeStrategy strategy(blockMesh::DEFAULT_MERGE);
-
-    if (args.found("merge-points"))
-    {
-        strategy = blockMesh::MERGE_POINTS;
-    }
+    const blockMesh::mergeStrategy strategy =
+    (
+        args.found("merge-points")
+      ? blockMesh::MERGE_POINTS
+      : blockMesh::DEFAULT_MERGE
+    );
 
     word regionName(polyMesh::defaultRegion);
     word regionPath;
@@ -216,7 +221,7 @@ int main(int argc, char *argv[])
     // Locate appropriate blockMeshDict
     #include "findBlockMeshDict.H"
 
-    blockMesh blocks(meshDict, regionName, strategy);
+    blockMesh blocks(meshDict, regionName, strategy, args.verbose());
 
     if (!blocks.valid())
     {
diff --git a/src/OpenFOAM/meshes/meshShapes/hexCell/hexCell.C b/src/OpenFOAM/meshes/meshShapes/hexCell/hexCell.C
index 15fd47bf3f79e85af6638929381bf4e2c02b3fbd..98cc0654c14bdc65b7f42ee1eeee688c9337b1f3 100644
--- a/src/OpenFOAM/meshes/meshShapes/hexCell/hexCell.C
+++ b/src/OpenFOAM/meshes/meshShapes/hexCell/hexCell.C
@@ -74,17 +74,16 @@ const Foam::faceList& Foam::hexCell::modelFaces()
 
     if (!ptr)
     {
-        ptr.reset(new Foam::faceList(6));
+        ptr.reset(new Foam::faceList(hexCell::nFaces(), Foam::face(4)));
 
-        for (label facei = 0; facei < 6; ++facei)
+        label facei = 0;
+        for (auto& f : *ptr)
         {
-            auto& f = (*ptr)[facei];
-
-            f.resize(4);
             f[0] = modelFaces_[facei][0];
             f[1] = modelFaces_[facei][1];
             f[2] = modelFaces_[facei][2];
             f[3] = modelFaces_[facei][3];
+            ++facei;
         }
     }
 
@@ -98,14 +97,14 @@ const Foam::edgeList& Foam::hexCell::modelEdges()
 
     if (!ptr)
     {
-        ptr.reset(new Foam::edgeList(12));
+        ptr.reset(new Foam::edgeList(hexCell::nEdges()));
 
-        for (label edgei = 0; edgei < 12; ++edgei)
+        label edgei = 0;
+        for (auto& e : *ptr)
         {
-            auto& e = (*ptr)[edgei];
-
-            e.first()  = modelEdges_[edgei][0];
-            e.second() = modelEdges_[edgei][1];
+            e[0] = modelEdges_[edgei][0];
+            e[1] = modelEdges_[edgei][1];
+            ++edgei;
         }
     }
 
@@ -115,39 +114,35 @@ const Foam::edgeList& Foam::hexCell::modelEdges()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-/// Foam::faceList Foam::hexCell::faces() const
-/// {
-///     Foam::faceList result(6);
-///
-///     for (label facei = 0; facei < 6; ++facei)
-///     {
-///         auto& f = result[facei];
-///
-///         f.resize(4);
-///         f[0] = (*this)[modelFaces_[facei][0]];
-///         f[1] = (*this)[modelFaces_[facei][1]];
-///         f[2] = (*this)[modelFaces_[facei][2]];
-///         f[3] = (*this)[modelFaces_[facei][3]];
-///     }
-///
-///     return result;
-/// }
-///
-///
-/// Foam::edgeList Foam::hexCell::edges() const
-/// {
-///     Foam::edgeList result(12);
-///
-///     for (label edgei = 0; edgei < 12; ++edgei)
-///     {
-///         auto& e = result[edgei];
-///
-///         e.first()  = (*this)[modelEdges_[edgei][0]],
-///         e.second() = (*this)[modelEdges_[edgei][1]]
-///     }
-///
-///     return result;
-/// }
+Foam::faceList Foam::hexCell::faces() const
+{
+    Foam::faceList theFaces(hexCell::nFaces(), Foam::face(4));
+
+    label facei = 0;
+    for (auto& f : theFaces)
+    {
+        copyFace(f, facei);
+        ++facei;
+    }
+
+    return theFaces;
+}
+
+
+Foam::edgeList Foam::hexCell::edges() const
+{
+    Foam::edgeList theEdges(hexCell::nEdges());
+
+    label edgei = 0;
+    for (auto& e : theEdges)
+    {
+        e[0] = (*this)[modelEdges_[edgei][0]];
+        e[1] = (*this)[modelEdges_[edgei][1]];
+        ++edgei;
+    }
+
+    return theEdges;
+}
 
 
 Foam::cellShape Foam::hexCell::shape(const bool doCollapse) const
diff --git a/src/OpenFOAM/meshes/meshShapes/hexCell/hexCell.H b/src/OpenFOAM/meshes/meshShapes/hexCell/hexCell.H
index fb6b9ddace20f8fcd6af17948cb4d1100c68a157..07817fa8ea3207737bc383a53adfa4c6b5fc00a9 100644
--- a/src/OpenFOAM/meshes/meshShapes/hexCell/hexCell.H
+++ b/src/OpenFOAM/meshes/meshShapes/hexCell/hexCell.H
@@ -71,6 +71,12 @@ class hexCell
         static const label modelEdges_[12][2];
 
 
+    // Private Member Functions
+
+        //- Copy vertices for given face - no checks
+        inline void copyFace(Foam::face& f, const label facei) const;
+
+
 public:
 
     // Constructors
@@ -139,6 +145,15 @@ public:
         //- Return i-th edge reversed
         inline Foam::edge reverseEdge(const label edgei) const;
 
+        //- Return list of cell faces [6]
+        Foam::faceList faces() const;
+
+        //- Return list of cell edges [12]
+        Foam::edgeList edges() const;
+
+        //- Cell centre - uses simple average of points
+        inline point centre(const UList<point>& meshPoints) const;
+
         //- The points corresponding to this shape
         inline pointField points(const UList<point>& meshPoints) const;
 
diff --git a/src/OpenFOAM/meshes/meshShapes/hexCell/hexCellI.H b/src/OpenFOAM/meshes/meshShapes/hexCell/hexCellI.H
index 4ec3870b1c7ceeba9413351496ae4304e1ab8127..bd2d319cbf99cd31f5816846c5d81f8e2322ff87 100644
--- a/src/OpenFOAM/meshes/meshShapes/hexCell/hexCellI.H
+++ b/src/OpenFOAM/meshes/meshShapes/hexCell/hexCellI.H
@@ -25,6 +25,17 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline void Foam::hexCell::copyFace(Foam::face& f, const label facei) const
+{
+    f[0] = (*this)[modelFaces_[facei][0]];
+    f[1] = (*this)[modelFaces_[facei][1]];
+    f[2] = (*this)[modelFaces_[facei][2]];
+    f[3] = (*this)[modelFaces_[facei][3]];
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 inline Foam::hexCell::hexCell()
@@ -77,7 +88,7 @@ inline Foam::hexCell::hexCell(Istream& is)
 inline Foam::face Foam::hexCell::face(const label facei) const
 {
     #ifdef FULLDEBUG
-    if (facei < 0 || facei >= 6)
+    if (facei < 0 || facei >= hexCell::nFaces())
     {
         FatalErrorInFunction
             << "Face index (" << facei << ") out of range 0..5\n"
@@ -86,10 +97,7 @@ inline Foam::face Foam::hexCell::face(const label facei) const
     #endif
 
     Foam::face f(4);
-    f[0] = (*this)[modelFaces_[facei][0]];
-    f[1] = (*this)[modelFaces_[facei][1]];
-    f[2] = (*this)[modelFaces_[facei][2]];
-    f[3] = (*this)[modelFaces_[facei][3]];
+    copyFace(f, facei);
 
     return f;
 }
@@ -98,7 +106,7 @@ inline Foam::face Foam::hexCell::face(const label facei) const
 inline Foam::edge Foam::hexCell::edge(const label edgei) const
 {
     #ifdef FULLDEBUG
-    if (edgei < 0 || edgei >= 12)
+    if (edgei < 0 || edgei >= hexCell::nEdges())
     {
         FatalErrorInFunction
             << "Edge index (" << edgei << ") out of range 0..11\n"
@@ -121,6 +129,27 @@ inline Foam::edge Foam::hexCell::reverseEdge(const label edgei) const
 }
 
 
+inline Foam::point Foam::hexCell::centre
+(
+    const UList<point>& meshPoints
+) const
+{
+    // Simple estimate of cell centre by averaging cell points
+    point cEst = Zero;
+    int npts = 0;
+    for (const label pointi : *this)
+    {
+        if (pointi >= 0)
+        {
+            cEst += meshPoints[pointi];
+            ++npts;
+        }
+    }
+
+    return (npts > 1 ? (cEst/scalar(npts)) : cEst);
+}
+
+
 inline Foam::pointField Foam::hexCell::points
 (
     const UList<point>& meshPoints
diff --git a/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.C b/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.C
index cdeb5f25d89d7bf89e197ca381b6dbb64896a7e9..d05982a22e8892f2da765dc2d6873f8d431351fd 100644
--- a/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.C
+++ b/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.C
@@ -67,16 +67,15 @@ const Foam::faceList& Foam::tetCell::modelFaces()
 
     if (!ptr)
     {
-        ptr.reset(new Foam::faceList(4));
+        ptr.reset(new Foam::faceList(tetCell::nFaces(), Foam::face(3)));
 
-        for (label facei = 0; facei < 4; ++facei)
+        label facei = 0;
+        for (auto& f : *ptr)
         {
-            auto& f = (*ptr)[facei];
-
-            f.resize(3);
             f[0] = modelFaces_[facei][0];
             f[1] = modelFaces_[facei][1];
             f[2] = modelFaces_[facei][2];
+            ++facei;
         }
     }
 
@@ -90,14 +89,14 @@ const Foam::edgeList& Foam::tetCell::modelEdges()
 
     if (!ptr)
     {
-        ptr.reset(new Foam::edgeList(6));
+        ptr.reset(new Foam::edgeList(tetCell::nEdges()));
 
-        for (label edgei = 0; edgei < 6; ++edgei)
+        label edgei = 0;
+        for (auto& e : *ptr)
         {
-            auto& e = (*ptr)[edgei];
-
-            e.first()  = modelEdges_[edgei][0];
-            e.second() = modelEdges_[edgei][1];
+            e[0] = modelEdges_[edgei][0];
+            e[1] = modelEdges_[edgei][1];
+            ++edgei;
         }
     }
 
@@ -107,35 +106,34 @@ const Foam::edgeList& Foam::tetCell::modelEdges()
 
 /// Foam::faceList Foam::tetCell::faces() const
 /// {
-///     Foam::faceList result(4);
+///     Foam::faceList theFaces(tetCell::nFaces(), Foam::face(3));
 ///
-///     for (label facei = 0; facei < 4; ++facei)
+///     label facei = 0;
+///     for (auto& f : theFaces)
 ///     {
-///         auto& f = result[facei];
-///
-///         f.resize(3);
 ///         f[0] = (*this)[modelFaces_[facei][0]];
 ///         f[1] = (*this)[modelFaces_[facei][1]];
 ///         f[2] = (*this)[modelFaces_[facei][2]];
+///         ++facei;
 ///     }
 ///
-///     return result;
+///     return theFaces;
 /// }
 ///
 ///
 /// Foam::edgeList Foam::tetCell::edges() const
 /// {
-///     Foam::edgeList result(6);
+///     Foam::edgeList theEdges(tetCell::nEdges());
 ///
-///     for (label edgei = 0; edgei < 6; ++edgei)
+///     label edgei = 0;
+///     for (auto& e : theEdges)
 ///     {
-///         auto& e = result[edgei];
-///
-///         e.first()  = (*this)[modelEdges_[edgei][0]],
-///         e.second() = (*this)[modelEdges_[edgei][1]]
+///         e[0] = (*this)[modelEdges_[edgei][0]];
+///         e[1] = (*this)[modelEdges_[edgei][1]];
+///         ++edgei;
 ///     }
 ///
-///     return result;
+///     return theEdges;
 /// }
 
 
diff --git a/src/OpenFOAM/meshes/meshShapes/tetCell/tetCellI.H b/src/OpenFOAM/meshes/meshShapes/tetCell/tetCellI.H
index 21761e7e9a35ba96881ef43dd3f43e031b6743fc..66f90767ff9b8de563cca114dbe4dd2668e984c2 100644
--- a/src/OpenFOAM/meshes/meshShapes/tetCell/tetCellI.H
+++ b/src/OpenFOAM/meshes/meshShapes/tetCell/tetCellI.H
@@ -93,7 +93,7 @@ inline Foam::tetCell::tetCell(Istream& is)
 inline Foam::triFace Foam::tetCell::face(const label facei) const
 {
     #ifdef FULLDEBUG
-    if (facei < 0 || facei >= 4)
+    if (facei < 0 || facei >= tetCell::nFaces())
     {
         FatalErrorInFunction
             << "Face index (" << facei << ") out of range 0..3\n"
@@ -117,7 +117,7 @@ inline Foam::label Foam::tetCell::edgeFace(const label edgei) const
     static const label edgeFaces[6] = {2, 3, 1, 0, 0, 1};
 
     #ifdef FULLDEBUG
-    if (edgei < 0 || edgei >= 6)
+    if (edgei < 0 || edgei >= tetCell::nEdges())
     {
         FatalErrorInFunction
             << "Edge index (" << edgei << ") out of range 0..5\n"
@@ -148,14 +148,14 @@ inline Foam::label Foam::tetCell::edgeAdjacentFace
     };
 
     #ifdef FULLDEBUG
-    if (facei < 0 || facei >= 4)
+    if (facei < 0 || facei >= tetCell::nFaces())
     {
         FatalErrorInFunction
             << "Face index (" << facei << ") out of range 0..3\n"
             << abort(FatalError);
     }
 
-    if (edgei < 0 || edgei >= 6)
+    if (edgei < 0 || edgei >= tetCell::nEdges())
     {
         FatalErrorInFunction
             << "Edge index (" << edgei << ") out of range 0..5\n"
@@ -170,7 +170,7 @@ inline Foam::label Foam::tetCell::edgeAdjacentFace
 inline Foam::edge Foam::tetCell::edge(const label edgei) const
 {
     #ifdef FULLDEBUG
-    if (edgei < 0 || edgei >= 6)
+    if (edgei < 0 || edgei >= tetCell::nEdges())
     {
         FatalErrorInFunction
             << "Edge index (" << edgei << ") out of range 0..5\n"
diff --git a/src/OpenFOAM/meshes/primitiveShapes/line/line.H b/src/OpenFOAM/meshes/primitiveShapes/line/line.H
index caf65ef76cdc551efbc2bc038a3c257890595c41..f2771a4bc0778779361a2e07f22fddc8a8f9a083 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/line/line.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/line/line.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -67,7 +67,7 @@ inline Ostream& operator<<(Ostream& os, const line<Point, PointRef>& l);
 template<class Point, class PointRef>
 class line
 {
-    // Private data
+    // Private Data
 
         //- First point
         PointRef a_;
@@ -81,7 +81,7 @@ public:
     // Constructors
 
         //- Construct from two points
-        inline line(const Point& start, const Point& end);
+        inline line(const Point& from, const Point& to);
 
         //- Construct from two points in the list of points
         //  The indices could be from edge etc.
@@ -100,19 +100,19 @@ public:
     // Access
 
         //- Return first point
-        inline PointRef first() const;
+        inline PointRef first() const noexcept;
 
         //- Return second (last) point
-        inline PointRef second() const;
+        inline PointRef second() const noexcept;
 
         //- Return last (second) point
-        inline PointRef last() const;
+        inline PointRef last() const noexcept;
 
         //- Return first point
-        inline PointRef start() const;
+        inline PointRef start() const noexcept;
 
         //- Return second (last) point
-        inline PointRef end() const;
+        inline PointRef end() const noexcept;
 
 
     // Properties
diff --git a/src/OpenFOAM/meshes/primitiveShapes/line/lineI.H b/src/OpenFOAM/meshes/primitiveShapes/line/lineI.H
index ec03aa43b1fcedfcb7370217b2563ffead932182..3ed1e17d457d94c225ff49a77c6d1ceffbf370b5 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/line/lineI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/line/lineI.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,10 +32,10 @@ License
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Point, class PointRef>
-inline Foam::line<Point, PointRef>::line(const Point& start, const Point& end)
+inline Foam::line<Point, PointRef>::line(const Point& from, const Point& to)
 :
-    a_(start),
-    b_(end)
+    a_(from),
+    b_(to)
 {}
 
 
@@ -61,34 +61,34 @@ inline Foam::line<Point, PointRef>::line(Istream& is)
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Point, class PointRef>
-inline PointRef Foam::line<Point, PointRef>::first() const
+inline PointRef Foam::line<Point, PointRef>::first() const noexcept
 {
     return a_;
 }
 
 
 template<class Point, class PointRef>
-inline PointRef Foam::line<Point, PointRef>::second() const
+inline PointRef Foam::line<Point, PointRef>::second() const noexcept
 {
     return b_;
 }
 
 
 template<class Point, class PointRef>
-inline PointRef Foam::line<Point, PointRef>::last() const
+inline PointRef Foam::line<Point, PointRef>::last() const noexcept
 {
     return b_;
 }
 
 
 template<class Point, class PointRef>
-inline PointRef Foam::line<Point, PointRef>::start() const
+inline PointRef Foam::line<Point, PointRef>::start() const noexcept
 {
     return first();
 }
 
 template<class Point, class PointRef>
-inline PointRef Foam::line<Point, PointRef>::end() const
+inline PointRef Foam::line<Point, PointRef>::end() const noexcept
 {
     return second();
 }
@@ -171,7 +171,7 @@ Foam::scalar Foam::line<Point, PointRef>::nearestDist
     Point c(edge.start() - start());
 
     Point crossab = a ^ b;
-    scalar magCrossSqr = magSqr(crossab);
+    const scalar magCrossSqr = Foam::magSqr(crossab);
 
     if (magCrossSqr > VSMALL)
     {
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
index acc1bf549be9a5b0d4e62b05033d33dd30e8a470..a5a66a23deb9613c1f640663cd2ad1f2ad7b26cc 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
@@ -210,6 +210,7 @@ Foam::blockDescriptor::blockDescriptor
     blockFaces_(faces),
     blockShape_(bshape),
     expand_(),
+    index_(-1),
     zoneName_(zoneName),
     curvedFaces_(-1),
     nCurvedFaces_(0)
@@ -241,6 +242,7 @@ Foam::blockDescriptor::blockDescriptor
     blockFaces_(faces),
     blockShape_(),
     expand_(),
+    index_(blockIndex),
     zoneName_(),
     curvedFaces_(-1),
     nCurvedFaces_(0)
@@ -383,7 +385,7 @@ void Foam::blockDescriptor::correctFacePoints
 {
     forAll(curvedFaces_, blockFacei)
     {
-        if (curvedFaces_[blockFacei] != -1)
+        if (curvedFaces_[blockFacei] >= 0)
         {
             blockFaces_[curvedFaces_[blockFacei]].project
             (
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
index 4af6734cd3c6ead18d4977409b4fca54b7296b94..9c04e352a6589aad5717df923a434e4e61a793fe 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
@@ -35,18 +35,18 @@ Description
     and face labels is shown below.  For vertex numbering in the sequence 0 to 7
     (block, centre): faces 0 (f0) and 1 are left and right, respectively; faces
     2 and 3 are front and back; and faces 4 and 5 are bottom and top:
+
     \verbatim
-                   7 ---- 6
-              f5   |\     |\   f3
-               |   | 4 ---- 5    \
-               |   3 |--- 2 |     \
-               |    \|     \|      f2
-              f4     0 ---- 1
-
-         Z         f0 ----- f1
-         |  Y
-         | /
-         O --- X
+                    7 ---- 6
+              f5    |\     :\    f3
+               |    | 4 ---- 5    \
+               |    3.|....2 |     \
+               |     \|     \|     f2
+              f4      0 ---- 1
+      Y  Z
+       \ |          f0 ------ f1
+        \|
+         o--- X
      \endverbatim
 
 SourceFiles
@@ -91,12 +91,15 @@ class blockDescriptor
         //- Reference to the list of curved faces
         const blockFaceList& blockFaces_;
 
-        //- Block shape
+        //- Block shape. Likely only hex is supportable
         cellShape blockShape_;
 
         //- Expansion ratios in all directions
         List<gradingDescriptors> expand_;
 
+        //- The block index in the originating list (-1 for unknown)
+        label index_;
+
         //- Name of the zone (empty word if none)
         word zoneName_;
 
@@ -122,8 +125,7 @@ class blockDescriptor
         (
             pointField& edgePoints,
             scalarList& edgeWeights,
-            const label start,
-            const label end,
+            const Foam::edge& cellModelEdge,
             const label nDiv,
             const gradingDescriptors& expand
         ) const;
@@ -179,7 +181,7 @@ public:
         //- Return the block shape
         inline const cellShape& blockShape() const noexcept;
 
-        //- Return the mesh density (number of cells) in the i,j,k directions
+        //- The mesh density (number of cells) in the i,j,k directions
         inline const labelVector& density() const noexcept;
 
         //- Expansion ratios in all directions
@@ -206,10 +208,10 @@ public:
             const label j
         ) const;
 
-        //- Return true if point i,j,k addresses a block vertex
+        //- True if point i,j,k addresses a block vertex
         inline bool vertex(const label i, const label j, const label k) const;
 
-        //- Return true if point i,j,k addresses a block edge
+        //- True if point i,j,k addresses a block edge
         inline bool edge(const label i, const label j, const label k) const;
 
         //- Calculate the points and weights for all edges.
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
index af615b13bbe5be8e7a7b94294d2a06a9797bb3f6..f296a4c8ef88a7f864602ba92ae806ceebedf569 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
@@ -29,14 +29,7 @@ License
 #include "blockDescriptor.H"
 #include "lineEdge.H"
 #include "lineDivide.H"
-
-// * * * * * * * * * * * * * * Local Data Members  * * * * * * * * * * * * * //
-
-// Warning.
-// Ordering of edges needs to be the same as hex cell shape model
-static const int hexEdge0[12] = { 0, 3, 7, 4,  0, 1, 5, 4,  0, 1, 2, 3 };
-static const int hexEdge1[12] = { 1, 2, 6, 5,  3, 2, 6, 7,  4, 5, 6, 7 };
-
+#include "hexCell.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -44,21 +37,35 @@ int Foam::blockDescriptor::calcEdgePointsWeights
 (
     pointField& edgePoints,
     scalarList& edgeWeights,
-    const label start,
-    const label end,
+    const Foam::edge& cellModelEdge,
     const label nDiv,
     const gradingDescriptors& expand
 ) const
 {
-    // Set reference to the list of labels defining the block
-    const labelList& blockLabels = blockShape_;
+    // The topological edge on the block
+    const Foam::edge thisEdge(blockShape_, cellModelEdge);
+
+    const bool isCollapsedEdge = !thisEdge.valid();
+
+    if (blockEdge::debug && isCollapsedEdge)
+    {
+        Info<< "Collapsed edge:" << thisEdge;
+        if (index_ >= 0)
+        {
+            Info << " block:" << index_;
+        }
+        Info<< " model edge:" << cellModelEdge << nl;
+    }
+
+    // FUTURE: skip point generation for collapsed edge
+
 
     // Set the edge points/weights
     // The edge is a straight-line if it is not in the list of blockEdges
 
     for (const blockEdge& cedge : blockEdges_)
     {
-        const int cmp = cedge.compare(blockLabels[start], blockLabels[end]);
+        const int cmp = cedge.compare(thisEdge);
 
         if (cmp > 0)
         {
@@ -104,7 +111,7 @@ int Foam::blockDescriptor::calcEdgePointsWeights
 
     lineDivide divEdge
     (
-        blockEdges::lineEdge(blockPoints, start, end),
+        blockEdges::lineEdge(blockPoints, cellModelEdge),
         nDiv,
         expand
     );
@@ -126,14 +133,13 @@ int Foam::blockDescriptor::edgesPointsWeights
 {
     int nCurved = 0;
 
-    for (label edgei = 0; edgei < 12; ++edgei)
+    for (label edgei = 0; edgei < 12; ++edgei)  //< hexCell::nEdges()
     {
         nCurved += calcEdgePointsWeights
         (
             edgesPoints[edgei],
             edgesWeights[edgei],
-            hexEdge0[edgei],
-            hexEdge1[edgei],
+            hexCell::modelEdges()[edgei],
 
             sizes()[edgei/4],   // 12 edges -> 3 components (x,y,z)
             expand_[edgei]
@@ -153,7 +159,7 @@ bool Foam::blockDescriptor::edgePointsWeights
     const gradingDescriptors& gd
 ) const
 {
-    if (edgei < 0 || edgei >= 12)
+    if (edgei < 0 || edgei >= 12)  //< hexCell::nEdges()
     {
         FatalErrorInFunction
             << "Edge label " << edgei
@@ -165,8 +171,8 @@ bool Foam::blockDescriptor::edgePointsWeights
     (
         edgePoints,
         edgeWeights,
-        hexEdge0[edgei],
-        hexEdge1[edgei],
+        hexCell::modelEdges()[edgei],
+
         nDiv,
         gd
     );
@@ -182,7 +188,7 @@ bool Foam::blockDescriptor::edgePointsWeights
     scalarList& edgeWeights
 ) const
 {
-    if (edgei < 0 || edgei >= 12)
+    if (edgei < 0 || edgei >= 12)  //< hexCell::nEdges()
     {
         FatalErrorInFunction
             << "Edge label " << edgei
@@ -194,8 +200,8 @@ bool Foam::blockDescriptor::edgePointsWeights
     (
         edgePoints,
         edgeWeights,
-        hexEdge0[edgei],
-        hexEdge1[edgei],
+        hexCell::modelEdges()[edgei],
+
         sizes()[edgei/4],   // 12 edges -> 3 components (x,y,z)
         expand_[edgei]
     );
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
index cc8c83cb2236738d34a88812bd75014ea96f9580..736fa5eff495a31f49e9e937595058db9be2ad3c 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
@@ -156,14 +156,14 @@ inline bool Foam::blockDescriptor::flatFaceOrEdge
     const label i, const label j, const label k
 ) const
 {
-    if (i == 0 && curvedFaces_[0] == -1) return true;
-    if (i == sizes().x() && curvedFaces_[1] == -1) return true;
+    if (i == 0 && curvedFaces_[0] < 0) return true;
+    if (i == sizes().x() && curvedFaces_[1] < 0) return true;
 
-    if (j == 0 && curvedFaces_[2] == -1) return true;
-    if (j == sizes().y() && curvedFaces_[3] == -1) return true;
+    if (j == 0 && curvedFaces_[2] < 0) return true;
+    if (j == sizes().y() && curvedFaces_[3] < 0) return true;
 
-    if (k == 0 && curvedFaces_[4] == -1) return true;
-    if (k == sizes().z() && curvedFaces_[5] == -1) return true;
+    if (k == 0 && curvedFaces_[4] < 0) return true;
+    if (k == sizes().z() && curvedFaces_[5] < 0) return true;
 
     return this->edge(i, j, k);
 }
diff --git a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C
index af5e4d89ed8f763dc6b31a3c70ace092f2889019..c9be346db8ec0f2b997a025ce906ba081f433c68 100644
--- a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C
+++ b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2014-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -53,24 +53,35 @@ namespace blockEdges
 Foam::blockEdges::BSplineEdge::BSplineEdge
 (
     const pointField& points,
-    const label start,
-    const label end,
+    const edge& fromTo,
     const pointField& internalPoints
 )
 :
-    blockEdge(points, start, end),
+    blockEdge(points, fromTo),
     BSpline
     (
-        polyLine::concat(points[start_], internalPoints, points[end_])
+        polyLine::concat(firstPoint(), internalPoints, lastPoint())
     )
 {}
 
 
+Foam::blockEdges::BSplineEdge::BSplineEdge
+(
+    const pointField& points,
+    const label from,
+    const label to,
+    const pointField& internalPoints
+)
+:
+    BSplineEdge(points, edge(from,to), internalPoints)
+{}
+
+
 Foam::blockEdges::BSplineEdge::BSplineEdge
 (
     const dictionary& dict,
     const label index,
-    const searchableSurfaces& geometry,
+    const searchableSurfaces&,
     const pointField& points,
     Istream& is
 )
@@ -78,7 +89,7 @@ Foam::blockEdges::BSplineEdge::BSplineEdge
     blockEdge(dict, index, points, is),
     BSpline
     (
-        polyLine::concat(points[start_], pointField(is), points[end_])
+        polyLine::concat(firstPoint(), pointField(is), lastPoint())
     )
 {
     token tok(is);
diff --git a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H
index 68eefdab87ec36b9c9cc18b7eea950a536c8edb0..193836853c1c35749ee9e89d6c92557c47bc3870 100644
--- a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H
+++ b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2014-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -78,8 +78,16 @@ public:
         BSplineEdge
         (
             const pointField& points,   //!< Referenced point field
-            const label start,      //!< Start point in referenced point field
-            const label end,        //!< End point in referenced point field
+            const edge& fromTo,         //!< Start/end in point field
+            const pointField& internalPoints
+        );
+
+        //- Construct from components
+        BSplineEdge
+        (
+            const pointField& points,   //!< Referenced point field
+            const label from,           //!< Start point in point field
+            const label to,             //!< End point in point field
             const pointField& internalPoints
         );
 
@@ -88,8 +96,8 @@ public:
         (
             const dictionary& dict,
             const label index,
-            const searchableSurfaces& geometry,
-            const pointField& points,   //!< Referenced point field
+            const searchableSurfaces&  /*unused*/,
+            const pointField& points,  //!< Referenced point field
             Istream& is
         );
 
diff --git a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C
index 4e5d33badd181cbc0dfdddcbbd2f0f57f2b4f1bd..035855e84438d0fe0d4cdfd0028d88b10248c948 100644
--- a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C
+++ b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.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.
@@ -193,41 +193,63 @@ Foam::blockEdges::arcEdge::arcEdge
 (
     const pointField& points,
     const point& origin,
-    const label start,
-    const label end
+    const edge& fromTo
 )
 :
-    blockEdge(points, start, end),
+    blockEdge(points, fromTo),
     radius_(0),
     angle_(0),
     cs_()
 {
-    calcFromCentre(points[start_], points[end_], origin);
+    calcFromCentre(firstPoint(), lastPoint(), origin);
 }
 
 
 Foam::blockEdges::arcEdge::arcEdge
 (
     const pointField& points,
-    const label start,
-    const label end,
+    const edge& fromTo,
     const point& midPoint
 )
 :
-    blockEdge(points, start, end),
+    blockEdge(points, fromTo),
     radius_(0),
     angle_(0),
     cs_()
 {
-    calcFromMidPoint(points[start_], points[end_], midPoint);
+    calcFromMidPoint(firstPoint(), lastPoint(), midPoint);
 }
 
 
+Foam::blockEdges::arcEdge::arcEdge
+(
+    const pointField& points,
+    const point& origin,
+    const label from,
+    const label to
+)
+:
+    arcEdge(points, origin, edge(from,to))
+{}
+
+
+Foam::blockEdges::arcEdge::arcEdge
+(
+    const pointField& points,
+    const label from,
+    const label to,
+    const point& midPoint
+)
+:
+    arcEdge(points, edge(from,to), midPoint)
+{}
+
+
 Foam::blockEdges::arcEdge::arcEdge
 (
     const dictionary& dict,
     const label index,
-    const searchableSurfaces& geometry,
+    const searchableSurfaces&,
     const pointField& points,
     Istream& is
 )
@@ -260,7 +282,7 @@ Foam::blockEdges::arcEdge::arcEdge
 
         is >> p;  // The origin (centre)
 
-        calcFromCentre(points_[start_], points_[end_], p, true, rMultiplier);
+        calcFromCentre(firstPoint(), lastPoint(), p, true, rMultiplier);
     }
     else
     {
@@ -268,7 +290,7 @@ Foam::blockEdges::arcEdge::arcEdge
 
         is >> p;  // A mid-point
 
-        calcFromMidPoint(points_[start_], points_[end_], p);
+        calcFromMidPoint(firstPoint(), lastPoint(), p);
     }
 
     if (debug)
@@ -295,11 +317,11 @@ Foam::point Foam::blockEdges::arcEdge::position(const scalar lambda) const
 
     if (lambda < SMALL)
     {
-        return points_[start_];
+        return firstPoint();
     }
     else if (lambda >= 1 - SMALL)
     {
-        return points_[end_];
+        return lastPoint();
     }
 
     return cs_.globalPosition(vector(radius_, (lambda*angle_), 0));
diff --git a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H
index 459b6a70f0d3f525c42ac06e8ec3d7b14b4f1d96..0a9d3c361fe827120ed15db1ecbf75b91f4bf531 100644
--- a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H
+++ b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.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.
@@ -136,18 +136,34 @@ public:
         arcEdge
         (
             const pointField& points,   //!< Referenced point field
-            const point& origin,    //!< The origin of the circle
-            const label start,      //!< Start point in referenced point field
-            const label end         //!< End point in referenced point field
+            const point& origin,        //!< The origin of the circle
+            const edge& fromTo          //!< Start/end in point field
         );
 
         //- Construct from components, using a point on the circumference
         arcEdge
         (
             const pointField& points,   //!< Referenced point field
-            const label start,      //!< Start point in referenced point field
-            const label end,        //!< End point in referenced point field
-            const point& midPoint   //!< A point on the circumference
+            const edge& fromTo,         //!< Start/end in point field
+            const point& midPoint       //!< A point on the circumference
+        );
+
+        //- Construct from components, given the origin of the circle
+        arcEdge
+        (
+            const pointField& points,   //!< Referenced point field
+            const point& origin,        //!< The origin of the circle
+            const label from,           //!< Start point in point field
+            const label to              //!< End point in point field
+        );
+
+        //- Construct from components, using a point on the circumference
+        arcEdge
+        (
+            const pointField& points,   //!< Referenced point field
+            const label from,           //!< Start point in point field
+            const label to,             //!< End point in point field
+            const point& midPoint       //!< A point on the circumference
         );
 
         //- Construct from Istream and point field.
@@ -157,8 +173,8 @@ public:
         (
             const dictionary& dict,
             const label index,
-            const searchableSurfaces& geometry, // unsed
-            const pointField& points,   //!< Referenced point field
+            const searchableSurfaces&  /*unused*/,
+            const pointField& points,  //!< Referenced point field
             Istream& is
         );
 
diff --git a/src/mesh/blockMesh/blockEdges/bezier/bezier.C b/src/mesh/blockMesh/blockEdges/bezier/bezier.C
index cac55d1ba3bfd3e246ae8b87dbc8238ebaa3a8e0..928f9e8917fcf921296d46f338e14f7397f596e6 100644
--- a/src/mesh/blockMesh/blockEdges/bezier/bezier.C
+++ b/src/mesh/blockMesh/blockEdges/bezier/bezier.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2019 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,21 +47,32 @@ namespace blockEdges
 Foam::blockEdges::bezier::bezier
 (
     const pointField& points,
-    const label start,
-    const label end,
+    const edge& fromTo,
     const pointField& control
 )
 :
-    blockEdge(points, start, end),
+    blockEdge(points, fromTo),
     control_(control)
 {}
 
 
+Foam::blockEdges::bezier::bezier
+(
+    const pointField& points,
+    const label from,
+    const label to,
+    const pointField& control
+)
+:
+    bezier(points, edge(from, to), control)
+{}
+
+
 Foam::blockEdges::bezier::bezier
 (
     const dictionary& dict,
     const label index,
-    const searchableSurfaces& geometry,
+    const searchableSurfaces&,
     const pointField& points,
     Istream& is
 )
@@ -69,7 +80,7 @@ Foam::blockEdges::bezier::bezier
     blockEdge(dict, index, points, is),
     control_
     (
-        polyLine::concat(points[start_], pointField(is), points[end_])
+        polyLine::concat(firstPoint(), pointField(is), lastPoint())
     )
 {}
 
diff --git a/src/mesh/blockMesh/blockEdges/bezier/bezier.H b/src/mesh/blockMesh/blockEdges/bezier/bezier.H
index 5feaeb4139205511df37f97fe69f8db89a726bea..4a7d2780f5fd32d216c188653d3993bc33745a30 100644
--- a/src/mesh/blockMesh/blockEdges/bezier/bezier.H
+++ b/src/mesh/blockMesh/blockEdges/bezier/bezier.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2019 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -90,8 +90,16 @@ public:
         bezier
         (
             const pointField& points,   //!< Referenced point field
-            const label start,      //!< Start point in referenced point field
-            const label end,        //!< End point in referenced point field
+            const edge& fromTo,         //!< Start/end in point field
+            const pointField& control   //!< The control points
+        );
+
+        //- Construct from components
+        bezier
+        (
+            const pointField& points,   //!< Referenced point field
+            const label from,           //!< Start point in point field
+            const label to,             //!< End point in point field
             const pointField& control   //!< The control points
         );
 
@@ -100,8 +108,8 @@ public:
         (
             const dictionary& dict,
             const label index,
-            const searchableSurfaces& geometry,
-            const pointField& points,   //!< Referenced point field
+            const searchableSurfaces&  /*unused*/,
+            const pointField& points,  //!< Referenced point field
             Istream& is
         );
 
diff --git a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C
index 871b064cbfb1b748cb3b576d80940df3a53ccc2d..81033557d2f6c4b172e1372120b0aba0152f2f93 100644
--- a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C
+++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C
@@ -44,13 +44,12 @@ namespace Foam
 Foam::blockEdge::blockEdge
 (
     const pointField& points,
-    const label start,
-    const label end
+    const edge& fromTo
 )
 :
     points_(points),
-    start_(start),
-    end_(end)
+    start_(fromTo.first()),
+    end_(fromTo.last())
 {}
 
 
@@ -109,13 +108,13 @@ Foam::autoPtr<Foam::blockEdge> Foam::blockEdge::New
 
 Foam::pointField Foam::blockEdge::appendEndPoints
 (
-    const pointField& pts,
-    const label start,
-    const label end,
+    const pointField& p,
+    const label from,
+    const label to,
     const pointField& intermediate
 )
 {
-    return pointField(polyLine::concat(pts[start], intermediate, pts[end]));
+    return pointField(polyLine::concat(p[from], intermediate, p[to]));
 }
 
 
diff --git a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H
index ae87a2fc5bdbcfa7332a9a25b98e441ef61c230f..1c6a51fb84b09d52ed5b4dcc849ae00ee9fef3f7 100644
--- a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H
+++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,6 +46,7 @@ SourceFiles
 #define blockEdge_H
 
 #include "searchableSurfaces.H"
+#include "edge.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,13 +70,15 @@ protected:
         //- The referenced point field
         const pointField& points_;
 
-        //- Index of the start point
+        //- Index of the first point
         const label start_;
 
-        //- Index of the end point
+        //- Index of the last point
         const label end_;
 
 
+protected:
+
     // Protected Member Functions
 
         //- Return a complete point field by appending the start/end points
@@ -83,12 +86,23 @@ protected:
         //  \deprecated(2020-10) use polyLine::concat
         static pointField appendEndPoints
         (
-            const pointField& points,   //!< Referenced point field
-            const label start,      //!< Start point in referenced point field
-            const label end,        //!< End point in referenced point field
+            const pointField& p,        //!< Referenced point field
+            const label from,           //!< Start point in point field
+            const label to,             //!< End point in point field
             const pointField& intermediate  //!< Intermediate points (knots)
         );
 
+        //- Construct from components
+        //  \deprecated(2021-11) use constructor with edge
+        blockEdge
+        (
+            const pointField& points,   //!< Referenced point field
+            const label from,           //!< Start point in point field
+            const label to              //!< End point in point field
+        )
+        :
+            blockEdge(points, edge(from,to))
+        {}
 
 public:
 
@@ -119,8 +133,7 @@ public:
         blockEdge
         (
             const pointField& points,   //!< Referenced point field
-            const label start,      //!< Start point in referenced point field
-            const label end         //!< End point in referenced point field
+            const edge& fromTo          //!< Start/end in point field
         );
 
         //- Construct from Istream and point field.
@@ -182,11 +195,21 @@ public:
 
     // Member Functions
 
-        //- Index of start point
-        inline label start() const;
+        //- True if first/last indices are unique and non-negative.
+        inline bool valid() const noexcept;
+
+        //- Index of start (first) point
+        inline label start() const noexcept;
+
+        //- Index of end (last) point
+        inline label end() const noexcept;
+
+        //- The location of the first point
+        inline const point& firstPoint() const;
+
+        //- The location of the last point
+        inline const point& lastPoint() const;
 
-        //- Index of end point
-        inline label end() const;
 
         //- Compare the given start/end points with this block edge
         //  Return:
@@ -209,6 +232,11 @@ public:
         //  - -1: same edge, but different orientation
         inline int compare(const label start, const label end) const;
 
+
+        //- The point position in the straight line
+        //  0 <= lambda <= 1
+        inline point linearPosition(const scalar lambda) const;
+
         //- The point position corresponding to the curve parameter
         //  0 <= lambda <= 1
         virtual point position(const scalar lambda) const = 0;
diff --git a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdgeI.H b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdgeI.H
index 6aab3064f06d052b2d49d07a294b6afbc9546c2a..2f781b8785162b5574ef4d1116b677017022a50a 100644
--- a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdgeI.H
+++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdgeI.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,18 +28,36 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::label Foam::blockEdge::start() const
+inline bool Foam::blockEdge::valid() const noexcept
+{
+    return (start_ != end_ && start_ >= 0 && end_ >= 0);
+}
+
+
+inline Foam::label Foam::blockEdge::start() const noexcept
 {
     return start_;
 }
 
 
-inline Foam::label Foam::blockEdge::end() const
+inline Foam::label Foam::blockEdge::end() const noexcept
 {
     return end_;
 }
 
 
+inline const Foam::point& Foam::blockEdge::firstPoint() const
+{
+    return points_[start_];
+}
+
+
+inline const Foam::point& Foam::blockEdge::lastPoint() const
+{
+    return points_[end_];
+}
+
+
 inline int Foam::blockEdge::compare(const label start, const label end) const
 {
     if (start_ == start && end_ == end)
@@ -66,4 +85,27 @@ inline int Foam::blockEdge::compare(const edge& e) const
 }
 
 
+inline Foam::point Foam::blockEdge::linearPosition(const scalar lambda) const
+{
+    #ifdef FULLDEBUG
+    if (lambda < -SMALL || lambda > 1 + SMALL)
+    {
+        InfoInFunction
+            << "Limit parameter to [0-1] range: " << lambda << nl;
+    }
+    #endif
+
+    if (lambda < SMALL)
+    {
+        return firstPoint();
+    }
+    else if (lambda >= 1 - SMALL)
+    {
+        return lastPoint();
+    }
+
+    return firstPoint() + lambda * (lastPoint() - firstPoint());
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C b/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C
index fde57357c6db325fbed45bf3a48750a8d6ca1016..0ff8bdc93d25f67eb09ed06603538e6822fdbb6c 100644
--- a/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C
+++ b/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -139,13 +140,13 @@ Foam::lineDivide::lineDivide
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::pointField& Foam::lineDivide::points() const
+const Foam::pointField& Foam::lineDivide::points() const noexcept
 {
     return points_;
 }
 
 
-const Foam::scalarList& Foam::lineDivide::lambdaDivisions() const
+const Foam::scalarList& Foam::lineDivide::lambdaDivisions() const noexcept
 {
     return divisions_;
 }
diff --git a/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.H b/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.H
index 1c58386a0ffc62093d3bd3144071a81e020f0caf..16c99fe0ecf524a531ce0c8515848b257892d11c 100644
--- a/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.H
+++ b/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.H
@@ -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.
@@ -78,10 +78,10 @@ public:
     // Member Functions
 
         //- The points
-        const pointField& points() const;
+        const pointField& points() const noexcept;
 
         //- The list of lambda values
-        const scalarList& lambdaDivisions() const;
+        const scalarList& lambdaDivisions() const noexcept;
 };
 
 
diff --git a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C
index 564240197c16a233008e27afa000687ee0114460..af50d632426725feeb527449d86158dcda9953ce 100644
--- a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C
+++ b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,11 +46,21 @@ namespace blockEdges
 Foam::blockEdges::lineEdge::lineEdge
 (
     const pointField& points,
-    const label start,
-    const label end
+    const edge& fromTo
 )
 :
-    blockEdge(points, start, end)
+    blockEdge(points, fromTo)
+{}
+
+
+Foam::blockEdges::lineEdge::lineEdge
+(
+    const pointField& points,
+    const label from,
+    const label to
+)
+:
+    blockEdge(points, from, to)
 {}
 
 
@@ -58,7 +68,7 @@ Foam::blockEdges::lineEdge::lineEdge
 (
     const dictionary& dict,
     const label index,
-    const searchableSurfaces& geometry,
+    const searchableSurfaces&,
     const pointField& points,
     Istream& is
 )
@@ -71,30 +81,13 @@ Foam::blockEdges::lineEdge::lineEdge
 
 Foam::point Foam::blockEdges::lineEdge::position(const scalar lambda) const
 {
-    #ifdef FULLDEBUG
-    if (lambda < -SMALL || lambda > 1 + SMALL)
-    {
-        InfoInFunction
-            << "Limit parameter to [0-1] range: " << lambda << nl;
-    }
-    #endif
-
-    if (lambda < SMALL)
-    {
-        return points_[start_];
-    }
-    else if (lambda >= 1 - SMALL)
-    {
-        return points_[end_];
-    }
-
-    return points_[start_] + lambda * (points_[end_] - points_[start_]);
+    return blockEdge::linearPosition(lambda);
 }
 
 
 Foam::scalar Foam::blockEdges::lineEdge::length() const
 {
-    return mag(points_[end_] - points_[start_]);
+    return Foam::mag(lastPoint() - firstPoint());
 }
 
 
diff --git a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H
index 1bbc9d67ce5b3d6f24bcba2aa4f27488830b3a33..b3c2b208dea0a43b9c933873c933ee994d829214 100644
--- a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H
+++ b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -67,8 +67,15 @@ public:
         lineEdge
         (
             const pointField& points,   //!< Referenced point field
-            const label start,      //!< Start point in referenced point field
-            const label end         //!< End point in referenced point field
+            const edge& fromTo          //!< Start/end in point field
+        );
+
+        //- Construct from components
+        lineEdge
+        (
+            const pointField& points,   //!< Referenced point field
+            const label from,           //!< Start point in point field
+            const label to              //!< End point in point field
         );
 
         //- Construct from Istream and point field.
@@ -76,8 +83,8 @@ public:
         (
             const dictionary& dict,
             const label index,
-            const searchableSurfaces& geometry,
-            const pointField& points,   //!< Referenced point field
+            const searchableSurfaces&  /*unused*/,
+            const pointField& points,  //!< Referenced point field
             Istream& is
         );
 
diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.C b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.C
index 42f6a74db691bebd9fac04b3fda427d3f7853c9c..e9ebf872e76b24ad51b3c5427ed3f2ef6faef89d 100644
--- a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.C
+++ b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.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.
@@ -81,9 +81,9 @@ void Foam::polyLine::calcParam()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::polyLine::polyLine(const pointField& ps, const bool)
+Foam::polyLine::polyLine(const pointField& p, const bool)
 :
-    points_(ps),
+    points_(p),
     lineLength_(0),
     param_()
 {
@@ -109,13 +109,13 @@ Foam::polyLine::polyLine
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::pointField& Foam::polyLine::points() const
+const Foam::pointField& Foam::polyLine::points() const noexcept
 {
     return points_;
 }
 
 
-Foam::label Foam::polyLine::nSegments() const
+Foam::label Foam::polyLine::nSegments() const noexcept
 {
     return points_.size()-1;
 }
@@ -207,7 +207,7 @@ Foam::point Foam::polyLine::position
 }
 
 
-Foam::scalar Foam::polyLine::length() const
+Foam::scalar Foam::polyLine::length() const noexcept
 {
     return lineLength_;
 }
diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.H b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.H
index b7bf40e8afc424b538418ee179d9649b35ba6c08..48ce09a90070036231744215a4f405c955d765e0 100644
--- a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.H
+++ b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.H
@@ -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.
@@ -115,10 +115,10 @@ public:
     // Member Functions
 
         //- Return const-access to the control-points
-        const pointField& points() const;
+        const pointField& points() const noexcept;
 
         //- The number of line segments
-        label nSegments() const;
+        label nSegments() const noexcept;
 
         //- The point position corresponding to the curve parameter
         //  0 <= lambda <= 1
@@ -129,7 +129,7 @@ public:
         point position(const label segment, const scalar) const;
 
         //- The length of the curve
-        scalar length() const;
+        scalar length() const noexcept;
 };
 
 
diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C
index 9a1b014fee0fe3517953d7648bd613515c10ea8d..393d589b6cb8ad171b1cf5a65b95cfd91c34b7b6 100644
--- a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C
+++ b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,14 +45,25 @@ namespace blockEdges
 
 Foam::blockEdges::polyLineEdge::polyLineEdge
 (
-    const pointField& ps,
-    const label start,
-    const label end,
+    const pointField& points,
+    const edge& fromTo,
     const pointField& intermediate
 )
 :
-    blockEdge(ps, start, end),
-    polyLine(ps[start_], intermediate, ps[end_])
+    blockEdge(points, fromTo),
+    polyLine(firstPoint(), intermediate, lastPoint())
+{}
+
+
+Foam::blockEdges::polyLineEdge::polyLineEdge
+(
+    const pointField& points,
+    const label from,
+    const label to,
+    const pointField& intermediate
+)
+:
+    polyLineEdge(points, edge(from, to), intermediate)
 {}
 
 
@@ -60,13 +71,13 @@ Foam::blockEdges::polyLineEdge::polyLineEdge
 (
     const dictionary& dict,
     const label index,
-    const searchableSurfaces& geometry,
-    const pointField& ps,
+    const searchableSurfaces&,
+    const pointField& points,
     Istream& is
 )
 :
-    blockEdge(dict, index, ps, is),
-    polyLine(ps[start_], pointField(is), ps[end_])
+    blockEdge(dict, index, points, is),
+    polyLine(firstPoint(), pointField(is), lastPoint())
 {}
 
 
diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H
index bcce9c738e2599aaafc3655caa3e49b1f79a4dbf..2fbadf3a29f2eedb3737e83265d0288d897620b3 100644
--- a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H
+++ b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -78,8 +78,16 @@ public:
         polyLineEdge
         (
             const pointField& points,       //!< Referenced point field
-            const label start,      //!< Start point in referenced point field
-            const label end,        //!< End point in referenced point field
+            const edge& fromTo,             //!< Start/end in point field
+            const pointField& intermediate  //!< The intermediate points
+        );
+
+        //- Construct from components
+        polyLineEdge
+        (
+            const pointField& points,       //!< Referenced point field
+            const label from,               //!< Start point in point field
+            const label to,                 //!< End point in point field
             const pointField& intermediate  //!< The intermediate points
         );
 
@@ -88,8 +96,8 @@ public:
         (
             const dictionary& dict,
             const label index,
-            const searchableSurfaces& geometry,
-            const pointField& points,       //!< Referenced point field
+            const searchableSurfaces&  /*unused*/,
+            const pointField& points,  //!< Referenced point field
             Istream& is
         );
 
diff --git a/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C
index da862fd7d705097792b145a0b9d56ae48207e906..6acef8ca8618f40290983b1588c508a04939a682 100644
--- a/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C
+++ b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C
@@ -115,14 +115,12 @@ Foam::blockEdges::projectCurveEdge::position(const scalarList& lambdas) const
     auto tpoints = tmp<pointField>::New(lambdas.size());
     auto& points = tpoints.ref();
 
-    const point& startPt = points_[start_];
-    const point& endPt = points_[end_];
-    const vector d = endPt-startPt;
+    const scalar distSqr = Foam::magSqr(lastPoint()-firstPoint());
 
     // Initial guess
     forAll(lambdas, i)
     {
-        points[i] = startPt+lambdas[i]*d;
+        points[i] = blockEdge::linearPosition(lambdas[i]);
     }
 
     // Use special interpolation to keep initial guess on same position on
@@ -142,7 +140,7 @@ Foam::blockEdges::projectCurveEdge::position(const scalarList& lambdas) const
                 points[0],
                 points.last(),
                 scalarField(lambdas),
-                scalarField(points.size(), magSqr(d)),
+                scalarField(points.size(), distSqr),
                 nearInfo
             );
             forAll(nearInfo, i)
@@ -179,7 +177,7 @@ Foam::blockEdges::projectCurveEdge::position(const scalarList& lambdas) const
                 geometry_,
                 surfaces_,
                 start,
-                scalarField(start.size(), magSqr(d)),
+                scalarField(start.size(), distSqr),
                 points,
                 constraints
             );
@@ -187,11 +185,11 @@ Foam::blockEdges::projectCurveEdge::position(const scalarList& lambdas) const
             // Reset start and end point
             if (lambdas[0] < SMALL)
             {
-                points[0] = startPt;
+                points[0] = firstPoint();
             }
             if (lambdas.last() > 1.0-SMALL)
             {
-                points.last() = endPt;
+                points.last() = lastPoint();
             }
 
             if (debugStr)
diff --git a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C
index 40ce7fe178d98ac24fbf12864008f54d7dca0101..96903d13e61591d8a517f35e64dfd89b876f355f 100644
--- a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C
+++ b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,7 +55,7 @@ void Foam::blockEdges::projectEdge::findNearest
 {
     if (surfaces_.size())
     {
-        const scalar distSqr = magSqr(points_[end_]-points_[start_]);
+        const scalar distSqr = Foam::magSqr(lastPoint()-firstPoint());
 
         pointField boundaryNear(1);
         List<pointConstraint> boundaryConstraint(1);
@@ -95,7 +95,7 @@ Foam::blockEdges::projectEdge::projectEdge
     geometry_(geometry)
 {
     wordList names(is);
-    surfaces_.setSize(names.size());
+    surfaces_.resize(names.size());
     forAll(names, i)
     {
         surfaces_[i] = geometry_.findSurfaceID(names[i]);
@@ -115,7 +115,7 @@ Foam::blockEdges::projectEdge::projectEdge
 Foam::point Foam::blockEdges::projectEdge::position(const scalar lambda) const
 {
     // Initial guess
-    const point start(points_[start_] + lambda*(points_[end_]-points_[start_]));
+    const point start(blockEdge::linearPosition(lambda));
 
     point near(start);
 
@@ -150,14 +150,12 @@ Foam::blockEdges::projectEdge::position(const scalarList& lambdas) const
     auto tpoints = tmp<pointField>::New(lambdas.size());
     auto& points = tpoints.ref();
 
-    const point& startPt = points_[start_];
-    const point& endPt = points_[end_];
-    const vector d = endPt-startPt;
+    const scalar distSqr = Foam::magSqr(lastPoint()-firstPoint());
 
     // Initial guess
     forAll(lambdas, i)
     {
-        points[i] = startPt+lambdas[i]*d;
+        points[i] = blockEdge::linearPosition(lambdas[i]);
     }
 
 
@@ -181,7 +179,7 @@ Foam::blockEdges::projectEdge::position(const scalarList& lambdas) const
                 geometry_,
                 surfaces_,
                 start,
-                scalarField(start.size(), magSqr(d)),
+                scalarField(start.size(), distSqr),
                 points,
                 constraints
             );
@@ -189,11 +187,11 @@ Foam::blockEdges::projectEdge::position(const scalarList& lambdas) const
             // Reset start and end point
             if (lambdas[0] < SMALL)
             {
-                points[0] = startPt;
+                points[0] = firstPoint();
             }
             if (lambdas.last() > 1.0-SMALL)
             {
-                points.last() = endPt;
+                points.last() = lastPoint();
             }
 
             if (debugStr)
diff --git a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C
index 5494d4ed884fe19464b05c0dd10a11bdc2f8a380..124a740362cc148269904965cff96e0dba7498b7 100644
--- a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C
+++ b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -53,24 +53,35 @@ namespace blockEdges
 Foam::blockEdges::splineEdge::splineEdge
 (
     const pointField& points,
-    const label start,
-    const label end,
+    const edge& fromTo,
     const pointField& internalPoints
 )
 :
-    blockEdge(points, start, end),
+    blockEdge(points, fromTo),
     CatmullRomSpline
     (
-        polyLine::concat(points[start_], internalPoints, points[end_])
+        polyLine::concat(firstPoint(), internalPoints, lastPoint())
     )
 {}
 
 
+Foam::blockEdges::splineEdge::splineEdge
+(
+    const pointField& points,
+    const label from,
+    const label to,
+    const pointField& internalPoints
+)
+:
+    splineEdge(points, edge(from,to), internalPoints)
+{}
+
+
 Foam::blockEdges::splineEdge::splineEdge
 (
     const dictionary& dict,
     const label index,
-    const searchableSurfaces& geometry,
+    const searchableSurfaces&,
     const pointField& points,
     Istream& is
 )
@@ -78,7 +89,7 @@ Foam::blockEdges::splineEdge::splineEdge
     blockEdge(dict, index, points, is),
     CatmullRomSpline
     (
-        polyLine::concat(points[start_], pointField(is), points[end_])
+        polyLine::concat(firstPoint(), pointField(is), lastPoint())
     )
 {
     token tok(is);
diff --git a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H
index 3ac595cbac7fc0b87e91ad72cbb475f743a62ba9..b138c7db21b64a4770ae72c4d77547b51e015d2f 100644
--- a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H
+++ b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -82,8 +82,16 @@ public:
         splineEdge
         (
             const pointField& points,   //!< Referenced point field
-            const label start,      //!< Start point in referenced point field
-            const label end,        //!< End point in referenced point field
+            const edge& fromTo,         //!< Start/end in point field
+            const pointField& internalPoints
+        );
+
+        //- Construct from components
+        splineEdge
+        (
+            const pointField& points,   //!< Referenced point field
+            const label from,           //!< Start point in point field
+            const label to,             //!< End point in point field
             const pointField& internalPoints
         );
 
@@ -92,8 +100,8 @@ public:
         (
             const dictionary& dict,
             const label index,
-            const searchableSurfaces& geometry,
-            const pointField& points,   //!< Referenced point field
+            const searchableSurfaces&  /*unused*/,
+            const pointField& points,  //!< Referenced point field
             Istream&
         );
 
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.C b/src/mesh/blockMesh/blockMesh/blockMesh.C
index 18321776522b13b19363805140f8445918aa68d0..f9777d83c35e8b66f02d99962bde1a9287fab9ad 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.C
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.C
@@ -230,11 +230,12 @@ Foam::blockMesh::blockMesh
 (
     const IOdictionary& dict,
     const word& regionName,
-    mergeStrategy strategy
+    mergeStrategy strategy,
+    int verbosity
 )
 :
     meshDict_(dict),
-    verbose_(meshDict_.getOrDefault("verbose", verboseOutput)),
+    verbose_(verbosity),
     checkFaceCorrespondence_
     (
         meshDict_.getOrDefault("checkFaceCorrespondence", true)
@@ -268,13 +269,32 @@ Foam::blockMesh::blockMesh
     transform_(),
     topologyPtr_(createTopology(meshDict_, regionName))
 {
-    // Command-line option has precedence over dictionary setting
+    // Command-line options have precedence over dictionary setting
+
+    if (!verbose_)
+    {
+        verbose_ = meshDict_.getOrDefault("verbose", verboseOutput);
+    }
 
     if (mergeStrategy_ == mergeStrategy::DEFAULT_MERGE)
     {
         strategyNames_.readIfPresent("mergeType", meshDict_, mergeStrategy_);
 
         // Warn about fairly obscure old "fastMerge" option?
+
+        if
+        (
+            mergeStrategy_ == mergeStrategy::DEFAULT_MERGE
+         && checkDegenerate()
+        )
+        {
+            Info<< nl
+                << "Detected collapsed blocks "
+                << "- using merge points instead of merge topology" << nl
+                << endl;
+
+            mergeStrategy_ = mergeStrategy::MERGE_POINTS;
+        }
     }
 
     if (mergeStrategy_ == mergeStrategy::MERGE_POINTS)
@@ -298,16 +318,16 @@ bool Foam::blockMesh::valid() const noexcept
 }
 
 
-bool Foam::blockMesh::verbose() const noexcept
+int Foam::blockMesh::verbose() const noexcept
 {
     return verbose_;
 }
 
 
-bool Foam::blockMesh::verbose(const bool on) noexcept
+int Foam::blockMesh::verbose(const int level) noexcept
 {
-    bool old(verbose_);
-    verbose_ = on;
+    int old(verbose_);
+    verbose_ = level;
     return old;
 }
 
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.H b/src/mesh/blockMesh/blockMesh/blockMesh.H
index df613d58063f7e963027ab9fd65b23d8ed57e04b..c8cb1c1e2939d36c2fb17231f9568ced35cfdd14 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.H
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.H
@@ -68,6 +68,7 @@ SourceFiles
 #define blockMesh_H
 
 #include "Enum.H"
+#include "Switch.H"
 #include "block.H"
 #include "PtrList.H"
 #include "cartesianCS.H"
@@ -138,8 +139,8 @@ private:
         //- Reference to mesh dictionary
         const IOdictionary& meshDict_;
 
-        //- Output verbosity
-        bool verbose_;
+        //- Output verbosity level
+        int verbose_;
 
         //- Check face consistency (defaults to true)
         bool checkFaceCorrespondence_;
@@ -228,6 +229,10 @@ private:
             const word& regionName
         );
 
+
+        //- Simple checks for collapsed hex cells
+        bool checkDegenerate() const;
+
         void check(const polyMesh& bm, const dictionary& dict) const;
 
         //- Determine merge info and final number of cells/points
@@ -272,7 +277,8 @@ public:
         (
             const IOdictionary& dict,
             const word& regionName = polyMesh::defaultRegion,
-            mergeStrategy strategy = mergeStrategy::DEFAULT_MERGE
+            mergeStrategy strategy = mergeStrategy::DEFAULT_MERGE,
+            int verbosity = 0   // 0: use static or dictionary value
         );
 
 
@@ -367,12 +373,12 @@ public:
 
     // Verbosity
 
-        //- Verbose output
-        bool verbose() const noexcept;
+        //- Output verbosity level
+        int verbose() const noexcept;
 
-        //- Enable/disable verbose output
-        //  \return old value
-        bool verbose(const bool on) noexcept;
+        //- Change the output verbosity level.
+        //  \return old level
+        int verbose(const int level) noexcept;
 
 
     // Mesh Generation
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshCheck.C b/src/mesh/blockMesh/blockMesh/blockMeshCheck.C
index 086c4ceb5c7a1e9de0f7e631e35ac3e97ef82fd2..922aaa331c070c62a4a20876fbcadc1b7902d6be 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshCheck.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshCheck.C
@@ -27,10 +27,43 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "blockMesh.H"
+#include "IOmanip.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
+bool Foam::blockMesh::checkDegenerate() const
+{
+    const blockList& blocks = *this;
+
+    for (const block& blk : blocks)
+    {
+        const cellShape& shape = blk.blockShape();
+
+        if (shape.model().index() == cellModel::HEX)
+        {
+            // Check for collapsed edges
+            // - limit to HEX only for now.
+            for (label edgei = 0; edgei < shape.nEdges(); ++edgei)
+            {
+                edge e(shape.edge(edgei));
+
+                if (!e.valid())
+                {
+                    return true;  // Looks like a collapsed edge
+                }
+            }
+        }
+    }
+
+    return false;
+}
+
+
+void Foam::blockMesh::check
+(
+    const polyMesh& bm,
+    const dictionary& dict
+) const
 {
     Info<< nl << "Check topology" << endl;
 
@@ -45,8 +78,8 @@ void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
             {
                 Info<< "    Curved edge ";
                 edges_[cej].write(Info, dict);
-                Info<< "    is a duplicate of curved edge " << edges_[cei]
-                    << endl;
+                Info<< "    is a duplicate of curved edge "
+                    << edges_[cei] << endl;
                 ok = false;
                 break;
             }
@@ -60,17 +93,15 @@ void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
     // repeated point labels
     const blockList& blocks = *this;
 
-    forAll(edges_, cei)
+    for (const blockEdge& curvEdge : edges_)
     {
         bool found = false;
 
-        forAll(blocks, blocki)
+        for (const block& blk : blocks)
         {
-            edgeList edges = blocks[blocki].blockShape().edges();
-
-            forAll(edges, ei)
+            for (const edge& blkEdge : blk.blockShape().edges())
             {
-                found = edges_[cei].compare(edges[ei][0], edges[ei][1]) != 0;
+                found = curvEdge.compare(blkEdge) != 0;
                 if (found) break;
             }
             if (found) break;
@@ -79,7 +110,7 @@ void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
         if (!found)
         {
             Info<< "    Curved edge ";
-            edges_[cei].write(Info, dict);
+            curvEdge.write(Info, dict);
             Info<< "    does not correspond to a block edge." << endl;
             ok = false;
         }
@@ -142,22 +173,22 @@ void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
         }
     }
 
+
     const pointField& points = bm.points();
     const cellList& cells = bm.cells();
     const polyPatchList& patches = bm.boundaryMesh();
 
     label nBoundaryFaces = 0;
-    forAll(cells, celli)
+    for (const cell& c : cells)
     {
-        nBoundaryFaces += cells[celli].nFaces();
+        nBoundaryFaces += c.nFaces();
     }
-
     nBoundaryFaces -= 2*bm.nInternalFaces();
 
     label nDefinedBoundaryFaces = 0;
-    forAll(patches, patchi)
+    for (const polyPatch& pp : patches)
     {
-        nDefinedBoundaryFaces += patches[patchi].size();
+        nDefinedBoundaryFaces += pp.size();
     }
 
 
@@ -184,22 +215,21 @@ void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
     }
 
 
-    forAll(patches, patchi)
+    for (const polyPatch& pp : patches)
     {
-        const faceList& Patch = patches[patchi];
-
-        forAll(Patch, patchFacei)
+        forAll(pp, patchFacei)
         {
-            const face& patchFace = Patch[patchFacei];
+            const face& patchFace = pp[patchFacei];
+
             bool patchFaceOK = false;
 
-            forAll(cells, celli)
+            for (const labelList& cellFaces : cells)
             {
-                const labelList& cellFaces = cells[celli];
-
-                forAll(cellFaces, cellFacei)
+                for (const label cellFacei : cellFaces)
                 {
-                    if (patchFace == faces[cellFaces[cellFacei]])
+                    const face& cellFace = faces[cellFacei];
+
+                    if (patchFace == cellFace)
                     {
                         patchFaceOK = true;
 
@@ -207,14 +237,14 @@ void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
                         (
                             (
                                 patchFace.areaNormal(points)
-                              & faces[cellFaces[cellFacei]].areaNormal(points)
+                              & cellFace.areaNormal(points)
                             ) < 0
                         )
                         {
                             Info<< tab << tab
                                 << "Face " << patchFacei
-                                << " of patch " << patchi
-                                << " (" << patches[patchi].name() << ")"
+                                << " of patch " << pp.index()
+                                << " (" << pp.name() << ')'
                                 << " points inwards"
                                 << endl;
 
@@ -228,8 +258,8 @@ void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
             {
                 Info<< tab << tab
                     << "Face " << patchFacei
-                    << " of patch " << patchi
-                    << " (" << patches[patchi].name() << ")"
+                    << " of patch " << pp.index()
+                    << " (" << pp.name() << ')'
                     << " does not match any block faces" << endl;
 
                 ok = false;
@@ -242,6 +272,42 @@ void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
         Info<< endl;
     }
 
+    // Report patch block/face correspondence
+    if (verbose_ > 1)
+    {
+        const labelList& own = bm.faceOwner();
+
+        Info.stream().setf(ios_base::left);
+
+        Info<< setw(20) << "patch" << "block/face" << nl
+            << setw(20) << "-----" << "----------" << nl;
+
+        for (const polyPatch& pp : patches)
+        {
+            Info<< setw(20) << pp.name();
+
+            label meshFacei = pp.start();
+
+            forAll(pp, bfacei)
+            {
+                const label celli = own[meshFacei];
+                const label cellFacei = cells[celli].find(meshFacei);
+
+                if (bfacei) Info<< token::SPACE;
+
+                Info<< token::BEGIN_LIST
+                    << celli << ' ' << cellFacei
+                    << token::END_LIST;
+
+                ++meshFacei;
+            }
+            Info<< nl;
+        }
+
+        Info<< setw(20) << "-----" << "----------" << nl
+            << nl;
+    }
+
     if (!ok)
     {
         FatalErrorInFunction