diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index dd8bfc5246..f99aaab857 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -62,7 +62,7 @@ Description
 #include "fvMeshTools.H"
 #include "profiling.H"
 #include "processorMeshes.H"
-#include "vtkSetWriter.H"
+#include "snappyVoxelMeshDriver.H"
 
 using namespace Foam;
 
@@ -325,7 +325,8 @@ autoPtr<refinementSurfaces> createRefinementSurfaces
             maxLevel,
             gapLevel,
             scalarField(nRegions, -GREAT),  //perpendicularAngle,
-            patchInfo
+            patchInfo,
+            false                           //dryRun
         )
     );
 
@@ -537,8 +538,61 @@ void extractSurface
 }
 
 
+label checkAlignment(const polyMesh& mesh, const scalar tol, Ostream& os)
+{
+    // Check all edges aligned with one of the coordinate axes
+    const faceList& faces = mesh.faces();
+    const pointField& points = mesh.points();
+
+    label nUnaligned = 0;
+
+    forAll(faces, facei)
+    {
+        const face& f = faces[facei];
+        forAll(f, fp)
+        {
+            label fp1 = f.fcIndex(fp);
+            const linePointRef e(edge(f[fp], f[fp1]).line(points));
+            const vector v(e.vec());
+            const scalar magV(mag(v));
+            if (magV > ROOTVSMALL)
+            {
+                for
+                (
+                    direction dir = 0;
+                    dir < pTraits<vector>::nComponents;
+                    ++dir
+                )
+                {
+                    const scalar s(mag(v[dir]));
+                    if (s > magV*tol && s < magV*(1-tol))
+                    {
+                        ++nUnaligned;
+                        break;
+                    }
+                }
+            }
+        }
+    }
+
+    reduce(nUnaligned, sumOp<label>());
+
+    if (nUnaligned)
+    {
+        os << "Initial mesh has " << nUnaligned
+            << " edges unaligned with any of the coordinate axes" << nl << endl;
+    }
+    return nUnaligned;
+}
+
+
 // Check writing tolerance before doing any serious work
-scalar getMergeDistance(const polyMesh& mesh, const scalar mergeTol)
+scalar getMergeDistance
+(
+    const polyMesh& mesh,
+    const scalar mergeTol,
+    const bool dryRun
+)
 {
     const boundBox& meshBb = mesh.bounds();
     scalar mergeDist = mergeTol * meshBb.mag();
@@ -550,7 +604,7 @@ scalar getMergeDistance(const polyMesh& mesh, const scalar mergeTol)
         << endl;
 
     // check writing tolerance
-    if (mesh.time().writeFormat() == IOstream::ASCII)
+    if (mesh.time().writeFormat() == IOstream::ASCII && !dryRun)
     {
         const scalar writeTol = std::pow
         (
@@ -690,6 +744,11 @@ int main(int argc, char *argv[])
         "checkGeometry",
         "Check all surface geometry for quality"
     );
+    argList::addBoolOption
+    (
+        "dry-run",
+        "Check case set-up only using a single time step"
+    );
     argList::addOption
     (
         "surfaceSimplify",
@@ -718,39 +777,17 @@ int main(int argc, char *argv[])
     const bool overwrite = args.found("overwrite");
     const bool checkGeometry = args.found("checkGeometry");
     const bool surfaceSimplify = args.found("surfaceSimplify");
+    const bool dryRun = args.optionFound("dry-run");
 
-    autoPtr<fvMesh> meshPtr;
-
+    if (dryRun)
     {
-        word regionName = fvMesh::defaultRegion;
-        if (args.readIfPresent("region", regionName))
-        {
-            Info<< "Create mesh " << regionName << " for time = "
-                << runTime.timeName() << nl << endl;
-        }
-        else
-        {
-            Info<< "Create mesh for time = "
-                << runTime.timeName() << nl << endl;
-        }
-
-        meshPtr.reset
-        (
-            new fvMesh
-            (
-                IOobject
-                (
-                    regionName,
-                    runTime.timeName(),
-                    runTime,
-                    IOobject::MUST_READ
-                )
-            )
-        );
+        Info<< "Operating in dry-run mode to detect set-up errors"
+            << nl << endl;
     }
 
-    fvMesh& mesh = meshPtr();
 
+
+    #include "createNamedMesh.H"
     Info<< "Read mesh in = "
         << runTime.cpuTimeIncrement() << " s" << endl;
 
@@ -758,6 +795,13 @@ int main(int argc, char *argv[])
     mesh.boundaryMesh().checkParallelSync(true);
     meshRefinement::checkCoupledFaceZones(mesh);
 
+    if (dryRun)
+    {
+        // Check if mesh aligned with cartesian axes
+        checkAlignment(mesh, 1e-6, Pout);   //FatalIOError);
+    }
+
+
 
     // Read meshing dictionary
     const word dictName("snappyHexMeshDict");
@@ -766,25 +810,36 @@ int main(int argc, char *argv[])
 
 
     // all surface geometry
-    const dictionary& geometryDict = meshDict.subDict("geometry");
+    const dictionary& geometryDict =
+        meshRefinement::subDict(meshDict, "geometry", dryRun);
 
     // refinement parameters
-    const dictionary& refineDict = meshDict.subDict("castellatedMeshControls");
+    const dictionary& refineDict =
+        meshRefinement::subDict(meshDict, "castellatedMeshControls", dryRun);
 
     // mesh motion and mesh quality parameters
-    const dictionary& motionDict = meshDict.subDict("meshQualityControls");
+    const dictionary& motionDict =
+        meshRefinement::subDict(meshDict, "meshQualityControls", dryRun);
 
     // snap-to-surface parameters
-    const dictionary& snapDict = meshDict.subDict("snapControls");
+    const dictionary& snapDict =
+        meshRefinement::subDict(meshDict, "snapControls", dryRun);
 
     // layer addition parameters
-    const dictionary& layerDict = meshDict.subDict("addLayersControls");
+    const dictionary& layerDict =
+        meshRefinement::subDict(meshDict, "addLayersControls", dryRun);
 
     // absolute merge distance
     const scalar mergeDist = getMergeDistance
     (
         mesh,
-        meshDict.get<scalar>("mergeTolerance")
+        meshRefinement::get<scalar>
+        (
+            meshDict,
+            "mergeTolerance",
+            dryRun
+        ),
+        dryRun
     );
 
     const bool keepPatches(meshDict.lookupOrDefault("keepPatches", false));
@@ -975,7 +1030,7 @@ int main(int argc, char *argv[])
         const scalar defaultCellSize =
             motionDict.get<scalar>("defaultCellSize");
 
-        const scalar initialCellSize = ::pow(meshPtr().V()[0], 1.0/3.0);
+        const scalar initialCellSize = ::pow(mesh.V()[0], 1.0/3.0);
 
         //Info<< "Wanted cell size  = " << defaultCellSize << endl;
         //Info<< "Current cell size = " << initialCellSize << endl;
@@ -1001,8 +1056,14 @@ int main(int argc, char *argv[])
             new refinementSurfaces
             (
                 allGeometry,
-                refineDict.subDict("refinementSurfaces"),
-                refineDict.lookupOrDefault("gapLevelIncrement", 0)
+                meshRefinement::subDict
+                (
+                    refineDict,
+                    "refinementSurfaces",
+                    dryRun
+                ),
+                refineDict.lookupOrDefault("gapLevelIncrement", 0),
+                dryRun
             )
         );
 
@@ -1017,6 +1078,8 @@ int main(int argc, char *argv[])
 
     if (checkGeometry)
     {
+        // Check geometry amongst itself (e.g. intersection, size differences)
+
         // Extract patchInfo
         List<wordList> patchTypes(allGeometry.size());
 
@@ -1035,7 +1098,14 @@ int main(int argc, char *argv[])
                 if (patchInfo.set(globalRegioni))
                 {
                     patchTypes[geomi][regioni] =
-                        patchInfo[globalRegioni].get<word>("type");
+                        meshRefinement::get<word>
+                        (
+                            patchInfo[globalRegioni],
+                            "type",
+                            dryRun,
+                            keyType::REGEX,
+                            word::null
+                        );
                 }
                 else
                 {
@@ -1058,10 +1128,49 @@ int main(int argc, char *argv[])
             true
         );
 
-        return 0;
+        if (!dryRun)
+        {
+            return 0;
+        }
     }
 
 
+    if (dryRun)
+    {
+        // Check geometry to mesh bounding box
+        Info<< "Checking for geometry size relative to mesh." << endl;
+        const boundBox& meshBb = mesh.bounds();
+        forAll(allGeometry, geomi)
+        {
+            const searchableSurface& s = allGeometry[geomi];
+            const boundBox& bb = s.bounds();
+
+            scalar ratio = bb.mag() / meshBb.mag();
+            if (ratio > 100.0 || ratio < 1.0/100.0)
+            {
+                Warning
+                    << "    " << allGeometry.names()[geomi]
+                    << " bounds differ from mesh"
+                    << " by more than a factor 100:" << nl
+                    << "        bounding box      : " << bb << nl
+                    << "        mesh bounding box : " << meshBb
+                    << endl;
+            }
+            if (!meshBb.contains(bb))
+            {
+                Warning
+                    << "    " << allGeometry.names()[geomi]
+                    << " bounds not fully contained in mesh" << nl
+                    << "        bounding box      : " << bb << nl
+                    << "        mesh bounding box : " << meshBb
+                    << endl;
+            }
+        }
+        Info<< endl;
+    }
+
+
+
 
     // Read refinement shells
     // ~~~~~~~~~~~~~~~~~~~~~~
@@ -1070,7 +1179,8 @@ int main(int argc, char *argv[])
     shellSurfaces shells
     (
         allGeometry,
-        refineDict.subDict("refinementRegions")
+        meshRefinement::subDict(refineDict, "refinementRegions", dryRun),
+        dryRun
     );
     Info<< "Read refinement shells in = "
         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
@@ -1093,7 +1203,7 @@ int main(int argc, char *argv[])
         Info<< "Reading limit shells." << endl;
     }
 
-    shellSurfaces limitShells(allGeometry, limitDict);
+    shellSurfaces limitShells(allGeometry, limitDict, dryRun);
 
     if (!limitDict.empty())
     {
@@ -1101,6 +1211,29 @@ int main(int argc, char *argv[])
             << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
     }
 
+    if (dryRun)
+    {
+        // Check for use of all geometry
+        const wordList& allGeomNames = allGeometry.names();
+
+        labelHashSet unusedGeometries(identity(allGeomNames.size()));
+        unusedGeometries.erase(surfaces.surfaces());
+        unusedGeometries.erase(shells.shells());
+        unusedGeometries.erase(limitShells.shells());
+
+        if (unusedGeometries.size())
+        {
+            IOWarningInFunction(geometryDict)
+                << "The following geometry entries are not used:" << nl;
+            for (const label geomi : unusedGeometries)
+            {
+                Info<< "    " << allGeomNames[geomi] << nl;
+            }
+            Info<< endl;
+        }
+    }
+
+
 
 
     // Read feature meshes
@@ -1110,12 +1243,33 @@ int main(int argc, char *argv[])
     refinementFeatures features
     (
         mesh,
-        refineDict.lookup("features")
+        meshRefinement::lookup(refineDict, "features", dryRun),
+        dryRun
     );
     Info<< "Read features in = "
         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 
 
+    if (dryRun)
+    {
+        // Check geometry to mesh bounding box
+        Info<< "Checking for line geometry size relative to surface geometry."
+            << endl;
+
+        OStringStream os;
+        bool hasErrors = features.checkSizes
+        (
+            100.0,  //const scalar maxRatio,
+            mesh.bounds(),
+            true,   //const bool report,
+            os //FatalIOError
+        );
+        if (hasErrors)
+        {
+            Warning<< os.str() << endl;
+        }
+    }
+
 
     // Refinement engine
     // ~~~~~~~~~~~~~~~~~
@@ -1134,10 +1288,17 @@ int main(int argc, char *argv[])
         surfaces,           // for surface intersection refinement
         features,           // for feature edges/point based refinement
         shells,             // for volume (inside/outside) refinement
-        limitShells         // limit of volume refinement
+        limitShells,        // limit of volume refinement
+        labelList(),        // initial faces to test
+        dryRun
     );
-    Info<< "Calculated surface intersections in = "
-        << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
+
+    if (!dryRun)
+    {
+        meshRefiner.updateIntersections(identity(mesh.nFaces()));
+        Info<< "Calculated surface intersections in = "
+            << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
+    }
 
     // Some stats
     meshRefiner.printMeshInfo(debugLevel, "Initial mesh");
@@ -1151,10 +1312,10 @@ int main(int argc, char *argv[])
 
 
     // Refinement parameters
-    const refinementParameters refineParams(refineDict);
+    const refinementParameters refineParams(refineDict, dryRun);
 
     // Snap parameters
-    const snapParameters snapParams(snapDict);
+    const snapParameters snapParams(snapDict, dryRun);
 
 
 
@@ -1464,9 +1625,46 @@ int main(int argc, char *argv[])
     // Now do the real work -refinement -snapping -layers
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    const bool wantRefine(meshDict.get<bool>("castellatedMesh"));
-    const bool wantSnap(meshDict.get<bool>("snap"));
-    const bool wantLayers(meshDict.get<bool>("addLayers"));
+    const bool wantRefine
+    (
+        meshRefinement::get<bool>(meshDict, "castellatedMesh", dryRun)
+    );
+    const bool wantSnap
+    (
+        meshRefinement::get<bool>(meshDict, "snap", dryRun)
+    );
+    const bool wantLayers
+    (
+        meshRefinement::get<bool>(meshDict, "addLayers", dryRun)
+    );
+
+    if (dryRun)
+    {
+        string errorMsg(FatalError.message());
+        string IOerrorMsg(FatalIOError.message());
+
+        if (errorMsg.size() || IOerrorMsg.size())
+        {
+            //errorMsg = "[dryRun] " + errorMsg;
+            //errorMsg.replaceAll("\n", "\n[dryRun] ");
+            //IOerrorMsg = "[dryRun] " + IOerrorMsg;
+            //IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
+
+            Warning
+                << nl
+                << "Missing/incorrect required dictionary entries:" << nl
+                << nl
+                << IOerrorMsg.c_str() << nl
+                << errorMsg.c_str() << nl << nl
+                << "Exiting dry-run" << nl << endl;
+
+            FatalError.clear();
+            FatalIOError.clear();
+
+            return 0;
+        }
+    }
+
 
     const bool mergePatchFaces
     (
@@ -1492,7 +1690,8 @@ int main(int argc, char *argv[])
             distributor,
             globalToMasterPatch,
             globalToSlavePatch,
-            setFormatter
+            setFormatter,
+            dryRun
         );
 
 
@@ -1518,13 +1717,16 @@ int main(int argc, char *argv[])
             fvMeshTools::removeEmptyPatches(mesh, true);
         }
 
-        writeMesh
-        (
-            "Refined mesh",
-            meshRefiner,
-            debugLevel,
-            meshRefinement::writeLevel()
-        );
+        if (!dryRun)
+        {
+            writeMesh
+            (
+                "Refined mesh",
+                meshRefiner,
+                debugLevel,
+                meshRefinement::writeLevel()
+            );
+        }
 
         Info<< "Mesh refined in = "
             << timer.cpuTimeIncrement() << " s." << endl;
@@ -1540,7 +1742,8 @@ int main(int argc, char *argv[])
         (
             meshRefiner,
             globalToMasterPatch,
-            globalToSlavePatch
+            globalToSlavePatch,
+            dryRun
         );
 
         if (!overwrite && !debugLevel)
