diff --git a/applications/utilities/mesh/generation/blockMesh/addCellZones.H b/applications/utilities/mesh/generation/blockMesh/addCellZones.H
index eca57f2ed4..7cec1b4541 100644
--- a/applications/utilities/mesh/generation/blockMesh/addCellZones.H
+++ b/applications/utilities/mesh/generation/blockMesh/addCellZones.H
@@ -1,87 +1 @@
-// Set any cellZones
-// Note cell labelling unaffected by previous mergePatchPairs
-
-{
-    const label nZones = blocks.numZonedBlocks();
-    if (nZones)
-    {
-        Info<< nl << "Adding cell zones" << endl;
-
-        // Map from zoneName to cellZone index
-        HashTable<label> zoneMap(nZones);
-
-        // Cells per zone.
-        List<DynamicList<label>> zoneCells(nZones);
-
-        // Running cell counter
-        label celli = 0;
-
-        // Largest zone so far
-        label freeZoneI = 0;
-
-        for (const block& b : blocks)
-        {
-            const word& zoneName = b.zoneName();
-            const label nCellsInBlock = b.cells().size();
-
-            if (zoneName.size())
-            {
-                const auto iter = zoneMap.cfind(zoneName);
-
-                label zonei;
-
-                if (iter.found())
-                {
-                    zonei = *iter;
-                }
-                else
-                {
-                    zonei = freeZoneI++;
-
-                    Info<< "    " << zonei << '\t' << zoneName << endl;
-
-                    zoneMap.insert(zoneName, zonei);
-                }
-
-
-                // Fill with cell ids
-
-                zoneCells[zonei].reserve
-                (
-                    zoneCells[zonei].size() + nCellsInBlock
-                );
-
-                const label endOfFill = celli + nCellsInBlock;
-
-                for (; celli < endOfFill; ++celli)
-                {
-                    zoneCells[zonei].append(celli);
-                }
-            }
-            else
-            {
-                celli += nCellsInBlock;
-            }
-        }
-
-        List<cellZone*> cz(zoneMap.size());
-        forAllConstIters(zoneMap, iter)
-        {
-            const word& zoneName = iter.key();
-            const label zonei = iter.val();
-
-            cz[zonei] = new cellZone
-            (
-                zoneName,
-                zoneCells[zonei].shrink(),
-                zonei,
-                mesh.cellZones()
-            );
-        }
-
-        mesh.pointZones().resize(0);
-        mesh.faceZones().resize(0);
-        mesh.cellZones().resize(0);
-        mesh.addZones(List<pointZone*>(), List<faceZone*>(), cz);
-    }
-}
+#warning File removed - left for old dependency check only
diff --git a/applications/utilities/mesh/generation/blockMesh/blockMesh.C b/applications/utilities/mesh/generation/blockMesh/blockMesh.C
index 0eec9f01b8..b510d91170 100644
--- a/applications/utilities/mesh/generation/blockMesh/blockMesh.C
+++ b/applications/utilities/mesh/generation/blockMesh/blockMesh.C
@@ -77,7 +77,6 @@ Usage
 #include "foamVtkInternalMeshWriter.H"
 #include "attachPolyTopoChanger.H"
 #include "polyTopoChange.H"
-#include "emptyPolyPatch.H"
 #include "cyclicPolyPatch.H"
 #include "cellSet.H"
 
@@ -85,7 +84,7 @@ Usage
 #include "OSspecific.H"
 #include "OFstream.H"
 
-#include "Pair.H"
+#include "wordPair.H"
 #include "slidingInterface.H"
 
 using namespace Foam;
@@ -122,7 +121,8 @@ int main(int argc, char *argv[])
     argList::addBoolOption
     (
         "write-obj",
-        "Write block edges and centres as obj files and exit"
+        "Write block edges and centres as obj files and exit",
+        true  // (old) mark as advanced option. -write-vtk is preferred
     );
     argList::addOptionCompat("write-obj", {"blockTopology", 1912});
 
@@ -135,9 +135,9 @@ int main(int argc, char *argv[])
     argList::addBoolOption
     (
         "merge-points",
-        "Geometric (point) merging instead of topological merging "
-        "(slower, fails with high-aspect cells. default for 1912 and earlier)",
-        true  // mark as an advanced option
+        "Geometric point merging instead of topological merging"
+        " [default for 1912 and earlier]."
+        // NOTE: " Slower, fails with high-aspect cells."
     );
     argList::addBoolOption
     (
@@ -164,6 +164,9 @@ int main(int argc, char *argv[])
     // Remove old files, unless disabled
     const bool removeOldFiles = !args.found("noClean");
 
+    // Write cellSets
+    const bool writeCellSets = args.found("sets");
+
     // Default merge (topology), unless otherwise specified
     blockMesh::mergeStrategy strategy(blockMesh::DEFAULT_MERGE);
 
@@ -329,56 +332,22 @@ int main(int argc, char *argv[])
     }
 
 
-    Info<< nl << "Creating polyMesh from blockMesh" << endl;
+    // Ensure we get information messages, even if turned off in dictionary
+    blocks.verbose(true);
 
-    polyMesh mesh
-    (
-        IOobject
+    autoPtr<polyMesh> meshPtr =
+        blocks.mesh
         (
-            regionName,
-            meshInstance,
-            runTime
-        ),
-        pointField(blocks.points()),  // Copy, could we re-use space?
-        blocks.cells(),
-        blocks.patches(),
-        blocks.patchNames(),
-        blocks.patchDicts(),
-        "defaultFaces",               // Default patch name
-        emptyPolyPatch::typeName      // Default patch type
-    );
+            IOobject(regionName, meshInstance, runTime)
+        );
 
+    polyMesh& mesh = *meshPtr;
 
-    // Handle merging of patch pairs. Dictionary entry "mergePatchPairs"
+    // Merge patch pairs (dictionary entry "mergePatchPairs")
     #include "mergePatchPairs.H"
 
-    // Set any cellZones
-    #include "addCellZones.H"
-
-
-    // Detect any cyclic patches and force re-ordering of the faces
-    {
-        bool hasCyclic = false;
-        for (const polyPatch& pp : mesh.boundaryMesh())
-        {
-            if (isA<cyclicPolyPatch>(pp))
-            {
-                hasCyclic = true;
-                break;
-            }
-        }
-
-        if (hasCyclic)
-        {
-            Info<< nl << "Detected cyclic patches; ordering boundary faces"
-                << endl;
-            const word oldInstance = mesh.instance();
-            polyTopoChange meshMod(mesh);
-            meshMod.changeMesh(mesh, false);
-            mesh.setInstance(oldInstance);
-        }
-    }
-
+    // Handle cyclic patches
+    #include "handleCyclicPatches.H"
 
     // Set the precision of the points data to 10
     IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
@@ -386,7 +355,7 @@ int main(int argc, char *argv[])
     Info<< nl << "Writing polyMesh with "
         << mesh.cellZones().size() << " cellZones";
 
-    if (args.found("sets") && !mesh.cellZones().empty())
+    if (writeCellSets && !mesh.cellZones().empty())
     {
         Info<< " (written as cellSets too)";
     }
@@ -400,7 +369,7 @@ int main(int argc, char *argv[])
             << exit(FatalError);
     }
 
