diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H b/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H
new file mode 100644
index 0000000000..3fb49e2982
--- /dev/null
+++ b/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H
@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::ijkAddressing
+
+Description
+    A simple i-j-k (row-major order) to linear addressing.
+
+SourceFiles
+    ijkAddressingI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ijkAddressing_H
+#define ijkAddressing_H
+
+#include "labelVector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class ijkAddressing Declaration
+\*---------------------------------------------------------------------------*/
+
+class ijkAddressing
+{
+    // Private Data
+
+        //- The number of cells in the i,j,k directions.
+        labelVector sizes_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct zero-size addressing
+        inline ijkAddressing();
+
+        //- Construct with addressing
+        inline explicit ijkAddressing(const labelVector& ijk);
+
+        //- Construct with addressing components
+        inline ijkAddressing(const label ni, const label nj, const label nk);
+
+
+    // Member Functions
+
+    // Access
+
+        //- Return the i,j,k addressing sizes
+        inline const labelVector& sizes() const;
+
+        //- Return the i,j,k addressing sizes for modification
+        inline labelVector& sizes();
+
+        //- Return the total i*j*k size
+        inline label size() const;
+
+        //- Addressing is considered empty if any component is zero
+        inline bool empty() const;
+
+        //- Reset to (0,0,0) sizing
+        inline void clear();
+
+        //- Change the sizing parameters
+        inline void reset(const label ni, const label nj, const label nk);
+
+        //- Linear addressing index (offset) for an i,j,k position.
+        inline label index(const label i, const label j, const label k) const;
+
+
+    // Checks
+
+        //- Check indices are within valid ni,nj,nk range.
+        //  Optionally allow an extra index for point addressing
+        inline void checkIndex
+        (
+            const label i,
+            const label j,
+            const label k,
+            const bool isPoint = false
+        ) const;
+
+        //- Check that sizes() are valid
+        inline void checkSizes() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "ijkAddressingI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H b/src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H
new file mode 100644
index 0000000000..3c8f442c8f
--- /dev/null
+++ b/src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H
@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::ijkAddressing::ijkAddressing()
+:
+    sizes_(0, 0, 0)
+{}
+
+
+inline Foam::ijkAddressing::ijkAddressing(const labelVector& ijk)
+:
+    sizes_(ijk)
+{
+    #ifdef FULLDEBUG
+    checkSizes();
+    #endif
+}
+
+
+inline Foam::ijkAddressing::ijkAddressing
+(
+    const label ni,
+    const label nj,
+    const label nk
+)
+:
+    sizes_(ni, nj, nk)
+{
+    #ifdef FULLDEBUG
+    checkSizes();
+    #endif
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::labelVector& Foam::ijkAddressing::sizes() const
+{
+    return sizes_;
+}
+
+
+inline Foam::labelVector& Foam::ijkAddressing::sizes()
+{
+    return sizes_;
+}
+
+
+inline Foam::label Foam::ijkAddressing::size() const
+{
+    // Could also use  cmptProduct(sizes_);
+    return (sizes_.x() * sizes_.y() * sizes_.z());
+}
+
+
+inline bool Foam::ijkAddressing::empty() const
+{
+    return (!sizes_.x() || !sizes_.y() || !sizes_.z());
+}
+
+
+inline void Foam::ijkAddressing::clear()
+{
+    sizes_ = Zero;
+}
+
+
+inline void Foam::ijkAddressing::reset
+(
+    const label ni,
+    const label nj,
+    const label nk
+)
+{
+    sizes_.x() = ni;
+    sizes_.y() = nj;
+    sizes_.z() = nk;
+
+    #ifdef FULLDEBUG
+    checkSizes();
+    #endif
+}
+
+
+inline Foam::label Foam::ijkAddressing::index
+(
+    const label i,
+    const label j,
+    const label k
+) const
+{
+    #ifdef FULLDEBUG
+    checkIndex(i, j, k);
+    #endif
+
+    return (i + (sizes_.x() * (j + (sizes_.y() * k))));
+}
+
+
+inline void Foam::ijkAddressing::checkIndex
+(
+    const label i,
+    const label j,
+    const label k,
+    const bool isPoint
+) const
+{
+    const label extra = (isPoint ? 1 : 0);
+
+    if (i < 0 || i >= (sizes_.x() + extra))
+    {
+        FatalErrorInFunction
+            << "The i-index " << i
+            << " is out of range [0," << (sizes_.x() + extra) << ']' << nl
+            << abort(FatalError);
+    }
+    if (j < 0 || j >= (sizes_.y() + extra))
+    {
+        FatalErrorInFunction
+            << "The j-index " << j
+            << " is out of range [0," << (sizes_.y() + extra) << ']' << nl
+            << abort(FatalError);
+    }
+    if (k < 0 || k >= (sizes_.z() + extra))
+    {
+        FatalErrorInFunction
+            << "The k-index " << k
+            << " is out of range [0," << (sizes_.z() + extra) << ']' << nl
+            << abort(FatalError);
+    }
+}
+
+
+inline void Foam::ijkAddressing::checkSizes() const
+{
+    if (sizes_.x() < 0)
+    {
+        FatalErrorInFunction
+            << "The i-size is negative" << nl
+            << abort(FatalError);
+    }
+    if (sizes_.y() < 0)
+    {
+        FatalErrorInFunction
+            << "The j-size is negative" << nl
+            << abort(FatalError);
+    }
+    if (sizes_.z() < 0)
+    {
+        FatalErrorInFunction
+            << "The k-size is negative" << nl
+            << abort(FatalError);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkMesh.H b/src/OpenFOAM/meshes/ijkMesh/ijkMesh.H
new file mode 100644
index 0000000000..976e3a6c9d
--- /dev/null
+++ b/src/OpenFOAM/meshes/ijkMesh/ijkMesh.H
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::ijkMesh
+
+Description
+    A simple i-j-k (row-major order) to linear addressing for a
+    rectilinear mesh. Since the underlying mesh is rectilinear, some
+    mesh-related sizing information can be derived directly from the
+    addressing information.
+
+SourceFiles
+    ijkMeshI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ijkMesh_H
+#define ijkMesh_H
+
+#include "ijkAddressing.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class ijkMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class ijkMesh
+:
+    public ijkAddressing
+{
+public:
+
+    // Constructors
+
+        //- Construct zero-sized
+        inline ijkMesh();
+
+        //- Construct with addressing
+        inline explicit ijkMesh(const labelVector& ijk);
+
+        //- Construct with addressing
+        inline ijkMesh(const label nx, const label ny, const label nz);
+
+
+    // Member Functions
+
+    // Access
+
+        //- The number of mesh points (nx+1)*(ny+1)*(nz+1) in the i-j-k mesh.
+        inline label nPoints() const;
+
+        //- The number of mesh cells (nx*ny*nz) in the i-j-k mesh
+        //  This is the same as the ijkAddressing::size()
+        inline label nCells() const;
+
+        //- The total number of mesh faces in the i-j-k mesh
+        inline label nFaces() const;
+
+        //- The number of internal faces in the i-j-k mesh
+        inline label nInternalFaces() const;
+
+        //- The number of boundary faces in the i-j-k mesh
+        inline label nBoundaryFaces() const;
+
+        //- The linear cell index for an i-j-k position - same as index()
+        inline label cellLabel
+        (
+            const label i,
+            const label j,
+            const label k
+        ) const;
+
+        //- The linear point index for an i-j-k position.
+        //  Addressable in the range (ni+1, nj+1, nk+1).
+        inline label pointLabel
+        (
+            const label i,
+            const label j,
+            const label k
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "ijkMeshI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H b/src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H
new file mode 100644
index 0000000000..cc89c2c520
--- /dev/null
+++ b/src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H
@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::ijkMesh::ijkMesh()
+:
+    ijkAddressing()
+{}
+
+
+inline Foam::ijkMesh::ijkMesh(const labelVector& ijk)
+:
+    ijkAddressing(ijk)
+{}
+
+
+inline Foam::ijkMesh::ijkMesh
+(
+    const label nx,
+    const label ny,
+    const label nz
+)
+:
+    ijkAddressing(nx, ny, nz)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::label Foam::ijkMesh::nPoints() const
+{
+    if (ijkAddressing::empty())
+    {
+        return 0;
+    }
+
+    const labelVector& n = ijkAddressing::sizes();
+
+    return ((n.x()+1) * (n.y()+1) * (n.z()+1));
+}
+
+
+inline Foam::label Foam::ijkMesh::nCells() const
+{
+    return ijkAddressing::size();
+}
+
+
+inline Foam::label Foam::ijkMesh::nFaces() const
+{
+    if (ijkAddressing::empty())
+    {
+        return 0;
+    }
+
+    const labelVector& n = ijkAddressing::sizes();
+
+    return
+    (
+        ((n.x()+1) * n.y() * n.z())
+      + ((n.y()+1) * n.z() * n.x())
+      + ((n.z()+1) * n.x() * n.y())
+    );
+}
+
+
+inline Foam::label Foam::ijkMesh::nInternalFaces() const
+{
+    if (ijkAddressing::empty())
+    {
+        return 0;
+    }
+
+    const labelVector& n = ijkAddressing::sizes();
+
+    return
+    (
+        ((n.x()-1) * n.y() * n.z())
+      + ((n.y()-1) * n.z() * n.x())
+      + ((n.z()-1) * n.x() * n.y())
+    );
+}
+
+
+inline Foam::label Foam::ijkMesh::nBoundaryFaces() const
+{
+    if (ijkAddressing::empty())
+    {
+        return 0;
+    }
+
+    const labelVector& n = ijkAddressing::sizes();
+
+    return
+    (
+        (2 * n.y() * n.z())
+      + (2 * n.z() * n.x())
+      + (2 * n.x() * n.y())
+    );
+}
+
+
+inline Foam::label Foam::ijkMesh::cellLabel
+(
+    const label i,
+    const label j,
+    const label k
+) const
+{
+    return ijkAddressing::index(i, j, k);
+}
+
+
+inline Foam::label Foam::ijkMesh::pointLabel
+(
+    const label i,
+    const label j,
+    const label k
+) const
+{
+    #ifdef FULLDEBUG
+    checkIndex(i, j, k, true);
+    #endif
+
+    const labelVector& n = sizes();
+
+    return (i + ((n.x()+1) * (j + (n.y()+1) * k)));
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
index a54d490946..4efc3094dd 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
@@ -133,16 +133,16 @@ Foam::blockDescriptor::blockDescriptor
     const pointField& vertices,
     const blockEdgeList& edges,
     const blockFaceList& faces,
-    const Vector<label>& density,
+    const labelVector& density,
     const UList<gradingDescriptors>& expand,
     const word& zoneName
 )
 :
+    ijkMesh(density),
     vertices_(vertices),
     blockEdges_(edges),
     blockFaces_(faces),
     blockShape_(bshape),
-    density_(density),
     expand_(expand),
     zoneName_(zoneName),
     curvedFaces_(-1),
@@ -173,10 +173,10 @@ Foam::blockDescriptor::blockDescriptor
     Istream& is
 )
 :
+    ijkMesh(),
     vertices_(vertices),
     blockEdges_(edges),
     blockFaces_(faces),
-    density_(0, 0, 0),
     expand_(12, gradingDescriptors()),
     zoneName_(),
     curvedFaces_(-1),
@@ -212,7 +212,7 @@ Foam::blockDescriptor::blockDescriptor
         // New-style: read a list of 3 values
         if (t.pToken() == token::BEGIN_LIST)
         {
-            is >> density_;
+            is >> ijkMesh::sizes();
         }
         else
         {
@@ -225,9 +225,9 @@ Foam::blockDescriptor::blockDescriptor
     else
     {
         // Old-style: read three labels
-        is  >> density_.x()
-            >> density_.y()
-            >> density_.z();
+        is  >> ijkMesh::sizes().x()
+            >> ijkMesh::sizes().y()
+            >> ijkMesh::sizes().z();
     }
 
     is >> t;
@@ -285,14 +285,13 @@ Foam::blockDescriptor::blockDescriptor
 Foam::FixedList<Foam::pointField, 6>
 Foam::blockDescriptor::facePoints(const pointField& points) const
 {
+    const label ni = sizes().x();
+    const label nj = sizes().y();
+    const label nk = sizes().z();
+
     // Caches points for curvature correction
     FixedList<pointField, 6> facePoints;
 
-    // Set local variables for mesh specification
-    const label ni = density_.x();
-    const label nj = density_.y();
-    const label nk = density_.z();
-
     facePoints[0].setSize((nj + 1)*(nk + 1));
     facePoints[1].setSize((nj + 1)*(nk + 1));
 
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
index 0b89d948d3..c3dc94f37f 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
@@ -57,6 +57,7 @@ SourceFiles
 #ifndef blockDescriptor_H
 #define blockDescriptor_H
 
+#include "ijkMesh.H"
 #include "cellShape.H"
 #include "pointField.H"
 #include "scalarList.H"
@@ -75,6 +76,8 @@ namespace Foam
 \*---------------------------------------------------------------------------*/
 
 class blockDescriptor
+:
+    public ijkMesh
 {
     // Private Data
 
@@ -90,9 +93,6 @@ class blockDescriptor
         //- Block shape
         cellShape blockShape_;
 
-        //- The number of cells in the i,j,k directions
-        Vector<label> density_;
-
         //- Expansion ratios in all directions
         List<gradingDescriptors> expand_;
 
@@ -140,7 +140,7 @@ public:
             const pointField& vertices,
             const blockEdgeList& edges,
             const blockFaceList& faces,
-            const Vector<label>& density,
+            const labelVector& density,
             const UList<gradingDescriptors>& expand,
             const word& zoneName = ""
         );
@@ -169,7 +169,7 @@ public:
         inline const cellShape& blockShape() const;
 
         //- Return the mesh density (number of cells) in the i,j,k directions
-        inline const Vector<label>& density() const;
+        inline const labelVector& density() const;
 
         //- Expansion ratios in all directions
         const List<gradingDescriptors>& grading() const;
@@ -177,21 +177,6 @@ public:
         //- Return the (optional) zone name
         inline const word& zoneName() const;
 
-        //- The number of meshed points described: (nx+1) * (ny+1) * (nz+1)
-        inline label nPoints() const;
-
-        //- The number of meshed cells described: (nx * ny * nz)
-        inline label nCells() const;
-
-        //- The number of meshed faces described
-        inline label nFaces() const;
-
-        //- The number of meshed internal faces described
-        inline label nInternalFaces() const;
-
-        //- The number of meshed internal faces described
-        inline label nBoundaryFaces() const;
-
         //- Curved-face labels for each block-face (-1 for flat faces)
         inline const FixedList<label, 6>& curvedFaces() const;
 
@@ -201,34 +186,11 @@ public:
         //- Return block point for local label i
         inline const point& blockPoint(const label i) const;
 
-        //- Check indices are within valid range ni,nj,nk
-        inline void checkIndex
-        (
-            const label i,
-            const label j,
-            const label k
-        ) const;
-
-        //- Cell label offset for a particular i,j,k position
-        inline label cellLabel
-        (
-            const label i,
-            const label j,
-            const label k
-        ) const;
-
-        //- Vertex label offset for a particular i,j,k position
-        inline label pointLabel
-        (
-            const label i,
-            const label j,
-            const label k
-        ) const;
-
         //- Face vertex label offset for a particular i,j,k position
+        //- on hex face (0-5)
         inline label facePointLabel
         (
-            const label facei,
+            const direction facei,
             const label i,
             const label j
         ) const;
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
index 011b7df0e0..9ff03ee636 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
@@ -116,24 +116,25 @@ Foam::label Foam::blockDescriptor::edgesPointsWeights
     scalarList (&edgeWeights)[12]
 ) const
 {
+    const label ni = sizes().x();
+    const label nj = sizes().y();
+    const label nk = sizes().z();
+
     label nCurvedEdges = 0;
 
     // X-direction
-    const label ni = density_.x();
     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
-    const label nj = density_.y();
     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
-    const label nk = density_.z();
     nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 8,  0, 4, nk);
     nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 9,  1, 5, nk);
     nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 10, 2, 6, nk);
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
index 1f3857c95b..593b594100 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
@@ -45,9 +45,9 @@ inline const Foam::cellShape& Foam::blockDescriptor::blockShape() const
 }
 
 
-inline const Foam::Vector<Foam::label>& Foam::blockDescriptor::density() const
+inline const Foam::labelVector& Foam::blockDescriptor::density() const
 {
-    return density_;
+    return ijkMesh::sizes();
 }
 
 
@@ -64,61 +64,6 @@ inline const Foam::word& Foam::blockDescriptor::zoneName() const
 }
 
 
-inline Foam::label Foam::blockDescriptor::nPoints() const
-{
-    return
-    (
-        (density_.x() + 1)
-      * (density_.y() + 1)
-      * (density_.z() + 1)
-    );
-}
-
-
-inline Foam::label Foam::blockDescriptor::nCells() const
-{
-    return
-    (
-        density_.x()
-      * density_.y()
-      * density_.z()
-    );
-}
-
-
-inline Foam::label Foam::blockDescriptor::nFaces() const
-{
-    return
-    (
-        ((density_.x() + 1) * density_.y() * density_.z())
-      + ((density_.y() + 1) * density_.z() * density_.x())
-      + ((density_.z() + 1) * density_.x() * density_.y())
-    );
-}
-
-
-inline Foam::label Foam::blockDescriptor::nInternalFaces() const
-{
-    return
-    (
-        ((density_.x() - 1) * density_.y() * density_.z())
-      + ((density_.y() - 1) * density_.z() * density_.x())
-      + ((density_.z() - 1) * density_.x() * density_.y())
-    );
-}
-
-
-inline Foam::label Foam::blockDescriptor::nBoundaryFaces() const
-{
-    return
-    (
-        (2 * density_.y() * density_.z())
-      + (2 * density_.z() * density_.x())
-      + (2 * density_.x() * density_.y())
-    );
-}
-
-
 inline const Foam::FixedList<Foam::label, 6>&
 Foam::blockDescriptor::curvedFaces() const
 {
@@ -138,63 +83,9 @@ inline const Foam::point& Foam::blockDescriptor::blockPoint(const label i) const
 }
 
 
-inline void Foam::blockDescriptor::checkIndex
-(
-    const label i, const label j, const label k
-) const
-{
-    if (i < 0 || i >= density_.x())
-    {
-        FatalErrorInFunction
-            << "index " << i << " out of range [0," << density_.x() << "]"
-            << abort(FatalError);
-    }
-    if (j < 0 || j >= density_.y())
-    {
-        FatalErrorInFunction
-            << "index " << j << " out of range [0," << density_.y() << "]"
-            << abort(FatalError);
-    }
-    if (k < 0 || k >= density_.z())
-    {
-        FatalErrorInFunction
-            << "index " << k << " out of range [0," << density_.z() << "]"
-            << abort(FatalError);
-    }
-}
-
-
-inline Foam::label Foam::blockDescriptor::cellLabel
-(
-    const label i,
-    const label j,
-    const label k
-) const
-{
-    #ifdef FULLDEBUG
-    checkIndex(i, j, k);
-    #endif
-    return (i + j*density_.x() + k*density_.x()*density_.y());
-}
-
-
-inline Foam::label Foam::blockDescriptor::pointLabel
-(
-    const label i,
-    const label j,
-    const label k
-) const
-{
-    #ifdef FULLDEBUG
-    checkIndex(i, j, k);
-    #endif
-    return (i + j*(density_.x()+1) + k*(density_.x()+1)*(density_.y()+1));
-}
-
-
 inline Foam::label Foam::blockDescriptor::facePointLabel
 (
-    const label facei,
+    const direction facei,
     const label i,
     const label j
 ) const
@@ -205,7 +96,7 @@ inline Foam::label Foam::blockDescriptor::facePointLabel
         return
         (
             i
-          + j*(density_.y() + 1)
+          + j*(sizes().y() + 1)
         );
     }
     else if (facei == 2 || facei == 3)
@@ -214,7 +105,7 @@ inline Foam::label Foam::blockDescriptor::facePointLabel
         return
         (
             i
-          + j*(density_.x() + 1)
+          + j*(sizes().x() + 1)
         );
     }
     else
