diff --git a/applications/utilities/mesh/generation/blockMesh/blockMeshOBJ.H b/applications/utilities/mesh/generation/blockMesh/blockMeshOBJ.H
index 83325bef651bae58d9ff82c077a560d86f95afb1..726ec450e8db5ef1d27f37edbf453b9299887971 100644
--- a/applications/utilities/mesh/generation/blockMesh/blockMeshOBJ.H
+++ b/applications/utilities/mesh/generation/blockMesh/blockMeshOBJ.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -16,7 +16,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    const polyMesh& topoMesh = blocks.topology();
+    refPtr<polyMesh> topoMeshPtr(blocks.topology(true));
+    const polyMesh& topoMesh = topoMeshPtr();
 
     // Write mesh as edges
     {
diff --git a/applications/utilities/mesh/generation/blockMesh/blockMeshVTK.H b/applications/utilities/mesh/generation/blockMesh/blockMeshVTK.H
index 431b0b357eed41a7536c1aa2363403c4974d8fc1..55c2554dece3f73a440b07f3931ec904057895d1 100644
--- a/applications/utilities/mesh/generation/blockMesh/blockMeshVTK.H
+++ b/applications/utilities/mesh/generation/blockMesh/blockMeshVTK.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -19,7 +19,9 @@ Description
     // Non-legacy and ASCII (mesh is small, want readable output)
     const vtk::outputOptions writeOpts = vtk::formatType::INLINE_ASCII;
 
-    const polyMesh& topoMesh = blocks.topology();
+    refPtr<polyMesh> topoMeshPtr(blocks.topology(true));
+    const polyMesh& topoMesh = topoMeshPtr();
+
     const vtk::vtuCells topoCells(topoMesh, writeOpts);
 
     vtk::internalMeshWriter writer
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
index 327d402c1c5d14cd71d24eb699074201e78fd7a0..acc1bf549be9a5b0d4e62b05033d33dd30e8a470 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,6 +31,67 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+bool Foam::blockDescriptor::assignGradings
+(
+    const UList<gradingDescriptors>& ratios
+)
+{
+    bool ok = true;
+
+    switch (ratios.size())
+    {
+        case 0:
+        {
+            expand_.resize(12);
+            expand_ = gradingDescriptors();
+            break;
+        }
+        case 1:
+        {
+            // Identical in x/y/z-directions
+            expand_.resize(12);
+            expand_ = ratios[0];
+            break;
+        }
+        case 3:
+        {
+            expand_.resize(12);
+
+            // x-direction
+            expand_[0]  = ratios[0];
+            expand_[1]  = ratios[0];
+            expand_[2]  = ratios[0];
+            expand_[3]  = ratios[0];
+
+            // y-direction
+            expand_[4]  = ratios[1];
+            expand_[5]  = ratios[1];
+            expand_[6]  = ratios[1];
+            expand_[7]  = ratios[1];
+
+            // z-direction
+            expand_[8]  = ratios[2];
+            expand_[9]  = ratios[2];
+            expand_[10] = ratios[2];
+            expand_[11] = ratios[2];
+            break;
+        }
+        case 12:
+        {
+            expand_ = ratios;
+            break;
+        }
+        default:
+        {
+            ok = false;
+            break;
+        }
+    }
+
+    return ok;
+}
+
+
 void Foam::blockDescriptor::check(const Istream& is)
 {
     for (const label pointi : blockShape_)
@@ -38,8 +99,8 @@ void Foam::blockDescriptor::check(const Istream& is)
         if (pointi < 0 || pointi >= vertices_.size())
         {
             FatalIOErrorInFunction(is)
-                << "Point label " << pointi
-                << " out of range 0.." << vertices_.size() - 1
+                << "Point label (" << pointi
+                << ") out of range 0.." << vertices_.size() - 1
                 << " in block " << *this
                 << exit(FatalIOError);
         }
@@ -100,7 +161,7 @@ void Foam::blockDescriptor::check(const Istream& is)
 }
 
 
-void Foam::blockDescriptor::findCurvedFaces()
+void Foam::blockDescriptor::findCurvedFaces(const label blockIndex)
 {
     const faceList shapeFaces(blockShape().faces());
 
@@ -108,17 +169,21 @@ void Foam::blockDescriptor::findCurvedFaces()
     {
         forAll(blockFaces_, facei)
         {
+            const face& f = blockFaces_[facei].vertices();
+
+            // Accept (<block> <face>) face description
             if
             (
-                face::sameVertices
                 (
-                    blockFaces_[facei].vertices(),
-                    shapeFaces[shapeFacei]
+                    f.size() == 2
+                 && f[0] == blockIndex
+                 && f[1] == shapeFacei
                 )
+             || face::sameVertices(f, shapeFaces[shapeFacei])
             )
             {
                 curvedFaces_[shapeFacei] = facei;
-                nCurvedFaces_++;
+                ++nCurvedFaces_;
                 break;
             }
         }
@@ -144,19 +209,15 @@ Foam::blockDescriptor::blockDescriptor
     blockEdges_(edges),
     blockFaces_(faces),
     blockShape_(bshape),
-    expand_(expand),
+    expand_(),
     zoneName_(zoneName),
     curvedFaces_(-1),
     nCurvedFaces_(0)
 {
-    if (expand_.empty())
-    {
-        expand_.resize(12, gradingDescriptors());
-    }
-    else if (expand_.size() != 12)
+    if (!assignGradings(expand))
     {
         FatalErrorInFunction
-            << "Unknown definition of expansion ratios"
+            << "Unknown definition of expansion ratios: " << expand
             << exit(FatalError);
     }
 
@@ -167,7 +228,7 @@ Foam::blockDescriptor::blockDescriptor
 Foam::blockDescriptor::blockDescriptor
 (
     const dictionary& dict,
-    const label index,
+    const label blockIndex,
     const pointField& vertices,
     const blockEdgeList& edges,
     const blockFaceList& faces,
@@ -178,7 +239,8 @@ Foam::blockDescriptor::blockDescriptor
     vertices_(vertices),
     blockEdges_(edges),
     blockFaces_(faces),
-    expand_(12, gradingDescriptors()),
+    blockShape_(),
+    expand_(),
     zoneName_(),
     curvedFaces_(-1),
     nCurvedFaces_(0)
@@ -218,7 +280,7 @@ Foam::blockDescriptor::blockDescriptor
         else
         {
             FatalIOErrorInFunction(is)
-                << "incorrect token while reading n, expected '(', found "
+                << "Incorrect token while reading n, expected '(', found "
                 << t.info()
                 << exit(FatalIOError);
         }
@@ -226,6 +288,10 @@ Foam::blockDescriptor::blockDescriptor
     else
     {
         // Old-style: read three labels
+        IOWarningInFunction(is)
+            << "Encountered old-style specification of mesh divisions"
+            << endl;
+
         is  >> ijkMesh::sizes().x()
             >> ijkMesh::sizes().y()
             >> ijkMesh::sizes().z();
@@ -237,47 +303,18 @@ Foam::blockDescriptor::blockDescriptor
         is.putBack(t);
     }
 
-    List<gradingDescriptors> expRatios(is);
+    List<gradingDescriptors> expand(is);
 
-    if (expRatios.size() == 1)
-    {
-        // Identical in x/y/z-directions
-        expand_ = expRatios[0];
-    }
-    else if (expRatios.size() == 3)
+    if (!assignGradings(expand))
     {
-        // x-direction
-        expand_[0]  = expRatios[0];
-        expand_[1]  = expRatios[0];
-        expand_[2]  = expRatios[0];
-        expand_[3]  = expRatios[0];
-
-        // y-direction
-        expand_[4]  = expRatios[1];
-        expand_[5]  = expRatios[1];
-        expand_[6]  = expRatios[1];
-        expand_[7]  = expRatios[1];
-
-        // z-direction
-        expand_[8]  = expRatios[2];
-        expand_[9]  = expRatios[2];
-        expand_[10] = expRatios[2];
-        expand_[11] = expRatios[2];
-    }
-    else if (expRatios.size() == 12)
-    {
-        expand_ = expRatios;
-    }
-    else
-    {
-        FatalIOErrorInFunction(is)
-            << "Unknown definition of expansion ratios: " << expRatios
-            << exit(FatalIOError);
+        FatalErrorInFunction
+            << "Unknown definition of expansion ratios: " << expand
+            << exit(FatalError);
     }
 
     check(is);
 
-    findCurvedFaces();
+    findCurvedFaces(blockIndex);
 }
 
 
@@ -403,7 +440,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const blockDescriptor& bd)
     }
 
     os  << ' '  << bd.density()
-        << " simpleGrading (";
+        << " grading (";
 
 
     const List<gradingDescriptors>& expand = bd.grading();
@@ -445,8 +482,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const blockDescriptor& bd)
         }
     }
 
-
-    os  << ")";
+    os  << ')';
 
     return os;
 }
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
index b168363fb80d9391fbacd49e4756697a47d5dfb6..4af6734cd3c6ead18d4977409b4fca54b7296b94 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
@@ -109,22 +109,26 @@ class blockDescriptor
 
     // Private Member Functions
 