@@ -1568,13 +1771,16 @@ int main(int argc, char *argv[])
             fvMeshTools::removeEmptyPatches(mesh, true);
         }
 
-        writeMesh
-        (
-            "Snapped mesh",
-            meshRefiner,
-            debugLevel,
-            meshRefinement::writeLevel()
-        );
+        if (!dryRun)
+        {
+            writeMesh
+            (
+                "Snapped mesh",
+                meshRefiner,
+                debugLevel,
+                meshRefinement::writeLevel()
+            );
+        }
 
         Info<< "Mesh snapped in = "
             << timer.cpuTimeIncrement() << " s." << endl;
@@ -1587,13 +1793,19 @@ int main(int argc, char *argv[])
         cpuTime timer;
 
         // Layer addition parameters
-        const layerParameters layerParams(layerDict, mesh.boundaryMesh());
+        const layerParameters layerParams
+        (
+            layerDict,
+            mesh.boundaryMesh(),
+            dryRun
+        );
 
         snappyLayerDriver layerDriver
         (
             meshRefiner,
             globalToMasterPatch,
-            globalToSlavePatch
+            globalToSlavePatch,
+            dryRun
         );
 
         // Use the maxLocalCells from the refinement parameters
@@ -1626,13 +1838,16 @@ int main(int argc, char *argv[])
             fvMeshTools::removeEmptyPatches(mesh, true);
         }
 
-        writeMesh
-        (
-            "Layer mesh",
-            meshRefiner,
-            debugLevel,
-            meshRefinement::writeLevel()
-        );
+        if (!dryRun)
+        {
+            writeMesh
+            (
+                "Layer mesh",
+                meshRefiner,
+                debugLevel,
+                meshRefinement::writeLevel()
+            );
+        }
 
         Info<< "Layers added in = "
             << timer.cpuTimeIncrement() << " s." << endl;
@@ -1647,7 +1862,7 @@ int main(int argc, char *argv[])
         // Check final mesh
         Info<< "Checking final mesh ..." << endl;
         faceSet wrongFaces(mesh, "wrongFaces", mesh.nFaces()/100);
-        motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
+        motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces, dryRun);
         const label nErrors = returnReduce
         (
             wrongFaces.size(),
@@ -1736,6 +1951,34 @@ int main(int argc, char *argv[])
     Info<< "Finished meshing in = "
         << runTime.elapsedCpuTime() << " s." << endl;
 
+
+    if (dryRun)
+    {
+        string errorMsg(FatalError.message());
+        string IOerrorMsg(FatalIOError.message());
+
+        if (errorMsg.size() || IOerrorMsg.size())
+        {
+            //errorMsg = "[dryRun] " + errorMsg;
+            //errorMsg.replaceAll("\n", "\n[dryRun] ");
+            //IOerrorMsg = "[dryRun] " + IOerrorMsg;
+            //IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
+
+            Perr<< nl
+                << "Missing/incorrect required dictionary entries:" << nl
+                << nl
+                << IOerrorMsg.c_str() << nl
+                << errorMsg.c_str() << nl << nl
+                << "Exiting dry-run" << nl << endl;
+
+            FatalError.clear();
+            FatalIOError.clear();
+
+            return 0;
+        }
+    }
+
+
     Info<< "End\n" << endl;
 
     return 0;
diff --git a/src/OpenFOAM/db/error/error.C b/src/OpenFOAM/db/error/error.C
index 4cc8b0a8fb..1ff7a96919 100644
--- a/src/OpenFOAM/db/error/error.C
+++ b/src/OpenFOAM/db/error/error.C
@@ -203,6 +203,12 @@ Foam::string Foam::error::message() const
 }
 
 
+void Foam::error::clear() const
+{
+    return messageStreamPtr_->reset();
+}
+
+
 void Foam::error::exit(const int errNo)
 {
     if (!throwExceptions_ && JobInfo::constructed)
diff --git a/src/OpenFOAM/db/error/error.H b/src/OpenFOAM/db/error/error.H
index 79dd7f3166..39d45d72e8 100644
--- a/src/OpenFOAM/db/error/error.H
+++ b/src/OpenFOAM/db/error/error.H
@@ -114,6 +114,9 @@ public:
 
         string message() const;
 
+        //- Clear any messages
+        void clear() const;
+
         inline const string& functionName() const
         {
             return functionName_;
diff --git a/src/dynamicMesh/motionSmoother/motionSmoother.C b/src/dynamicMesh/motionSmoother/motionSmoother.C
index d30ee35732..9c4a300e16 100644
--- a/src/dynamicMesh/motionSmoother/motionSmoother.C
+++ b/src/dynamicMesh/motionSmoother/motionSmoother.C
@@ -41,7 +41,8 @@ Foam::motionSmoother::motionSmoother
     pointMesh& pMesh,
     indirectPrimitivePatch& pp,
     const labelList& adaptPatchIDs,
-    const dictionary& paramDict
+    const dictionary& paramDict,
+    const bool dryRun
 )
 :
     motionSmootherData(pMesh),
@@ -54,7 +55,8 @@ Foam::motionSmoother::motionSmoother
         motionSmootherData::scale_,
         motionSmootherData::oldPoints_,
         adaptPatchIDs,
-        paramDict
+        paramDict,
+        dryRun
     )
 {}
 
@@ -65,7 +67,8 @@ Foam::motionSmoother::motionSmoother
     indirectPrimitivePatch& pp,
     const labelList& adaptPatchIDs,
     const pointVectorField& displacement,
-    const dictionary& paramDict
+    const dictionary& paramDict,
+    const bool dryRun
 )
 :
     motionSmootherData(displacement),
@@ -78,7 +81,8 @@ Foam::motionSmoother::motionSmoother
         motionSmootherData::scale_,
         motionSmootherData::oldPoints_,
         adaptPatchIDs,
-        paramDict
+        paramDict,
+        dryRun
     )
 {}
 
diff --git a/src/dynamicMesh/motionSmoother/motionSmoother.H b/src/dynamicMesh/motionSmoother/motionSmoother.H
index 2b45646004..31c29bb8b5 100644
--- a/src/dynamicMesh/motionSmoother/motionSmoother.H
+++ b/src/dynamicMesh/motionSmoother/motionSmoother.H
@@ -107,7 +107,8 @@ public:
             pointMesh& pMesh,
             indirectPrimitivePatch& pp,         //!< 'outside' points
             const labelList& adaptPatchIDs,     //!< patches forming 'outside'
-            const dictionary& paramDict
+            const dictionary& paramDict,
+            const bool dryRun = false
         );
 
         //- Construct from mesh, patches to work on and smoothing parameters and
@@ -118,7 +119,8 @@ public:
             indirectPrimitivePatch& pp,         //!< 'outside' points
             const labelList& adaptPatchIDs,     //!< patches forming 'outside'
             const pointVectorField& displacement,
-            const dictionary& paramDict
+            const dictionary& paramDict,
+            const bool dryRun = false
         );
 };
 
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
index ee8a1137e4..35b5c1cdf7 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
@@ -328,7 +328,8 @@ Foam::motionSmootherAlgo::motionSmootherAlgo
     pointScalarField& scale,
     pointField& oldPoints,
     const labelList& adaptPatchIDs,
-    const dictionary& paramDict
+    const dictionary& paramDict,
+    const bool dryRun
 )
 :
     mesh_(mesh),
@@ -339,6 +340,7 @@ Foam::motionSmootherAlgo::motionSmootherAlgo
     oldPoints_(oldPoints),
     adaptPatchIDs_(adaptPatchIDs),
     paramDict_(paramDict),
+    dryRun_(dryRun),
     isInternalPoint_(mesh_.nPoints(), true)
 {
     updateMesh();
@@ -856,9 +858,14 @@ bool Foam::motionSmootherAlgo::scaleMesh
         }
     }
 
-    const scalar errorReduction = paramDict.get<scalar>("errorReduction");
-    const label nSmoothScale = paramDict.get<label>("nSmoothScale");
-
+    const scalar errorReduction = get<scalar>
+    (
+        paramDict, "errorReduction", dryRun_, keyType::REGEX_RECURSIVE
+    );
+    const label nSmoothScale = get<label>
+    (
+        paramDict, "nSmoothScale", dryRun_, keyType::REGEX_RECURSIVE
+    );
 
     // Note: displacement_ should already be synced already from setDisplacement
     // but just to make sure.
@@ -885,7 +892,16 @@ bool Foam::motionSmootherAlgo::scaleMesh
 
     // Check. Returns parallel number of incorrect faces.
     faceSet wrongFaces(mesh_, "wrongFaces", mesh_.nFaces()/100+100);
-    checkMesh(false, mesh_, meshQualityDict, checkFaces, baffles, wrongFaces);
+    checkMesh
+    (
+        false,
+        mesh_,
+        meshQualityDict,
+        checkFaces,
+        baffles,
+        wrongFaces,
+        dryRun_
+    );
 
     if (returnReduce(wrongFaces.size(), sumOp<label>()) <= nAllowableErrors)
     {
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
index 742962e5bb..cabc041af2 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
@@ -156,9 +156,12 @@ class motionSmootherAlgo
             // displacement on.
             const labelList adaptPatchIDs_;
 
-            // Smoothing and checking parameters
+            //- Smoothing and checking parameters
             dictionary paramDict_;
 
+            //- In test/dry-run mode?
+            const bool dryRun_;
+
             //- Is mesh point on boundary or not
             bitSet isInternalPoint_;
 
@@ -303,7 +306,8 @@ public:
             pointScalarField& scale,
             pointField& oldPoints,
             const labelList& adaptPatchIDs,     // patches forming 'outside'
-            const dictionary& paramDict
+            const dictionary& paramDict,
+            const bool dryRun = false
         );
 
 
@@ -442,7 +446,8 @@ public:
                 const bool report,
                 const polyMesh& mesh,
                 const dictionary& dict,
-                labelHashSet& wrongFaces
+                labelHashSet& wrongFaces,
+                const bool dryRun = false
             );
 
             //- Check (subset of mesh) with mesh settings in dict.
@@ -454,7 +459,8 @@ public:
                 const polyMesh& mesh,
                 const dictionary& dict,
                 const labelList& checkFaces,
-                labelHashSet& wrongFaces
+                labelHashSet& wrongFaces,
+                const bool dryRun = false
             );
 
             //- Check (subset of mesh including baffles) with mesh settings
@@ -467,7 +473,8 @@ public:
                 const dictionary& dict,
                 const labelList& checkFaces,
                 const List<labelPair>& baffles,
-                labelHashSet& wrongFaces
+                labelHashSet& wrongFaces,
+                const bool dryRun = false
             );
 
             //- Check part of mesh with mesh settings in dict.
@@ -480,7 +487,8 @@ public:
                 const polyMeshGeometry&,
                 const pointField&,
                 const labelList& checkFaces,
-                labelHashSet& wrongFaces
+                labelHashSet& wrongFaces,
+                const bool dryRun = false
             );
 
             //- Check part of mesh including baffles with mesh settings in dict.
@@ -494,7 +502,8 @@ public:
                 const pointField&,
                 const labelList& checkFaces,
                 const List<labelPair>& baffles,
-                labelHashSet& wrongFaces
+                labelHashSet& wrongFaces,
+                const bool dryRun = false
             );
 
             // Helper functions to manipulate displacement vector.
@@ -508,6 +517,28 @@ public:
                     const scalarField& edgeWeight,
                     GeometricField<Type, pointPatchField, pointMesh>& newFld
                 ) const;
+
+            //- Wrapper around dictionary::get which does not exit
+            template<class Type>
+            static Type get
+            (
+                const dictionary& dict,
+                const word& keyword,
+                const bool noExit,
+                enum keyType::option matchOpt,
+                const Type& defaultValue = Zero
+            );
+//
+//            //- Wrapper around dictionary::get which does not exit
+//            template<class Type>
+//            static Type get
+//            (
+//                const dictionary& dict,
+//                const word& keyword,
+//                const bool noExit,
+//                bool recursive,
+//                bool patternMatch
+//            );
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
index 361fbe216e..003067656f 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
@@ -35,7 +35,8 @@ bool Foam::motionSmootherAlgo::checkMesh
     const polyMesh& mesh,
     const dictionary& dict,
     const labelList& checkFaces,
-    labelHashSet& wrongFaces
+    labelHashSet& wrongFaces,
+    const bool dryRun
 )
 {
     List<labelPair> emptyBaffles;
@@ -46,7 +47,8 @@ bool Foam::motionSmootherAlgo::checkMesh
         dict,
         checkFaces,
         emptyBaffles,
-        wrongFaces
+        wrongFaces,
+        dryRun
     );
 }
 
@@ -57,54 +59,60 @@ bool Foam::motionSmootherAlgo::checkMesh
     const dictionary& dict,
     const labelList& checkFaces,
     const List<labelPair>& baffles,