@@ -223,7 +114,7 @@ inline Foam::label Foam::blockDescriptor::facePointLabel
         return
         (
             i
-          + j*(density_.x() + 1)
+          + j*(sizes().x() + 1)
         );
     }
 }
@@ -234,9 +125,9 @@ inline bool Foam::blockDescriptor::vertex
     const label i, const label j, const label k
 ) const
 {
-    bool iEnd = (i == 0 || i == density_.x());
-    bool jEnd = (j == 0 || j == density_.y());
-    bool kEnd = (k == 0 || k == density_.z());
+    bool iEnd = (i == 0 || i == sizes().x());
+    bool jEnd = (j == 0 || j == sizes().y());
+    bool kEnd = (k == 0 || k == sizes().z());
 
     return (iEnd && jEnd && kEnd);
 }
@@ -247,9 +138,9 @@ inline bool Foam::blockDescriptor::edge
     const label i, const label j, const label k
 ) const
 {
-    bool iEnd = (i == 0 || i == density_.x());
-    bool jEnd = (j == 0 || j == density_.y());
-    bool kEnd = (k == 0 || k == density_.z());
+    bool iEnd = (i == 0 || i == sizes().x());
+    bool jEnd = (j == 0 || j == sizes().y());
+    bool kEnd = (k == 0 || k == sizes().z());
 
     return (iEnd && jEnd) || (iEnd && kEnd) || (jEnd && kEnd);
 }