+        //- Assign edge grading.
+        //  \return false for unsupported specification
+        bool assignGradings(const UList<gradingDescriptors>& ratios);
+
         //- Check block has outward-pointing faces
         void check(const Istream& is);
 
-        //- Calculate the points and weights for the specified edge.
-        //  Return the number of curved edges
-        label edgePointsWeights
+        //- Calculate the points and weights for the specified hex edge.
+        //  Return the number of curved edges (0-1)
+        int calcEdgePointsWeights
         (
-            pointField (&edgePoints)[12],
-            scalarList (&edgeWeights)[12],
-            const label edgei,
+            pointField& edgePoints,
+            scalarList& edgeWeights,
             const label start,
             const label end,
-            const label dim
+            const label nDiv,
+            const gradingDescriptors& expand
         ) const;
 
-        void findCurvedFaces();
+        void findCurvedFaces(const label blockIndex = -1);
 
 public:
 
@@ -156,7 +160,7 @@ public:
         blockDescriptor
         (
             const dictionary& dict,
-            const label index,
+            const label blockIndex,
             const pointField& vertices,
             const blockEdgeList& edges,
             const blockFaceList& faces,
@@ -167,28 +171,28 @@ public:
     // Member Functions
 
         //- Reference to point field defining the block mesh
-        inline const pointField& vertices() const;
+        inline const pointField& vertices() const noexcept;
 
         //- Return reference to the list of curved faces
-        inline const blockFaceList& blockFaces() const;
+        inline const blockFaceList& blockFaces() const noexcept;
 
         //- Return the block shape
-        inline const cellShape& blockShape() const;
+        inline const cellShape& blockShape() const noexcept;
 
         //- Return the mesh density (number of cells) in the i,j,k directions
-        inline const labelVector& density() const;
+        inline const labelVector& density() const noexcept;
 
         //- Expansion ratios in all directions
-        const List<gradingDescriptors>& grading() const;
+        inline const List<gradingDescriptors>& grading() const noexcept;
 
         //- Return the (optional) zone name
-        inline const word& zoneName() const;
+        inline const word& zoneName() const noexcept;
 
         //- Curved-face labels for each block-face (-1 for flat faces)
-        inline const FixedList<label, 6>& curvedFaces() const;
+        inline const FixedList<label, 6>& curvedFaces() const noexcept;
 
         //- Number of curved faces in this block
-        inline label nCurvedFaces() const;
+        inline label nCurvedFaces() const noexcept;
 
         //- Return block point for local label i
         inline const point& blockPoint(const label i) const;
@@ -209,11 +213,32 @@ public:
         inline bool edge(const label i, const label j, const label k) const;
 
         //- Calculate the points and weights for all edges.
-        //  Return the number of curved edges
-        label edgesPointsWeights
+        //  \return the number of curved edges (0-12)
+        int edgesPointsWeights
+        (
+            pointField (&edgesPoints)[12],
+            scalarList (&edgesWeights)[12]
+        ) const;
+
+        //- Calculate points and weights for specified edge,
+        //- using the specified number of divisions and grading
+        //  \return True if the edge is curved
+        bool edgePointsWeights
         (
-            pointField (&edgePoints)[12],
-            scalarList (&edgeWeights)[12]
+            const label edgei,
+            pointField& edgePoints,
+            scalarList& edgeWeights,
+            const label nDiv,
+            const gradingDescriptors& gd = gradingDescriptors()
+        ) const;
+
+        //- Calculate points and weights for specified edge.
+        //  \return True if the edge is curved
+        bool edgePointsWeights
+        (
+            const label edgei,
+            pointField& edgePoints,
+            scalarList& edgeWeights
         ) const;
 
         //- Return true if point i,j,k addresses a block flat face or edge
@@ -235,7 +260,7 @@ public:
         static void write(Ostream&, const label blocki, const dictionary&);
 
         //- Return info proxy
-        inline InfoProxy<blockDescriptor> info() const
+        InfoProxy<blockDescriptor> info() const
         {
             return *this;
         }
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
index f3ab13497c63b978e2d1bd4a781aa3bad5398572..af615b13bbe5be8e7a7b94294d2a06a9797bb3f6 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,24 +30,29 @@ License
 #include "lineEdge.H"
 #include "lineDivide.H"
 
+// * * * * * * * * * * * * * * Local Data Members  * * * * * * * * * * * * * //
+
+// Warning.
+// Ordering of edges needs to be the same as hex cell shape model
+static const int hexEdge0[12] = { 0, 3, 7, 4,  0, 1, 5, 4,  0, 1, 2, 3 };
+static const int hexEdge1[12] = { 1, 2, 6, 5,  3, 2, 6, 7,  4, 5, 6, 7 };
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::label Foam::blockDescriptor::edgePointsWeights
+int Foam::blockDescriptor::calcEdgePointsWeights
 (
-    pointField (&edgePoints)[12],
-    scalarList (&edgeWeights)[12],
-    const label edgei,
+    pointField& edgePoints,
+    scalarList& edgeWeights,
     const label start,
     const label end,
-    const label nDiv
+    const label nDiv,
+    const gradingDescriptors& expand
 ) const
 {
     // Set reference to the list of labels defining the block
     const labelList& blockLabels = blockShape_;
 
-    // Get list of points for this block
-    const pointField blockPoints = blockShape_.points(vertices_);
-
     // Set the edge points/weights
     // The edge is a straight-line if it is not in the list of blockEdges
 
@@ -55,55 +60,57 @@ Foam::label Foam::blockDescriptor::edgePointsWeights
     {
         const int cmp = cedge.compare(blockLabels[start], blockLabels[end]);
 
-        if (cmp)
+        if (cmp > 0)
         {
-            if (cmp > 0)
-            {
-                // Curve has the same orientation
+            // Curve has the same orientation
 
-                // Divide the line
-                const lineDivide divEdge(cedge, nDiv, expand_[edgei]);
+            // Divide the line
+            const lineDivide divEdge(cedge, nDiv, expand);
 
-                edgePoints[edgei]  = divEdge.points();
-                edgeWeights[edgei] = divEdge.lambdaDivisions();
-            }
-            else
-            {
-                // Curve has the opposite orientation
+            edgePoints  = divEdge.points();
+            edgeWeights = divEdge.lambdaDivisions();
+
+            return 1;  // Found curved-edge: done
+        }
+        else if (cmp < 0)
+        {
+            // Curve has the opposite orientation
 
-                // Divide the line
-                const lineDivide divEdge(cedge, nDiv, expand_[edgei].inv());
+            // Divide the line
+            const lineDivide divEdge(cedge, nDiv, expand.inv());
 
-                const pointField& p = divEdge.points();
-                const scalarList& d = divEdge.lambdaDivisions();
+            const pointField& p = divEdge.points();
+            const scalarList& d = divEdge.lambdaDivisions();
 
-                edgePoints[edgei].setSize(p.size());
-                edgeWeights[edgei].setSize(d.size());
+            edgePoints.resize(p.size());
+            edgeWeights.resize(d.size());
 
-                label pn = p.size() - 1;
-                forAll(p, pi)
-                {
-                    edgePoints[edgei][pi]  = p[pn - pi];
-                    edgeWeights[edgei][pi] = 1 - d[pn - pi];
-                }
+            // Copy in reverse order
+            const label pn = (p.size() - 1);
+            forAll(p, pi)
+            {
+                edgePoints[pi]  = p[pn - pi];
+                edgeWeights[pi] = 1 - d[pn - pi];
             }
 
-            // Found curved-edge: done
-            return 1;
+            return 1;  // Found curved-edge: done
         }
     }
 
-
     // Not curved-edge: divide the edge as a straight line
+
+    // Get list of points for this block
+    const pointField blockPoints(blockShape_.points(vertices_));
+
     lineDivide divEdge
     (
         blockEdges::lineEdge(blockPoints, start, end),
         nDiv,
-        expand_[edgei]
+        expand
     );
 
-    edgePoints[edgei]  = divEdge.points();
-    edgeWeights[edgei] = divEdge.lambdaDivisions();
+    edgePoints = divEdge.points();
+    edgeWeights = divEdge.lambdaDivisions();
 
     return 0;
 }
@@ -111,37 +118,89 @@ Foam::label Foam::blockDescriptor::edgePointsWeights
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::label Foam::blockDescriptor::edgesPointsWeights
+int Foam::blockDescriptor::edgesPointsWeights
+(
+    pointField (&edgesPoints)[12],
+    scalarList (&edgesWeights)[12]
+) const
+{
+    int nCurved = 0;
+
+    for (label edgei = 0; edgei < 12; ++edgei)
+    {
+        nCurved += calcEdgePointsWeights
+        (
+            edgesPoints[edgei],
+            edgesWeights[edgei],
+            hexEdge0[edgei],
+            hexEdge1[edgei],
+
+            sizes()[edgei/4],   // 12 edges -> 3 components (x,y,z)
+            expand_[edgei]
+        );
+    }
+
+    return nCurved;
+}
+
+
+bool Foam::blockDescriptor::edgePointsWeights
 (
-    pointField (&edgePoints)[12],
-    scalarList (&edgeWeights)[12]
+    const label edgei,
+    pointField& edgePoints,
+    scalarList& edgeWeights,
+    const label nDiv,
+    const gradingDescriptors& gd
 ) const
 {
-    const label ni = sizes().x();
-    const label nj = sizes().y();
-    const label nk = sizes().z();
-
-    label nCurvedEdges = 0;
-
-    // X-direction
-    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 0,  0, 1, ni);
-    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 1,  3, 2, ni);
-    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 2,  7, 6, ni);
-    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 3,  4, 5, ni);
-
-    // Y-direction
-    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 4,  0, 3, nj);
-    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 5,  1, 2, nj);
-    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 6,  5, 6, nj);
-    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 7,  4, 7, nj);
-
-    // Z-direction
-    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 8,  0, 4, nk);
-    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 9,  1, 5, nk);
-    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 10, 2, 6, nk);
-    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 11, 3, 7, nk);
-
-    return nCurvedEdges;
+    if (edgei < 0 || edgei >= 12)
+    {
+        FatalErrorInFunction
+            << "Edge label " << edgei
+            << " out of range 0..11"
+            << exit(FatalError);
+    }
+
+    const int nCurved = calcEdgePointsWeights
+    (
+        edgePoints,
+        edgeWeights,
+        hexEdge0[edgei],
+        hexEdge1[edgei],
+        nDiv,
+        gd
+    );
+
+    return nCurved;
+}
+
+
+bool Foam::blockDescriptor::edgePointsWeights
+(
+    const label edgei,
+    pointField& edgePoints,
+    scalarList& edgeWeights
+) const
+{
+    if (edgei < 0 || edgei >= 12)
+    {
+        FatalErrorInFunction
+            << "Edge label " << edgei
+            << " out of range 0..11"
+            << exit(FatalError);
+    }
+
+    const int nCurved = calcEdgePointsWeights
+    (
+        edgePoints,
+        edgeWeights,
+        hexEdge0[edgei],
+        hexEdge1[edgei],
+        sizes()[edgei/4],   // 12 edges -> 3 components (x,y,z)
+        expand_[edgei]
+    );
+
+    return nCurved;
 }
 
 
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
index 88c3d66dcdd08788ce8ae6b5592c2889c4e4581f..cc8c83cb2236738d34a88812bd75014ea96f9580 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,51 +28,55 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const Foam::pointField& Foam::blockDescriptor::vertices() const
+inline const Foam::pointField&
+Foam::blockDescriptor::vertices() const noexcept
 {
     return vertices_;
 }
 
 
