diff --git a/applications/test/PDRblockMesh/Make/files b/applications/test/PDRblockMesh/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..5ce1f00cc10e6cda50315a74e4e2602df1b27bc4
--- /dev/null
+++ b/applications/test/PDRblockMesh/Make/files
@@ -0,0 +1,3 @@
+Test-PDRblockMesh.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-PDRblockMesh
diff --git a/applications/test/PDRblockMesh/Make/options b/applications/test/PDRblockMesh/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..447e8d18fe81f46eb47f592fdf7b63ff9ded2265
--- /dev/null
+++ b/applications/test/PDRblockMesh/Make/options
@@ -0,0 +1,6 @@
+EXE_INC = \
+    -DFULLDEBUG \
+    -I$(LIB_SRC)/mesh/blockMesh/lnInclude
+
+EXE_LIBS = \
+    -lblockMesh
diff --git a/applications/test/PDRblockMesh/Test-PDRblockMesh.C b/applications/test/PDRblockMesh/Test-PDRblockMesh.C
new file mode 100644
index 0000000000000000000000000000000000000000..c505c069771f3c044b256707927eeca30f109163
--- /dev/null
+++ b/applications/test/PDRblockMesh/Test-PDRblockMesh.C
@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Description
+    Test accessors for PDRblock
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "PDRblock.H"
+
+using namespace Foam;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    argList::noFunctionObjects();
+    argList::addOption
+    (
+        "dict",
+        "file",
+        "Alternative dictionary for the PDRblockMesh description"
+    );
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    const word dictName("PDRblockMeshDict");
+
+    #include "setSystemRunTimeDictionaryIO.H"
+
+    Info<< "Reading " << dictName << nl << endl;
+
+    IOdictionary meshDict(dictIO);
+
+    PDRblock mesh(meshDict, true);
+
+    // Write summary
+    Info<< "----------------" << nl
+        << "Mesh Information" << nl
+        << "----------------" << nl
+        << "  " << "bounds: " << mesh.bounds() << nl
+        << "  " << "nPoints: " << mesh.nPoints()
+        << "  " << (mesh.sizes() + labelVector::one) << nl
+        << "  " << "nCells: " << mesh.nCells()
+        << "  " <<  mesh.sizes() << nl
+        << "  " << "nFaces: " << mesh.nFaces() << nl
+        << "  " << "nInternalFaces: " << mesh.nInternalFaces() << nl;
+
+    Info<< "----------------" << nl
+        << "Addressing" << nl
+        << "----------------" << nl
+        << "  " << "volume: " << mesh.bounds().volume() << nl
+        << "  " << "avg-vol: "<< (mesh.bounds().volume() / mesh.nCells()) << nl;
+
+
+    Info<< nl << "Check sizes: " << mesh.bounds().span() << nl;
+
+    for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+    {
+        // Three different ways of getting the same information
+        // not all are equally efficient
+
+        scalar totalWidth = 0;
+
+        // Using location
+        forAll(mesh.grid()[cmpt], i)
+        {
+            totalWidth += mesh.grid()[cmpt].width(i);
+        }
+
+        Info<< vector::componentNames[cmpt] << "   min/max "
+            << minMax(mesh.grid()[cmpt])
+            << ' ' << mesh.grid()[cmpt].first()
+            << ' ' << mesh.grid()[cmpt].last()
+            << "  " << mesh.grid()[cmpt].size() << " cells"  << nl;
+
+
+        // Use global (i,j,k) accessors, with mid-mesh for other directions
+
+        const label nx = mesh.sizes().x() / (cmpt == vector::X ? 1 : 2);
+        const label ny = mesh.sizes().y() / (cmpt == vector::Y ? 1 : 2);
+        const label nz = mesh.sizes().z() / (cmpt == vector::Z ? 1 : 2);
+
+
+        scalar totalSpan = 0;
+        scalar totalDelta = 0;
+
+        switch (cmpt)
+        {
+            case vector::X:
+            {
+                for (label i=0; i < nx; ++i)
+                {
+                    totalSpan += mesh.span(i, ny, nz).x();
+                    totalDelta += mesh.dx(i);
+                }
+            }
+            break;
+
+            case vector::Y:
+            {
+                for (label j=0; j < ny; ++j)
+                {
+                    totalSpan += mesh.span(nx, j, nz).y();
+                    totalDelta += mesh.dy(j);
+                }
+            }
+            break;
+
+            case vector::Z:
+            {
+                for (label k=0; k < nz; ++k)
+                {
+                    totalSpan += mesh.span(nx, ny, k).z();
+                    totalDelta += mesh.dz(k);
+                }
+            }
+            break;
+        }
+
+        Info<< "    width = " << totalWidth << " delta = " << totalDelta
+            << " span = " << totalSpan << nl;
+    }
+
+    {
+        const labelVector mid = mesh.sizes() / 2;
+
+        Info<< nl << "Mid-mesh information at " << mid << nl
+            << "  centre = " << mesh.C(mid) << nl
+            << "  volume = " << mesh.V(mid) << nl
+            << "  length = " << mesh.width(mid) << nl;
+    }
+
+    Info<< nl;
+
+    // Fatal with FULLDEBUG
+    {
+        const bool throwingError = FatalError.throwExceptions();
+
+        const label nx = mesh.sizes().x();
+        const label ny = mesh.sizes().y();
+        const label nz = mesh.sizes().z();
+
+        // Here we get error in x-y-z order
+        try
+        {
+            mesh.checkIndex(nx, ny, nz);
+        }
+        catch (const Foam::error& err)
+        {
+            Info<< nl << "Expected: Caught accesor error\n    "
+                << err.message().c_str() << nl;
+        }
+
+        // No order for the error since it is called parameter-wise
+        try
+        {
+            Info<< "centre at mesh(nx, ny, nz) = "
+                << mesh.C(mesh.sizes()) << nl;
+        }
+        catch (const Foam::error& err)
+        {
+            Info<< nl << "Expected: Caught accesor error\n    "
+                << err.message().c_str() << nl;
+        }
+
+        // Sizing error
+        try
+        {
+            mesh.checkSizes(labelVector(nx+1, ny, nz));
+        }
+        catch (const Foam::error& err)
+        {
+            Info<< nl << "Expected: Caught sizing error\n    "
+                << err.message().c_str() << nl;
+        }
+
+        try
+        {
+            mesh.checkSizes(labelMax/4);
+        }
+        catch (const Foam::error& err)
+        {
+            Info<< nl << "Expected: Caught sizing error\n    "
+                << err.message().c_str() << nl;
+        }
+
+        FatalError.throwExceptions(throwingError);
+    }
+
+
+    Info<< "\nEnd\n" << nl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/PDRblockMesh/box0/system/PDRblockMeshDict b/applications/test/PDRblockMesh/box0/system/PDRblockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..2a04fb12f6bc36aeba6640bd2e3645f9ef96a542
--- /dev/null
+++ b/applications/test/PDRblockMesh/box0/system/PDRblockMeshDict
@@ -0,0 +1,78 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      PDRblockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Point scaling (optional)
+scale 1.0;
+
+x
+{
+    points  ( -13.28 -0.10 6.0 19.19 );
+    nCells  ( 10 12 10 );
+    ratios  ( -5.16 1 5.16 );
+}
+
+y
+{
+    points  ( -12.98 0 5.50 18.48 );
+    nCells  ( 10 11 10 );
+    ratios  ( -5.16 1 5.16 );
+}
+
+z
+{
+    points  ( 0.00 4.80 17.26 );
+    nCells  ( 10 10 );
+    ratios  ( 1 5.16 );
+}
+
+
+defaultPatch
+{
+    name    walls;
+    type    wall;
+}
+
+
+// Faces: 0=x-min, 1=x-max, 2=y-min, 3=y-max, 4=z-min, 5=z-max
+
+boundary
+(
+    outer
+    {
+        type    patch;
+        faces   ( 0 1 2 3 5 );
+    }
+
+    mergingFaces
+    {
+        type    wall;
+        faces   ();
+    }
+
+    blockedFaces
+    {
+        type    wall;
+        faces   ();
+    }
+
+    ground
+    {
+        type    wall;
+        faces   ( 4 );
+    }
+);
+
+// ************************************************************************* //
diff --git a/applications/test/PDRblockMesh/box0/system/controlDict b/applications/test/PDRblockMesh/box0/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..a5a945dc8d0886e0443395f54d2efac7ffda5c65
--- /dev/null
+++ b/applications/test/PDRblockMesh/box0/system/controlDict
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+libs            ("libblockMesh.so");
+
+application     PDRblockMesh;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0;
+
+deltaT          0;
+
+writeControl    timeStep;
+
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  8;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //
diff --git a/applications/test/PDRblockMesh/box0/system/fvSchemes b/applications/test/PDRblockMesh/box0/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..7405277d729114929bb44f27b4b9cc6dfae9de0d
--- /dev/null
+++ b/applications/test/PDRblockMesh/box0/system/fvSchemes
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{}
+
+gradSchemes
+{}
+
+divSchemes
+{}
+
+laplacianSchemes
+{}
+
+interpolationSchemes
+{}
+
+snGradSchemes
+{}
+
+
+// ************************************************************************* //
diff --git a/applications/test/PDRblockMesh/box0/system/fvSolution b/applications/test/PDRblockMesh/box0/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..41d5821731215142aa28599c2f8fddc0c04bd836
--- /dev/null
+++ b/applications/test/PDRblockMesh/box0/system/fvSolution
@@ -0,0 +1,18 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/PDRblockMesh/Make/files b/applications/utilities/mesh/generation/PDRblockMesh/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..c945d7414e3bb0e4747339192b6a25c5afa39e1f
--- /dev/null
+++ b/applications/utilities/mesh/generation/PDRblockMesh/Make/files
@@ -0,0 +1,3 @@
+PDRblockMesh.C
+
+EXE = $(FOAM_APPBIN)/PDRblockMesh
diff --git a/applications/utilities/mesh/generation/PDRblockMesh/Make/options b/applications/utilities/mesh/generation/PDRblockMesh/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..b4f99e5295f698bac85899f3f379fe8b894ca1ef
--- /dev/null
+++ b/applications/utilities/mesh/generation/PDRblockMesh/Make/options
@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/mesh/blockMesh/lnInclude
+
+EXE_LIBS = \
+    -lblockMesh
diff --git a/applications/utilities/mesh/generation/PDRblockMesh/PDRblockMesh.C b/applications/utilities/mesh/generation/PDRblockMesh/PDRblockMesh.C
new file mode 100644
index 0000000000000000000000000000000000000000..60904946fc0c11ed195af91d78245d306a5d8160
--- /dev/null
+++ b/applications/utilities/mesh/generation/PDRblockMesh/PDRblockMesh.C
@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Application
+    PDRblockMesh
+
+Group
+    grpMeshGenerationUtilities
+
+Description
+    A specialized single-block mesh generator for a rectilinear mesh
+    in x-y-z.
+
+    Uses the mesh description found in
+      - \c system/PDRblockMeshDict
+
+Usage
+    \b PDRblockMesh [OPTION]
+
+    Options:
+      - \par -dict \<filename\>
+        Alternative dictionary for the mesh description.
+
+      - \par -noClean
+        Do not remove any existing polyMesh/ directory or files
+
+      - \par -time
+        Write resulting mesh to a time directory (instead of constant)
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "polyMesh.H"
+#include "PDRblock.H"
+#include "Time.H"
+#include "IOdictionary.H"
+#include "OSspecific.H"
+#include "OFstream.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "A block mesh generator for a rectilinear mesh in x-y-z.\n"
+        "  The ordering of vertex and face labels within a block as shown "
+        "below.\n"
+        "  For the local vertex numbering in the sequence 0 to 7:\n"
+        "    Faces 0, 1  ==  x-min, x-max.\n"
+        "    Faces 2, 3  ==  y-min, y-max.\n"
+        "    Faces 4, 5  ==  z-min, z-max.\n"
+        "\n"
+        "                        7 ---- 6\n"
+        "                 f5     |\\     |\\     f3\n"
+        "                 |      | 4 ---- 5     \\\n"
+        "                 |      3 |--- 2 |      \\\n"
+        "                 |       \\|     \\|      f2\n"
+        "                 f4       0 ---- 1\n"
+        "    Y  Z\n"
+        "     \\ |                f0 ------ f1\n"
+        "      \\|\n"
+        "       O--- X\n"
+    );
+
+    argList::noParallel();
+    argList::noFunctionObjects();
+
+    argList::addBoolOption
+    (
+        "noClean",
+        "Do not remove any existing polyMesh/ directory or files"
+    );
+    argList::addOption
+    (
+        "dict",
+        "file",
+        "Alternative dictionary for the PDRblockMesh description"
+    );
+    argList::addOption
+    (
+        "time",
+        "time",
+        "Specify a time to write mesh to (default: constant)"
+    );
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    // Instance for resulting mesh
+    bool useTime = false;
+    word meshInstance(runTime.constant());
+
+    if
+    (
+        args.readIfPresent("time", meshInstance)
+     && runTime.constant() != meshInstance
+    )
+    {
+        // Verify that the value is actually good
+        scalar timeValue;
+
+        useTime = readScalar(meshInstance, timeValue);
+        if (!useTime)
+        {
+            FatalErrorInFunction
+                << "Bad input value: " << meshInstance
+                << "Should be a scalar or 'constant'"
+                << nl << endl
+                << exit(FatalError);
+        }
+    }
+
+
+    const word dictName("PDRblockMeshDict");
+
+    #include "setSystemRunTimeDictionaryIO.H"
+
+    IOdictionary meshDict(dictIO);
+
+    Info<< "Creating PDRblockMesh from "
+        << runTime.relativePath(dictIO.objectPath()) << endl;
+
+    PDRblock blkMesh(meshDict, true);
+
+
+    // Instance for resulting mesh
+    if (useTime)
+    {
+        Info<< "Writing polyMesh to " << meshInstance << nl << endl;
+
+        // Make sure that the time is seen to be the current time.
+        // This is the logic inside regIOobject that resets the instance
+        // to the current time before writing
+        runTime.setTime(instant(meshInstance), 0);
+    }
+
+    if (!args.found("noClean"))
+    {
+        const fileName polyMeshPath
+        (
+            runTime.path()/meshInstance/polyMesh::meshSubDir
+        );
+
+        if (exists(polyMeshPath))
+        {
+            Info<< "Deleting polyMesh directory "
+                << runTime.relativePath(polyMeshPath) << endl;
+            rmDir(polyMeshPath);
+        }
+    }
+
+
+    Info<< nl << "Creating polyMesh from PDRblockMesh" << endl;
+
+    auto meshPtr = blkMesh.mesh
+    (
+        IOobject
+        (
+            polyMesh::defaultRegion,
+            meshInstance,
+            runTime,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        )
+    );
+
+
+    const polyMesh& mesh = *meshPtr;
+
+    // Set the precision of the points data to 10
+    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
+
+    Info<< nl << "Writing polyMesh with "
+        << mesh.cellZones().size() << " cellZones" << endl;
+
+    mesh.removeFiles();
+    if (!mesh.write())
+    {
+        FatalErrorInFunction
+            << "Failed writing polyMesh."
+            << exit(FatalError);
+    }
+
+    // Mesh summary
+    {
+        Info<< "----------------" << nl
+            << "Mesh Information" << nl
+            << "----------------" << nl
+            << "  " << "boundingBox: " << boundBox(mesh.points()) << nl
+            << "  " << "nPoints: " << mesh.nPoints() << nl
+            << "  " << "nCells: " << mesh.nCells() << nl
+            << "  " << "nFaces: " << mesh.nFaces() << nl
+            << "  " << "nInternalFaces: " << mesh.nInternalFaces() << nl;
+
+        Info<< "----------------" << nl
+            << "Patches" << nl
+            << "----------------" << nl;
+
+        for (const polyPatch& p : mesh.boundaryMesh())
+        {
+            Info<< "  " << "patch " << p.index()
+                << " (start: " << p.start()
+                << " size: " << p.size()
+                << ") name: " << p.name()
+                << nl;
+        }
+    }
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/blockMesh/blockMesh.C b/applications/utilities/mesh/generation/blockMesh/blockMesh.C
index 34fa38c9fd2e8dacabab4916226ad5fae35a8863..6498fab43450d36a1a446efd15ca2cb817507eb3 100644
--- a/applications/utilities/mesh/generation/blockMesh/blockMesh.C
+++ b/applications/utilities/mesh/generation/blockMesh/blockMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011, 2016-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -109,6 +109,8 @@ int main(int argc, char *argv[])
     );
 
     argList::noParallel();