@@ -261,17 +152,13 @@ inline bool Foam::blockDescriptor::flatFaceOrEdge
 ) const
 {
     if (i == 0 && curvedFaces_[0] == -1) return true;
-    if (i == density_.x() && curvedFaces_[1] == -1) return true;
+    if (i == sizes().x() && curvedFaces_[1] == -1) return true;
     if (j == 0 && curvedFaces_[2] == -1) return true;
-    if (j == density_.y() && curvedFaces_[3] == -1) return true;
+    if (j == sizes().y() && curvedFaces_[3] == -1) return true;
     if (k == 0 && curvedFaces_[4] == -1) return true;
-    if (k == density_.z() && curvedFaces_[5] == -1) return true;
+    if (k == sizes().z() && curvedFaces_[5] == -1) return true;
 
-    bool iEnd = (i == 0 || i == density_.x());
-    bool jEnd = (j == 0 || j == density_.y());
-    bool kEnd = (k == 0 || k == density_.z());
-
-    return (iEnd && jEnd) || (iEnd && kEnd) || (jEnd && kEnd);
+    return this->edge(i, j, k);
 }
 
 
diff --git a/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C b/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C
index b774f55386..b26bff6645 100644
--- a/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C
+++ b/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C
@@ -28,11 +28,11 @@ License
 #include "lineDivide.H"
 #include "blockEdge.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    //- Calculate the geometric expension factor from the expansion ratio