-inline const Foam::blockFaceList& Foam::blockDescriptor::blockFaces() const
+inline const Foam::blockFaceList&
+Foam::blockDescriptor::blockFaces() const noexcept
 {
     return blockFaces_;
 }
 
 
-inline const Foam::cellShape& Foam::blockDescriptor::blockShape() const
+inline const Foam::cellShape&
+Foam::blockDescriptor::blockShape() const noexcept
 {
     return blockShape_;
 }
 
 
-inline const Foam::labelVector& Foam::blockDescriptor::density() const
+inline const Foam::labelVector&
+Foam::blockDescriptor::density() const noexcept
 {
     return ijkMesh::sizes();
 }
 
 
 inline const Foam::List<Foam::gradingDescriptors>&
-Foam::blockDescriptor::grading() const
+Foam::blockDescriptor::grading() const noexcept
 {
     return expand_;
 }
 
 
-inline const Foam::word& Foam::blockDescriptor::zoneName() const
+inline const Foam::word& Foam::blockDescriptor::zoneName() const noexcept
 {
     return zoneName_;
 }
 
 
 inline const Foam::FixedList<Foam::label, 6>&
-Foam::blockDescriptor::curvedFaces() const
+Foam::blockDescriptor::curvedFaces() const noexcept
 {
     return curvedFaces_;
 }
 
 
-inline Foam::label Foam::blockDescriptor::nCurvedFaces() const
+inline Foam::label Foam::blockDescriptor::nCurvedFaces() const noexcept
 {
     return nCurvedFaces_;
 }