-    labelHashSet& wrongFaces
+    labelHashSet& wrongFaces,
+    const bool dryRun
 )
 {
     const scalar maxNonOrtho
     (
-        dict.get<scalar>("maxNonOrtho", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "maxNonOrtho", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minVol
     (
-        dict.get<scalar>("minVol", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minVol", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minTetQuality
     (
-        dict.get<scalar>("minTetQuality", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minTetQuality", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar maxConcave
     (
-        dict.get<scalar>("maxConcave", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "maxConcave", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minArea
     (
-        dict.get<scalar>("minArea", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minArea", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar maxIntSkew
     (
-        dict.get<scalar>("maxInternalSkewness", keyType::REGEX_RECURSIVE)
+        get<scalar>
+        (
+            dict, "maxInternalSkewness", dryRun, keyType::REGEX_RECURSIVE
+        )
     );
     const scalar maxBounSkew
     (
-        dict.get<scalar>("maxBoundarySkewness", keyType::REGEX_RECURSIVE)
+        get<scalar>
+        (
+            dict, "maxBoundarySkewness", dryRun, keyType::REGEX_RECURSIVE
+        )
     );
     const scalar minWeight
     (
-        dict.get<scalar>("minFaceWeight", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minFaceWeight", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minVolRatio
     (
-        dict.get<scalar>("minVolRatio", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minVolRatio", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minTwist
     (
-        dict.get<scalar>("minTwist", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minTwist", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minTriangleTwist
     (
-        dict.get<scalar>("minTriangleTwist", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minTriangleTwist", dryRun, keyType::REGEX_RECURSIVE)
     );
-
     const scalar minFaceFlatness
     (
         dict.lookupOrDefault<scalar>
@@ -114,8 +122,37 @@ bool Foam::motionSmootherAlgo::checkMesh
     );
     const scalar minDet
     (
-        dict.get<scalar>("minDeterminant", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minDeterminant", dryRun, keyType::REGEX_RECURSIVE)
     );
+
+
+    if (dryRun)
+    {
+        string errorMsg(FatalError.message());
+        string IOerrorMsg(FatalIOError.message());
+
+        if (errorMsg.size() || IOerrorMsg.size())
+        {
+            //errorMsg = "[dryRun] " + errorMsg;
+            //errorMsg.replaceAll("\n", "\n[dryRun] ");
+            //IOerrorMsg = "[dryRun] " + IOerrorMsg;
+            //IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
+
+            IOWarningInFunction(dict)
+                << nl
+                << "Missing/incorrect required dictionary entries:" << nl
+                << nl
+                << IOerrorMsg.c_str() << nl
+                << errorMsg.c_str() << nl << nl
+                << "Exiting dry-run" << nl << endl;
+
+            FatalError.clear();
+            FatalIOError.clear();
+        }
+        return false;
+    }
+
+
     label nWrongFaces = 0;
 
     Info<< "Checking faces in error :" << endl;
@@ -422,7 +459,8 @@ bool Foam::motionSmootherAlgo::checkMesh
     const bool report,
     const polyMesh& mesh,
     const dictionary& dict,
-    labelHashSet& wrongFaces
+    labelHashSet& wrongFaces,
+    const bool dryRun
 )
 {
     return checkMesh
@@ -431,7 +469,8 @@ bool Foam::motionSmootherAlgo::checkMesh
         mesh,
         dict,
         identity(mesh.nFaces()),
-        wrongFaces
+        wrongFaces,
+        dryRun
     );
 }
 
@@ -442,7 +481,8 @@ bool Foam::motionSmootherAlgo::checkMesh
     const polyMeshGeometry& meshGeom,
     const pointField& points,
     const labelList& checkFaces,
-    labelHashSet& wrongFaces
+    labelHashSet& wrongFaces,
+    const bool dryRun
 )
 {
     List<labelPair> emptyBaffles;
@@ -455,7 +495,8 @@ bool Foam::motionSmootherAlgo::checkMesh
         points,
         checkFaces,
         emptyBaffles,
-        wrongFaces
+        wrongFaces,
+        dryRun
      );
 }
 
@@ -468,64 +509,92 @@ bool Foam::motionSmootherAlgo::checkMesh
     const pointField& points,
     const labelList& checkFaces,
     const List<labelPair>& baffles,
-    labelHashSet& wrongFaces
+    labelHashSet& wrongFaces,
+    const bool dryRun
 )
 {
     const scalar maxNonOrtho
     (
-        dict.get<scalar>("maxNonOrtho", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "maxNonOrtho", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minVol
     (
-        dict.get<scalar>("minVol", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minVol", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minTetQuality
     (
-        dict.get<scalar>("minTetQuality", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minTetQuality", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar maxConcave
     (
-        dict.get<scalar>("maxConcave", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "maxConcave", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minArea
     (
-        dict.get<scalar>("minArea", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minArea", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar maxIntSkew
     (
-        dict.get<scalar>("maxInternalSkewness", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "maxInternalSkewness", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar maxBounSkew
     (
-        dict.get<scalar>("maxBoundarySkewness", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "maxBoundarySkewness", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minWeight
     (
-        dict.get<scalar>("minFaceWeight", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minFaceWeight", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minVolRatio
     (
-        dict.get<scalar>("minVolRatio", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minVolRatio", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minTwist
     (
-        dict.get<scalar>("minTwist", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minTwist", dryRun, keyType::REGEX_RECURSIVE)
     );
     const scalar minTriangleTwist
     (
-        dict.get<scalar>("minTriangleTwist", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minTriangleTwist", dryRun, keyType::REGEX_RECURSIVE)
     );
-    const scalar minFaceFlatness
+    scalar minFaceFlatness = -1.0;
+    dict.readIfPresent
     (
-        dict.lookupOrDefault<scalar>
-        (
-            "minFaceFlatness", -1, keyType::REGEX_RECURSIVE
-        )
+        "minFaceFlatness",
+        minFaceFlatness,
+        keyType::REGEX_RECURSIVE
     );
     const scalar minDet
     (
-        dict.get<scalar>("minDeterminant", keyType::REGEX_RECURSIVE)
+        get<scalar>(dict, "minDeterminant", dryRun, keyType::REGEX_RECURSIVE)
     );
+
+    if (dryRun)
+    {
+        string errorMsg(FatalError.message());
+        string IOerrorMsg(FatalIOError.message());
+
+        if (errorMsg.size() || IOerrorMsg.size())
+        {
+            errorMsg = "[dryRun] " + errorMsg;
+            errorMsg.replaceAll("\n", "\n[dryRun] ");
+            IOerrorMsg = "[dryRun] " + IOerrorMsg;
+            IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
+
+            Perr<< nl
+                << "Missing/incorrect required dictionary entries:" << nl
+                << nl
+                << IOerrorMsg.c_str() << nl
+                << errorMsg.c_str() << nl << nl
+                << "Exiting dry-run" << nl << endl;
+
+            FatalError.clear();
+            FatalIOError.clear();
+        }
+        return false;
+    }
+
+
     label nWrongFaces = 0;
 
     Info<< "Checking faces in error :" << endl;
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
index a228c61bb7..eea46f3ee0 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
@@ -287,4 +287,35 @@ void Foam::motionSmootherAlgo::testSyncField
 }
 
 
+template<class Type>
+Type Foam::motionSmootherAlgo::get
+(
+    const dictionary& dict,
+    const word& keyword,
+    const bool noExit,
+    enum keyType::option matchOpt,
+    const Type& defaultValue
+)
+{
+    Type val(defaultValue);
+    
+    if
+    (
+       !dict.readEntry
+        (
+            keyword,
+            val,
+            matchOpt,
+            !noExit
+        )
+    )
+    {
+        FatalIOError
+            << "Entry '" << keyword << "' not found in dictionary "
+            << dict.name() << endl;
+    }
+    return val;
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/snappyHexMesh/Make/files b/src/mesh/snappyHexMesh/Make/files
index 49a981e4af..00c46b066e 100644
--- a/src/mesh/snappyHexMesh/Make/files
+++ b/src/mesh/snappyHexMesh/Make/files
@@ -2,6 +2,7 @@ snappyHexMeshDriver/snappyLayerDriver.C
 snappyHexMeshDriver/snappySnapDriver.C
 snappyHexMeshDriver/snappySnapDriverFeature.C
 snappyHexMeshDriver/snappyRefineDriver.C
+snappyHexMeshDriver/snappyVoxelMeshDriver.C
 
 snappyHexMeshDriver/layerParameters/layerParameters.C
 snappyHexMeshDriver/refinementParameters/refinementParameters.C
diff --git a/src/mesh/snappyHexMesh/Make/options b/src/mesh/snappyHexMesh/Make/options
index bff1be5f77..2cf98acf6b 100644
--- a/src/mesh/snappyHexMesh/Make/options
+++ b/src/mesh/snappyHexMesh/Make/options
@@ -1,10 +1,12 @@
 EXE_INC = \
+    /* -g -DFULLDEBUG -O0 */ \
     -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/overset/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
@@ -14,4 +16,5 @@ LIB_LIBS = \
     -llagrangian \
     -lmeshTools \
     -lfvMotionSolvers \
+    -loverset \
     -ldistributed
diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C
index 5fd4aa2af6..e4fa24a874 100644
--- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C
+++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C
@@ -109,10 +109,11 @@ Foam::displacementMotionSolverMeshMover::displacementMotionSolverMeshMover
 (
     const dictionary& dict,
     const List<labelPair>& baffles,
-    pointVectorField& pointDisplacement
+    pointVectorField& pointDisplacement,
+    const bool dryRun
 )
 :
-    externalDisplacementMeshMover(dict, baffles, pointDisplacement),
+    externalDisplacementMeshMover(dict, baffles, pointDisplacement, dryRun),
 
     solverPtr_
     (
@@ -178,7 +179,8 @@ Foam::displacementMotionSolverMeshMover::displacementMotionSolverMeshMover
         scale_,
         oldPoints_,
         adaptPatchIDs_,
-        dict
+        dict,
+        dryRun
     ),
 
     fieldSmoother_(mesh())
@@ -205,12 +207,15 @@ bool Foam::displacementMotionSolverMeshMover::move
     // Note that this has to update the pointDisplacement boundary conditions
     // as well, not just the internal field.
     {
-        const label nSmoothPatchThickness
+        const label nSmoothPatchThickness = meshRefinement::get<label>
         (
-            moveDict.get<label>("nSmoothThickness")
+            moveDict, "nSmoothThickness", dryRun_, keyType::REGEX
         );
 
-        const word minThicknessName(moveDict.get<word>("minThicknessName"));
+        const word minThicknessName = meshRefinement::get<word>
+        (
+            moveDict, "minThicknessName", dryRun_, keyType::REGEX, word::null
+        );
 
         scalarField zeroMinThickness;
 
diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H
index b4df2ef5c4..6462958f99 100644
--- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H
+++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H
@@ -112,7 +112,8 @@ public:
         (
             const dictionary& dict,
             const List<labelPair>& baffles,
-            pointVectorField& pointDisplacemen
+            pointVectorField& pointDisplacement,
+            const bool dryRun
         );
 
 
diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C
index 43ee2771ad..f7fa6213ee 100644
--- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C
+++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C
@@ -117,11 +117,13 @@ Foam::externalDisplacementMeshMover::externalDisplacementMeshMover
 (
     const dictionary& dict,
     const List<labelPair>& baffles,
-    pointVectorField& pointDisplacement
+    pointVectorField& pointDisplacement,
+    const bool dryRun
 )
 :
     baffles_(baffles),
-    pointDisplacement_(pointDisplacement)
+    pointDisplacement_(pointDisplacement),
+    dryRun_(dryRun)
 {}
 
 
@@ -133,7 +135,8 @@ Foam::externalDisplacementMeshMover::New
     const word& type,
     const dictionary& dict,
     const List<labelPair>& baffles,
-    pointVectorField& pointDisplacement
+    pointVectorField& pointDisplacement,
+    const bool dryRun
 )
 {
     Info<< "Selecting externalDisplacementMeshMover " << type << endl;
@@ -152,7 +155,7 @@ Foam::externalDisplacementMeshMover::New
 
     return autoPtr<externalDisplacementMeshMover>
     (
-        cstrIter()(dict, baffles, pointDisplacement)
+        cstrIter()(dict, baffles, pointDisplacement, dryRun)
     );
 }
 
diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H
index f9d65ef825..9656145f64 100644
--- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H
+++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H
@@ -66,6 +66,9 @@ protected:
         //- Reference to point motion field
         pointVectorField& pointDisplacement_;
 
+        //- In dry-run mode?
+        const bool dryRun_;
+
 
     // Protected Member functions
 
@@ -110,9 +113,10 @@ public:
             (
                 const dictionary& dict,
                 const List<labelPair>& baffles,
-                pointVectorField& pointDisplacement
+                pointVectorField& pointDisplacement,
+                const bool dryRun
             ),
-            (dict, baffles, pointDisplacement)
+            (dict, baffles, pointDisplacement, dryRun)
         );
 
 
@@ -124,7 +128,8 @@ public:
         (
             const dictionary& dict,
             const List<labelPair>& baffles,
-            pointVectorField& pointDisplacement
+            pointVectorField& pointDisplacement,
+            const bool dryRun
         );
 
 
@@ -136,7 +141,8 @@ public:
             const word& type,
             const dictionary& dict,
             const List<labelPair>& baffles,
-            pointVectorField& pointDisplacement
+            pointVectorField& pointDisplacement,
+            const bool dryRun = false
         );
 
 
diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
index 523b2a8bc5..6b044a8700 100644
--- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
+++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
@@ -138,38 +138,55 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
     //- Smooth surface normals
     const label nSmoothSurfaceNormals
     (
-        coeffDict.get<label>("nSmoothSurfaceNormals")
+        meshRefinement::get<label>
+        (
+            coeffDict,
+            "nSmoothSurfaceNormals",
+            dryRun_
+        )
     );
 
     // Note: parameter name changed
     // "minMedianAxisAngle" -> "minMedialAxisAngle" (DEC-2013)
     // but not previously reported.
-    scalar minMedialAxisAngleCos = Foam::cos
+    scalar minMedialAxisAngle(Zero);
+    if
     (
-        degToRad
+       !coeffDict.readCompat
         (
-            coeffDict.getCompat<scalar>
-            (
-                "minMedialAxisAngle", {{ "minMedianAxisAngle", 1712 }}
-            )
+            "minMedialAxisAngle",
+            {{ "minMedianAxisAngle", 1712 }},
+            minMedialAxisAngle,
+            keyType::REGEX,
+            !dryRun_
         )
-    );
+    )
+    {
+        FatalIOError
+            << "Entry '" << "minMedialAxisAngle"
+            << "' not found in dictionary " << coeffDict.name() << endl;
+    }
+
+    const scalar minMedialAxisAngleCos(Foam::cos(degToRad(minMedialAxisAngle)));
 
     //- Feature angle when to stop adding layers
-    const scalar featureAngle = coeffDict.get<scalar>("featureAngle");
+    const scalar featureAngle
+    (
+        meshRefinement::get<scalar>(coeffDict, "featureAngle", dryRun_)
+    );
 
     //- When to slip along wall
     const scalar slipFeatureAngle =
-        coeffDict.lookupOrDefault<scalar>
-        (
-            "slipFeatureAngle",
-            0.5*featureAngle
-        );
+    (
+        coeffDict.found("slipFeatureAngle")
+      ? coeffDict.get<scalar>("slipFeatureAngle")
+      : 0.5*featureAngle
+    );
 
     //- Smooth internal normals
     const label nSmoothNormals
     (
-        coeffDict.get<label>("nSmoothNormals")
+        meshRefinement::get<label>(coeffDict, "nSmoothNormals", dryRun_)
     );
 
     //- Number of edges walking out
@@ -1222,10 +1239,11 @@ Foam::medialAxisMeshMover::medialAxisMeshMover
 (
     const dictionary& dict,
     const List<labelPair>& baffles,
-    pointVectorField& pointDisplacement
+    pointVectorField& pointDisplacement,
+    const bool dryRun
 )
 :
-    externalDisplacementMeshMover(dict, baffles, pointDisplacement),
+    externalDisplacementMeshMover(dict, baffles, pointDisplacement, dryRun),
     adaptPatchIDs_(getFixedValueBCs(pointDisplacement)),
     adaptPatchPtr_(getPatch(mesh(), adaptPatchIDs_)),
     scale_
@@ -1251,7 +1269,8 @@ Foam::medialAxisMeshMover::medialAxisMeshMover
         scale_,
         oldPoints_,
         adaptPatchIDs_,
-        dict
+        dict,
+        dryRun
     ),
     fieldSmoother_(mesh()),
     dispVec_
@@ -1341,23 +1360,15 @@ void Foam::medialAxisMeshMover::calculateDisplacement
     // ~~~~~~~~~~~~~
 
     //- (lambda-mu) smoothing of internal displacement
-    const label nSmoothDisplacement =  coeffDict.lookupOrDefault
-    (
-        "nSmoothDisplacement",
-        0
-    );
+    const label nSmoothDisplacement =
+        coeffDict.lookupOrDefault("nSmoothDisplacement", 0);
 
     //- Layer thickness too big
-    const scalar maxThicknessToMedialRatio
-    (
-        coeffDict.get<scalar>("maxThicknessToMedialRatio")
-    );
+    const scalar maxThicknessToMedialRatio =
+        coeffDict.get<scalar>("maxThicknessToMedialRatio");
 
     //- Feature angle when to stop adding layers
-    const scalar featureAngle
-    (
-        coeffDict.get<scalar>("featureAngle")
-    );
+    const scalar featureAngle = coeffDict.get<scalar>("featureAngle");
 
     //- Stop layer growth where mesh wraps around sharp edge
     scalar layerTerminationAngle = coeffDict.lookupOrDefault<scalar>
@@ -1368,10 +1379,8 @@ void Foam::medialAxisMeshMover::calculateDisplacement
     scalar minCosLayerTermination = Foam::cos(degToRad(layerTerminationAngle));
 
     //- Smoothing wanted patch thickness
-    const label nSmoothPatchThickness
-    (
-        coeffDict.get<label>("nSmoothThickness")
-    );
+    const label nSmoothPatchThickness =
+        coeffDict.get<label>("nSmoothThickness");
 
     //- Number of edges walking out
     const label nMedialAxisIter = coeffDict.lookupOrDefault<label>
diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
index 8f3ea990f4..f47a99c19a 100644
--- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
+++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
@@ -207,7 +207,8 @@ public:
         (
             const dictionary& dict,
             const List<labelPair>& baffles,
-            pointVectorField& pointDisplacement
+            pointVectorField& pointDisplacement,
+            const bool dryRun
         );
 
 
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
index 404042246b..effc9cb935 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
@@ -487,7 +487,11 @@ void Foam::meshRefinement::checkData()
             neiBoundaryFc
         );
     }
+
     // Check meshRefinement
+    const labelList& surfIndex = surfaceIndex();
+
+
     {
         // Get boundary face centre and level. Coupled aware.
         labelList neiLevel(nBnd);
@@ -549,14 +553,14 @@ void Foam::meshRefinement::checkData()
         // Check
         forAll(surfaceHit, facei)
         {
-            if (surfaceIndex_[facei] != surfaceHit[facei])
+            if (surfIndex[facei] != surfaceHit[facei])
             {
                 if (mesh_.isInternalFace(facei))
                 {
                     WarningInFunction
                         << "Internal face:" << facei
                         << " fc:" << mesh_.faceCentres()[facei]
-                        << " cached surfaceIndex_:" << surfaceIndex_[facei]
+                        << " cached surfaceIndex_:" << surfIndex[facei]
                         << " current:" << surfaceHit[facei]
                         << " ownCc:"
                         << mesh_.cellCentres()[mesh_.faceOwner()[facei]]
@@ -566,14 +570,14 @@ void Foam::meshRefinement::checkData()
                 }
                 else if
                 (
-                    surfaceIndex_[facei]
+                    surfIndex[facei]
                  != neiHit[facei-mesh_.nInternalFaces()]
                 )
                 {
                     WarningInFunction
                         << "Boundary face:" << facei
                         << " fc:" << mesh_.faceCentres()[facei]
-                        << " cached surfaceIndex_:" << surfaceIndex_[facei]
+                        << " cached surfaceIndex_:" << surfIndex[facei]
                         << " current:" << surfaceHit[facei]
                         << " ownCc:"
                         << mesh_.cellCentres()[mesh_.faceOwner()[facei]]
@@ -996,7 +1000,8 @@ Foam::label Foam::meshRefinement::splitFacesUndo
             false,  // report
             mesh_,
             motionDict,
-            errorFaces
+            errorFaces,
+            dryRun_
         );
         if (!hasErrors)
         {
@@ -1208,7 +1213,9 @@ Foam::meshRefinement::meshRefinement
     const refinementSurfaces& surfaces,
     const refinementFeatures& features,
     const shellSurfaces& shells,
-    const shellSurfaces& limitShells
+    const shellSurfaces& limitShells,
+    const labelUList& checkFaces,
+    const bool dryRun
 )
 :
     mesh_(mesh),
@@ -1219,6 +1226,7 @@ Foam::meshRefinement::meshRefinement
     features_(features),
     shells_(shells),
     limitShells_(limitShells),
+    dryRun_(dryRun),
     meshCutter_
     (
         mesh,
@@ -1241,12 +1249,35 @@ Foam::meshRefinement::meshRefinement
     userFaceData_(0)
 {
     // recalculate intersections for all faces
-    updateIntersections(identity(mesh_.nFaces()));
+    updateIntersections(checkFaces);
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+const Foam::labelList& Foam::meshRefinement::surfaceIndex() const
+{
+    if (surfaceIndex_.size() != mesh_.nFaces())
+    {
+        const_cast<meshRefinement&>(*this).updateIntersections
+        (
+            identity(mesh_.nFaces())
+        );
+    }
+    return surfaceIndex_;
+}
+
+
+Foam::labelList& Foam::meshRefinement::surfaceIndex()
+{
+    if (surfaceIndex_.size() != mesh_.nFaces())
+    {
+        updateIntersections(identity(mesh_.nFaces()));
+    }
+    return surfaceIndex_;
+}
+
+
 Foam::label Foam::meshRefinement::countHits() const
 {
     // Stats on edges to test. Count proc faces only once.
@@ -1254,9 +1285,11 @@ Foam::label Foam::meshRefinement::countHits() const
 
     label nHits = 0;
 
-    forAll(surfaceIndex_, facei)
+    const labelList& surfIndex = surfaceIndex();
+
+    forAll(surfIndex, facei)
     {
-        if (surfaceIndex_[facei] >= 0 && isMasterFace.test(facei))
+        if (surfIndex[facei] >= 0 && isMasterFace.test(facei))
         {
             ++nHits;
         }
@@ -1520,9 +1553,11 @@ Foam::labelList Foam::meshRefinement::intersectedFaces() const
 {
     label nBoundaryFaces = 0;
 
-    forAll(surfaceIndex_, facei)
+    const labelList& surfIndex = surfaceIndex();
+
+    forAll(surfIndex, facei)
     {
-        if (surfaceIndex_[facei] != -1)
+        if (surfIndex[facei] != -1)
         {
             nBoundaryFaces++;
         }
@@ -1531,9 +1566,9 @@ Foam::labelList Foam::meshRefinement::intersectedFaces() const
     labelList surfaceFaces(nBoundaryFaces);
     nBoundaryFaces = 0;
 
-    forAll(surfaceIndex_, facei)
+    forAll(surfIndex, facei)
     {
-        if (surfaceIndex_[facei] != -1)
+        if (surfIndex[facei] != -1)
         {
             surfaceFaces[nBoundaryFaces++] = facei;
         }
@@ -1550,9 +1585,11 @@ Foam::labelList Foam::meshRefinement::intersectedPoints() const
     bitSet isBoundaryPoint(mesh_.nPoints());
     label nBoundaryPoints = 0;
 
-    forAll(surfaceIndex_, facei)
+    const labelList& surfIndex = surfaceIndex();
+
+    forAll(surfIndex, facei)
     {
-        if (surfaceIndex_[facei] != -1)
+        if (surfIndex[facei] != -1)
         {
             const face& f = faces[facei];
 
@@ -3094,6 +3131,9 @@ void Foam::meshRefinement::write
     if (writeFlags && !(writeFlags & NOWRITEREFINEMENT))
     {
         meshCutter_.write();
+
+        // Force calculation before writing
+        (void)surfaceIndex();
         surfaceIndex_.write();
     }
 
@@ -3157,4 +3197,76 @@ void Foam::meshRefinement::writeLevel(const writeType flags)
 //}
 
 
+const Foam::dictionary& Foam::meshRefinement::subDict
+(
+    const dictionary& dict,
+    const word& keyword,
+    const bool noExit
+)
+{
+    if (noExit)
+    {
+        // Find non-recursive with patterns
+        const dictionary::const_searcher finder
+        (
+            dict.csearch
+            (
+                keyword,
+                keyType::REGEX
+            )
+        );
+
+        if (!finder.found())
+        {
+            FatalIOErrorInFunction(dict)
+                << "Entry '" << keyword << "' not found in dictionary "
+                << dict.name();
+            return dictionary::null;
+        }
+        else
+        {
+            return finder.dict();
+        }
+    }
+    else
+    {
+        return dict.subDict(keyword);
+    }
+}
+
+
+Foam::ITstream& Foam::meshRefinement::lookup
+(
+    const dictionary& dict,
+    const word& keyword,
+    const bool noExit
+)
+{
+    if (noExit)
+    {
+        const dictionary::const_searcher finder
+        (
+            dict.csearch(keyword, keyType::REGEX)
+        );
+
+        if (!finder.found())
+        {
+            FatalIOErrorInFunction(dict)
+                << "Entry '" << keyword << "' not found in dictionary "
+                << dict.name();
+            // Fake entry
+            return dict.first()->stream();
+        }
+        else
+        {
+            return finder.ref().stream();
+        }
+    }
+    else
+    {
+        return dict.lookup(keyword);
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
index f67d0722ec..30dcbd8878 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
@@ -165,6 +165,9 @@ private:
         //- All limit-refinement interaction
         const shellSurfaces& limitShells_;
 
+        //- Are we operating in test mode?
+        const bool dryRun_;
+
         //- Refinement engine
         hexRef8 meshCutter_;
 
@@ -819,7 +822,9 @@ public:
             const refinementSurfaces&,
             const refinementFeatures&,
             const shellSurfaces&,   // omnidirectional refinement
-            const shellSurfaces&    // limit refinement
+            const shellSurfaces&,   // limit refinement
+            const labelUList& checkFaces,   // initial faces to check
+            const bool dryRun
         );
 
 
@@ -885,15 +890,9 @@ public:
             }
 
             //- Per start-end edge the index of the surface hit
-            const labelList& surfaceIndex() const
-            {
-                return surfaceIndex_;
-            }
+            const labelList& surfaceIndex() const;
 
-            labelList& surfaceIndex()
-            {
-                return surfaceIndex_;
-            }
+            labelList& surfaceIndex();
 
             //- For faces originating from processor faces store the original
             //  patch
@@ -1531,6 +1530,33 @@ public:
                 const EnumContainer& namedEnum,
                 const wordList& words
             );
+
+            //- Wrapper around dictionary::get which does not exit
+            template<class Type>
+            static Type get
+            (
+                const dictionary& dict,
+                const word& keyword,
+                const bool noExit,
+                enum keyType::option matchOpt = keyType::REGEX,
+                const Type& defaultValue = Zero
+            );
+
+            //- Wrapper around dictionary::subDict which does not exit
+            static const dictionary& subDict
+            (
+                const dictionary& dict,
+                const word& keyword,
+                const bool noExit
+            );
+
+            //- Wrapper around dictionary::lookup which does not exit
+            static ITstream& lookup
+            (
+                const dictionary& dict,
+                const word& keyword,
+                const bool noExit
+            );
 };
 
 
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C
index 01da7090da..4982bba7f9 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C
@@ -1378,7 +1378,7 @@ Foam::label Foam::meshRefinement::markInternalGapRefinement
                 << mesh_.globalData().nTotalFaces() << endl;
 
 
-            transportData::trackData td(surfaceIndex_);
+            transportData::trackData td(surfaceIndex());
 
             FaceCellWave<transportData, transportData::trackData> deltaCalc
             (
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C
index a7cc9e377e..14c293e2fb 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C
@@ -414,7 +414,8 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
                 false,  // report
                 mesh_,
                 motionDict,
-                errorFaces
+                errorFaces,
+                dryRun_
             );
 
             //if (checkEdgeConnectivity)
@@ -898,7 +899,8 @@ Foam::label Foam::meshRefinement::mergeEdgesUndo
                 false,  // report
                 mesh_,
                 motionDict,
-                errorFaces
+                errorFaces,
+                dryRun_
             );
             //if (checkEdgeConnectivity)
             //{
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C
index 2476b0922e..b0eb7cc6eb 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -633,8 +633,8 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
     if (checkCollapse)
     {
-        motionDict.readEntry("minArea", minArea);
-        motionDict.readEntry("maxNonOrtho", maxNonOrtho);
+        minArea = get<scalar>(motionDict, "minArea", dryRun_);
+        maxNonOrtho = get<scalar>(motionDict, "maxNonOrtho", dryRun_);
 
         Info<< "markFacesOnProblemCells :"
             << " Deleting all-anchor surface cells only if"
@@ -1123,7 +1123,14 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
         // Check initial mesh
         Info<< "Checking initial mesh ..." << endl;
         labelHashSet wrongFaces(mesh_.nFaces()/100);
-        motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
+        motionSmoother::checkMesh
+        (
+            false,
+            mesh_,
+            motionDict,
+            wrongFaces,
+            dryRun_
+        );
         const label nInitErrors = returnReduce
         (
             wrongFaces.size(),
@@ -1235,7 +1242,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
             //    nWrongFaces = nNewWrongFaces;
             //}
 
-            scalar minArea(motionDict.get<scalar>("minArea"));
+            scalar minArea(get<scalar>(motionDict, "minArea", dryRun_));
             if (minArea > -SMALL)
             {
                 polyMeshGeometry::checkFaceArea
@@ -1262,7 +1269,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
                 nWrongFaces = nNewWrongFaces;
             }
 
-            scalar minDet(motionDict.get<scalar>("minDeterminant"));
+            scalar minDet(get<scalar>(motionDict, "minDeterminant", dryRun_));
             if (minDet > -1)
             {
                 polyMeshGeometry::checkCellDeterminant
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
index 1315caa1c7..5ed47ac1ef 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
@@ -913,9 +913,11 @@ Foam::labelList Foam::meshRefinement::getRefineCandidateFaces
 
     label nTest = 0;
 
-    forAll(surfaceIndex_, faceI)
+    const labelList& surfIndex = surfaceIndex();
+
+    forAll(surfIndex, faceI)
     {
-        if (surfaceIndex_[faceI] != -1)
+        if (surfIndex[faceI] != -1)
         {
             label own = mesh_.faceOwner()[faceI];
 
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C
index ff7c829fcf..2d7c6a3680 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C
@@ -325,4 +325,35 @@ void Foam::meshRefinement::weightedSum
 }
 
 
+template<class Type>
+Type Foam::meshRefinement::get
+(
+    const dictionary& dict,
+    const word& keyword,
+    const bool noExit,
+    enum keyType::option matchOpt,
+    const Type& defaultValue
+)
+{
+    Type val(defaultValue);
+
+    if
+    (
+       !dict.readEntry
+        (
+            keyword,
+            val,
+            matchOpt,
+            !noExit
+        )
+    )
+    {
+        FatalIOError
+            << "Entry '" << keyword << "' not found in dictionary "
+            << dict.name() << endl;
+    }
+    return val;
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
index 86d322c192..d16b9eb0df 100644
--- a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
+++ b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
@@ -28,6 +28,7 @@ License
 #include "Tuple2.H"
 #include "DynamicField.H"
 #include "featureEdgeMesh.H"
+#include "meshRefinement.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -41,7 +42,17 @@ void Foam::refinementFeatures::read
     {
         const dictionary& dict = featDicts[featI];
 
-        fileName featFileName(dict.lookup("file"));
+        fileName featFileName
+        (
+            meshRefinement::get<fileName>
+            (
+                dict,
+                "file",
+                dryRun_,
+                keyType::REGEX,
+                fileName::null
+            )
+        );
 
 
         // Try reading extendedEdgeMesh first
@@ -233,7 +244,19 @@ void Foam::refinementFeatures::read
         else
         {
             // Look up 'level' for single level
-            levels_[featI] = labelList(1, dict.get<label>("level"));
+            levels_[featI] =
+                labelList
+                (
+                    1,
+                    meshRefinement::get<label>
+                    (
+                        dict,
+                        "level",
+                        dryRun_,
+                        keyType::REGEX,
+                        0
+                    )
+                );
             distances_[featI] = scalarField(1, Zero);
         }
 
@@ -439,14 +462,16 @@ Foam::refinementFeatures::regionEdgeTrees() const
 Foam::refinementFeatures::refinementFeatures
 (
     const objectRegistry& io,
-    const PtrList<dictionary>& featDicts
+    const PtrList<dictionary>& featDicts,
+    const bool dryRun
 )
 :
     PtrList<extendedFeatureEdgeMesh>(featDicts.size()),
     distances_(featDicts.size()),
     levels_(featDicts.size()),
     edgeTrees_(featDicts.size()),
-    pointTrees_(featDicts.size())
+    pointTrees_(featDicts.size()),
+    dryRun_(dryRun)
 {
     // Read features
     read(io, featDicts);
@@ -531,6 +556,76 @@ Foam::refinementFeatures::refinementFeatures
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::refinementFeatures::checkSizes
+(
+    const scalar maxRatio,
+    const boundBox& meshBb,
+    const bool report,
+    Ostream& os
+) const
+{
+    if (report)
+    {
+        os<< "Checking for size." << endl;
+    }
+
+    bool hasError = false;
+
+    forAll(*this, i)
+    {
+        const extendedFeatureEdgeMesh& em = operator[](i);
+        const boundBox bb(em.points(), true);
+
+        for (label j = i+1; j < size(); j++)
+        {
+            const extendedFeatureEdgeMesh& em2 = operator[](j);
+            const boundBox bb2(em2.points(), true);
+
+            scalar ratio = bb.mag()/bb2.mag();
+
+            if (ratio > maxRatio || ratio < 1.0/maxRatio)
+            {
+                hasError = true;
+
+                if (report)
+                {
+                    os  << "    " << em.name()
+                        << " bounds differ from " << em2.name()
+                        << " by more than a factor 100:" << nl
+                        << "        bounding box : " << bb << nl
+                        << "        bounding box : " << bb2
+                        << endl;
+                }
+            }
+        }
+    }
+
+    forAll(*this, i)
+    {
+        const extendedFeatureEdgeMesh& em = operator[](i);
+        const boundBox bb(em.points(), true);
+        if (!meshBb.contains(bb))
+        {
+            if (report)
+            {
+                os  << "    " << em.name()
+                    << " bounds not fully contained in mesh"<< nl
+                    << "        bounding box      : " << bb << nl
+                    << "        mesh bounding box : " << meshBb
+                    << endl;
+            }
+        }
+    }
+
+    if (report)
+    {
+        os<< endl;
+    }
+
+    return returnReduce(hasError, orOp<bool>());
+}
+
+
 void Foam::refinementFeatures::findNearestEdge
 (
     const pointField& samples,
diff --git a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H
index 6147c81bcf..b8b46ac8f3 100644
--- a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H
+++ b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H
@@ -67,6 +67,9 @@ class refinementFeatures
         //- Features points
         PtrList<indexedOctree<treeDataPoint>> pointTrees_;
 
+        //- Dry-run mode?
+        const bool dryRun_;
+
         //- Region edge trees (demand driven)
         mutable autoPtr<PtrList<indexedOctree<treeDataEdge>>>
             regionEdgeTreesPtr_;
@@ -112,7 +115,8 @@ public:
         refinementFeatures
         (
             const objectRegistry& io,
-            const PtrList<dictionary>& featDicts
+            const PtrList<dictionary>& featDicts,
+            const bool dryRun = false
         );
 
 
@@ -189,6 +193,16 @@ public:
                 labelList& maxLevel
             ) const;
 
+        // Other
+
+            //- Check sizes - return true if error and stream to os
+            bool checkSizes
+            (
+                const scalar maxRatio,
+                const boundBox& meshBb,
+                const bool report,
+                Ostream& os
+            ) const;
 };
 
 
diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
index c300ceaa40..61d5c9607c 100644
--- a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -32,6 +32,8 @@ License
 #include "searchableSurfacesQueries.H"
 #include "UPtrList.H"
 #include "volumeType.H"
+// For dictionary::get wrapper
+#include "meshRefinement.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -124,14 +126,16 @@ Foam::refinementSurfaces::refinementSurfaces
 (
     const searchableSurfaces& allGeometry,
     const dictionary& surfacesDict,
-    const label gapLevelIncrement
+    const label gapLevelIncrement,
+    const bool dryRun
 )
 :
     allGeometry_(allGeometry),
     surfaces_(surfacesDict.size()),
     names_(surfacesDict.size()),
     surfZones_(surfacesDict.size()),
-    regionOffset_(surfacesDict.size())
+    regionOffset_(surfacesDict.size()),
+    dryRun_(dryRun)
 {
     // Wildcard specification : loop over all surface, all regions
     // and try to find a match.
@@ -195,7 +199,18 @@ Foam::refinementSurfaces::refinementSurfaces
             names_[surfI] = geomName;
             surfaces_[surfI] = geomI;
 
-            const labelPair refLevel(dict.lookup("level"));
+            const labelPair refLevel
+            (
+                meshRefinement::get<labelPair>
+                (
+                    dict,
+                    "level",
+                    dryRun_,
+                    keyType::REGEX,
+                    labelPair(0, 0)
+                )
+            );
+
             globalMinLevel[surfI] = refLevel[0];
             globalMaxLevel[surfI] = refLevel[1];
             globalLevelIncr[surfI] = dict.lookupOrDefault
@@ -281,7 +296,18 @@ Foam::refinementSurfaces::refinementSurfaces
                             regionNames[regionI]
                         );
 
-                        const labelPair refLevel(regionDict.lookup("level"));
+                        const labelPair refLevel
+                        (
+                            meshRefinement::get<labelPair>
+                            (
+                                regionDict,
+                                "level",
+                                dryRun_,
+                                keyType::REGEX,
+                                labelPair(0, 0)
+                            )
+                        );
+
 
                         regionMinLevel[surfI].insert(regionI, refLevel[0]);
                         regionMaxLevel[surfI].insert(regionI, refLevel[1]);
@@ -476,7 +502,8 @@ Foam::refinementSurfaces::refinementSurfaces
     const labelList& maxLevel,
     const labelList& gapLevel,
     const scalarField& perpendicularAngle,
-    PtrList<dictionary>& patchInfo
+    PtrList<dictionary>& patchInfo,
+    const bool dryRun
 )
 :
     allGeometry_(allGeometry),
@@ -488,7 +515,8 @@ Foam::refinementSurfaces::refinementSurfaces
     maxLevel_(maxLevel),
     gapLevel_(gapLevel),
     perpendicularAngle_(perpendicularAngle),
-    patchInfo_(patchInfo.size())
+    patchInfo_(patchInfo.size()),
+    dryRun_(dryRun)
 {
     forAll(patchInfo_, pI)
     {
diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H
index d3939c763e..c716f7ef5a 100644
--- a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -98,6 +98,9 @@ class refinementSurfaces
         //- From global region number to patchType
         PtrList<dictionary> patchInfo_;
 
+        //- Are we operating in test mode?
+        const bool dryRun_;
+
 
     // Private Member Functions
 
@@ -127,7 +130,8 @@ public:
         (
             const searchableSurfaces& allGeometry,
             const dictionary&,
-            const label gapLevelIncrement
+            const label gapLevelIncrement,
+            const bool dryRun
         );
 
         //- Construct from components
@@ -142,7 +146,8 @@ public:
             const labelList& maxLevel,
             const labelList& gapLevel,
             const scalarField& perpendicularAngle,
-            PtrList<dictionary>& patchInfo
+            PtrList<dictionary>& patchInfo,
+            const bool dryRun
         );
 
 
diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
index 5739217fe2..35c83377c5 100644
--- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
+++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
@@ -569,10 +569,12 @@ void Foam::shellSurfaces::findLevel
 Foam::shellSurfaces::shellSurfaces
 (
     const searchableSurfaces& allGeometry,
-    const dictionary& shellsDict
+    const dictionary& shellsDict,
+    const bool dryRun
 )
 :
-    allGeometry_(allGeometry)
+    allGeometry_(allGeometry),
+    dryRun_(dryRun)
 {
     // Wilcard specification : loop over all surfaces and try to find a match.
 
diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H
index 2390a98e95..f5bc6f7b63 100644
--- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H
+++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H
@@ -74,6 +74,9 @@ private:
         //- Reference to all geometry.
         const searchableSurfaces& allGeometry_;
 
+        //- Dry-run mode?
+        const bool dryRun_;
+
         //- Indices of surfaces that are shells
         labelList shells_;
 
@@ -173,7 +176,8 @@ public:
         shellSurfaces
         (
             const searchableSurfaces& allGeometry,
-            const dictionary& shellsDict
+            const dictionary& shellsDict,
+            const bool dryRun = false
         );
 
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.C
index 4e5f08bfb5..00a76b59aa 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.C
@@ -96,12 +96,13 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
 Foam::layerParameters::layerParameters
 (
     const dictionary& dict,
-    const polyBoundaryMesh& boundaryMesh
+    const polyBoundaryMesh& boundaryMesh,
+    const bool dryRun
 )
 :
     dict_(dict),
     numLayers_(boundaryMesh.size(), -1),
-    relativeSizes_(dict.lookup("relativeSizes")),
+    relativeSizes_(meshRefinement::get<bool>(dict, "relativeSizes", dryRun)),
     layerSpec_(ILLEGAL),
     firstLayerThickness_(boundaryMesh.size(), -123),
     finalLayerThickness_(boundaryMesh.size(), -123),
@@ -110,9 +111,9 @@ Foam::layerParameters::layerParameters
     minThickness_
     (
         boundaryMesh.size(),
-        dict.get<scalar>("minThickness")
+        meshRefinement::get<scalar>(dict, "minThickness", dryRun)
     ),
-    featureAngle_(dict.get<scalar>("featureAngle")),
+    featureAngle_(meshRefinement::get<scalar>(dict, "featureAngle", dryRun)),
     mergePatchFacesAngle_
     (
         dict.lookupOrDefault<scalar>
@@ -125,16 +126,16 @@ Foam::layerParameters::layerParameters
     (
         dict.lookupOrDefault("concaveAngle", defaultConcaveAngle)
     ),
-    nGrow_(dict.get<label>("nGrow")),
+    nGrow_(meshRefinement::get<label>(dict, "nGrow", dryRun)),
     maxFaceThicknessRatio_
     (
-        dict.get<scalar>("maxFaceThicknessRatio")
+        meshRefinement::get<scalar>(dict, "maxFaceThicknessRatio", dryRun)
     ),
     nBufferCellsNoExtrude_
     (
-        dict.get<label>("nBufferCellsNoExtrude")
+        meshRefinement::get<label>(dict, "nBufferCellsNoExtrude", dryRun)
     ),
-    nLayerIter_(dict.get<label>("nLayerIter")),
+    nLayerIter_(meshRefinement::get<label>(dict, "nLayerIter", dryRun)),
     nRelaxedIter_(labelMax),
     additionalReporting_(dict.lookupOrDefault("additionalReporting", false)),
     meshShrinker_
@@ -144,7 +145,8 @@ Foam::layerParameters::layerParameters
             "meshShrinker",
             medialAxisMeshMover::typeName
         )
-    )
+    ),
+    dryRun_(dryRun)
 {
     // Detect layer specification mode
 
@@ -255,7 +257,12 @@ Foam::layerParameters::layerParameters
     }
 
 
-    const dictionary& layersDict = dict.subDict("layers");
+    const dictionary& layersDict = meshRefinement::subDict
+    (
+        dict,
+        "layers",
+        dryRun
+    );
 
     for (const entry& dEntry : layersDict)
     {
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.H
index 34a9e813ec..e5d586e6a6 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.H
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.H
@@ -132,6 +132,9 @@ private:
 
         word meshShrinker_;
 
+        //- In dry-run mode?
+        const bool dryRun_;
+
 
     // Private Member Functions
 
@@ -155,7 +158,12 @@ public:
     // Constructors
 
         //- Construct from dictionary
-        layerParameters(const dictionary& dict, const polyBoundaryMesh&);
+        layerParameters
+        (
+            const dictionary& dict,
+            const polyBoundaryMesh&,
+            const bool dryRun = false
+        );
 
 
     // Member Functions
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C
index c63345f993..ab84e69831 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C
@@ -29,14 +29,28 @@ License
 #include "globalIndex.H"
 #include "Tuple2.H"
 #include "wallPolyPatch.H"
+#include "meshRefinement.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::refinementParameters::refinementParameters(const dictionary& dict)
+Foam::refinementParameters::refinementParameters
+(
+    const dictionary& dict,
+    const bool dryRun
+)
 :
-    maxGlobalCells_(dict.get<label>("maxGlobalCells")),
-    maxLocalCells_(dict.get<label>("maxLocalCells")),
-    minRefineCells_(dict.get<label>("minRefinementCells")),
+    maxGlobalCells_
+    (
+        meshRefinement::get<label>(dict, "maxGlobalCells", dryRun)
+    ),
+    maxLocalCells_
+    (
+        meshRefinement::get<label>(dict, "maxLocalCells", dryRun)
+    ),
+    minRefineCells_
+    (
+        meshRefinement::get<label>(dict, "minRefinementCells", dryRun)
+    ),
     planarAngle_
     (
         dict.lookupOrDefault
@@ -45,7 +59,10 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
             dict.get<scalar>("resolveFeatureAngle")
         )
     ),
-    nBufferLayers_(dict.get<label>("nCellsBetweenLevels")),
+    nBufferLayers_
+    (
+        meshRefinement::get<label>(dict, "nCellsBetweenLevels", dryRun)
+    ),
     locationsOutsideMesh_
     (
         dict.lookupOrDefault
@@ -55,7 +72,10 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
         )
     ),
     faceZoneControls_(dict.subOrEmptyDict("faceZoneControls")),
-    allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
+    allowFreeStandingZoneFaces_
+    (
+        meshRefinement::get<bool>(dict, "allowFreeStandingZoneFaces", dryRun)
+    ),
     useTopologicalSnapDetection_
     (
         dict.lookupOrDefault("useTopologicalSnapDetection", true)
@@ -69,9 +89,11 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
     (
         dict.lookupOrDefault<Switch>("interfaceRefine", true)
     ),
-    nErodeCellZone_(dict.lookupOrDefault<label>("nCellZoneErodeIter", 0))
+    nErodeCellZone_(dict.lookupOrDefault<label>("nCellZoneErodeIter", 0)),
+    dryRun_(dryRun)
 {
     point locationInMesh;
+    List<Tuple2<point, word>> pointsToZone;
     if (dict.readIfPresent("locationInMesh", locationInMesh))
     {
         locationsInMesh_.append(locationInMesh);
@@ -84,7 +106,7 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
                 << exit(FatalIOError);
         }
     }
-    else
+    else if (dict.readIfPresent("locationsInMesh", pointsToZone))
     {
         List<Tuple2<point, word>> pointsToZone(dict.lookup("locationsInMesh"));
         label nZones = locationsInMesh_.size();
@@ -102,9 +124,18 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
             nZones++;
         }
     }
+    else
+    {
+        IOWarningInFunction(dict)
+            << "No 'locationInMesh' or 'locationsInMesh' provided"
+            << endl;
+    }
 
 
-    const scalar featAngle(dict.get<scalar>("resolveFeatureAngle"));
+    const scalar featAngle
+    (
+        meshRefinement::get<scalar>(dict, "resolveFeatureAngle", dryRun)
+    );
 
     if (featAngle < 0 || featAngle > 180)
     {
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.H
index a6475c00d4..29112d02ef 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.H
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.H
@@ -70,7 +70,7 @@ class refinementParameters
         scalar curvature_;
 
         //- Planarity criterion
-        scalar planarAngle_;
+        const scalar planarAngle_;
 
         //- Number of layers between different refinement levels
         const label nBufferLayers_;
@@ -93,20 +93,22 @@ class refinementParameters
 
         //- FaceZone faces allowed which have owner and neighbour in same
         //  cellZone?
-        Switch allowFreeStandingZoneFaces_;
+        const Switch allowFreeStandingZoneFaces_;
 
         //- Use old topology based problem-cell removal (cells with 8 points
         //  on surface)
-        Switch useTopologicalSnapDetection_;
+        const Switch useTopologicalSnapDetection_;
 
         //- Allowed load unbalance
-        scalar maxLoadUnbalance_;
+        const scalar maxLoadUnbalance_;
 
-        Switch handleSnapProblems_;
+        const Switch handleSnapProblems_;
 
-        Switch interfaceRefine_;
+        const Switch interfaceRefine_;
 
-        label nErodeCellZone_;
+        const label nErodeCellZone_;
+
+        const bool dryRun_;
 
 
     // Private Member Functions
@@ -123,7 +125,7 @@ public:
     // Constructors
 
         //- Construct from dictionary - new syntax
-        refinementParameters(const dictionary& dict);
+        refinementParameters(const dictionary& dict, const bool dryRun = false);
 
 
     // Member Functions
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.C
index 3f19065202..7c9693dc59 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.C
@@ -24,17 +24,30 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "snapParameters.H"
+#include "meshRefinement.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from dictionary
-Foam::snapParameters::snapParameters(const dictionary& dict)
+Foam::snapParameters::snapParameters(const dictionary& dict, const bool dryRun)
 :
-    nSmoothPatch_(dict.get<label>("nSmoothPatch")),
+    nSmoothPatch_
+    (
+        meshRefinement::get<label>(dict, "nSmoothPatch", dryRun)
+    ),
     nSmoothInternal_(dict.lookupOrDefault("nSmoothInternal", 0)),
-    snapTol_(dict.get<scalar>("tolerance")),
-    nSmoothDispl_(dict.get<label>("nSolveIter")),
-    nSnap_(dict.get<label>("nRelaxIter")),
+    snapTol_
+    (
+        meshRefinement::get<scalar>(dict, "tolerance", dryRun)
+    ),
+    nSmoothDispl_
+    (
+        meshRefinement::get<label>(dict, "nSolveIter", dryRun)
+    ),
+    nSnap_
+    (
+        meshRefinement::get<label>(dict, "nRelaxIter", dryRun)
+    ),
     nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1)),
     explicitFeatureSnap_(dict.lookupOrDefault("explicitFeatureSnap", true)),
     implicitFeatureSnap_(dict.lookupOrDefault("implicitFeatureSnap", false)),
@@ -60,7 +73,8 @@ Foam::snapParameters::snapParameters(const dictionary& dict)
         dict.lookupOrDefault("nFaceSplitInterval", labelMin)
     ),
     concaveAngle_(dict.lookupOrDefault("concaveAngle", 45)),
-    minAreaRatio_(dict.lookupOrDefault("minAreaRatio", 0.3))
+    minAreaRatio_(dict.lookupOrDefault("minAreaRatio", 0.3)),
+    dryRun_(dryRun)
 {}
 
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.H
index 18aca4d8c7..1f72e3b8bd 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.H
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.H
@@ -86,13 +86,16 @@ class snapParameters
 
 
         //- How often needs face splitting be run
-        label nFaceSplitInterval_;
+        const label nFaceSplitInterval_;
 
         //- When is angle too concave too split
-        scalar concaveAngle_;
+        const scalar concaveAngle_;
 
         //- When is face-split not sufficiently diagonal
-        scalar minAreaRatio_;
+        const scalar minAreaRatio_;
+
+        //- in dry-run mode
+        const bool dryRun_;
 
 
     // Private Member Functions
@@ -109,7 +112,7 @@ public:
     // Constructors
 
         //- Construct from dictionary
-        snapParameters(const dictionary& dict);
+        snapParameters(const dictionary& dict, const bool dryRun = false);
 
 
     // Member Functions
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
index ebae51274f..f794e54437 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
@@ -2648,7 +2648,8 @@ Foam::label Foam::snappyLayerDriver::checkAndUnmark
         meshQualityDict,
         identity(newMesh.nFaces()),
         baffles,
-        wrongFaces
+        wrongFaces,
+        false           // dryRun_
     );
     Info<< "Detected " << returnReduce(wrongFaces.size(), sumOp<label>())
         << " illegal faces"
@@ -3215,12 +3216,14 @@ Foam::snappyLayerDriver::snappyLayerDriver
 (
     meshRefinement& meshRefiner,
     const labelList& globalToMasterPatch,
-    const labelList& globalToSlavePatch
+    const labelList& globalToSlavePatch,
+    const bool dryRun
 )
 :
     meshRefiner_(meshRefiner),
     globalToMasterPatch_(globalToMasterPatch),
-    globalToSlavePatch_(globalToSlavePatch)
+    globalToSlavePatch_(globalToSlavePatch),
+    dryRun_(dryRun)
 {}
 
 
@@ -3894,6 +3897,40 @@ void Foam::snappyLayerDriver::addLayers
                 internalBaffles,
                 displacement
             );
+
+
+            if (dryRun_)
+            {
+                string errorMsg(FatalError.message());
+                string IOerrorMsg(FatalIOError.message());
+
+                if (errorMsg.size() || IOerrorMsg.size())
+                {
+                    //errorMsg = "[dryRun] " + errorMsg;
+                    //errorMsg.replaceAll("\n", "\n[dryRun] ");
+                    //IOerrorMsg = "[dryRun] " + IOerrorMsg;
+                    //IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
+
+                    IOWarningInFunction(combinedDict)
+                        << nl
+                        << "Missing/incorrect required dictionary entries:"
+                        << nl << nl
+                        << IOerrorMsg.c_str() << nl
+                        << errorMsg.c_str() << nl << nl
+                        << "Exiting dry-run" << nl << endl;
+
+                    if (Pstream::parRun())
+                    {
+                        Perr<< "\nFOAM parallel run exiting\n" << endl;
+                        Pstream::exit(0);
+                    }
+                    else
+                    {
+                        Perr<< "\nFOAM exiting\n" << endl;
+                        ::exit(0);
+                    }
+                }
+            }
         }
 
 
@@ -4691,7 +4728,7 @@ void Foam::snappyLayerDriver::doLayers
         // Check initial mesh
         Info<< "Checking initial mesh ..." << endl;
         labelHashSet wrongFaces(mesh.nFaces()/100);
-        motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
+        motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces, dryRun_);
         const label nInitErrors = returnReduce
         (
             wrongFaces.size(),
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H
index 9c663cf70b..c0df3db8a2 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H
@@ -106,6 +106,9 @@ private:
         //- From surface region to patch
         const labelList globalToSlavePatch_;
 
+        //- Are we operating in test mode?
+        const bool dryRun_;
+
 
     // Private Member Functions
 
@@ -618,7 +621,8 @@ public:
         (
             meshRefinement& meshRefiner,
             const labelList& globalToMasterPatch,
-            const labelList& globalToSlavePatch
+            const labelList& globalToSlavePatch,
+            const bool dryRun = false
         );
 
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
index 8bd94a084e..f05b6fe659 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
@@ -43,6 +43,7 @@ License
 #include "searchableSurfaces.H"
 #include "fvMeshSubset.H"
 #include "interpolationTable.H"
+#include "snappyVoxelMeshDriver.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -61,7 +62,8 @@ Foam::snappyRefineDriver::snappyRefineDriver
     fvMeshDistribute& distributor,
     const labelUList& globalToMasterPatch,
     const labelUList& globalToSlavePatch,
-    const writer<scalar>& setFormatter
+    const writer<scalar>& setFormatter,
+    const bool dryRun
 )
 :
     meshRefiner_(meshRefiner),
@@ -69,7 +71,8 @@ Foam::snappyRefineDriver::snappyRefineDriver
     distributor_(distributor),
     globalToMasterPatch_(globalToMasterPatch),
     globalToSlavePatch_(globalToSlavePatch),
-    setFormatter_(setFormatter)
+    setFormatter_(setFormatter),
+    dryRun_(dryRun)
 {}
 
 
@@ -82,6 +85,11 @@ Foam::label Foam::snappyRefineDriver::featureEdgeRefine
     const label minRefine
 )
 {
+    if (dryRun_)
+    {
+        return 0;
+    }
+
     if (refineParams.minRefineCells() == -1)
     {
         // Special setting to be able to restart shm on meshes with inconsistent
@@ -199,6 +207,11 @@ Foam::label Foam::snappyRefineDriver::smallFeatureRefine
     const label maxIter
 )
 {
+    if (dryRun_)
+    {
+        return 0;
+    }
+
     if (refineParams.minRefineCells() == -1)
     {
         // Special setting to be able to restart shm on meshes with inconsistent
@@ -209,7 +222,6 @@ Foam::label Foam::snappyRefineDriver::smallFeatureRefine
     addProfiling(feature, "snappyHexMesh::refine::smallFeature");
     const fvMesh& mesh = meshRefiner_.mesh();
 
-
     label iter = 0;
 
     // See if any surface has an extendedGapLevel
@@ -329,6 +341,11 @@ Foam::label Foam::snappyRefineDriver::surfaceOnlyRefine
     const label maxIter
 )
 {
+    if (dryRun_)
+    {
+        return 0;
+    }
+
     if (refineParams.minRefineCells() == -1)
     {
         // Special setting to be able to restart shm on meshes with inconsistent
@@ -460,6 +477,11 @@ Foam::label Foam::snappyRefineDriver::gapOnlyRefine
     const label maxIter
 )
 {
+    if (dryRun_)
+    {
+        return 0;
+    }
+
     if (refineParams.minRefineCells() == -1)
     {
         // Special setting to be able to restart shm on meshes with inconsistent
@@ -706,6 +728,11 @@ Foam::label Foam::snappyRefineDriver::bigGapOnlyRefine
         return 0;
     }
 
+    if (dryRun_)
+    {
+        return 0;
+    }
+
     const fvMesh& mesh = meshRefiner_.mesh();
 
     label iter = 0;
@@ -855,6 +882,11 @@ Foam::label Foam::snappyRefineDriver::danglingCellRefine
         return 0;
     }
 
+    if (dryRun_)
+    {
+        return 0;
+    }
+
     addProfiling(dangling, "snappyHexMesh::refine::danglingCell");
     const fvMesh& mesh = meshRefiner_.mesh();
 
@@ -1008,6 +1040,11 @@ Foam::label Foam::snappyRefineDriver::refinementInterfaceRefine
         return 0;
     }
 
+    if (dryRun_)
+    {
+        return 0;
+    }
+
     addProfiling(interface, "snappyHexMesh::refine::transition");
     const fvMesh& mesh = meshRefiner_.mesh();
 
@@ -1355,6 +1392,11 @@ void Foam::snappyRefineDriver::removeInsideCells
         return;
     }
 
+    if (dryRun_)
+    {
+        return;
+    }
+
     Info<< nl
         << "Removing mesh beyond surface intersections" << nl
         << "------------------------------------------" << nl
@@ -1420,6 +1462,11 @@ Foam::label Foam::snappyRefineDriver::shellRefine
     const label maxIter
 )
 {
+    if (dryRun_)
+    {
+        return 0;
+    }
+
     if (refineParams.minRefineCells() == -1)
     {
         // Special setting to be able to restart shm on meshes with inconsistent
@@ -1601,6 +1648,11 @@ Foam::label Foam::snappyRefineDriver::directionalShellRefine
     const label maxIter
 )
 {
+    if (dryRun_)
+    {
+        return 0;
+    }
+
     addProfiling(shell, "snappyHexMesh::refine::directionalShell");
     const fvMesh& mesh = meshRefiner_.mesh();
     const shellSurfaces& shells = meshRefiner_.shells();
@@ -2323,6 +2375,11 @@ void Foam::snappyRefineDriver::baffleAndSplitMesh
     const dictionary& motionDict
 )
 {
+    if (dryRun_)
+    {
+        return;
+    }
+
     addProfiling(split, "snappyHexMesh::refine::splitting");
     Info<< nl
         << "Splitting mesh at surface intersections" << nl
@@ -2454,6 +2511,11 @@ void Foam::snappyRefineDriver::splitAndMergeBaffles
     const dictionary& motionDict
 )
 {
+    if (dryRun_)
+    {
+        return;
+    }
+
     Info<< nl
         << "Handling cells with snap problems" << nl
         << "---------------------------------" << nl
@@ -2656,6 +2718,11 @@ void Foam::snappyRefineDriver::mergePatchFaces
     const dictionary& motionDict
 )
 {
+    if (dryRun_)
+    {
+        return;
+    }
+
     addProfiling(merge, "snappyHexMesh::refine::merge");
     Info<< nl
         << "Merge refined boundary faces" << nl
@@ -2719,9 +2786,29 @@ void Foam::snappyRefineDriver::doRefine
 
     const fvMesh& mesh = meshRefiner_.mesh();
 
+
     // Check that all the keep points are inside the mesh.
     refineParams.findCells(true, mesh, refineParams.locationsInMesh());
 
+    // Check that all the keep points are inside the mesh.
+    if (dryRun_)
+    {
+        refineParams.findCells(true, mesh, refineParams.locationsOutsideMesh());
+    }
+
+    // Estimate cell sizes
+    if (dryRun_)
+    {
+        snappyVoxelMeshDriver voxelDriver
+        (
+            meshRefiner_,
+            globalToMasterPatch_,
+            globalToSlavePatch_
+        );
+        voxelDriver.doRefine(refineParams);
+    }
+
+
     // Refine around feature edges
     featureEdgeRefine
     (
@@ -2892,7 +2979,7 @@ void Foam::snappyRefineDriver::doRefine
     }
 
 
-    if (Pstream::parRun())
+    if (!dryRun_ && Pstream::parRun())
     {
         Info<< nl
             << "Doing final balancing" << nl
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H
index bab0d95ed7..fdb9ca7dc0 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H
@@ -39,6 +39,8 @@ SourceFiles
 #include "scalarField.H"
 #include "Tuple2.H"
 #include "writer.H"
+#include "DynamicList.H"
+#include "labelVector.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,9 +82,27 @@ class snappyRefineDriver
         //- How to write lines. Used e.g. when writing leak-paths
         const writer<scalar>& setFormatter_;
 
+        //- Are we operating in test mode?
+        const bool dryRun_;
+
 
     // Private Member Functions
 
+        void addNeighbours
+        (
+            const labelVector& n,
+            const labelList& cellLevel,
+            const labelVector& voxel,
+            const label voxeli,
+            DynamicList<labelVector>& front
+        ) const;
+
+        //- Rough estimate of cell size and cell count
+        void estimateCellSizeAndCount
+        (
+            const refinementParameters& refineParams
+        ) const;
+
         //- Refine all cells pierced by explicit feature edges
         label featureEdgeRefine
         (
@@ -225,7 +245,8 @@ public:
             fvMeshDistribute& distributor,
             const labelUList& globalToMasterPatch,
             const labelUList& globalToSlavePatch,
-            const writer<scalar>& setFormatter
+            const writer<scalar>& setFormatter,
+            const bool dryRun = false
         );
 
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
index b90d4b3327..c7a7f043b9 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
@@ -757,12 +757,14 @@ Foam::snappySnapDriver::snappySnapDriver
 (
     meshRefinement& meshRefiner,
     const labelList& globalToMasterPatch,
-    const labelList& globalToSlavePatch
+    const labelList& globalToSlavePatch,
+    const bool dryRun
 )
 :
     meshRefiner_(meshRefiner),
     globalToMasterPatch_(globalToMasterPatch),
-    globalToSlavePatch_(globalToSlavePatch)
+    globalToSlavePatch_(globalToSlavePatch),
+    dryRun_(dryRun)
 {}
 
 
@@ -2606,7 +2608,11 @@ void Foam::snappySnapDriver::doSnap
     if (snapParams.nFeatureSnap() > 0)
     {
         doFeatures = true;
-        nFeatIter = snapParams.nFeatureSnap();
+
+        if (!dryRun_)
+        {
+            nFeatIter = snapParams.nFeatureSnap();
+        }
 
         Info<< "Snapping to features in " << nFeatIter
             << " iterations ..." << endl;
@@ -2652,7 +2658,8 @@ void Foam::snappySnapDriver::doSnap
                     pointMesh::New(mesh),
                     adaptPatchIDs
                 ),
-                motionDict
+                motionDict,
+                dryRun_
             )
         );
 
@@ -2660,7 +2667,7 @@ void Foam::snappySnapDriver::doSnap
         // Check initial mesh
         Info<< "Checking initial mesh ..." << endl;
         labelHashSet wrongFaces(mesh.nFaces()/100);
-        motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
+        motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces, dryRun_);
         const label nInitErrors = returnReduce
         (
             wrongFaces.size(),
@@ -2934,7 +2941,8 @@ void Foam::snappySnapDriver::doSnap
                             pointMesh::New(mesh),
                             adaptPatchIDs
                         ),
-                        motionDict
+                        motionDict,
+                        dryRun_
                     )
                 );
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H
index 7d7f9c4056..42783d354c 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H
@@ -67,6 +67,9 @@ class snappySnapDriver
         //- From global surface region to slave side patch
         const labelList globalToSlavePatch_;
 
+        //- Are we operating in test mode?
+        const bool dryRun_;
+
 
     // Private Member Functions
 
@@ -670,7 +673,8 @@ public:
         (
             meshRefinement& meshRefiner,
             const labelList& globalToMasterPatch,
-            const labelList& globalToSlavePatch
+            const labelList& globalToSlavePatch,
+            const bool dryRun = false
         );
 
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.C
new file mode 100644
index 0000000000..126444be2b
--- /dev/null
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.C
@@ -0,0 +1,594 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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 "snappyVoxelMeshDriver.H"
+#include "meshRefinement.H"
+#include "fvMesh.H"
+#include "Time.H"
+#include "refinementParameters.H"
+#include "refinementSurfaces.H"
+#include "refinementFeatures.H"
+#include "shellSurfaces.H"
+#include "searchableSurfaces.H"
+#include "voxelMeshSearch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(snappyVoxelMeshDriver, 0);
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::snappyVoxelMeshDriver::addNeighbours
+(
+    const labelList& cellLevel,
+    const labelVector& voxel,
+    const label voxeli,
+    DynamicList<labelVector>& front
+) const
+{
+    const labelVector off(voxelMeshSearch::offset(n_));
+
+    for (direction dir = 0; dir < 3; ++dir)
+    {
+        if (voxel[dir] > 0)
+        {
+            labelVector left(voxel);
+            left[dir] -= 1;
+            if (cellLevel[voxeli-off[dir]] == -1)
+            {
+                front.append(left);
+            }
+        }
+        if (voxel[dir] < n_[dir]-1)
+        {
+            labelVector right(voxel);
+            right[dir] += 1;
+            if (cellLevel[voxeli+off[dir]] == -1)
+            {
+                front.append(right);
+            }
+        }
+    }
+}
+
+
+// Insert cell level for all volume refinement
+Foam::tmp<Foam::pointField> Foam::snappyVoxelMeshDriver::voxelCentres() const
+{
+    tmp<pointField> tcc(tmp<pointField>::New(n_.x()*n_.y()*n_.z()));
+    pointField& cc = tcc.ref();
+
+    const labelVector off(voxelMeshSearch::offset(n_));
+    label voxeli = voxelMeshSearch::index(n_, labelVector(0, 0, 0));
+    for (label k = 0; k < n_[2]; k++)
+    {
+        const label start1 = voxeli;
+        for (label j = 0; j < n_[1]; j++)
+        {
+            const label start0 = voxeli;
+            for (label i = 0; i < n_[0]; i++)
+            {
+                const labelVector voxel(i, j, k);
+                cc[voxeli] = voxelMeshSearch::centre(bb_, n_, voxel);
+                voxeli += off[0];
+            }
+            voxeli = start0 + off[1];
+        }
+        voxeli = start1 + off[2];
+    }
+    return tcc;
+}
+
+
+void Foam::snappyVoxelMeshDriver::isInside
+(
+    const pointField& cc,
+    boolList& isVoxelInMesh
+) const
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+    isVoxelInMesh.setSize(cc.size());
+    if (isVoxelInMesh.size() < mesh.globalData().nTotalCells())
+    {
+        forAll(cc, voxeli)
+        {
+            const label celli = mesh.findCell
+            (
+                cc[voxeli],
+                polyMesh::FACE_PLANES
+            );
+            isVoxelInMesh[voxeli] = (celli != -1);
+        }
+        Pstream::listCombineGather(isVoxelInMesh, orEqOp<bool>());
+    }
+    else
+    {
+        const cellList& cells = mesh.cells();
+        const faceList& faces = mesh.faces();
+        const pointField& points = mesh.points();
+
+        for (label celli = 0; celli < mesh.nCells(); celli++)
+        {
+            const cell& cFaces = cells[celli];
+            boundBox cellBb(boundBox::invertedBox);
+            forAll(cFaces, cFacei)
+            {
+                const face& f = faces[cFaces[cFacei]];
+                forAll(f, fp)
+                {
+                    cellBb.add(points[f[fp]]);
+                }
+            }
+            voxelMeshSearch::fill
+            (
+                isVoxelInMesh,
+                bb_,
+                n_,
+                cellBb,
+                1,
+                orEqOp<bool>()
+            );
+        }
+        Pstream::listCombineGather(isVoxelInMesh, orEqOp<bool>());
+    }
+}
+
+
+void Foam::snappyVoxelMeshDriver::markSurfaceRefinement
+(
+    labelList& voxelLevel,
+    labelList& globalRegion
+) const
+{
+    // Insert cell level for all refinementSurfaces
+
+    const refinementSurfaces& s = meshRefiner_.surfaces();
+    forAll(s.surfaces(), surfi)
+    {
+        label geomi = s.surfaces()[surfi];
+        const searchableSurface& geom = s.geometry()[geomi];
+        //Pout<< "Geometry:" << s.names()[surfi] << endl;
+        if (isA<triSurface>(geom))
+        {
+            const triSurface& ts = refCast<const triSurface>(geom);
+            const pointField& points = ts.points();
+
+            forAll(ts, trii)
+            {
+                label regioni = ts[trii].region();
+                label globalRegioni = s.regionOffset()[surfi]+regioni;
+                const boundBox triBb(points, ts[trii], false);
+
+                // Fill cellLevel
+                label level = s.minLevel()[globalRegioni];
+                voxelMeshSearch::fill
+                (
+                    voxelLevel,
+                    bb_,
+                    n_,
+                    triBb,
+                    level,
+                    maxEqOp<label>()
+                );
+                voxelMeshSearch::fill
+                (
+                    globalRegion,
+                    bb_,
+                    n_,
+                    triBb,
+                    globalRegioni,
+                    maxEqOp<label>()
+                );
+            }
+        }
+        // else: maybe do intersection tests?
+    }
+}
+
+
+void Foam::snappyVoxelMeshDriver::findVoxels
+(
+    const labelList& voxelLevel,
+    const pointField& locationsOutsideMesh,
+    labelList& voxels
+) const
+{
+    voxels.setSize(locationsOutsideMesh.size());
+    voxels = -1;
+    forAll(locationsOutsideMesh, loci)
+    {
+        const point& pt = locationsOutsideMesh[loci];
+        label voxeli = voxelMeshSearch::index(bb_, n_, pt, false);
+
+        if (voxeli == -1 || voxelLevel[voxeli] == labelMax)
+        {
+            WarningInFunction << "Location outside mesh "
+                << pt << " is outside mesh with bounding box "
+                << bb_ << endl;
+        }
+        else
+        {
+            voxels[loci] = voxeli;
+        }
+    }
+}
+
+
+void Foam::snappyVoxelMeshDriver::floodFill
+(
+    const label startVoxeli,
+    const label newLevel,
+    labelList& voxelLevel
+) const
+{
+    DynamicList<labelVector> front;
+    front.append(voxelMeshSearch::index3(n_, startVoxeli));
+
+    DynamicList<labelVector> newFront;
+    while (true)
+    {
+        newFront.clear();
+        for (const auto& voxel : front)
+        {
+            label voxeli = voxelMeshSearch::index(n_, voxel);
+            if (voxelLevel[voxeli] == -1)
+            {
+                voxelLevel[voxeli] = 0;
+                addNeighbours
+                (
+                    voxelLevel,
+                    voxel,
+                    voxeli,
+                    newFront
+                );
+            }
+        }
+
+        if (newFront.empty())
+        {
+            break;
+        }
+        front.transfer(newFront);
+    }
+}
+
+
+void Foam::snappyVoxelMeshDriver::max
+(
+    const labelList& maxLevel,
+    labelList& voxelLevel
+) const
+{
+    // Mark voxels with level
+    const labelVector off(voxelMeshSearch::offset(n_));
+
+    label voxeli = voxelMeshSearch::index(n_, labelVector(0, 0, 0));
+    for (label k = 0; k < n_[2]; k++)
+    {
+        const label start1 = voxeli;
+        for (label j = 0; j < n_[1]; j++)
+        {
+            const label start0 = voxeli;
+            for (label i = 0; i < n_[0]; i++)
+            {
+                voxelLevel[voxeli] = Foam::max
+                (
+                    voxelLevel[voxeli],
+                    maxLevel[voxeli]
+                );
+                voxeli += off[0];
+            }
+            voxeli = start0 + off[1];
+        }
+        voxeli = start1 + off[2];
+    }
+}
+
+
+Foam::labelList Foam::snappyVoxelMeshDriver::count
+(
+    const labelList& voxelLevel
+) const
+{
+
+    label maxLevel = 0;
+    for (const auto level : voxelLevel)
+    {
+        if (level != labelMax)
+        {
+            maxLevel = Foam::max(maxLevel, level);
+        }
+    }
+    labelList count(maxLevel+1, 0);
+
+    const labelVector off(voxelMeshSearch::offset(n_));
+
+    label voxeli = voxelMeshSearch::index(n_, labelVector(0, 0, 0));
+    for (label k = 0; k < n_[2]; k++)
+    {
+        const label start1 = voxeli;
+        for (label j = 0; j < n_[1]; j++)
+        {
+            const label start0 = voxeli;
+            for (label i = 0; i < n_[0]; i++)
+            {
+                label level = voxelLevel[voxeli];
+
+                if (level != -1 && level != labelMax)
+                {
+                    ++count[level];
+                }
+                voxeli += off[0];
+            }
+            voxeli = start0 + off[1];
+        }
+        voxeli = start1 + off[2];
+    }
+
+    return count;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::snappyVoxelMeshDriver::snappyVoxelMeshDriver
+(
+    meshRefinement& meshRefiner,
+    const labelUList& globalToMasterPatch,
+    const labelUList& globalToSlavePatch
+)
+:
+    meshRefiner_(meshRefiner),
+    globalToMasterPatch_(globalToMasterPatch),
+    globalToSlavePatch_(globalToSlavePatch),
+    bb_(meshRefiner_.mesh().bounds())
+{
+    label maxLevel = labelMin;
+
+    // Feature refinement
+    const labelListList& featLevels = meshRefiner_.features().levels();
+    forAll(featLevels, feati)
+    {
+        maxLevel = Foam::max(maxLevel, Foam::max(featLevels[feati]));
+    }
+
+    // Surface refinement
+    const labelList& surfaceLevels = meshRefiner_.surfaces().maxLevel();
+    maxLevel = Foam::max(maxLevel, Foam::max(surfaceLevels));
+
+    // Shell refinement
+    maxLevel = Foam::max(maxLevel, meshRefiner_.shells().maxLevel());
+
+    const scalar level0Len = meshRefiner_.meshCutter().level0EdgeLength();
+
+    Info<< "Cell size estimate :" << nl
+        << "    Maximum : " << level0Len << nl
+        << "    Minimum : " << level0Len/pow(2.0, maxLevel) << nl
+        << endl;
+
+
+    // Define voxel mesh with similar dimensions as mesh
+    const vector meshSpan(bb_.span());
+    n_ = labelVector
+    (
+        round(meshSpan.x()/level0Len),
+        round(meshSpan.y()/level0Len),
+        round(meshSpan.z()/level0Len)
+    );
+    label nTot = n_.x()*n_.y()*n_.z();
+    while (nTot < 1000000)  //1048576)
+    {
+        n_ *= 2;
+        nTot = n_.x()*n_.y()*n_.z();
+    }
+
+    Info<< "Voxellating initial mesh : " << n_ << nl << endl;
+
+    tmp<pointField> tcc(voxelCentres());
+    const pointField& cc = tcc();
+
+    Info<< "Voxel refinement :" << nl
+        << "    Initial                      : (" << nTot << ')' << endl;
+
+    boolList isVoxelInMesh;
+    isInside(cc, isVoxelInMesh);
+
+    if (Pstream::master())
+    {
+        voxelLevel_.setSize(nTot, -1);
+        globalRegion_.setSize(nTot, -1);
+
+        // Remove cells outside initial mesh
+        forAll(isVoxelInMesh, voxeli)
+        {
+            if (!isVoxelInMesh[voxeli])
+            {
+                voxelLevel_[voxeli] = labelMax;
+                globalRegion_[voxeli] = -1;
+            }
+        }
+
+        //if (debug)
+        {
+            Info<< "    After removing outside cells : " << count(voxelLevel_)
+                << endl;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::snappyVoxelMeshDriver::doRefine
+(
+    const refinementParameters& refineParams
+)
+{
+    const scalar level0Len = meshRefiner_.meshCutter().level0EdgeLength();
+
+    tmp<pointField> tcc(voxelCentres());
+    const pointField& cc = tcc();
+
+    boolList isVoxelInMesh;
+    isInside(cc, isVoxelInMesh);
+
+    if (Pstream::master())
+    {
+        // Mark voxels containing (parts of) triangles
+        markSurfaceRefinement(voxelLevel_, globalRegion_);
+
+        //if (debug)
+        {
+            Info<< "    After surface refinement     : " << count(voxelLevel_)
+                << endl;
+        }
+
+
+        // Find outside locations (and their current refinement level)
+        const pointField& outsidePoints = refineParams.locationsOutsideMesh();
+        labelList outsideMeshVoxels;
+        findVoxels
+        (
+            voxelLevel_,
+            outsidePoints,
+            outsideMeshVoxels
+        );
+        labelList outsideOldLevel(outsideMeshVoxels.size(), -1);
+        forAll(outsideMeshVoxels, loci)
+        {
+            label voxeli = outsideMeshVoxels[loci];
+            if (voxeli >= 0)
+            {
+                outsideOldLevel[loci] = voxelLevel_[outsideMeshVoxels[loci]];
+                if (outsideOldLevel[loci] >= 0)
+                {
+                    WarningInFunction << "Location outside mesh "
+                        << outsidePoints[loci]
+                        << " is inside mesh or close to surface" << endl;
+                }
+            }
+        }
+
+
+        // Find inside locations
+        labelList insideMeshVoxels;
+        findVoxels
+        (
+            voxelLevel_,
+            refineParams.locationsInMesh(),
+            insideMeshVoxels
+        );
+
+        forAll(insideMeshVoxels, loci)
+        {
+            label voxeli = insideMeshVoxels[loci];
+            if (voxeli != -1)
+            {
+                if (voxelLevel_[voxeli] != -1)
+                {
+                    WarningInFunction << "Location inside mesh "
+                        << refineParams.locationsInMesh()[loci]
+                        << " is marked as a surface voxel " << voxeli
+                        << " with cell level " << voxelLevel_[voxeli] << endl;
+                }
+                else
+                {
+                    // Flood-fill out from voxel
+                    floodFill(voxeli, 0, voxelLevel_);
+                }
+            }
+        }
+
+        //if (debug)
+        {
+            Info<< "    After keeping inside voxels  : " << count(voxelLevel_)
+                << endl;
+        }
+
+
+        // Re-check the outside locations to see if they have been bled into
+        {
+            forAll(outsideMeshVoxels, loci)
+            {
+                label voxeli = outsideMeshVoxels[loci];
+                if (voxeli >= 0 && voxelLevel_[voxeli] != outsideOldLevel[loci])
+                {
+                    WarningInFunction << "Location outside mesh "
+                        << outsidePoints[loci]
+                        << " is reachable from an inside location" << nl
+                        << "Either your locations are too close to the"
+                        << " geometry or there might be a leak in the"
+                        << " geometry" << endl;
+                }
+            }
+        }
+
+
+        // Shell refinement : find ccs inside higher shells
+        labelList maxLevel;
+        meshRefiner_.shells().findHigherLevel(cc, voxelLevel_, maxLevel);
+
+        // Assign max of maxLevel and voxelLevel
+        max(maxLevel, voxelLevel_);
+
+        // Determine number of levels
+        const labelList levelCounts(count(voxelLevel_));
+
+        //if (debug)
+        {
+            Info<< "    After shell refinement       : " << levelCounts << endl;
+        }
+
+
+        const vector meshSpan(bb_.span());
+        const vector voxel0Size
+        (
+            meshSpan[0]/n_[0],
+            meshSpan[1]/n_[1],
+            meshSpan[2]/n_[2]
+        );
+        label cellCount = 0;
+        forAll(levelCounts, leveli)
+        {
+            const scalar s = level0Len/pow(2.0, leveli);
+            const scalar nCellsPerVoxel
+            (
+                voxel0Size[0]/s
+               *voxel0Size[1]/s
+               *voxel0Size[2]/s
+            );
+            cellCount += levelCounts[leveli]*nCellsPerVoxel;
+        }
+        Info<< "Estimated cell count : " << cellCount << endl;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.H
new file mode 100644
index 0000000000..e241bfcd85
--- /dev/null
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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::snappyVoxelMeshDriver
+
+Description
+    Equivalent of snappyRefineDriver but operating on a voxel mesh.
+
+    Used to estimate cell size count from refinement (currently).
+
+SourceFiles
+    snappyVoxelMeshDriver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef snappyVoxelMeshDriver_H
+#define snappyVoxelMeshDriver_H
+
+#include "labelList.H"
+#include "DynamicList.H"
+#include "labelVector.H"
+#include "boundBox.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declarations
+class refinementParameters;
+class layerParameters;
+class meshRefinement;
+
+/*---------------------------------------------------------------------------*\
+                           Class snappyVoxelMeshDriver Declaration
+\*---------------------------------------------------------------------------*/
+
+class snappyVoxelMeshDriver
+{
+    // Private data
+
+        //- Mesh+surface
+        meshRefinement& meshRefiner_;
+
+        //- From surface region to patch
+        const labelList globalToMasterPatch_;
+
+        //- From surface region to patch
+        const labelList globalToSlavePatch_;
+
+        //- Bounding box
+        boundBox bb_;
+
+        //- Number of voxels
+        labelVector n_;
+
+        //- Per voxel the refinement level
+        labelList voxelLevel_;
+
+        //- Per insersected voxel the originating surface
+        labelList globalRegion_;
+
+
+    // Private Member Functions
+
+        void addNeighbours
+        (
+            const labelList& cellLevel,
+            const labelVector& voxel,
+            const label voxeli,
+            DynamicList<labelVector>& front
+        ) const;
+
+        tmp<pointField> voxelCentres() const;
+
+        void isInside
+        (
+            const pointField& voxelCentres,
+            boolList& isVoxelInMesh
+        ) const;
+
+        void markSurfaceRefinement
+        (
+            labelList& voxelLevel,
+            labelList& globalRegion
+        ) const;
+
+        void findVoxels
+        (
+            const labelList& voxelLevel,
+            const pointField& locationsOutsideMesh,
+            labelList& voxels
+        ) const;
+
+        void floodFill
+        (
+            const label voxeli,
+            const label newLevel,
+            labelList& voxelLevel
+        ) const;
+
+        void max
+        (
+            const labelList& maxLevel,
+            labelList& voxelLevel
+        ) const;
+
+        labelList count(const labelList& voxelLevel) const;
+
+
+       //- No copy construct
+        snappyVoxelMeshDriver(const snappyVoxelMeshDriver&) = delete;
+
+        //- No copy assignment
+        void operator=(const snappyVoxelMeshDriver&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    ClassName("snappyVoxelMeshDriver");
+
+
+    // Constructors
+
+        //- Construct from components
+        snappyVoxelMeshDriver
+        (
+            meshRefinement& meshRefiner,
+            const labelUList& globalToMasterPatch,
+            const labelUList& globalToSlavePatch
+        );
+
+
+    // Member Functions
+
+        void doRefine(const refinementParameters& refineParams);
+
+        //void doLayers(const layerParameters& layerParams);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurfaces/searchableSurfaces/searchableSurfaces.C b/src/meshTools/searchableSurfaces/searchableSurfaces/searchableSurfaces.C
index 144373d15d..91ccdc8d3d 100644
--- a/src/meshTools/searchableSurfaces/searchableSurfaces/searchableSurfaces.C
+++ b/src/meshTools/searchableSurfaces/searchableSurfaces/searchableSurfaces.C
@@ -48,18 +48,51 @@ bool Foam::searchableSurfaces::connected
     const pointIndexHit& hit
 )
 {
-    const triFace& localFace = s.localFaces()[hit.index()];
     const edge& e = s.edges()[edgeI];
+    const labelList& mp = s.meshPoints();
+    const edge meshE(mp[e[0]], mp[e[1]]);
 
-    forAll(localFace, i)
+    const triFace& f = s[hit.index()];
+
+    forAll(f, i)
     {
-        if (e.otherVertex(localFace[i]) != -1)
+        if (meshE.otherVertex(f[i]) != -1)
         {
             return true;
         }
     }
 
-    return false;
+    // Account for triangle intersection routine going wrong for
+    // lines in same plane as triangle. Tbd!
+
+    vector eVec(meshE.vec(s.points()));
+    scalar magEVec(mag(eVec));
+    if (magEVec > ROOTVSMALL)
+    {
+        vector n(f.areaNormal(s.points()));
+        scalar magArea(mag(n));
+        if (magArea > ROOTVSMALL)
+        {
+            n /= magArea;
+            if (mag(n&(eVec/magEVec)) < SMALL)
+            {
+                // Bad intersection
+                return true;
+            }
+            else
+            {
+                return false;
+            }
+        }
+        else
+        {
+            return true;
+        }
+    }
+    else
+    {
+        return true;
+    }
 }
 
 
@@ -619,20 +652,20 @@ bool Foam::searchableSurfaces::checkIntersection
             {
                 List<pointIndexHit> hits;
 
-                if (i == j)
-                {
-                    // Slightly shorten the edges to prevent finding lots of
-                    // intersections. The fast triangle intersection routine
-                    // has problems with rays passing through a point of the
-                    // triangle.
-                    vectorField d
-                    (
-                        max(tolerance, 10*s0.tolerance())
-                       *(end-start)
-                    );
-                    start += d;
-                    end -= d;
-                }
+                //if (i == j)
+                //{
+                //    // Slightly shorten the edges to prevent finding lots of
+                //    // intersections. The fast triangle intersection routine
+                //    // has problems with rays passing through a point of the
+                //    // triangle. Now done in 'connected' routine. Tbd.
+                //    vectorField d
+                //    (
+                //        max(tolerance, 10*s0.tolerance())
+                //       *(end-start)
+                //    );
+                //    start += d;
+                //    end -= d;
+                //}
 
                 operator[](j).findLineAny(start, end, hits);
 
@@ -654,6 +687,9 @@ bool Foam::searchableSurfaces::checkIntersection
                     }
                 }
 
+                // tdb. What about distributedTriSurfaceMesh?
+                //reduce(nHits, sumOp<label>());
+
                 if (nHits > 0)
                 {
                     if (report)
diff --git a/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.C b/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.C
index d1e23e278e..2b6e23baf8 100644
--- a/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.C
+++ b/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.C
@@ -42,6 +42,43 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+Foam::labelVector Foam::voxelMeshSearch::offset
+(
+    const labelVector& nDivs
+)
+{
+    return labelVector(1, nDivs.x(), nDivs.x()*nDivs.y());
+}
+
+
+Foam::label Foam::voxelMeshSearch::index
+(
+    const labelVector& nDivs,
+    const labelVector& voxel
+)
+{
+    return voxel.x()+voxel.y()*nDivs.x()+voxel.z()*nDivs.x()*nDivs.y();
+}
+
+
+Foam::labelVector Foam::voxelMeshSearch::index3
+(
+    const labelVector& nDivs,
+    label voxeli
+)
+{
+    const label nxy = nDivs.x()*nDivs.y();
+
+    labelVector voxel;
+    voxel.z() = voxeli/nxy;
+    voxeli = voxeli % nxy;
+    voxel.y() = voxeli/nDivs.x();
+    voxel.x() = voxeli%nDivs.x();
+
+    return voxel;
+}
+
+
 Foam::labelVector Foam::voxelMeshSearch::index3
 (
     const boundBox& bb,
@@ -98,7 +135,7 @@ Foam::label Foam::voxelMeshSearch::index
         return -1;
     }
 
-    return v[0] + g[0]*v[1] + g[1]*g.y()*v[2];
+    return index(g, v);
 }
 
 
diff --git a/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.H b/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.H
index d7ce1e543f..691c3c73cc 100644
--- a/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.H
+++ b/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.H
@@ -122,6 +122,26 @@ public:
 
         //Voxel helper functions
 
+            //- Voxel indices to combined voxel index
+            static label index
+            (
+                const labelVector& nDivs,
+                const labelVector& voxel
+            );
+
+            //- Change in combined voxel index for change in components
+            static labelVector offset
+            (
+                const labelVector& nDivs
+            );
+
+            //- Combined voxel index to individual indices
+            static labelVector index3
+            (
+                const labelVector& nDivs,
+                const label voxeli
+            );
+
             //- Coordinate to voxel indices
             static labelVector index3
             (
@@ -130,6 +150,14 @@ public:
                 const point& p
             );
 
+            //- Voxel index to voxel centre
+            static point centre
+            (
+                const boundBox& bb,
+                const labelVector& nDivs,
+                const labelVector& voxel
+            );
+
             //- Coordinate to combined voxel index. If clip makes sure
             //  components are all inside. If not clip returns -1 if outside bb.
             static label index
@@ -140,14 +168,6 @@ public:
                 const bool clip
             );
 
-            //- Voxel index to voxel centre
-            static point centre
-            (
-                const boundBox& bb,
-                const labelVector& nDivs,
-                const labelVector& voxel
-            );
-
             //- Fill voxels indicated by bounding box
             template<class Container, class Type>
             static void fill
@@ -159,6 +179,18 @@ public:
                 const Type val
             );
 
+            //- Fill voxels indicated by bounding box
+            template<class Container, class Type, class CombineOp>
+            static void fill
+            (
+                Container& elems,
+                const boundBox& bb,
+                const labelVector& nDivs,
+                const boundBox& subBb,
+                const Type val,
+                const CombineOp& cop = eqOp<Type>()
+            );
+
             //- Check if any voxel inside bounding box is set to val or
             //  not set to val (isNot = true)
             template<class Container, class Type>
diff --git a/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearchTemplates.C b/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearchTemplates.C
index e748d23281..adcac1dad3 100644
--- a/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearchTemplates.C
+++ b/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearchTemplates.C
@@ -52,15 +52,67 @@ void Foam::voxelMeshSearch::fill
         minIds[cmpt] = max(minIds[cmpt], 0);
     }
 
-    for (label i = minIds[0]; i <= maxIds[0]; i++)
+    const labelVector off(offset(nDivs));
+    label voxeli = index(nDivs, minIds);
+    for (label k = minIds[2]; k <= maxIds[2]; k++)
     {
+        const label start1 = voxeli;
         for (label j = minIds[1]; j <= maxIds[1]; j++)
         {
-            for (label k = minIds[2]; k <= maxIds[2]; k++)
+            const label start0 = voxeli;
+            for (label i = minIds[0]; i <= maxIds[0]; i++)
             {
-                elems[i+j*nDivs.x()+k*nDivs.x()*nDivs.y()] = val;
+                elems[voxeli] = val;
+                voxeli += off[0];
             }
+            voxeli = start0 + off[1];
         }
+        voxeli = start1 + off[2];
+    }
+}
+
+
+template<class Container, class Type, class CombineOp>
+void Foam::voxelMeshSearch::fill
+(
+    Container& elems,
+    const boundBox& bb,
+    const labelVector& nDivs,
+    const boundBox& subBb,
+    const Type val,
+    const CombineOp& cop
+)
+{
+    labelVector minIds(index3(bb, nDivs, subBb.min()));
+    labelVector maxIds(index3(bb, nDivs, subBb.max()));
+
+    for (direction cmpt = 0; cmpt < 3; cmpt++)
+    {
+        if (maxIds[cmpt] < 0 || minIds[cmpt] >= nDivs[cmpt])
+        {
+            return;
+        }
+        // Clip
+        maxIds[cmpt] = min(maxIds[cmpt], nDivs[cmpt]-1);
+        minIds[cmpt] = max(minIds[cmpt], 0);
+    }
+
+    const labelVector off(offset(nDivs));
+    label voxeli = index(nDivs, minIds);
+    for (label k = minIds[2]; k <= maxIds[2]; k++)
+    {
+        const label start1 = voxeli;
+        for (label j = minIds[1]; j <= maxIds[1]; j++)
+        {
+            const label start0 = voxeli;
+            for (label i = minIds[0]; i <= maxIds[0]; i++)
+            {
+                cop(elems[voxeli], val);
+                voxeli += off[0];
+            }
+            voxeli = start0 + off[1];
+        }
+        voxeli = start1 + off[2];
     }
 }
 
@@ -100,19 +152,26 @@ bool Foam::voxelMeshSearch::overlaps
     }
 
 
-    for (label i = minIds[0]; i <= maxIds[0]; i++)
+    const labelVector off(offset(nDivs));
+    label voxeli = index(nDivs, minIds);
+    for (label k = minIds[2]; k <= maxIds[2]; k++)
     {
+        const label start1 = voxeli;
         for (label j = minIds[1]; j <= maxIds[1]; j++)
         {
-            for (label k = minIds[2]; k <= maxIds[2]; k++)
+            const label start0 = voxeli;
+            for (label i = minIds[0]; i <= maxIds[0]; i++)
             {
-                const Type elemVal = elems[i+j*nDivs.x()+k*nDivs.x()*nDivs.y()];
+                const Type elemVal = elems[voxeli];
                 if (isNot != (elemVal == val))
                 {
                     return true;
                 }
+                voxeli += off[0];
             }
+            voxeli = start0 + off[1];
         }
+        voxeli = start1 + off[2];
     }
     return false;
 }
@@ -136,19 +195,26 @@ void Foam::voxelMeshSearch::write
             << exit(FatalError);
     }
 
-    for (label i = 0; i < nDivs[0]; i++)
+    const labelVector off(offset(nDivs));
+    label voxeli = index(nDivs, labelVector(0, 0, 0));
+    for (label k = 0; k < nDivs[2]; k++)
     {
+        const label start1 = voxeli;
         for (label j = 0; j < nDivs[1]; j++)
         {
-            for (label k = 0; k < nDivs[2]; k++)
+            const label start0 = voxeli;
+            for (label i = 0; i < nDivs[0]; i++)
             {
-                const Type elemVal = elems[i+j*nDivs.x()+k*nDivs.x()*nDivs.y()];
+                const Type& elemVal = elems[voxeli];
                 if (isNot != (elemVal == val))
                 {
                     os.write(centre(bb, nDivs, labelVector(i, j, k)));
                 }
+                voxeli += off[0];
             }
+            voxeli = start0 + off[1];
         }
+        voxeli = start1 + off[2];
     }
 }