+    //- Calculate the geometric expansion factor from the expansion ratio
     inline scalar calcGexp(const scalar expRatio, const label nDiv)
     {
         return nDiv > 1 ? pow(expRatio, 1.0/(nDiv - 1)) : 0.0;
@@ -96,7 +96,7 @@ Foam::lineDivide::lineDivide
             label secnEnd = secnStart + secnDiv;
 
             // Calculate the spacing
-            if (expRatio == 1.0)
+            if (equal(expRatio, 1))
             {
                 for (label i = secnStart; i < secnEnd; i++)
                 {
diff --git a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C
index 57a7b40e42..1c9619df0a 100644
--- a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C
+++ b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2016 OpenFOAM Foundation
@@ -74,7 +74,7 @@ Foam::label Foam::blockFaces::projectFace::index
 (
     const labelPair& n,
     const labelPair& coord
-) const
+)
 {
     return coord.first() + coord.second()*n.first();
 }
@@ -175,24 +175,24 @@ void Foam::blockFaces::projectFace::project
         case 0:
         case 1:
         {
-            n.first() = desc.density()[1] + 1;
-            n.second() = desc.density()[2] + 1;
+            n.first() = desc.density().y() + 1;
+            n.second() = desc.density().z() + 1;
         }
         break;
 
         case 2:
         case 3:
         {
-            n.first() = desc.density()[0] + 1;
-            n.second() = desc.density()[2] + 1;
+            n.first() = desc.density().x() + 1;
+            n.second() = desc.density().z() + 1;
         }
         break;
 
         case 4:
         case 5:
         {
-            n.first() = desc.density()[0] + 1;
-            n.second() = desc.density()[1] + 1;
+            n.first() = desc.density().x() + 1;
+            n.second() = desc.density().y() + 1;
         }
         break;
     }
diff --git a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H
index 1216204d11..b91d914684 100644
--- a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H
+++ b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H
@@ -70,11 +70,7 @@ class projectFace
         ) const;
 
         //- Convert i,j to single index