@@ -154,8 +158,10 @@ inline bool Foam::blockDescriptor::flatFaceOrEdge
 {
     if (i == 0 && curvedFaces_[0] == -1) return true;
     if (i == sizes().x() && curvedFaces_[1] == -1) return true;
+
     if (j == 0 && curvedFaces_[2] == -1) return true;
     if (j == sizes().y() && curvedFaces_[3] == -1) return true;
+
     if (k == 0 && curvedFaces_[4] == -1) return true;
     if (k == sizes().z() && curvedFaces_[5] == -1) return true;
 
diff --git a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H
index 850353274ab2842f6401d2e2c2e008613ce4a6b3..3e87fe570a51a6d10681bc0b6ddc47b9f6d106b9 100644
--- a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H
+++ b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H
@@ -45,7 +45,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
+// Forward Declarations
 class blockDescriptor;
 class blockFace;
 
@@ -91,7 +91,7 @@ public:
     // Constructors
 
         //- Construct from face vertices
-        blockFace(const face& vertices);
+        explicit blockFace(const face& vertices);
 
         //- Construct from Istream
         blockFace
@@ -163,7 +163,7 @@ public:
         void write(Ostream&, const dictionary&) const;
 
 
-    // Ostream operator
+    // Ostream Operator
 
         friend Ostream& operator<<(Ostream&, const blockFace&);
 };
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.C b/src/mesh/blockMesh/blockMesh/blockMesh.C
index ef0c90f1dac7efad967dc0dba1eb3e97bddf0046..18321776522b13b19363805140f8445918aa68d0 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.C
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,6 +27,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "blockMesh.H"
+#include "transform.H"
 #include "Time.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -47,6 +48,182 @@ Foam::blockMesh::strategyNames_
 });
 
 
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Process dictionary entry with single scalar or vector quantity
+// - return 0 if scaling is not needed. Eg, Not found or is unity
+// - return 1 for uniform scaling
+// - return 3 for non-uniform scaling
+
+int readScaling(const entry *eptr, vector& scale)
+{
+    int nCmpt = 0;
+
+    if (!eptr)
+    {
+        return nCmpt;
+    }
+
+    ITstream& is = eptr->stream();
+
+    if (is.peek().isNumber())
+    {
+        // scalar value
+        scalar val;
+        is >> val;
+
+        if ((val > 0) && !equal(val, 1))
+        {
+            // Uniform scaling
+            nCmpt = 1;
+            scale = vector::uniform(val);
+        }
+    }
+    else
+    {
+        // vector value
+        is >> scale;
+
+        bool nonUnity = false;
+        for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+        {
+            if (scale[cmpt] <= 0)
+            {
+                scale[cmpt] = 1;
+            }
+            else if (!equal(scale[cmpt], 1))
+            {
+                nonUnity = true;
+            }
+        }
+
+        if (nonUnity)
+        {
+            if (equal(scale.x(), scale.y()) && equal(scale.x(), scale.z()))
+            {
+                // Uniform scaling
+                nCmpt = 1;
+            }
+            else
+            {
+                nCmpt = 3;
+            }
+        }
+    }
+
+    eptr->checkITstream(is);
+
+    return nCmpt;
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::blockMesh::readPointTransforms(const dictionary& dict)
+{
+    transformType_ = transformTypes::NO_TRANSFORM;
+
+    // Optional cartesian coordinate system transform, since JUL-2021
+    if (const dictionary* dictptr = dict.findDict("transform"))
+    {
+        transform_ = coordSystem::cartesian(*dictptr);
+
+        constexpr scalar tol = VSMALL;
+
+        // Check for non-zero origin
+        const vector& o = transform_.origin();
+        if
+        (
+            (mag(o.x()) > tol)
+         || (mag(o.y()) > tol)
+         || (mag(o.z()) > tol)
+        )
+        {
+            transformType_ |= transformTypes::TRANSLATION;
+        }
+
+        // Check for non-identity rotation
+        const tensor& r = transform_.R();
+        if
+        (
+            (mag(r.xx() - 1) > tol)
+         || (mag(r.xy()) > tol)
+         || (mag(r.xz()) > tol)
+
+         || (mag(r.yx()) > tol)
+         || (mag(r.yy() - 1) > tol)
+         || (mag(r.yz()) > tol)
+
+         || (mag(r.zx()) > tol)
+         || (mag(r.zy()) > tol)
+         || (mag(r.zz() - 1) > tol)
+        )
+        {
+            transformType_ |= transformTypes::ROTATION;
+        }
+    }
+    else
+    {
+        transform_.clear();
+    }
+
+
+    // Optional 'prescale' factor.
+    {
+        prescaling_ = vector::uniform(1);
+
+        const int scaleType = readScaling
+        (
+            dict.findEntry("prescale", keyType::LITERAL),
+            prescaling_
+        );
+
+        if (scaleType == 1)
+        {
+            transformType_ |= transformTypes::PRESCALING;
+        }
+        else if (scaleType == 3)
+        {
+            transformType_ |= transformTypes::PRESCALING3;
+        }
+    }
+
+
+    // Optional 'scale' factor. Was 'convertToMeters' until OCT-2008
+    {
+        scaling_ = vector::uniform(1);
+
+        const int scaleType = readScaling
+        (
+            // Mark as changed from 2010 onwards
+            dict.findCompat
+            (
+                "scale",
+                {{"convertToMeters", 1012}},
+                keyType::LITERAL
+            ),
+            scaling_
+        );
+
+        if (scaleType == 1)
+        {
+            transformType_ |= transformTypes::SCALING;
+        }
+        else if (scaleType == 3)
+        {
+            transformType_ |= transformTypes::SCALING3;
+        }
+    }
+
+    return bool(transformType_);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::blockMesh::blockMesh
@@ -62,6 +239,8 @@ Foam::blockMesh::blockMesh
     (
         meshDict_.getOrDefault("checkFaceCorrespondence", true)
     ),
+    mergeStrategy_(strategy),
+    transformType_(transformTypes::NO_TRANSFORM),
     geometry_
     (
         IOobject
@@ -78,25 +257,27 @@ Foam::blockMesh::blockMesh
       : dictionary(),
         true
     ),
-    scaleFactor_(1),
     blockVertices_
     (
         meshDict_.lookup("vertices"),
         blockVertex::iNew(meshDict_, geometry_)
     ),
     vertices_(Foam::vertices(blockVertices_)),
+    prescaling_(vector::uniform(1)),
+    scaling_(vector::uniform(1)),
+    transform_(),
     topologyPtr_(createTopology(meshDict_, regionName))
 {
     // Command-line option has precedence over dictionary setting
 
-    if (strategy == mergeStrategy::DEFAULT_MERGE)
+    if (mergeStrategy_ == mergeStrategy::DEFAULT_MERGE)
     {
-        strategyNames_.readIfPresent("mergeType", meshDict_, strategy);
+        strategyNames_.readIfPresent("mergeType", meshDict_, mergeStrategy_);
 
         // Warn about fairly obscure old "fastMerge" option?
     }
 
-    if (strategy == mergeStrategy::MERGE_POINTS)
+    if (mergeStrategy_ == mergeStrategy::MERGE_POINTS)
     {
         // MERGE_POINTS
         calcGeometricalMerge();
@@ -131,35 +312,38 @@ bool Foam::blockMesh::verbose(const bool on) noexcept
 }
 
 
-const Foam::pointField& Foam::blockMesh::vertices() const
+const Foam::pointField& Foam::blockMesh::vertices() const noexcept
 {
     return vertices_;
 }
 
 
-const Foam::polyMesh& Foam::blockMesh::topology() const
+Foam::tmp<Foam::pointField>
+Foam::blockMesh::vertices(bool applyTransform) const
 {
-    if (!topologyPtr_)
+    if (applyTransform && hasPointTransforms())
     {
-        FatalErrorInFunction
-            << "topologyPtr_ not allocated"
-            << exit(FatalError);
+        auto tpts = tmp<pointField>::New(vertices_);
+
+        inplacePointTransforms(tpts.ref());
+
+        return tpts;
     }
 
-    return *topologyPtr_;
+    return vertices_;
 }
 
 
 Foam::PtrList<Foam::dictionary> Foam::blockMesh::patchDicts() const
 {
-    const polyPatchList& patchTopologies = topology().boundaryMesh();
+    const polyPatchList& topoPatches = topology().boundaryMesh();
 
-    PtrList<dictionary> patchDicts(patchTopologies.size());
+    PtrList<dictionary> patchDicts(topoPatches.size());
 
-    forAll(patchTopologies, patchi)
+    forAll(topoPatches, patchi)
     {
         OStringStream os;
-        patchTopologies[patchi].write(os);
+        topoPatches[patchi].write(os);
         IStringStream is(os.str());
         patchDicts.set(patchi, new dictionary(is));
     }
@@ -167,12 +351,6 @@ Foam::PtrList<Foam::dictionary> Foam::blockMesh::patchDicts() const
 }
 
 
-Foam::scalar Foam::blockMesh::scaleFactor() const
-{
-    return scaleFactor_;
-}
-
-
 const Foam::pointField& Foam::blockMesh::points() const
 {
     if (points_.empty())
@@ -242,4 +420,100 @@ Foam::label Foam::blockMesh::numZonedBlocks() const
 }
 
 
+bool Foam::blockMesh::hasPointTransforms() const noexcept
+{
+    return bool(transformType_);
+}
+
+
+bool Foam::blockMesh::inplacePointTransforms(pointField& pts) const
+{
+    if (!transformType_)
+    {
+        return false;
+    }
+
+    if (transformType_ & transformTypes::PRESCALING)
+    {
+        for (point& p : pts)
+        {
+            p *= prescaling_.x();
+        }
+    }
+    else if (transformType_ & transformTypes::PRESCALING3)
+    {
+        for (point& p : pts)
+        {
+            p = cmptMultiply(p, prescaling_);
+        }
+    }
+
+    if (transformType_ & transformTypes::ROTATION)
+    {
+        const tensor rot(transform_.R());
+
+        if (transformType_ & transformTypes::TRANSLATION)
+        {
+            const point origin(transform_.origin());
+
+            for (point& p : pts)
+            {
+                p = Foam::transform(rot, p) + origin;
+            }
+        }
+        else
+        {
+            for (point& p : pts)
+            {
+                p = Foam::transform(rot, p);
+            }
+        }
+    }
+    else if (transformType_ & transformTypes::TRANSLATION)
+    {
+        const point origin(transform_.origin());
+
+        for (point& p : pts)
+        {
+            p += origin;
+        }
+    }
+
+    if (transformType_ & transformTypes::SCALING)
+    {
+        for (point& p : pts)
+        {
+            p *= scaling_.x();
+        }
+    }
+    else if (transformType_ & transformTypes::SCALING3)
+    {
+        for (point& p : pts)
+        {
+            p = cmptMultiply(p, scaling_);
+        }
+    }
+
+    return true;
+}
+
+
+Foam::tmp<Foam::pointField>
+Foam::blockMesh::globalPosition(const pointField& localPoints) const
+{
+    if (hasPointTransforms())
+    {
+        auto tpts = tmp<pointField>::New(localPoints);
+
+        inplacePointTransforms(tpts.ref());
+
+        return tpts;
+    }
+    else
+    {
+        return localPoints;
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.H b/src/mesh/blockMesh/blockMesh/blockMesh.H
index ec55f4a20fe33c1c9f0532842aea214e0b2c62b0..df613d58063f7e963027ab9fd65b23d8ed57e04b 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.H
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,7 +33,9 @@ Description
     Dictionary controls
     \table
         Property    | Description                           | Required | Default
-        scale       | Point scaling                         | no  | 1.0
+        prescale    | Point scaling before transform        | no  | 1.0
+        scale       | Point scaling after transform         | no  | 1.0
+        transform   | Point transform (origin, rotation)    | no  |
         vertices    |                                       | yes |
         blocks      |                                       | yes |
         edges       |                                       | no  |
@@ -48,6 +50,9 @@ Description
     \endtable
 
 Note
+    The \c prescale and \c scale can be a single scalar or a vector of
+    values.
+
     The vertices, cells and patches for filling the blocks are demand-driven.
 
 SourceFiles
@@ -65,6 +70,7 @@ SourceFiles
 #include "Enum.H"
 #include "block.H"
 #include "PtrList.H"
+#include "cartesianCS.H"
 #include "searchableSurfaces.H"
 #include "polyMesh.H"
 #include "IOdictionary.H"
@@ -112,6 +118,21 @@ private:
         static const Enum<mergeStrategy> strategyNames_;
 
 
+    // Data Types
+
+        //- Point transformations. Internal usage
+        enum transformTypes : unsigned
+        {
+            NO_TRANSFORM = 0,   //!< No transformations
+            ROTATION = 0x1,     //!< Local coordinate system rotation
+            TRANSLATION = 0x2,  //!< Local coordinate system translation
+            PRESCALING = 0x4,   //!< Uniform scale before rotate/translate
+            PRESCALING3 = 0x8,  //!< Non-uniform scale before rotate/translate
+            SCALING = 0x10,     //!< Uniform scale after rotate/translate
+            SCALING3 = 0x20     //!< Non-uniform scale after rotate/translate
+        };
+
+
     // Private Data
 
         //- Reference to mesh dictionary
@@ -120,15 +141,18 @@ private:
         //- Output verbosity
         bool verbose_;
 
-        //- Switch checking face consistency (defaults to true)
+        //- Check face consistency (defaults to true)
         bool checkFaceCorrespondence_;
 
+        //- The mesh merging strategy
+        mergeStrategy mergeStrategy_;
+
+        //- Types of point transformations requested
+        unsigned transformType_;
+
         //- Optional searchable geometry to project face-points to
         searchableSurfaces geometry_;
 
-        //- The scaling factor to convert to metres
-        scalar scaleFactor_;
-
         //- The list of block vertices
         blockVertexList blockVertices_;
 
@@ -141,6 +165,15 @@ private:
         //- The list of curved faces
         blockFaceList faces_;
 
+        //- The scaling factor before rotate/translate
+        vector prescaling_;
+
+        //- The scaling factor after rotate/translate
+        vector scaling_;
+
+        //- Local coordinate system transformation
+        coordSystem::cartesian transform_;
+
         //- The blocks themselves (the topology) as a polyMesh
         autoPtr<polyMesh> topologyPtr_;
 
@@ -150,12 +183,12 @@ private:
         //- The sum of all cells in each block
         label nCells_;
 
-        //- The point offset added to each block
-        labelList blockOffsets_;
-
         //- The merge points information
         labelList mergeList_;
 
+        //- The point offset added to each block. Offsets into mergeList_
+        labelList blockOffsets_;
+
         mutable pointField points_;
 
         mutable cellShapeList cells_;
@@ -165,14 +198,9 @@ private:
 
     // Private Member Functions
 
-        template<class Source>
-        void checkPatchLabels
-        (
-            const Source& source,
-            const word& patchName,
-            const pointField& points,
-            faceList& patchShapes
-        ) const;
+        //- Get scaling and/or coordinate transforms
+        //  \return True if scaling and/or transformations are needed
+        bool readPointTransforms(const dictionary& dict);
 
         void readPatches
         (
@@ -191,7 +219,8 @@ private:
             PtrList<dictionary>& patchDicts
         );
 
-        void createCellShapes(cellShapeList& tmpBlockCells);
+        //- Topology blocks as cell shapes
+        cellShapeList getBlockShapes() const;
 
         autoPtr<polyMesh> createTopology
         (
@@ -253,63 +282,87 @@ public:
 
     // Member Functions
 
-        // Access
+    // General Access, Description
 
-            //- Access to input dictionary
-            const dictionary& meshDict() const
-            {
-                return meshDict_;
-            }
+        //- Access to input dictionary
+        const dictionary& meshDict() const noexcept
+        {
+            return meshDict_;
+        }
 
-            //- Optional searchable geometry to project face-points to
-            const searchableSurfaces& geometry() const
-            {
-                return geometry_;
-            }
+        //- Optional searchable geometry to project face-points to
+        const searchableSurfaces& geometry() const noexcept
+        {
+            return geometry_;
+        }
+
+        //- The curved edges
+        const blockEdgeList& edges() const noexcept
+        {
+            return edges_;
+        }
+
+        //- The curved faces
+        const blockFaceList& faces() const noexcept
+        {
+            return faces_;
+        }
+
+        //- True if the blockMesh topology exists
+        bool valid() const noexcept;
+
+        //- Return the patch names
+        wordList patchNames() const;
+
+        //- Number of blocks with specified zones
+        label numZonedBlocks() const;
 
-            //- True if the blockMesh topology exists
-            bool valid() const noexcept;
 
-            //- Reference to point field defining the blockMesh.
-            //  These points are \b not scaled by scaleFactor
-            const pointField& vertices() const;
+    // Point Transformations
 
-            //- Return the blockMesh topology as a polyMesh
-            const polyMesh& topology() const;
+        //- True if scaling and/or transformations are needed
+        bool hasPointTransforms() const noexcept;
 
-            //- Return the curved edges
-            const blockEdgeList& edges() const
-            {
-                return edges_;
-            }
+        //- Apply coordinate transforms and scaling
+        bool inplacePointTransforms(pointField& pts) const;
 
-            //- Return the curved faces
-            const blockFaceList& faces() const
-            {
-                return faces_;
-            }
+        //- Apply coordinate transforms and scaling
+        tmp<pointField> globalPosition(const pointField& localPoints) const;
 
-            //- The scaling factor used to convert to metres
-            scalar scaleFactor() const;
 
-            //- The points for the entire mesh.
-            //  These points \b are scaled by scaleFactor
-            const pointField& points() const;
+    // Block Topology
 
-            //- Return cell shapes list
-            const cellShapeList& cells() const;
+        //- Reference to point field defining the blocks,
+        //- these points are \b unscaled and \b non-transformed
+        const pointField& vertices() const noexcept;
 
-            //- Return the patch face lists
-            const faceListList& patches() const;
+        //- Point field defining the blocks,
+        //- optionally transformed and scaled
+        tmp<pointField> vertices(bool applyTransform) const;
 
-            //- Get patch information from the topology mesh
-            PtrList<dictionary> patchDicts() const;
+        //- The blockMesh topology as a polyMesh
+        //- \b unscaled and \b non-transformed
+        const polyMesh& topology() const;
 
-            //- Return patch names
-            wordList patchNames() const;
+        //- The blockMesh topology as a polyMesh
+        //- optionally transformed and scaled
+        refPtr<polyMesh> topology(bool applyTransform) const;
 
-            //- Number of blocks with specified zones
-            label numZonedBlocks() const;
+
+    // Detailed Mesh
+
+        //- The points for the entire mesh.
+        //- These points \b are scaled and transformed
+        const pointField& points() const;
+
+        //- Return cell shapes list
+        const cellShapeList& cells() const;
+
+        //- Return the patch face lists
+        const faceListList& patches() const;
+
+        //- Patch information from the topology mesh
+        PtrList<dictionary> patchDicts() const;
 
 
     // Verbosity
@@ -326,6 +379,16 @@ public:
 
         //- Create polyMesh, with cell zones
         autoPtr<polyMesh> mesh(const IOobject& io) const;
+
+
+    // Housekeeping
+
+        //- Old (v2106 and earlier) uniform scaling factor
+        //  \deprecated use inplacePointTransforms or globalPosition instead
+        scalar scaleFactor() const
+        {
+            return scaling_.x();
+        }
 };
 
 
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshCheck.C b/src/mesh/blockMesh/blockMesh/blockMeshCheck.C
index e38c73ec6ef638a8099dbbeb12b33b6dbc33f31b..086c4ceb5c7a1e9de0f7e631e35ac3e97ef82fd2 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshCheck.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshCheck.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -105,20 +106,37 @@ void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
     }
 
     // Check curved-face/block-face correspondence
-    forAll(faces_, cfi)
+    for (const blockFace& cface : faces_)
     {
+        const face& cf = cface.vertices();
+
         bool found = false;
 
-        forAll(faces, fi)
+        if (cf.size() == 2)
         {
-            found = faces_[cfi].compare(faces[fi]) != 0;
-            if (found) break;
+            const label bi = cf[0];
+            const label fi = cf[1];
+
+            found =
+            (
+                bi >= 0 && bi < blocks.size()
+             && fi >= 0 && fi < blocks[bi].blockShape().nFaces()
+            );
+        }
+
+        if (!found)
+        {
+            for (const face& bf : faces)
+            {
+                found = cface.compare(bf) != 0;
+                if (found) break;
+            }
         }
 
         if (!found)
         {
             Info<< "    Curved face ";
-            faces_[cfi].write(Info, dict);
+            cface.write(Info, dict);
             Info<< "    does not correspond to a block face." << endl;
             ok = false;
         }
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshCreate.C b/src/mesh/blockMesh/blockMesh/blockMeshCreate.C
index 3ab59809e61bd0f27d35cce5f3ab81cf783709be..0835309ce4029382dbf3bb0ba768742f63c70db5 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshCreate.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshCreate.C
@@ -36,64 +36,85 @@ void Foam::blockMesh::createPoints() const
 {
     const blockList& blocks = *this;
 
+    const vector scaleTot
+    (
+        prescaling_.x() * scaling_.x(),
+        prescaling_.y() * scaling_.y(),
+        prescaling_.z() * scaling_.z()
+    );
+
     if (verbose_)
     {
-        Info<< "Creating points with scale " << scaleFactor_ << endl;
+        Info<< "Creating points with scale " << scaleTot << endl;
     }
 
-    points_.setSize(nPoints_);
+    points_.resize(nPoints_);
 
     forAll(blocks, blocki)
     {
         const pointField& blockPoints = blocks[blocki].points();
 
+        const labelSubList pointAddr
+        (
+            mergeList_,
+            blockPoints.size(),
+            blockOffsets_[blocki]
+        );
+
+        UIndirectList<point>(points_, pointAddr) = blockPoints;
+
         if (verbose_)
         {
-            const label nx = blocks[blocki].density().x();
-            const label ny = blocks[blocki].density().y();
-            const label nz = blocks[blocki].density().z();
-
-            label v0 = blocks[blocki].pointLabel(0, 0, 0);
-            label vi1 = blocks[blocki].pointLabel(1, 0, 0);
-            scalar diStart = mag(blockPoints[vi1] - blockPoints[v0]);
-
-            label vinM1 = blocks[blocki].pointLabel(nx-1, 0, 0);
-            label vin = blocks[blocki].pointLabel(nx, 0, 0);
-            scalar diFinal = mag(blockPoints[vin] - blockPoints[vinM1]);
-
-            label vj1 = blocks[blocki].pointLabel(0, 1, 0);
-            scalar djStart = mag(blockPoints[vj1] - blockPoints[v0]);
-            label vjnM1 = blocks[blocki].pointLabel(0, ny-1, 0);
-            label vjn = blocks[blocki].pointLabel(0, ny, 0);
-            scalar djFinal = mag(blockPoints[vjn] - blockPoints[vjnM1]);
-
-            label vk1 = blocks[blocki].pointLabel(0, 0, 1);
-            scalar dkStart = mag(blockPoints[vk1] - blockPoints[v0]);
-            label vknM1 = blocks[blocki].pointLabel(0, 0, nz-1);
-            label vkn = blocks[blocki].pointLabel(0, 0, nz);
-            scalar dkFinal = mag(blockPoints[vkn] - blockPoints[vknM1]);
-
-            Info<< "    Block " << blocki << " cell size :" << nl
-                << "        i : " << scaleFactor_*diStart << " .. "
-                << scaleFactor_*diFinal << nl
-                << "        j : " << scaleFactor_*djStart << " .. "
-                << scaleFactor_*djFinal << nl
-                << "        k : " << scaleFactor_*dkStart << " .. "
-                << scaleFactor_*dkFinal << nl
-                << endl;
-        }
+            Info<< "    Block " << blocki << " cell size :" << nl;
 
-        forAll(blockPoints, blockPointi)
-        {
-            points_
-            [
-                mergeList_
-                [
-                    blockOffsets_[blocki] + blockPointi
-                ]
-            ] = scaleFactor_ * blockPoints[blockPointi];
+            const label v0 = blocks[blocki].pointLabel(0, 0, 0);
+
+            {
+                const label nx = blocks[blocki].density().x();
+                const label v1 = blocks[blocki].pointLabel(1, 0, 0);
+                const label vn = blocks[blocki].pointLabel(nx, 0, 0);
+                const label vn1 = blocks[blocki].pointLabel(nx-1, 0, 0);
+
+                const scalar cwBeg = mag(blockPoints[v1] - blockPoints[v0]);
+                const scalar cwEnd = mag(blockPoints[vn] - blockPoints[vn1]);
+
+                Info<< "        i : "
+                    << cwBeg*scaleTot.x() << " .. "
+                    << cwEnd*scaleTot.x() << nl;
+            }
+
+            {
+                const label ny = blocks[blocki].density().y();
+                const label v1 = blocks[blocki].pointLabel(0, 1, 0);
+                const label vn = blocks[blocki].pointLabel(0, ny, 0);
+                const label vn1 = blocks[blocki].pointLabel(0, ny-1, 0);
+
+                const scalar cwBeg = mag(blockPoints[v1] - blockPoints[v0]);
+                const scalar cwEnd = mag(blockPoints[vn] - blockPoints[vn1]);
+
+                Info<< "        j : "
+                    << cwBeg*scaleTot.y() << " .. "
+                    << cwEnd*scaleTot.y() << nl;
+            }
+
+            {
+                const label nz = blocks[blocki].density().z();
+                const label v1 = blocks[blocki].pointLabel(0, 0, 1);
+                const label vn = blocks[blocki].pointLabel(0, 0, nz);
+                const label vn1 = blocks[blocki].pointLabel(0, 0, nz-1);
+
+                const scalar cwBeg = mag(blockPoints[v1] - blockPoints[v0]);
+                const scalar cwEnd = mag(blockPoints[vn] - blockPoints[vn1]);
+
+                Info<< "        k : "
+                    << cwBeg*scaleTot.z() << " .. "
+                    << cwEnd*scaleTot.z() << nl;
+            }
+            Info<< endl;
         }
     }
+
+    inplacePointTransforms(points_);
 }
 
 
@@ -107,24 +128,22 @@ void Foam::blockMesh::createCells() const
         Info<< "Creating cells" << endl;
     }
 
-    cells_.setSize(nCells_);
+    cells_.resize(nCells_);
 
     label celli = 0;
 
+    labelList cellPoints(8);  // Hex cells - 8 points
+
     forAll(blocks, blocki)
     {
-        const auto& blockCells = blocks[blocki].cells();
-
-        forAll(blockCells, blockCelli)
+        for (const hexCell& blockCell : blocks[blocki].cells())
         {
-            labelList cellPoints(blockCells[blockCelli].size());
-
             forAll(cellPoints, cellPointi)
             {
                 cellPoints[cellPointi] =
                     mergeList_
                     [
-                        blockCells[blockCelli][cellPointi]
+                        blockCell[cellPointi]
                       + blockOffsets_[blocki]
                     ];
             }
@@ -170,7 +189,7 @@ Foam::faceList Foam::blockMesh::createPatchFaces
 
 
     faceList patchFaces(nFaces);
-    face quadFace(4);
+    face quad(4);
     label faceLabel = 0;
 
     forAll(patchTopologyFaces, patchTopologyFaceLabel)
@@ -195,7 +214,7 @@ Foam::faceList Foam::blockMesh::createPatchFaces
                     // Lookup the face points
                     // and collapse duplicate point labels
 
-                    quadFace[0] =
+                    quad[0] =
                         mergeList_
                         [
                             blockPatchFaces[blockFaceLabel][0]
@@ -211,33 +230,33 @@ Foam::faceList Foam::blockMesh::createPatchFaces
                         facePointLabel++
                     )
                     {
-                        quadFace[nUnique] =
+                        quad[nUnique] =
                             mergeList_
                             [
                                 blockPatchFaces[blockFaceLabel][facePointLabel]
                               + blockOffsets_[blocki]
                             ];
 
-                        if (quadFace[nUnique] != quadFace[nUnique-1])
+                        if (quad[nUnique] != quad[nUnique-1])
                         {
                             nUnique++;
                         }
                     }
 
-                    if (quadFace[nUnique-1] == quadFace[0])
+                    if (quad[nUnique-1] == quad[0])
                     {
                         nUnique--;
                     }
 
                     if (nUnique == 4)
                     {
-                        patchFaces[faceLabel++] = quadFace;
+                        patchFaces[faceLabel++] = quad;
                     }
                     else if (nUnique == 3)
                     {
                         patchFaces[faceLabel++] = face
                         (
-                            labelList::subList(quadFace, 3)
+                            labelList::subList(quad, 3)
                         );
                     }
                     // else the face has collapsed to an edge or point
@@ -246,7 +265,7 @@ Foam::faceList Foam::blockMesh::createPatchFaces
         }
     }
 
-    patchFaces.setSize(faceLabel);
+    patchFaces.resize(faceLabel);
 
     return patchFaces;
 }
@@ -261,7 +280,7 @@ void Foam::blockMesh::createPatches() const
         Info<< "Creating patches" << endl;
     }
 
-    patches_.setSize(topoPatches.size());
+    patches_.resize(topoPatches.size());
 
     forAll(topoPatches, patchi)
     {
@@ -272,6 +291,66 @@ void Foam::blockMesh::createPatches() const
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+const Foam::polyMesh& Foam::blockMesh::topology() const
+{
+    if (!topologyPtr_)
+    {
+        FatalErrorInFunction
+            << "topology not allocated"
+            << exit(FatalError);
+    }
+
+    return *topologyPtr_;
+}
+
+
+Foam::refPtr<Foam::polyMesh>
+Foam::blockMesh::topology(bool applyTransform) const
+{
+    const polyMesh& origTopo = topology();
+
+    if (applyTransform && hasPointTransforms())
+    {
+        // Copy mesh components
+
+        IOobject newIO(origTopo, IOobject::NO_READ, IOobject::NO_WRITE);
+        newIO.registerObject(false);
+
+        pointField newPoints(origTopo.points());
+        inplacePointTransforms(newPoints);
+
+        auto ttopoMesh = refPtr<polyMesh>::New
+        (
+            newIO,
+            std::move(newPoints),
+            faceList(origTopo.faces()),
+            labelList(origTopo.faceOwner()),
+            labelList(origTopo.faceNeighbour())
+        );
+        auto& topoMesh = ttopoMesh.ref();
+
+        // Patches
+        const polyBoundaryMesh& pbmOld = origTopo.boundaryMesh();
+        const polyBoundaryMesh& pbmNew = topoMesh.boundaryMesh();
+
+        PtrList<polyPatch> newPatches(pbmOld.size());
+
+        forAll(pbmOld, patchi)
+        {
+            newPatches.set(patchi, pbmOld[patchi].clone(pbmNew));
+        }
+
+        topoMesh.addPatches(newPatches);
+
+        return ttopoMesh;
+    }
+    else
+    {
+        return origTopo;
+    }
+}
+
+
 Foam::autoPtr<Foam::polyMesh>
 Foam::blockMesh::mesh(const IOobject& io) const
 {
@@ -285,7 +364,7 @@ Foam::blockMesh::mesh(const IOobject& io) const
     auto meshPtr = autoPtr<polyMesh>::New
     (
         io,
-        pointField(blkMesh.points()),   // Copy, could we re-use space?
+        pointField(blkMesh.points()),   // Use a copy of the block points
         blkMesh.cells(),
         blkMesh.patches(),
         blkMesh.patchNames(),
@@ -381,9 +460,9 @@ Foam::blockMesh::mesh(const IOobject& io) const
             );
         }
 
-        pmesh.pointZones().resize(0);
-        pmesh.faceZones().resize(0);
-        pmesh.cellZones().resize(0);
+        pmesh.pointZones().clear();
+        pmesh.faceZones().clear();
+        pmesh.cellZones().clear();
         pmesh.addZones(List<pointZone*>(), List<faceZone*>(), cz);
     }
 
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshMergeGeometrical.C b/src/mesh/blockMesh/blockMesh/blockMeshMergeGeometrical.C
index 7071290a989da07627ebab513f1eb09dbec5a95f..cfb733bc89acca64ebe7a19e66b58ed65988d6c0 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshMergeGeometrical.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshMergeGeometrical.C
@@ -39,7 +39,7 @@ void Foam::blockMesh::calcGeometricalMerge()
         Info<< "Creating block offsets" << endl;
     }
 
-    blockOffsets_.setSize(blocks.size());
+    blockOffsets_.resize(blocks.size());
 
     nPoints_ = 0;
     nCells_  = 0;
@@ -58,21 +58,28 @@ void Foam::blockMesh::calcGeometricalMerge()
         Info<< "Creating merge list (geometric search).." << flush;
     }
 