-    if (args.found("sets"))
+    if (writeCellSets)
     {
         for (const cellZone& cz : mesh.cellZones())
         {
diff --git a/applications/utilities/mesh/generation/blockMesh/handleCyclicPatches.H b/applications/utilities/mesh/generation/blockMesh/handleCyclicPatches.H
new file mode 100644
index 0000000000..eb79d02de2
--- /dev/null
+++ b/applications/utilities/mesh/generation/blockMesh/handleCyclicPatches.H
@@ -0,0 +1,22 @@
+// Detect any cyclic patches and force re-ordering of the faces
+{
+    bool hasCyclic = false;
+    for (const polyPatch& pp : mesh.boundaryMesh())
+    {
+        if (isA<cyclicPolyPatch>(pp))
+        {
+            hasCyclic = true;
+            break;
+        }
+    }
+
+    if (hasCyclic)
+    {
+        Info<< nl << "Detected cyclic patches; ordering boundary faces" << endl;
+
+        const word oldInstance = mesh.instance();
+        polyTopoChange meshMod(mesh);
+        meshMod.changeMesh(mesh, false);
+        mesh.setInstance(oldInstance);
+    }
+}
diff --git a/applications/utilities/mesh/generation/blockMesh/mergePatchPairs.H b/applications/utilities/mesh/generation/blockMesh/mergePatchPairs.H
index 2c94bcc0fa..be551d01f6 100644
--- a/applications/utilities/mesh/generation/blockMesh/mergePatchPairs.H
+++ b/applications/utilities/mesh/generation/blockMesh/mergePatchPairs.H
@@ -1,6 +1,6 @@
 // Handle merging of patch pairs
 {
-    List<Pair<word>> mergePatchPairs;
+    wordPairList mergePatchPairs;
 
     // Read in a list of dictionaries for the merge patch pairs
     if
@@ -14,19 +14,19 @@
         // Create and add point and face zones and mesh modifiers
         List<pointZone*> pz(mergePatchPairs.size());
         List<faceZone*> fz(3*mergePatchPairs.size());
-        List<cellZone*> cz(0);
+        List<cellZone*> cz;
 
-        forAll(mergePatchPairs, pairI)
+        forAll(mergePatchPairs, pairi)
         {
             const word mergeName
             (
-                mergePatchPairs[pairI].first()
-              + mergePatchPairs[pairI].second()
-              + name(pairI)
+                mergePatchPairs[pairi].first()
+              + mergePatchPairs[pairi].second()
+              + name(pairi)
             );
 
             // An empty zone for cut points
-            pz[pairI] = new pointZone
+            pz[pairi] = new pointZone
             (
                 mergeName + "CutPointZone",
                 0,
@@ -34,11 +34,11 @@
             );
 
             // Master patch
-            const word masterPatchName(mergePatchPairs[pairI].first());
+            const word masterPatchName(mergePatchPairs[pairi].first());
             const polyPatch& masterPatch =
                 mesh.boundaryMesh()[masterPatchName];
 
-            fz[3*pairI] = new faceZone
+            fz[3*pairi] = new faceZone
             (
                 mergeName + "MasterZone",
                 identity(masterPatch.size(), masterPatch.start()),
@@ -48,11 +48,11 @@
             );
 
             // Slave patch
-            const word slavePatchName(mergePatchPairs[pairI].second());
+            const word slavePatchName(mergePatchPairs[pairi].second());
             const polyPatch& slavePatch =
                 mesh.boundaryMesh()[slavePatchName];
 
-            fz[3*pairI + 1] = new faceZone
+            fz[3*pairi + 1] = new faceZone
             (
                 mergeName + "SlaveZone",
                 identity(slavePatch.size(), slavePatch.start()),
@@ -62,7 +62,7 @@
             );
 
             // An empty zone for cut faces
-            fz[3*pairI + 2] = new faceZone
+            fz[3*pairi + 2] = new faceZone
             (
                 mergeName + "CutFaceZone",
                 2,
@@ -78,30 +78,30 @@
         attachPolyTopoChanger polyMeshAttacher(mesh);
         polyMeshAttacher.setSize(mergePatchPairs.size());
 
-        forAll(mergePatchPairs, pairI)
+        forAll(mergePatchPairs, pairi)
         {
             const word mergeName
             (
-                mergePatchPairs[pairI].first()
-              + mergePatchPairs[pairI].second()
-              + name(pairI)
+                mergePatchPairs[pairi].first()
+              + mergePatchPairs[pairi].second()
+              + name(pairi)
             );
 
             // Add the sliding interface mesh modifier
             polyMeshAttacher.set
             (
-                pairI,
+                pairi,
                 new slidingInterface
                 (
-                    "couple" + name(pairI),
-                    pairI,
+                    "couple" + name(pairi),
+                    pairi,
                     polyMeshAttacher,
                     mergeName + "MasterZone",
                     mergeName + "SlaveZone",
                     mergeName + "CutPointZone",
                     mergeName + "CutFaceZone",
-                    mergePatchPairs[pairI].first(),
-                    mergePatchPairs[pairI].second(),
+                    mergePatchPairs[pairi].first(),
+                    mergePatchPairs[pairi].second(),
                     slidingInterface::INTEGRAL, // always integral
                     false,
                     intersection::VISIBLE
diff --git a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.C b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.C
index 39d3b99b8c..29957a0ceb 100644
--- a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.C
+++ b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.C
@@ -134,7 +134,7 @@ Foam::point Foam::BSpline::position
 Foam::scalar Foam::BSpline::length() const
 {
     NotImplemented;
-    return 1.0;
+    return 1;
 }
 
 
diff --git a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.H b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.H
index e9da49bcdc..00f77df5b0 100644
--- a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.H
+++ b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,21 +78,12 @@ class BSpline
 :
     public polyLine
 {
-    // Private Member Functions
-
-        //- No copy construct
-        BSpline(const BSpline&) = delete;
-
-        //- No copy assignment
-        void operator=(const BSpline&) = delete;
-
-
 public:
 
     // Constructors
 
         //- Construct from components
-        BSpline
+        explicit BSpline
         (
             const pointField& knots,
             const bool notImplementedClosed = false
@@ -100,15 +92,16 @@ public:
 
     // Member Functions
 
-        //- Return the point position corresponding to the curve parameter
+        //- The point position corresponding to the global curve parameter
         //  0 <= lambda <= 1
         point position(const scalar lambda) const;
 
-        //- Return the point position corresponding to the local parameter
-        //  0 <= lambda <= 1 on the given segment
+        //- The point position corresponding to the local lambda (0-1)
+        //- on the given segment
         point position(const label segment, const scalar lambda) const;
 
-        //- Return the length of the curve
+        //- The length of the curve
+        //  \note NotImplemented
         scalar length() const;
 };
 
diff --git a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C
index 81c3ab24ee..af5e4d89ed 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 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,9 +27,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "BSplineEdge.H"
+#include "polyLine.H"
 #include "addToRunTimeSelectionTable.H"
 
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
@@ -59,7 +59,10 @@ Foam::blockEdges::BSplineEdge::BSplineEdge
 )
 :
     blockEdge(points, start, end),
-    BSpline(appendEndPoints(points, start, end, internalPoints))
+    BSpline
+    (
+        polyLine::concat(points[start_], internalPoints, points[end_])
+    )
 {}
 
 
@@ -73,13 +76,16 @@ Foam::blockEdges::BSplineEdge::BSplineEdge
 )
 :
     blockEdge(dict, index, points, is),
-    BSpline(appendEndPoints(points, start_, end_, pointField(is)))
+    BSpline
+    (
+        polyLine::concat(points[start_], pointField(is), points[end_])
+    )
 {
-    token t(is);
-    is.putBack(t);
+    token tok(is);
+    is.putBack(tok);
 
-    // discard unused start/end tangents
-    if (t == token::BEGIN_LIST)
+    // Discard unused start/end tangents
+    if (tok == token::BEGIN_LIST)
     {
         vector tangent0Ignored(is);
         vector tangent1Ignored(is);
diff --git a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H
index e04434c593..68eefdab87 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 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,20 +77,20 @@ public:
         //- Construct from components
         BSplineEdge
         (
-            const pointField&,
-            const label start,
-            const label end,
+            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& internalPoints
         );
 
-        //- Construct from Istream, setting pointsList
+        //- Construct from Istream and point field.
         BSplineEdge
         (
             const dictionary& dict,
             const label index,
             const searchableSurfaces& geometry,
-            const pointField&,
-            Istream&
+            const pointField& points,   //!< Referenced point field
+            Istream& is
         );
 
 
@@ -100,11 +100,12 @@ public:
 
     // Member Functions
 
-        //- Return the point position corresponding to the curve parameter
+        //- The point position corresponding to the curve parameter
         //  0 <= lambda <= 1
         virtual point position(const scalar) const;
 
-        //- Return the length of the spline curve (not implemented)
+        //- The length of the spline curve
+        //  \note NotImplemented
         virtual scalar length() const;
 };
 
diff --git a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C
index 7c614b8829..698fec2a8d 100644
--- a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C
+++ b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,10 +44,15 @@ namespace blockEdges
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::coordSystem::cylindrical Foam::blockEdges::arcEdge::calcAngle()
+void Foam::blockEdges::arcEdge::calcFromMidPoint
+(
+    const point& p1,
+    const point& p3,
+    const point& p2
+)
 {
-    const vector a = p2_ - p1_;
-    const vector b = p3_ - p1_;
+    const vector a = p2 - p1;
+    const vector b = p3 - p1;
 
     // Find centre of arcEdge
     const scalar asqr = a & a;
@@ -55,7 +61,7 @@ Foam::coordSystem::cylindrical Foam::blockEdges::arcEdge::calcAngle()
 
     const scalar denom = asqr*bsqr - adotb*adotb;
 
-    if (mag(denom) < VSMALL)
+    if (mag(denom) < ROOTVSMALL)
     {
         FatalErrorInFunction
             << denom
@@ -64,44 +70,120 @@ Foam::coordSystem::cylindrical Foam::blockEdges::arcEdge::calcAngle()
 
     const scalar fact = 0.5*(bsqr - adotb)/denom;
 
-    point centre = 0.5*a + fact*((a ^ b) ^ a);
+    const point centre = p1 + 0.5*a + fact*((a ^ b) ^ a);
 
-    centre += p1_;
+    // Position vectors from centre
+    const vector r1(p1 - centre);
+    const vector r2(p2 - centre);
+    const vector r3(p3 - centre);
 
-    // Find position vectors w.r.t. the arcEdge centre
-    const vector r1(p1_ - centre);
-    const vector r2(p2_ - centre);
-    const vector r3(p3_ - centre);
+    const scalar mag1(mag(r1));
+    const scalar mag3(mag(r3));
 
-    // Find angle (in degrees)
-    angle_ = radToDeg(acos((r3 & r1)/(mag(r3) * mag(r1))));
+    vector arcAxis(r1 ^ r3);
+
+    // The radius from r1 and from r3 will be identical
+    radius_ = mag(r3);
+
+
+    // Determine the angle
+    angle_ = acos((r1 & r3)/(mag1*mag3));
 
     // Check if the vectors define an exterior or an interior arcEdge
     if (((r1 ^ r2) & (r1 ^ r3)) < 0.0)
     {
-        angle_ = 360.0 - angle_;
+        angle_ = constant::mathematical::twoPi - angle_;
     }
 
-    vector arcAxis;
-
-    if (angle_ <= 180.0)
+    if (angle_ <= constant::mathematical::pi)
     {
-        arcAxis = r1 ^ r3;
-
-        if (mag(arcAxis)/(mag(r1)*mag(r3)) < 0.001)
+        if (mag(arcAxis)/(mag1*mag3) < 0.001)
         {
             arcAxis = r1 ^ r2;
         }
     }
     else
     {
-        arcAxis = r3 ^ r1;
+        arcAxis = -arcAxis;
     }
 
-    radius_ = mag(r3);
+    // Corresponding local cylindrical coordinate system
+    cs_ = coordSystem::cylindrical(centre, arcAxis, r1);
+}
 
-    // The corresponding local cylindrical coordinate system (radians)
-    return coordSystem::cylindrical("arc", centre, arcAxis, r1);
+
+void Foam::blockEdges::arcEdge::calcFromCentre
+(
+    const point& p1,
+    const point& p3,
+    const point& centre,
+    bool adjustCentre,
+    scalar rMultiplier
+)
+{
+    // Position vectors from centre
+    const vector r1(p1 - centre);
+    const vector r3(p3 - centre);
+
+    const scalar mag1(mag(r1));
+    const scalar mag3(mag(r3));
+
+    const vector chord(p3 - p1);
+
+    const vector arcAxis(r1 ^ r3);
+
+    // The average radius
+    radius_ = 0.5*(mag1 + mag3);
+
+    // The included angle
+    angle_ = acos((r1 & r3)/(mag1*mag3));
+
+    // TODO? check for 180 degrees (co-linear points)?
+
+    bool needsAdjust = false;
+
+    if (adjustCentre)
+    {
+        needsAdjust = !equal(mag1, mag3);
+
+        if (!equal(rMultiplier, 1))
+        {
+            // The min radius is constrained by the chord,
+            // otherwise bad things will happen.
+
+            needsAdjust = true;
+            radius_ *= rMultiplier;
+            radius_ = max(radius_, (1.001*0.5*mag(chord)));
+        }
+    }
+
+    if (needsAdjust)
+    {
+        // The centre is not equidistant to p1 and p3.
+        // Use the chord and the arcAxis to determine the vector to
+        // the midpoint of the chord and adjust the centre along this
+        // line.
+
+        const point newCentre =
+        (
+            (0.5 * (p3 + p1))                   // mid-chord point
+          + sqrt(sqr(radius_) - 0.25 * magSqr(chord))
+          * normalised(arcAxis ^ chord)         // mid-chord -> centre
+        );
+
+        //// Info<< nl << "Adjust centre. r1=" << mag1 << " r3=" << mag3
+        ////     << " radius=" << radius_ << nl
+        ////     << "angle=" << radToDeg(angle_) << ' '
+        ////     << coordSystem::cylindrical(centre, arcAxis, r1) << nl;
+
+        // Recalculate - do attempt to readjust
+        calcFromCentre(p1, p3, newCentre, false);
+    }
+    else
+    {
+        // Corresponding local cylindrical coordinate system
+        cs_ = coordSystem::cylindrical(centre, arcAxis, r1);
+    }
 }
 
 
@@ -110,17 +192,35 @@ Foam::coordSystem::cylindrical Foam::blockEdges::arcEdge::calcAngle()
 Foam::blockEdges::arcEdge::arcEdge
 (
     const pointField& points,
+    const point& origin,
     const label start,
-    const label end,
-    const point& pMid
+    const label end
 )
 :
     blockEdge(points, start, end),
-    p1_(points_[start_]),
-    p2_(pMid),
-    p3_(points_[end_]),
-    cs_(calcAngle())
-{}
+    radius_(0),
+    angle_(0),
+    cs_()
+{
+    calcFromCentre(points[start_], points[end_], origin);
+}
+
+
+Foam::blockEdges::arcEdge::arcEdge
+(
+    const pointField& points,
+    const label start,
+    const label end,
+    const point& midPoint
+)
+:
+    blockEdge(points, start, end),
+    radius_(0),
+    angle_(0),
+    cs_()
+{
+    calcFromMidPoint(points[start_], points[end_], midPoint);
+}
 
 
 Foam::blockEdges::arcEdge::arcEdge
@@ -133,41 +233,82 @@ Foam::blockEdges::arcEdge::arcEdge
 )
 :
     blockEdge(dict, index, points, is),
-    p1_(points_[start_]),
-    p2_(is),
-    p3_(points_[end_]),
-    cs_(calcAngle())
-{}
+    radius_(0),
+    angle_(0),
+    cs_()
+{
+    point p;
+
+    token tok(is);
+    if (tok.isWord())
+    {
+        // Can be
+        //   - origin (0 0 0)
+        //   - origin 1.2 (0 0 0)
+
+        scalar rMultiplier = 1;
+
+        is >> tok;
+        if (tok.isNumber())
+        {
+            rMultiplier = tok.number();
+        }
+        else
+        {
+            is.putBack(tok);
+        }
+
+        is >> p;  // The origin (centre)
+
+        calcFromCentre(points_[start_], points_[end_], p, true, rMultiplier);
+    }
+    else
+    {
+        is.putBack(tok);
+
+        is >> p;  // A mid-point
+
+        calcFromMidPoint(points_[start_], points_[end_], p);
+    }
+
+    // Debug information
+    #if 0
+    Info<< "arc " << start_ << ' ' << end_
+        << ' ' << position(0.5) << ' ' << cs_
+        // << " radius=" << radius_ << " angle=" << radToDeg(angle_)
+        << nl;
+    #endif
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::point Foam::blockEdges::arcEdge::position(const scalar lambda) const
 {
+    #ifdef FULLDEBUG
     if (lambda < -SMALL || lambda > 1 + SMALL)
     {
-        FatalErrorInFunction
-            << "Parameter out of range, lambda = " << lambda
-            << abort(FatalError);
+        WarningInFunction
+            << "Parameter out of range, lambda = " << lambda << nl;
     }
+    #endif
 
     if (lambda < SMALL)
     {
-        return p1_;
+        return points_[start_];
     }
-    else if (lambda > 1 - SMALL)
+    else if (lambda >= 1 - SMALL)
     {
-        return p3_;
+        return points_[end_];
     }
 
-    // The angle is degrees, the coordinate system in radians
-    return cs_.globalPosition(vector(radius_, degToRad(lambda*angle_), 0));
+    return cs_.globalPosition(vector(radius_, (lambda*angle_), 0));
 }
 
 
-Foam::scalar Foam::blockEdges::arcEdge::length() const
+Foam::scalar Foam::blockEdges::arcEdge::length() const noexcept
 {
-    return degToRad(radius_*angle_);
+    return (radius_*angle_);
 }
 
 
diff --git a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H
index 8db6fcaa97..459b6a70f0 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 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,32 @@ Class
     Foam::blockEdges::arcEdge
 
 Description
-    Defines the arcEdge of a circle in terms of 3 points on its circumference
+    A blockEdge defined as an arc of a circle.
+
+    The arc is normally defined by its endpoints and a point on
+    its circumference, typically a midpoint. For example,
+    \verbatim
+    points ((1 0 0) (0 1 0));
+
+    arc 0 1 (0.707107 0.707107 0);
+    \endverbatim
+    The arc can enclose an angle greater than 0 and less than 360 degrees.
+
+    The arc will frequently enclose an angle less than 180 degrees.
+    For the case, it is possible to define the arc by its endpoints and its
+    centre (origin) point. For example,
+    \verbatim
+    arc 0 1 origin (0 0 0);
+    \endverbatim
+    When defined in the way, any discrepancy in the arc radius for the
+    endpoints is resolved by adjusting the origin to ensure that the average
+    radius is satisfied.
+
+    It is also possible to define a \em flatness factor as a multiplier of
+    the calculated radius. For example,
+    \verbatim
+    arc 0 1 origin 1.1 (0 0 0);
+    \endverbatim
 
 SourceFiles
     arcEdge.C
@@ -56,26 +81,41 @@ class arcEdge
 :
     public blockEdge
 {
-    // Private data
-
-        // Begin, mid, end points
-        point p1_, p2_, p3_;
-
-        //- The arc angle (in degrees)
-        scalar angle_;
+    // Private Data
 
         //- The arc radius
         scalar radius_;
 
+        //- The arc angle (radians)
+        scalar angle_;
+
         //- The local cylindrical coordinate system
         coordSystem::cylindrical cs_;
 
 
     // Private Member Functions
 
-        //- Calculate the angle, radius and axis
-        //  \return the cylindrical coordinate system
-        coordSystem::cylindrical calcAngle();
+        //- Calculate angle, radius, cylindrical coordinate system
+        //- from end points and the given point on the circumference
+        void calcFromMidPoint
+        (
+            const point& p1,        //!< Start point
+            const point& p3,        //!< End point
+            const point& p2         //!< Point on circumference
+        );
+
+        //- Calculate angle, radius, cylindrical coordinate system
+        //- from end points and the given origin.
+        //  Optionally adjust centre to accommodate deviations in the
+        //  effective radius to the end-points
+        void calcFromCentre
+        (
+            const point& p1,        //!< Start point
+            const point& p3,        //!< End point
+            const point& centre,    //!< Centre
+            bool adjustCentre = false,  //!< Adjust centre
+            scalar rMultiplier = 1      //!< Adjust radius by this factor
+        );
 
         //- No copy construct
         arcEdge(const arcEdge&) = delete;
@@ -92,23 +132,34 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from components, given the origin of the circle
         arcEdge
         (
-            const pointField& points,
-            const label start,
-            const label end,
-            const point& pMid
+            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
         );
 
-        //- Construct from Istream setting pointsList
+        //- 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
+        );
+
+        //- Construct from Istream and point field.
+        //  The Istream can specify either a point on the circumference,
+        //  or with a tag to specify the origin.
         arcEdge
         (
             const dictionary& dict,
             const label index,
-            const searchableSurfaces& geometry,
-            const pointField& points,
-            Istream&
+            const searchableSurfaces& geometry, // unsed
+            const pointField& points,   //!< Referenced point field
+            Istream& is
         );
 
 
@@ -122,7 +173,7 @@ public:
         point position(const scalar lambda) const;
 
         //- The length of the curve
-        scalar length() const;
+        scalar length() const noexcept;
 };
 
 
diff --git a/src/mesh/blockMesh/blockEdges/bezier/bezier.C b/src/mesh/blockMesh/blockEdges/bezier/bezier.C
index b1b884d32a..cac55d1ba3 100644
--- a/src/mesh/blockMesh/blockEdges/bezier/bezier.C
+++ b/src/mesh/blockMesh/blockEdges/bezier/bezier.C
@@ -26,6 +26,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "bezier.H"
+#include "polyLine.H"
 #include "SubList.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -66,7 +67,10 @@ Foam::blockEdges::bezier::bezier
 )
 :
     blockEdge(dict, index, points, is),
-    control_(appendEndPoints(points, start_, end_, pointField(is)))
+    control_
+    (
+        polyLine::concat(points[start_], pointField(is), points[end_])
+    )
 {}
 
 
@@ -94,7 +98,7 @@ Foam::point Foam::blockEdges::bezier::position(const scalar lambda) const
 Foam::scalar Foam::blockEdges::bezier::length() const
 {
     NotImplemented;
-    return 1.0;
+    return 1;
 }
 
 
diff --git a/src/mesh/blockMesh/blockEdges/bezier/bezier.H b/src/mesh/blockMesh/blockEdges/bezier/bezier.H
index b1f14dab27..5feaeb4139 100644
--- a/src/mesh/blockMesh/blockEdges/bezier/bezier.H
+++ b/src/mesh/blockMesh/blockEdges/bezier/bezier.H
@@ -63,7 +63,7 @@ class bezier
 :
     public blockEdge
 {
-private:
+    // Private Data
 
         //- Control points
         pointField control_;
@@ -89,20 +89,20 @@ public:
         //- Construct from components
         bezier
         (
-            const pointField& points,
-            const label start,
-            const label end,
-            const pointField& control
+            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& control   //!< The control points
         );
 
-        //- Construct from Istream
+        //- Construct from Istream and point field.
         bezier
         (
             const dictionary& dict,
             const label index,
             const searchableSurfaces& geometry,
-            const pointField&,
-            Istream&
+            const pointField& points,   //!< Referenced point field
+            Istream& is
         );
 
 
@@ -117,6 +117,7 @@ public:
         point position(const scalar lambda) const;
 
         //- Return the length of the curve
+        //  \not NotImplemented
         scalar length() const;
 };
 
diff --git a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C
index 01621e08c8..a49880b121 100644
--- a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C
+++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +28,7 @@ License
 
 #include "blockEdge.H"
 #include "blockVertex.H"
+#include "polyLine.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -104,32 +105,22 @@ Foam::autoPtr<Foam::blockEdge> Foam::blockEdge::New
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
 Foam::pointField Foam::blockEdge::appendEndPoints
 (
-    const pointField& points,
+    const pointField& pts,
     const label start,
     const label end,
-    const pointField& otherKnots
+    const pointField& intermediate
 )
 {
-    pointField allKnots(otherKnots.size() + 2);
-
-    // Start/end knots
-    allKnots[0] = points[start];
-    allKnots[otherKnots.size() + 1] = points[end];
-
-    // Intermediate knots
-    forAll(otherKnots, knotI)
-    {
-        allKnots[knotI+1] = otherKnots[knotI];
-    }
-
-    return allKnots;
+    return pointField(polyLine::concat(pts[start], intermediate, pts[end]));
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 Foam::tmp<Foam::pointField>
 Foam::blockEdge::position(const scalarList& lambdas) const
 {
diff --git a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H
index 30d09f00f1..dd48cba834 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 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,6 +24,12 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
+Namespace
+    Foam::blockEdges
+
+Description
+    A namespace for various blockEdge types.
+
 Class
     Foam::blockEdge
 
@@ -46,13 +52,10 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
-
+// Forward Declarations
 class blockEdge;
-
 Ostream& operator<<(Ostream&, const blockEdge&);
 
-
 /*---------------------------------------------------------------------------*\
                          Class blockEdge Declaration
 \*---------------------------------------------------------------------------*/
@@ -61,24 +64,29 @@ class blockEdge
 {
 protected:
 
-    // Protected data
+    // Protected Data
 
+        //- The referenced point field
         const pointField& points_;
 
+        //- Index of the start point
         const label start_;
+
+        //- Index of the end point
         const label end_;
 
 
     // Protected Member Functions
 
         //- Return a complete point field by appending the start/end points
-        //  to the given list
+        //- to the given list
+        //  \deprecated(2020-10) use polyLine::concat
         static pointField appendEndPoints
         (
-            const pointField&,
-            const label start,
-            const label end,
-            const pointField& otherKnots
+            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& intermediate  //!< Intermediate points (knots)
         );
 
 
@@ -110,18 +118,18 @@ public:
         //- Construct from components
         blockEdge
         (
-            const pointField& points,
-            const label start,
-            const label end
+            const pointField& points,   //!< Referenced point field
+            const label start,      //!< Start point in referenced point field
+            const label end         //!< End point in referenced point field
         );
 
-        //- Construct from Istream setting pointsList
+        //- Construct from Istream and point field.
         blockEdge
         (
             const dictionary& dict,
             const label index,
-            const pointField&,
-            Istream&
+            const pointField& points,   //!< Referenced point field
+            Istream& is
         );
 
         //- Clone function
@@ -133,8 +141,8 @@ public:
             const dictionary& dict,
             const label index,
             const searchableSurfaces& geometry,
-            const pointField&,
-            Istream&
+            const pointField& points,
+            Istream& is
         );
 
         //- Class used for the read-construction of
@@ -174,10 +182,10 @@ public:
 
     // Member Functions
 
-        //- Return label of start point
+        //- Index of start point
         inline label start() const;
 
-        //- Return label of end point
+        //- Index of end point
         inline label end() const;
 
         //- Compare the given start and end points with this curve
@@ -201,22 +209,22 @@ public:
         //  - -1: same edge, but different orientation
         inline int compare(const label start, const label end) const;
 
-        //- Return the point position corresponding to the curve parameter
+        //- The point position corresponding to the curve parameter
         //  0 <= lambda <= 1
         virtual point position(const scalar) const = 0;
 
-        //- Return the point positions corresponding to the curve parameters
+        //- The point positions corresponding to the curve parameters
         //  0 <= lambda <= 1
         virtual tmp<pointField> position(const scalarList&) const;
 
-        //- Return the length of the curve
+        //- The length of the curve
         virtual scalar length() const = 0;
 
         //- Write edge with variable backsubstitution
         void write(Ostream&, const dictionary&) const;
 
 
-    // Ostream operator
+    // Ostream Operator
 
         friend Ostream& operator<<(Ostream&, const blockEdge&);
 };
diff --git a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdgeI.H b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdgeI.H
index a781be5f96..6aab3064f0 100644
--- a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdgeI.H
+++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdgeI.H
@@ -49,10 +49,8 @@ inline int Foam::blockEdge::compare(const label start, const label end) const
     {
         return -1;
     }
-    else
-    {
-        return 0;
-    }
+
+    return 0;
 }
 
 
diff --git a/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.H b/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.H
index acd035602d..1c58386a0f 100644
--- a/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.H
+++ b/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,6 +47,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class blockEdge;
 
 /*---------------------------------------------------------------------------*\
@@ -54,7 +56,7 @@ class blockEdge;
 
 class lineDivide
 {
-    // Private data
+    // Private Data
 
         pointField points_;
 
@@ -67,18 +69,18 @@ public:
         //- Construct from components
         lineDivide
         (
-            const blockEdge&,
-            const label ndiv,
+            const blockEdge& cedge,
+            const label nDiv,
             const gradingDescriptors& gd = gradingDescriptors()
         );
 
 
     // Member Functions
 
-        //- Return the points
+        //- The points
         const pointField& points() const;
 
-        //- Return the list of lambda values
+        //- The list of lambda values
         const scalarList& lambdaDivisions() const;
 };
 
diff --git a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C
index 713d9bec15..203f16220a 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 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -71,11 +71,21 @@ Foam::blockEdges::lineEdge::lineEdge
 
 Foam::point Foam::blockEdges::lineEdge::position(const scalar lambda) const
 {
-    if (lambda < -SMALL || lambda > 1+SMALL)
+    #ifdef FULLDEBUG
+    if (lambda < -SMALL || lambda > 1 + SMALL)
     {
-        FatalErrorInFunction
-            << "Parameter out of range, lambda = " << lambda
-            << abort(FatalError);
+        WarningInFunction
+            << "Parameter out of range, lambda = " << 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_]);
diff --git a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H
index b44b3ab14a..1bbc9d67ce 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 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,12 +51,10 @@ namespace blockEdges
                           Class lineEdge Declaration
 \*---------------------------------------------------------------------------*/
 
-
 class lineEdge
 :
     public blockEdge
 {
-
 public:
 
     //- Runtime type information
@@ -66,15 +64,20 @@ public:
     // Constructors
 
         //- Construct from components
-        lineEdge(const pointField&, const label start, const label end);
+        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
+        );
 
-        //- Construct from Istream with a pointField
+        //- Construct from Istream and point field.
         lineEdge
         (
             const dictionary& dict,
             const label index,
             const searchableSurfaces& geometry,
-            const pointField&,
+            const pointField& points,   //!< Referenced point field
             Istream& is
         );
 
@@ -85,11 +88,11 @@ public:
 
     // Member Functions
 
-        //- Return the point position corresponding to the curve parameter
+        //- The point position corresponding to the curve parameter
         //  0 <= lambda <= 1
         point position(const scalar) const;
 
-        //- Return the length of the curve
+        //- The length of the curve
         scalar length() const;
 };
 
diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.C b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.C
index 0703696698..42f6a74db6 100644
--- a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.C
+++ b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,16 +27,41 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "polyLine.H"
+#include "SubList.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::tmp<Foam::pointField> Foam::polyLine::concat
+(
+    const point& p0,
+    const pointField& intermediate,
+    const point& p1
+)
+{
+    auto tresult = tmp<pointField>::New(intermediate.size() + 2);
+    auto& result = tresult.ref();
+
+    // Intermediate points (knots)
+    SubList<point>(result, intermediate.size(), 1) = intermediate;
+
+    // Start/end points (knots)
+    result.first() = p0;
+    result.last() = p1;
+
+    return tresult;
+}
+
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 void Foam::polyLine::calcParam()
 {
-    param_.setSize(points_.size());
+    lineLength_ = 0;
+    param_.resize(points_.size());
 
     if (param_.size())
     {
-        param_[0] = 0.0;
+        param_[0] = 0;
 
         for (label i=1; i < param_.size(); i++)
         {
@@ -48,23 +74,34 @@ void Foam::polyLine::calcParam()
         {
             param_[i] /= lineLength_;
         }
-        param_.last() = 1.0;
-    }
-    else
-    {
-        lineLength_ = 0.0;
+        param_.last() = 1;
     }
 }
 
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::polyLine::polyLine(const pointField& ps, const bool)
 :
     points_(ps),
-    lineLength_(0.0),
-    param_(0)
+    lineLength_(0),
+    param_()
+{
+    calcParam();
+}
+
+
+Foam::polyLine::polyLine
+(
+    const point& start,
+    const pointField& intermediate,
+    const point& end,
+    const bool
+)
+:
+    points_(polyLine::concat(start, intermediate, end)),
+    lineLength_(0),
+    param_()
 {
     calcParam();
 }
@@ -101,19 +138,20 @@ Foam::label Foam::polyLine::localParameter(scalar& lambda) const
     // Search table of cumulative distances to find which line-segment
     // we are on.
     // Check the upper bound.
+    // Too small to bother with a binary search...
 
-    label segmentI = 1;
-    while (param_[segmentI] < lambda)
+    label segment = 1;
+    while (param_[segment] < lambda)
     {
-        segmentI++;
+        ++segment;
     }
-    segmentI--;   // We want the corresponding lower bound
+    --segment;   // We want the corresponding lower bound
 
     // The local parameter [0-1] on this line segment
     lambda =
-        (lambda - param_[segmentI])/(param_[segmentI+1] - param_[segmentI]);
+        (lambda - param_[segment])/(param_[segment+1] - param_[segment]);
 
-    return segmentI;
+    return segment;
 }
 
 
@@ -151,8 +189,8 @@ Foam::point Foam::polyLine::position
         return points_.last();
     }
 
-    const point& p0 = points()[segment];
-    const point& p1 = points()[segment+1];
+    const point& p0 = points_[segment];
+    const point& p1 = points_[segment+1];
 
     // Special cases - no calculation needed
     if (mu <= 0.0)
@@ -163,11 +201,9 @@ Foam::point Foam::polyLine::position
     {
         return p1;
     }
-    else
-    {
-        // Linear interpolation
-        return points_[segment] + mu*(p1 - p0);
-    }
+
+    // Linear interpolation
+    return points_[segment] + mu*(p1 - p0);
 }
 
 
diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.H b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.H
index e720812201..b7bf40e8af 100644
--- a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.H
+++ b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,26 +53,16 @@ namespace Foam
                           Class polyLine Declaration
 \*---------------------------------------------------------------------------*/
 
-
 class polyLine
 {
-    // Private Member Functions
-
-        //- No copy construct
-        polyLine(const polyLine&) = delete;
-
-        //- No copy assignment
-        void operator=(const polyLine&) = delete;
-
-
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- The control points or ends of each segments
         pointField points_;
 
-        //- The real line length
+        //- The real (total) line length
         scalar lineLength_;
 
         //- The rational (0-1) cumulative parameter value for each point
@@ -81,11 +72,11 @@ protected:
     // Protected Member Functions
 
         //- Precalculate the rational cumulative parameter value
-        //  and the line-length
+        //- and the line-length
         void calcParam();
 
         //- Return the line segment and the local parameter [0..1]
-        //  corresponding to the global lambda [0..1]
+        //- corresponding to the global lambda [0..1]
         label localParameter(scalar& lambda) const;
 
 
@@ -94,30 +85,50 @@ public:
     // Constructors
 
         //- Construct from components
+        explicit polyLine
+        (
+            const pointField& points,   //!< The poly-line points
+            const bool notImplementedClosed = false
+        );
+
+        //- Construct from begin, intermediate, end points
         polyLine
         (
-            const pointField&,
+            const point& start,             //!< The start point
+            const pointField& intermediate, //!< The intermediate points
+            const point& end,               //!< The end point
             const bool notImplementedClosed = false
         );
 
 
+    // Static Member Functions
+
+        //- Concatenate begin, intermediate and end points
+        static tmp<pointField> concat
+        (
+            const point& start,             //!< The start point
+            const pointField& intermediate, //!< The intermediate points
+            const point& end                //!< The end point
+        );
+
+
     // Member Functions
 
         //- Return const-access to the control-points
         const pointField& points() const;
 
-        //- Return the number of line segments
+        //- The number of line segments
         label nSegments() const;
 
-        //- Return the point position corresponding to the curve parameter
+        //- The point position corresponding to the curve parameter
         //  0 <= lambda <= 1
         point position(const scalar) const;
 
-        //- Return the point position corresponding to the local parameter
+        //- The point position corresponding to the local parameter
         //  0 <= lambda <= 1 on the given segment
         point position(const label segment, const scalar) const;
 
-        //- Return the length of the curve
+        //- The length of the curve
         scalar length() const;
 };
 
diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C
index c350a82322..9a1b014fee 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 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,11 +48,11 @@ Foam::blockEdges::polyLineEdge::polyLineEdge
     const pointField& ps,
     const label start,
     const label end,
-    const pointField& otherPoints
+    const pointField& intermediate
 )
 :
     blockEdge(ps, start, end),
-    polyLine(appendEndPoints(ps, start_, end_, otherPoints))
+    polyLine(ps[start_], intermediate, ps[end_])
 {}
 
 
@@ -66,7 +66,7 @@ Foam::blockEdges::polyLineEdge::polyLineEdge
 )
 :
     blockEdge(dict, index, ps, is),
-    polyLine(appendEndPoints(ps, start_, end_, pointField(is)))
+    polyLine(ps[start_], pointField(is), ps[end_])
 {}
 
 
diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H
index 9eb6690143..bcce9c738e 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 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,19 +77,19 @@ public:
         //- Construct from components
         polyLineEdge
         (
-            const pointField&,
-            const label start,
-            const label end,
-            const pointField& otherPoints
+            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& intermediate  //!< The intermediate points
         );
 
-        //- Construct from Istream
+        //- Construct from Istream and point field.
         polyLineEdge
         (
             const dictionary& dict,
             const label index,
             const searchableSurfaces& geometry,
-            const pointField&,
+            const pointField& points,       //!< Referenced point field
             Istream& is
         );
 
@@ -100,11 +100,11 @@ public:
 
     // Member Functions
 
-        //- Return the point position corresponding to the curve parameter
+        //- The point position corresponding to the curve parameter
         //  0 <= lambda <= 1
         point position(const scalar lambda) const;
 
-        //- Return the length of the curve
+        //- The length of the curve
         scalar length() const;
 };
 
diff --git a/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C
index 08fa1186a5..f6f182207a 100644
--- a/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C
+++ b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C
@@ -59,7 +59,7 @@ Foam::projectCurveEdge::projectCurveEdge
     geometry_(geometry)
 {
     wordList names(is);
-    surfaces_.setSize(names.size());
+    surfaces_.resize(names.size());
     forAll(names, i)
     {
         surfaces_[i] = geometry_.findSurfaceID(names[i]);
@@ -83,6 +83,14 @@ Foam::projectCurveEdge::projectCurveEdge
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::point
+Foam::projectCurveEdge::position(const scalar) const
+{
+    NotImplemented;
+    return point::max;
+}
+
+
 Foam::tmp<Foam::pointField>
 Foam::projectCurveEdge::position(const scalarList& lambdas) const
 {
@@ -101,8 +109,8 @@ Foam::projectCurveEdge::position(const scalarList& lambdas) const
     }
 
 
-    tmp<pointField> tpoints(new pointField(lambdas.size()));
-    pointField& points = tpoints.ref();
+    auto tpoints = tmp<pointField>::New(lambdas.size());
+    auto& points = tpoints.ref();
 
     const point& startPt = points_[start_];
     const point& endPt = points_[end_];
@@ -149,10 +157,11 @@ Foam::projectCurveEdge::position(const scalarList& lambdas) const
 
 
     // Upper limit for number of iterations
-    const label maxIter = 10;
+    constexpr label maxIter = 10;
+
     // Residual tolerance
-    const scalar relTol = 0.1;
-    const scalar absTol = 1e-4;
+    constexpr scalar relTol = 0.1;
+    constexpr scalar absTol = 1e-4;
 
     scalar initialResidual = 0.0;
 
@@ -258,4 +267,11 @@ Foam::projectCurveEdge::position(const scalarList& lambdas) const
 }
 
 
+Foam::scalar Foam::projectCurveEdge::length() const
+{
+    NotImplemented;
+    return 1;
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.H b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.H
index 1e6ca665ee..7f939b9d58 100644
--- a/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.H
+++ b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,6 +46,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class pointConstraint;
 
 /*---------------------------------------------------------------------------*\
@@ -81,13 +82,13 @@ public:
 
     // Constructors
 
-        //- Construct from Istream setting pointsList
+        //- Construct from Istream and point field.
         projectCurveEdge
         (
             const dictionary& dict,
             const label index,
             const searchableSurfaces& geometry,
-            const pointField& points,
+            const pointField& points,   //!< Referenced point field
             Istream& is
         );
 
@@ -98,24 +99,18 @@ public:
 
     // Member Functions
 
-        //- Return the point positions corresponding to the curve parameters
+        //- The point position corresponding to the curve parameter
         //  0 <= lambda <= 1
-        virtual point position(const scalar) const
-        {
-            NotImplemented;
-            return point::max;
-        }
+        //  \note NotImplemented
+        virtual point position(const scalar) const;
 
-        //- Return the point positions corresponding to the curve parameters
+        //- The point positions corresponding to the curve parameters
         //  0 <= lambda <= 1
         virtual tmp<pointField> position(const scalarList&) const;
 
-        //- Return the length of the curve
-        virtual scalar length() const
-        {
-            NotImplemented;
-            return 1;
-        }
+        //- The length of the curve
+        //  \note NotImplemented
+        virtual scalar length() const;
 };
 
 
diff --git a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C
index 8259bc77ed..16b99b1bbf 100644
--- a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C
+++ b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C
@@ -26,13 +26,12 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "searchableSurfacesQueries.H"
 #include "projectEdge.H"
-#include "unitConversion.H"
-#include "addToRunTimeSelectionTable.H"
 #include "pointConstraint.H"
+#include "searchableSurfacesQueries.H"
 #include "OBJstream.H"
 #include "linearInterpolationWeights.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -146,8 +145,8 @@ Foam::projectEdge::position(const scalarList& lambdas) const
     }
 
 
-    tmp<pointField> tpoints(new pointField(lambdas.size()));
-    pointField& points = tpoints.ref();
+    auto tpoints = tmp<pointField>::New(lambdas.size());
+    auto& points = tpoints.ref();
 
     const point& startPt = points_[start_];
     const point& endPt = points_[end_];
@@ -161,10 +160,11 @@ Foam::projectEdge::position(const scalarList& lambdas) const
 
 
     // Upper limit for number of iterations
-    const label maxIter = 10;
+    constexpr label maxIter = 10;
+
     // Residual tolerance
-    const scalar relTol = 0.1;
-    const scalar absTol = 1e-4;
+    constexpr scalar relTol = 0.1;
+    constexpr scalar absTol = 1e-4;
 
     scalar initialResidual = 0.0;
 
@@ -270,4 +270,11 @@ Foam::projectEdge::position(const scalarList& lambdas) const
 }
 
 
+Foam::scalar Foam::projectEdge::length() const
+{
+    NotImplemented;
+    return 1;
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H
index 1c3d180d7f..44c9d143d8 100644
--- a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H
+++ b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,6 +46,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class pointConstraint;
 
 /*---------------------------------------------------------------------------*\
@@ -84,13 +85,13 @@ public:
 
     // Constructors
 
-        //- Construct from Istream setting pointsList
+        //- Construct from Istream and point field.
         projectEdge
         (
             const dictionary& dict,
             const label index,
             const searchableSurfaces& geometry,
-            const pointField& points,
+            const pointField& points,   //!< Referenced point field
             Istream& is
         );
 
@@ -101,20 +102,17 @@ public:
 
     // Member Functions
 
-        //- Return the point positions corresponding to the curve parameters
+        //- The point position corresponding to the curve parameter
         //  0 <= lambda <= 1
         virtual point position(const scalar) const;
 
-        //- Return the point positions corresponding to the curve parameters
+        //- The point positions corresponding to the curve parameters
         //  0 <= lambda <= 1
         virtual tmp<pointField> position(const scalarList&) const;
 
-        //- Return the length of the curve
-        virtual scalar length() const
-        {
-            NotImplemented;
-            return 1;
-        }
+        //- The length of the edge
+        //  \note NotImplemented
+        virtual scalar length() const;
 };
 
 
diff --git a/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.C b/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.C
index a8ee8c1d71..bfa0aaf27b 100644
--- a/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.C
+++ b/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.C
@@ -134,7 +134,7 @@ Foam::point Foam::CatmullRomSpline::position
 Foam::scalar Foam::CatmullRomSpline::length() const
 {
     NotImplemented;
-    return 1.0;
+    return 1;
 }
 
 
diff --git a/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.H b/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.H
index c78803d9d3..27093f9cbd 100644
--- a/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.H
+++ b/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -79,21 +80,12 @@ class CatmullRomSpline
 :
     public polyLine
 {
-    // Private Member Functions
-
-        //- No copy construct
-        CatmullRomSpline(const CatmullRomSpline&) = delete;
-
-        //- No copy assignment
-        void operator=(const CatmullRomSpline&) = delete;
-
-
 public:
 
     // Constructors
 
         //- Construct from components
-        CatmullRomSpline
+        explicit CatmullRomSpline
         (
             const pointField& knots,
             const bool notImplementedClosed = false
@@ -102,15 +94,16 @@ public:
 
     // Member Functions
 
-        //- Return the point position corresponding to the curve parameter
+        //- The point position corresponding to the curve parameter
         //  0 <= lambda <= 1
         point position(const scalar lambda) const;
 
-        //- Return the point position corresponding to the local parameter
+        //- The point position corresponding to the local parameter
         //  0 <= lambda <= 1 on the given segment
         point position(const label segment, const scalar lambda) const;
 
-        //- Return the length of the curve
+        //- The length of the curve
+        //  \note NotImplemented
         scalar length() const;
 };
 
diff --git a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C
index 048d19878b..5494d4ed88 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 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,9 +27,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "splineEdge.H"
+#include "polyLine.H"
 #include "addToRunTimeSelectionTable.H"
 
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
@@ -59,7 +59,10 @@ Foam::blockEdges::splineEdge::splineEdge
 )
 :
     blockEdge(points, start, end),
-    CatmullRomSpline(appendEndPoints(points, start, end, internalPoints))
+    CatmullRomSpline
+    (
+        polyLine::concat(points[start_], internalPoints, points[end_])
+    )
 {}
 
 
@@ -73,13 +76,16 @@ Foam::blockEdges::splineEdge::splineEdge
 )
 :
     blockEdge(dict, index, points, is),
-    CatmullRomSpline(appendEndPoints(points, start_, end_, pointField(is)))
+    CatmullRomSpline
+    (
+        polyLine::concat(points[start_], pointField(is), points[end_])
+    )
 {
-    token t(is);
-    is.putBack(t);
+    token tok(is);
+    is.putBack(tok);
 
-    // discard unused start/end tangents
-    if (t == token::BEGIN_LIST)
+    // Discard unused start/end tangents
+    if (tok == token::BEGIN_LIST)
     {
         vector tangent0Ignored(is);
         vector tangent1Ignored(is);
diff --git a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H
index a5519ec883..3ac595cbac 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 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,6 +30,10 @@ Class
 Description
     A blockEdge interface for Catmull-Rom splines.
 
+See also
+    BSpline
+    CatmullRomSpline
+
 SourceFiles
     splineEdge.C
 
@@ -77,19 +81,19 @@ public:
         //- Construct from components
         splineEdge
         (
-            const pointField&,
-            const label start,
-            const label end,
+            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& internalPoints
         );
 
-        //- Construct from Istream, setting pointsList
+        //- Construct from Istream and point field.
         splineEdge
         (
             const dictionary& dict,
             const label index,
             const searchableSurfaces& geometry,
-            const pointField&,
+            const pointField& points,   //!< Referenced point field
             Istream&
         );
 
@@ -100,11 +104,12 @@ public:
 
     // Member Functions
 
-        //- Return the point position corresponding to the curve parameter
+        //- The point position corresponding to the curve parameter
         //  0 <= lambda <= 1
         virtual point position(const scalar) const;
 
-        //- Return the length of the spline curve (not implemented)
+        //- The length of the spline curve
+        //  \note NotImplemented
         virtual scalar length() const;
 };
 
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.C b/src/mesh/blockMesh/blockMesh/blockMesh.C
index 3df87de8d2..04db49d80c 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.C
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.C
@@ -109,15 +109,23 @@ Foam::blockMesh::blockMesh
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::blockMesh::valid() const
+bool Foam::blockMesh::valid() const noexcept
 {
     return bool(topologyPtr_);
 }
 
 
-void Foam::blockMesh::verbose(const bool on)
+bool Foam::blockMesh::verbose() const noexcept
 {
+    return verboseOutput;
+}
+
+
+bool Foam::blockMesh::verbose(const bool on) noexcept
+{
+    bool old(verboseOutput);
     verboseOutput = on;
+    return old;
 }
 
 
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.H b/src/mesh/blockMesh/blockMesh/blockMesh.H
index 4227b3de06..d08a475585 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.H
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.H
@@ -226,10 +226,10 @@ public:
 
         //- Construct from IOdictionary for given region
         //  Default is topological merging.
-        blockMesh
+        explicit blockMesh
         (
             const IOdictionary& dict,
-            const word& regionName,
+            const word& regionName = polyMesh::defaultRegion,
             mergeStrategy strategy = mergeStrategy::DEFAULT_MERGE
         );
 
@@ -255,10 +255,10 @@ public:
             }
 
             //- True if the blockMesh topology exists
-            bool valid() const;
+            bool valid() const noexcept;
 
-            //- Reference to point field defining the blockMesh
-            //  these points have not been scaled by scaleFactor
+            //- Reference to point field defining the blockMesh.
+            //  These points are \b not scaled by scaleFactor
             const pointField& vertices() const;
 
             //- Return the blockMesh topology as a polyMesh
@@ -279,8 +279,8 @@ public:
             //- The scaling factor used to convert to metres
             scalar scaleFactor() const;
 
-            //- The points for the entire mesh
-            //  these points have been scaled by scaleFactor
+            //- The points for the entire mesh.
+            //  These points \b are scaled by scaleFactor
             const pointField& points() const;
 
             //- Return cell shapes list
@@ -299,16 +299,26 @@ public:
             label numZonedBlocks() const;
 
 
-        // Edit
+    // Verbosity
 
-            //- Enable/disable verbose information about the progress
-            void verbose(const bool on=true);
+        //- Verbose information?
+        bool verbose() const noexcept;
+
+        //- Enable/disable verbose information about the progress
+        //  \return old value
+        bool verbose(const bool on) noexcept;
 
 
-        // Write
+    // Mesh Generation
 
-            //- Writes edges of blockMesh in OBJ format.
-            void writeTopology(Ostream&) const;
+        //- Create polyMesh, with cell zones
+        autoPtr<polyMesh> mesh(const IOobject& io) const;
+
+
+    // Write
+
+        //- Writes edges of blockMesh in OBJ format.
+        void writeTopology(Ostream& os) const;
 };
 
 
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshCreate.C b/src/mesh/blockMesh/blockMesh/blockMeshCreate.C
index 4e8dad0cc6..8d93b69b87 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshCreate.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshCreate.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +28,7 @@ License
 
 #include "blockMesh.H"
 #include "cellModel.H"
+#include "emptyPolyPatch.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -270,4 +271,129 @@ void Foam::blockMesh::createPatches() const
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::polyMesh>
+Foam::blockMesh::mesh(const IOobject& io) const
+{
+    const blockMesh& blkMesh = *this;
+
+    if (verboseOutput)
+    {
+        Info<< nl << "Creating polyMesh from blockMesh" << endl;
+    }
+
+    auto meshPtr = autoPtr<polyMesh>::New
+    (
+        io,
+        pointField(blkMesh.points()),   // Copy, could we re-use space?
+        blkMesh.cells(),
+        blkMesh.patches(),
+        blkMesh.patchNames(),
+        blkMesh.patchDicts(),
+        "defaultFaces",                 // Default patch name
+        emptyPolyPatch::typeName        // Default patch type
+    );
+
+
+    // Set any cellZones
+    const label nZones = blkMesh.numZonedBlocks();
+
+    if (nZones)
+    {
+        polyMesh& pmesh = *meshPtr;
+
+        if (verboseOutput)
+        {
+            Info<< "Adding cell zones" << endl;
+        }
+
+        // Map from zoneName to cellZone index
+        HashTable<label> zoneMap(2*nZones);
+
+        // Cells per zone
+        List<DynamicList<label>> zoneCells(nZones);
+
+        // Running cell counter
+        label celli = 0;
+
+        // Largest zone so far
+        label freeZonei = 0;
+
+        for (const block& b : blkMesh)
+        {
+            const word& zoneName = b.zoneName();
+            const label nCellsInBlock = b.cells().size();
+
+            if (zoneName.size())
+            {
+                const auto iter = zoneMap.cfind(zoneName);
+
+                label zonei = freeZonei;
+
+                if (iter.found())
+                {
+                    zonei = *iter;
+                }
+                else
+                {
+                    zoneMap.insert(zoneName, zonei);
+                    ++freeZonei;
+
+                    if (verboseOutput)
+                    {
+                        Info<< "    " << zonei << '\t' << zoneName << endl;
+                    }
+                }
+
+
+                // Fill with cell ids
+
+                zoneCells[zonei].reserve
+                (
+                    zoneCells[zonei].size() + nCellsInBlock
+                );
+
+                const label endOfFill = celli + nCellsInBlock;
+
+                for (; celli < endOfFill; ++celli)
+                {
+                    zoneCells[zonei].append(celli);
+                }
+            }
+            else
+            {
+                celli += nCellsInBlock;
+            }
+        }
+
+        List<cellZone*> cz(zoneMap.size());
+        forAllConstIters(zoneMap, iter)
+        {
+            const word& zoneName = iter.key();
+            const label zonei = iter.val();
+
+            cz[zonei] = new cellZone
+            (
+                zoneName,
+                zoneCells[zonei].shrink(),
+                zonei,
+                pmesh.cellZones()
+            );
+        }
+
+        pmesh.pointZones().resize(0);
+        pmesh.faceZones().resize(0);
+        pmesh.cellZones().resize(0);
+        pmesh.addZones(List<pointZone*>(), List<faceZone*>(), cz);
+    }
+
+
+    // Merge patch pairs, cyclic must be done elsewhere
+    // - requires libdynamicMesh
+
+    return meshPtr;
+}
+
+
 // ************************************************************************* //
diff --git a/tutorials/basic/potentialFoam/cylinder/system/blockMeshDict b/tutorials/basic/potentialFoam/cylinder/system/blockMeshDict
index ac3995c8dc..ceb0865644 100644
--- a/tutorials/basic/potentialFoam/cylinder/system/blockMeshDict
+++ b/tutorials/basic/potentialFoam/cylinder/system/blockMeshDict
@@ -16,6 +16,10 @@ FoamFile
 
 scale   1;
 
+// Front/back locations
+zmin -0.5;
+zmax 0.5;
+
 vertices #codeStream
 {
     codeInclude
@@ -25,34 +29,39 @@ vertices #codeStream
 
     code
     #{
-        pointField points(19);
-        points[0]  = point(0.5, 0, -0.5);
-        points[1]  = point(1, 0, -0.5);
-        points[2]  = point(2, 0, -0.5);
-        points[3]  = point(2, 0.707107, -0.5);
-        points[4]  = point(0.707107, 0.707107, -0.5);
-        points[5]  = point(0.353553, 0.353553, -0.5);
-        points[6]  = point(2, 2, -0.5);
-        points[7]  = point(0.707107, 2, -0.5);
-        points[8]  = point(0, 2, -0.5);
-        points[9]  = point(0, 1, -0.5);
-        points[10] = point(0, 0.5, -0.5);
-        points[11] = point(-0.5, 0, -0.5);
-        points[12] = point(-1, 0, -0.5);
-        points[13] = point(-2, 0, -0.5);
-        points[14] = point(-2, 0.707107, -0.5);
-        points[15] = point(-0.707107, 0.707107, -0.5);
-        points[16] = point(-0.353553, 0.353553, -0.5);
-        points[17] = point(-2, 2, -0.5);
-        points[18] = point(-0.707107, 2, -0.5);
+        // sin(45), cos(45)
+        const scalar sqrt05 = sqrt(0.5);
 
-        // Duplicate z points
-        label sz = points.size();
-        points.setSize(2*sz);
-        for (label i = 0; i < sz; i++)
+        pointField points
+        ({
+            /* 0*/ {0.5, 0, $zmin},
+            /* 1*/ {1, 0, $zmin},
+            /* 2*/ {2, 0, $zmin},
+            /* 3*/ {2, sqrt05, $zmin},
+            /* 4*/ {sqrt05, sqrt05, $zmin},
+            /* 5*/ {sqrt05/2, sqrt05/2, $zmin},
+            /* 6*/ {2, 2, $zmin},
+            /* 7*/ {sqrt05, 2, $zmin},
+            /* 8*/ {0, 2, $zmin},
+            /* 9*/ {0, 1, $zmin},
+            /*10*/ {0, 0.5, $zmin},
+            /*11*/ {-0.5, 0, $zmin},
+            /*12*/ {-1, 0, $zmin},
+            /*13*/ {-2, 0, $zmin},
+            /*14*/ {-2, sqrt05, $zmin},
+            /*15*/ {-sqrt05, sqrt05, $zmin},
+            /*16*/ {-sqrt05/2, sqrt05/2, $zmin},
+            /*17*/ {-2, 2, $zmin},
+            /*18*/ {-sqrt05, 2, $zmin}
+        });
+
+        // Duplicate z points for zmax
+        const label sz = points.size();
+        points.resize(2*sz);
+        for (label i = 0; i < sz; ++i)
         {
             const point& pt = points[i];
-            points[i+sz] = point(pt.x(), pt.y(), -pt.z());
+            points[i + sz] = point(pt.x(), pt.y(), $zmax);
         }
 
         os  << points;
@@ -76,22 +85,24 @@ blocks
 
 edges
 (
-    arc 0 5 (0.469846 0.17101 -0.5)
-    arc 5 10 (0.17101 0.469846 -0.5)
-    arc 1 4 (0.939693 0.34202 -0.5)
-    arc 4 9 (0.34202 0.939693 -0.5)
-    arc 19 24 (0.469846 0.17101 0.5)
-    arc 24 29 (0.17101 0.469846 0.5)
-    arc 20 23 (0.939693 0.34202 0.5)
-    arc 23 28 (0.34202 0.939693 0.5)
-    arc 11 16 (-0.469846 0.17101 -0.5)
-    arc 16 10 (-0.17101 0.469846 -0.5)
-    arc 12 15 (-0.939693 0.34202 -0.5)
-    arc 15 9 (-0.34202 0.939693 -0.5)
-    arc 30 35 (-0.469846 0.17101 0.5)
-    arc 35 29 (-0.17101 0.469846 0.5)
-    arc 31 34 (-0.939693 0.34202 0.5)
-    arc 34 28 (-0.34202 0.939693 0.5)
+    // Inner cylinder
+    arc  0  5 origin (0 0 $zmin)
+    arc  5 10 origin (0 0 $zmin)
+    arc  1  4 origin (0 0 $zmin)
+    arc  4  9 origin (0 0 $zmin)
+    arc 19 24 origin (0 0 $zmax)
+    arc 24 29 origin (0 0 $zmax)
+    arc 20 23 origin (0 0 $zmax)
+    arc 23 28 origin (0 0 $zmax)
+    // Intermediate cylinder
+    arc 11 16 origin (0 0 $zmin)
+    arc 16 10 origin (0 0 $zmin)
+    arc 12 15 origin (0 0 $zmin)
+    arc 15  9 origin (0 0 $zmin)
+    arc 30 35 origin (0 0 $zmax)
+    arc 35 29 origin (0 0 $zmax)
+    arc 31 34 origin (0 0 $zmax)
+    arc 34 28 origin (0 0 $zmax)
 );
 
 boundary
diff --git a/tutorials/combustion/XiDyMFoam/oscillatingCylinder/system/blockMeshDict b/tutorials/combustion/XiDyMFoam/oscillatingCylinder/system/blockMeshDict
index cc11d554db..ee3b75f0be 100644
--- a/tutorials/combustion/XiDyMFoam/oscillatingCylinder/system/blockMeshDict
+++ b/tutorials/combustion/XiDyMFoam/oscillatingCylinder/system/blockMeshDict
@@ -73,21 +73,22 @@ blocks
 
 edges
 (
-    arc  8  9 (-3.46410161513775 2 0)
-    arc  9 10 (-4                0 0)
-    arc 10 11 (-3.46410161513775 2 0)
-    arc 11 13 ( 0                4 0)
-    arc 13 12 ( 4                0 0)
-    arc 12  8 ( 0               -4 0)
+    arc  8  9 origin (0 0 0)
+    arc  9 10 origin (0 0 0)
+    arc 10 11 origin (0 0 0)
+    arc 11 13 origin (0 0 0)
+    arc 13 12 origin (0 0 0)
+    arc 12  8 origin (0 0 0)
 
-    arc 26 27 (-3.46410161513775 2 1)
-    arc 27 28 (-4                0 1)
-    arc 28 29 (-3.46410161513775 2 1)
-    arc 29 31 ( 0                4 1)
-    arc 31 30 ( 4                0 1)
-    arc 30 26 ( 0               -4 1)
+    arc 26 27 origin (0 0 1)
+    arc 27 28 origin (0 0 1)
+    arc 28 29 origin (0 0 1)
+    arc 29 31 origin (0 0 1)
+    arc 31 30 origin (0 0 1)
+    arc 30 26 origin (0 0 1)
 );
 
+
 defaultPatch
 {
     name    frontAndBack;
diff --git a/tutorials/combustion/fireFoam/LES/compartmentFire/system/blockMeshDict.m4 b/tutorials/combustion/fireFoam/LES/compartmentFire/system/blockMeshDict.m4
index 30801ad7b5..e54dee9a38 100644
--- a/tutorials/combustion/fireFoam/LES/compartmentFire/system/blockMeshDict.m4
+++ b/tutorials/combustion/fireFoam/LES/compartmentFire/system/blockMeshDict.m4
@@ -153,15 +153,15 @@ blocks
 
 edges
 (
-    arc 4 5 (0 0 -R_burner)
-    arc 5 6 (R_burner 0 0)
-    arc 6 7 (0 0 R_burner)
-    arc 7 4 (-R_burner 0 0)
+    arc 4 5 origin (0 0 0)
+    arc 5 6 origin (0 0 0)
+    arc 6 7 origin (0 0 0)
+    arc 7 4 origin (0 0 0)
 
-    arc 24 25 (0 H_box -R_burner)
-    arc 25 26 (R_burner H_box 0)
-    arc 26 27 (0 H_box R_burner)
-    arc 27 24 (-R_burner H_box 0)
+    arc 24 25 origin (0 H_box 0)
+    arc 25 26 origin (0 H_box 0)
+    arc 26 27 origin (0 H_box 0)
+    arc 27 24 origin (0 H_box 0)
 );
 
 boundary
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/externalCoupledSquareBendLiq/system/blockMeshDict b/tutorials/compressible/rhoPimpleFoam/RAS/externalCoupledSquareBendLiq/system/blockMeshDict
index d638986f2e..4aa031189d 100644
--- a/tutorials/compressible/rhoPimpleFoam/RAS/externalCoupledSquareBendLiq/system/blockMeshDict
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/externalCoupledSquareBendLiq/system/blockMeshDict
@@ -16,36 +16,39 @@ FoamFile
 
 scale   0.001;
 
+zmin -25;
+zmax 25;
+
 vertices
 (
-    // front-plane: z = +25mm
+    // front-plane:
     // inlet region
-    (   -50  25   25)           // pt 0
-    (     0  25   25)           // pt 1
-    (   -50  75   25)           // pt 2
-    (     0  75   25)           // pt 3
+    (   -50  25   $zmax)    // pt 0
+    (     0  25   $zmax)    // pt 1
+    (   -50  75   $zmax)    // pt 2
+    (     0  75   $zmax)    // pt 3
     // outlet region
-    (  -500 -75   25)           // pt 4
-    (     0 -75   25)           // pt 5
-    (  -500 -25   25)           // pt 6
-    (     0 -25   25)           // pt 7
+    (  -500 -75   $zmax)    // pt 4
+    (     0 -75   $zmax)    // pt 5
+    (  -500 -25   $zmax)    // pt 6
+    (     0 -25   $zmax)    // pt 7
     // bend mid-points
-    (    25   0   25)           // pt 8
-    (    75   0   25)           // pt 9
-    // back-plane: z = -25mm
+    (    25   0   $zmax)    // pt 8
+    (    75   0   $zmax)    // pt 9
+    // back-plane:
     // inlet region
-    (   -50  25   -25)          // pt 0 + 10
-    (     0  25   -25)          // pt 1 + 10
-    (   -50  75   -25)          // pt 2 + 10
-    (     0  75   -25)          // pt 3 + 10
+    (   -50  25   $zmin)    // pt 0 + 10
+    (     0  25   $zmin)    // pt 1 + 10
+    (   -50  75   $zmin)    // pt 2 + 10
+    (     0  75   $zmin)    // pt 3 + 10
     // outlet region
-    (  -500 -75   -25)          // pt 4 + 10
-    (     0 -75   -25)          // pt 5 + 10
-    (  -500 -25   -25)          // pt 7 + 10
-    (     0 -25   -25)          // pt 8 + 10
+    (  -500 -75   $zmin)    // pt 4 + 10
+    (     0 -75   $zmin)    // pt 5 + 10
+    (  -500 -25   $zmin)    // pt 7 + 10
+    (     0 -25   $zmin)    // pt 8 + 10
     // bend mid-points
-    (    25   0   -25)          // pt 8 + 10
-    (    75   0   -25)          // pt 9 + 10
+    (    25   0   $zmin)    // pt 8 + 10
+    (    75   0   $zmin)    // pt 9 + 10
 );
 
 blocks
@@ -59,16 +62,16 @@ blocks
 
 edges
 (
-   // block 2
-   arc  1  8  ( 17.678  17.678  25)
-   arc 11 18  ( 17.678  17.678 -25)
-   arc  3  9  ( 53.033  53.033  25)
-   arc 13 19  ( 53.033  53.033 -25)
-   // block 3
-   arc  7  8  ( 17.678  -17.678  25)
-   arc 17 18  ( 17.678  -17.678 -25)
-   arc  5  9  ( 53.033  -53.033  25)
-   arc 15 19  ( 53.033  -53.033 -25)
+    // block 2
+    arc  1  8  origin (0 0 $zmax)
+    arc  3  9  origin (0 0 $zmax)
+    arc 11 18  origin (0 0 $zmin)
+    arc 13 19  origin (0 0 $zmin)
+    // block 3
+    arc  7  8  origin (0 0 $zmax)
+    arc  5  9  origin (0 0 $zmax)
+    arc 17 18  origin (0 0 $zmin)
+    arc 15 19  origin (0 0 $zmin)
 );
 
 boundary
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict
index 165f4e9b23..366520c4bb 100644
--- a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict
@@ -16,36 +16,39 @@ FoamFile
 
 scale   0.001;
 
+zmin -25;
+zmax 25;
+
 vertices
 (
-    // front-plane: z = +25mm
+    // front-plane:
     // inlet region
-    (   -50  25   25)           // pt 0
-    (     0  25   25)           // pt 1
-    (   -50  75   25)           // pt 2
-    (     0  75   25)           // pt 3
+    (   -50  25   $zmax)    // pt 0
+    (     0  25   $zmax)    // pt 1
+    (   -50  75   $zmax)    // pt 2
+    (     0  75   $zmax)    // pt 3
     // outlet region
-    (  -500 -75   25)           // pt 4
-    (     0 -75   25)           // pt 5
-    (  -500 -25   25)           // pt 6
-    (     0 -25   25)           // pt 7
+    (  -500 -75   $zmax)    // pt 4
+    (     0 -75   $zmax)    // pt 5
+    (  -500 -25   $zmax)    // pt 6
+    (     0 -25   $zmax)    // pt 7
     // bend mid-points
-    (    25   0   25)           // pt 8
-    (    75   0   25)           // pt 9
-    // back-plane: z = -25mm
+    (    25   0   $zmax)    // pt 8
+    (    75   0   $zmax)    // pt 9
+    // back-plane:
     // inlet region
-    (   -50  25   -25)          // pt 0 + 10
-    (     0  25   -25)          // pt 1 + 10
-    (   -50  75   -25)          // pt 2 + 10
-    (     0  75   -25)          // pt 3 + 10
+    (   -50  25   $zmin)    // pt 0 + 10
+    (     0  25   $zmin)    // pt 1 + 10
+    (   -50  75   $zmin)    // pt 2 + 10
+    (     0  75   $zmin)    // pt 3 + 10
     // outlet region
-    (  -500 -75   -25)          // pt 4 + 10
-    (     0 -75   -25)          // pt 5 + 10
-    (  -500 -25   -25)          // pt 7 + 10
-    (     0 -25   -25)          // pt 8 + 10
+    (  -500 -75   $zmin)    // pt 4 + 10
+    (     0 -75   $zmin)    // pt 5 + 10
+    (  -500 -25   $zmin)    // pt 7 + 10
+    (     0 -25   $zmin)    // pt 8 + 10
     // bend mid-points
-    (    25   0   -25)          // pt 8 + 10
-    (    75   0   -25)          // pt 9 + 10
+    (    25   0   $zmin)    // pt 8 + 10
+    (    75   0   $zmin)    // pt 9 + 10
 );
 
 blocks
@@ -60,15 +63,15 @@ blocks
 edges
 (
     // block 2
-    arc  1  8   (17.678  17.678  25)
-    arc 11 18   (17.678  17.678 -25)
-    arc  3  9   (53.033  53.033  25)
-    arc 13 19   (53.033  53.033 -25)
+    arc  1  8  origin (0 0 $zmax)
+    arc  3  9  origin (0 0 $zmax)
+    arc 11 18  origin (0 0 $zmin)
+    arc 13 19  origin (0 0 $zmin)
     // block 3
-    arc  7  8   (17.678  -17.678  25)
-    arc 17 18   (17.678  -17.678 -25)
-    arc  5  9   (53.033  -53.033  25)
-    arc 15 19   (53.033  -53.033 -25)
+    arc  7  8  origin (0 0 $zmax)
+    arc  5  9  origin (0 0 $zmax)
+    arc 17 18  origin (0 0 $zmin)
+    arc 15 19  origin (0 0 $zmin)
 );
 
 boundary
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/blockMeshDict b/tutorials/compressible/rhoSimpleFoam/squareBend/system/blockMeshDict
index 165f4e9b23..366520c4bb 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/blockMeshDict
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/blockMeshDict
@@ -16,36 +16,39 @@ FoamFile
 
 scale   0.001;
 
+zmin -25;
+zmax 25;
+
 vertices
 (
-    // front-plane: z = +25mm
+    // front-plane:
     // inlet region
-    (   -50  25   25)           // pt 0
-    (     0  25   25)           // pt 1
-    (   -50  75   25)           // pt 2
-    (     0  75   25)           // pt 3
+    (   -50  25   $zmax)    // pt 0
+    (     0  25   $zmax)    // pt 1
+    (   -50  75   $zmax)    // pt 2
+    (     0  75   $zmax)    // pt 3
     // outlet region
-    (  -500 -75   25)           // pt 4
-    (     0 -75   25)           // pt 5
-    (  -500 -25   25)           // pt 6
-    (     0 -25   25)           // pt 7
+    (  -500 -75   $zmax)    // pt 4
+    (     0 -75   $zmax)    // pt 5
+    (  -500 -25   $zmax)    // pt 6
+    (     0 -25   $zmax)    // pt 7
     // bend mid-points
-    (    25   0   25)           // pt 8
-    (    75   0   25)           // pt 9
-    // back-plane: z = -25mm
+    (    25   0   $zmax)    // pt 8
+    (    75   0   $zmax)    // pt 9
+    // back-plane:
     // inlet region
-    (   -50  25   -25)          // pt 0 + 10
-    (     0  25   -25)          // pt 1 + 10
-    (   -50  75   -25)          // pt 2 + 10
-    (     0  75   -25)          // pt 3 + 10
+    (   -50  25   $zmin)    // pt 0 + 10
+    (     0  25   $zmin)    // pt 1 + 10
+    (   -50  75   $zmin)    // pt 2 + 10
+    (     0  75   $zmin)    // pt 3 + 10
     // outlet region
-    (  -500 -75   -25)          // pt 4 + 10
-    (     0 -75   -25)          // pt 5 + 10
-    (  -500 -25   -25)          // pt 7 + 10
-    (     0 -25   -25)          // pt 8 + 10
+    (  -500 -75   $zmin)    // pt 4 + 10
+    (     0 -75   $zmin)    // pt 5 + 10
+    (  -500 -25   $zmin)    // pt 7 + 10
+    (     0 -25   $zmin)    // pt 8 + 10
     // bend mid-points
-    (    25   0   -25)          // pt 8 + 10
-    (    75   0   -25)          // pt 9 + 10
+    (    25   0   $zmin)    // pt 8 + 10
+    (    75   0   $zmin)    // pt 9 + 10
 );
 
 blocks
@@ -60,15 +63,15 @@ blocks
 edges
 (
     // block 2
-    arc  1  8   (17.678  17.678  25)
-    arc 11 18   (17.678  17.678 -25)
-    arc  3  9   (53.033  53.033  25)
-    arc 13 19   (53.033  53.033 -25)
+    arc  1  8  origin (0 0 $zmax)
+    arc  3  9  origin (0 0 $zmax)
+    arc 11 18  origin (0 0 $zmin)
+    arc 13 19  origin (0 0 $zmin)
     // block 3
-    arc  7  8   (17.678  -17.678  25)
-    arc 17 18   (17.678  -17.678 -25)
-    arc  5  9   (53.033  -53.033  25)
-    arc 15 19   (53.033  -53.033 -25)
+    arc  7  8  origin (0 0 $zmax)
+    arc  5  9  origin (0 0 $zmax)
+    arc 17 18  origin (0 0 $zmin)
+    arc 15 19  origin (0 0 $zmin)
 );
 
 boundary
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/blockMeshDict b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/blockMeshDict
index 165f4e9b23..366520c4bb 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/blockMeshDict
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/blockMeshDict
@@ -16,36 +16,39 @@ FoamFile
 
 scale   0.001;
 
+zmin -25;
+zmax 25;
+
 vertices
 (
-    // front-plane: z = +25mm
+    // front-plane:
     // inlet region
-    (   -50  25   25)           // pt 0
-    (     0  25   25)           // pt 1
-    (   -50  75   25)           // pt 2
-    (     0  75   25)           // pt 3
+    (   -50  25   $zmax)    // pt 0
+    (     0  25   $zmax)    // pt 1
+    (   -50  75   $zmax)    // pt 2
+    (     0  75   $zmax)    // pt 3
     // outlet region
-    (  -500 -75   25)           // pt 4
-    (     0 -75   25)           // pt 5
-    (  -500 -25   25)           // pt 6
-    (     0 -25   25)           // pt 7
+    (  -500 -75   $zmax)    // pt 4
+    (     0 -75   $zmax)    // pt 5
+    (  -500 -25   $zmax)    // pt 6
+    (     0 -25   $zmax)    // pt 7
     // bend mid-points
-    (    25   0   25)           // pt 8
-    (    75   0   25)           // pt 9
-    // back-plane: z = -25mm
+    (    25   0   $zmax)    // pt 8
+    (    75   0   $zmax)    // pt 9
+    // back-plane:
     // inlet region
-    (   -50  25   -25)          // pt 0 + 10
-    (     0  25   -25)          // pt 1 + 10
-    (   -50  75   -25)          // pt 2 + 10
-    (     0  75   -25)          // pt 3 + 10
+    (   -50  25   $zmin)    // pt 0 + 10
+    (     0  25   $zmin)    // pt 1 + 10
+    (   -50  75   $zmin)    // pt 2 + 10
+    (     0  75   $zmin)    // pt 3 + 10
     // outlet region
-    (  -500 -75   -25)          // pt 4 + 10
-    (     0 -75   -25)          // pt 5 + 10
-    (  -500 -25   -25)          // pt 7 + 10
-    (     0 -25   -25)          // pt 8 + 10
+    (  -500 -75   $zmin)    // pt 4 + 10
+    (     0 -75   $zmin)    // pt 5 + 10
+    (  -500 -25   $zmin)    // pt 7 + 10
+    (     0 -25   $zmin)    // pt 8 + 10
     // bend mid-points
-    (    25   0   -25)          // pt 8 + 10
-    (    75   0   -25)          // pt 9 + 10
+    (    25   0   $zmin)    // pt 8 + 10
+    (    75   0   $zmin)    // pt 9 + 10
 );
 
 blocks
@@ -60,15 +63,15 @@ blocks
 edges
 (
     // block 2
-    arc  1  8   (17.678  17.678  25)
-    arc 11 18   (17.678  17.678 -25)
-    arc  3  9   (53.033  53.033  25)
-    arc 13 19   (53.033  53.033 -25)
+    arc  1  8  origin (0 0 $zmax)
+    arc  3  9  origin (0 0 $zmax)
+    arc 11 18  origin (0 0 $zmin)
+    arc 13 19  origin (0 0 $zmin)
     // block 3
-    arc  7  8   (17.678  -17.678  25)
-    arc 17 18   (17.678  -17.678 -25)
-    arc  5  9   (53.033  -53.033  25)
-    arc 15 19   (53.033  -53.033 -25)
+    arc  7  8  origin (0 0 $zmax)
+    arc  5  9  origin (0 0 $zmax)
+    arc 17 18  origin (0 0 $zmin)
+    arc 15 19  origin (0 0 $zmin)
 );
 
 boundary
diff --git a/tutorials/electromagnetics/electrostaticFoam/chargedWire/system/blockMeshDict b/tutorials/electromagnetics/electrostaticFoam/chargedWire/system/blockMeshDict
index 0b40b5bfcf..57c872d8a0 100644
--- a/tutorials/electromagnetics/electrostaticFoam/chargedWire/system/blockMeshDict
+++ b/tutorials/electromagnetics/electrostaticFoam/chargedWire/system/blockMeshDict
@@ -53,14 +53,14 @@ blocks
 
 edges
 (
-    arc 0 5 (0.00092387 0.00038268 0)
-    arc 5 10 (0.00038268 0.00092387 0)
-    arc 1 4 (0.0351074 0.0145419 0)
-    arc 4 9 (0.0145419 0.0351074 0)
-    arc 11 16 (0.00092387 0.00038268 0.5)
-    arc 16 21 (0.00038268 0.00092387 0.5)
-    arc 12 15 (0.0351074 0.0145419 0.5)
-    arc 15 20 (0.0145419 0.0351074 0.5)
+    arc 0 5 origin (0 0 0)
+    arc 5 10 origin (0 0 0)
+    arc 1 4 origin (0 0 0)
+    arc 4 9 origin (0 0 0)
+    arc 11 16 origin (0 0 0.5)
+    arc 16 21 origin (0 0 0.5)
+    arc 12 15 origin (0 0 0.5)
+    arc 15 20 origin (0 0 0.5)
 );
 
 boundary
diff --git a/tutorials/finiteArea/liquidFilmFoam/cylinder/system/blockMeshDict b/tutorials/finiteArea/liquidFilmFoam/cylinder/system/blockMeshDict
index aededbfc0f..6b56721ed0 100644
--- a/tutorials/finiteArea/liquidFilmFoam/cylinder/system/blockMeshDict
+++ b/tutorials/finiteArea/liquidFilmFoam/cylinder/system/blockMeshDict
@@ -14,49 +14,54 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Almost identical to basic/potentialFoam/cylinder
+
 scale   0.016;
 
+// Front/back locations
+zmin 0;
+zmax 0.5;
+
 vertices
 (
-    (0.5 0 0)
-    (1 0 0)
-    (2 0 0)
-    (2 0.707107 0)
-    (0.707107 0.707107 0)
-    (0.353553 0.353553 0)
-    (2 2 0)
-    (0.707107 2 0)
-    (0 2 0)
-    (0 1 0)
-    (0 0.5 0)
-    (-0.5 0 0)
-    (-1 0 0)
-    (-2 0 0)
-    (-2 0.707107 0)
-    (-0.707107 0.707107 0)
-    (-0.353553 0.353553 0)
-    (-2 2 0)
-    (-0.707107 2 0)
-
-    (0.5 0 0.5)
-    (1 0 0.5)
-    (2 0 0.5)
-    (2 0.707107 0.5)
-    (0.707107 0.707107 0.5)
-    (0.353553 0.353553 0.5)
-    (2 2 0.5)
-    (0.707107 2 0.5)
-    (0 2 0.5)
-    (0 1 0.5)
-    (0 0.5 0.5)
-    (-0.5 0 0.5)
-    (-1 0 0.5)
-    (-2 0 0.5)
-    (-2 0.707107 0.5)
-    (-0.707107 0.707107 0.5)
-    (-0.353553 0.353553 0.5)
-    (-2 2 0.5)
-    (-0.707107 2 0.5)
+    (0.5 0 $zmin)
+    (1 0 $zmin)
+    (2 0 $zmin)
+    (2 0.707107 $zmin)
+    (0.707107 0.707107 $zmin)
+    (0.353553 0.353553 $zmin)
+    (2 2 $zmin)
+    (0.707107 2 $zmin)
+    (0 2 $zmin)
+    (0 1 $zmin)
+    (0 0.5 $zmin)
+    (-0.5 0 $zmin)
+    (-1 0 $zmin)
+    (-2 0 $zmin)
+    (-2 0.707107 $zmin)
+    (-0.707107 0.707107 $zmin)
+    (-0.353553 0.353553 $zmin)
+    (-2 2 $zmin)
+    (-0.707107 2 $zmin)
+    (0.5 0 $zmax)
+    (1 0 $zmax)
+    (2 0 $zmax)
+    (2 0.707107 $zmax)
+    (0.707107 0.707107 $zmax)
+    (0.353553 0.353553 $zmax)
+    (2 2 $zmax)
+    (0.707107 2 $zmax)
+    (0 2 $zmax)
+    (0 1 $zmax)
+    (0 0.5 $zmax)
+    (-0.5 0 $zmax)
+    (-1 0 $zmax)
+    (-2 0 $zmax)
+    (-2 0.707107 $zmax)
+    (-0.707107 0.707107 $zmax)
+    (-0.353553 0.353553 $zmax)
+    (-2 2 $zmax)
+    (-0.707107 2 $zmax)
 );
 
 blocks
@@ -75,22 +80,24 @@ blocks
 
 edges
 (
-    arc 0 5 (0.469846 0.17101 0)
-    arc 5 10 (0.17101 0.469846 0)
-    arc 1 4 (0.939693 0.34202 0)
-    arc 4 9 (0.34202 0.939693 0)
-    arc 19 24 (0.469846 0.17101 0.5)
-    arc 24 29 (0.17101 0.469846 0.5)
-    arc 20 23 (0.939693 0.34202 0.5)
-    arc 23 28 (0.34202 0.939693 0.5)
-    arc 11 16 (-0.469846 0.17101 0)
-    arc 16 10 (-0.17101 0.469846 0)
-    arc 12 15 (-0.939693 0.34202 0)
-    arc 15 9 (-0.34202 0.939693 0)
-    arc 30 35 (-0.469846 0.17101 0.5)
-    arc 35 29 (-0.17101 0.469846 0.5)
-    arc 31 34 (-0.939693 0.34202 0.5)
-    arc 34 28 (-0.34202 0.939693 0.5)
+    // Inner cylinder
+    arc  0  5 origin (0 0 $zmin)
+    arc  5 10 origin (0 0 $zmin)
+    arc  1  4 origin (0 0 $zmin)
+    arc  4  9 origin (0 0 $zmin)
+    arc 19 24 origin (0 0 $zmax)
+    arc 24 29 origin (0 0 $zmax)
+    arc 20 23 origin (0 0 $zmax)
+    arc 23 28 origin (0 0 $zmax)
+    // Intermediate cylinder
+    arc 11 16 origin (0 0 $zmin)
+    arc 16 10 origin (0 0 $zmin)
+    arc 12 15 origin (0 0 $zmin)
+    arc 15  9 origin (0 0 $zmin)
+    arc 30 35 origin (0 0 $zmax)
+    arc 35 29 origin (0 0 $zmax)
+    arc 31 34 origin (0 0 $zmax)
+    arc 34 28 origin (0 0 $zmax)
 );
 
 patches
diff --git a/tutorials/incompressible/SRFSimpleFoam/mixer/system/blockMeshDict b/tutorials/incompressible/SRFSimpleFoam/mixer/system/blockMeshDict
index 56729cdf23..578ecbf576 100644
--- a/tutorials/incompressible/SRFSimpleFoam/mixer/system/blockMeshDict
+++ b/tutorials/incompressible/SRFSimpleFoam/mixer/system/blockMeshDict
@@ -17,57 +17,60 @@ FoamFile
 
 scale   0.1;
 
+zmin 0;
+zmax 2;
+
 vertices
 (
-      (       0.500        0.000        0.000)
-      (       0.369        0.338        0.000)
-      (       0.338        0.369        0.000)
-      (       0.000        0.500        0.000)
-      (       0.737        0.676        0.000)
-      (       0.074        0.068        0.000)
-      (       0.676        0.737        0.000)
-      (       0.068        0.074        0.000)
-      (       0.000        1.000        0.000)
-      (       1.000        0.000        0.000)
-      (       0.100        0.000        0.000)
-      (       0.000        0.100        0.000)
-      (       0.500        0.000        2.000)
-      (       0.369        0.338        2.000)
-      (       0.338        0.369        2.000)
-      (       0.000        0.500        2.000)
-      (       0.737        0.676        2.000)
-      (       0.074        0.068        2.000)
-      (       0.676        0.737        2.000)
-      (       0.068        0.074        2.000)
-      (       0.000        1.000        2.000)
-      (       1.000        0.000        2.000)
-      (       0.100        0.000        2.000)
-      (       0.000        0.100        2.000)
+    (0.500        0.000    $zmin)
+    (0.369        0.338    $zmin)
+    (0.338        0.369    $zmin)
+    (0.000        0.500    $zmin)
+    (0.737        0.676    $zmin)
+    (0.074        0.068    $zmin)
+    (0.676        0.737    $zmin)
+    (0.068        0.074    $zmin)
+    (0.000        1.000    $zmin)
+    (1.000        0.000    $zmin)
+    (0.100        0.000    $zmin)
+    (0.000        0.100    $zmin)
+    (0.500        0.000    $zmax)
+    (0.369        0.338    $zmax)
+    (0.338        0.369    $zmax)
+    (0.000        0.500    $zmax)
+    (0.737        0.676    $zmax)
+    (0.074        0.068    $zmax)
+    (0.676        0.737    $zmax)
+    (0.068        0.074    $zmax)
+    (0.000        1.000    $zmax)
+    (1.000        0.000    $zmax)
+    (0.100        0.000    $zmax)
+    (0.000        0.100    $zmax)
 );
 
 blocks
 (
-      hex (1 0 9 4 13 12 21 16) (10 20 40) simpleGrading (1 1 1)
-      hex (2 1 4 6 14 13 16 18) (2 20 40) simpleGrading (1 1 1)
-      hex (3 2 6 8 15 14 18 20) (10 20 40) simpleGrading (1 1 1)
-      hex (5 10 0 1 17 22 12 13) (10 20 40) simpleGrading (1 1 1)
-      hex (11 7 2 3 23 19 14 15) (10 20 40) simpleGrading (1 1 1)
+    hex (1 0 9 4 13 12 21 16) (10 20 40) simpleGrading (1 1 1)
+    hex (2 1 4 6 14 13 16 18) (2 20 40) simpleGrading (1 1 1)
+    hex (3 2 6 8 15 14 18 20) (10 20 40) simpleGrading (1 1 1)
+    hex (5 10 0 1 17 22 12 13) (10 20 40) simpleGrading (1 1 1)
+    hex (11 7 2 3 23 19 14 15) (10 20 40) simpleGrading (1 1 1)
 );
 
 edges
 (
-      arc  0 1 (       0.470        0.171        0.000 )
-      arc  12 13 (       0.470        0.171        2.000 )
-      arc  2 3 (       0.171        0.470        0.000 )
-      arc  14 15 (       0.171        0.470        2.000 )
-      arc  9 4 (       0.940        0.342        0.000 )
-      arc  21 16 (       0.940        0.342        2.000 )
-      arc  5 10 (       0.094        0.034        0.000 )
-      arc  17 22 (       0.094        0.034        2.000 )
-      arc  6 8 (       0.342        0.940        0.000 )
-      arc  18 20 (       0.342        0.940        2.000 )
-      arc  11 7 (       0.034        0.094        0.000 )
-      arc  23 19 (       0.034        0.094        2.000 )
+    arc   0  1 origin (0 0 $zmin)
+    arc  12 13 origin (0 0 $zmax)
+    arc   2  3 origin (0 0 $zmin)
+    arc  14 15 origin (0 0 $zmax)
+    arc   9  4 origin (0 0 $zmin)
+    arc  21 16 origin (0 0 $zmax)
+    arc   5 10 origin (0 0 $zmin)
+    arc  17 22 origin (0 0 $zmax)
+    arc   6  8 origin (0 0 $zmin)
+    arc  18 20 origin (0 0 $zmax)
+    arc  11  7 origin (0 0 $zmin)
+    arc  23 19 origin (0 0 $zmax)
 );
 
 boundary
diff --git a/tutorials/incompressible/nonNewtonianIcoFoam/offsetCylinder/system/blockMeshDict b/tutorials/incompressible/nonNewtonianIcoFoam/offsetCylinder/system/blockMeshDict
index 2ec0e7a55a..5ec15e9252 100644
--- a/tutorials/incompressible/nonNewtonianIcoFoam/offsetCylinder/system/blockMeshDict
+++ b/tutorials/incompressible/nonNewtonianIcoFoam/offsetCylinder/system/blockMeshDict
@@ -16,72 +16,75 @@ FoamFile
 
 scale   1;
 
+zmin -1;
+zmax 1;
+
 vertices
 (
-    (1 0 -1)
-    (1.4 0 -1)
-    (5 0 -1)
-    (5 0.989949 -1)
-    (0.989949 0.989949 -1)
-    (0.707107 0.707107 -1)
-    (5 2.5 -1)
-    (0.989949 2.5 -1)
-    (0 2.5 -1)
-    (0 1.4 -1)
-    (0 1 -1)
-    (-1 0 -1)
-    (-1.4 0 -1)
-    (-5 0 -1)
-    (-5 0.989949 -1)
-    (-0.989949 0.989949 -1)
-    (-0.707107 0.707107 -1)
-    (-5 2.5 -1)
-    (-0.989949 2.5 -1)
-    (5 -0.989949 -1)
-    (0.989949 -0.989949 -1)
-    (0.707107 -0.707107 -1)
-    (5 -1.5 -1)
-    (0.989949 -1.5 -1)
-    (0 -1.5 -1)
-    (0 -1.4 -1)
-    (0 -1 -1)
-    (-5 -0.989949 -1)
-    (-0.989949 -0.989949 -1)
-    (-0.707107 -0.707107 -1)
-    (-5 -1.5 -1)
-    (-0.989949 -1.5 -1)
-    (1 0 1)
-    (1.4 0 1)
-    (5 0 1)
-    (5 0.989949 1)
-    (0.989949 0.989949 1)
-    (0.707107 0.707107 1)
-    (5 2.5 1)
-    (0.989949 2.5 1)
-    (0 2.5 1)
-    (0 1.4 1)
-    (0 1 1)
-    (-1 0 1)
-    (-1.4 0 1)
-    (-5 0 1)
-    (-5 0.989949 1)
-    (-0.989949 0.989949 1)
-    (-0.707107 0.707107 1)
-    (-5 2.5 1)
-    (-0.989949 2.5 1)
-    (5 -0.989949 1)
-    (0.989949 -0.989949 1)
-    (0.707107 -0.707107 1)
-    (5 -1.5 1)
-    (0.989949 -1.5 1)
-    (0 -1.5 1)
-    (0 -1.4 1)
-    (0 -1 1)
-    (-5 -0.989949 1)
-    (-0.989949 -0.989949 1)
-    (-0.707107 -0.707107 1)
-    (-5 -1.5 1)
-    (-0.989949 -1.5 1)
+    (1 0 $zmin)
+    (1.4 0 $zmin)
+    (5 0 $zmin)
+    (5 0.989949 $zmin)
+    (0.989949 0.989949 $zmin)
+    (0.707107 0.707107 $zmin)
+    (5 2.5 $zmin)
+    (0.989949 2.5 $zmin)
+    (0 2.5 $zmin)
+    (0 1.4 $zmin)
+    (0 1 $zmin)
+    (-1 0 $zmin)
+    (-1.4 0 $zmin)
+    (-5 0 $zmin)
+    (-5 0.989949 $zmin)
+    (-0.989949 0.989949 $zmin)
+    (-0.707107 0.707107 $zmin)
+    (-5 2.5 $zmin)
+    (-0.989949 2.5 $zmin)
+    (5 -0.989949 $zmin)
+    (0.989949 -0.989949 $zmin)
+    (0.707107 -0.707107 $zmin)
+    (5 -1.5 $zmin)
+    (0.989949 -1.5 $zmin)
+    (0 -1.5 $zmin)
+    (0 -1.4 $zmin)
+    (0 -1 $zmin)
+    (-5 -0.989949 $zmin)
+    (-0.989949 -0.989949 $zmin)
+    (-0.707107 -0.707107 $zmin)
+    (-5 -1.5 $zmin)
+    (-0.989949 -1.5 $zmin)
+    (1 0 $zmax)
+    (1.4 0 $zmax)
+    (5 0 $zmax)
+    (5 0.989949 $zmax)
+    (0.989949 0.989949 $zmax)
+    (0.707107 0.707107 $zmax)
+    (5 2.5 $zmax)
+    (0.989949 2.5 $zmax)
+    (0 2.5 $zmax)
+    (0 1.4 $zmax)
+    (0 1 $zmax)
+    (-1 0 $zmax)
+    (-1.4 0 $zmax)
+    (-5 0 $zmax)
+    (-5 0.989949 $zmax)
+    (-0.989949 0.989949 $zmax)
+    (-0.707107 0.707107 $zmax)
+    (-5 2.5 $zmax)
+    (-0.989949 2.5 $zmax)
+    (5 -0.989949 $zmax)
+    (0.989949 -0.989949 $zmax)
+    (0.707107 -0.707107 $zmax)
+    (5 -1.5 $zmax)
+    (0.989949 -1.5 $zmax)
+    (0 -1.5 $zmax)
+    (0 -1.4 $zmax)
+    (0 -1 $zmax)
+    (-5 -0.989949 $zmax)
+    (-0.989949 -0.989949 $zmax)
+    (-0.707107 -0.707107 $zmax)
+    (-5 -1.5 $zmax)
+    (-0.989949 -1.5 $zmax)
 );
 
 blocks
@@ -110,38 +113,38 @@ blocks
 
 edges
 (
-    arc 5 0 (0.92388 0.382683 -1)
-    arc 5 10 (0.382683 0.923879 -1)
-    arc 1 4 (1.29343 0.535757 -1)
-    arc 4 9 (0.535757 1.29343 -1)
-    arc 32 37 (0.92388 0.382683 1)
-    arc 37 42 (0.382683 0.923879 1)
-    arc 33 36 (1.29343 0.535757 1)
-    arc 36 41 (0.535757 1.29343 1)
-    arc 11 16 (-0.923879 0.382683 -1)
-    arc 16 10 (-0.382683 0.923879 -1)
-    arc 12 15 (-1.29343 0.535757 -1)
-    arc 15 9 (-0.535757 1.29343 -1)
-    arc 43 48 (-0.923879 0.382683 1)
-    arc 48 42 (-0.382683 0.923879 1)
-    arc 44 47 (-1.29343 0.535757 1)
-    arc 47 41 (-0.535757 1.29343 1)
-    arc 0 21 (0.923879 -0.382683 -1)
-    arc 21 26 (0.382683 -0.923879 -1)
-    arc 1 20 (1.29343 -0.535757 -1)
-    arc 20 25 (0.535757 -1.29343 -1)
-    arc 32 53 (0.923879 -0.382683 1)
-    arc 53 58 (0.382683 -0.923879 1)
-    arc 33 52 (1.29343 -0.535757 1)
-    arc 52 57 (0.535757 -1.29343 1)
-    arc 11 29 (-0.923879 -0.382683 -1)
-    arc 29 26 (-0.382683 -0.923879 -1)
-    arc 12 28 (-1.29343 -0.535757 -1)
-    arc 28 25 (-0.535757 -1.29343 -1)
-    arc 43 61 (-0.923879 -0.382683 1)
-    arc 61 58 (-0.382683 -0.923879 1)
-    arc 44 60 (-1.29343 -0.535757 1)
-    arc 60 57 (-0.535757 -1.29343 1)
+    arc 5 0 origin (0 0 $zmin)
+    arc 5 10 origin (0 0 $zmin)
+    arc 1 4 origin (0 0 $zmin)
+    arc 4 9 origin (0 0 $zmin)
+    arc 32 37 origin (0 0 $zmax)
+    arc 37 42 origin (0 0 $zmax)
+    arc 33 36 origin (0 0 $zmax)
+    arc 36 41 origin (0 0 $zmax)
+    arc 11 16 origin (0 0 $zmin)
+    arc 16 10 origin (0 0 $zmin)
+    arc 12 15 origin (0 0 $zmin)
+    arc 15 9 origin (0 0 $zmin)
+    arc 43 48 origin (0 0 $zmax)
+    arc 48 42 origin (0 0 $zmax)
+    arc 44 47 origin (0 0 $zmax)
+    arc 47 41 origin (0 0 $zmax)
+    arc 0 21 origin (0 0 $zmin)
+    arc 21 26 origin (0 0 $zmin)
+    arc 1 20 origin (0 0 $zmin)
+    arc 20 25 origin (0 0 $zmin)
+    arc 32 53 origin (0 0 $zmax)
+    arc 53 58 origin (0 0 $zmax)
+    arc 33 52 origin (0 0 $zmax)
+    arc 52 57 origin (0 0 $zmax)
+    arc 11 29 origin (0 0 $zmin)
+    arc 29 26 origin (0 0 $zmin)
+    arc 12 28 origin (0 0 $zmin)
+    arc 28 25 origin (0 0 $zmin)
+    arc 43 61 origin (0 0 $zmax)
+    arc 61 58 origin (0 0 $zmax)
+    arc 44 60 origin (0 0 $zmax)
+    arc 60 57 origin (0 0 $zmax)
 );
 
 boundary
diff --git a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.main b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.main
index f6c2adc8f7..223dc50e5a 100644
--- a/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.main
+++ b/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.main
@@ -23,18 +23,6 @@ x3        #eval{ 10.0*$R };
 xOutlet   #eval{ 18.6667*$R };
 xInlet    #eval{ -10.125*$R };
 
-RsinPi8   #eval{ $R*sin(0.125*pi()) };
-RsinPi8n  #eval{ -$RsinPi8 };
-RcosPi8   #eval{ $R*cos(0.125*pi()) };
-RcosPi8n  #eval{ -$RcosPi8 };
-RsinPi4   #eval{ $R*sin(0.25*pi()) };
-
-x2sinPi8  #eval{ $x2*sin(0.125*pi()) };
-x2sinPi8n #eval{ -$x2sinPi8 };
-x2cosPi8  #eval{ $x2*cos(0.125*pi()) };
-x2cosPi8n #eval{ -$x2cosPi8 };
-x2sinPi4  #eval{ $x2*sin(0.25*pi()) };
-
 z0        -0.0075;
 z1        0.0075;
 nz        1;
@@ -49,15 +37,18 @@ vertices #codeStream
 
     code
     #{
+        // sin(45), cos(45)
+        const scalar sqrt05 = sqrt(0.5);
+
         pointField points(19);
         points[0] = point($R, 0, $z0);
         points[1] = point($x2, 0, $z0);
         points[2] = point($x3, 0, $z0);
-        points[3] = point($x3, $x2sinPi4, $z0);
-        points[4] = point($x2sinPi4, $x2sinPi4, $z0);
-        points[5] = point($RsinPi4, $RsinPi4, $z0);
+        points[3] = point($x3, $x2*sqrt05, $z0);
+        points[4] = point($x2*sqrt05, $x2*sqrt05, $z0);
+        points[5] = point($R*sqrt05, $R*sqrt05, $z0);
         points[6] = point($x3, $x3, $z0);
-        points[7] = point($x2sinPi4, $x3, $z0);
+        points[7] = point($x2*sqrt05, $x3, $z0);
 
         // Mirror +x points to -x side
         points[11] = point(-points[0].x(), points[0].y(), points[0].z());
@@ -76,7 +67,7 @@ vertices #codeStream
 
         // Mirror -z points to +z side
         label sz = points.size();
-        points.setSize(2*sz);
+        points.resize(2*sz);
         for (label i = 0; i < sz; ++i)
         {
             const point& pt = points[i];
@@ -86,7 +77,7 @@ vertices #codeStream
         // Add an inner cylinder
         sz = points.size();
         label nAdd = 6;
-        points.setSize(sz + nAdd);
+        points.resize(sz + nAdd);
 
         // Points within the inner cylinder
         points[sz] = point(0, 0, $z0);
@@ -102,7 +93,7 @@ vertices #codeStream
 
         // Mirror -z points to +z side
         sz = points.size();
-        points.setSize(sz + nAdd);
+        points.resize(sz + nAdd);
         for (label i = 0; i < nAdd; ++i)
         {
             const point& pt = points[i+sz-nAdd];
@@ -112,21 +103,21 @@ vertices #codeStream
         // Add downstream and upstream blocks
         sz = points.size();
         nAdd = 6;
-        points.setSize(sz + nAdd);
+        points.resize(sz + nAdd);
 
         // Points on outlet
         points[sz] = point($xOutlet, 0, $z0);
         points[sz + 1] = point($xOutlet, $x3, $z0);
-        points[sz + 4] = point($xOutlet, $x2sinPi4, $z0);
+        points[sz + 4] = point($xOutlet, $x2*sqrt05, $z0);
 
         // Points on inlet
         points[sz + 2] = point($xInlet, 0, $z0);
         points[sz + 3] = point($xInlet, $x3, $z0);
-        points[sz + 5] = point($xInlet, $x2sinPi4, $z0);
+        points[sz + 5] = point($xInlet, $x2*sqrt05, $z0);
 
         // Mirror -z points to +z side
         sz = points.size();
-        points.setSize(sz + nAdd);
+        points.resize(sz + nAdd);
         for (label i = 0; i < nAdd; ++i)
         {
             const point& pt = points[i + sz - nAdd];
@@ -158,22 +149,22 @@ blocks
 
 edges
 (
-    arc  0  5 ($RcosPi8   $RsinPi8  $z0)
-    arc  5 10 ($RsinPi8   $RcosPi8  $z0)
-    arc  1  4 ($x2cosPi8  $x2sinPi8 $z0)
-    arc  4  9 ($x2sinPi8  $x2cosPi8 $z0)
-    arc 19 24 ($RcosPi8   $RsinPi8  $z1)
-    arc 24 29 ($RsinPi8   $RcosPi8  $z1)
-    arc 20 23 ($x2cosPi8  $x2sinPi8 $z1)
-    arc 23 28 ($x2sinPi8  $x2cosPi8 $z1)
-    arc 11 16 ($RcosPi8n  $RsinPi8  $z0)
-    arc 16 10 ($RsinPi8n  $RcosPi8  $z0)
-    arc 12 15 ($x2cosPi8n $x2sinPi8 $z0)
-    arc 15  9 ($x2sinPi8n $x2cosPi8 $z0)
-    arc 30 35 ($RcosPi8n  $RsinPi8  $z1)
-    arc 35 29 ($RsinPi8n  $RcosPi8  $z1)
-    arc 31 34 ($x2cosPi8n $x2sinPi8 $z1)
-    arc 34 28 ($x2sinPi8n $x2cosPi8 $z1)
+    arc  0  5 origin (0 0 $z0)
+    arc  5 10 origin (0 0 $z0)
+    arc  1  4 origin (0 0 $z0)
+    arc  4  9 origin (0 0 $z0)
+    arc 19 24 origin (0 0 $z1)
+    arc 24 29 origin (0 0 $z1)
+    arc 20 23 origin (0 0 $z1)
+    arc 23 28 origin (0 0 $z1)
+    arc 11 16 origin (0 0 $z0)
+    arc 16 10 origin (0 0 $z0)
+    arc 12 15 origin (0 0 $z0)
+    arc 15  9 origin (0 0 $z0)
+    arc 30 35 origin (0 0 $z1)
+    arc 35 29 origin (0 0 $z1)
+    arc 31 34 origin (0 0 $z1)
+    arc 34 28 origin (0 0 $z1)
 );
 
 boundary
diff --git a/tutorials/incompressible/simpleFoam/bump2D/system/blockMeshDict b/tutorials/incompressible/simpleFoam/bump2D/system/blockMeshDict
index 890fc02425..5d2a59b5cd 100644
--- a/tutorials/incompressible/simpleFoam/bump2D/system/blockMeshDict
+++ b/tutorials/incompressible/simpleFoam/bump2D/system/blockMeshDict
@@ -74,14 +74,14 @@ edges #codeStream
 
     code
     #{
-        const scalar xMin = 0.3;
-        const scalar xMax = 1.2;
-        const label nPoints = 100;
-        const scalar dx = (xMax - xMin)/scalar(nPoints - 1);
+        constexpr scalar xMin = 0.3;
+        constexpr scalar xMax = 1.2;
+        constexpr label nPoints = 100;
+        constexpr scalar dx = (xMax - xMin)/scalar(nPoints - 1);
+        constexpr scalar pi = constant::mathematical::pi;
 
         os  << "(" << nl << "spline 2 3" << nl;
         pointField profile(nPoints);
-        const scalar pi = constant::mathematical::pi;
 
         for (label i = 0; i < nPoints; ++i)
         {
diff --git a/tutorials/incompressible/simpleFoam/pipeCyclic/system/blockMeshDict b/tutorials/incompressible/simpleFoam/pipeCyclic/system/blockMeshDict
index bdbbf798dd..b97babdd44 100644
--- a/tutorials/incompressible/simpleFoam/pipeCyclic/system/blockMeshDict
+++ b/tutorials/incompressible/simpleFoam/pipeCyclic/system/blockMeshDict
@@ -25,34 +25,33 @@ halfAngle 45.0;
 radius  0.5;
 
 y               #eval{ $radius*sin(degToRad($halfAngle)) };
-minY            #eval{ -1*$y };
 z               #eval{ $radius*cos(degToRad($halfAngle)) };
+minY            #eval{ -1*$y };
 minZ            #eval{ -1*$z };
 
 vertices
 (
-    (0.0    0.0 0)      //0
-    (10     0.0 0)
-    (10     0.0 0)      //2
-    (0.0    0.0 0)
+    (0  0 0)        //0
+    (10 0 0)        //1
+    (10 0 0)        //2
+    (0  0 0)        //3
 
-    (0.0    $minY $z)   //4
-    (10     $minY $z)
-    (10     $y $z)      //6
-    (0.0    $y $z)
+    (0  $minY $z)   //4
+    (10 $minY $z)   //5
+    (10 $y $z)      //6
+    (0  $y $z)      //7
 
 );
 
 blocks
 (
-    // inlet block
-    hex (0 1 2 3  4 5 6 7) (50 5 5) simpleGrading (1 1 1)
+    hex (0 1 2 3 4 5 6 7) (50 5 5) simpleGrading (1 1 1)
 );
 
 edges
 (
-    arc 4 7 (0 0 $radius)
-    arc 5 6 (10 0 $radius)
+    arc 4 7 origin (0 0 0)
+    arc 5 6 origin (10 0 0)
 );
 
 boundary
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/blockMeshDict b/tutorials/incompressible/simpleFoam/squareBend/system/blockMeshDict
index 165f4e9b23..366520c4bb 100644
--- a/tutorials/incompressible/simpleFoam/squareBend/system/blockMeshDict
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/blockMeshDict
@@ -16,36 +16,39 @@ FoamFile
 
 scale   0.001;
 
+zmin -25;
+zmax 25;
+
 vertices
 (
-    // front-plane: z = +25mm
+    // front-plane:
     // inlet region
-    (   -50  25   25)           // pt 0
-    (     0  25   25)           // pt 1
-    (   -50  75   25)           // pt 2
-    (     0  75   25)           // pt 3
+    (   -50  25   $zmax)    // pt 0
+    (     0  25   $zmax)    // pt 1
+    (   -50  75   $zmax)    // pt 2
+    (     0  75   $zmax)    // pt 3
     // outlet region
-    (  -500 -75   25)           // pt 4
-    (     0 -75   25)           // pt 5
-    (  -500 -25   25)           // pt 6
-    (     0 -25   25)           // pt 7
+    (  -500 -75   $zmax)    // pt 4
+    (     0 -75   $zmax)    // pt 5
+    (  -500 -25   $zmax)    // pt 6
+    (     0 -25   $zmax)    // pt 7
     // bend mid-points
-    (    25   0   25)           // pt 8
-    (    75   0   25)           // pt 9
-    // back-plane: z = -25mm
+    (    25   0   $zmax)    // pt 8
+    (    75   0   $zmax)    // pt 9
+    // back-plane:
     // inlet region
-    (   -50  25   -25)          // pt 0 + 10
-    (     0  25   -25)          // pt 1 + 10
-    (   -50  75   -25)          // pt 2 + 10
-    (     0  75   -25)          // pt 3 + 10
+    (   -50  25   $zmin)    // pt 0 + 10
+    (     0  25   $zmin)    // pt 1 + 10
+    (   -50  75   $zmin)    // pt 2 + 10
+    (     0  75   $zmin)    // pt 3 + 10
     // outlet region
-    (  -500 -75   -25)          // pt 4 + 10
-    (     0 -75   -25)          // pt 5 + 10
-    (  -500 -25   -25)          // pt 7 + 10
-    (     0 -25   -25)          // pt 8 + 10
+    (  -500 -75   $zmin)    // pt 4 + 10
+    (     0 -75   $zmin)    // pt 5 + 10
+    (  -500 -25   $zmin)    // pt 7 + 10
+    (     0 -25   $zmin)    // pt 8 + 10
     // bend mid-points
-    (    25   0   -25)          // pt 8 + 10
-    (    75   0   -25)          // pt 9 + 10
+    (    25   0   $zmin)    // pt 8 + 10
+    (    75   0   $zmin)    // pt 9 + 10
 );
 
 blocks
@@ -60,15 +63,15 @@ blocks
 edges
 (
     // block 2
-    arc  1  8   (17.678  17.678  25)
-    arc 11 18   (17.678  17.678 -25)
-    arc  3  9   (53.033  53.033  25)
-    arc 13 19   (53.033  53.033 -25)
+    arc  1  8  origin (0 0 $zmax)
+    arc  3  9  origin (0 0 $zmax)
+    arc 11 18  origin (0 0 $zmin)
+    arc 13 19  origin (0 0 $zmin)
     // block 3
-    arc  7  8   (17.678  -17.678  25)
-    arc 17 18   (17.678  -17.678 -25)
-    arc  5  9   (53.033  -53.033  25)
-    arc 15 19   (53.033  -53.033 -25)
+    arc  7  8  origin (0 0 $zmax)
+    arc  5  9  origin (0 0 $zmax)
+    arc 17 18  origin (0 0 $zmin)
+    arc 15 19  origin (0 0 $zmin)
 );
 
 boundary
diff --git a/tutorials/lagrangian/reactingParcelFoam/cylinder/system/blockMeshDict b/tutorials/lagrangian/reactingParcelFoam/cylinder/system/blockMeshDict
index 7445c24428..2161b74b1c 100644
--- a/tutorials/lagrangian/reactingParcelFoam/cylinder/system/blockMeshDict
+++ b/tutorials/lagrangian/reactingParcelFoam/cylinder/system/blockMeshDict
@@ -14,48 +14,54 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Almost identical to basic/potentialFoam/cylinder
+
 scale   1;
 
+// Front/back locations
+zmin -0.5;
+zmax 0.5;
+
 vertices
 (
-    (0.5 0 -0.5)
-    (1 0 -0.5)
-    (2 0 -0.5)
-    (2 0.707107 -0.5)
-    (0.707107 0.707107 -0.5)
-    (0.353553 0.353553 -0.5)
-    (2 2 -0.5)
-    (0.707107 2 -0.5)
-    (0 2 -0.5)
-    (0 1 -0.5)
-    (0 0.5 -0.5)
-    (-0.5 0 -0.5)
-    (-1 0 -0.5)
-    (-2 0 -0.5)
-    (-2 0.707107 -0.5)
-    (-0.707107 0.707107 -0.5)
-    (-0.353553 0.353553 -0.5)
-    (-2 2 -0.5)
-    (-0.707107 2 -0.5)
-    (0.5 0 0.5)
-    (1 0 0.5)
-    (2 0 0.5)
-    (2 0.707107 0.5)
-    (0.707107 0.707107 0.5)
-    (0.353553 0.353553 0.5)
-    (2 2 0.5)
-    (0.707107 2 0.5)
-    (0 2 0.5)
-    (0 1 0.5)
-    (0 0.5 0.5)
-    (-0.5 0 0.5)
-    (-1 0 0.5)
-    (-2 0 0.5)
-    (-2 0.707107 0.5)
-    (-0.707107 0.707107 0.5)
-    (-0.353553 0.353553 0.5)
-    (-2 2 0.5)
-    (-0.707107 2 0.5)
+    (0.5 0 $zmin)
+    (1 0 $zmin)
+    (2 0 $zmin)
+    (2 0.707107 $zmin)
+    (0.707107 0.707107 $zmin)
+    (0.353553 0.353553 $zmin)
+    (2 2 $zmin)
+    (0.707107 2 $zmin)
+    (0 2 $zmin)
+    (0 1 $zmin)
+    (0 0.5 $zmin)
+    (-0.5 0 $zmin)
+    (-1 0 $zmin)
+    (-2 0 $zmin)
+    (-2 0.707107 $zmin)
+    (-0.707107 0.707107 $zmin)
+    (-0.353553 0.353553 $zmin)
+    (-2 2 $zmin)
+    (-0.707107 2 $zmin)
+    (0.5 0 $zmax)
+    (1 0 $zmax)
+    (2 0 $zmax)
+    (2 0.707107 $zmax)
+    (0.707107 0.707107 $zmax)
+    (0.353553 0.353553 $zmax)
+    (2 2 $zmax)
+    (0.707107 2 $zmax)
+    (0 2 $zmax)
+    (0 1 $zmax)
+    (0 0.5 $zmax)
+    (-0.5 0 $zmax)
+    (-1 0 $zmax)
+    (-2 0 $zmax)
+    (-2 0.707107 $zmax)
+    (-0.707107 0.707107 $zmax)
+    (-0.353553 0.353553 $zmax)
+    (-2 2 $zmax)
+    (-0.707107 2 $zmax)
 );
 
 blocks
@@ -74,22 +80,24 @@ blocks
 
 edges
 (
-    arc 0 5 (0.469846 0.17101 -0.5)
-    arc 5 10 (0.17101 0.469846 -0.5)
-    arc 1 4 (0.939693 0.34202 -0.5)
-    arc 4 9 (0.34202 0.939693 -0.5)
-    arc 19 24 (0.469846 0.17101 0.5)
-    arc 24 29 (0.17101 0.469846 0.5)
-    arc 20 23 (0.939693 0.34202 0.5)
-    arc 23 28 (0.34202 0.939693 0.5)
-    arc 11 16 (-0.469846 0.17101 -0.5)
-    arc 16 10 (-0.17101 0.469846 -0.5)
-    arc 12 15 (-0.939693 0.34202 -0.5)
-    arc 15 9 (-0.34202 0.939693 -0.5)
-    arc 30 35 (-0.469846 0.17101 0.5)
-    arc 35 29 (-0.17101 0.469846 0.5)
-    arc 31 34 (-0.939693 0.34202 0.5)
-    arc 34 28 (-0.34202 0.939693 0.5)
+    // Inner cylinder
+    arc  0  5 origin (0 0 $zmin)
+    arc  5 10 origin (0 0 $zmin)
+    arc  1  4 origin (0 0 $zmin)
+    arc  4  9 origin (0 0 $zmin)
+    arc 19 24 origin (0 0 $zmax)
+    arc 24 29 origin (0 0 $zmax)
+    arc 20 23 origin (0 0 $zmax)
+    arc 23 28 origin (0 0 $zmax)
+    // Intermediate cylinder
+    arc 11 16 origin (0 0 $zmin)
+    arc 16 10 origin (0 0 $zmin)
+    arc 12 15 origin (0 0 $zmin)
+    arc 15  9 origin (0 0 $zmin)
+    arc 30 35 origin (0 0 $zmax)
+    arc 35 29 origin (0 0 $zmax)
+    arc 31 34 origin (0 0 $zmax)
+    arc 34 28 origin (0 0 $zmax)
 );
 
 defaultPatch
diff --git a/tutorials/mesh/blockMesh/sphere/system/blockMeshDict b/tutorials/mesh/blockMesh/sphere/system/blockMeshDict
index 59d5a948f8..191c1f6a49 100644
--- a/tutorials/mesh/blockMesh/sphere/system/blockMeshDict
+++ b/tutorials/mesh/blockMesh/sphere/system/blockMeshDict
@@ -16,65 +16,71 @@ FoamFile
 
 scale   1;
 
+// Geometric parameters
+outerRadius 1;
+
+// Divisions in x/y/z directions. Can be unequal.
+nx   10;
+ny   $nx;
+nz   $nx;
+
 geometry
 {
     sphere
     {
         type   sphere;
         origin (0 0 0);
-        radius 1;
+        radius $outerRadius;
     }
 }
 
-v 0.5773502;
-mv -0.5773502;
-
-a 0.7071067;
-ma -0.7071067;
+// Box sizes
+vo   #eval{ sqrt($outerRadius/3) };
+mvo  #eval{ -$vo };
 
 vertices
 (
-    ($mv $mv $mv)
-    ( $v $mv $mv)
-    ( $v  $v $mv)
-    ($mv  $v $mv)
-    ($mv $mv  $v)
-    ( $v $mv  $v)
-    ( $v  $v  $v)
-    ($mv  $v  $v)
+    ($mvo $mvo $mvo)
+    ( $vo $mvo $mvo)
+    ( $vo  $vo $mvo)
+    ($mvo  $vo $mvo)
+    ($mvo $mvo  $vo)
+    ( $vo $mvo  $vo)
+    ( $vo  $vo  $vo)
+    ($mvo  $vo  $vo)
 );
 
 blocks
 (
-    hex (0 1 2 3 4 5 6 7) (10 10 10) simpleGrading (1 1 1)
+    hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) grading (1 1 1)
 );
 
 edges
 (
-    arc 0 1 (0 $ma $ma)
-    arc 2 3 (0 $a $ma)
-    arc 6 7 (0 $a $a)
-    arc 4 5 (0 $ma $a)
+    arc 0 1 origin (0 0 0)
+    arc 2 3 origin (0 0 0)
+    arc 6 7 origin (0 0 0)
+    arc 4 5 origin (0 0 0)
 
-    arc 0 3 ($ma 0 $ma)
-    arc 1 2 ($a 0 $ma)
-    arc 5 6 ($a 0 $a)
-    arc 4 7 ($ma 0 $a)
+    arc 0 3 origin (0 0 0)
+    arc 1 2 origin (0 0 0)
+    arc 5 6 origin (0 0 0)
+    arc 4 7 origin (0 0 0)
 
-    arc 0 4 ($ma $ma 0)
-    arc 1 5 ($a $ma 0)
-    arc 2 6 ($a $a 0)
-    arc 3 7 ($ma $a 0)
+    arc 0 4 origin (0 0 0)
+    arc 1 5 origin (0 0 0)
+    arc 2 6 origin (0 0 0)
+    arc 3 7 origin (0 0 0)
 );
 
 faces
 (
-    project (0 4 7 3) sphere
-    project (2 6 5 1) sphere
-    project (1 5 4 0) sphere
-    project (3 7 6 2) sphere
-    project (0 3 2 1) sphere
-    project (4 5 6 7) sphere
+    project (0 4 7 3) sphere  // x-min
+    project (1 2 6 5) sphere  // x-max
+    project (0 1 5 4) sphere  // y-min
+    project (3 7 6 2) sphere  // y-max
+    project (0 3 2 1) sphere  // z-min
+    project (4 5 6 7) sphere  // z-max
 );
 
 boundary
@@ -84,12 +90,12 @@ boundary
         type wall;
         faces
         (
-            (0 4 7 3)
-            (2 6 5 1)
-            (1 5 4 0)
-            (3 7 6 2)
-            (0 3 2 1)
-            (4 5 6 7)
+            (0 4 7 3)  // x-min
+            (1 2 6 5)  // x-max
+            (0 1 5 4)  // y-min
+            (3 7 6 2)  // y-max
+            (0 3 2 1)  // z-min
+            (4 5 6 7)  // z-max
         );
     }
 );
diff --git a/tutorials/mesh/blockMesh/sphere7/system/blockMeshDict b/tutorials/mesh/blockMesh/sphere7/system/blockMeshDict
index 3e7c71adee..6d19d4a536 100644
--- a/tutorials/mesh/blockMesh/sphere7/system/blockMeshDict
+++ b/tutorials/mesh/blockMesh/sphere7/system/blockMeshDict
@@ -17,103 +17,115 @@ FoamFile
 scale   1;
 verbose no;
 
+// Geometric parameters
+outerRadius 1;
+innerRatio  0.7;
+
+// Divisions in x/y/z and radial directions. Can be unequal.
+nx   10;
+ny   $nx;
+nz   $nx;
+nr   6;
+
 geometry
 {
     sphere
     {
         type   sphere;
         origin (0 0 0);
-        radius 1;
+        radius $outerRadius;
     }
 }
 
+// Outer box sizes
+vo   #eval{ sqrt($outerRadius/3) };
+mvo  #eval{ -$vo };
 
-n    10;
-
-v    0.5773502;
-mv  -0.5773502;
-vh   0.2886751;
-mvh -0.2886751;
-
-a    0.7071067;
-ma  -0.7071067;
-ah   0.3464101;
-mah -0.3464101;
+// Inner box sizes - % of overall dimension
+vi   #eval{ $innerRatio*$vo };
+mvi  #eval{ $innerRatio*$mvo };
 
 vertices
 (
-    ($mvh $mvh $mvh)
-    ( $vh $mvh $mvh)
-    ( $vh  $vh $mvh)
-    ($mvh  $vh $mvh)
-    ($mvh $mvh  $vh)
-    ( $vh $mvh  $vh)
-    ( $vh  $vh  $vh)
-    ($mvh  $vh  $vh)
+    // Inner block
+    ($mvi $mvi $mvi)
+    ( $vi $mvi $mvi)
+    ( $vi  $vi $mvi)
+    ($mvi  $vi $mvi)
+    ($mvi $mvi  $vi)
+    ( $vi $mvi  $vi)
+    ( $vi  $vi  $vi)
+    ($mvi  $vi  $vi)
 
-    ($mv $mv $mv)
-    ( $v $mv $mv)
-    ( $v  $v $mv)
-    ($mv  $v $mv)
-    ($mv $mv  $v)
-    ( $v $mv  $v)
-    ( $v  $v  $v)
-    ($mv  $v  $v)
+    // Outer blocks
+    ($mvo $mvo $mvo)
+    ( $vo $mvo $mvo)
+    ( $vo  $vo $mvo)
+    ($mvo  $vo $mvo)
+    ($mvo $mvo  $vo)
+    ( $vo $mvo  $vo)
+    ( $vo  $vo  $vo)
+    ($mvo  $vo  $vo)
 );
 
 blocks
 (
-    hex ( 0  1  2  3  4  5  6  7) ($n $n $n) simpleGrading (1 1 1)
-    hex ( 9  8 12 13  1  0  4  5) ($n $n $n) simpleGrading (1 1 1)
-    hex (10  9 13 14  2  1  5  6) ($n $n $n) simpleGrading (1 1 1)
-    hex (11 10 14 15  3  2  6  7) ($n $n $n) simpleGrading (1 1 1)
-    hex ( 8 11 15 12  0  3  7  4) ($n $n $n) simpleGrading (1 1 1)
-    hex ( 8  9 10 11  0  1  2  3) ($n $n $n) simpleGrading (1 1 1)
-    hex (13 12 15 14  5  4  7  6) ($n $n $n) simpleGrading (1 1 1)
+    hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) grading (1 1 1)  // Inner block
+
+    // Outer blocks
+    hex ( 8  0  3 11 12  4  7 15) ($nr $ny $nz) grading (1 1 1)  // x-min
+    hex ( 1  9 10  2  5 13 14  6) ($nr $ny $nz) grading (1 1 1)  // x-max
+    hex ( 8  9  1  0 12 13  5  4) ($nx $nr $nz) grading (1 1 1)  // y-min
+    hex ( 3  2 10 11  7  6 14 15) ($nx $nr $nz) grading (1 1 1)  // y-max
+    hex ( 8  9 10 11  0  1  2  3) ($nx $ny $nr) grading (1 1 1)  // z-min
+    hex ( 4  5  6  7 12 13 14 15) ($nx $ny $nr) grading (1 1 1)  // z-max
 );
 
+
 edges
 (
-    arc  8  9 (0 $ma $ma)
-    arc 10 11 (0 $a $ma)
-    arc 14 15 (0 $a $a)
-    arc 12 13 (0 $ma $a)
+    // Outer blocks
+    arc  8  9 origin (0 0 0)
+    arc 10 11 origin (0 0 0)
+    arc 14 15 origin (0 0 0)
+    arc 12 13 origin (0 0 0)
 
-    arc  8 11 ($ma 0 $ma)
-    arc  9 10 ($a 0 $ma)
-    arc 13 14 ($a 0 $a)
-    arc 12 15 ($ma 0 $a)
+    arc  8 11 origin (0 0 0)
+    arc  9 10 origin (0 0 0)
+    arc 13 14 origin (0 0 0)
+    arc 12 15 origin (0 0 0)
 
-    arc  8 12 ($ma $ma 0)
-    arc  9 13 ($a $ma 0)
-    arc 10 14 ($a $a 0)
-    arc 11 15 ($ma $a 0)
+    arc  8 12 origin (0 0 0)
+    arc  9 13 origin (0 0 0)
+    arc 10 14 origin (0 0 0)
+    arc 11 15 origin (0 0 0)
 
+    // Inner block - flattened by factor 1.1
+    arc 0 1 origin 1.1 (0 0 0)
+    arc 2 3 origin 1.1 (0 0 0)
+    arc 6 7 origin 1.1 (0 0 0)
+    arc 4 5 origin 1.1 (0 0 0)
 
-    arc  0  1 (0 $mah $mah)
-    arc  2  3 (0 $ah $mah)
-    arc  6  7 (0 $ah $ah)
-    arc  4  5 (0 $mah $ah)
+    arc 0 3 origin 1.1 (0 0 0)
+    arc 1 2 origin 1.1 (0 0 0)
+    arc 5 6 origin 1.1 (0 0 0)
+    arc 4 7 origin 1.1 (0 0 0)
 
-    arc  0  3 ($mah 0 $mah)
-    arc  1  2 ($ah 0 $mah)
-    arc  5  6 ($ah 0 $ah)
-    arc  4  7 ($mah 0 $ah)
-
-    arc  0  4 ($mah $mah 0)
-    arc  1  5 ($ah $mah 0)
-    arc  2  6 ($ah $ah 0)
-    arc  3  7 ($mah $ah 0)
+    arc 0 4 origin 1.1 (0 0 0)
+    arc 1 5 origin 1.1 (0 0 0)
+    arc 2 6 origin 1.1 (0 0 0)
+    arc 3 7 origin 1.1 (0 0 0)
 );
 
 faces
 (
-    project ( 8 12 15 11) sphere
-    project (10 14 13  9) sphere
-    project ( 9 13 12  8) sphere
-    project (11 15 14 10) sphere
-    project ( 8 11 10  9) sphere
-    project (12 13 14 15) sphere
+    // Outer blocks
+    project ( 8 12 15 11) sphere  // x-min
+    project ( 9 10 14 13) sphere  // x-max
+    project ( 8  9 13 12) sphere  // y-min
+    project (11 15 14 10) sphere  // y-max
+    project ( 8 11 10  9) sphere  // z-min
+    project (12 13 14 15) sphere  // z-max
 );
 
 boundary
@@ -123,12 +135,12 @@ boundary
         type wall;
         faces
         (
-            ( 8 12 15 11)
-            (10 14 13  9)
-            ( 9 13 12  8)
-            (11 15 14 10)
-            ( 8 11 10  9)
-            (12 13 14 15)
+            ( 8 12 15 11)  // x-min
+            ( 9 10 14 13)  // x-max
+            ( 8  9 13 12)  // y-min
+            (11 15 14 10)  // y-max
+            ( 8 11 10  9)  // z-min
+            (12 13 14 15)  // z-max
         );
     }
 );