-        label index
-        (
-            const labelPair& n,
-            const labelPair& coord
-        ) const;
+        static label index(const labelPair& n, const labelPair& coord);
 
         //- Calulate lambdas (but unnormalised)
         void calcLambdas
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshCreate.C b/src/mesh/blockMesh/blockMesh/blockMeshCreate.C
index bfe61d6961..3e99d7f859 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshCreate.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshCreate.C
@@ -47,26 +47,28 @@ void Foam::blockMesh::createPoints() const
 
         if (verboseOutput)
         {
-            const Vector<label>& density = blocks[blocki].density();
+            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(density.x()-1, 0, 0);
-            label vin = blocks[blocki].pointLabel(density.x(), 0, 0);
+            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, density.y()-1, 0);
-            label vjn = blocks[blocki].pointLabel(0, density.y(), 0);
+            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, density.z()-1);
-            label vkn = blocks[blocki].pointLabel(0, 0, density.z());
+            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
diff --git a/src/mesh/blockMesh/blocks/block/block.C b/src/mesh/blockMesh/blocks/block/block.C
index 8bd8f373fa..0e88b8458d 100644
--- a/src/mesh/blockMesh/blocks/block/block.C
+++ b/src/mesh/blockMesh/blocks/block/block.C
@@ -43,7 +43,7 @@ Foam::block::block
     const pointField& vertices,
     const blockEdgeList& edges,
     const blockFaceList& faces,