-    // set unused to -1
-    mergeList_.setSize(nPoints_);
+    // Set unused to -1
+    mergeList_.resize(nPoints_);
     mergeList_ = -1;
 
 
-    const pointField& blockPoints = topology().points();
-    const cellList& blockCells = topology().cells();
-    const faceList& blockFaces = topology().faces();
-    const labelList& faceOwnerBlocks = topology().faceOwner();
+    // Block mesh topology
+    const polyMesh& topoMesh = topology();
 
-    // For efficiency, create merge pairs in the first pass
-    labelListListList glueMergePairs(blockFaces.size());
+    const pointField& blockPoints = topoMesh.points();
+    const cellList& blockCells = topoMesh.cells();
+    const faceList& blockFaces = topoMesh.faces();
+    const labelList& faceOwnerBlocks = topoMesh.faceOwner();
+    const labelList& faceNeighbourBlocks = topoMesh.faceNeighbour();
 
-    const labelList& faceNeighbourBlocks = topology().faceNeighbour();
+    const faceList::subList blockinternalFaces
+    (
+        blockFaces,
+        topoMesh.nInternalFaces()
+    );
 
+    // For efficiency, create merge pairs in the first pass
+    labelListListList glueMergePairs(blockFaces.size());
 
     forAll(blockFaces, blockFaceLabel)
     {
@@ -179,7 +186,7 @@ void Foam::blockMesh::calcGeometricalMerge()
         sqrMergeTol /= 10.0;
 
 
-        if (topology().isInternalFace(blockFaceLabel))
+        if (topoMesh.isInternalFace(blockFaceLabel))
         {
         label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
         const pointField& blockNpoints = blocks[blockNlabel].points();
@@ -306,12 +313,6 @@ void Foam::blockMesh::calcGeometricalMerge()
     }
 
 
-    const faceList::subList blockinternalFaces
-    (
-        blockFaces,
-        topology().nInternalFaces()
-    );
-
     bool changedPointMerge = false;
     label nPasses = 0;
 
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshMergeTopological.C b/src/mesh/blockMesh/blockMesh/blockMeshMergeTopological.C
index 15c3e5186866824d65a083987cec19e114c2d05a..19d2ffc55c395216f6feb70d99043162428b1acd 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshMergeTopological.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshMergeTopological.C
@@ -313,7 +313,7 @@ void Foam::blockMesh::calcTopologicalMerge()
         Info<< "Creating block offsets" << endl;
     }
 