diff --git a/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/blockMeshDict b/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/blockMeshDict
index 27624af930..6f7108118f 100644
--- a/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/blockMeshDict
+++ b/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/blockMeshDict
@@ -17,68 +17,81 @@ FoamFile
 scale   1;
 verbose no;
 
+// Geometric parameters
+outerRadius 1;
+innerRatio  0.5;
+
+// Divisions in x/y/z and radial directions. Can be unequal.
+nx   10;
+ny   $nx;
+nz   $nx;
+nr   6;
+
 geometry
 {
     sphere
     {
         type   sphere;
         origin (0 0 0);
-        radius 1;
+        radius $outerRadius;
     }
 
     innerSphere
     {
         type   sphere;
         origin (0 0 0);
-        radius 0.5;
+        radius #eval{ $innerRatio * $outerRadius };
     }
 }
 
-n    10;
 
-// Rough approximation of corner points
-v    0.5773502;
-mv  -0.5773502;
-vh   0.2886751;
-mvh -0.2886751;
+// Outer box sizes (approximate)
+vo   #eval{ sqrt($outerRadius/3) };
+mvo  #eval{ -$vo };
+
+// Inner box sizes - % of overall dimension
+vi   #eval{ $innerRatio*$vo };
+mvi  #eval{ $innerRatio*$mvo };
 
 vertices
 (
     // Inner block points
-    project ($mvh $mvh $mvh) (innerSphere)
-    project ( $vh $mvh $mvh) (innerSphere)
-    project ( $vh  $vh $mvh) (innerSphere)
-    project ($mvh  $vh $mvh) (innerSphere)
-    project ($mvh $mvh  $vh) (innerSphere)
-    project ( $vh $mvh  $vh) (innerSphere)
-    project ( $vh  $vh  $vh) (innerSphere)
-    project ($mvh  $vh  $vh) (innerSphere)
+    project ($mvi $mvi $mvi) (innerSphere)
+    project ( $vi $mvi $mvi) (innerSphere)
+    project ( $vi  $vi $mvi) (innerSphere)
+    project ($mvi  $vi $mvi) (innerSphere)
+    project ($mvi $mvi  $vi) (innerSphere)
+    project ( $vi $mvi  $vi) (innerSphere)
+    project ( $vi  $vi  $vi) (innerSphere)
+    project ($mvi  $vi  $vi) (innerSphere)
 
     // Outer block points
-    project ($mv $mv $mv) (sphere)
-    project ( $v $mv $mv) (sphere)
-    project ( $v  $v $mv) (sphere)
-    project ($mv  $v $mv) (sphere)
-    project ($mv $mv  $v) (sphere)
-    project ( $v $mv  $v) (sphere)
-    project ( $v  $v  $v) (sphere)
-    project ($mv  $v  $v) (sphere)
+    project ($mvo $mvo $mvo) (sphere)
+    project ( $vo $mvo $mvo) (sphere)
+    project ( $vo  $vo $mvo) (sphere)
+    project ($mvo  $vo $mvo) (sphere)
+    project ($mvo $mvo  $vo) (sphere)
+    project ( $vo $mvo  $vo) (sphere)
+    project ( $vo  $vo  $vo) (sphere)
+    project ($mvo  $vo  $vo) (sphere)
 );
 
 blocks
 (
-    hex ( 0  1  2  3  4  5  6  7) ($n $n $n) simpleGrading (1 1 1)
-    hex ( 9  8 12 13  1  0  4  5) ($n $n $n) simpleGrading (1 1 1)
-    hex (10  9 13 14  2  1  5  6) ($n $n $n) simpleGrading (1 1 1)
-    hex (11 10 14 15  3  2  6  7) ($n $n $n) simpleGrading (1 1 1)
-    hex ( 8 11 15 12  0  3  7  4) ($n $n $n) simpleGrading (1 1 1)
-    hex ( 8  9 10 11  0  1  2  3) ($n $n $n) simpleGrading (1 1 1)
-    hex (13 12 15 14  5  4  7  6) ($n $n $n) simpleGrading (1 1 1)
+    hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) grading (1 1 1)  // Inner block
+
+    // Outer blocks
+    hex ( 8  0  3 11 12  4  7 15) ($nr $ny $nz) grading (1 1 1)  // x-min
+    hex ( 1  9 10  2  5 13 14  6) ($nr $ny $nz) grading (1 1 1)  // x-max
+    hex ( 8  9  1  0 12 13  5  4) ($nx $nr $nz) grading (1 1 1)  // y-min
+    hex ( 3  2 10 11  7  6 14 15) ($nx $nr $nz) grading (1 1 1)  // y-max
+    hex ( 8  9 10 11  0  1  2  3) ($nx $ny $nr) grading (1 1 1)  // z-min
+    hex ( 4  5  6  7 12 13 14 15) ($nx $ny $nr) grading (1 1 1)  // z-max
 );
 
 edges
 (
-    // Outer edges
+    // Outer blocks
     project  8  9 (sphere)
     project 10 11 (sphere)
     project 14 15 (sphere)
@@ -94,32 +107,32 @@ edges
     project 10 14 (sphere)
     project 11 15 (sphere)
 
+    // Inner block
+    project 0 1 (innerSphere)
+    project 2 3 (innerSphere)
+    project 6 7 (innerSphere)
+    project 4 5 (innerSphere)
 
-    // Inner edges
-    project  0  1 (innerSphere)
-    project  2  3 (innerSphere)
-    project  6  7 (innerSphere)
-    project  4  5 (innerSphere)
+    project 0 3 (innerSphere)
+    project 1 2 (innerSphere)
+    project 5 6 (innerSphere)
+    project 4 7 (innerSphere)
 
-    project  0  3 (innerSphere)
-    project  1  2 (innerSphere)
-    project  5  6 (innerSphere)
-    project  4  7 (innerSphere)
-
-    project  0  4 (innerSphere)
-    project  1  5 (innerSphere)
-    project  2  6 (innerSphere)
-    project  3  7 (innerSphere)
+    project 0 4 (innerSphere)
+    project 1 5 (innerSphere)
+    project 2 6 (innerSphere)
+    project 3 7 (innerSphere)
 );
 
 faces
 (
-    project ( 8 12 15 11) sphere
-    project (10 14 13  9) sphere
-    project ( 9 13 12  8) sphere
-    project (11 15 14 10) sphere
-    project ( 8 11 10  9) sphere
-    project (12 13 14 15) sphere
+    // Outer blocks
+    project ( 8 12 15 11) sphere  // x-min
+    project ( 9 10 14 13) sphere  // x-max
+    project ( 8  9 13 12) sphere  // y-min
+    project (11 15 14 10) sphere  // y-max
+    project ( 8 11 10  9) sphere  // z-min
+    project (12 13 14 15) sphere  // z-max
 );
 
 boundary