-    const Vector<label>& density,
+    const labelVector& density,
     const UList<gradingDescriptors>& expand,
     const word& zoneName
 )
@@ -52,7 +52,11 @@ Foam::block::block
     points_(),
     blockCells_(),
     blockPatches_()
-{}
+{
+    // Always need points, and demand-driven data leaves dangling addressing?
+    createPoints();
+    createBoundary();
+}
 
 
 Foam::block::block
@@ -69,7 +73,11 @@ Foam::block::block
     points_(),
     blockCells_(),
     blockPatches_()
-{}
+{
+    // Always need points, and demand-driven data leaves dangling addressing?
+    createPoints();
+    createBoundary();
+}
 
 
 Foam::block::block(const blockDescriptor& blockDesc)
@@ -78,7 +86,11 @@ Foam::block::block(const blockDescriptor& blockDesc)
     points_(),
     blockCells_(),
     blockPatches_()
-{}
+{
+    // Always need points, and demand-driven data leaves dangling addressing?
+    createPoints();
+    createBoundary();
+}
 
 
 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
diff --git a/src/mesh/blockMesh/blocks/block/block.H b/src/mesh/blockMesh/blocks/block/block.H
index 7721b9c98b..f693a6f840 100644
--- a/src/mesh/blockMesh/blocks/block/block.H
+++ b/src/mesh/blockMesh/blocks/block/block.H
@@ -31,8 +31,7 @@ Description
     cells in each direction and an expansion ratio.
 
 Note