-    blockOffsets_.setSize(blocks.size());
+    blockOffsets_.resize(blocks.size());
 
     nPoints_ = 0;
     nCells_  = 0;
@@ -335,11 +335,13 @@ void Foam::blockMesh::calcTopologicalMerge()
     mergeList_.setSize(nPoints_, -1);
 
     // Block mesh topology
-    const pointField& topoPoints = topology().points();
-    const cellList& topoCells = topology().cells();
-    const faceList& topoFaces = topology().faces();
-    const labelList& topoFaceOwn = topology().faceOwner();
-    const labelList& topoFaceNei = topology().faceNeighbour();
+    const polyMesh& topoMesh = topology();
+
+    const pointField& topoPoints = topoMesh.points();
+    const cellList& topoCells = topoMesh.cells();
+    const faceList& topoFaces = topoMesh.faces();
+    const labelList& topoFaceOwn = topoMesh.faceOwner();
+    const labelList& topoFaceNei = topoMesh.faceNeighbour();
 
     // Topological merging is only necessary for internal block faces
     // Note edge and face collapse may apply to boundary faces
@@ -347,7 +349,7 @@ void Foam::blockMesh::calcTopologicalMerge()
     const faceList::subList topoInternalFaces
     (
         topoFaces,
-        topology().nInternalFaces()
+        topoMesh.nInternalFaces()
     );
 
     List<Pair<label>> mergeBlockP(topoInternalFaces.size());
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshTopology.C b/src/mesh/blockMesh/blockMesh/blockMeshTopology.C
index 5e3a01fa33c5cb56d2fe7c69d164ae1e3f5ebdc3..c264e96701cd7d57d437da6c20e7fd57c34c807f 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshTopology.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshTopology.C
@@ -33,17 +33,27 @@ License
 #include "emptyPolyPatch.H"
 #include "cyclicPolyPatch.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Replace (<block> <face>) face description