+    argList::noFunctionObjects();
+
     argList::addBoolOption
     (
         "blockTopology",
@@ -134,7 +136,7 @@ int main(int argc, char *argv[])
     (
         "time",
         "time",
-        "specify a time to write mesh to"
+        "Specify a time to write mesh to (default: constant)"
     );
 
     #include "addRegionOption.H"
@@ -234,15 +236,15 @@ int main(int argc, char *argv[])
     {
         Info<< "Writing polyMesh to " << meshInstance << nl << endl;
 
-        // Make sure that the time is seen to be the current time. This
-        // is the logic inside regIOobject which resets the instance to the
-        // current time before writing
+        // Make sure that the time is seen to be the current time.
+        // This is the logic inside regIOobject that resets the instance
+        // to the current time before writing
         runTime.setTime(instant(meshInstance), 0);
     }
 
     if (!args.found("noClean"))
     {
-        fileName polyMeshPath
+        const fileName polyMeshPath
         (
             runTime.path()/meshInstance/regionPath/polyMesh::meshSubDir
         );
@@ -251,14 +253,14 @@ int main(int argc, char *argv[])
         {
             if (exists(polyMeshPath/dictName))
             {
-                Info<< "Not deleting polyMesh directory " << nl
-                    << "    " << polyMeshPath << nl
+                Info<< "Not deleting polyMesh directory "
+                    << runTime.relativePath(polyMeshPath) << nl
                     << "    because it contains " << dictName << endl;
             }
             else
             {
-                Info<< "Deleting polyMesh directory" << nl
-                    << "    " << polyMeshPath << endl;
+                Info<< "Deleting polyMesh directory "
+                    << runTime.relativePath(polyMeshPath) << endl;
                 rmDir(polyMeshPath);
             }
         }
@@ -267,9 +269,6 @@ int main(int argc, char *argv[])
 
     Info<< nl << "Creating polyMesh from blockMesh" << endl;
 
-    word defaultFacesName = "defaultFaces";
-    word defaultFacesType = emptyPolyPatch::typeName;
-
     polyMesh mesh
     (
         IOobject
@@ -283,8 +282,8 @@ int main(int argc, char *argv[])
         blocks.patches(),
         blocks.patchNames(),
         blocks.patchDicts(),
-        defaultFacesName,
-        defaultFacesType
+        "defaultFaces",               // Default patch name
+        emptyPolyPatch::typeName      // Default patch type
     );
 
 
@@ -297,9 +296,8 @@ int main(int argc, char *argv[])
 
     // Detect any cyclic patches and force re-ordering of the faces
     {
-        const polyPatchList& patches = mesh.boundaryMesh();
         bool hasCyclic = false;
-        for (const polyPatch& pp : patches)
+        for (const polyPatch& pp : mesh.boundaryMesh())
         {
             if (isA<cyclicPolyPatch>(pp))
             {
@@ -350,8 +348,6 @@ int main(int argc, char *argv[])
 
     // Write summary
     {
-        const polyPatchList& patches = mesh.boundaryMesh();
-
         Info<< "----------------" << nl
             << "Mesh Information" << nl
             << "----------------" << nl
@@ -365,7 +361,7 @@ int main(int argc, char *argv[])
             << "Patches" << nl
             << "----------------" << nl;
 
-        for (const polyPatch& p : patches)
+        for (const polyPatch& p : mesh.boundaryMesh())
         {
             Info<< "  " << "patch " << p.index()
                 << " (start: " << p.start()
diff --git a/applications/utilities/mesh/generation/blockMesh/findBlockMeshDict.H b/applications/utilities/mesh/generation/blockMesh/findBlockMeshDict.H
index c49b5e561dc72caf465c91312b666774a51fa83d..b05bda22c82e3261f7b69076b65d4231492524ff 100644
--- a/applications/utilities/mesh/generation/blockMesh/findBlockMeshDict.H
+++ b/applications/utilities/mesh/generation/blockMesh/findBlockMeshDict.H
@@ -24,7 +24,7 @@ autoPtr<IOdictionary> meshDictPtr;
         )
     )
     {
-        // Dictionary present constant polyMesh directory (old-style)
+        // Dictionary present in constant polyMesh directory (old-style)
 
         dictPath =
             runTime.constant()
@@ -63,8 +63,8 @@ autoPtr<IOdictionary> meshDictPtr;
             << exit(FatalError);
     }
 
-    Info<< "Creating block mesh from\n    "
-        << meshDictIO.objectPath() << endl;
+    Info<< "Creating block mesh from "
+        << runTime.relativePath(meshDictIO.objectPath()) << endl;
 
     meshDictPtr = autoPtr<IOdictionary>::New(meshDictIO);
 }
diff --git a/etc/config.sh/completion_cache b/etc/config.sh/completion_cache
index 4af50d7b6f5b31ee77b882dc0884abca7f09f6de..b105c588a324fd5dfba32417f6c73e82508d9c74 100644
--- a/etc/config.sh/completion_cache
+++ b/etc/config.sh/completion_cache
@@ -18,7 +18,7 @@ _of_complete_cache_[ansysToFoam]="-case -fileHandler -scale | -noFunctionObjects
 _of_complete_cache_[applyBoundaryLayer]="-Cbl -case -decomposeParDict -fileHandler -region -ybl | -noFunctionObjects -parallel -write-nut -doc -doc-source -help"
 _of_complete_cache_[attachMesh]="-case -fileHandler | -overwrite -doc -doc-source -help"
 _of_complete_cache_[autoPatch]="-case -fileHandler | -overwrite -doc -doc-source -help"
-_of_complete_cache_[blockMesh]="-case -dict -fileHandler -region -time | -blockTopology -noClean -noFunctionObjects -sets -doc -doc-source -help"
+_of_complete_cache_[blockMesh]="-case -dict -fileHandler -region -time | -blockTopology -noClean -sets -doc -doc-source -help"
 _of_complete_cache_[boundaryFoam]="-case -fileHandler | -dry-run -dry-run-write -listFunctionObjects -listFvOptions -listRegisteredSwitches -listScalarBCs -listSwitches -listTurbulenceModels -listUnsetSwitches -listVectorBCs -noFunctionObjects -doc -doc-source -help"
 _of_complete_cache_[boxTurb]="-case -fileHandler | -noFunctionObjects -doc -doc-source -help"
 _of_complete_cache_[buoyantBoussinesqPimpleFoam]="-case -decomposeParDict -fileHandler | -dry-run -dry-run-write -listFunctionObjects -listFvOptions -listRegisteredSwitches -listScalarBCs -listSwitches -listTurbulenceModels -listUnsetSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -doc -doc-source -help"
@@ -97,8 +97,7 @@ _of_complete_cache_[foamToGMV]="-case -decomposeParDict -fileHandler | -noFuncti
 _of_complete_cache_[foamToStarMesh]="-case -fileHandler -scale -time | -constant -latestTime -noBnd -noFunctionObjects -noZero -doc -doc-source -help"
 _of_complete_cache_[foamToSurface]="-case -fileHandler -scale -time | -constant -latestTime -noFunctionObjects -noZero -tri -doc -doc-source -help"
 _of_complete_cache_[foamToTetDualMesh]="-case -decomposeParDict -fileHandler -time | -constant -latestTime -noFunctionObjects -noZero -overwrite -parallel -doc -doc-source -help"
-_of_complete_cache_[foamToVTK]="-case -cellSet -cellZone -decomposeParDict -excludePatches -faceSet -faceZones -fields -fileHandler -name -patches -pointSet -region -regions -time | -allRegions -ascii -constant -finiteAreaFields -latestTime -legacy -nearCellValue -no-boundary -no-internal -no-lagrangian -no-point-data -noFunctionObjects -noZero -one-boundary -overwrite -parallel -poly-decomp -surfaceFields -doc -doc-source -help"
-_of_complete_cache_[foamToVTU]="-case -cellSet -decomposeParDict -excludePatches -faceSet -fields -fileHandler -region -time | -allPatches -ascii -constant -inline -latestTime -legacy -listFunctionObjects -listRegisteredSwitches -listScalarBCs -listSwitches -listUnsetSwitches -listVectorBCs -nearCellValue -noFaceZones -noFunctionObjects -noInternal -noLagrangian -noLinks -noPointZones -noZero -parallel -poly -surfaceFields -useTimeName -doc -doc-source -help"
+_of_complete_cache_[foamToVTK]="-case -cellSet -cellZone -decomposeParDict -excludePatches -faceSet -faceZones -fields -fileHandler -name -patches -pointSet -region -regions -time | -allRegions -ascii -constant -finiteAreaFields -latestTime -legacy -nearCellValue -no-boundary -no-internal -no-lagrangian -no-point-data -noFunctionObjects -noZero -one-boundary -overwrite -parallel -poly-decomp -processor-fields -surfaceFields -with-point-ids -doc -doc-source -help"
 _of_complete_cache_[foamUpgradeCyclics]="-case -decomposeParDict -fileHandler -region -time | -constant -dry-run -enableFunctionEntries -latestTime -noFunctionObjects -noZero -parallel -doc -doc-source -help"
 _of_complete_cache_[foamyHexMesh]="-case -decomposeParDict -fileHandler | -checkGeometry -conformationOnly -parallel -doc -doc-source -help"
 _of_complete_cache_[foamyQuadMesh]="-case -fileHandler -pointsFile | -noFunctionObjects -overwrite -doc -doc-source -help"
@@ -138,7 +137,7 @@ _of_complete_cache_[mergeOrSplitBaffles]="-case -decomposeParDict -dict -fileHan
 _of_complete_cache_[mergeSurfacePatches]="-case -fileHandler -output -patchIdRange -patchIds -patchNames | -keep -noFunctionObjects -doc -doc-source -help"
 _of_complete_cache_[meshToFPMA]="-case -decomposeParDict -fileHandler | -noFunctionObjects -parallel -doc -doc-source -help"
 _of_complete_cache_[mhdFoam]="-case -decomposeParDict -fileHandler | -dry-run -dry-run-write -listFunctionObjects -listRegisteredSwitches -listScalarBCs -listSwitches -listUnsetSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -doc -doc-source -help"
-_of_complete_cache_[mirrorMesh]="-case -dict -decomposeParDict -fileHandler | -noFunctionObjects -overwrite -parallel -doc -doc-source -help"
+_of_complete_cache_[mirrorMesh]="-case -decomposeParDict -dict -fileHandler | -noFunctionObjects -overwrite -parallel -doc -doc-source -help"
 _of_complete_cache_[mixtureAdiabaticFlameT]="-case -fileHandler | -doc -doc-source -help"
 _of_complete_cache_[modifyMesh]="-case -decomposeParDict -dict -fileHandler | -overwrite -parallel -doc -doc-source -help"
 _of_complete_cache_[moveDynamicMesh]="-case -decomposeParDict -fileHandler -region | -checkAMI -noFunctionObjects -parallel -doc -doc-source -help"
@@ -166,6 +165,7 @@ _of_complete_cache_[particleTracks]="-case -decomposeParDict -fileHandler -regio
 _of_complete_cache_[patchesToSubsets]="-case -fileHandler | -noFunctionObjects -doc -doc-source -help"
 _of_complete_cache_[patchSummary]="-case -decomposeParDict -fileHandler -region -time | -constant -expand -latestTime -noFunctionObjects -noZero -parallel -doc -doc-source -help"
 _of_complete_cache_[pdfPlot]="-case -decomposeParDict -fileHandler | -noFunctionObjects -parallel -doc -doc-source -help"
+_of_complete_cache_[PDRblockMesh]="-case -dict -fileHandler -time | -noClean -doc -doc-source -help"
 _of_complete_cache_[PDRFoam]="-case -decomposeParDict -fileHandler | -dry-run -dry-run-write -listFunctionObjects -listFvOptions -listRegisteredSwitches -listScalarBCs -listSwitches -listTurbulenceModels -listUnsetSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -doc -doc-source -help"
 _of_complete_cache_[PDRMesh]="-case -decomposeParDict -fileHandler | -overwrite -parallel -doc -doc-source -help"
 _of_complete_cache_[pimpleFoam]="-case -decomposeParDict -fileHandler | -dry-run -dry-run-write -listFunctionObjects -listFvOptions -listRegisteredSwitches -listScalarBCs -listSwitches -listTurbulenceModels -listUnsetSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -doc -doc-source -help"
@@ -221,7 +221,7 @@ _of_complete_cache_[simpleSprayFoam]="-case -decomposeParDict -fileHandler | -dr
 _of_complete_cache_[singleCellMesh]="-case -decomposeParDict -fileHandler -time | -constant -latestTime -noFunctionObjects -noZero -parallel -doc -doc-source -help"
 _of_complete_cache_[slopeMesh]="-case -decomposeParDict -fileHandler | -noFunctionObjects -parallel -doc -doc-source -help"
 _of_complete_cache_[smapToFoam]="-case -fileHandler | -noFunctionObjects -doc -doc-source -help"
-_of_complete_cache_[snappyHexMesh]="-case -decomposeParDict -dict -fileHandler -outFile -patches -region -surfaceSimplify | -checkGeometry -overwrite -parallel -profiling -doc -doc-source -help"
+_of_complete_cache_[snappyHexMesh]="-case -decomposeParDict -dict -fileHandler -outFile -patches -region -surfaceSimplify | -checkGeometry -dry-run -overwrite -parallel -profiling -doc -doc-source -help"
 _of_complete_cache_[snappyRefineMesh]="-case -fileHandler | -noFunctionObjects -doc -doc-source -help"
 _of_complete_cache_[solidDisplacementFoam]="-case -decomposeParDict -fileHandler | -dry-run -dry-run-write -listFunctionObjects -listRegisteredSwitches -listScalarBCs -listSwitches -listUnsetSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -doc -doc-source -help"
 _of_complete_cache_[solidEquilibriumDisplacementFoam]="-case -decomposeParDict -fileHandler | -dry-run -dry-run-write -listFunctionObjects -listRegisteredSwitches -listScalarBCs -listSwitches -listUnsetSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -doc -doc-source -help"
@@ -275,7 +275,6 @@ _of_complete_cache_[surfaceToPatch]="-case -faceSet -fileHandler -tol | -noFunct
 _of_complete_cache_[surfaceTransformPoints]="-case -fileHandler -origin -rollPitchYaw -rotate -rotate-angle -scale -translate -yawPitchRoll | -noFunctionObjects -doc -doc-source -help"
 _of_complete_cache_[surfactantFoam]="-case -decomposeParDict -fileHandler | -listFunctionObjects -listRegisteredSwitches -listScalarBCs -listSwitches -listUnsetSwitches -listVectorBCs -noFunctionObjects -parallel -doc -doc-source -help"
 _of_complete_cache_[temporalInterpolate]="-case -decomposeParDict -divisions -fields -fileHandler -interpolationType -region -time | -constant -latestTime -noZero -parallel -doc -doc-source -help"
-_of_complete_cache_[Test-decomposePar]="-case -decomposeParDict -domains -fileHandler -method -region -time | -allRegions -constant -latestTime -listFunctionObjects -listRegisteredSwitches -listScalarBCs -listSwitches -listUnsetSwitches -listVectorBCs -noFunctionObjects -noZero -verbose -doc -doc-source -help"
 _of_complete_cache_[tetgenToFoam]="-case -decomposeParDict -fileHandler | -noFaceFile -noFunctionObjects -parallel -doc -doc-source -help"
 _of_complete_cache_[tetMesh]="-case -decomposeParDict -fileHandler | -noFunctionObjects -parallel -doc -doc-source -help"
 _of_complete_cache_[thermoFoam]="-case -decomposeParDict -fileHandler | -dry-run -dry-run-write -listFunctionObjects -listFvOptions -listRegisteredSwitches -listScalarBCs -listSwitches -listTurbulenceModels -listUnsetSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -doc -doc-source -help"
diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H b/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H
index 3fb49e29820adda73feb08cdea58f887976e0180..feffeb6c17db9c62d1c99118bb0787d6b99e6b3b 100644
--- a/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H
+++ b/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H
@@ -90,24 +90,45 @@ public:
         //- Change the sizing parameters
         inline void reset(const label ni, const label nj, const label nk);
 
+        //- Change the sizing parameters
+        inline void reset(const labelVector& newSizes);
+
         //- Linear addressing index (offset) for an i,j,k position.
         inline label index(const label i, const label j, const label k) const;
 
+        //- Linear addressing index (offset) for an i,j,k position.
+        inline label index(const labelVector& ijk) const;
+
 
     // Checks
 
-        //- Check indices are within valid ni,nj,nk range.
+        //- Check indices are within 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 bool allowExtra = false
+        ) const;
+
+        //- Check indices are within ni,nj,nk range.
+        //  Optionally allow an extra index for point addressing
+        inline void checkIndex
+        (
+            const labelVector& ijk,
+            const bool allowExtra = false
         ) const;
 
-        //- Check that sizes() are valid
+        //- Check that all components of sizes() are non-negative
         inline void checkSizes() const;
+
+        //- Check that all components of sizes() match
+        inline void checkSizes(const labelVector& other) const;
+
+        //- Check that the total size matches
+        inline void checkSizes(const label nTotal) const;
+
 };
 
 
diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H b/src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H
index 3c8f442c8f81f9a322eda929121ccef5129b6ac7..9c526d7214154fbbbd40472bc807a60b7e53f786 100644
--- a/src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H
+++ b/src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H
@@ -72,7 +72,7 @@ inline Foam::labelVector& Foam::ijkAddressing::sizes()
 
 inline Foam::label Foam::ijkAddressing::size() const
 {
-    // Could also use  cmptProduct(sizes_);
+    // Could also use cmptProduct(sizes_);
     return (sizes_.x() * sizes_.y() * sizes_.z());
 }
 
@@ -106,6 +106,16 @@ inline void Foam::ijkAddressing::reset
 }
 
 
+inline void Foam::ijkAddressing::reset(const labelVector& newSizes)
+{
+    sizes_ = newSizes;
+
+    #ifdef FULLDEBUG
+    checkSizes();
+    #endif
+}
+
+
 inline Foam::label Foam::ijkAddressing::index
 (
     const label i,
@@ -121,15 +131,25 @@ inline Foam::label Foam::ijkAddressing::index
 }
 
 
+inline Foam::label Foam::ijkAddressing::index(const labelVector& ijk) const
+{
+    #ifdef FULLDEBUG
+    checkIndex(ijk);
+    #endif
+
+    return (ijk.x() + (sizes_.x() * (ijk.y() + (sizes_.y() * ijk.z()))));
+}
+
+
 inline void Foam::ijkAddressing::checkIndex
 (
     const label i,
     const label j,
     const label k,
-    const bool isPoint
+    const bool allowExtra
 ) const
 {
-    const label extra = (isPoint ? 1 : 0);
+    const label extra = (allowExtra ? 1 : 0);
 
     if (i < 0 || i >= (sizes_.x() + extra))
     {
@@ -155,6 +175,16 @@ inline void Foam::ijkAddressing::checkIndex
 }
 
 
+inline void Foam::ijkAddressing::checkIndex
+(
+    const labelVector& ijk,
+    const bool allowExtra
+) const
+{
+    checkIndex(ijk.x(), ijk.y(), ijk.z(), allowExtra);
+}
+
+
 inline void Foam::ijkAddressing::checkSizes() const
 {
     if (sizes_.x() < 0)
@@ -178,4 +208,28 @@ inline void Foam::ijkAddressing::checkSizes() const
 }
 
 
+inline void Foam::ijkAddressing::checkSizes(const labelVector& other) const
+{
+    if (sizes_ != other)
+    {
+        FatalErrorInFunction
+            << "The i-j-k sizes are different. "
+            << sizes_ << " vs. " << other << nl
+            << abort(FatalError);
+    }
+}
+
+
+inline void Foam::ijkAddressing::checkSizes(const label nTotal) const
+{
+    if (size() != nTotal)
+    {
+        FatalErrorInFunction
+            << "The total size is different. "
+            << size() << " vs. " << nTotal << nl
+            << abort(FatalError);
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkMesh.H b/src/OpenFOAM/meshes/ijkMesh/ijkMesh.H
index 976e3a6c9dac3a5f5b301eb57c4a3ed614241c34..783ddbfb4bc93b6d6e3bc56099457fce98f96e89 100644
--- a/src/OpenFOAM/meshes/ijkMesh/ijkMesh.H
+++ b/src/OpenFOAM/meshes/ijkMesh/ijkMesh.H
@@ -87,6 +87,10 @@ public:
         //- The number of boundary faces in the i-j-k mesh
         inline label nBoundaryFaces() const;
 
+        //- The number of boundary faces on the box face
+        //  Face: 0=x-min, 1=x-max, 2=y-min, 3=y-max, 4=z-min, 5=z-max
+        inline label nBoundaryFaces(const direction shapeFacei) const;
+
         //- The linear cell index for an i-j-k position - same as index()
         inline label cellLabel
         (
@@ -95,6 +99,9 @@ public:
             const label k
         ) const;
 
+        //- The linear cell index for an i-j-k position - same as index()
+        inline label cellLabel(const labelVector& ijk) const;
+
         //- The linear point index for an i-j-k position.
         //  Addressable in the range (ni+1, nj+1, nk+1).
         inline label pointLabel
@@ -103,6 +110,10 @@ public:
             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 labelVector& ijk) const;
 };
 
 
diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H b/src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H
index cc89c2c520e3dc4bb9b0b3f60cb49be4eaedcd32..da73e24be321aa7535de6b784de6b8ec29bb3cf7 100644
--- a/src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H
+++ b/src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H
@@ -123,6 +123,43 @@ inline Foam::label Foam::ijkMesh::nBoundaryFaces() const
 }
 
 
+inline Foam::label Foam::ijkMesh::nBoundaryFaces
+(
+    const direction shapeFacei
+) const
+{
+    if (ijkAddressing::empty())
+    {
+        return 0;
+    }
+
+    const labelVector& n = ijkAddressing::sizes();
+
+    switch (shapeFacei)
+    {
+        // Face 0,1 == x-min, x-max
+        case 0:
+        case 1:
+            return n.y()*n.z();
+            break;
+
+        // Face 2, 3 == y-min, y-max
+        case 2:
+        case 3:
+            return n.z()*n.x();
+            break;
+
+        // Face 4, 5 == z-min, z-max
+        case 4:
+        case 5:
+            return n.x()*n.y();
+            break;
+    }
+
+    return 0;
+}
+
+
 inline Foam::label Foam::ijkMesh::cellLabel
 (
     const label i,
@@ -134,6 +171,12 @@ inline Foam::label Foam::ijkMesh::cellLabel
 }
 
 
+inline Foam::label Foam::ijkMesh::cellLabel(const labelVector& ijk) const
+{
+    return ijkAddressing::index(ijk);
+}
+
+
 inline Foam::label Foam::ijkMesh::pointLabel
 (
     const label i,
@@ -151,4 +194,16 @@ inline Foam::label Foam::ijkMesh::pointLabel
 }
 
 
+Foam::label Foam::ijkMesh::pointLabel(const labelVector& ijk) const
+{
+    #ifdef FULLDEBUG
+    checkIndex(ijk, true);
+    #endif
+
+    const labelVector& n = sizes();
+
+    return (ijk.x() + ((n.x()+1) * (ijk.y() + (n.y()+1) * ijk.z())));
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/blockMesh/Make/files b/src/mesh/blockMesh/Make/files
index 7d567453f887fd53f1b77e2bb2192f5e723e7aaf..84473bfc0f232e12373adf983eb1a3f4c3bdfac8 100644
--- a/src/mesh/blockMesh/Make/files
+++ b/src/mesh/blockMesh/Make/files
@@ -39,4 +39,7 @@ blockMesh/blockMeshMergeFast.C
 
 blockMeshTools/blockMeshTools.C
 
+PDRblockMesh/PDRblock.C
+PDRblockMesh/PDRblockCreate.C
+
 LIB = $(FOAM_LIBBIN)/libblockMesh
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblock.C b/src/mesh/blockMesh/PDRblockMesh/PDRblock.C
new file mode 100644
index 0000000000000000000000000000000000000000..0ae9bab0bb2d8218d35b3f9f431208b03127030f
--- /dev/null
+++ b/src/mesh/blockMesh/PDRblockMesh/PDRblock.C
@@ -0,0 +1,440 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRblock.H"
+#include "emptyPolyPatch.H"
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    //- Calculate geometric expansion factor from expansion ratio
+    inline scalar calcGexp(const scalar expRatio, const label nDiv)
+    {
+        return nDiv > 1 ? pow(expRatio, 1.0/(nDiv - 1)) : 0.0;
+    }
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::PDRblock::readGridControl
+(
+    const direction cmpt,
+    const dictionary& dict,
+    const scalar scaleFactor
+)
+{
+    //- The begin/end nodes for each segment
+    scalarList knots;
+
+    // The number of division per segment
+    labelList divisions;
+
+    // The expansion ratio per segment
+    scalarList expansion;
+
+    if (verbose_)
+    {
+        Info<< "Reading grid control for "
+            << vector::componentNames[cmpt] << " direction" << nl;
+    }
+
+    // Points
+    // ~~~~~~
+
+    dict.readEntry("points", knots);
+
+    const label nSegments = (knots.size()-1);
+
+    if (nSegments < 1)
+    {
+        FatalIOErrorInFunction(dict)
+            << "Must define at least two control points:"
+            << " in " << vector::componentNames[cmpt]
+            << " direction" << nl
+            << exit(FatalIOError);
+    }
+
+    // Apply scaling
+    if (scaleFactor > 0)
+    {
+        for (scalar& pt : knots)
+        {
+            pt *= scaleFactor;
+        }
+    }
+
+    // Do points increase monotonically?
+    {
+        const scalar& minVal = knots.first();
+
+        for (label pointi = 1; pointi < knots.size(); ++pointi)
+        {
+            if (knots[pointi] <= minVal)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "Points are not monotonically increasing:"
+                    << " in " << vector::componentNames[cmpt]
+                    << " direction" << nl
+                    << exit(FatalIOError);
+            }
+        }
+    }
+
+
+    if (verbose_)
+    {
+        Info<< "    points : " << flatOutput(knots) << nl;
+    }
+
+
+    // Divisions
+    // ~~~~~~~~~
+
+    dict.readEntry("nCells", divisions);
+
+    label nTotalDivs = 0;
+    for (const label ndiv : divisions)
+    {
+        nTotalDivs += ndiv;
+
+        if (ndiv < 1)
+        {
+            FatalIOErrorInFunction(dict)
+                << "Negative or zero divisions:"
+                << " in " << vector::componentNames[cmpt]
+                << " direction" << nl
+                << exit(FatalIOError);
+        }
+    }
+
+    if (divisions.size() != nSegments)
+    {
+        FatalIOErrorInFunction(dict)
+            << "Mismatch in number of divisions and number of points:"
+            << " in " << vector::componentNames[cmpt]
+            << " direction" << nl
+            << exit(FatalIOError);
+    }
+    else if (!nTotalDivs)
+    {
+        FatalIOErrorInFunction(dict)
+            << "No grid divisions at all:"
+            << " in " << vector::componentNames[cmpt]
+            << " direction" << nl
+            << exit(FatalIOError);
+    }
+
+
+    if (verbose_)
+    {
+        Info<< "    nCells : " << flatOutput(divisions) << nl;
+    }
+
+
+    // Expansion ratios
+    // ~~~~~~~~~~~~~~~~
+
+    dict.readIfPresent("ratios", expansion);
+
+    // Correct expansion ratios - negative is the same as inverse.
+    for (scalar& expRatio : expansion)
+    {
+        if (expRatio < 0)
+        {
+            expRatio = 1.0/(-expRatio);
+        }
+    }
+
+    if (expansion.size() && expansion.size() != nSegments)
+    {
+        FatalIOErrorInFunction(dict)
+            << "Mismatch in number of expansion ratios and divisions:"
+            << " in " << vector::componentNames[cmpt]
+            << " direction" << nl
+            << exit(FatalIOError);
+    }
+
+    if (expansion.empty())
+    {
+        expansion.resize(nSegments, scalar(1));
+
+        if (verbose_)
+        {
+            Info<< "Warning: 'ratios' not specified, using constant spacing"
+                << nl;
+        }
+    }
+    else
+    {
+        if (verbose_)
+        {
+            Info<< "    ratios : " << flatOutput(expansion) << nl;
+        }
+    }
+
+
+
+    // Define the grid points
+
+    scalarList& gridPoint = grid_[cmpt];
+    gridPoint.resize(nTotalDivs+1);
+
+    label start = 0;
+
+    for (label segmenti=0; segmenti < nSegments; ++segmenti)
+    {
+        const label nDiv = divisions[segmenti];
+        const scalar expRatio = expansion[segmenti];
+
+        SubList<scalar> subPoint(gridPoint, nDiv+1, start);
+        start += nDiv;
+
+        subPoint[0] = knots[segmenti];
+        subPoint[nDiv] = knots[segmenti+1];
+
+        const scalar dist = (subPoint.last() - subPoint.first());
+
+        if (equal(expRatio, 1.0))
+        {
+            for (label i=1; i < nDiv; ++i)
+            {
+                subPoint[i] = (subPoint[0] + (dist * i)/nDiv);
+            }
+        }
+        else
+        {
+            // Geometric expansion factor from the expansion ratio
+            const scalar expFact = calcGexp(expRatio, nDiv);
+
+            for (label i=1; i < nDiv; ++i)
+            {
+                subPoint[i] =
+                (
+                    subPoint[0]
+                  + dist * (1.0 - pow(expFact, i))/(1.0 - pow(expFact, nDiv))
+                );
+            }
+        }
+    }
+}
+
+
+void Foam::PDRblock::readBoundary(const dictionary& dict)
+{
+    if (verbose_)
+    {
+        Info<< "Reading boundary entries" << nl;
+    }
+
+    PtrList<entry> patchEntries;
+
+    if (dict.found("boundary"))
+    {
+        PtrList<entry> newEntries(dict.lookup("boundary"));
+        patchEntries.transfer(newEntries);
+    }
+
+    patches_.clear();
+    patches_.resize(patchEntries.size());
+
+    // Hex cell has 6 sides
+    const labelRange validFaceId(0, 6);
+
+    // Track which sides have been handled
+    labelHashSet handled;
+
+    forAll(patchEntries, patchi)
+    {
+        if (!patchEntries.set(patchi))
+        {
+            continue;
+        }
+
+        const entry& e = patchEntries[patchi];
+
+        if (!e.isDict())
+        {
+            FatalIOErrorInFunction(dict)
+                << "Entry " << e
+                << " in boundary section is not a valid dictionary."
+                << exit(FatalIOError);
+        }
+
+        const dictionary& patchDict = e.dict();
+
+        const word& patchName = e.keyword();
+
+        const word patchType(patchDict.get<word>("type"));
+
+        labelList faceLabels(patchDict.get<labelList>("faces"));
+
+        labelHashSet patchFaces(faceLabels);
+
+        if (faceLabels.size() != patchFaces.size())
+        {
+            WarningInFunction
+                << "Patch: " << patchName
+                << ": Ignoring duplicate face ids" << nl;
+        }
+
+        label nErased = patchFaces.filterKeys(validFaceId);
+        if (nErased)
+        {
+            WarningInFunction
+                << "Patch: " << patchName << ": Ignoring " << nErased
+                << " faces with invalid identifiers" << nl;
+        }
+
+        nErased = patchFaces.erase(handled);
+        if (nErased)
+        {
+            WarningInFunction
+                << "Patch: " << patchName << ": Ignoring " << nErased
+                << " faces, which were already handled" << nl;
+        }
+
+        // Mark these faces as having been handled
+        handled += patchFaces;
+
+
+        // Save information for later access during mesh creation.
+
+        patches_.set(patchi, new boundaryEntry());
+
+        boundaryEntry& bentry = patches_[patchi];
+
+        bentry.name_ = patchName;
+        bentry.type_ = patchType;
+        bentry.size_ = 0;
+        bentry.faces_ = patchFaces.sortedToc();
+
+        // Count patch faces
+        for (const label shapeFacei : bentry.faces_)
+        {
+            bentry.size_ += nBoundaryFaces(shapeFacei);
+        }
+    }
+
+
+    // Deal with unhandled faces
+
+    labelHashSet missed(identity(6));
+    missed.erase(handled);
+
+    if (missed.size())
+    {
+        patches_.append(new boundaryEntry());
+
+        boundaryEntry& bentry = patches_.last();
+
+        bentry.name_ = "defaultFaces";
+        bentry.type_ = emptyPolyPatch::typeName;
+        bentry.size_ = 0;
+        bentry.faces_ = missed.sortedToc();
+
+        // Count patch faces
+        for (const label shapeFacei : bentry.faces_)
+        {
+            bentry.size_ += nBoundaryFaces(shapeFacei);
+        }
+
+
+        // Use name/type from "defaultPatch" entry if available
+        // - be generous with error handling
+
+        const dictionary* defaultEntry = dict.findDict("defaultPatch");
+        if (defaultEntry)
+        {
+            defaultEntry->readIfPresent("name", bentry.name_);
+            defaultEntry->readIfPresent("type", bentry.type_);
+        }
+    }
+
+    if (verbose_)
+    {
+        label patchi = 0;
+        for (const boundaryEntry& bentry : patches_)
+        {
+            Info<< "    patch: " << patchi
+                << " (size: " << bentry.size_
+                << " type: " << bentry.type_
+                << ") name: " << bentry.name_
+                << " faces: " << flatOutput(bentry.faces_) << nl;
+
+            ++patchi;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::PDRblock::PDRblock(const dictionary& dict, bool verbose)
+:
+    PDRblock()
+{
+    verbose_ = verbose;
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::PDRblock::read(const dictionary& dict)
+{
+    // Mark no scaling with invalid value
+    const scalar scaleFactor = dict.lookupOrDefault<scalar>("scale", -1);
+
+    // Grid controls
+    readGridControl(0, dict.subDict("x"), scaleFactor);
+    readGridControl(1, dict.subDict("y"), scaleFactor);
+    readGridControl(2, dict.subDict("z"), scaleFactor);
+
+    // Adjust i-j-k addressing
+    sizes().x() = grid_.x().nCells();
+    sizes().y() = grid_.y().nCells();
+    sizes().z() = grid_.z().nCells();
+
+    // Adjust boundBox
+    bounds_.min().x() = grid_.x().first();
+    bounds_.min().y() = grid_.y().first();
+    bounds_.min().z() = grid_.z().first();
+
+    bounds_.max().x() = grid_.x().last();
+    bounds_.max().y() = grid_.y().last();
+    bounds_.max().z() = grid_.z().last();
+
+
+    // Boundaries
+    readBoundary(dict);
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblock.H b/src/mesh/blockMesh/PDRblockMesh/PDRblock.H
new file mode 100644
index 0000000000000000000000000000000000000000..419d23f0d1d8ed168076c62db4a22c5c804d33fc
--- /dev/null
+++ b/src/mesh/blockMesh/PDRblockMesh/PDRblock.H
@@ -0,0 +1,289 @@
+/*---------------------------------------------------------------------------* \
+  =========                 |
+  \\      /  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::PDRblock
+
+Description
+    A single block x-y-z rectilinear mesh addressable as i,j,k with
+    simplified creation. Some of the input is similar to blockMeshDict,
+    but since this specialization is for a single-block that is aligned
+    with the x-y-z directions, it provides a different means of specifying
+    the mesh.
+
+    Dictionary controls
+    \table
+        Property    | Description                          | Required | Default
+        x           | X-direction grid specification       | yes |
+        y           | Y-direction grid specification       | yes |
+        z           | Z-direction grid specification       | yes |
+        scale       | Point scaling                        | no  | 1.0
+        boundary    | Boundary patches                     | yes |
+        defaultPatch | Default patch specification         | no  |
+    \endtable
+
+    Grid coordinate controls
+    \table
+        Property| Description                          | Required | Default
+        points  | Locations defining the mesh segment  | yes |
+        nCells  | Divisions per mesh segment           | yes |
+        ratios  | Expansion ratio (end/start) per segment | no | uniform
+    \endtable
+
+    The expansion ratios are defined as per blockMesh and represent the ratio
+    of end-size / start-size for the section. A negative value is trapped
+    and treated as its reciprocal.
+
+SourceFiles
+    PDRblockI.H
+    PDRblock.C
+    PDRblockCreate.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRblock_H
+#define PDRblock_H
+
+#include "ijkMesh.H"
+#include "boundBox.H"
+#include "pointField.H"
+#include "faceList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class IOobject;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class PDRblock Declaration
+\*---------------------------------------------------------------------------*/
+
+class PDRblock
+:
+    public ijkMesh
+{
+
+public:
+
+    // Public Classes
+
+        //- Grid locations in an axis direction.
+        //  The number of points is one larger than the number of elements
+        //  it represents
+        struct location
+        :
+            public scalarList
+        {
+            //- The number of cells in this direction.
+            inline label nCells() const;
+
+            //- The number of points in this direction.
+            inline label nPoints() const;
+
+            //- The number of divisions (nCells) in this direction.
+            inline label size() const;
+
+            //- Check that element index is within valid range.
+            inline void checkIndex(const label i) const;
+
+            //- Cell size at element position.
+            inline scalar width(const label i) const;
+
+            //- Cell centre at element position.
+            inline scalar C(const label i) const;
+        };
+
+
+private:
+
+    // Private Classes
+
+        //- Extracted patch settings
+        struct boundaryEntry
+        {
+            //- The patch name
+            word name_;
+
+            //- The patch type
+            word type_;
+
+            //- The patch size
+            label size_;
+
+            //- The associated block face ids [0,5]
+            labelList faces_;
+        };
+
+
+    // Private Data
+
+        //- Verbosity
+        bool verbose_;
+
+        //- The grid points in all (i,j,k / x,y,z) directions.
+        Vector<location> grid_;
+
+        //- The mesh bounding box
+        boundBox bounds_;
+
+        //- The boundary patch information
+        PtrList<boundaryEntry> patches_;
+
+
+    // Private Member Functions
+
+        //- Read and define grid points in given direction
+        void readGridControl
+        (
+            const direction cmpt,
+            const dictionary& dict,
+            const scalar scaleFactor = -1
+        );
+
+        //- Read "boundary" information
+        void readBoundary(const dictionary& dict);
+
+
+        //- Populate point field for the block
+        void createPoints(pointField& pts) const;
+
+        //- Add internal faces to lists.
+        //  Lists must be properly sized!
+        //  \return the number of faces added
+        label addInternalFaces
+        (
+            faceList::iterator& faceIter,
+            labelList::iterator& ownIter,
+            labelList::iterator& neiIter
+        ) const;
+
+
+        //- Add boundary faces for the shape face to lists
+        //  Lists must be properly sized!
+        //  \return the number of faces added
+        label addBoundaryFaces
+        (
+            const direction shapeFacei,
+            faceList::iterator& faceIter,
+            labelList::iterator& ownIter
+        ) const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct zero-size
+        inline PDRblock();
+
+        //- Construct from dictionary
+        explicit PDRblock(const dictionary& dict, bool verbose=false);
+
+
+    // Member Functions
+
+        //- Read dictionary
+        bool read(const dictionary& dict);
+
+
+    // Access
+
+        //- The grid point locations in the i,j,k (x,y,z) directions.
+        inline const Vector<location>& grid() const;
+
+
+    // Mesh information
+
+        //- The mesh bounding box
+        inline const boundBox& bounds() const;
+
+        //- Cell size in x-direction at i position.
+        inline scalar dx(const label i) const;
+
+        //- Cell size in x-direction at i position.
+        inline scalar dx(const labelVector& ijk) const;
+
+        //- Cell size in y-direction at j position.
+        inline scalar dy(const label j) const;
+
+        //- Cell size in y-direction at j position.
+        inline scalar dy(const labelVector& ijk) const;
+
+        //- Cell size in z-direction at k position.
+        inline scalar dz(const label k) const;
+
+        //- Cell size in z-direction at k position.
+        inline scalar dz(const labelVector& ijk) const;
+
+        //- Cell dimensions at i,j,k position.
+        inline vector span(const label i, const label j, const label k) const;
+
+        //- Cell dimensions at i,j,k position.
+        inline vector span(const labelVector& ijk) const;
+
+        //- Cell centre at i,j,k position.
+        inline point C(const label i, const label j, const label k) const;
+
+        //- Cell centre at i,j,k position.
+        inline point C(const labelVector& ijk) const;
+
+        //- Cell volume at i,j,k position.
+        inline scalar V(const label i, const label j, const label k) const;
+
+        //- Cell volume at i,j,k position.
+        inline scalar V(const labelVector& ijk) const;
+
+        //- Characteristic cell size at i,j,k position.
+        //  This is the cubic root of the volume
+        inline scalar width(const label i, const label j, const label k) const;
+
+        //- Characteristic cell size at i,j,k position.
+        //  This is the cubic root of the volume
+        inline scalar width(const labelVector& ijk) const;
+
+
+    // Mesh generation
+
+        //- Create polyMesh for grid definition and patch information
+        autoPtr<polyMesh> mesh(const IOobject& io) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "PDRblockI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockCreate.C b/src/mesh/blockMesh/PDRblockMesh/PDRblockCreate.C
new file mode 100644
index 0000000000000000000000000000000000000000..d50e94f78bdcfce5bf7a6f5ee25c250af3ce7c14
--- /dev/null
+++ b/src/mesh/blockMesh/PDRblockMesh/PDRblockCreate.C
@@ -0,0 +1,372 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRblock.H"
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::PDRblock::createPoints(pointField& pts) const
+{
+    const label ni = sizes().x();
+    const label nj = sizes().y();
+    const label nk = sizes().z();
+
+    pts.resize(nPoints());
+
+    for (label k=0; k<=nk; ++k)
+    {
+        for (label j=0; j<=nj; ++j)
+        {
+            for (label i=0; i<=ni; ++i)
+            {
+                point& pt = pts[pointLabel(i,j,k)];
+
+                pt.x() = grid_.x()[i];
+                pt.y() = grid_.y()[j];
+                pt.z() = grid_.z()[k];
+            }
+        }
+    }
+}
+
+
+Foam::label Foam::PDRblock::addInternalFaces
+(
+    faceList::iterator& faceIter,
+    labelList::iterator& ownIter,
+    labelList::iterator& neiIter
+) const
+{
+    const label ni = sizes().x();
+    const label nj = sizes().y();
+    const label nk = sizes().z();
+
+    const labelList::iterator firstIter = ownIter;
+
+    for (label k=0; k<nk; ++k)
+    {
+        for (label j=0; j<nj; ++j)
+        {
+            for (label i=0; i<ni; ++i)
+            {
+                const label celli = cellLabel(i, j, k);
+
+                // Local Face 1 == x-max
+                if (i < ni-1)
+                {
+                    auto& f = *faceIter;
+                    ++faceIter;
+                    f.resize(4);
+
+                    f[0] = pointLabel(i+1, j,   k);
+                    f[1] = pointLabel(i+1, j+1, k);
+                    f[2] = pointLabel(i+1, j+1, k+1);
+                    f[3] = pointLabel(i+1, j,   k+1);
+
+                    *ownIter = celli;
+                    *neiIter = cellLabel(i+1, j, k);
+
+                    ++ownIter;
+                    ++neiIter;
+                }
+
+                // Local Face 3 == y-max
+                if (j < nj-1)
+                {
+                    auto& f = *faceIter;
+                    ++faceIter;
+                    f.resize(4);
+
+                    f[0] = pointLabel(i,   j+1, k);
+                    f[1] = pointLabel(i,   j+1, k+1);
+                    f[2] = pointLabel(i+1, j+1, k+1);
+                    f[3] = pointLabel(i+1, j+1, k);
+
+                    *ownIter = celli;
+                    *neiIter = cellLabel(i, j+1, k);
+
+                    ++ownIter;
+                    ++neiIter;
+                }
+
+                // Local Face 5 == z-max
+                if (k < nk-1)
+                {
+                    auto& f = *faceIter;
+                    ++faceIter;
+                    f.resize(4);
+
+                    f[0] = pointLabel(i,   j,   k+1);
+                    f[1] = pointLabel(i+1, j,   k+1);
+                    f[2] = pointLabel(i+1, j+1, k+1);
+                    f[3] = pointLabel(i,   j+1, k+1);
+
+                    *ownIter = celli;
+                    *neiIter = cellLabel(i, j, k+1);
+
+                    ++ownIter;
+                    ++neiIter;
+                }
+            }
+        }
+    }
+
+    // Return the number of faces added
+    return (ownIter - firstIter);
+}
+
+
+Foam::label Foam::PDRblock::addBoundaryFaces
+(
+    const direction shapeFacei,
+    faceList::iterator& faceIter,
+    labelList::iterator& ownIter
+) const
+{
+    const label ni = sizes().x();
+    const label nj = sizes().y();
+    const label nk = sizes().z();
+
+    const labelList::iterator firstIter = ownIter;
+
+    switch (shapeFacei)
+    {
+        // Face 0 == x-min
+        case 0:
+        {
+            for (label k=0; k<nk; ++k)
+            {
+                for (label j=0; j<nj; ++j)
+                {
+                    auto& f = *faceIter;
+                    ++faceIter;
+                    f.resize(4);
+
+                    f[0] = pointLabel(0, j,   k);
+                    f[1] = pointLabel(0, j,   k+1);
+                    f[2] = pointLabel(0, j+1, k+1);
+                    f[3] = pointLabel(0, j+1, k);
+
+                    *ownIter = cellLabel(0, j, k);
+                    ++ownIter;
+                }
+            }
+        }
+        break;
+
+        // Face 1 == x-max
+        case 1:
+        {
+            for (label k=0; k<nk; ++k)
+            {
+                for (label j=0; j<nj; ++j)
+                {
+                    auto& f = *faceIter;
+                    ++faceIter;
+                    f.resize(4);
+
+                    f[0] = pointLabel(ni, j,   k);
+                    f[1] = pointLabel(ni, j+1, k);
+                    f[2] = pointLabel(ni, j+1, k+1);
+                    f[3] = pointLabel(ni, j,   k+1);
+
+                    *ownIter = cellLabel(ni-1, j, k);
+                    ++ownIter;
+                }
+            }
+        }
+        break;
+
+        // Face 2 == y-min
+        case 2:
+        {
+            for (label i=0; i<ni; ++i)
+            {
+                for (label k=0; k<nk; ++k)
+                {
+                    auto& f = *faceIter;
+                    ++faceIter;
+                    f.resize(4);
+
+                    f[0] = pointLabel(i,   0, k);
+                    f[1] = pointLabel(i+1, 0, k);
+                    f[2] = pointLabel(i+1, 0, k+1);
+                    f[3] = pointLabel(i,   0, k+1);
+
+                    *ownIter = cellLabel(i, 0, k);
+                    ++ownIter;
+                }
+            }
+        }
+        break;
+
+        // Face 3 == y-max
+        case 3:
+        {
+            for (label i=0; i<ni; ++i)
+            {
+                for (label k=0; k<nk; ++k)
+                {
+                    auto& f = *faceIter;
+                    ++faceIter;
+                    f.resize(4);
+
+                    f[0] = pointLabel(i,   nj, k);
+                    f[1] = pointLabel(i,   nj, k+1);
+                    f[2] = pointLabel(i+1, nj, k+1);
+                    f[3] = pointLabel(i+1, nj, k);
+
+                    *ownIter = cellLabel(i, nj-1, k);
+                    ++ownIter;
+                }
+            }
+        }
+        break;
+
+        // Face 4 == z-min
+        case 4:
+        {
+            for (label i=0; i<ni; ++i)
+            {
+                for (label j=0; j<nj; ++j)
+                {
+                    auto& f = *faceIter;
+                    ++faceIter;
+                    f.resize(4);
+
+                    f[0] = pointLabel(i,   j,   0);
+                    f[1] = pointLabel(i,   j+1, 0);
+                    f[2] = pointLabel(i+1, j+1, 0);
+                    f[3] = pointLabel(i+1, j,   0);
+
+                    *ownIter = cellLabel(i, j, 0);
+                    ++ownIter;
+                }
+            }
+        }
+        break;
+
+        // Face 5 == z-max
+        case 5:
+        {
+            for (label i=0; i<ni; ++i)
+            {
+                for (label j=0; j<nj; ++j)
+                {
+                    auto& f = *faceIter;
+                    ++faceIter;
+                    f.resize(4);
+
+                    f[0] = pointLabel(i,   j,   nk);
+                    f[1] = pointLabel(i+1, j,   nk);
+                    f[2] = pointLabel(i+1, j+1, nk);
+                    f[3] = pointLabel(i,   j+1, nk);
+
+                    *ownIter = cellLabel(i, j, nk-1);
+                    ++ownIter;
+                }
+            }
+        }
+        break;
+    }
+
+    // Return the number of faces added
+    return (ownIter - firstIter);
+}
+
+
+Foam::autoPtr<Foam::polyMesh> Foam::PDRblock::mesh(const IOobject& io) const
+{
+    pointField pts(nPoints());
+
+    faceList faces(nFaces());
+    labelList own(nFaces());
+    labelList nei(nInternalFaces());
+
+    auto faceIter = faces.begin();
+    auto ownIter  = own.begin();
+    auto neiIter  = nei.begin();
+
+    createPoints(pts);
+
+    addInternalFaces(faceIter, ownIter, neiIter);
+
+    // After readBoundary() we have a complete list of patches
+    // without any conflicts, and the correct size per-patch
+
+    // Add boundary faces and adjust patch sizes
+    for (const boundaryEntry& bentry : patches_)
+    {
+        for (const label shapeFacei : bentry.faces_)
+        {
+            addBoundaryFaces(shapeFacei, faceIter, ownIter);
+        }
+    }
+
+    auto meshPtr = autoPtr<polyMesh>::New
+    (
+        io,
+        std::move(pts),
+        std::move(faces),
+        std::move(own),
+        std::move(nei)
+    );
+
+    PtrList<polyPatch> patches(patches_.size());
+
+    label startFace = nInternalFaces();
+
+    label patchi = 0;
+
+    for (const boundaryEntry& bentry : patches_)
+    {
+        patches.set
+        (
+            patchi,
+            polyPatch::New
+            (
+                bentry.type_,
+                bentry.name_,
+                bentry.size_,
+                startFace,
+                patchi,  // index
+                meshPtr->boundaryMesh()
+            )
+        );
+
+        // physicalType?
+
+        startFace += bentry.size_;
+        ++patchi;
+    }
+
+    meshPtr->addPatches(patches);
+
+    return meshPtr;
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H b/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H
new file mode 100644
index 0000000000000000000000000000000000000000..64becbd4b8631fe704239d1d51ca8ac191d7102b
--- /dev/null
+++ b/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H
@@ -0,0 +1,213 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::PDRblock::PDRblock()
+:
+    ijkMesh(),
+    verbose_(false),
+    grid_(),
+    bounds_(),
+    patches_()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::label Foam::PDRblock::location::nCells() const
+{
+    return (scalarList::size()-1);
+}
+
+
+inline Foam::label Foam::PDRblock::location::nPoints() const
+{
+    return scalarList::size();
+}
+
+
+inline Foam::label Foam::PDRblock::location::size() const
+{
+    return (scalarList::size()-1);
+}
+
+
+inline void Foam::PDRblock::location::checkIndex(const label i) const
+{
+    if (i < 0 || i >= size())
+    {
+        FatalErrorInFunction
+            << "The index " << i
+            << " is out of range [0," << size() << ']' << nl
+            << abort(FatalError);
+    }
+}
+
+
+inline Foam::scalar Foam::PDRblock::location::width(const label i) const
+{
+    #ifdef FULLDEBUG
+    checkIndex(i);
+    #endif
+
+    return (operator[](i+1) - operator[](i));
+}
+
+
+inline Foam::scalar Foam::PDRblock::location::C(const label i) const
+{
+    #ifdef FULLDEBUG
+    checkIndex(i);
+    #endif
+
+    return (operator[](i+1) + operator[](i));
+}
+
+
+inline const Foam::Vector<Foam::PDRblock::location>&
+Foam::PDRblock::grid() const
+{
+    return grid_;
+}
+
+
+inline const Foam::boundBox& Foam::PDRblock::bounds() const
+{
+    return bounds_;
+}
+
+
+inline Foam::scalar Foam::PDRblock::dx(const label i) const
+{
+    return grid_.x().width(i);
+}
+
+
+inline Foam::scalar Foam::PDRblock::dx(const labelVector& ijk) const
+{
+    return grid_.x().width(ijk.x());
+}
+
+
+inline Foam::scalar Foam::PDRblock::dy(const label j) const
+{
+    return grid_.y().width(j);
+}
+
+
+inline Foam::scalar Foam::PDRblock::dy(const labelVector& ijk) const
+{
+    return grid_.y().width(ijk.y());
+}
+
+
+inline Foam::scalar Foam::PDRblock::dz(const label k) const
+{
+    return grid_.z().width(k);
+}
+
+
+inline Foam::scalar Foam::PDRblock::dz(const labelVector& ijk) const
+{
+    return grid_.z().width(ijk.z());
+}
+
+
+inline Foam::vector Foam::PDRblock::span
+(
+    const label i,
+    const label j,
+    const label k
+) const
+{
+    return vector(dx(i), dy(j), dz(k));
+}
+
+
+inline Foam::vector Foam::PDRblock::span(const labelVector& ijk) const
+{
+    return vector(dx(ijk), dy(ijk), dz(ijk));
+}
+
+
+inline Foam::point Foam::PDRblock::C
+(
+    const label i,
+    const label j,
+    const label k
+) const
+{
+    return point(grid_.x().C(i), grid_.y().C(j), grid_.z().C(k));
+}
+
+
+inline Foam::point Foam::PDRblock::C(const labelVector& ijk) const
+{
+    return
+        point
+        (
+            grid_.x().C(ijk.x()),
+            grid_.y().C(ijk.y()),
+            grid_.z().C(ijk.z())
+        );
+}
+
+
+inline Foam::scalar Foam::PDRblock::V
+(
+    const label i,
+    const label j,
+    const label k
+) const
+{
+    return dx(i)*dy(j)*dz(k);
+}
+
+
+inline Foam::scalar Foam::PDRblock::V(const labelVector& ijk) const
+{
+    return dx(ijk.x())*dy(ijk.y())*dz(ijk.z());
+}
+
+
+inline Foam::scalar Foam::PDRblock::width
+(
+    const label i,
+    const label j,
+    const label k
+) const
+{
+    return Foam::cbrt(V(i, j, k));
+}
+
+
+inline Foam::scalar Foam::PDRblock::width(const labelVector& ijk) const
+{
+    return Foam::cbrt(V(ijk));
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/PDRblockMesh/box0/Allrun b/tutorials/mesh/PDRblockMesh/box0/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..dc4cece02f99480c6245dd4860d047e3f14af1b7
--- /dev/null
+++ b/tutorials/mesh/PDRblockMesh/box0/Allrun
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1                        # Run from this directory
+. $WM_PROJECT_DIR/bin/tools/RunFunctions    # Tutorial run functions
+
+runApplication PDRblockMesh
+runApplication checkMesh
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/mesh/PDRblockMesh/box0/system/PDRblockMeshDict b/tutorials/mesh/PDRblockMesh/box0/system/PDRblockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..2a04fb12f6bc36aeba6640bd2e3645f9ef96a542
--- /dev/null
+++ b/tutorials/mesh/PDRblockMesh/box0/system/PDRblockMeshDict
@@ -0,0 +1,78 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      PDRblockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Point scaling (optional)
+scale 1.0;
+
+x
+{
+    points  ( -13.28 -0.10 6.0 19.19 );
+    nCells  ( 10 12 10 );
+    ratios  ( -5.16 1 5.16 );
+}
+
+y
+{
+    points  ( -12.98 0 5.50 18.48 );
+    nCells  ( 10 11 10 );
+    ratios  ( -5.16 1 5.16 );
+}
+
+z
+{
+    points  ( 0.00 4.80 17.26 );
+    nCells  ( 10 10 );
+    ratios  ( 1 5.16 );
+}
+
+
+defaultPatch
+{
+    name    walls;
+    type    wall;
+}
+
+
+// Faces: 0=x-min, 1=x-max, 2=y-min, 3=y-max, 4=z-min, 5=z-max
+
+boundary
+(
+    outer
+    {
+        type    patch;
+        faces   ( 0 1 2 3 5 );
+    }
+
+    mergingFaces
+    {
+        type    wall;
+        faces   ();
+    }
+
+    blockedFaces
+    {
+        type    wall;
+        faces   ();
+    }
+
+    ground
+    {
+        type    wall;
+        faces   ( 4 );
+    }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/PDRblockMesh/box0/system/controlDict b/tutorials/mesh/PDRblockMesh/box0/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..a5a945dc8d0886e0443395f54d2efac7ffda5c65
--- /dev/null
+++ b/tutorials/mesh/PDRblockMesh/box0/system/controlDict
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+libs            ("libblockMesh.so");
+
+application     PDRblockMesh;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0;
+
+deltaT          0;
+
+writeControl    timeStep;
+
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  8;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/PDRblockMesh/box0/system/fvSchemes b/tutorials/mesh/PDRblockMesh/box0/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..7405277d729114929bb44f27b4b9cc6dfae9de0d
--- /dev/null
+++ b/tutorials/mesh/PDRblockMesh/box0/system/fvSchemes
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{}
+
+gradSchemes
+{}
+
+divSchemes
+{}
+
+laplacianSchemes
+{}
+
+interpolationSchemes
+{}
+
+snGradSchemes
+{}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/PDRblockMesh/box0/system/fvSolution b/tutorials/mesh/PDRblockMesh/box0/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..41d5821731215142aa28599c2f8fddc0c04bd836
--- /dev/null
+++ b/tutorials/mesh/PDRblockMesh/box0/system/fvSolution
@@ -0,0 +1,18 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ************************************************************************* //