-    The vertices, cells, boundary patches for filling the block are
-    demand-driven.
+    The cells for filling the block are demand-driven.
 
 SourceFiles
     block.C
@@ -83,11 +82,12 @@ class block
         //- Create boundary patch faces for the block
         void createBoundary();
 
-        //- Add boundary faces to output list at iterator location
+        //- Add boundary faces for the shape face to the output list at
+        //- the iterator location
         template<class OutputIterator>
         OutputIterator addBoundaryFaces
         (
-            const int shapeFacei,
+            const direction shapeFacei,
             OutputIterator iter
         ) const;
 
@@ -133,7 +133,7 @@ public:
             const pointField& vertices,
             const blockEdgeList& edges,
             const blockFaceList& faces,
-            const Vector<label>& density,
+            const labelVector& density,
             const UList<gradingDescriptors>& expand,
             const word& zoneName = ""
         );
diff --git a/src/mesh/blockMesh/blocks/block/blockCreate.C b/src/mesh/blockMesh/blocks/block/blockCreate.C
index 0dd8bbccc0..dd8e20d73c 100644
--- a/src/mesh/blockMesh/blocks/block/blockCreate.C
+++ b/src/mesh/blockMesh/blocks/block/blockCreate.C
@@ -383,7 +383,7 @@ void Foam::block::createCells()
 template<class OutputIterator>
 OutputIterator Foam::block::addBoundaryFaces
 (
-    const int shapeFacei,
+    const direction shapeFacei,
     OutputIterator iter
 ) const
 {
@@ -524,7 +524,7 @@ void Foam::block::createBoundary()
     const label county = (density().z() * density().x());
     const label countz = (density().x() * density().y());
 
-    int patchi = 0;
+    direction patchi = 0;
 
     // 0 == x-min
     blockPatches_[patchi].resize(countx);
diff --git a/src/mesh/blockMesh/blocks/block/blockI.H b/src/mesh/blockMesh/blocks/block/blockI.H
index 51623c451b..514266c249 100644
--- a/src/mesh/blockMesh/blocks/block/blockI.H
+++ b/src/mesh/blockMesh/blocks/block/blockI.H
@@ -29,11 +29,6 @@ License
 
 inline const Foam::pointField& Foam::block::points() const
 {
-    if (points_.empty())
-    {
-        const_cast<block&>(*this).createPoints();
-    }
-
     return points_;
 }
 
@@ -53,11 +48,6 @@ Foam::block::cells() const
 inline const Foam::FixedList<Foam::List<Foam::FixedList<Foam::label, 4>>, 6>&
 Foam::block::boundaryPatches() const
 {
-    if (blockPatches_.empty())
-    {
-        const_cast<block&>(*this).createBoundary();
-    }
-
     return blockPatches_;
 }
 
diff --git a/src/mesh/blockMesh/gradingDescriptor/gradingDescriptor.C b/src/mesh/blockMesh/gradingDescriptor/gradingDescriptor.C
index 3ba189b571..2bde851f3b 100644
--- a/src/mesh/blockMesh/gradingDescriptor/gradingDescriptor.C
+++ b/src/mesh/blockMesh/gradingDescriptor/gradingDescriptor.C
@@ -49,7 +49,12 @@ Foam::gradingDescriptor::gradingDescriptor
     blockFraction_(blockFraction),
     nDivFraction_(nDivFraction),
     expansionRatio_(expansionRatio)
-{}
+{
+    if (expansionRatio_ < 0)
+    {
+        expansionRatio_ = 1.0/(-expansionRatio_);
+    }
+}
 
 
 Foam::gradingDescriptor::gradingDescriptor
@@ -60,7 +65,12 @@ Foam::gradingDescriptor::gradingDescriptor
     blockFraction_(1),
     nDivFraction_(1),
     expansionRatio_(expansionRatio)
-{}
+{
+    if (expansionRatio_ < 0)
+    {
+        expansionRatio_ = 1.0/(-expansionRatio_);
+    }
+}
 
 
 Foam::gradingDescriptor::gradingDescriptor(Istream& is)
diff --git a/src/mesh/blockMesh/gradingDescriptor/gradingDescriptor.H b/src/mesh/blockMesh/gradingDescriptor/gradingDescriptor.H
index 9dbc0d37f8..0cd0dce991 100644
--- a/src/mesh/blockMesh/gradingDescriptor/gradingDescriptor.H
+++ b/src/mesh/blockMesh/gradingDescriptor/gradingDescriptor.H
@@ -29,14 +29,17 @@ Class
 Description
     Handles the specification for grading within a section of a block
 
-    blockFraction: the fraction of the block the section occupies
+    The grading specification is handled is controlled by three parameters:
 
-    nDivFraction: the fraction of the divisions of the block allocated to
+    - blockFraction: the fraction of the block the section occupies
+
+    - nDivFraction: the fraction of the divisions of the block allocated to
         the section
 
-    expansionRatio: the expansions ratio for the grading with the section of
-        block defined as the ratio of the size of the division at either and
-        of the section.
+    - expansionRatio:
+       the expansion ratio for the grading with the section of
+       block defined as the ratio of end-size / start-size for the section.
+       A negative value is trapped and treated as its inverse.
 
 SourceFiles
     gradingDescriptor.C
@@ -53,12 +56,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declarations
 class Istream;
 class Ostream;
-
-// Forward declaration of friend functions and operators
 class gradingDescriptor;
 class gradingDescriptors;
+
 Istream& operator>>(Istream&, gradingDescriptor&);
 Ostream& operator<<(Ostream&, const gradingDescriptor&);