+// with the corresponding block face
+// - otherwise check point labels are in range
 
 template<class Source>
-void Foam::blockMesh::checkPatchLabels
+static void rewritePatchLabels
 (
     const Source& source,
     const word& patchName,
-    const pointField& points,
+    const PtrList<block>& blocks,
+    const label nPoints,
     faceList& patchFaces
-) const
+)
 {
+    const label nBlocks = blocks.size();
+
     forAll(patchFaces, facei)
     {
         face& f = patchFaces[facei];
@@ -55,40 +65,41 @@ void Foam::blockMesh::checkPatchLabels
             const label bi = f[0];
             const label fi = f[1];
 
-            if (bi >= size())
+            if (bi < 0 || bi >= nBlocks)
             {
                 FatalIOErrorInFunction(source)
-                    << "Block index out of range for patch face " << f << nl
-                    << "    Number of blocks = " << size()
-                    << ", index = " << f[0] << nl
+                    << "Block index out of range."
+                    << " Patch face (" << bi << ' ' << fi << ")\n"
+                    << "    Number of blocks = " << nBlocks
+                    << ", block index = " << bi << nl
                     << "    on patch " << patchName << ", face " << facei
                     << exit(FatalIOError);
             }
-            else if (fi >= operator[](bi).blockShape().faces().size())
+            else if (fi >= blocks[bi].blockShape().nFaces())
             {
                 FatalIOErrorInFunction(source)
-                    << "Block face index out of range for patch face " << f
-                    << nl
+                    << "Block face index out of range."
+                    << " Patch face (" << bi << ' ' << fi << ")\n"
                     << "    Number of block faces = "
-                    << operator[](bi).blockShape().faces().size()
-                    << ", index = " << f[1] << nl
+                    << blocks[bi].blockShape().nFaces()
+                    << ", face index = " << fi << nl
                     << "    on patch " << patchName << ", face " << facei
                     << exit(FatalIOError);
             }
             else
             {
-                f = operator[](bi).blockShape().faces()[fi];
+                f = blocks[bi].blockShape().face(fi);
             }
         }
         else
         {
             for (const label pointi : f)
             {
-                if (pointi < 0 || pointi >= points.size())
+                if (pointi < 0 || pointi >= nPoints)
                 {
                     FatalIOErrorInFunction(source)
                         << "Point label " << pointi
-                        << " out of range 0.." << points.size() - 1 << nl
+                        << " out of range 0.." << (nPoints - 1) << nl
                         << "    on patch " << patchName << ", face " << facei
                         << exit(FatalIOError);
                 }
@@ -97,6 +108,10 @@ void Foam::blockMesh::checkPatchLabels
     }
 }
 
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::blockMesh::readPatches
 (
@@ -112,26 +127,20 @@ void Foam::blockMesh::readPatches
     varDict.merge(meshDescription.subOrEmptyDict("namedBlocks"));
 
 
-    ITstream& patchStream(meshDescription.lookup("patches"));
+    ITstream& patchStream = meshDescription.lookup("patches");
 
     // Read number of patches in mesh
     label nPatches = 0;
 
-    token firstToken(patchStream);
-
-    if (firstToken.isLabel())
+    if (patchStream.peek().isLabel())
     {
-        nPatches = firstToken.labelToken();
+        patchStream >> nPatches;
 
         tmpBlocksPatches.setSize(nPatches);
         patchNames.setSize(nPatches);
         patchTypes.setSize(nPatches);
         nbrPatchNames.setSize(nPatches);
     }
-    else
-    {
-        patchStream.putBack(firstToken);
-    }
 
     // Read beginning of blocks
     patchStream.readBegin("patches");
@@ -175,15 +184,16 @@ void Foam::blockMesh::readPatches
             }
         }
 
-        checkPatchLabels
+        rewritePatchLabels
         (
             patchStream,
             patchNames[nPatches],
-            vertices_,
+            *this,
+            vertices_.size(),
             tmpBlocksPatches[nPatches]
         );
 
-        nPatches++;
+        ++nPatches;
 
 
         // Split old style cyclics
@@ -298,67 +308,57 @@ void Foam::blockMesh::readBoundary
         );
 
 
-        checkPatchLabels
+        rewritePatchLabels
         (
             patchInfo.dict(),
             patchNames[patchi],
-            vertices_,
+            *this,
+            vertices_.size(),
             tmpBlocksPatches[patchi]
         );
     }
 }
 
 