@@ -129,12 +142,12 @@ boundary
         type wall;
         faces
         (
-            ( 8 12 15 11)
-            (10 14 13  9)
-            ( 9 13 12  8)
-            (11 15 14 10)
-            ( 8 11 10  9)
-            (12 13 14 15)
+            ( 8 12 15 11)  // x-min
+            ( 9 10 14 13)  // x-max
+            ( 8  9 13 12)  // y-min
+            (11 15 14 10)  // y-max
+            ( 8 11 10  9)  // z-min
+            (12 13 14 15)  // z-max
         );
     }
 );
diff --git a/tutorials/mesh/refineMesh/refineFieldDirs/system/blockMeshDict b/tutorials/mesh/refineMesh/refineFieldDirs/system/blockMeshDict
index a7266f14e1..f97a1d3387 100644
--- a/tutorials/mesh/refineMesh/refineFieldDirs/system/blockMeshDict
+++ b/tutorials/mesh/refineMesh/refineFieldDirs/system/blockMeshDict
@@ -1,9 +1,15 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
 FoamFile
 {
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
     object      blockMeshDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -25,13 +31,10 @@ scale   1;
                          D
 */
 
-Ax  0.438912;    Ay  0.000000;
-Bx  18.28800;    By  0.000000;
-Cx  0.310358;    Cy  -0.310358;
-Dx  12.931569;   Dy  -12.931569;
-
-Hx  0.405502;    Hy  -0.167964;
-Gx  16.895909;   Gy  -6.998515;
+Ax  0.438912;    Ay  0;
+Bx  18.28800;    By  0;
+Cx  0.310358;    Cy  #eval{-$Cx};
+Dx  12.931569;   Dy  #eval{-$Dx};
 
 Z_DB_low  -0.88706;
 Z_AC_low  -1.07;
@@ -51,16 +54,15 @@ vertices
 
 blocks
 (
-   hex (0 1 2 3 4 5 6 7)    (47 10 4) simpleGrading (41.6669 1 1)
+    hex (0 1 2 3 4 5 6 7) (47 10 4) simpleGrading (41.6669 1 1)
 );
 
 edges
 (
-  arc 0 3 ($Hx $Hy $Z_AC_low)
-  arc 4 7 ($Hx $Hy $Z_high)
-
-  arc 1 2 ($Gx $Gy $Z_DB_low)
-  arc 5 6 ($Gx $Gy $Z_high)
+    arc 0 3 origin (0 0 $Z_AC_low)
+    arc 4 7 origin (0 0 $Z_high)
+    arc 1 2 origin (0 0 $Z_DB_low)
+    arc 5 6 origin (0 0 $Z_high)
 );
 
 patches
diff --git a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/system/blockMeshDict b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/system/blockMeshDict
index bce421f9a9..b194cb8d50 100644
--- a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/system/blockMeshDict
+++ b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/system/blockMeshDict
@@ -53,14 +53,14 @@ blocks
 
 edges
 (
-    arc 0 5 (0.469846 0.17101 0)
-    arc 5 10 (0.17101 0.469846 0)
-    arc 1 4 (0.939693 0.34202 0)
-    arc 4 9 (0.34202 0.939693 0)
-    arc 11 16 (0.469846 0.17101 0.5)
-    arc 16 21 (0.17101 0.469846 0.5)
-    arc 12 15 (0.939693 0.34202 0.5)
-    arc 15 20 (0.34202 0.939693 0.5)
+    arc  0  5 origin (0 0 0)
+    arc  5 10 origin (0 0 0)
+    arc  1  4 origin (0 0 0)
+    arc  4  9 origin (0 0 0)
+    arc 11 16 origin (0 0 0.5)
+    arc 16 21 origin (0 0 0.5)
+    arc 12 15 origin (0 0 0.5)
+    arc 15 20 origin (0 0 0.5)
 );
 
 boundary