-void Foam::blockMesh::createCellShapes
-(
-    cellShapeList& tmpBlockCells
-)
+Foam::cellShapeList Foam::blockMesh::getBlockShapes() const
 {
     const blockMesh& blocks = *this;
 
-    tmpBlockCells.setSize(blocks.size());
+    cellShapeList shapes(blocks.size());
+
     forAll(blocks, blocki)
     {
-        tmpBlockCells[blocki] = blocks[blocki].blockShape();
+        shapes[blocki] = blocks[blocki].blockShape();
     }
+
+    return shapes;
 }
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
+Foam::autoPtr<Foam::polyMesh>
+Foam::blockMesh::createTopology
 (
     const IOdictionary& meshDescription,
     const word& regionName
 )
 {
-    blockList& blocks = *this;
-
     word defaultPatchName = "defaultFaces";
     word defaultPatchType = emptyPolyPatch::typeName;
 
     // Read the names/types for the unassigned patch faces
     // this is a bit heavy handed (and ugly), but there is currently
     // no easy way to rename polyMesh patches subsequently
-    if (const dictionary* dictPtr = meshDescription.findDict("defaultPatch"))
+    if (const dictionary* dictptr = meshDescription.findDict("defaultPatch"))
     {
-        dictPtr->readIfPresent("name", defaultPatchName);
-        dictPtr->readIfPresent("type", defaultPatchType);
+        dictptr->readIfPresent("name", defaultPatchName);
+        dictptr->readIfPresent("type", defaultPatchType);
     }
 
-    // Optional 'scale' factor. Was 'convertToMeters' until OCT-2008
-    meshDescription.readIfPresentCompat
-    (
-        "scale",
-        {{"convertToMeters", 1012}},  // Mark as changed from 2010 onwards
-        scaleFactor_
-    );
 
-    // Treat (scale <= 0) as scaling == 1 (no scaling).
-    if (scaleFactor_ <= 0)
-    {
-        scaleFactor_ = 1.0;
-    }
+    // Scaling, transformations
+    readPointTransforms(meshDescription);
 
     // Read the block edges
     if (meshDescription.found("edges"))
@@ -376,9 +376,13 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
 
         edges_.transfer(edges);
     }
-    else if (verbose_)
+    else
     {
-        Info<< "No non-linear block edges defined" << endl;
+        edges_.clear();
+        if (verbose_)
+        {
+            Info<< "No non-linear block edges defined" << endl;
+        }
     }
 
 
@@ -398,9 +402,13 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
 
         faces_.transfer(faces);
     }
-    else if (verbose_)
+    else
     {
-        Info<< "No non-planar block faces defined" << endl;
+        faces_.clear();
+        if (verbose_)
+        {
+            Info<< "No non-planar block faces defined" << endl;
+        }
     }
 
 
@@ -420,21 +428,25 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
     }
 
 
-    autoPtr<polyMesh> blockMeshPtr;
+    // Create patch information
+
+    faceListList tmpBlocksPatches;
+    wordList patchNames;
+    PtrList<dictionary> patchDicts;
 
-    // Create the patches
 
     if (verbose_)
     {
-        Info<< "Creating topology patches" << endl;
+        Info<< nl << "Creating topology patches - ";
     }
 
     if (meshDescription.found("patches"))
     {
-        Info<< nl << "Reading patches section" << endl;
+        if (verbose_)
+        {
+            Info<< "from patches section" << endl;
+        }
 
-        faceListList tmpBlocksPatches;
-        wordList patchNames;
         wordList patchTypes;
         wordList nbrPatchNames;
 
@@ -447,16 +459,13 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
             nbrPatchNames
         );
 
-        Info<< nl << "Creating block mesh topology" << endl;
-
-        cellShapeList tmpBlockCells(blocks.size());
-        createCellShapes(tmpBlockCells);
-
-
-        Info<< nl << "Reading physicalType from existing boundary file" << endl;
+        if (verbose_)
+        {
+            Info<< nl
+                << "Reading physicalType from existing boundary file" << endl;
+        }
 
-        PtrList<dictionary> patchDicts(patchNames.size());
-        word defaultFacesType;
+        patchDicts.resize(patchNames.size());
 
         preservePatchTypes
         (
@@ -495,38 +504,19 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
                     << exit(FatalIOError);
             }
 
-            // Override neighbourpatch name
-            if (nbrPatchNames[patchi] != word::null)
+            // Override neighbour patch name
+            if (!nbrPatchNames[patchi].empty())
             {
                 dict.set("neighbourPatch", nbrPatchNames[patchi]);
             }
         }
-
-        blockMeshPtr = autoPtr<polyMesh>::New
-        (
-            IOobject
-            (
-                regionName,
-                meshDescription.time().constant(),
-                meshDescription.time(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            pointField(vertices_),   // Copy these points, do NOT move
-            tmpBlockCells,
-            tmpBlocksPatches,
-            patchNames,
-            patchDicts,
-            defaultPatchName,
-            defaultPatchType
-        );
     }
     else if (meshDescription.found("boundary"))
     {
-        wordList patchNames;
-        faceListList tmpBlocksPatches;
-        PtrList<dictionary> patchDicts;
+        if (verbose_)
+        {
+            Info<< "from boundary section" << endl;
+        }
 
         readBoundary
         (
@@ -535,32 +525,45 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
             tmpBlocksPatches,
             patchDicts
         );
+    }
+    else
+    {
+        if (verbose_)
+        {
+            Info<< "with default boundary only!!" << endl;
+        }
+    }
 
-        Info<< nl << "Creating block mesh topology" << endl;
 
-        cellShapeList tmpBlockCells(blocks.size());
-        createCellShapes(tmpBlockCells);
+    if (verbose_)
+    {
+        Info<< nl << "Creating block mesh topology";
+        if (hasPointTransforms())
+        {
+            Info<< " - scaling/transform applied later";
+        }
+        Info<< endl;
+    }
 
-        blockMeshPtr = autoPtr<polyMesh>::New
+    auto blockMeshPtr = autoPtr<polyMesh>::New
+    (
+        IOobject
         (
-            IOobject
-            (
-                regionName,
-                meshDescription.time().constant(),
-                meshDescription.time(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            pointField(vertices_),   // Copy these points, do NOT move
-            tmpBlockCells,
-            tmpBlocksPatches,
-            patchNames,
-            patchDicts,
-            defaultPatchName,
-            defaultPatchType
-        );
-    }
+            regionName,
+            meshDescription.time().constant(),
+            meshDescription.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        pointField(vertices_),   // Use a copy of vertices
+        getBlockShapes(),
+        tmpBlocksPatches,
+        patchNames,
+        patchDicts,
+        defaultPatchName,
+        defaultPatchType
+    );
 
     check(*blockMeshPtr, meshDescription);
 
diff --git a/src/mesh/blockMesh/blocks/block/block.H b/src/mesh/blockMesh/blocks/block/block.H
index 6bac8acac0d0d69e16597da0cefa04c637654d94..2d2c2e6ecec17180437b58ea51e0b0c734da88c4 100644
--- a/src/mesh/blockMesh/blocks/block/block.H
+++ b/src/mesh/blockMesh/blocks/block/block.H
@@ -215,14 +215,14 @@ public:
     // Access
 
         //- The points for filling the block
-        inline const pointField& points() const;
+        inline const pointField& points() const noexcept;
 
         //- The hex cells for filling the block
         inline const List<hexCell>& cells() const;
 
         //- The boundary patch faces for the block
         inline const FixedList<List<FixedList<label, 4>>, 6>&
-        boundaryPatches() const;
+        boundaryPatches() const noexcept;
 
 
     // Mesh Components
diff --git a/src/mesh/blockMesh/blocks/block/blockCreate.C b/src/mesh/blockMesh/blocks/block/blockCreate.C
index dbc5b7f5c704f528ebd0aefe0af1c298dd45e01a..175b9331187aad9c4272778471aaebc5e4c4d83c 100644
--- a/src/mesh/blockMesh/blocks/block/blockCreate.C
+++ b/src/mesh/blockMesh/blocks/block/blockCreate.C
@@ -65,7 +65,7 @@ void Foam::block::createPoints()
     // List of edge point and weighting factors
     pointField p[12];
     scalarList w[12];
-    label nCurvedEdges = edgesPointsWeights(p, w);
+    const int nCurvedEdges = edgesPointsWeights(p, w);
 
     points_.resize(nPoints());
 
diff --git a/src/mesh/blockMesh/blocks/block/blockI.H b/src/mesh/blockMesh/blocks/block/blockI.H
index e8aea8b5ca82517693507c81e7478135385164fe..f7fa311f6cf40a7cd90cf17119a72378d839dc84 100644
--- a/src/mesh/blockMesh/blocks/block/blockI.H
+++ b/src/mesh/blockMesh/blocks/block/blockI.H
@@ -28,7 +28,7 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const Foam::pointField& Foam::block::points() const
+inline const Foam::pointField& Foam::block::points() const noexcept
 {
     return points_;
 }
@@ -46,7 +46,7 @@ inline const Foam::List<Foam::hexCell>& Foam::block::cells() const
 
 
 inline const Foam::FixedList<Foam::List<Foam::FixedList<Foam::label, 4>>, 6>&
-Foam::block::boundaryPatches() const
+Foam::block::boundaryPatches() const noexcept
 {
     return blockPatches_;
 }