diff --git a/applications/utilities/mesh/advanced/runDynamicMesh/Make/files b/applications/utilities/mesh/advanced/runDynamicMesh/Make/files
new file mode 100644
index 000000000..55ca1fb38
--- /dev/null
+++ b/applications/utilities/mesh/advanced/runDynamicMesh/Make/files
@@ -0,0 +1,3 @@
+runDynamicMesh.C
+
+EXE = $(FOAM_APPBIN)/runDynamicMesh
diff --git a/applications/utilities/mesh/advanced/runDynamicMesh/Make/options b/applications/utilities/mesh/advanced/runDynamicMesh/Make/options
new file mode 100644
index 000000000..9bb9fdd78
--- /dev/null
+++ b/applications/utilities/mesh/advanced/runDynamicMesh/Make/options
@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -ldynamicFvMesh \
+    -ldynamicMesh \
+    -lmeshTools \
+    -lfiniteVolume \
+    -L$(MESQUITE_LIB_DIR) -lmesquite
diff --git a/applications/utilities/mesh/advanced/runDynamicMesh/createFields.H b/applications/utilities/mesh/advanced/runDynamicMesh/createFields.H
new file mode 100644
index 000000000..979389f96
--- /dev/null
+++ b/applications/utilities/mesh/advanced/runDynamicMesh/createFields.H
@@ -0,0 +1,13 @@
+    Info<< "Reading field alpha\n" << endl;
+    volScalarField alpha
+    (
+        IOobject
+        (
+            "alpha",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
diff --git a/applications/utilities/mesh/advanced/runDynamicMesh/runDynamicMesh.C b/applications/utilities/mesh/advanced/runDynamicMesh/runDynamicMesh.C
new file mode 100644
index 000000000..b4d50f036
--- /dev/null
+++ b/applications/utilities/mesh/advanced/runDynamicMesh/runDynamicMesh.C
@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    runDynamicMesh
+
+Description
+    Creates dynamicFvMesh and updates it in a time loop. Reads indicator field
+    alpha used for dynamic mesh refinement to trigger refinement in certain
+    areas.
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "simpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createDynamicFvMesh.H"
+
+    simpleControl simple(mesh);
+
+#   include "createFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        bool meshChanged = mesh.update();
+
+        // Write only if mesh has changed
+        if (meshChanged)
+        {
+            runTime.write();
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/Make/files b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/Make/files
index c8d0cc571..abcfb8261 100644
--- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/Make/files
+++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/Make/files
@@ -4,6 +4,7 @@ linearNormal/linearNormal.C
 linearRadial/linearRadial.C
 sigmaRadial/sigmaRadial.C
 wedge/wedge.C
+gradedNormal/gradedNormal.C
 
 LIB = $(FOAM_LIBBIN)/libextrudeModel
 
diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/Make/options b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/Make/options
index 3c57befe4..afc486fed 100644
--- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/Make/options
+++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/Make/options
@@ -1,8 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude
 
 EXE_LIBS = \
     -lmeshTools \
-    -ldynamicMesh
-
+    -ldynamicMesh \
+    -lODE
diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/gradedNormal/gradedNormal.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/gradedNormal/gradedNormal.C
new file mode 100644
index 000000000..37bb95212
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/gradedNormal/gradedNormal.C
@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "gradedNormal.H"
+#include "addToRunTimeSelectionTable.H"
+#include "BisectionRoot.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace extrudeModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(gradedNormal, 0);
+
+addToRunTimeSelectionTable(extrudeModel, gradedNormal, dictionary);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+gradedNormal::gradedNormal(const dictionary& dict)
+:
+    extrudeModel(typeName, dict),
+    thickness_(readScalar(coeffDict_.lookup("thickness"))),
+    delta0_(readScalar(coeffDict_.lookup("initialCellLength"))),
+    expansionFactor_(1.0)
+{
+    // Sanity checks
+    if (thickness_ <= SMALL)
+    {
+        FatalErrorIn("gradedNormal(const dictionary&)")
+            << "thickness should be positive: " << thickness_
+            << exit(FatalError);
+    }
+
+    if (delta0_ <= SMALL)
+    {
+        FatalErrorIn("gradedNormal(const dictionary&)")
+            << "initialCellLength should be positive: " << delta0_
+            << exit(FatalError);
+    }
+
+    const scalar maxExpFactor =
+        coeffDict_.lookupOrDefault<scalar>("maxExpansionFactor", 3.0);
+
+    if (maxExpFactor <= SMALL)
+    {
+        FatalErrorIn("gradedNormal(const dictionary&)")
+            << "maxExpansionFactor should be positive: " << maxExpFactor
+            << exit(FatalError);
+    }
+
+    const scalar bisectionTol =
+        coeffDict_.lookupOrDefault<scalar>("bisectionTol", 1e-5);
+
+    if (bisectionTol <= SMALL)
+    {
+        FatalErrorIn("gradedNormal(const dictionary&)")
+            << "bisectionTolerance should be positive: " << bisectionTol
+            << exit(FatalError);
+    }
+
+    // Create expansion factor equation represented as a function object
+    expansionFactorEqn eqn(*this);
+
+    // Calculate the expansionFactor using the bisection algorithm with the
+    // default tolerance of 1e-5
+    BisectionRoot<expansionFactorEqn> rootFinder
+    (
+        eqn,
+        bisectionTol
+    );
+
+    // Search for the root in [0, 3], where upper bound 3 is default
+    expansionFactor_ = rootFinder.root
+    (
+        0.0,
+        maxExpFactor
+    );
+
+    // Report the result
+    Info<< "Calculated expansion factor: " << expansionFactor_ << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+gradedNormal::~gradedNormal()
+{}
+
+
+// * * * * * * * * * * * * * * * * Operators * * * * * * * * * * * * * * * * //
+
+point gradedNormal::operator()
+(
+    const point& surfacePoint,
+    const vector& surfaceNormal,
+    const label layer
+) const
+{
+    scalar d = 0.0;
+
+    for (label i = 0; i < layer; ++i)
+    {
+        d += delta0_*pow(expansionFactor_, i);
+    }
+
+    return surfacePoint + d*surfaceNormal;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace extrudeModels
+} // End namespace Foam
+
+// ************************************************************************* //
+
diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/gradedNormal/gradedNormal.H b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/gradedNormal/gradedNormal.H
new file mode 100644
index 000000000..bd2127bf0
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/gradedNormal/gradedNormal.H
@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::extrudeModels::gradedNormal
+
+Description
+    Extrudes by transforming points normal to the surface. Input parameters are:
+    1. Extrusion thickness (in meters),
+    2. Initial delta (in meters),
+    3. Number of layers.
+
+    Expansion factor is calculated numerically using bisection algorithm.
+
+Author
+    Vuko Vukcevic. FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef gradedNormal_H
+#define gradedNormal_H
+
+#include "point.H"
+#include "extrudeModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace extrudeModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class gradedNormal Declaration
+\*---------------------------------------------------------------------------*/
+
+class gradedNormal
+:
+    public extrudeModel
+{
+    // Private data
+
+        //- Layer thickness
+        scalar thickness_;
+
+        //- Initial delta (cell size at the extruded patch)
+        scalar delta0_;
+
+        //- Expansion factor
+        scalar expansionFactor_;
+
+
+    // Expansion factor equation (functor class used by BisectionRoot)
+    class expansionFactorEqn
+    {
+        // Private data
+
+            //- Const reference to underlying gradedNormal model
+            const gradedNormal& gnm_;
+
+
+    public:
+
+        //- Construct given gradedNormal model
+        explicit expansionFactorEqn(const gradedNormal& gnm)
+        :
+            gnm_(gnm)
+        {}
+
+
+        //- Function call operator
+        scalar operator()(const scalar& x) const
+        {
+            scalar result = gnm_.thickness();
+
+            for (label i = 0; i <= gnm_.nLayers(); ++i)
+            {
+                result -= gnm_.delta0()*pow(x, i);
+            }
+
+            return result;
+        }
+    };
+
+
+public:
+
+    //- Runtime type information
+    TypeName("gradedNormal");
+
+    // Constructors
+
+        //- Construct from components
+        gradedNormal(const dictionary& dict);
+
+
+    //- Destructor
+    ~gradedNormal();
+
+
+    // Access functions
+
+        //- Return const reference to thickness
+        const scalar& thickness() const
+        {
+            return thickness_;
+        }
+
+
+        //- Return const reference to initial delta
+        const scalar& delta0() const
+        {
+            return delta0_;
+        }
+
+
+    // Member Operators
+
+        point operator()
+        (
+            const point& surfacePoint,
+            const vector& surfaceNormal,
+            const label layer
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace extrudeModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
index 68c984865..015761d1b 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
@@ -103,6 +103,37 @@ bool domainDecomposition::writeDecomposition()
     label maxProcFaces = 0;
 
 
+    // Note: get cellLevel and pointLevel. Avoid checking whether they exist or
+    // not by hand. If they don't exist, simpy assume that the level is 0
+    const labelIOList globalCellLevel
+    (
+        IOobject
+        (
+            "cellLevel",
+            this->facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        labelList(nCells(), 0)
+    );
+
+    const labelIOList globalPointLevel
+    (
+        IOobject
+        (
+            "pointLevel",
+            this->facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        labelList(nPoints(), 0)
+    );
+
+
     // Write out the meshes
     for (label procI = 0; procI < nProcs_; procI++)
     {
@@ -615,6 +646,39 @@ bool domainDecomposition::writeDecomposition()
             procBoundaryAddressing_[procI]
         );
         boundaryProcAddressing.write();
+
+        // Create and write cellLevel and pointLevel information
+        const unallocLabelList& cellMap = cellProcAddressing;
+        labelIOList procCellLevel
+        (
+            IOobject
+            (
+                "cellLevel",
+                procMesh.facesInstance(),
+                procMesh.meshSubDir,
+                procMesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            labelList(globalCellLevel, cellMap)
+        );
+        procCellLevel.write();
+
+        const unallocLabelList& pointMap = pointProcAddressing;
+        labelIOList procPointLevel
+        (
+            IOobject
+            (
+                "pointLevel",
+                procMesh.facesInstance(),
+                procMesh.meshSubDir,
+                procMesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            labelList(globalPointLevel, pointMap)
+        );
+        procPointLevel.write();
     }
 
     Info<< nl
diff --git a/src/dynamicMesh/dynamicFvMesh/Make/files b/src/dynamicMesh/dynamicFvMesh/Make/files
index 46aee8b5e..1e0932d17 100644
--- a/src/dynamicMesh/dynamicFvMesh/Make/files
+++ b/src/dynamicMesh/dynamicFvMesh/Make/files
@@ -10,6 +10,7 @@ dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C
 subsetMotionSolverFvMesh/subsetMotionSolverFvMesh.C
 dynamicInkJetFvMesh/dynamicInkJetFvMesh.C
 dynamicRefineFvMesh/dynamicRefineFvMesh.C
+dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.C
 
 mixerGgiFvMesh/mixerGgiFvMesh.C
 turboFvMesh/turboFvMesh.C
diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicMeshDict b/src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicMeshDict
new file mode 100644
index 000000000..b25b10d6e
--- /dev/null
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicMeshDict
@@ -0,0 +1,77 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     4.0                                |
+|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//dynamicFvMeshLib   "libtopoChangerFvMesh.so";
+dynamicFvMesh      dynamicRefinePolyFvMesh;
+//staticFvMesh;
+
+mixerFvMeshCoeffs
+{
+    coordinateSystem
+    {
+        type            cylindrical;
+        origin          (0 0 0);
+        axis            (0 0 1);
+        direction       (1 0 0);
+    }
+
+    rpm             10;
+
+    slider
+    {
+        inside      insideSlider;
+        outside     outsideSlider;
+    }
+}
+
+// Refinement
+dynamicRefineFvMeshCoeffs
+{
+    // Refine every refineInterval timesteps
+    refineInterval 3;
+
+    // Maximum refinement level (starts from 0)
+    maxRefinement 2;
+
+    // Maximum cell limit (approximate)
+    maxCells 1000000;
+
+    // volScalarField to base refinement on
+    field gamma;
+
+    // Which cells to un/refine: based on point values (simple averaging).
+    // - refine pointCells of point value inbetween minLevel..maxLevel
+    // - unrefine pointCells that are within nBufferLayers of points marked
+    //   for refinement.
+    minLevel 0.01;
+    maxLevel 0.99;
+    nBufferLayers 1;
+
+    // Newly introduced patch points optionally get projected onto a surface
+    //projectSurfaces ("fixedWalls4.stl");
+    //projectPatches (fixedWalls);
+    // Maximum project distance
+    //projectDistance 1;
+
+    // Fluxes to adapt. For newly created faces or split faces the flux
+    // gets estimated from an interpolated volVectorField ('velocity')
+    // First is name of the flux to adapt, second is velocity that will
+    // be interpolated and inner-producted with the face area vector.
+    correctFluxes ((phi U));
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.C b/src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.C
new file mode 100644
index 000000000..b7d966e01
--- /dev/null
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.C
@@ -0,0 +1,1035 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "dynamicRefinePolyFvMesh.H"
+#include "addToRunTimeSelectionTable.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "fvCFD.H"
+#include "syncTools.H"
+#include "pointFields.H"
+#include "directTopoChange.H"
+#include "cellSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(dynamicRefinePolyFvMesh, 0);
+
+addToRunTimeSelectionTable(dynamicFvMesh, dynamicRefinePolyFvMesh, IOobject);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+label dynamicRefinePolyFvMesh::count
+(
+    const PackedBoolList& l,
+    const unsigned int val
+)
+{
+    label n = 0;
+    forAll(l, i)
+    {
+        if (l.get(i) == val)
+        {
+            n++;
+        }
+
+        // debug also serves to get-around Clang compiler trying to optimise
+        // out this forAll loop under O3 optimisation
+
+        if (debug)
+        {
+            Info<< "n=" << n << endl;
+        }
+    }
+
+    return n;
+}
+
+
+void dynamicRefinePolyFvMesh::readDict()
+{
+    dictionary refineDict
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "dynamicMeshDict",
+                time().constant(),
+                *this,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).subDict(typeName + "Coeffs")
+    );
+
+    correctFluxes_ = List<Pair<word> >(refineDict.lookup("correctFluxes"));
+
+    dumpLevel_ = Switch(refineDict.lookup("dumpLevel"));
+}
+
+
+// Refines cells, maps fields and recalculates (an approximate) flux
+autoPtr<mapPolyMesh> dynamicRefinePolyFvMesh::refine
+(
+    const labelList& cellsToRefine
+)
+{
+    // Mesh changing engine.
+    directTopoChange meshMod(*this);
+
+    // Play refinement commands into mesh changer.
+    meshCutter_.setRefinement(cellsToRefine, meshMod);
+
+    // Create mesh (with inflation), return map from old to new mesh.
+    //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
+
+    Info<< "Refined from "
+        << returnReduce(map().nOldCells(), sumOp<label>())
+        << " to " << globalData().nTotalCells() << " cells." << endl;
+
+    if (debug)
+    {
+        // Check map.
+        for (label faceI = 0; faceI < nInternalFaces(); faceI++)
+        {
+            label oldFaceI = map().faceMap()[faceI];
+
+            if (oldFaceI >= nInternalFaces())
+            {
+                FatalErrorIn("dynamicRefinePolyFvMesh::refine(const labelList&)")
+                    << "New internal face:" << faceI
+                    << " fc:" << faceCentres()[faceI]
+                    << " originates from boundary oldFace:" << oldFaceI
+                    << abort(FatalError);
+            }
+        }
+    }
+
+
+    // Update fields
+    updateMesh(map);
+
+    // Move mesh
+    /*
+    pointField newPoints;
+    if (map().hasMotionPoints())
+    {
+        newPoints = map().preMotionPoints();
+    }
+    else
+    {
+        newPoints = points();
+    }
+    movePoints(newPoints);
+    */
+
+    // Correct the flux for modified/added faces. All the faces which only
+    // have been renumbered will already have been handled by the mapping.
+    {
+        const labelList& faceMap = map().faceMap();
+        const labelList& reverseFaceMap = map().reverseFaceMap();
+
+        // Storage for any master faces. These will be the original faces on the
+        // coarse cell that gets split into four for hex cell and quad
+        // face. Rather, the master face gets modified and three faces get added
+        // from the master.
+        labelHashSet masterFaces(4*cellsToRefine.size());
+
+        forAll(faceMap, faceI)
+        {
+            label oldFaceI = faceMap[faceI];
+
+            if (oldFaceI >= 0)
+            {
+                label masterFaceI = reverseFaceMap[oldFaceI];
+
+                if (masterFaceI < 0)
+                {
+                    FatalErrorIn
+                    (
+                        "dynamicRefinePolyFvMesh::refine(const labelList&)"
+                    )   << "Problem: should not have removed faces"
+                        << " when refining."
+                        << nl << "face:" << faceI << abort(FatalError);
+                }
+                else if (masterFaceI != faceI)
+                {
+                    masterFaces.insert(masterFaceI);
+                }
+            }
+        }
+        if (debug)
+        {
+            Info<< "Found " << returnReduce(masterFaces.size(), sumOp<label>())
+                << " split faces " << endl;
+        }
+
+        forAll(correctFluxes_, i)
+        {
+            if (debug)
+            {
+                Info<< "Mapping flux " << correctFluxes_[i][0]
+                    << " using interpolated flux " << correctFluxes_[i][1]
+                    << endl;
+            }
+            surfaceScalarField& phi = const_cast<surfaceScalarField&>
+            (
+                lookupObject<surfaceScalarField>(correctFluxes_[i][0])
+            );
+            surfaceScalarField phiU =
+                fvc::interpolate
+                (
+                    lookupObject<volVectorField>(correctFluxes_[i][1])
+                )
+              & Sf();
+
+            // Recalculate new internal faces.
+            for (label faceI = 0; faceI < nInternalFaces(); faceI++)
+            {
+                label oldFaceI = faceMap[faceI];
+
+                if (oldFaceI == -1)
+                {
+                    // Inflated/appended
+                    phi[faceI] = phiU[faceI];
+                }
+                else if (reverseFaceMap[oldFaceI] != faceI)
+                {
+                    // face-from-masterface
+                    phi[faceI] = phiU[faceI];
+                }
+            }
+
+            // Recalculate new boundary faces.
+            forAll(phi.boundaryField(), patchI)
+            {
+                fvsPatchScalarField& patchPhi = phi.boundaryField()[patchI];
+                const fvsPatchScalarField& patchPhiU =
+                    phiU.boundaryField()[patchI];
+
+                label faceI = patchPhi.patch().patch().start();
+
+                forAll(patchPhi, i)
+                {
+                    label oldFaceI = faceMap[faceI];
+
+                    if (oldFaceI == -1)
+                    {
+                        // Inflated/appended
+                        patchPhi[i] = patchPhiU[i];
+                    }
+                    else if (reverseFaceMap[oldFaceI] != faceI)
+                    {
+                        // face-from-masterface
+                        patchPhi[i] = patchPhiU[i];
+                    }
+
+                    faceI++;
+                }
+            }
+
+            // Update master faces
+            forAllConstIter(labelHashSet, masterFaces, iter)
+            {
+                label faceI = iter.key();
+
+                if (isInternalFace(faceI))
+                {
+                    phi[faceI] = phiU[faceI];
+                }
+                else
+                {
+                    label patchI = boundaryMesh().whichPatch(faceI);
+                    label i = faceI - boundaryMesh()[patchI].start();
+
+                    const fvsPatchScalarField& patchPhiU =
+                        phiU.boundaryField()[patchI];
+
+                    fvsPatchScalarField& patchPhi =
+                        phi.boundaryField()[patchI];
+
+                    patchPhi[i] = patchPhiU[i];
+                }
+            }
+        }
+    }
+
+
+
+    // Update numbering of cells/vertices.
+    meshCutter_.updateMesh(map);
+
+    // Debug: Check refinement levels (across faces only)
+    meshCutter_.checkRefinementLevels(-1, labelList(0));
+
+    return map;
+}
+
+
+// Combines previously split cells, maps fields and recalculates
+// (an approximate) flux
+autoPtr<mapPolyMesh> dynamicRefinePolyFvMesh::unrefine
+(
+    const labelList& splitPoints
+)
+{
+    directTopoChange meshMod(*this);
+
+    // Play refinement commands into mesh changer.
+    meshCutter_.setUnrefinement(splitPoints, meshMod);
+
+
+    // Save information on faces that will be combined
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Find the faceMidPoints on cells to be combined.
+    // for each face resulting of split of face into four store the
+    // midpoint
+    Map<label> faceToSplitPoint(3*splitPoints.size());
+
+    {
+        forAll(splitPoints, i)
+        {
+            label pointI = splitPoints[i];
+
+            const labelList& pEdges = pointEdges()[pointI];
+
+            forAll(pEdges, j)
+            {
+                label otherPointI = edges()[pEdges[j]].otherVertex(pointI);
+
+                const labelList& pFaces = pointFaces()[otherPointI];
+
+                forAll(pFaces, pFaceI)
+                {
+                    faceToSplitPoint.insert(pFaces[pFaceI], otherPointI);
+                }
+            }
+        }
+    }
+
+
+    // Change mesh and generate map.
+    //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
+
+    Info<< "Unrefined from "
+        << returnReduce(map().nOldCells(), sumOp<label>())
+        << " to " << globalData().nTotalCells() << " cells."
+        << endl;
+
+    // Update fields
+    updateMesh(map);
+
+
+    // Move mesh
+    /*
+    pointField newPoints;
+    if (map().hasMotionPoints())
+    {
+        newPoints = map().preMotionPoints();
+    }
+    else
+    {
+        newPoints = points();
+    }
+    movePoints(newPoints);
+    */
+
+    // Correct the flux for modified faces.
+    {
+        const labelList& reversePointMap = map().reversePointMap();
+        const labelList& reverseFaceMap = map().reverseFaceMap();
+
+        forAll(correctFluxes_, i)
+        {
+            if (debug)
+            {
+                Info<< "Mapping flux " << correctFluxes_[i][0]
+                    << " using interpolated flux " << correctFluxes_[i][1]
+                    << endl;
+            }
+            surfaceScalarField& phi = const_cast<surfaceScalarField&>
+            (
+                lookupObject<surfaceScalarField>(correctFluxes_[i][0])
+            );
+            surfaceScalarField phiU =
+                fvc::interpolate
+                (
+                    lookupObject<volVectorField>(correctFluxes_[i][1])
+                )
+              & Sf();
+
+            forAllConstIter(Map<label>, faceToSplitPoint, iter)
+            {
+                label oldFaceI = iter.key();
+                label oldPointI = iter();
+
+                if (reversePointMap[oldPointI] < 0)
+                {
+                    // midpoint was removed. See if face still exists.
+                    label faceI = reverseFaceMap[oldFaceI];
+
+                    if (faceI >= 0)
+                    {
+                        if (isInternalFace(faceI))
+                        {
+                            phi[faceI] = phiU[faceI];
+                        }
+                        else
+                        {
+                            label patchI = boundaryMesh().whichPatch(faceI);
+                            label i = faceI - boundaryMesh()[patchI].start();
+
+                            const fvsPatchScalarField& patchPhiU =
+                                phiU.boundaryField()[patchI];
+
+                            fvsPatchScalarField& patchPhi =
+                                phi.boundaryField()[patchI];
+
+                            patchPhi[i] = patchPhiU[i];
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
+    // Update numbering of cells/vertices.
+    meshCutter_.updateMesh(map);
+
+    // Debug: Check refinement levels (across faces only)
+    meshCutter_.checkRefinementLevels(-1, labelList(0));
+
+    return map;
+}
+
+
+// Get max of connected point
+scalarField dynamicRefinePolyFvMesh::maxPointField(const scalarField& pFld) const
+{
+    scalarField vFld(nCells(), -GREAT);
+
+    forAll(pointCells(), pointI)
+    {
+        const labelList& pCells = pointCells()[pointI];
+
+        forAll(pCells, i)
+        {
+            vFld[pCells[i]] = max(vFld[pCells[i]], pFld[pointI]);
+        }
+    }
+    return vFld;
+}
+
+
+// Get min of connected cell
+scalarField dynamicRefinePolyFvMesh::minCellField(const volScalarField& vFld) const
+{
+    scalarField pFld(nPoints(), GREAT);
+
+    forAll(pointCells(), pointI)
+    {
+        const labelList& pCells = pointCells()[pointI];
+
+        forAll(pCells, i)
+        {
+            pFld[pointI] = min(pFld[pointI], vFld[pCells[i]]);
+        }
+    }
+    return pFld;
+}
+
+
+// Simple (non-parallel) interpolation by averaging.
+scalarField dynamicRefinePolyFvMesh::cellToPoint(const scalarField& vFld) const
+{
+    scalarField pFld(nPoints());
+
+    forAll(pointCells(), pointI)
+    {
+        const labelList& pCells = pointCells()[pointI];
+
+        scalar sum = 0.0;
+        forAll(pCells, i)
+        {
+            sum += vFld[pCells[i]];
+        }
+        pFld[pointI] = sum/pCells.size();
+    }
+    return pFld;
+}
+
+
+// Calculate error. Is < 0 or distance from inbetween levels
+scalarField dynamicRefinePolyFvMesh::error
+(
+    const scalarField& fld,
+    const scalar minLevel,
+    const scalar maxLevel
+) const
+{
+    const scalar halfLevel = 0.5*(minLevel + maxLevel);
+
+    scalarField c(fld.size(), -1);
+
+    forAll(fld, i)
+    {
+        if (fld[i] >= minLevel && fld[i] < maxLevel)
+        {
+            c[i] = mag(fld[i] - halfLevel);
+        }
+    }
+    return c;
+}
+
+
+void dynamicRefinePolyFvMesh::selectRefineCandidates
+(
+    const scalar lowerRefineLevel,
+    const scalar upperRefineLevel,
+    const scalarField& vFld,
+    PackedBoolList& candidateCell
+) const
+{
+    // Get error per cell. Is -1 (not to be refined) to >0 (to be refined,
+    // higher more desirable to be refined).
+    scalarField cellError
+    (
+        maxPointField
+        (
+            error
+            (
+                cellToPoint(vFld),
+                lowerRefineLevel,
+                upperRefineLevel
+            )
+        )
+    );
+
+    // Mark cells that are candidates for refinement.
+    forAll(cellError, cellI)
+    {
+        if (cellError[cellI] > 0)
+        {
+            candidateCell.set(cellI, 1);
+        }
+    }
+}
+
+
+labelList dynamicRefinePolyFvMesh::selectRefineCells
+(
+    const label maxCells,
+    const label maxRefinement,
+    const PackedBoolList& candidateCell
+) const
+{
+    // Note: Assuming predominantly hex mesh, i.e. every refined cell causes 7
+    // extra cells. This might lead to slight over shoot of maximum number of
+    // cells.
+    label nTotToRefine = (maxCells - globalData().nTotalCells()) / 7;
+
+    const labelList& cellLevel = meshCutter_.cellLevel();
+
+    // Count current selection
+    label nCandidates = returnReduce(count(candidateCell, 1), sumOp<label>());
+
+    // Collect all cells
+    dynamicLabelList candidates(nCells());
+
+    if (nCandidates < nTotToRefine)
+    {
+        // We won't exceed the maximum number of cells. Collect all candidate
+        // cells that have refinement level smaller than max level
+        forAll(candidateCell, cellI)
+        {
+            if
+            (
+                cellLevel[cellI] < maxRefinement
+             && candidateCell.get(cellI) == 1
+            )
+            {
+                candidates.append(cellI);
+            }
+        }
+    }
+    else
+    {
+        // We will exceed the maximum number of cells. Simply truncate the
+        // list. It is possible to prefer certain cells based on error field.
+        for (label level = 0; level < maxRefinement; level++)
+        {
+            forAll(candidateCell, cellI)
+            {
+                if
+                (
+                    cellLevel[cellI] == level
+                 && candidateCell.get(cellI) == 1
+                )
+                {
+                    candidates.append(cellI);
+                }
+            }
+
+            if (returnReduce(candidates.size(), sumOp<label>()) > nTotToRefine)
+            {
+                break;
+            }
+        }
+    }
+
+    // Guarantee 2:1 refinement after refinement
+    labelList consistentSet
+    (
+        meshCutter_.consistentRefinement
+        (
+            candidates.shrink(),
+            true               // Add to set to guarantee 2:1
+        )
+    );
+
+    Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
+        << " cells for refinement out of " << globalData().nTotalCells()
+        << "." << endl;
+
+    return consistentSet;
+}
+
+
+labelList dynamicRefinePolyFvMesh::selectUnrefinePoints
+(
+    const scalar unrefineLevel,
+    const PackedBoolList& markedCell,
+    const scalarField& pFld
+) const
+{
+    // All points that can be unrefined
+    const labelList splitPoints(meshCutter_.getSplitPoints());
+
+    dynamicLabelList newSplitPoints(splitPoints.size());
+
+    forAll(splitPoints, i)
+    {
+        label pointI = splitPoints[i];
+
+        if (pFld[pointI] < unrefineLevel)
+        {
+            // Check that all cells are not marked
+            const labelList& pCells = pointCells()[pointI];
+
+            bool hasMarked = false;
+
+            forAll(pCells, pCellI)
+            {
+                if (markedCell.get(pCells[pCellI]) == 1)
+                {
+                    hasMarked = true;
+                    break;
+                }
+            }
+
+            if (!hasMarked)
+            {
+                newSplitPoints.append(pointI);
+            }
+        }
+    }
+
+
+    newSplitPoints.shrink();
+
+    // Guarantee 2:1 refinement after unrefinement
+    labelList consistentSet
+    (
+        meshCutter_.consistentUnrefinement
+        (
+            newSplitPoints,
+            false
+        )
+    );
+    Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
+        << " split points out of a possible "
+        << returnReduce(splitPoints.size(), sumOp<label>())
+        << "." << endl;
+
+    return consistentSet;
+}
+
+
+void dynamicRefinePolyFvMesh::extendMarkedCells(PackedBoolList& markedCell) const
+{
+    // Mark faces using any marked cell
+    boolList markedFace(nFaces(), false);
+
+    forAll(markedCell, cellI)
+    {
+        if (markedCell.get(cellI) == 1)
+        {
+            const cell& cFaces = cells()[cellI];
+
+            forAll(cFaces, i)
+            {
+                markedFace[cFaces[i]] = true;
+            }
+        }
+    }
+
+    syncTools::syncFaceList(*this, markedFace, orEqOp<bool>(), false);
+
+    // Update cells using any markedFace
+    for (label faceI = 0; faceI < nInternalFaces(); faceI++)
+    {
+        if (markedFace[faceI])
+        {
+            markedCell.set(faceOwner()[faceI], 1);
+            markedCell.set(faceNeighbour()[faceI], 1);
+        }
+    }
+    for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
+    {
+        if (markedFace[faceI])
+        {
+            markedCell.set(faceOwner()[faceI], 1);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+dynamicRefinePolyFvMesh::dynamicRefinePolyFvMesh(const IOobject& io)
+:
+    dynamicFvMesh(io),
+    singleMotionUpdate_
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "dynamicMeshDict",
+                time().constant(),
+                *this,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).subDict(typeName + "Coeffs")
+        .lookupOrDefault<Switch>("singleMotionUpdate", true)
+    ),
+    curTimeIndex_(-1),
+    meshCutter_(*this),
+    dumpLevel_(false),
+    nRefinementIterations_(0)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+dynamicRefinePolyFvMesh::~dynamicRefinePolyFvMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool dynamicRefinePolyFvMesh::update()
+{
+    // Handling multiple calls in a single time step
+    if
+    (
+        singleMotionUpdate_
+     && curTimeIndex_ == this->time().timeIndex()
+    )
+    {
+        // This is not the first call to update, simply return false
+        return false;
+    }
+
+    // Update local time index
+    curTimeIndex_ = this->time().timeIndex();
+
+    // Re-read dictionary. Choosen since usually -small so trivial amount
+    // of time compared to actual refinement. Also very useful to be able
+    // to modify on-the-fly.
+    dictionary refineDict
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "dynamicMeshDict",
+                time().constant(),
+                *this,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).subDict(typeName + "Coeffs")
+    );
+
+    label refineInterval = readLabel(refineDict.lookup("refineInterval"));
+
+    bool hasChanged = false;
+
+    if (refineInterval == 0)
+    {
+        changing(hasChanged);
+
+        return false;
+    }
+    else if (refineInterval < 0)
+    {
+        FatalErrorIn("dynamicRefinePolyFvMesh::update()")
+            << "Illegal refineInterval " << refineInterval << nl
+            << "The refineInterval setting in the dynamicMeshDict should"
+            << " be >= 1." << nl
+            << exit(FatalError);
+    }
+
+
+
+
+    // Note: cannot refine at time 0 since no V0 present since mesh not
+    //       moved yet.
+
+    if (time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0)
+    {
+        label maxCells = readLabel(refineDict.lookup("maxCells"));
+
+        if (maxCells <= 0)
+        {
+            FatalErrorIn("dynamicRefinePolyFvMesh::update()")
+                << "Illegal maximum number of cells " << maxCells << nl
+                << "The maxCells setting in the dynamicMeshDict should"
+                << " be > 0." << nl
+                << exit(FatalError);
+        }
+
+        label maxRefinement = readLabel(refineDict.lookup("maxRefinement"));
+
+        if (maxRefinement <= 0)
+        {
+            FatalErrorIn("dynamicRefinePolyFvMesh::update()")
+                << "Illegal maximum refinement level " << maxRefinement << nl
+                << "The maxCells setting in the dynamicMeshDict should"
+                << " be > 0." << nl
+                << exit(FatalError);
+        }
+
+        const word fieldName(refineDict.lookup("field"));
+
+        const volScalarField& vFld = lookupObject<volScalarField>(fieldName);
+
+        const scalar lowerRefineLevel =
+            readScalar(refineDict.lookup("lowerRefineLevel"));
+        const scalar upperRefineLevel =
+            readScalar(refineDict.lookup("upperRefineLevel"));
+        const scalar unrefineLevel =
+            readScalar(refineDict.lookup("unrefineLevel"));
+        const label nBufferLayers =
+            readLabel(refineDict.lookup("nBufferLayers"));
+
+        // Cells marked for refinement or otherwise protected from unrefinement.
+        PackedBoolList refineCell(nCells());
+
+        if (globalData().nTotalCells() < maxCells)
+        {
+            // Determine candidates for refinement (looking at field only)
+            selectRefineCandidates
+            (
+                lowerRefineLevel,
+                upperRefineLevel,
+                vFld,
+                refineCell
+            );
+
+            // Select subset of candidates. Take into account max allowable
+            // cells, refinement level, protected cells.
+            labelList cellsToRefine
+            (
+                selectRefineCells
+                (
+                    maxCells,
+                    maxRefinement,
+                    refineCell
+                )
+            );
+
+            label nCellsToRefine = returnReduce
+            (
+                cellsToRefine.size(), sumOp<label>()
+            );
+
+            if (nCellsToRefine > 0)
+            {
+                // Refine/update mesh and map fields
+                autoPtr<mapPolyMesh> map = refine(cellsToRefine);
+
+                // Update refineCell. Note that some of the marked ones have
+                // not been refined due to constraints.
+                {
+                    const labelList& cellMap = map().cellMap();
+                    const labelList& reverseCellMap = map().reverseCellMap();
+
+                    PackedBoolList newRefineCell(cellMap.size());
+
+                    forAll(cellMap, cellI)
+                    {
+                        label oldCellI = cellMap[cellI];
+
+                        if (oldCellI < 0)
+                        {
+                            newRefineCell.set(cellI, 1);
+                        }
+                        else if (reverseCellMap[oldCellI] != cellI)
+                        {
+                            newRefineCell.set(cellI, 1);
+                        }
+                        else
+                        {
+                            newRefineCell.set(cellI, refineCell.get(oldCellI));
+                        }
+                    }
+                    refineCell.transfer(newRefineCell);
+                }
+
+                // Extend with a buffer layer to prevent neighbouring points
+                // being unrefined.
+                for (label i = 0; i < nBufferLayers; i++)
+                {
+                    extendMarkedCells(refineCell);
+                }
+
+                hasChanged = true;
+            }
+        }
+
+
+        {
+            // Select unrefineable points that are not marked in refineCell
+            labelList pointsToUnrefine
+            (
+                selectUnrefinePoints
+                (
+                    unrefineLevel,
+                    refineCell,
+                    minCellField(vFld)
+                )
+            );
+
+            label nSplitPoints = returnReduce
+            (
+                pointsToUnrefine.size(),
+                sumOp<label>()
+            );
+
+            if (nSplitPoints > 0)
+            {
+                // Refine/update mesh
+                unrefine(pointsToUnrefine);
+
+                hasChanged = true;
+            }
+        }
+
+
+        if ((nRefinementIterations_ % 10) == 0)
+        {
+            // Compact refinement history occassionally (how often?).
+            // Unrefinement causes holes in the polyRefinementHistory.
+            const_cast<polyRefinementHistory&>
+            (
+                meshCutter().history()
+            ).compact();
+        }
+        nRefinementIterations_++;
+    }
+
+    changing(hasChanged);
+
+    return hasChanged;
+}
+
+
+bool dynamicRefinePolyFvMesh::writeObject
+(
+    IOstream::streamFormat fmt,
+    IOstream::versionNumber ver,
+    IOstream::compressionType cmp
+) const
+{
+    // Force refinement data to go to the current time directory.
+    const_cast<polyRef&>(meshCutter_).setInstance(time().timeName());
+
+    bool writeOk =
+        dynamicFvMesh::writeObjects(fmt, ver, cmp)
+     && meshCutter_.write();
+
+    if (dumpLevel_)
+    {
+        volScalarField scalarCellLevel
+        (
+            IOobject
+            (
+                "cellLevel",
+                time().timeName(),
+                *this,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            *this,
+            dimensionedScalar("level", dimless, 0)
+        );
+
+        const labelList& cellLevel = meshCutter_.cellLevel();
+
+        forAll(cellLevel, cellI)
+        {
+            scalarCellLevel[cellI] = cellLevel[cellI];
+        }
+
+        writeOk = writeOk && scalarCellLevel.write();
+    }
+
+    return writeOk;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.H b/src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.H
new file mode 100644
index 000000000..b86029edb
--- /dev/null
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.H
@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::dynamicRefinePolyFvMesh
+
+Description
+    A fvMesh with built-in refinement of arbitrary polyhedral cells.
+
+    Determines which cells to refine/unrefine and does all in update().
+
+SourceFiles
+    dynamicRefinePolyFvMesh.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved
+
+Notes
+    Generalisation of dynamicRefineMesh for polyhedral cells
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dynamicRefinePolyFvMesh_H
+#define dynamicRefinePolyFvMesh_H
+
+#include "dynamicFvMesh.H"
+#include "polyRef.H"
+#include "PackedBoolList.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class dynamicRefinePolyFvMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class dynamicRefinePolyFvMesh
+:
+    public dynamicFvMesh
+{
+
+    // Private data
+
+        // Helper variables to enable switching between a single and multiple
+        // mesh motion updates within a time step (if update() is called more
+        // than once in a single time step)
+
+            //- Switch for single motion update (true by default)
+            Switch singleMotionUpdate_;
+
+            //- Helper varaible: current time index
+            label curTimeIndex_;
+
+
+protected:
+
+        //- Mesh cutting engine
+        polyRef meshCutter_;
+
+        //- Dump cellLevel for postprocessing
+        Switch dumpLevel_;
+
+        //- Fluxes to map
+        List<Pair<word> > correctFluxes_;
+
+        //- Number of refinement/unrefinement steps done so far.
+        label nRefinementIterations_;
+
+
+    // Private Member Functions
+
+        //- Count set/unset elements in packedlist.
+        static label count(const PackedBoolList&, const unsigned int);
+
+        //- Read the projection parameters from dictionary
+        void readDict();
+
+
+        //- Refine cells. Update mesh and fields.
+        autoPtr<mapPolyMesh> refine(const labelList&);
+
+        //- Unrefine cells. Gets passed in centre points of cells to combine.
+        autoPtr<mapPolyMesh> unrefine(const labelList&);
+
+
+        // Selection of cells to un/refine
+
+            //- Get per cell max of connected point
+            scalarField maxPointField(const scalarField&) const;
+
+            //- Get point min of connected cell
+            scalarField minCellField(const volScalarField&) const;
+
+            //- Simple (non-parallel) interpolation by averaging
+            scalarField cellToPoint(const scalarField& vFld) const;
+
+            //- Calculate error (= -1 by default or distance from inbetween
+            //  levels
+            scalarField error
+            (
+                const scalarField& fld,
+                const scalar minLevel,
+                const scalar maxLevel
+            ) const;
+
+            //- Select candidate cells for refinement
+            virtual void selectRefineCandidates
+            (
+                const scalar lowerRefineLevel,
+                const scalar upperRefineLevel,
+                const scalarField& vFld,
+                PackedBoolList& candidateCell
+            ) const;
+
+            //- Subset candidate cells for refinement
+            virtual labelList selectRefineCells
+            (
+                const label maxCells,
+                const label maxRefinement,
+                const PackedBoolList& candidateCell
+            ) const;
+
+            //- Select points that can be unrefined
+            virtual labelList selectUnrefinePoints
+            (
+                const scalar unrefineLevel,
+                const PackedBoolList& markedCell,
+                const scalarField& pFld
+            ) const;
+
+            //- Extend markedCell with cell-face-cell
+            void extendMarkedCells(PackedBoolList& markedCell) const;
+
+
+private:
+
+        //- Disallow default bitwise copy construct
+        dynamicRefinePolyFvMesh(const dynamicRefinePolyFvMesh&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const dynamicRefinePolyFvMesh&);
+
+public:
+
+    //- Runtime type information
+    TypeName("dynamicRefinePolyFvMesh");
+
+
+    // Constructors
+
+        //- Construct from IOobject
+        explicit dynamicRefinePolyFvMesh(const IOobject& io);
+
+
+    // Destructor
+
+        virtual ~dynamicRefinePolyFvMesh();
+
+
+    // Member Functions
+
+        //- Direct access to the refinement engine
+        const polyRef& meshCutter() const
+        {
+            return meshCutter_;
+        }
+
+        //- Update the mesh for both mesh motion and topology change
+        virtual bool update();
+
+
+    // Writing
+
+        //- Write using given format, version and compression
+        virtual bool writeObject
+        (
+            IOstream::streamFormat fmt,
+            IOstream::versionNumber ver,
+            IOstream::compressionType cmp
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/dynamicMesh/Make/files b/src/dynamicMesh/dynamicMesh/Make/files
index 0ec6268cc..3305c6e7e 100644
--- a/src/dynamicMesh/dynamicMesh/Make/files
+++ b/src/dynamicMesh/dynamicMesh/Make/files
@@ -52,6 +52,7 @@ motionSolver/motionSolver.C
 refinementData/refinementData.C
 refinementData/refinementDistanceData.C
 refinementData/refinementHistory.C
+refinementData/polyRefinementHistory.C
 
 directTopoChange/directTopoChange/directTopoChange.C
 
@@ -60,6 +61,7 @@ $(directActions)/addPatchCellLayer.C
 $(directActions)/edgeCollapser.C
 $(directActions)/faceCollapser.C
 $(directActions)/hexRef8.C
+$(directActions)/polyRef.C
 $(directActions)/removeCells.C
 $(directActions)/removeFaces.C
 $(directActions)/removePoints.C
diff --git a/src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/polyRef.C b/src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/polyRef.C
new file mode 100644
index 000000000..9c38356ce
--- /dev/null
+++ b/src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/polyRef.C
@@ -0,0 +1,5473 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyRef.H"
+
+#include "polyMesh.H"
+#include "directTopoChange.H"
+#include "meshTools.H"
+#include "polyAddFace.H"
+#include "polyAddPoint.H"
+#include "polyAddCell.H"
+#include "polyModifyFace.H"
+#include "syncTools.H"
+#include "faceSet.H"
+#include "cellSet.H"
+#include "pointSet.H"
+#include "OFstream.H"
+#include "foamTime.H"
+#include "FaceCellWave.H"
+#include "mapDistributePolyMesh.H"
+#include "refinementData.H"
+#include "refinementDistanceData.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(polyRef, 0);
+
+    //- Reduction class. If x and y are not equal assign value.
+    template< int value >
+    class ifEqEqOp
+    {
+        public:
+        void operator()( label& x, const label& y ) const
+        {
+            x = (x==y) ? x:value;
+        }
+    };
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::polyRef::reorder
+(
+    const labelList& map,
+    const label len,
+    const label null,
+    labelList& elems
+)
+{
+    labelList newElems(len, null);
+
+    forAll(elems, i)
+    {
+        label newI = map[i];
+
+        if (newI >= len)
+        {
+            FatalErrorIn("polyRef::reorder(..)") << abort(FatalError);
+        }
+
+        if (newI >= 0)
+        {
+            newElems[newI] = elems[i];
+        }
+    }
+
+    elems.transfer(newElems);
+}
+
+
+void Foam::polyRef::getFaceInfo
+(
+    const label faceI,
+    label& patchID,
+    label& zoneID,
+    label& zoneFlip
+) const
+{
+    patchID = -1;
+
+    if (!mesh_.isInternalFace(faceI))
+    {
+        patchID = mesh_.boundaryMesh().whichPatch(faceI);
+    }
+
+    zoneID = mesh_.faceZones().whichZone(faceI);
+
+    zoneFlip = false;
+
+    if (zoneID >= 0)
+    {
+        const faceZone& fZone = mesh_.faceZones()[zoneID];
+
+        zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
+    }
+}
+
+
+// Adds a face on top of existing faceI.
+Foam::label Foam::polyRef::addFace
+(
+    directTopoChange& meshMod,
+    const label faceI,
+    const face& newFace,
+    const label own,
+    const label nei
+) const
+{
+    label patchID, zoneID, zoneFlip;
+
+    getFaceInfo(faceI, patchID, zoneID, zoneFlip);
+
+    label newFaceI = -1;
+
+    if ((nei == -1) || (own < nei))
+    {
+        // Ordering ok
+        newFaceI = meshMod.setAction
+        (
+            polyAddFace
+            (
+                newFace,                    // face
+                own,                        // owner
+                nei,                        // neighbour
+                -1,                         // master point
+                -1,                         // master edge
+                faceI,                      // master face for addition
+                false,                      // flux flip
+                patchID,                    // patch for face
+                zoneID,                     // zone for face
+                zoneFlip                    // face zone flip
+            )
+        );
+    }
+    else
+    {
+        // Reverse owner/neighbour
+        newFaceI = meshMod.setAction
+        (
+            polyAddFace
+            (
+                newFace.reverseFace(),      // face
+                nei,                        // owner
+                own,                        // neighbour
+                -1,                         // master point
+                -1,                         // master edge
+                faceI,                      // master face for addition
+                false,                      // flux flip
+                patchID,                    // patch for face
+                zoneID,                     // zone for face
+                zoneFlip                    // face zone flip
+            )
+        );
+    }
+
+    return newFaceI;
+}
+
+
+// Adds an internal face from an edge. Assumes orientation correct.
+// Problem is that the face is between four new vertices. So what do we provide
+// as master? The only existing mesh item we have is the edge we have split.
+// Have to be careful in only using it if it has internal faces since otherwise
+// polyMeshMorph will complain (because it cannot generate a sensible mapping
+// for the face)
+Foam::label Foam::polyRef::addInternalFace
+(
+    directTopoChange& meshMod,
+    const label meshFaceI,
+    const label meshPointI,
+    const face& newFace,
+    const label own,
+    const label nei
+) const
+{
+    if (mesh_.isInternalFace(meshFaceI))
+    {
+        return meshMod.setAction
+        (
+            polyAddFace
+            (
+                newFace,                    // face
+                own,                        // owner
+                nei,                        // neighbour
+                -1,                         // master point
+                -1,                         // master edge
+                meshFaceI,                  // master face for addition
+                false,                      // flux flip
+                -1,                         // patch for face
+                -1,                         // zone for face
+                false                       // face zone flip
+            )
+        );
+    }
+    else
+    {
+        // Two choices:
+        // - append (i.e. create out of nothing - will not be mapped)
+        //   problem: field does not get mapped.
+        // - inflate from point.
+        //   problem: does interpolative mapping which constructs full
+        //   volPointInterpolation!
+
+        // For now create out of nothing
+
+        // Map from zero face to achieve some initialisation and
+        // avoid uninitialised memory in mapping
+        // HJ, 4/Dec/2013
+        return meshMod.setAction
+        (
+            polyAddFace
+            (
+                newFace,                    // face
+                own,                        // owner
+                nei,                        // neighbour
+                -1,                         // master point
+                -1,                         // master edge
+                0,                          // master face for addition
+                false,                      // flux flip
+                -1,                         // patch for face
+                -1,                         // zone for face
+                false                       // face zone flip
+            )
+        );
+
+
+        ////- Inflate-from-point:
+        //// Check if point has any internal faces we can use.
+        //label masterPointI = -1;
+        //
+        //const labelList& pFaces = mesh_.pointFaces()[meshPointI];
+        //
+        //forAll(pFaces, i)
+        //{
+        //    if (mesh_.isInternalFace(pFaces[i]))
+        //    {
+        //        // meshPoint uses internal faces so ok to inflate from it
+        //        masterPointI = meshPointI;
+        //
+        //        break;
+        //    }
+        //}
+        //
+        //return meshMod.setAction
+        //(
+        //    polyAddFace
+        //    (
+        //        newFace,                    // face
+        //        own,                        // owner
+        //        nei,                        // neighbour
+        //        masterPointI,               // master point
+        //        -1,                         // master edge
+        //        -1,                         // master face for addition
+        //        false,                      // flux flip
+        //        -1,                         // patch for face
+        //        -1,                         // zone for face
+        //        false                       // face zone flip
+        //    )
+        //);
+    }
+}
+
+
+// Modifies existing faceI for either new owner/neighbour or new face points.
+void Foam::polyRef::modFace
+(
+    directTopoChange& meshMod,
+    const label faceI,
+    const face& newFace,
+    const label own,
+    const label nei
+) const
+{
+    label patchID, zoneID, zoneFlip;
+
+    getFaceInfo(faceI, patchID, zoneID, zoneFlip);
+
+    if
+    (
+        (own != mesh_.faceOwner()[faceI])
+     || (
+            mesh_.isInternalFace(faceI)
+         && (nei != mesh_.faceNeighbour()[faceI])
+        )
+     || (newFace != mesh_.faces()[faceI])
+    )
+    {
+        if ((nei == -1) || (own < nei))
+        {
+            // Ordering ok
+            meshMod.setAction
+            (
+                polyModifyFace
+                (
+                    newFace,            // modified face
+                    faceI,              // label of face being modified
+                    own,                // owner
+                    nei,                // neighbour
+                    false,              // face flip
+                    patchID,            // patch for face
+                    false,              // remove from zone
+                    zoneID,             // zone for face
+                    zoneFlip            // face flip in zone
+                )
+            );
+        }
+        else
+        {
+            // Reverse owner/neighbour
+            meshMod.setAction
+            (
+                polyModifyFace
+                (
+                    newFace.reverseFace(),  // modified face
+                    faceI,                  // label of face being modified
+                    nei,                    // owner
+                    own,                    // neighbour
+                    false,                  // face flip
+                    patchID,                // patch for face
+                    false,                  // remove from zone
+                    zoneID,                 // zone for face
+                    zoneFlip                // face flip in zone
+                )
+            );
+        }
+    }
+}
+
+
+// Bit complex way to determine the unrefined edge length.
+Foam::scalar Foam::polyRef::getLevel0EdgeLength() const
+{
+    if (cellLevel_.size() != mesh_.nCells())
+    {
+        FatalErrorIn("scalar polyRef::getLevel0EdgeLength() const")
+            << "Number of cells in mesh:" << mesh_.nCells()
+            << " does not equal size of cellLevel:" << cellLevel_.size()
+            << endl
+            << "This might be because of a restart with inconsistent cellLevel."
+            << abort(FatalError);
+    }
+
+    // Determine minimum edge length per refinement level
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    const scalar GREAT2 = sqr(GREAT);
+
+    label nLevels = gMax(cellLevel_) + 1;
+
+    scalarField typEdgeLenSqr(nLevels, GREAT2);
+
+    // 1. Look only at edges surrounded by cellLevel cells only.
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    {
+        // Per edge the cellLevel of connected cells. -1 if not set,
+        // labelMax if different levels, otherwise levels of connected cells.
+        labelList edgeLevel(mesh_.nEdges(), -1);
+
+        forAll(cellLevel_, cellI)
+        {
+            const label cLevel = cellLevel_[cellI];
+
+            const labelList& cEdges = mesh_.cellEdges()[cellI];
+
+            forAll(cEdges, i)
+            {
+                label edgeI = cEdges[i];
+
+                if (edgeLevel[edgeI] == -1)
+                {
+                    edgeLevel[edgeI] = cLevel;
+                }
+                else if (edgeLevel[edgeI] == labelMax)
+                {
+                    // Already marked as on different cellLevels
+                }
+                else if (edgeLevel[edgeI] != cLevel)
+                {
+                    edgeLevel[edgeI] = labelMax;
+                }
+            }
+        }
+
+        // Make sure that edges with different levels on different processors
+        // are also marked. Do the same test (edgeLevel != cLevel) on coupled
+        // edges.
+        syncTools::syncEdgeList
+        (
+            mesh_,
+            edgeLevel,
+            ifEqEqOp<labelMax>(),
+            labelMin,
+            false               // no separation
+        );
+
+        // Now use the edgeLevel with a valid value to determine the
+        // length per level.
+        forAll(edgeLevel, edgeI)
+        {
+            const label eLevel = edgeLevel[edgeI];
+
+            if (eLevel >= 0 && eLevel < labelMax)
+            {
+                const edge& e = mesh_.edges()[edgeI];
+
+                scalar edgeLenSqr = magSqr(e.vec(mesh_.points()));
+
+                typEdgeLenSqr[eLevel] = min(typEdgeLenSqr[eLevel], edgeLenSqr);
+            }
+        }
+    }
+
+    // Get the minimum per level over all processors. Note minimum so if
+    // cells are not cubic we use the smallest edge side.
+    Pstream::listCombineGather(typEdgeLenSqr, minEqOp<scalar>());
+    Pstream::listCombineScatter(typEdgeLenSqr);
+
+    if (debug)
+    {
+        Pout<< "polyRef::getLevel0EdgeLength() :"
+            << " After phase1: Edgelengths (squared) per refinementlevel:"
+            << typEdgeLenSqr << endl;
+    }
+
+
+    // 2. For any levels where we haven't determined a valid length yet
+    //    use any surrounding cell level. Here we use the max so we don't
+    //    pick up levels between celllevel and higher celllevel (will have
+    //    edges sized according to highest celllevel)
+    //    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    scalarField maxEdgeLenSqr(nLevels, -GREAT2);
+
+    forAll(cellLevel_, cellI)
+    {
+        const label cLevel = cellLevel_[cellI];
+
+        const labelList& cEdges = mesh_.cellEdges()[cellI];
+
+        forAll(cEdges, i)
+        {
+            const edge& e = mesh_.edges()[cEdges[i]];
+
+            scalar edgeLenSqr = magSqr(e.vec(mesh_.points()));
+
+            maxEdgeLenSqr[cLevel] = max(maxEdgeLenSqr[cLevel], edgeLenSqr);
+        }
+    }
+
+    Pstream::listCombineGather(maxEdgeLenSqr, maxEqOp<scalar>());
+    Pstream::listCombineScatter(maxEdgeLenSqr);
+
+    if (debug)
+    {
+        Pout<< "polyRef::getLevel0EdgeLength() :"
+            << " Crappy Edgelengths (squared) per refinementlevel:"
+            << maxEdgeLenSqr << endl;
+    }
+
+
+    // 3. Combine the two sets of lengths
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    forAll(typEdgeLenSqr, levelI)
+    {
+        if (typEdgeLenSqr[levelI] == GREAT2 && maxEdgeLenSqr[levelI] >= 0)
+        {
+            typEdgeLenSqr[levelI] = maxEdgeLenSqr[levelI];
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "polyRef::getLevel0EdgeLength() :"
+            << " Final Edgelengths (squared) per refinementlevel:"
+            << typEdgeLenSqr << endl;
+    }
+
+    // Find lowest level present
+    scalar level0Size = -1;
+
+    forAll(typEdgeLenSqr, levelI)
+    {
+        scalar lenSqr = typEdgeLenSqr[levelI];
+
+        if (lenSqr < GREAT2)
+        {
+            level0Size = Foam::sqrt(lenSqr)*(1<<levelI);
+
+            if (debug)
+            {
+                Pout<< "polyRef::getLevel0EdgeLength() :"
+                    << " For level:" << levelI
+                    << " have edgeLen:" << Foam::sqrt(lenSqr)
+                    << " with equivalent level0 len:" << level0Size
+                    << endl;
+            }
+            break;
+        }
+    }
+
+    if (level0Size == -1)
+    {
+        FatalErrorIn("polyRef::getLevel0EdgeLength()")
+            << "Problem : typEdgeLenSqr:" << typEdgeLenSqr << abort(FatalError);
+    }
+
+    return level0Size;
+}
+
+
+// Check whether pointI is an anchor on cellI.
+// If it is not check whether any other point on the face is an anchor cell.
+Foam::label Foam::polyRef::getAnchorCell
+(
+    const labelListList& cellAnchorPoints,
+    const labelListList& cellAddedCells,
+    const label cellI,
+    const label faceI,
+    const label pointI
+) const
+{
+    if (cellAnchorPoints[cellI].size() > 0)
+    {
+        label index = findIndex(cellAnchorPoints[cellI], pointI);
+
+        if (index != -1)
+        {
+            return cellAddedCells[cellI][index];
+        }
+
+
+        // pointI is not an anchor cell.
+        // Maybe we are already a refined face so check all the face
+        // vertices.
+        const face& f = mesh_.faces()[faceI];
+
+        forAll(f, fp)
+        {
+            label index = findIndex(cellAnchorPoints[cellI], f[fp]);
+
+            if (index != -1)
+            {
+                return cellAddedCells[cellI][index];
+            }
+        }
+
+        // Problem.
+
+        // Pick up points of the cell
+        const labelList cPoints(cellPoints(cellI));
+
+        Perr<< "cell: " << cellI << ", points: " << endl;
+        forAll(cPoints, i)
+        {
+            label pointI = cPoints[i];
+
+            Perr<< "    " << pointI << " coord: " << mesh_.points()[pointI]
+                << nl;
+        }
+
+        Perr<< "cell: " << cellI << " anchorPoints: " << cellAnchorPoints[cellI]
+            << endl;
+
+        FatalErrorIn("polyRef::getAnchorCell(..)")
+            << "Could not find point " << pointI
+            << " in the anchorPoints for cell " << cellI << endl
+            << "Does your original mesh obey the 2:1 constraint and"
+            << " did you use consistentRefinement to make your cells to refine"
+            << " obey this constraint as well?"
+            << abort(FatalError);
+
+        return -1;
+    }
+    else
+    {
+        return cellI;
+    }
+}
+
+
+// Get new owner and neighbour
+void Foam::polyRef::getFaceNeighbours
+(
+    const labelListList& cellAnchorPoints,
+    const labelListList& cellAddedCells,
+    const label faceI,
+    const label pointI,
+
+    label& own,
+    label& nei
+) const
+{
+    // Is owner split?
+    own = getAnchorCell
+    (
+        cellAnchorPoints,
+        cellAddedCells,
+        mesh_.faceOwner()[faceI],
+        faceI,
+        pointI
+    );
+
+    if (mesh_.isInternalFace(faceI))
+    {
+        nei = getAnchorCell
+        (
+            cellAnchorPoints,
+            cellAddedCells,
+            mesh_.faceNeighbour()[faceI],
+            faceI,
+            pointI
+        );
+    }
+    else
+    {
+        nei = -1;
+    }
+}
+
+
+// Get point with the lowest pointLevel
+Foam::label Foam::polyRef::findMinLevel(const labelList& f) const
+{
+    label minLevel = labelMax;
+    label minFp = -1;
+
+    forAll(f, fp)
+    {
+        const label level = pointLevel_[f[fp]];
+
+        if (level < minLevel)
+        {
+            minLevel = level;
+            minFp = fp;
+        }
+    }
+
+    return minFp;
+}
+
+
+// Get point with the highest pointLevel
+Foam::label Foam::polyRef::findMaxLevel(const labelList& f) const
+{
+    label maxLevel = labelMin;
+    label maxFp = -1;
+
+    forAll(f, fp)
+    {
+        const label level = pointLevel_[f[fp]];
+
+        if (level > maxLevel)
+        {
+            maxLevel = level;
+            maxFp = fp;
+        }
+    }
+
+    return maxFp;
+}
+
+
+Foam::label Foam::polyRef::countAnchors
+(
+    const labelList& f,
+    const label anchorLevel
+) const
+{
+    label nAnchors = 0;
+
+    forAll(f, fp)
+    {
+        if (pointLevel_[f[fp]] <= anchorLevel)
+        {
+            nAnchors++;
+        }
+    }
+    return nAnchors;
+}
+
+
+// Find point with certain pointLevel. Skip any higher levels.
+Foam::label Foam::polyRef::findLevel
+(
+    const face& f,
+    const label startFp,
+    const bool searchForward,
+    const label wantedLevel
+) const
+{
+    label fp = startFp;
+
+    forAll(f, i)
+    {
+        label pointI = f[fp];
+
+        if (pointLevel_[pointI] < wantedLevel)
+        {
+            FatalErrorIn("polyRef::findLevel")
+                << "face:" << f
+                << " level:" << IndirectList<label>(pointLevel_, f)()
+                << " startFp:" << startFp
+                << " wantedLevel:" << wantedLevel
+                << abort(FatalError);
+        }
+        else if (pointLevel_[pointI] == wantedLevel)
+        {
+            return fp;
+        }
+
+        if (searchForward)
+        {
+            fp = f.fcIndex(fp);
+        }
+        else
+        {
+            fp = f.rcIndex(fp);
+        }
+    }
+
+    FatalErrorIn("polyRef::findLevel")
+        << "face:" << f
+        << " level:" << IndirectList<label>(pointLevel_, f)()
+        << " startFp:" << startFp
+        << " wantedLevel:" << wantedLevel
+        << abort(FatalError);
+
+    return -1;
+}
+
+
+// Get cell level such that the face has the most points that are <= level (at
+// least three points)
+Foam::label Foam::polyRef::getAnchorLevel(const label faceI) const
+{
+    const face& f = mesh_.faces()[faceI];
+
+    // Get number of points in this face
+    const label nPoints = f.size();
+
+    // Get unique levels that can be found in this face
+    labelHashSet uniqueLevels(nPoints);
+    forAll(f, fp)
+    {
+        uniqueLevels.insert(pointLevel_[f[fp]]);
+    }
+
+    // Get sorted levels (in increasing order)
+    const labelList allLevels = uniqueLevels.sortedToc();
+
+    // Count number of points per level
+    labelList nPointsPerLevel(allLevels.size(), 0);
+
+    // Loop through all points
+    forAll(f, fp)
+    {
+        const label curPointLevel = pointLevel_[f[fp]];
+
+        // Loop through all levels
+        forAll(allLevels, levelI)
+        {
+            // Get current level
+            const label curLevel = allLevels[levelI];
+
+            if (curPointLevel == curLevel)
+            {
+                // Note: storing nPoints per level based on list index, not
+                // actual level
+                ++nPointsPerLevel[levelI];
+
+                break;
+            }
+        }
+    }
+
+    if (debug)
+    {
+        // Print out number of points per level
+        Info<< "nPointsPerLevel: " << nPointsPerLevel << endl;
+    }
+
+    // Find level with maximum number of points
+    const label maxLevelI = findMax(nPointsPerLevel);
+
+    // Check whether there are at least three points with the level <= of this
+    // level with maximum number of occurences
+    label nPointsBelow = 0;
+    for (label belowLevelI = 0; belowLevelI <= maxLevelI; ++belowLevelI)
+    {
+        nPointsBelow += nPointsPerLevel[belowLevelI];
+    }
+
+    if (nPointsBelow < 3)
+    {
+        // For some reason, the minimum requirement that at least three points
+        // have level below the anchor level is not met
+        WarningIn("label polyRef::getAnchorLevel(const label faceI) const")
+            << "Found less than three points below an anchor level."
+            << nl
+            << "Face: " << faceI
+            << ", anchor level: " << allLevels[maxLevelI]
+            << ", n points below level: " << nPointsBelow
+            << ", n total points of this face: " << nPoints
+            << endl;
+    }
+
+    return allLevels[maxLevelI];
+}
+
+
+void Foam::polyRef::checkInternalOrientation
+(
+    directTopoChange& meshMod,
+    const label cellI,
+    const label faceI,
+    const point& ownPt,
+    const point& neiPt,
+    const face& newFace
+)
+{
+    const face compactFace(identity(newFace.size()));
+    const pointField compactPoints
+    (
+        IndirectList<point>(meshMod.points(), newFace)()
+    );
+
+    const vector n(compactFace.normal(compactPoints));
+
+    const vector dir(neiPt - ownPt);
+
+    // Check orientation error
+    if ((dir & n) < 0)
+    {
+        FatalErrorIn("checkInternalOrientation(..)")
+            << "cell:" << cellI << " old face:" << faceI
+            << " newFace:" << newFace << endl
+            << " coords:" << compactPoints
+            << " ownPt:" << ownPt
+            << " neiPt:" << neiPt
+            << abort(FatalError);
+    }
+
+    // Note: report significant non-orthogonality error
+    const scalar severeNonOrthogonalityThreshold =
+      ::cos
+        (
+            primitiveMesh::nonOrthThreshold_()/180.0*mathematicalConstant::pi
+        );
+
+    const vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
+
+    // Note: normal vector already normalised
+    const scalar dDotN = (fcToOwn & n)/(mag(fcToOwn) + VSMALL);
+
+    if (dDotN > severeNonOrthogonalityThreshold)
+    {
+        WarningIn("checkInternalOrientation(..)")
+            << "Detected severely non-orthogonal face with non-orthogonality: "
+            << ::acos(dDotN)/mathematicalConstant::pi*180.0
+            << "cell:" << cellI << " old face:" << faceI
+            << " newFace:" << newFace << endl
+            << " coords:" << compactPoints
+            << " ownPt:" << ownPt
+            << " neiPt:" << neiPt
+            << " dDotN:" << dDotN
+            << endl;
+    }
+}
+
+
+void Foam::polyRef::checkBoundaryOrientation
+(
+    directTopoChange& meshMod,
+    const label cellI,
+    const label faceI,
+    const point& ownPt,
+    const point& boundaryPt,
+    const face& newFace
+)
+{
+    const face compactFace(identity(newFace.size()));
+    const pointField compactPoints
+    (
+        IndirectList<point>(meshMod.points(), newFace)()
+    );
+
+    const vector n(compactFace.normal(compactPoints));
+
+    const vector dir(boundaryPt - ownPt);
+
+    // Check orientation error
+    if ((dir & n) < 0)
+    {
+        FatalErrorIn("checkBoundaryOrientation(..)")
+            << "cell:" << cellI << " old face:" << faceI
+            << " newFace:" << newFace
+            << " coords:" << compactPoints
+            << " ownPt:" << ownPt
+            << " boundaryPt:" << boundaryPt
+            << abort(FatalError);
+    }
+
+    // Note: report significant non-orthogonality error
+    const scalar severeNonOrthogonalityThreshold =
+      ::cos
+        (
+            primitiveMesh::nonOrthThreshold_()/180.0*mathematicalConstant::pi
+        );
+
+    const vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
+
+    // Note: normal vector already normalised
+    const scalar dDotN = (fcToOwn & n)/(mag(fcToOwn) + VSMALL);
+
+    if (dDotN > severeNonOrthogonalityThreshold)
+    {
+        WarningIn("checkBoundaryOrientation(..)")
+            << "Detected severely non-orthogonal face with non-orthogonality: "
+            << ::acos(dDotN)/mathematicalConstant::pi*180.0
+            << "cell:" << cellI << " old face:" << faceI
+            << " newFace:" << newFace
+            << " coords:" << compactPoints
+            << " ownPt:" << ownPt
+            << " boundaryPt:" << boundaryPt
+            << " dDotN:" << dDotN
+            << endl;
+    }
+}
+
+
+// If p0 and p1 are existing vertices check if edge is split and insert
+// splitPoint.
+void Foam::polyRef::insertEdgeSplit
+(
+    const labelList& edgeMidPoint,
+    const label p0,
+    const label p1,
+    dynamicLabelList& verts
+) const
+{
+    if (p0 < mesh_.nPoints() && p1 < mesh_.nPoints())
+    {
+        label edgeI = meshTools::findEdge(mesh_, p0, p1);
+
+        if (edgeI != -1 && edgeMidPoint[edgeI] != -1)
+        {
+            verts.append(edgeMidPoint[edgeI]);
+        }
+    }
+}
+
+
+// Internal faces are one per edge between anchor points. So one per midPoint
+// between the anchor points. Here we store the information on the midPoint
+// and if we have enough information:
+// - two anchors
+// - two face mid points
+// we add the face. Note that this routine can get called anywhere from
+// two times (two unrefined faces) to four times (two refined faces) so
+// the first call that adds the information creates the face.
+Foam::label Foam::polyRef::storeMidPointInfo
+(
+    const labelListList& cellAnchorPoints,
+    const labelListList& cellAddedCells,
+    const labelList& cellMidPoint,
+    const labelList& edgeMidPoint,
+    const label cellI,
+    const label faceI,
+    const bool faceOrder,
+    const label edgeMidPointI,
+    const label anchorPointI,
+    const label faceMidPointI,
+
+    Map<edge>& midPointToAnchors,
+    Map<edge>& midPointToFaceMids,
+    directTopoChange& meshMod
+) const
+{
+    // See if need to store anchors.
+
+    bool changed = false;
+    bool haveTwoAnchors = false;
+
+    Map<edge>::iterator edgeMidFnd = midPointToAnchors.find(edgeMidPointI);
+
+    if (edgeMidFnd == midPointToAnchors.end())
+    {
+        midPointToAnchors.insert(edgeMidPointI, edge(anchorPointI, -1));
+    }
+    else
+    {
+        edge& e = edgeMidFnd();
+
+        if (anchorPointI != e[0])
+        {
+            if (e[1] == -1)
+            {
+                e[1] = anchorPointI;
+                changed = true;
+            }
+        }
+
+        if (e[0] != -1 && e[1] != -1)
+        {
+            haveTwoAnchors = true;
+        }
+    }
+
+    bool haveTwoFaceMids = false;
+
+    Map<edge>::iterator faceMidFnd = midPointToFaceMids.find(edgeMidPointI);
+
+    if (faceMidFnd == midPointToFaceMids.end())
+    {
+        midPointToFaceMids.insert(edgeMidPointI, edge(faceMidPointI, -1));
+    }
+    else
+    {
+        edge& e = faceMidFnd();
+
+        if (faceMidPointI != e[0])
+        {
+            if (e[1] == -1)
+            {
+                e[1] = faceMidPointI;
+                changed = true;
+            }
+        }
+
+        if (e[0] != -1 && e[1] != -1)
+        {
+            haveTwoFaceMids = true;
+        }
+    }
+
+    // Check if this call of storeMidPointInfo is the one that completed all
+    // the nessecary information.
+
+    if (changed && haveTwoAnchors && haveTwoFaceMids)
+    {
+        const edge& anchors = midPointToAnchors[edgeMidPointI];
+        const edge& faceMids = midPointToFaceMids[edgeMidPointI];
+
+        label otherFaceMidPointI = faceMids.otherVertex(faceMidPointI);
+
+        // Create face consistent with anchorI being the owner.
+        // Note that the edges between the edge mid point and the face mids
+        // might be marked for splitting. Note that these edge splits cannot
+        // be between cellMid and face mids.
+
+        dynamicLabelList newFaceVerts(4);
+        if (faceOrder == (mesh_.faceOwner()[faceI] == cellI))
+        {
+            newFaceVerts.append(faceMidPointI);
+
+            // Check & insert edge split if any
+            insertEdgeSplit
+            (
+                edgeMidPoint,
+                faceMidPointI,  // edge between faceMid and
+                edgeMidPointI,  // edgeMid
+                newFaceVerts
+            );
+
+            newFaceVerts.append(edgeMidPointI);
+
+            insertEdgeSplit
+            (
+                edgeMidPoint,
+                edgeMidPointI,
+                otherFaceMidPointI,
+                newFaceVerts
+            );
+
+            newFaceVerts.append(otherFaceMidPointI);
+            newFaceVerts.append(cellMidPoint[cellI]);
+        }
+        else
+        {
+            newFaceVerts.append(otherFaceMidPointI);
+
+            insertEdgeSplit
+            (
+                edgeMidPoint,
+                otherFaceMidPointI,
+                edgeMidPointI,
+                newFaceVerts
+            );
+
+            newFaceVerts.append(edgeMidPointI);
+
+            insertEdgeSplit
+            (
+                edgeMidPoint,
+                edgeMidPointI,
+                faceMidPointI,
+                newFaceVerts
+            );
+
+            newFaceVerts.append(faceMidPointI);
+            newFaceVerts.append(cellMidPoint[cellI]);
+        }
+
+        face newFace;
+        newFace.transfer(newFaceVerts.shrink());
+        newFaceVerts.clear();
+
+        label anchorCell0 = getAnchorCell
+        (
+            cellAnchorPoints,
+            cellAddedCells,
+            cellI,
+            faceI,
+            anchorPointI
+        );
+        label anchorCell1 = getAnchorCell
+        (
+            cellAnchorPoints,
+            cellAddedCells,
+            cellI,
+            faceI,
+            anchors.otherVertex(anchorPointI)
+        );
+
+        // Get mesh points
+        const pointField& meshPoints = mesh_.points();
+
+        label own, nei;
+        point ownPt, neiPt;
+
+        if (anchorCell0 < anchorCell1)
+        {
+            own = anchorCell0;
+            nei = anchorCell1;
+
+            ownPt = meshPoints[anchorPointI];
+            neiPt = meshPoints[anchors.otherVertex(anchorPointI)];
+
+        }
+        else
+        {
+            own = anchorCell1;
+            nei = anchorCell0;
+            newFace = newFace.reverseFace();
+
+            ownPt = meshPoints[anchors.otherVertex(anchorPointI)];
+            neiPt = meshPoints[anchorPointI];
+        }
+
+        if (debug)
+        {
+            point ownPt, neiPt;
+
+            if (anchorCell0 < anchorCell1)
+            {
+                ownPt = meshPoints[anchorPointI];
+                neiPt = meshPoints[anchors.otherVertex(anchorPointI)];
+            }
+            else
+            {
+                ownPt = meshPoints[anchors.otherVertex(anchorPointI)];
+                neiPt = meshPoints[anchorPointI];
+            }
+
+            checkInternalOrientation
+            (
+                meshMod,
+                cellI,
+                faceI,
+                ownPt,
+                neiPt,
+                newFace
+            );
+        }
+
+        return addInternalFace
+        (
+            meshMod,
+            faceI,
+            anchorPointI,
+            newFace,
+            own,
+            nei
+        );
+    }
+    else
+    {
+        return -1;
+    }
+}
+
+
+// Creates all internal faces for cellI.
+void Foam::polyRef::createInternalFaces
+(
+    const labelListList& cellAnchorPoints,
+    const labelListList& cellAddedCells,
+    const labelList& cellMidPoint,
+    const labelList& faceMidPoint,
+    const labelList& faceAnchorLevel,
+    const labelList& edgeMidPoint,
+    const label cellI,
+
+    directTopoChange& meshMod
+) const
+{
+    // Find in every face the cellLevel + 1 points (from edge subdivision)
+    // and the anchor points.
+
+    const cell& cFaces = mesh_.cells()[cellI];
+    const label cLevel = cellLevel_[cellI];
+
+    // From edge mid to anchor points
+    Map<edge> midPointToAnchors(24);
+    // From edge mid to face mids
+    Map<edge> midPointToFaceMids(24);
+
+    // Running count of number of internal faces added so far.
+    label nFacesAdded = 0;
+
+    forAll(cFaces, i)
+    {
+        label faceI = cFaces[i];
+
+        const face& f = mesh_.faces()[faceI];
+        const labelList& fEdges = mesh_.faceEdges()[faceI];
+
+        // We are on the cellI side of face f. The face will have 1 or n
+        // cLevel points
+
+        label faceMidPointI = -1;
+
+        const label nAnchors = countAnchors(f, cLevel);
+
+        if (nAnchors == 1)
+        {
+            // Only one anchor point. So the other side of the face has already
+            // been split using cLevel + 1 and cLevel + 2 points.
+
+            // Find the one anchor.
+            label anchorFp = -1;
+
+            forAll(f, fp)
+            {
+                if (pointLevel_[f[fp]] <= cLevel)
+                {
+                    anchorFp = fp;
+                    break;
+                }
+            }
+
+            // Now the face mid point is the second cLevel + 1 point
+            label edgeMid = findLevel(f, f.fcIndex(anchorFp), true, cLevel + 1);
+            label faceMid = findLevel(f, f.fcIndex(edgeMid), true, cLevel + 1);
+
+            faceMidPointI = f[faceMid];
+        }
+        else
+        {
+            // There is no face middle yet but the face will be marked for
+            // splitting.
+
+            faceMidPointI = faceMidPoint[faceI];
+        }
+
+
+        // Now loop over all the anchors (might be just one) and store
+        // the edge mids connected to it. storeMidPointInfo will collect
+        // all the info and combine it all.
+
+        forAll(f, fp0)
+        {
+            const label point0 = f[fp0];
+
+            if (pointLevel_[point0] <= cLevel)
+            {
+                // Anchor.
+
+                // Walk forward
+                // ~~~~~~~~~~~~
+                // to cLevel + 1 or edgeMidPoint of this level.
+
+
+                label edgeMidPointI = -1;
+
+                const label fp1 = f.fcIndex(fp0);
+
+                if (pointLevel_[f[fp1]] <= cLevel)
+                {
+                    // Anchor. Edge will be split.
+                    const label edgeI = fEdges[fp0];
+
+                    edgeMidPointI = edgeMidPoint[edgeI];
+
+                    if (edgeMidPointI == -1)
+                    {
+                        const labelList cPoints(cellPoints(cellI));
+
+                        FatalErrorIn("createInternalFaces(..)")
+                            << "cell:" << cellI << " cLevel:" << cLevel
+                            << " cell points:" << cPoints
+                            << " pointLevel:"
+                            << IndirectList<label>(pointLevel_, cPoints)()
+                            << " face:" << faceI
+                            << " f:" << f
+                            << " pointLevel:"
+                            << IndirectList<label>(pointLevel_, f)()
+                            << " faceAnchorLevel:" << faceAnchorLevel[faceI]
+                            << " faceMidPoint:" << faceMidPoint[faceI]
+                            << " faceMidPointI:" << faceMidPointI
+                            << " fp:" << fp0
+                            << abort(FatalError);
+                    }
+                }
+                else
+                {
+                    // Search forward in face to clevel + 1
+                    const label edgeMid = findLevel(f, fp1, true, cLevel + 1);
+
+                    edgeMidPointI = f[edgeMid];
+                }
+
+                label newFaceI = storeMidPointInfo
+                (
+                    cellAnchorPoints,
+                    cellAddedCells,
+                    cellMidPoint,
+                    edgeMidPoint,
+
+                    cellI,
+                    faceI,
+                    true,                   // mid point after anchor
+                    edgeMidPointI,          // edgemid
+                    point0,                 // anchor
+                    faceMidPointI,
+
+                    midPointToAnchors,
+                    midPointToFaceMids,
+                    meshMod
+                );
+
+                if (newFaceI != -1)
+                {
+                    ++nFacesAdded;
+                }
+
+
+
+                // Walk backward
+                // ~~~~~~~~~~~~~
+
+                label fpMin1 = f.rcIndex(fp0);
+
+                if (pointLevel_[f[fpMin1]] <= cLevel)
+                {
+                    // Anchor. Edge will be split.
+                    label edgeI = fEdges[fpMin1];
+
+                    edgeMidPointI = edgeMidPoint[edgeI];
+
+                    if (edgeMidPointI == -1)
+                    {
+                        const labelList cPoints(cellPoints(cellI));
+
+                        FatalErrorIn("createInternalFaces(..)")
+                            << "cell:" << cellI << " cLevel:" << cLevel
+                            << " cell points:" << cPoints
+                            << " pointLevel:"
+                            << IndirectList<label>(pointLevel_, cPoints)()
+                            << " face:" << faceI
+                            << " f:" << f
+                            << " pointLevel:"
+                            << IndirectList<label>(pointLevel_, f)()
+                            << " faceAnchorLevel:" << faceAnchorLevel[faceI]
+                            << " faceMidPoint:" << faceMidPoint[faceI]
+                            << " faceMidPointI:" << faceMidPointI
+                            << " fp:" << fp0
+                            << abort(FatalError);
+                    }
+                }
+                else
+                {
+                    // Search back to clevel + 1
+                    label edgeMid = findLevel(f, fpMin1, false, cLevel + 1);
+
+                    edgeMidPointI = f[edgeMid];
+                }
+
+                newFaceI = storeMidPointInfo
+                (
+                    cellAnchorPoints,
+                    cellAddedCells,
+                    cellMidPoint,
+                    edgeMidPoint,
+
+                    cellI,
+                    faceI,
+                    false,                  // mid point before anchor
+                    edgeMidPointI,          // edgemid
+                    point0,                 // anchor
+                    faceMidPointI,
+
+                    midPointToAnchors,
+                    midPointToFaceMids,
+                    meshMod
+                );
+
+                if (newFaceI != -1)
+                {
+                    ++nFacesAdded;
+                }
+            } // End for this anchor point
+        } // End for all face points
+    } // End for all cell faces
+}
+
+
+void Foam::polyRef::walkFaceToMid
+(
+    const labelList& edgeMidPoint,
+    const label cLevel,
+    const label faceI,
+    const label startFp,
+    dynamicLabelList& faceVerts
+) const
+{
+    const face& f = mesh_.faces()[faceI];
+    const labelList& fEdges = mesh_.faceEdges()[faceI];
+
+    label fp = startFp;
+
+    // Starting from fp store all (1 or 2) vertices until where the face
+    // gets split
+
+    while (true)
+    {
+        if (edgeMidPoint[fEdges[fp]] >= 0)
+        {
+            faceVerts.append(edgeMidPoint[fEdges[fp]]);
+        }
+
+        fp = f.fcIndex(fp);
+
+        if (pointLevel_[f[fp]] <= cLevel)
+        {
+            // Next anchor. Have already append split point on edge in code
+            // above.
+            return;
+        }
+        else if (pointLevel_[f[fp]] == cLevel + 1)
+        {
+            // Mid level
+            faceVerts.append(f[fp]);
+
+            return;
+        }
+        else if (pointLevel_[f[fp]] == cLevel + 2)
+        {
+            // Store and continue to cLevel + 1.
+            faceVerts.append(f[fp]);
+        }
+    }
+}
+
+
+// Same as walkFaceToMid but now walk back.
+void Foam::polyRef::walkFaceFromMid
+(
+    const labelList& edgeMidPoint,
+    const label cLevel,
+    const label faceI,
+    const label startFp,
+    dynamicLabelList& faceVerts
+) const
+{
+    const face& f = mesh_.faces()[faceI];
+    const labelList& fEdges = mesh_.faceEdges()[faceI];
+
+    label fp = f.rcIndex(startFp);
+
+    while (true)
+    {
+        if (pointLevel_[f[fp]] <= cLevel)
+        {
+            // anchor.
+            break;
+        }
+        else if (pointLevel_[f[fp]] == cLevel + 1)
+        {
+            // Mid level
+            faceVerts.append(f[fp]);
+            break;
+        }
+        else if (pointLevel_[f[fp]] == cLevel + 2)
+        {
+            // Continue to cLevel+1.
+        }
+        fp = f.rcIndex(fp);
+    }
+
+    // Store
+    while (true)
+    {
+        if (edgeMidPoint[fEdges[fp]] >= 0)
+        {
+            faceVerts.append(edgeMidPoint[fEdges[fp]]);
+        }
+
+        fp = f.fcIndex(fp);
+
+        if (fp == startFp)
+        {
+            break;
+        }
+        faceVerts.append(f[fp]);
+    }
+}
+
+
+// Updates refineCell (cells marked for refinement) so across all faces
+// there will be 2:1 consistency after refinement.
+Foam::label Foam::polyRef::faceConsistentRefinement
+(
+    const bool maxSet,
+    PackedList<1>& refineCell
+) const
+{
+    label nChanged = 0;
+
+    // Internal faces.
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        label own = mesh_.faceOwner()[faceI];
+        label ownLevel = cellLevel_[own] + refineCell.get(own);
+
+        label nei = mesh_.faceNeighbour()[faceI];
+        label neiLevel = cellLevel_[nei] + refineCell.get(nei);
+
+        if (ownLevel > (neiLevel+1))
+        {
+            if (maxSet)
+            {
+                refineCell.set(nei, 1);
+            }
+            else
+            {
+                refineCell.set(own, 0);
+            }
+            nChanged++;
+        }
+        else if (neiLevel > (ownLevel+1))
+        {
+            if (maxSet)
+            {
+                refineCell.set(own, 1);
+            }
+            else
+            {
+                refineCell.set(nei, 0);
+            }
+            nChanged++;
+        }
+    }
+
+
+    // Coupled faces. Swap owner level to get neighbouring cell level.
+    // (only boundary faces of neiLevel used)
+    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+
+    forAll(neiLevel, i)
+    {
+        label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
+
+        neiLevel[i] = cellLevel_[own] + refineCell.get(own);
+    }
+
+    // Swap to neighbour
+    syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
+
+    // Now we have neighbour value see which cells need refinement
+    forAll(neiLevel, i)
+    {
+        label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
+        label ownLevel = cellLevel_[own] + refineCell.get(own);
+
+        if (ownLevel > (neiLevel[i]+1))
+        {
+            if (!maxSet)
+            {
+                refineCell.set(own, 0);
+                nChanged++;
+            }
+        }
+        else if (neiLevel[i] > (ownLevel+1))
+        {
+            if (maxSet)
+            {
+                refineCell.set(own, 1);
+                nChanged++;
+            }
+        }
+    }
+
+    return nChanged;
+}
+
+
+// Debug: check if wanted refinement is compatible with 2:1
+void Foam::polyRef::checkWantedRefinementLevels
+(
+    const labelList& cellsToRefine
+) const
+{
+    PackedList<1> refineCell(mesh_.nCells(), 0);
+    forAll(cellsToRefine, i)
+    {
+        refineCell.set(cellsToRefine[i], 1);
+    }
+
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        label own = mesh_.faceOwner()[faceI];
+        label ownLevel = cellLevel_[own] + refineCell.get(own);
+
+        label nei = mesh_.faceNeighbour()[faceI];
+        label neiLevel = cellLevel_[nei] + refineCell.get(nei);
+
+        if (mag(ownLevel-neiLevel) > 1)
+        {
+            FatalErrorIn
+            (
+                "polyRef::checkWantedRefinementLevels(const labelList&)"
+            )   << "cell:" << own
+                << " current level:" << cellLevel_[own]
+                << " level after refinement:" << ownLevel
+                << nl
+                << "neighbour cell:" << nei
+                << " current level:" << cellLevel_[nei]
+                << " level after refinement:" << neiLevel
+                << nl
+                << "which does not satisfy 2:1 constraints anymore."
+                << abort(FatalError);
+        }
+    }
+
+    // Coupled faces. Swap owner level to get neighbouring cell level.
+    // (only boundary faces of neiLevel used)
+    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+
+    forAll(neiLevel, i)
+    {
+        label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
+
+        neiLevel[i] = cellLevel_[own] + refineCell.get(own);
+    }
+
+    // Swap to neighbour
+    syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
+
+    // Now we have neighbour value see which cells need refinement
+    forAll(neiLevel, i)
+    {
+        label faceI = i + mesh_.nInternalFaces();
+
+        label own = mesh_.faceOwner()[faceI];
+        label ownLevel = cellLevel_[own] + refineCell.get(own);
+
+        if (mag(ownLevel - neiLevel[i]) > 1)
+        {
+            label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+
+            FatalErrorIn
+            (
+                "polyRef::checkWantedRefinementLevels(const labelList&)"
+            )   << "Cell level does not satisfy 2:1 constraint."
+                << " On coupled face "
+                << faceI
+                << " on patch " << patchI << " "
+                << mesh_.boundaryMesh()[patchI].name()
+                << " owner cell " << own
+                << " current level:" << cellLevel_[own]
+                << " level after refinement:" << ownLevel
+                << nl
+                << " (coupled) neighbour cell will get refinement "
+                << neiLevel[i]
+                << abort(FatalError);
+        }
+    }
+}
+
+
+// Set instance for mesh files
+void Foam::polyRef::setInstance(const fileName& inst)
+{
+    if (debug)
+    {
+        Pout<< "polyRef::setInstance(const fileName& inst) : "
+            << "Resetting file instance to " << inst << endl;
+    }
+
+    cellLevel_.instance() = inst;
+    pointLevel_.instance() = inst;
+    history_.instance() = inst;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from mesh, read refinement data
+Foam::polyRef::polyRef(const polyMesh& mesh)
+:
+    mesh_(mesh),
+    cellLevel_
+    (
+        IOobject
+        (
+            "cellLevel",
+            mesh_.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        labelList(mesh_.nCells(), 0)
+    ),
+    pointLevel_
+    (
+        IOobject
+        (
+            "pointLevel",
+            mesh_.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        labelList(mesh_.nPoints(), 0)
+    ),
+    level0Edge_(getLevel0EdgeLength()),
+    history_
+    (
+        IOobject
+        (
+            "polyRefinementHistory",
+            mesh_.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_.nCells()    // All cells visible if could not be read
+    ),
+    faceRemover_(mesh_, GREAT),     // merge boundary faces wherever possible
+    savedPointLevel_(0),
+    savedCellLevel_(0)
+{
+    if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
+    {
+        FatalErrorIn
+        (
+            "polyRef::polyRef(const polyMesh&)"
+        )   << "History enabled but number of visible cells in it "
+            << history_.visibleCells().size()
+            << " is not equal to the number of cells in the mesh "
+            << mesh_.nCells()
+            << abort(FatalError);
+    }
+
+    if
+    (
+        cellLevel_.size() != mesh_.nCells()
+     || pointLevel_.size() != mesh_.nPoints()
+    )
+    {
+        FatalErrorIn
+        (
+            "polyRef::polyRef(const polyMesh&)"
+        )   << "Restarted from inconsistent cellLevel or pointLevel files."
+            << endl
+            << "Number of cells in mesh:" << mesh_.nCells()
+            << " does not equal size of cellLevel:" << cellLevel_.size() << nl
+            << "Number of points in mesh:" << mesh_.nPoints()
+            << " does not equal size of pointLevel:" << pointLevel_.size()
+            << abort(FatalError);
+    }
+
+
+    // Check refinement levels for consistency
+    checkRefinementLevels(-1, labelList(0));
+
+
+    // Check initial mesh for consistency
+
+    //if (debug)
+    {
+        checkMesh();
+    }
+}
+
+
+// Construct from components
+Foam::polyRef::polyRef
+(
+    const polyMesh& mesh,
+    const labelList& cellLevel,
+    const labelList& pointLevel,
+    const polyRefinementHistory& history
+)
+:
+    mesh_(mesh),
+    cellLevel_
+    (
+        IOobject
+        (
+            "cellLevel",
+            mesh_.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        cellLevel
+    ),
+    pointLevel_
+    (
+        IOobject
+        (
+            "pointLevel",
+            mesh_.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pointLevel
+    ),
+    level0Edge_(getLevel0EdgeLength()),
+    history_
+    (
+        IOobject
+        (
+            "polyRefinementHistory",
+            mesh_.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        history
+    ),
+    faceRemover_(mesh_, GREAT),     // merge boundary faces wherever possible
+    savedPointLevel_(0),
+    savedCellLevel_(0)
+{
+    if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
+    {
+        FatalErrorIn
+        (
+            "polyRef::polyRef(const polyMesh&, const labelList&"
+            ", const labelList&, const polyRefinementHistory&)"
+        )   << "History enabled but number of visible cells in it "
+            << history_.visibleCells().size()
+            << " is not equal to the number of cells in the mesh "
+            << mesh_.nCells() << abort(FatalError);
+    }
+
+    if
+    (
+        cellLevel_.size() != mesh_.nCells()
+     || pointLevel_.size() != mesh_.nPoints()
+    )
+    {
+        FatalErrorIn
+        (
+            "polyRef::polyRef(const polyMesh&, const labelList&"
+            ", const labelList&, const polyRefinementHistory&)"
+        )   << "Incorrect cellLevel or pointLevel size." << endl
+            << "Number of cells in mesh:" << mesh_.nCells()
+            << " does not equal size of cellLevel:" << cellLevel_.size() << endl
+            << "Number of points in mesh:" << mesh_.nPoints()
+            << " does not equal size of pointLevel:" << pointLevel_.size()
+            << abort(FatalError);
+    }
+
+    // Check refinement levels for consistency
+    checkRefinementLevels(-1, labelList(0));
+
+
+    // Check initial mesh for consistency
+
+    //if (debug)
+    {
+        checkMesh();
+    }
+}
+
+
+// Construct from components
+Foam::polyRef::polyRef
+(
+    const polyMesh& mesh,
+    const labelList& cellLevel,
+    const labelList& pointLevel
+)
+:
+    mesh_(mesh),
+    cellLevel_
+    (
+        IOobject
+        (
+            "cellLevel",
+            mesh_.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        cellLevel
+    ),
+    pointLevel_
+    (
+        IOobject
+        (
+            "pointLevel",
+            mesh_.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pointLevel
+    ),
+    level0Edge_(getLevel0EdgeLength()),
+    history_
+    (
+        IOobject
+        (
+            "polyRefinementHistory",
+            mesh_.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        List<polyRefinementHistory::splitPolyCell>(0),
+        labelList(0)
+    ),
+    faceRemover_(mesh_, GREAT),     // merge boundary faces wherever possible
+    savedPointLevel_(0),
+    savedCellLevel_(0)
+{
+    if
+    (
+        cellLevel_.size() != mesh_.nCells()
+     || pointLevel_.size() != mesh_.nPoints()
+    )
+    {
+        FatalErrorIn
+        (
+            "polyRef::polyRef(const polyMesh&, const labelList&"
+            ", const labelList&)"
+        )   << "Incorrect cellLevel or pointLevel size." << endl
+            << "Number of cells in mesh:" << mesh_.nCells()
+            << " does not equal size of cellLevel:" << cellLevel_.size() << endl
+            << "Number of points in mesh:" << mesh_.nPoints()
+            << " does not equal size of pointLevel:" << pointLevel_.size()
+            << abort(FatalError);
+    }
+
+    // Check refinement levels for consistency
+    checkRefinementLevels(-1, labelList(0));
+
+    // Check initial mesh for consistency
+
+    //if (debug)
+    {
+        checkMesh();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//- Get points of a cell (without using cellPoints addressing)
+Foam::labelList Foam::polyRef::cellPoints(const label cellI) const
+{
+    // Pick up points of the cell
+    const cell& cFaces = mesh_.cells()[cellI];
+
+    labelHashSet cPoints(4*cFaces.size());
+
+    forAll(cFaces, i)
+    {
+        const face& f = mesh_.faces()[cFaces[i]];
+
+        forAll(f, fp)
+        {
+            cPoints.insert(f[fp]);
+        }
+    }
+    return cPoints.toc();
+}
+
+
+Foam::labelList Foam::polyRef::consistentRefinement
+(
+    const labelList& cellsToRefine,
+    const bool maxSet
+) const
+{
+    // Loop, modifying cellsToRefine, until no more changes to due to 2:1
+    // conflicts.
+    // maxSet = false : unselect cells to refine
+    // maxSet = true  : select cells to refine
+
+    // Go to straight boolList.
+    PackedList<1> refineCell(mesh_.nCells(), 0);
+    forAll(cellsToRefine, i)
+    {
+        refineCell.set(cellsToRefine[i], 1);
+    }
+
+    while (true)
+    {
+        label nChanged = faceConsistentRefinement(maxSet, refineCell);
+
+        reduce(nChanged, sumOp<label>());
+
+        if (debug)
+        {
+            Pout<< "polyRef::consistentRefinement : Changed " << nChanged
+                << " refinement levels due to 2:1 conflicts."
+                << endl;
+        }
+
+        if (nChanged == 0)
+        {
+            break;
+        }
+    }
+
+
+    // Convert back to labelList.
+    label nRefined = 0;
+
+    forAll(refineCell, cellI)
+    {
+        if (refineCell.get(cellI) == 1)
+        {
+            nRefined++;
+        }
+    }
+
+    labelList newCellsToRefine(nRefined);
+    nRefined = 0;
+
+    forAll(refineCell, cellI)
+    {
+        if (refineCell.get(cellI) == 1)
+        {
+            newCellsToRefine[nRefined++] = cellI;
+        }
+    }
+
+    if (debug)
+    {
+        checkWantedRefinementLevels(newCellsToRefine);
+    }
+
+    return newCellsToRefine;
+}
+
+
+// Given a list of cells to refine determine additional cells to refine
+// such that the overall refinement:
+// - satisfies maxFaceDiff (e.g. 2:1) across neighbouring faces
+// - satisfies maxPointDiff (e.g. 4:1) across selected point connected
+//   cells. This is used to ensure that e.g. cells on the surface are not
+//   point connected to cells which are 8 times smaller.
+Foam::labelList Foam::polyRef::consistentSlowRefinement
+(
+    const label maxFaceDiff,
+    const labelList& cellsToRefine,
+    const labelList& facesToCheck,
+    const label maxPointDiff,
+    const labelList& pointsToCheck
+) const
+{
+    const labelList& faceOwner = mesh_.faceOwner();
+    const labelList& faceNeighbour = mesh_.faceNeighbour();
+
+
+    if (maxFaceDiff <= 0)
+    {
+        FatalErrorIn
+        (
+            "polyRef::consistentSlowRefinement"
+            "(const label, const labelList&, const labelList&"
+            ", const label, const labelList&)"
+        )   << "Illegal maxFaceDiff " << maxFaceDiff << nl
+            << "Value should be >= 1" << exit(FatalError);
+    }
+
+
+    // Bit tricky. Say we want a distance of three cells between two
+    // consecutive refinement levels. This is done by using FaceCellWave to
+    // transport out the new refinement level. It gets decremented by one
+    // every cell it crosses so if we initialize it to maxFaceDiff
+    // we will get a field everywhere that tells us whether an unselected cell
+    // needs refining as well.
+
+
+    // Initial information about (distance to) cellLevel on all cells
+    List<refinementData> allCellInfo(mesh_.nCells());
+
+    // Initial information about (distance to) cellLevel on all faces
+    List<refinementData> allFaceInfo(mesh_.nFaces());
+
+    forAll(allCellInfo, cellI)
+    {
+        // maxFaceDiff since refinementData counts both
+        // faces and cells.
+        allCellInfo[cellI] = refinementData
+        (
+            maxFaceDiff*(cellLevel_[cellI]+1),// when cell is to be refined
+            maxFaceDiff*cellLevel_[cellI]     // current level
+        );
+    }
+
+    // Cells to be refined will have cellLevel+1
+    forAll(cellsToRefine, i)
+    {
+        label cellI = cellsToRefine[i];
+
+        allCellInfo[cellI].count() = allCellInfo[cellI].refinementCount();
+    }
+
+
+    // Labels of seed faces
+    dynamicLabelList seedFaces(mesh_.nFaces()/100);
+    // refinementLevel data on seed faces
+    DynamicList<refinementData> seedFacesInfo(mesh_.nFaces()/100);
+
+
+    // Additional buffer layer thickness by changing initial count. Usually
+    // this happens on boundary faces. Bit tricky. Use allFaceInfo to mark
+    // off thus marked faces so they're skipped in the next loop.
+    forAll(facesToCheck, i)
+    {
+        label faceI = facesToCheck[i];
+
+        if (allFaceInfo[faceI].valid())
+        {
+            // Can only occur if face has already gone through loop below.
+            FatalErrorIn
+            (
+                "polyRef::consistentSlowRefinement"
+                "(const label, const labelList&, const labelList&"
+                ", const label, const labelList&)"
+            )   << "Argument facesToCheck seems to have duplicate entries!"
+                << endl
+                << "face:" << faceI << " occurs at positions "
+                << findIndices(facesToCheck, faceI)
+                << abort(FatalError);
+        }
+
+
+        const refinementData& ownData = allCellInfo[faceOwner[faceI]];
+
+        if (mesh_.isInternalFace(faceI))
+        {
+            // Seed face if neighbouring cell (after possible refinement)
+            // will be refined one more than the current owner or neighbour.
+
+            const refinementData& neiData = allCellInfo[faceNeighbour[faceI]];
+
+            label faceCount;
+            label faceRefineCount;
+            if (neiData.count() > ownData.count())
+            {
+                faceCount = neiData.count() + maxFaceDiff;
+                faceRefineCount = faceCount + maxFaceDiff;
+            }
+            else
+            {
+                faceCount = ownData.count() + maxFaceDiff;
+                faceRefineCount = faceCount + maxFaceDiff;
+            }
+
+            seedFaces.append(faceI);
+            seedFacesInfo.append
+            (
+                refinementData
+                (
+                    faceRefineCount,
+                    faceCount
+                )
+            );
+            allFaceInfo[faceI] = seedFacesInfo[seedFacesInfo.size()-1];
+        }
+        else
+        {
+            label faceCount = ownData.count() + maxFaceDiff;
+            label faceRefineCount = faceCount + maxFaceDiff;
+
+            seedFaces.append(faceI);
+            seedFacesInfo.append
+            (
+                refinementData
+                (
+                    faceRefineCount,
+                    faceCount
+                )
+            );
+            allFaceInfo[faceI] = seedFacesInfo[seedFacesInfo.size()-1];
+        }
+    }
+
+
+    // Just seed with all faces inbetween different refinement levels for now
+    // (alternatively only seed faces on cellsToRefine but that gives problems
+    //  if no cells to refine)
+    forAll(faceNeighbour, faceI)
+    {
+        // Check if face already handled in loop above
+        if (!allFaceInfo[faceI].valid())
+        {
+            label own = faceOwner[faceI];
+            label nei = faceNeighbour[faceI];
+
+            // Seed face with transported data from highest cell.
+
+            if (allCellInfo[own].count() > allCellInfo[nei].count())
+            {
+                allFaceInfo[faceI].updateFace
+                (
+                    mesh_,
+                    faceI,
+                    own,
+                    allCellInfo[own],
+                    FaceCellWave<refinementData>::propagationTol()
+                );
+                seedFaces.append(faceI);
+                seedFacesInfo.append(allFaceInfo[faceI]);
+            }
+            else if (allCellInfo[own].count() < allCellInfo[nei].count())
+            {
+                allFaceInfo[faceI].updateFace
+                (
+                    mesh_,
+                    faceI,
+                    nei,
+                    allCellInfo[nei],
+                    FaceCellWave<refinementData>::propagationTol()
+                );
+                seedFaces.append(faceI);
+                seedFacesInfo.append(allFaceInfo[faceI]);
+            }
+        }
+    }
+
+    // Seed all boundary faces with owner value. This is to make sure that
+    // they are visited (probably only important for coupled faces since
+    // these need to be visited from both sides)
+    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+    {
+        // Check if face already handled in loop above
+        if (!allFaceInfo[faceI].valid())
+        {
+            label own = faceOwner[faceI];
+
+            // Seed face with transported data from owner.
+            refinementData faceData;
+            faceData.updateFace
+            (
+                mesh_,
+                faceI,
+                own,
+                allCellInfo[own],
+                FaceCellWave<refinementData>::propagationTol()
+            );
+            seedFaces.append(faceI);
+            seedFacesInfo.append(faceData);
+        }
+    }
+
+
+    // face-cell-face transport engine
+    FaceCellWave<refinementData> levelCalc
+    (
+        mesh_,
+        allFaceInfo,
+        allCellInfo
+    );
+
+    while(true)
+    {
+        if (debug)
+        {
+            Pout<< "polyRef::consistentSlowRefinement : Seeded "
+                << seedFaces.size() << " faces between cells with different"
+                << " refinement level." << endl;
+        }
+
+        // Set seed faces
+        levelCalc.setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink());
+        seedFaces.clear();
+        seedFacesInfo.clear();
+
+        // Iterate until no change. Now 2:1 face difference should be satisfied
+        levelCalc.iterate(mesh_.globalData().nTotalFaces());
+
+
+        // Now check point-connected cells (face-connected cells already ok):
+        // - get per point max of connected cells
+        // - sync across coupled points
+        // - check cells against above point max
+
+        if (maxPointDiff == -1)
+        {
+            // No need to do any point checking.
+            break;
+        }
+
+        // Determine per point the max cell level. (done as count, not
+        // as cell level purely for ease)
+        labelList maxPointCount(mesh_.nPoints(), 0);
+
+        const labelListList& pointCells = mesh_.pointCells();
+
+        forAll(pointCells, pointI)
+        {
+            label& pLevel = maxPointCount[pointI];
+
+            const labelList& pCells = pointCells[pointI];
+
+            forAll(pCells, i)
+            {
+                pLevel = max(pLevel, allCellInfo[pCells[i]].count());
+            }
+        }
+
+        // Sync maxPointCount to neighbour
+        syncTools::syncPointList
+        (
+            mesh_,
+            maxPointCount,
+            maxEqOp<label>(),
+            labelMin,           // null value
+            false               // no separation
+        );
+
+        // Update allFaceInfo from maxPointCount for all points to check
+        // (usually on boundary faces)
+
+        // Per face the new refinement data
+        Map<refinementData> changedFacesInfo(pointsToCheck.size());
+
+        forAll(pointsToCheck, i)
+        {
+            label pointI = pointsToCheck[i];
+
+            // Loop over all cells using the point and check whether their
+            // refinement level is much less than the maximum.
+
+            const labelList& pCells = pointCells[pointI];
+
+            forAll(pCells, pCellI)
+            {
+                label cellI = pCells[pCellI];
+
+                refinementData& cellInfo = allCellInfo[cellI];
+
+                if
+                (
+                   !cellInfo.isRefined()
+                 && (
+                        maxPointCount[pointI]
+                      > cellInfo.count() + maxFaceDiff*maxPointDiff
+                    )
+                )
+                {
+                    // Mark cell for refinement
+                    cellInfo.count() = cellInfo.refinementCount();
+
+                    // Insert faces of cell as seed faces.
+                    const cell& cFaces = mesh_.cells()[cellI];
+
+                    forAll(cFaces, cFaceI)
+                    {
+                        label faceI = cFaces[cFaceI];
+
+                        refinementData faceData;
+                        faceData.updateFace
+                        (
+                            mesh_,
+                            faceI,
+                            cellI,
+                            cellInfo,
+                            FaceCellWave<refinementData>::propagationTol()
+                        );
+
+                        if (faceData.count() > allFaceInfo[faceI].count())
+                        {
+                            changedFacesInfo.insert(faceI, faceData);
+                        }
+                    }
+                }
+            }
+        }
+
+        label nChanged = changedFacesInfo.size();
+        reduce(nChanged, sumOp<label>());
+
+        if (nChanged == 0)
+        {
+            break;
+        }
+
+
+        // Transfer into seedFaces, seedFacesInfo
+        seedFaces.setSize(changedFacesInfo.size());
+        seedFacesInfo.setSize(changedFacesInfo.size());
+
+        forAllConstIter(Map<refinementData>, changedFacesInfo, iter)
+        {
+            seedFaces.append(iter.key());
+            seedFacesInfo.append(iter());
+        }
+    }
+
+
+    if (debug)
+    {
+        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+        {
+            label own = mesh_.faceOwner()[faceI];
+            label ownLevel =
+                cellLevel_[own]
+              + (allCellInfo[own].isRefined() ? 1 : 0);
+
+            label nei = mesh_.faceNeighbour()[faceI];
+            label neiLevel =
+                cellLevel_[nei]
+              + (allCellInfo[nei].isRefined() ? 1 : 0);
+
+            if (mag(ownLevel-neiLevel) > 1)
+            {
+                FatalErrorIn
+                (
+                    "polyRef::consistentSlowRefinement"
+                )   << "cell:" << own
+                    << " current level:" << cellLevel_[own]
+                    << " current refData:" << allCellInfo[own]
+                    << " level after refinement:" << ownLevel
+                    << nl
+                    << "neighbour cell:" << nei
+                    << " current level:" << cellLevel_[nei]
+                    << " current refData:" << allCellInfo[nei]
+                    << " level after refinement:" << neiLevel
+                    << nl
+                    << "which does not satisfy 2:1 constraints anymore." << nl
+                    << "face:" << faceI << " faceRefData:" << allFaceInfo[faceI]
+                    << abort(FatalError);
+            }
+        }
+
+
+        // Coupled faces. Swap owner level to get neighbouring cell level.
+        // (only boundary faces of neiLevel used)
+
+        labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+        labelList neiCount(mesh_.nFaces()-mesh_.nInternalFaces());
+        labelList neiRefCount(mesh_.nFaces()-mesh_.nInternalFaces());
+
+        forAll(neiLevel, i)
+        {
+            label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
+            neiLevel[i] = cellLevel_[own];
+            neiCount[i] = allCellInfo[own].count();
+            neiRefCount[i] = allCellInfo[own].refinementCount();
+        }
+
+        // Swap to neighbour
+        syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
+        syncTools::swapBoundaryFaceList(mesh_, neiCount, false);
+        syncTools::swapBoundaryFaceList(mesh_, neiRefCount, false);
+
+        // Now we have neighbour value see which cells need refinement
+        forAll(neiLevel, i)
+        {
+            label faceI = i+mesh_.nInternalFaces();
+
+            label own = mesh_.faceOwner()[faceI];
+            label ownLevel =
+                cellLevel_[own]
+              + (allCellInfo[own].isRefined() ? 1 : 0);
+
+            label nbrLevel =
+                neiLevel[i]
+              + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0);
+
+            if (mag(ownLevel - nbrLevel) > 1)
+            {
+                label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+
+                FatalErrorIn
+                (
+                    "polyRef::consistentSlowRefinement"
+                    "(const label, const labelList&, const labelList&"
+                    ", const label, const labelList&)"
+                )   << "Celllevel does not satisfy 2:1 constraint."
+                    << " On coupled face "
+                    << faceI
+                    << " refData:" << allFaceInfo[faceI]
+                    << " on patch " << patchI << " "
+                    << mesh_.boundaryMesh()[patchI].name() << nl
+                    << "owner cell " << own
+                    << " current level:" << cellLevel_[own]
+                    << " current count:" << allCellInfo[own].count()
+                    << " current refCount:"
+                    << allCellInfo[own].refinementCount()
+                    << " level after refinement:" << ownLevel
+                    << nl
+                    << "(coupled) neighbour cell"
+                    << " has current level:" << neiLevel[i]
+                    << " current count:" << neiCount[i]
+                    << " current refCount:" << neiRefCount[i]
+                    << " level after refinement:" << nbrLevel
+                    << abort(FatalError);
+            }
+        }
+    }
+
+    // Convert back to labelList of cells to refine.
+
+    label nRefined = 0;
+
+    forAll(allCellInfo, cellI)
+    {
+        if (allCellInfo[cellI].isRefined())
+        {
+            nRefined++;
+        }
+    }
+
+    // Updated list of cells to refine
+    labelList newCellsToRefine(nRefined);
+    nRefined = 0;
+
+    forAll(allCellInfo, cellI)
+    {
+        if (allCellInfo[cellI].isRefined())
+        {
+            newCellsToRefine[nRefined++] = cellI;
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "polyRef::consistentSlowRefinement : From "
+            << cellsToRefine.size() << " to " << newCellsToRefine.size()
+            << " cells to refine." << endl;
+    }
+
+    return newCellsToRefine;
+}
+
+
+Foam::labelList Foam::polyRef::consistentSlowRefinement2
+(
+    const label maxFaceDiff,
+    const labelList& cellsToRefine,
+    const labelList& facesToCheck
+) const
+{
+    const labelList& faceOwner = mesh_.faceOwner();
+    const labelList& faceNeighbour = mesh_.faceNeighbour();
+
+    if (maxFaceDiff <= 0)
+    {
+        FatalErrorIn
+        (
+            "polyRef::consistentSlowRefinement2"
+            "(const label, const labelList&, const labelList&)"
+        )   << "Illegal maxFaceDiff " << maxFaceDiff << nl
+            << "Value should be >= 1" << exit(FatalError);
+    }
+
+    const scalar level0Size = 2*maxFaceDiff*level0Edge_;
+
+
+    // Bit tricky. Say we want a distance of three cells between two
+    // consecutive refinement levels. This is done by using FaceCellWave to
+    // transport out the 'refinement shell'. Anything inside the refinement
+    // shell (given by a distance) gets marked for refinement.
+
+    // Initial information about (distance to) cellLevel on all cells
+    List<refinementDistanceData> allCellInfo(mesh_.nCells());
+
+    // Initial information about (distance to) cellLevel on all faces
+    List<refinementDistanceData> allFaceInfo(mesh_.nFaces());
+
+
+    // Mark cells with wanted refinement level
+    forAll(cellsToRefine, i)
+    {
+        label cellI = cellsToRefine[i];
+
+        allCellInfo[cellI] = refinementDistanceData
+        (
+            level0Size,
+            mesh_.cellCentres()[cellI],
+            cellLevel_[cellI]+1             // wanted refinement
+        );
+    }
+    // Mark all others with existing refinement level
+    forAll(allCellInfo, cellI)
+    {
+        if (!allCellInfo[cellI].valid())
+        {
+            allCellInfo[cellI] = refinementDistanceData
+            (
+                level0Size,
+                mesh_.cellCentres()[cellI],
+                cellLevel_[cellI]           // wanted refinement
+            );
+        }
+    }
+
+
+    // Labels of seed faces
+    dynamicLabelList seedFaces(mesh_.nFaces()/100);
+    // refinementLevel data on seed faces
+    DynamicList<refinementDistanceData> seedFacesInfo(mesh_.nFaces()/100);
+
+
+    const pointField& cc = mesh_.cellCentres();
+
+    forAll(facesToCheck, i)
+    {
+        label faceI = facesToCheck[i];
+
+        if (allFaceInfo[faceI].valid())
+        {
+            // Can only occur if face has already gone through loop below.
+            FatalErrorIn
+            (
+                "polyRef::consistentSlowRefinement2"
+                "(const label, const labelList&, const labelList&)"
+            )   << "Argument facesToCheck seems to have duplicate entries!"
+                << endl
+                << "face:" << faceI << " occurs at positions "
+                << findIndices(facesToCheck, faceI)
+                << abort(FatalError);
+        }
+
+        label own = faceOwner[faceI];
+
+        label ownLevel =
+        (
+            allCellInfo[own].valid()
+          ? allCellInfo[own].originLevel()
+          : cellLevel_[own]
+        );
+
+        if (!mesh_.isInternalFace(faceI))
+        {
+            // Do as if boundary face would have neighbour with one higher
+            // refinement level.
+            const point& fc = mesh_.faceCentres()[faceI];
+
+            refinementDistanceData neiData
+            (
+                level0Size,
+                2*fc - cc[own],    // est'd cell centre
+                ownLevel+1
+            );
+
+            allFaceInfo[faceI].updateFace
+            (
+                mesh_,
+                faceI,
+                own,        // not used (should be nei)
+                neiData,
+                FaceCellWave<refinementDistanceData>::propagationTol()
+            );
+        }
+        else
+        {
+            label nei = faceNeighbour[faceI];
+
+            label neiLevel =
+            (
+                allCellInfo[nei].valid()
+              ? allCellInfo[nei].originLevel()
+              : cellLevel_[nei]
+            );
+
+            if (ownLevel == neiLevel)
+            {
+                // Fake as if nei>own or own>nei (whichever one 'wins')
+                allFaceInfo[faceI].updateFace
+                (
+                    mesh_,
+                    faceI,
+                    nei,
+                    refinementDistanceData(level0Size, cc[nei], neiLevel+1),
+                    FaceCellWave<refinementDistanceData>::propagationTol()
+                );
+                allFaceInfo[faceI].updateFace
+                (
+                    mesh_,
+                    faceI,
+                    own,
+                    refinementDistanceData(level0Size, cc[own], ownLevel+1),
+                    FaceCellWave<refinementDistanceData>::propagationTol()
+                );
+            }
+            else
+            {
+                // Difference in level anyway.
+                allFaceInfo[faceI].updateFace
+                (
+                    mesh_,
+                    faceI,
+                    nei,
+                    refinementDistanceData(level0Size, cc[nei], neiLevel),
+                    FaceCellWave<refinementDistanceData>::propagationTol()
+                );
+                allFaceInfo[faceI].updateFace
+                (
+                    mesh_,
+                    faceI,
+                    own,
+                    refinementDistanceData(level0Size, cc[own], ownLevel),
+                    FaceCellWave<refinementDistanceData>::propagationTol()
+                );
+            }
+        }
+        seedFaces.append(faceI);
+        seedFacesInfo.append(allFaceInfo[faceI]);
+    }
+
+
+    // Create some initial seeds to start walking from. This is only if there
+    // are no facesToCheck.
+    // Just seed with all faces inbetween different refinement levels for now
+    forAll(faceNeighbour, faceI)
+    {
+        // Check if face already handled in loop above
+        if (!allFaceInfo[faceI].valid())
+        {
+            label own = faceOwner[faceI];
+
+            label ownLevel =
+            (
+                allCellInfo[own].valid()
+              ? allCellInfo[own].originLevel()
+              : cellLevel_[own]
+            );
+
+            label nei = faceNeighbour[faceI];
+
+            label neiLevel =
+            (
+                allCellInfo[nei].valid()
+              ? allCellInfo[nei].originLevel()
+              : cellLevel_[nei]
+            );
+
+            if (ownLevel > neiLevel)
+            {
+                // Set face to owner data. (since face not yet would be copy)
+                seedFaces.append(faceI);
+                allFaceInfo[faceI].updateFace
+                (
+                    mesh_,
+                    faceI,
+                    own,
+                    refinementDistanceData(level0Size, cc[own], ownLevel),
+                    FaceCellWave<refinementDistanceData>::propagationTol()
+                );
+                seedFacesInfo.append(allFaceInfo[faceI]);
+            }
+            else if (neiLevel > ownLevel)
+            {
+                seedFaces.append(faceI);
+                allFaceInfo[faceI].updateFace
+                (
+                    mesh_,
+                    faceI,
+                    nei,
+                    refinementDistanceData(level0Size, cc[nei], neiLevel),
+                    FaceCellWave<refinementDistanceData>::propagationTol()
+                );
+                seedFacesInfo.append(allFaceInfo[faceI]);
+            }
+        }
+    }
+
+    seedFaces.shrink();
+    seedFacesInfo.shrink();
+
+    // face-cell-face transport engine
+    FaceCellWave<refinementDistanceData> levelCalc
+    (
+        mesh_,
+        seedFaces,
+        seedFacesInfo,
+        allFaceInfo,
+        allCellInfo,
+        mesh_.globalData().nTotalCells()
+    );
+
+
+    //if (debug)
+    //{
+    //    // Dump wanted level
+    //    volScalarField wantedLevel
+    //    (
+    //        IOobject
+    //        (
+    //            "wantedLevel",
+    //            fMesh.time().timeName(),
+    //            fMesh,
+    //            IOobject::NO_READ,
+    //            IOobject::AUTO_WRITE,
+    //            false
+    //        ),
+    //        fMesh,
+    //        dimensionedScalar("zero", dimless, 0)
+    //    );
+    //
+    //    forAll(wantedLevel, cellI)
+    //    {
+    //        wantedLevel[cellI] = allCellInfo[cellI].wantedLevel(cc[cellI]);
+    //    }
+    //
+    //    Pout<< "Writing " << wantedLevel.objectPath() << endl;
+    //    wantedLevel.write();
+    //}
+
+
+    // Convert back to labelList of cells to refine.
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // 1. Force original refinement cells to be picked up by setting the
+    // originLevel of input cells to be a very large level (but within range
+    // of 1<< shift inside refinementDistanceData::wantedLevel)
+    forAll(cellsToRefine, i)
+    {
+        label cellI = cellsToRefine[i];
+
+        allCellInfo[cellI].originLevel() = sizeof(label)*8-2;
+        allCellInfo[cellI].origin() = cc[cellI];
+    }
+
+    // 2. Extend to 2:1. For non-cube cells the scalar distance does not work
+    // so make sure it at least provides 2:1.
+    PackedList<1> refineCell(mesh_.nCells(), 0);
+    forAll(allCellInfo, cellI)
+    {
+        label wanted = allCellInfo[cellI].wantedLevel(cc[cellI]);
+
+        if (wanted > cellLevel_[cellI]+1)
+        {
+            refineCell.set(cellI, 1);
+        }
+    }
+
+    while (true)
+    {
+        label nChanged = faceConsistentRefinement(true, refineCell);
+
+        reduce(nChanged, sumOp<label>());
+
+        if (debug)
+        {
+            Pout<< "polyRef::consistentSlowRefinement2 : Changed " << nChanged
+                << " refinement levels due to 2:1 conflicts."
+                << endl;
+        }
+
+        if (nChanged == 0)
+        {
+            break;
+        }
+    }
+
+    // 3. Convert back to labelList.
+    label nRefined = 0;
+
+    forAll(refineCell, cellI)
+    {
+        if (refineCell.get(cellI) == 1)
+        {
+            nRefined++;
+        }
+    }
+
+    labelList newCellsToRefine(nRefined);
+    nRefined = 0;
+
+    forAll(refineCell, cellI)
+    {
+        if (refineCell.get(cellI) == 1)
+        {
+            newCellsToRefine[nRefined++] = cellI;
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "polyRef::consistentSlowRefinement2 : From "
+            << cellsToRefine.size() << " to " << newCellsToRefine.size()
+            << " cells to refine." << endl;
+
+        // Check that newCellsToRefine obeys at least 2:1.
+
+        {
+            cellSet cellsIn
+            (
+                mesh_,
+                "cellsToRefineIn",
+                labelHashSet(cellsToRefine)
+            );
+            Pout<< "polyRef::consistentSlowRefinement2 : writing "
+                << cellsIn.size() << " to cellSet "
+                << cellsIn.objectPath() << endl;
+            cellsIn.write();
+        }
+        {
+            cellSet cellsOut
+            (
+                mesh_,
+                "cellsToRefineOut",
+                labelHashSet(newCellsToRefine)
+            );
+            Pout<< "polyRef::consistentSlowRefinement2 : writing "
+                << cellsOut.size() << " to cellSet "
+                << cellsOut.objectPath() << endl;
+            cellsOut.write();
+        }
+
+        // Extend to 2:1
+        PackedList<1> refineCell(mesh_.nCells(), 0);
+        forAll(newCellsToRefine, i)
+        {
+            refineCell.set(newCellsToRefine[i], 1);
+        }
+        const PackedList<1> savedRefineCell(refineCell);
+
+        label nChanged = faceConsistentRefinement(true, refineCell);
+
+        {
+            cellSet cellsOut2
+            (
+                mesh_, "cellsToRefineOut2", newCellsToRefine.size()
+            );
+            forAll(refineCell, cellI)
+            {
+                if (refineCell.get(cellI) == 1)
+                {
+                    cellsOut2.insert(cellI);
+                }
+            }
+            Pout<< "polyRef::consistentSlowRefinement2 : writing "
+                << cellsOut2.size() << " to cellSet "
+                << cellsOut2.objectPath() << endl;
+            cellsOut2.write();
+        }
+
+        if (nChanged > 0)
+        {
+            forAll(refineCell, cellI)
+            {
+                if
+                (
+                    refineCell.get(cellI) == 1
+                 && savedRefineCell.get(cellI) == 0
+                )
+                {
+                    FatalErrorIn
+                    (
+                        "polyRef::consistentSlowRefinement2"
+                        "(const label, const labelList&, const labelList&)"
+                    )   << "Cell:" << cellI << " cc:"
+                        << mesh_.cellCentres()[cellI]
+                        << " was not marked for refinement but does not obey"
+                        << " 2:1 constraints."
+                        << abort(FatalError);
+                }
+            }
+        }
+    }
+
+    return newCellsToRefine;
+}
+
+
+// Top level driver to insert topo changes to do all refinement.
+Foam::labelListList Foam::polyRef::setRefinement
+(
+    const labelList& cellLabels,
+    directTopoChange& meshMod
+)
+{
+    if (debug)
+    {
+        Pout<< "polyRef::setRefinement :"
+            << " Checking initial mesh just to make sure" << endl;
+
+        checkMesh();
+        checkRefinementLevels(-1, labelList(0));
+    }
+
+    // Clear any saved point/cell data.
+    savedPointLevel_.clear();
+    savedCellLevel_.clear();
+
+
+    // New point/cell level. Copy of pointLevel for existing points.
+    dynamicLabelList newCellLevel(cellLevel_.size());
+    forAll(cellLevel_, cellI)
+    {
+        newCellLevel.append(cellLevel_[cellI]);
+    }
+    dynamicLabelList newPointLevel(pointLevel_.size());
+    forAll(pointLevel_, pointI)
+    {
+        newPointLevel.append(pointLevel_[pointI]);
+    }
+
+
+    if (debug)
+    {
+        Pout<< "polyRef::setRefinement :"
+            << " Allocating " << cellLabels.size() << " cell midpoints."
+            << endl;
+    }
+
+    // Get necessary mesh data
+    const faceList& meshFaces = mesh_.faces();
+    const cellList& meshCells = mesh_.cells();
+
+    // Mid point per refined cell.
+    // -1 : not refined
+    // >=0: label of mid point.
+    labelList cellMidPoint(mesh_.nCells(), -1);
+
+    forAll(cellLabels, i)
+    {
+        label cellI = cellLabels[i];
+
+        label anchorPointI = meshFaces[meshCells[cellI][0]][0];
+
+        cellMidPoint[cellI] = meshMod.setAction
+        (
+            polyAddPoint
+            (
+                mesh_.cellCentres()[cellI],     // point
+                anchorPointI,                   // master point
+                -1,                             // zone for point
+                true                            // supports a cell
+            )
+        );
+
+        newPointLevel(cellMidPoint[cellI]) = cellLevel_[cellI] + 1;
+    }
+
+
+    if (debug)
+    {
+        cellSet splitCells(mesh_, "splitCells", cellLabels.size());
+
+        forAll(cellMidPoint, cellI)
+        {
+            if (cellMidPoint[cellI] >= 0)
+            {
+                splitCells.insert(cellI);
+            }
+        }
+
+        Pout<< "polyRef::setRefinement : Dumping " << splitCells.size()
+            << " cells to split to cellSet " << splitCells.objectPath()
+            << endl;
+
+        splitCells.write();
+    }
+
+
+
+    // Split edges
+    // ~~~~~~~~~~~
+
+    if (debug)
+    {
+        Pout<< "polyRef::setRefinement :"
+            << " Allocating edge midpoints."
+            << endl;
+    }
+
+    // Unrefined edges are ones between cellLevel or lower points.
+    // If any cell using this edge gets split then the edge needs to be split.
+
+    // Get necessary mesh data
+    const labelListList& meshCellEdges = mesh_.cellEdges();
+    const edgeList& meshEdges = mesh_.edges();
+
+    // -1  : no need to split edge
+    // >=0 : label of introduced mid point
+    labelList edgeMidPoint(mesh_.nEdges(), -1);
+
+    // Note: Loop over cells to be refined or edges?
+    forAll(cellMidPoint, cellI)
+    {
+        if (cellMidPoint[cellI] >= 0)
+        {
+            const labelList& cEdges = meshCellEdges[cellI];
+
+            forAll(cEdges, i)
+            {
+                label edgeI = cEdges[i];
+
+                const edge& e = meshEdges[edgeI];
+
+                if
+                (
+                    pointLevel_[e[0]] <= cellLevel_[cellI]
+                 && pointLevel_[e[1]] <= cellLevel_[cellI]
+                )
+                {
+                    edgeMidPoint[edgeI] = 12345;    // mark need for splitting
+                }
+            }
+        }
+    }
+
+    // Synchronize edgeMidPoint across coupled patches. Take max so that
+    // any split takes precedence.
+    syncTools::syncEdgeList
+    (
+        mesh_,
+        edgeMidPoint,
+        maxEqOp<label>(),
+        labelMin,
+        false               // no separation
+    );
+
+
+    // Introduce edge points
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    // Get necessary mesh data
+    const pointField& meshPoints = mesh_.points();
+
+    {
+        // Phase 1: calculate midpoints and sync.
+        // This needs doing for if people do not write binary and we slowly
+        // get differences.
+
+        pointField edgeMids(mesh_.nEdges(), point(-GREAT, -GREAT, -GREAT));
+
+        forAll(edgeMidPoint, edgeI)
+        {
+            if (edgeMidPoint[edgeI] >= 0)
+            {
+                // Edge marked to be split.
+                edgeMids[edgeI] = meshEdges[edgeI].centre(meshPoints);
+            }
+        }
+
+        syncTools::syncEdgeList
+        (
+            mesh_,
+            edgeMids,
+            maxEqOp<vector>(),
+            point(-GREAT, -GREAT, -GREAT),
+            true               // apply separation
+        );
+
+
+        // Phase 2: introduce points at the synced locations.
+        forAll(edgeMidPoint, edgeI)
+        {
+            if (edgeMidPoint[edgeI] >= 0)
+            {
+                // Edge marked to be split. Replace edgeMidPoint with actual
+                // point label.
+
+                const edge& e = meshEdges[edgeI];
+
+                edgeMidPoint[edgeI] = meshMod.setAction
+                (
+                    polyAddPoint
+                    (
+                        edgeMids[edgeI],            // point
+                        e[0],                       // master point
+                        -1,                         // zone for point
+                        true                        // supports a cell
+                    )
+                );
+
+                // Update level: take maximum of the two points in the original
+                // edge and increment
+                newPointLevel(edgeMidPoint[edgeI]) =
+                    max
+                    (
+                        pointLevel_[e[0]],
+                        pointLevel_[e[1]]
+                    )
+                  + 1;
+            }
+        }
+    }
+
+    if (debug)
+    {
+        OFstream str(mesh_.time().path()/"edgeMidPoint.obj");
+
+        forAll(edgeMidPoint, edgeI)
+        {
+            if (edgeMidPoint[edgeI] >= 0)
+            {
+                const edge& e = meshEdges[edgeI];
+
+                meshTools::writeOBJ(str, e.centre(meshPoints));
+            }
+        }
+
+        Pout<< "polyRef::setRefinement :"
+            << " Dumping edge centres to split to file " << str.name() << endl;
+    }
+
+
+    // Calculate face level
+    // ~~~~~~~~~~~~~~~~~~~~
+    // (after splitting)
+
+    if (debug)
+    {
+        Pout<< "polyRef::setRefinement :"
+            << " Allocating face midpoints."
+            << endl;
+    }
+
+    // Face anchor level. There are guaranteed 4 points with level
+    // <= anchorLevel. These are the corner points.
+    labelList faceAnchorLevel(mesh_.nFaces());
+
+    for (label faceI = 0; faceI < mesh_.nFaces(); faceI++)
+    {
+        faceAnchorLevel[faceI] = getAnchorLevel(faceI);
+    }
+
+    // -1  : no need to split face
+    // >=0 : label of introduced mid point
+    labelList faceMidPoint(mesh_.nFaces(), -1);
+
+    // Get necessary mesh data
+    const labelList& meshFaceOwner = mesh_.faceOwner();
+    const labelList& meshFaceNeighbour = mesh_.faceNeighbour();
+
+    // Internal faces: look at cells on both sides. Uniquely determined since
+    // face itself guaranteed to be same level as most refined neighbour.
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        // Note: no need to check whether the face has valid anchor level since
+        // all faces can be split
+
+        const label own = meshFaceOwner[faceI];
+        const label ownLevel = cellLevel_[own];
+        const label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
+
+        const label nei = meshFaceNeighbour[faceI];
+        const label neiLevel = cellLevel_[nei];
+        const label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0);
+
+        if
+        (
+            newOwnLevel > faceAnchorLevel[faceI]
+         || newNeiLevel > faceAnchorLevel[faceI]
+        )
+        {
+            faceMidPoint[faceI] = 12345;    // mark to be split
+        }
+    }
+
+    // Coupled patches handled like internal faces except now all information
+    // from neighbour comes from across processor.
+    // Boundary faces are more complicated since the boundary face can
+    // be more refined than its owner (or neighbour for coupled patches)
+    // (does not happen if refining/unrefining only, but does e.g. when
+    // refinining and subsetting)
+
+    {
+        labelList newNeiLevel(mesh_.nFaces() - mesh_.nInternalFaces());
+
+        forAll(newNeiLevel, i)
+        {
+            const label own = meshFaceOwner[i + mesh_.nInternalFaces()];
+            const label ownLevel = cellLevel_[own];
+            const label newOwnLevel =
+                ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
+
+            newNeiLevel[i] = newOwnLevel;
+        }
+
+        // Swap the list which now contains data from the other side
+        syncTools::swapBoundaryFaceList(mesh_, newNeiLevel, false);
+
+        forAll(newNeiLevel, i)
+        {
+            const label faceI = i + mesh_.nInternalFaces();
+
+            // Note: no need to check whether the face has valid anchor level
+            // since all faces can be split
+            const label own = meshFaceOwner[faceI];
+            const label ownLevel = cellLevel_[own];
+            const label newOwnLevel =
+                ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
+
+            if
+            (
+                newOwnLevel > faceAnchorLevel[faceI]
+             || newNeiLevel[i] > faceAnchorLevel[faceI]
+            )
+            {
+                faceMidPoint[faceI] = 12345;    // mark to be split
+            }
+        }
+    }
+
+
+    // Synchronize faceMidPoint across coupled patches. (logical or)
+    syncTools::syncFaceList
+    (
+        mesh_,
+        faceMidPoint,
+        maxEqOp<label>(),
+        false
+    );
+
+
+
+    // Introduce face points
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    // Get face centres
+    const vectorField& meshFaceCentres = mesh_.faceCentres();
+
+    {
+        // Phase 1: determine mid points and sync. See comment for edgeMids
+        // above
+        pointField bFaceMids
+        (
+            mesh_.nFaces() - mesh_.nInternalFaces(),
+            point(-GREAT, -GREAT, -GREAT)
+        );
+
+        forAll(bFaceMids, i)
+        {
+            const label faceI = i + mesh_.nInternalFaces();
+
+            if (faceMidPoint[faceI] >= 0)
+            {
+                bFaceMids[i] = meshFaceCentres[faceI];
+            }
+        }
+
+        syncTools::syncBoundaryFaceList
+        (
+            mesh_,
+            bFaceMids,
+            maxEqOp<vector>(),
+            true               // apply separation
+        );
+
+        forAll(faceMidPoint, faceI)
+        {
+            if (faceMidPoint[faceI] >= 0)
+            {
+                // Face marked to be split. Replace faceMidPoint with actual
+                // point label.
+
+                const face& f = meshFaces[faceI];
+
+                faceMidPoint[faceI] = meshMod.setAction
+                (
+                    polyAddPoint
+                    (
+                        (
+                            faceI < mesh_.nInternalFaces()
+                          ? meshFaceCentres[faceI]
+                          : bFaceMids[faceI - mesh_.nInternalFaces()]
+                        ),                          // point
+                        f[0],                       // master point
+                        -1,                         // zone for point
+                        true                        // supports a cell
+                    )
+                );
+
+                // Determine the level of the corner points and midpoint will
+                // be one higher.
+                newPointLevel(faceMidPoint[faceI]) = faceAnchorLevel[faceI] + 1;
+            }
+        }
+    }
+
+    if (debug)
+    {
+        faceSet splitFaces(mesh_, "splitFaces", cellLabels.size());
+
+        forAll(faceMidPoint, faceI)
+        {
+            if (faceMidPoint[faceI] >= 0)
+            {
+                splitFaces.insert(faceI);
+            }
+        }
+
+        Pout<< "polyRef::setRefinement : Dumping " << splitFaces.size()
+            << " faces to split to faceSet " << splitFaces.objectPath() << endl;
+
+        splitFaces.write();
+    }
+
+
+    // Information complete
+    // ~~~~~~~~~~~~~~~~~~~~
+    // At this point we have all the information we need. We should no
+    // longer reference the cellLabels to refine. All the information is:
+    // - cellMidPoint >= 0 : cell needs to be split
+    // - faceMidPoint >= 0 : face needs to be split
+    // - edgeMidPoint >= 0 : edge needs to be split
+
+
+
+    // Get the corner/anchor points
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    if (debug)
+    {
+        Pout<< "polyRef::setRefinement :"
+            << " Finding cell anchorPoints (8 per cell)"
+            << endl;
+    }
+
+    // Get anchor points for each cell: points that have the same or lower
+    // refinement level as the cell.
+    List<dynamicLabelList> cellAnchorPointsDynamic(mesh_.nCells());
+
+    {
+        // Loop through all cells
+        forAll(cellMidPoint, cellI)
+        {
+            if (cellMidPoint[cellI] >= 0)
+            {
+                // Set capacity for this cell (that's supposed to be cut) to 16
+                // anchor points to prevent excessive resizing
+                cellAnchorPointsDynamic[cellI].setCapacity(16);
+            }
+        }
+
+        // Loop through all points
+        forAll(pointLevel_, pointI)
+        {
+            const labelList& pCells = mesh_.pointCells()[pointI];
+
+            // Loop through all cells sharing this point
+            forAll(pCells, pCellI)
+            {
+                const label cellI = pCells[pCellI];
+
+                if
+                (
+                    cellMidPoint[cellI] >= 0
+                 && pointLevel_[pointI] <= cellLevel_[cellI]
+                )
+                {
+                    cellAnchorPointsDynamic[cellI].append(pointI);
+                }
+            }
+        }
+
+        // Loop through all cells and check whether at least 4 anchor points
+        // have been found (minimum requirement for a tet cell)
+        forAll(cellMidPoint, cellI)
+        {
+            if (cellMidPoint[cellI] >= 0)
+            {
+                if (cellAnchorPointsDynamic[cellI].size() <= 3)
+                {
+                    const labelList cPoints(cellPoints(cellI));
+
+                    FatalErrorIn
+                    (
+                        "polyRef::setRefinement(const labelList&"
+                        ", directTopoChange&)"
+                    )   << "cell " << cellI
+                        << " of level " << cellLevel_[cellI]
+                        << " does not seem to have enough points of "
+                        << " lower level" << endl
+                        << "cellPoints:" << cPoints << endl
+                        << "pointLevels:"
+                        << IndirectList<label>(pointLevel_, cPoints)() << endl
+                        << abort(FatalError);
+                }
+            }
+        }
+    }
+
+    // Collect cellAnchorPoints into a List<labelList> instead of
+    // List<dynamicList>
+    labelListList cellAnchorPoints(mesh_.nCells());
+
+    forAll(cellAnchorPointsDynamic, cellI)
+    {
+        // Tranfer the dynamic list for each cell into an ordinary list
+        cellAnchorPoints[cellI].transfer(cellAnchorPointsDynamic[cellI]);
+    }
+
+    // Add the cells
+    // ~~~~~~~~~~~~~
+
+    if (debug)
+    {
+        Pout<< "polyRef::setRefinement :"
+            << " Adding cells (1 per anchorPoint)"
+            << endl;
+    }
+
+    // We should have exactly n new cells per each split cells, where n is the
+    // number of anchor points in a cell
+    labelListList cellAddedCells(mesh_.nCells());
+
+    forAll(cellAnchorPoints, cellI)
+    {
+        // Check whether this is a split cell
+        if (cellMidPoint[cellI] >= 0)
+        {
+            // Get cell anchors
+            const labelList& cAnchors = cellAnchorPoints[cellI];
+
+            // Set the total number of added cells
+            labelList& cAdded = cellAddedCells[cellI];
+            cAdded.setSize(cAnchors.size());
+
+            // Original cell at 0
+            cAdded[0] = cellI;
+
+            // Update cell level
+            newCellLevel[cellI] = cellLevel_[cellI] + 1;
+
+            for (label i = 1; i < cAdded.size(); ++i)
+            {
+                cAdded[i] = meshMod.setAction
+                (
+                    polyAddCell
+                    (
+                        -1,                                 // master point
+                        -1,                                 // master edge
+                        -1,                                 // master face
+                        cellI,                              // master cell
+                        mesh_.cellZones().whichZone(cellI)  // zone for cell
+                    )
+                );
+
+                // Update level for this cell
+                newCellLevel(cAdded[i]) = cellLevel_[cellI] + 1;
+            }
+        }
+    }
+
+
+    // Faces
+    // ~~~~~
+    // 1. existing faces that get split (into n faces where n is the number of
+    //    points or edges)
+    // 2. existing faces that do not get split but only edges get split
+    // 3. existing faces that do not get split but get new owner/neighbour
+    // 4. new internal faces inside split cells.
+
+    if (debug)
+    {
+        Pout<< "polyRef::setRefinement :"
+            << " Marking faces to be handled"
+            << endl;
+    }
+
+    // Get all affected faces.
+    PackedList<1> affectedFace(mesh_.nFaces(), 0);
+
+    // Get mesh edges
+    const labelListList& meshEdgeFaces = mesh_.edgeFaces();
+
+    {
+        // All faces of a cell that's being split
+        forAll(cellMidPoint, cellI)
+        {
+            if (cellMidPoint[cellI] >= 0)
+            {
+                const cell& cFaces = meshCells[cellI];
+
+                forAll(cFaces, i)
+                {
+                    affectedFace.set(cFaces[i], 1);
+                }
+            }
+        }
+
+        // All faces that are being split
+        forAll(faceMidPoint, faceI)
+        {
+            if (faceMidPoint[faceI] >= 0)
+            {
+                affectedFace.set(faceI, 1);
+            }
+        }
+
+        // Both faces of an edge that are being split
+        forAll(edgeMidPoint, edgeI)
+        {
+            if (edgeMidPoint[edgeI] >= 0)
+            {
+                const labelList& eFaces = meshEdgeFaces[edgeI];
+
+                forAll(eFaces, i)
+                {
+                    affectedFace.set(eFaces[i], 1);
+                }
+            }
+        }
+    }
+
+
+    // 1. Faces that get split
+    // ~~~~~~~~~~~~~~~~~~~~~~~
+
+    if (debug)
+    {
+        Pout<< "polyRef::setRefinement : Splitting faces" << endl;
+    }
+
+    forAll(faceMidPoint, faceI)
+    {
+        if (faceMidPoint[faceI] >= 0 && affectedFace.get(faceI) == 1)
+        {
+            // Face needs to be split and hasn't yet been done in some way
+            // (affectedFace - is impossible since this is first change but
+            // just for completeness)
+
+            const face& f = meshFaces[faceI];
+
+            // Has original faceI been used (n - 1 faces added, original gets
+            // modified). n is the number of points/edges of a face
+            bool modifiedFace = false;
+
+            label anchorLevel = faceAnchorLevel[faceI];
+
+            // New face always has four points/edges
+            face newFace(4);
+
+            // Loop through all points of original face
+            forAll(f, fp)
+            {
+                const label pointI = f[fp];
+
+                if (pointLevel_[pointI] <= anchorLevel)
+                {
+                    // point is anchor. Start collecting face.
+
+                    // Create a dynamic storage for face vertices and append the
+                    // first (anchor) point
+                    dynamicLabelList faceVerts(4);
+                    faceVerts.append(pointI);
+
+                    // Walk forward to mid point.
+                    // - if next is +2 midpoint is +1
+                    // - if next is +1 it is midpoint
+                    // - if next is +0 there has to be edgeMidPoint
+
+                    // Appends all points from this point to face mid point
+                    walkFaceToMid
+                    (
+                        edgeMidPoint,
+                        anchorLevel,
+                        faceI,
+                        fp,
+                        faceVerts
+                    );
+
+                    // Append face mid point
+                    faceVerts.append(faceMidPoint[faceI]);
+
+                    // Append all points from face mid point to starting point
+                    walkFaceFromMid
+                    (
+                        edgeMidPoint,
+                        anchorLevel,
+                        faceI,
+                        fp,
+                        faceVerts
+                    );
+
+                    // Transfer dynamic list to a face (ordinary list)
+                    newFace.transfer(faceVerts.shrink());
+                    faceVerts.clear();
+
+                    if (debug)
+                    {
+                        Pout<< "Split face: " << faceI << ", verts: " << f
+                            << ", into quad: " << newFace << endl;
+                    }
+
+                    // Get new owner/neighbour
+                    label own, nei;
+                    getFaceNeighbours
+                    (
+                        cellAnchorPoints,
+                        cellAddedCells,
+                        faceI,
+                        pointI,          // Anchor point
+
+                        own,
+                        nei
+                    );
+
+
+                    if (debug)
+                    {
+                        // Get mesh cell centres
+                        const vectorField& meshCellCentres =
+                            mesh_.cellCentres();
+
+                        if (mesh_.isInternalFace(faceI))
+                        {
+                            label oldOwn = meshFaceOwner[faceI];
+                            label oldNei = meshFaceNeighbour[faceI];
+
+                            checkInternalOrientation
+                            (
+                                meshMod,
+                                oldOwn,
+                                faceI,
+                                meshCellCentres[oldOwn],
+                                meshCellCentres[oldNei],
+                                newFace
+                            );
+                        }
+                        else
+                        {
+                            label oldOwn = meshFaceOwner[faceI];
+
+                            checkBoundaryOrientation
+                            (
+                                meshMod,
+                                oldOwn,
+                                faceI,
+                                meshCellCentres[oldOwn],
+                                meshFaceCentres[faceI],
+                                newFace
+                            );
+                        }
+                    } // End debug
+
+
+                    if (!modifiedFace)
+                    {
+                        modifiedFace = true;
+
+                        modFace(meshMod, faceI, newFace, own, nei);
+                    }
+                    else
+                    {
+                        addFace(meshMod, faceI, newFace, own, nei);
+                    }
+                } // End point anchor chech
+            } // End for all points
+
+            // Mark face as having been handled
+            affectedFace.set(faceI, 0);
+        }
+    }
+
+
+    // 2. faces that do not get split but use edges that get split
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    if (debug)
+    {
+        Pout<< "polyRef::setRefinement :"
+            << " Adding edge splits to unsplit faces"
+            << endl;
+    }
+
+    // Get face edges
+    const labelListList& meshFaceEdges = mesh_.faceEdges();
+
+    forAll(edgeMidPoint, edgeI)
+    {
+        if (edgeMidPoint[edgeI] >= 0)
+        {
+            // Split edge. Check that face not already handled above.
+
+            const labelList& eFaces = meshEdgeFaces[edgeI];
+
+            forAll(eFaces, i)
+            {
+                const label faceI = eFaces[i];
+
+                if (faceMidPoint[faceI] < 0 && affectedFace.get(faceI) == 1)
+                {
+                    // Unsplit face that hasn't been handled. Add edge splits to
+                    // face.
+
+                    const face& f = meshFaces[faceI];
+                    const labelList& fEdges = meshFaceEdges[faceI];
+
+                    dynamicLabelList newFaceVerts(f.size());
+
+                    // Append all original points and all edge mid points
+                    forAll(f, fp)
+                    {
+                        newFaceVerts.append(f[fp]);
+
+                        const label edgeI = fEdges[fp];
+
+                        if (edgeMidPoint[edgeI] >= 0)
+                        {
+                            newFaceVerts.append(edgeMidPoint[edgeI]);
+                        }
+                    }
+
+                    face newFace;
+                    newFace.transfer(newFaceVerts.shrink());
+
+
+                    // The point with the lowest level should be an anchor
+                    // point of the neighbouring cells.
+                    const label anchorFp = findMinLevel(f);
+
+                    label own, nei;
+                    getFaceNeighbours
+                    (
+                        cellAnchorPoints,
+                        cellAddedCells,
+                        faceI,
+                        f[anchorFp],          // Anchor point
+
+                        own,
+                        nei
+                    );
+
+
+                    if (debug)
+                    {
+                        // Get mesh cell centres
+                        const vectorField& meshCellCentres =
+                            mesh_.cellCentres();
+
+                        if (mesh_.isInternalFace(faceI))
+                        {
+                            label oldOwn = meshFaceOwner[faceI];
+                            label oldNei = meshFaceNeighbour[faceI];
+
+                            checkInternalOrientation
+                            (
+                                meshMod,
+                                oldOwn,
+                                faceI,
+                                meshCellCentres[oldOwn],
+                                meshCellCentres[oldNei],
+                                newFace
+                            );
+                        }
+                        else
+                        {
+                            label oldOwn = meshFaceOwner[faceI];
+
+                            checkBoundaryOrientation
+                            (
+                                meshMod,
+                                oldOwn,
+                                faceI,
+                                meshCellCentres[oldOwn],
+                                meshFaceCentres[faceI],
+                                newFace
+                            );
+                        }
+                    } // End debug
+
+                    modFace(meshMod, faceI, newFace, own, nei);
+
+                    // Mark face as having been handled
+                    affectedFace.set(faceI, 0);
+                } // End if unsplit, unhandled face
+            } // End for all edge faces
+        } // End if edge has been cut
+    } // End for all edges
+
+
+    // 3. faces that do not get split but whose owner/neighbour change
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    if (debug)
+    {
+        Pout<< "polyRef::setRefinement :"
+            << " Changing owner/neighbour for otherwise unaffected faces"
+            << endl;
+    }
+
+    forAll(affectedFace, faceI)
+    {
+        if (affectedFace.get(faceI) == 1)
+        {
+            const face& f = meshFaces[faceI];
+
+            // The point with the lowest level should be an anchor
+            // point of the neighbouring cells.
+            label anchorFp = findMinLevel(f);
+
+            label own, nei;
+            getFaceNeighbours
+            (
+                cellAnchorPoints,
+                cellAddedCells,
+                faceI,
+                f[anchorFp],          // Anchor point
+
+                own,
+                nei
+            );
+
+            modFace(meshMod, faceI, f, own, nei);
+
+            // Mark face as having been handled
+            affectedFace.set(faceI, 0);
+        }
+    }
+
+
+    // 4. new internal faces inside split cells.
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+    // This is the hard one. We have to find the splitting points between
+    // the anchor points. But the edges between the anchor points might have
+    // been split (into two,three or four edges).
+
+    if (debug)
+    {
+        Pout<< "polyRef::setRefinement :"
+            << " Create new internal faces for split cells"
+            << endl;
+    }
+
+    forAll(cellMidPoint, cellI)
+    {
+        if (cellMidPoint[cellI] >= 0)
+        {
+            createInternalFaces
+            (
+                cellAnchorPoints,
+                cellAddedCells,
+                cellMidPoint,
+                faceMidPoint,
+                faceAnchorLevel,
+                edgeMidPoint,
+                cellI,
+                meshMod
+            );
+        }
+    }
+
+    // Extend pointLevels and cellLevels for the new cells. Could also be done
+    // in updateMesh but saves passing cellAddedCells out of this routine.
+
+    // Check
+    if (debug)
+    {
+        label minPointI = labelMax;
+        label maxPointI = labelMin;
+
+        forAll(cellMidPoint, cellI)
+        {
+            if (cellMidPoint[cellI] >= 0)
+            {
+                minPointI = min(minPointI, cellMidPoint[cellI]);
+                maxPointI = max(maxPointI, cellMidPoint[cellI]);
+            }
+        }
+        forAll(faceMidPoint, faceI)
+        {
+            if (faceMidPoint[faceI] >= 0)
+            {
+                minPointI = min(minPointI, faceMidPoint[faceI]);
+                maxPointI = max(maxPointI, faceMidPoint[faceI]);
+            }
+        }
+        forAll(edgeMidPoint, edgeI)
+        {
+            if (edgeMidPoint[edgeI] >= 0)
+            {
+                minPointI = min(minPointI, edgeMidPoint[edgeI]);
+                maxPointI = max(maxPointI, edgeMidPoint[edgeI]);
+            }
+        }
+
+        if (minPointI != labelMax && minPointI != mesh_.nPoints())
+        {
+            FatalErrorIn("polyRef::setRefinement(..)")
+                << "Added point labels not consecutive to existing mesh points."
+                << nl
+                << "mesh_.nPoints():" << mesh_.nPoints()
+                << " minPointI:" << minPointI
+                << " maxPointI:" << maxPointI
+                << abort(FatalError);
+        }
+    }
+
+    pointLevel_.transfer(newPointLevel.shrink());
+    newPointLevel.clear();
+    cellLevel_.transfer(newCellLevel.shrink());
+    newCellLevel.clear();
+
+    // Mark files as changed
+    setInstance(mesh_.facesInstance());
+
+
+    // Update the live split cells tree.
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // New unrefinement structure
+    if (history_.active())
+    {
+        if (debug)
+        {
+            Pout<< "polyRef::setRefinement :"
+                << " Updating refinement history to " << cellLevel_.size()
+                << " cells" << endl;
+        }
+
+        // Extend refinement history for new cells
+        history_.resize(cellLevel_.size());
+
+        forAll(cellAddedCells, cellI)
+        {
+            const labelList& addedCells = cellAddedCells[cellI];
+
+            if (addedCells.size() > 0)
+            {
+                // Cell was split.
+                history_.storeSplit(cellI, addedCells);
+            }
+        }
+    }
+
+    // Compact cellAddedCells.
+
+    labelListList refinedCells(cellLabels.size());
+
+    forAll(cellLabels, i)
+    {
+        label cellI = cellLabels[i];
+
+        refinedCells[i].transfer(cellAddedCells[cellI]);
+    }
+
+    return refinedCells;
+}
+
+
+void Foam::polyRef::storeData
+(
+    const labelList& pointsToStore,
+    const labelList& facesToStore,
+    const labelList& cellsToStore
+)
+{
+    savedPointLevel_.resize(2*pointsToStore.size());
+    forAll(pointsToStore, i)
+    {
+        label pointI = pointsToStore[i];
+        savedPointLevel_.insert(pointI, pointLevel_[pointI]);
+    }
+
+    savedCellLevel_.resize(2*cellsToStore.size());
+    forAll(cellsToStore, i)
+    {
+        label cellI = cellsToStore[i];
+        savedCellLevel_.insert(cellI, cellLevel_[cellI]);
+    }
+}
+
+
+// Gets called after the mesh change. setRefinement will already have made
+// sure the pointLevel_ and cellLevel_ are the size of the new mesh so we
+// only need to account for reordering.
+void Foam::polyRef::updateMesh(const mapPolyMesh& map)
+{
+    Map<label> dummyMap(0);
+
+    updateMesh(map, dummyMap, dummyMap, dummyMap);
+}
+
+
+// Gets called after the mesh change. setRefinement will already have made
+// sure the pointLevel_ and cellLevel_ are the size of the new mesh so we
+// only need to account for reordering.
+void Foam::polyRef::updateMesh
+(
+    const mapPolyMesh& map,
+    const Map<label>& pointsToRestore,
+    const Map<label>& facesToRestore,
+    const Map<label>& cellsToRestore
+)
+{
+    // Update celllevel
+    if (debug)
+    {
+        Pout<< "polyRef::updateMesh :"
+            << " Updating various lists"
+            << endl;
+    }
+
+    {
+        const labelList& reverseCellMap = map.reverseCellMap();
+
+        if (debug)
+        {
+            Pout<< "polyRef::updateMesh :"
+                << " reverseCellMap:" << map.reverseCellMap().size()
+                << " cellMap:" << map.cellMap().size()
+                << " nCells:" << mesh_.nCells()
+                << " nOldCells:" << map.nOldCells()
+                << " cellLevel_:" << cellLevel_.size()
+                << " reversePointMap:" << map.reversePointMap().size()
+                << " pointMap:" << map.pointMap().size()
+                << " nPoints:" << mesh_.nPoints()
+                << " nOldPoints:" << map.nOldPoints()
+                << " pointLevel_:" << pointLevel_.size()
+                << endl;
+        }
+
+        if (reverseCellMap.size() == cellLevel_.size())
+        {
+            // Assume it is after polyRef that this routine is called.
+            // Just account for reordering. We cannot use cellMap since
+            // then cells created from cells would get cellLevel_ of
+            // cell they were created from.
+            reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_);
+        }
+        else
+        {
+            // Map data
+            const labelList& cellMap = map.cellMap();
+
+            labelList newCellLevel(cellMap.size());
+            forAll(cellMap, newCellI)
+            {
+                label oldCellI = cellMap[newCellI];
+
+                if (oldCellI == -1)
+                {
+                    //FatalErrorIn("polyRef::updateMesh(const mapPolyMesh&)")
+                    //    << "Problem : cell " << newCellI
+                    //    << " at " << mesh_.cellCentres()[newCellI]
+                    //    << " does not originate from another cell"
+                    //    << " (i.e. is inflated)." << nl
+                    //    << "Hence we cannot determine the new cellLevel"
+                    //    << " for it." << abort(FatalError);
+                    newCellLevel[newCellI] = -1;
+                }
+                else
+                {
+                    newCellLevel[newCellI] = cellLevel_[oldCellI];
+                }
+            }
+            cellLevel_.transfer(newCellLevel);
+        }
+
+        // See if any cells to restore. This will be for some new cells
+        // the corresponding old cell.
+        forAllConstIter(Map<label>, cellsToRestore, iter)
+        {
+            label newCellI = iter.key();
+            label storedCellI = iter();
+
+            Map<label>::iterator fnd = savedCellLevel_.find(storedCellI);
+
+            if (fnd == savedCellLevel_.end())
+            {
+                FatalErrorIn("polyRef::updateMesh(const mapPolyMesh&)")
+                    << "Problem : trying to restore old value for new cell "
+                    << newCellI << " but cannot find old cell " << storedCellI
+                    << " in map of stored values " << savedCellLevel_
+                    << abort(FatalError);
+            }
+            cellLevel_[newCellI] = fnd();
+        }
+
+        //if (findIndex(cellLevel_, -1) != -1)
+        //{
+        //    WarningIn("polyRef::updateMesh(const mapPolyMesh&)")
+        //        << "Problem : "
+        //        << "cellLevel_ contains illegal value -1 after mapping
+        //        << " at cell " << findIndex(cellLevel_, -1) << endl
+        //        << "This means that another program has inflated cells"
+        //        << " (created cells out-of-nothing) and hence we don't know"
+        //        << " their cell level. Continuing with illegal value."
+        //        << abort(FatalError);
+        //}
+    }
+
+
+    // Update pointlevel
+    {
+        const labelList& reversePointMap = map.reversePointMap();
+
+        if (reversePointMap.size() == pointLevel_.size())
+        {
+            // Assume it is after polyRef that this routine is called.
+            reorder(reversePointMap, mesh_.nPoints(), -1,  pointLevel_);
+        }
+        else
+        {
+            // Map data
+            const labelList& pointMap = map.pointMap();
+
+            labelList newPointLevel(pointMap.size());
+
+            forAll(pointMap, newPointI)
+            {
+                label oldPointI = pointMap[newPointI];
+
+                if (oldPointI == -1)
+                {
+                    //FatalErrorIn("polyRef::updateMesh(const mapPolyMesh&)")
+                    //    << "Problem : point " << newPointI
+                    //    << " at " << mesh_.points()[newPointI]
+                    //    << " does not originate from another point"
+                    //    << " (i.e. is inflated)." << nl
+                    //    << "Hence we cannot determine the new pointLevel"
+                    //    << " for it." << abort(FatalError);
+                    newPointLevel[newPointI] = -1;
+                }
+                else
+                {
+                    newPointLevel[newPointI] = pointLevel_[oldPointI];
+                }
+            }
+            pointLevel_.transfer(newPointLevel);
+        }
+
+        // See if any points to restore. This will be for some new points
+        // the corresponding old point (the one from the call to storeData)
+        forAllConstIter(Map<label>, pointsToRestore, iter)
+        {
+            label newPointI = iter.key();
+            label storedPointI = iter();
+
+            Map<label>::iterator fnd = savedPointLevel_.find(storedPointI);
+
+            if (fnd == savedPointLevel_.end())
+            {
+                FatalErrorIn("polyRef::updateMesh(const mapPolyMesh&)")
+                    << "Problem : trying to restore old value for new point "
+                    << newPointI << " but cannot find old point "
+                    << storedPointI
+                    << " in map of stored values " << savedPointLevel_
+                    << abort(FatalError);
+            }
+            pointLevel_[newPointI] = fnd();
+        }
+
+        //if (findIndex(pointLevel_, -1) != -1)
+        //{
+        //    WarningIn("polyRef::updateMesh(const mapPolyMesh&)")
+        //        << "Problem : "
+        //        << "pointLevel_ contains illegal value -1 after mapping"
+        //        << " at point" << findIndex(pointLevel_, -1) << endl
+        //        << "This means that another program has inflated points"
+        //        << " (created points out-of-nothing) and hence we don't know"
+        //        << " their point level. Continuing with illegal value."
+        //        //<< abort(FatalError);
+        //}
+    }
+
+    // Update refinement tree
+    if (history_.active())
+    {
+        history_.updateMesh(map);
+    }
+
+    // Mark files as changed
+    setInstance(mesh_.facesInstance());
+
+    // Update face removal engine
+    faceRemover_.updateMesh(map);
+}
+
+
+// Gets called after mesh subsetting. Maps are from new back to old.
+void Foam::polyRef::subset
+(
+    const labelList& pointMap,
+    const labelList& faceMap,
+    const labelList& cellMap
+)
+{
+    // Update celllevel
+    if (debug)
+    {
+        Pout<< "polyRef::subset :"
+            << " Updating various lists"
+            << endl;
+    }
+
+    if (history_.active())
+    {
+        WarningIn
+        (
+            "polyRef::subset(const labelList&, const labelList&"
+            ", const labelList&)"
+        )   << "Subsetting will not work in combination with unrefinement."
+            << nl
+            << "Proceed at your own risk." << endl;
+    }
+
+
+    // Update celllevel
+    {
+        labelList newCellLevel(cellMap.size());
+
+        forAll(cellMap, newCellI)
+        {
+            newCellLevel[newCellI] = cellLevel_[cellMap[newCellI]];
+        }
+
+        cellLevel_.transfer(newCellLevel);
+
+        if (findIndex(cellLevel_, -1) != -1)
+        {
+            FatalErrorIn("polyRef::subset(..)")
+                << "Problem : "
+                << "cellLevel_ contains illegal value -1 after mapping:"
+                << cellLevel_
+                << abort(FatalError);
+        }
+    }
+
+    // Update pointlevel
+    {
+        labelList newPointLevel(pointMap.size());
+
+        forAll(pointMap, newPointI)
+        {
+            newPointLevel[newPointI] = pointLevel_[pointMap[newPointI]];
+        }
+
+        pointLevel_.transfer(newPointLevel);
+
+        if (findIndex(pointLevel_, -1) != -1)
+        {
+            FatalErrorIn("polyRef::subset(..)")
+                << "Problem : "
+                << "pointLevel_ contains illegal value -1 after mapping:"
+                << pointLevel_
+                << abort(FatalError);
+        }
+    }
+
+    // Update refinement tree
+    if (history_.active())
+    {
+        history_.subset(pointMap, faceMap, cellMap);
+    }
+
+    // Mark files as changed
+    setInstance(mesh_.facesInstance());
+
+    // Nothing needs doing to faceRemover.
+    //faceRemover_.subset(pointMap, faceMap, cellMap);
+}
+
+
+// Gets called after the mesh distribution
+void Foam::polyRef::distribute(const mapDistributePolyMesh& map)
+{
+    if (debug)
+    {
+        Pout<< "polyRef::distribute :"
+            << " Updating various lists"
+            << endl;
+    }
+
+    // Update celllevel
+    map.distributeCellData(cellLevel_);
+    // Update pointlevel
+    map.distributePointData(pointLevel_);
+
+    // Update refinement tree
+    if (history_.active())
+    {
+        history_.distribute(map);
+    }
+
+    // Update face removal engine
+    faceRemover_.distribute(map);
+}
+
+
+void Foam::polyRef::checkMesh() const
+{
+    const boundBox& meshBb = mesh_.globalData().bb();
+    const scalar smallDim = 1E-6*mag(meshBb.max() - meshBb.min());
+
+    if (debug)
+    {
+        Pout<< "polyRef::checkMesh : Using matching tolerance smallDim:"
+            << smallDim << endl;
+    }
+
+    // Check owner on coupled faces.
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // There should be only one coupled face between two cells. Why? Since
+    // otherwise mesh redistribution might cause multiple faces between two
+    // cells
+    {
+        labelList nei(mesh_.nFaces()-mesh_.nInternalFaces());
+        forAll(nei, i)
+        {
+            nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
+        }
+
+        // Replace data on coupled patches with their neighbour ones.
+        syncTools::swapBoundaryFaceList(mesh_, nei, false);
+
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+        forAll(patches, patchI)
+        {
+            const polyPatch& pp = patches[patchI];
+
+            if (pp.coupled())
+            {
+                // Check how many faces between owner and neighbour. Should
+                // be only one.
+                HashTable<label, labelPair, labelPair::Hash<> >
+                    cellToFace(2*pp.size());
+
+                label faceI = pp.start();
+
+                forAll(pp, i)
+                {
+                    label own = mesh_.faceOwner()[faceI];
+                    label bFaceI = faceI-mesh_.nInternalFaces();
+
+                    if (!cellToFace.insert(labelPair(own, nei[bFaceI]), faceI))
+                    {
+                        FatalErrorIn("polyRef::checkMesh()")
+                            << "Faces do not seem to be correct across coupled"
+                            << " boundaries" << endl
+                            << "Coupled face " << faceI
+                            << " between owner " << own
+                            << " on patch " << pp.name()
+                            << " and coupled neighbour " << nei[bFaceI]
+                            << " has two faces connected to it:"
+                            << faceI << " and "
+                            << cellToFace[labelPair(own, nei[bFaceI])]
+                            << abort(FatalError);
+                    }
+
+                    faceI++;
+                }
+            }
+        }
+    }
+
+    // Check face areas.
+    // ~~~~~~~~~~~~~~~~~
+
+    // Get necessary mesh data
+    const faceList& meshFaces = mesh_.faces();
+
+    {
+        scalarField neiFaceAreas(mesh_.nFaces()-mesh_.nInternalFaces());
+        forAll(neiFaceAreas, i)
+        {
+            neiFaceAreas[i] = mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
+        }
+
+        // Replace data on coupled patches with their neighbour ones.
+        syncTools::swapBoundaryFaceList(mesh_, neiFaceAreas, false);
+
+        forAll(neiFaceAreas, i)
+        {
+            label faceI = i+mesh_.nInternalFaces();
+
+            const scalar magArea = mag(mesh_.faceAreas()[faceI]);
+
+            if (mag(magArea - neiFaceAreas[i]) > smallDim)
+            {
+                const face& f = meshFaces[faceI];
+                label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+
+                FatalErrorIn("polyRef::checkMesh()")
+                    << "Faces do not seem to be correct across coupled"
+                    << " boundaries" << endl
+                    << "Coupled face " << faceI
+                    << " on patch " << patchI
+                    << " " << mesh_.boundaryMesh()[patchI].name()
+                    << " coords:" << IndirectList<point>(mesh_.points(), f)()
+                    << " has face area:" << magArea
+                    << " (coupled) neighbour face area differs:"
+                    << neiFaceAreas[i]
+                    << " to within tolerance " << smallDim
+                    << abort(FatalError);
+            }
+        }
+    }
+
+
+    // Check number of points on faces.
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    {
+        labelList nVerts(mesh_.nFaces()-mesh_.nInternalFaces());
+
+        forAll(nVerts, i)
+        {
+            nVerts[i] = meshFaces[i+mesh_.nInternalFaces()].size();
+        }
+
+        // Replace data on coupled patches with their neighbour ones.
+        syncTools::swapBoundaryFaceList(mesh_, nVerts, false);
+
+        forAll(nVerts, i)
+        {
+            label faceI = i+mesh_.nInternalFaces();
+
+            const face& f = meshFaces[faceI];
+
+            if (f.size() != nVerts[i])
+            {
+                label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+
+                FatalErrorIn("polyRef::checkMesh()")
+                    << "Faces do not seem to be correct across coupled"
+                    << " boundaries" << endl
+                    << "Coupled face " << faceI
+                    << " on patch " << patchI
+                    << " " << mesh_.boundaryMesh()[patchI].name()
+                    << " coords:" << IndirectList<point>(mesh_.points(), f)()
+                    << " has size:" << f.size()
+                    << " (coupled) neighbour face has size:"
+                    << nVerts[i]
+                    << abort(FatalError);
+            }
+        }
+    }
+
+
+    // Check points of face
+    // ~~~~~~~~~~~~~~~~~~~~
+    {
+        // Anchor points.
+        pointField anchorPoints(mesh_.nFaces()-mesh_.nInternalFaces());
+
+        forAll(anchorPoints, i)
+        {
+            label faceI = i+mesh_.nInternalFaces();
+            const point& fc = mesh_.faceCentres()[faceI];
+            const face& f = meshFaces[faceI];
+            const vector anchorVec(mesh_.points()[f[0]] - fc);
+
+            anchorPoints[i] = anchorVec;
+        }
+
+        // Replace data on coupled patches with their neighbour ones. Apply
+        // rotation transformation (but not separation since is relative vector
+        // to point on same face.
+        syncTools::swapBoundaryFaceList(mesh_, anchorPoints, false);
+
+        forAll(anchorPoints, i)
+        {
+            label faceI = i+mesh_.nInternalFaces();
+            const point& fc = mesh_.faceCentres()[faceI];
+            const face& f = meshFaces[faceI];
+            const vector anchorVec(mesh_.points()[f[0]] - fc);
+
+            if (mag(anchorVec - anchorPoints[i]) > smallDim)
+            {
+                label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+
+                FatalErrorIn("polyRef::checkMesh()")
+                    << "Faces do not seem to be correct across coupled"
+                    << " boundaries" << endl
+                    << "Coupled face " << faceI
+                    << " on patch " << patchI
+                    << " " << mesh_.boundaryMesh()[patchI].name()
+                    << " coords:" << IndirectList<point>(mesh_.points(), f)()
+                    << " has anchor vector:" << anchorVec
+                    << " (coupled) neighbour face anchor vector differs:"
+                    << anchorPoints[i]
+                    << " to within tolerance " << smallDim
+                    << abort(FatalError);
+            }
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "polyRef::checkMesh : Returning" << endl;
+    }
+}
+
+
+void Foam::polyRef::checkRefinementLevels
+(
+    const label maxPointDiff,
+    const labelList& pointsToCheck
+) const
+{
+    if (debug)
+    {
+        Pout<< "polyRef::checkRefinementLevels :"
+            << " Checking 2:1 refinement level" << endl;
+    }
+
+    if
+    (
+        cellLevel_.size() != mesh_.nCells()
+     || pointLevel_.size() != mesh_.nPoints()
+    )
+    {
+        FatalErrorIn("polyRef::checkRefinementLevels(const label)")
+            << "cellLevel size should be number of cells"
+            << " and pointLevel size should be number of points."<< nl
+            << "cellLevel:" << cellLevel_.size()
+            << " mesh.nCells():" << mesh_.nCells() << nl
+            << "pointLevel:" << pointLevel_.size()
+            << " mesh.nPoints():" << mesh_.nPoints()
+            << abort(FatalError);
+    }
+
+
+    // Check 2:1 consistency.
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    {
+        // Internal faces.
+        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+        {
+            label own = mesh_.faceOwner()[faceI];
+            label nei = mesh_.faceNeighbour()[faceI];
+
+            if (mag(cellLevel_[own] - cellLevel_[nei]) > 1)
+            {
+                FatalErrorIn
+                (
+                    "polyRef::checkRefinementLevels(const label)"
+                )   << "Celllevel does not satisfy 2:1 constraint." << nl
+                    << "On face " << faceI << " owner cell " << own
+                    << " has refinement " << cellLevel_[own]
+                    << " neighbour cell " << nei << " has refinement "
+                    << cellLevel_[nei]
+                    << abort(FatalError);
+            }
+        }
+
+        // Coupled faces. Get neighbouring value
+        labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+
+        forAll(neiLevel, i)
+        {
+            label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
+
+            neiLevel[i] = cellLevel_[own];
+        }
+
+        // No separation
+        syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
+
+        forAll(neiLevel, i)
+        {
+            label faceI = i+mesh_.nInternalFaces();
+
+            label own = mesh_.faceOwner()[faceI];
+
+            if (mag(cellLevel_[own] - neiLevel[i]) > 1)
+            {
+                label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+
+                FatalErrorIn
+                (
+                    "polyRef::checkRefinementLevels(const label)"
+                )   << "Celllevel does not satisfy 2:1 constraint."
+                    << " On coupled face " << faceI
+                    << " on patch " << patchI << " "
+                    << mesh_.boundaryMesh()[patchI].name()
+                    << " owner cell " << own << " has refinement "
+                    << cellLevel_[own]
+                    << " (coupled) neighbour cell has refinement "
+                    << neiLevel[i]
+                    << abort(FatalError);
+            }
+        }
+    }
+
+
+    // Check pointLevel is synchronized
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    {
+        labelList syncPointLevel(pointLevel_);
+
+        // Get min level
+        syncTools::syncPointList
+        (
+            mesh_,
+            syncPointLevel,
+            minEqOp<label>(),
+            labelMax,
+            false               // no separation
+        );
+
+
+        forAll(syncPointLevel, pointI)
+        {
+            if (pointLevel_[pointI] != syncPointLevel[pointI])
+            {
+                FatalErrorIn
+                (
+                    "polyRef::checkRefinementLevels(const label)"
+                )   << "PointLevel is not consistent across coupled patches."
+                    << endl
+                    << "point:" << pointI << " coord:" << mesh_.points()[pointI]
+                    << " has level " << pointLevel_[pointI]
+                    << " whereas the coupled point has level "
+                    << syncPointLevel[pointI]
+                    << abort(FatalError);
+            }
+        }
+    }
+
+
+    // Check 2:1 across points (instead of faces)
+    if (maxPointDiff != -1)
+    {
+        const labelListList& pointCells = mesh_.pointCells();
+
+        // Determine per point the max cell level.
+        labelList maxPointLevel(mesh_.nPoints(), 0);
+
+        forAll(pointCells, pointI)
+        {
+            const labelList& pCells = pointCells[pointI];
+
+            label& pLevel = maxPointLevel[pointI];
+
+            forAll(pCells, i)
+            {
+                pLevel = max(pLevel, cellLevel_[pCells[i]]);
+            }
+        }
+
+        // Sync maxPointLevel to neighbour
+        syncTools::syncPointList
+        (
+            mesh_,
+            maxPointLevel,
+            maxEqOp<label>(),
+            labelMin,           // null value
+            false               // no separation
+        );
+
+        // Check 2:1 across boundary points
+        forAll(pointsToCheck, i)
+        {
+            label pointI = pointsToCheck[i];
+
+            const labelList& pCells = pointCells[pointI];
+
+            forAll(pCells, i)
+            {
+                label cellI = pCells[i];
+
+                if
+                (
+                    mag(cellLevel_[cellI]-maxPointLevel[pointI])
+                  > maxPointDiff
+                )
+                {
+                    FatalErrorIn
+                    (
+                        "polyRef::checkRefinementLevels(const label)"
+                    )   << "Too big a difference between"
+                        << " point-connected cells." << nl
+                        << "cell:" << cellI
+                        << " cellLevel:" << cellLevel_[cellI]
+                        << " uses point:" << pointI
+                        << " coord:" << mesh_.points()[pointI]
+                        << " which is also used by a cell with level:"
+                        << maxPointLevel[pointI]
+                        << abort(FatalError);
+                }
+            }
+        }
+    }
+}
+
+
+
+//
+// Unrefinement
+// ~~~~~~~~~~~~
+//
+
+
+Foam::labelList Foam::polyRef::getSplitPoints() const
+{
+    if (debug)
+    {
+        checkRefinementLevels(-1, labelList(0));
+    }
+
+    if (debug)
+    {
+        Pout<< "polyRef::getSplitPoints :"
+            << " Calculating unrefineable points" << endl;
+    }
+
+
+    if (!history_.active())
+    {
+        FatalErrorIn("polyRef::getSplitPoints()")
+            << "Only call if constructed with history capability"
+            << abort(FatalError);
+    }
+
+    // Master cell
+    // -1 undetermined
+    // -2 certainly not split point
+    // >= label of master cell
+    labelList splitMaster(mesh_.nPoints(), -1);
+    labelList splitMasterLevel(mesh_.nPoints(), 0);
+
+    // Unmark all with different master cells
+    const labelList& visibleCells = history_.visibleCells();
+
+    forAll(visibleCells, cellI)
+    {
+        const labelList cPoints(cellPoints(cellI));
+
+        if (visibleCells[cellI] != -1 && history_.parentIndex(cellI) >= 0)
+        {
+            label parentIndex = history_.parentIndex(cellI);
+
+            // Check same master.
+            forAll(cPoints, i)
+            {
+                label pointI = cPoints[i];
+
+                label masterCellI = splitMaster[pointI];
+
+                if (masterCellI == -1)
+                {
+                    // First time visit of point. Store parent cell and
+                    // level of the parent cell (with respect to cellI). This
+                    // is additional guarantee that we're referring to the
+                    // same master at the same refinement level.
+
+                    splitMaster[pointI] = parentIndex;
+                    splitMasterLevel[pointI] = cellLevel_[cellI] - 1;
+                }
+                else if (masterCellI == -2)
+                {
+                    // Already decided that point is not splitPoint
+                }
+                else if
+                (
+                    (masterCellI != parentIndex)
+                 || (splitMasterLevel[pointI] != cellLevel_[cellI] - 1)
+                )
+                {
+                    // Different masters so point is on two refinement
+                    // patterns
+                    splitMaster[pointI] = -2;
+                }
+            }
+        }
+        else
+        {
+            // Either not visible or is unrefined cell
+            forAll(cPoints, i)
+            {
+                label pointI = cPoints[i];
+
+                splitMaster[pointI] = -2;
+            }
+        }
+    }
+
+    // Unmark boundary faces
+    for
+    (
+        label faceI = mesh_.nInternalFaces();
+        faceI < mesh_.nFaces();
+        faceI++
+    )
+    {
+        const face& f = mesh_.faces()[faceI];
+
+        forAll(f, fp)
+        {
+            splitMaster[f[fp]] = -2;
+        }
+    }
+
+
+    // Collect into labelList
+
+    label nSplitPoints = 0;
+
+    forAll(splitMaster, pointI)
+    {
+        if (splitMaster[pointI] >= 0)
+        {
+            nSplitPoints++;
+        }
+    }
+
+    labelList splitPoints(nSplitPoints);
+    nSplitPoints = 0;
+
+    forAll(splitMaster, pointI)
+    {
+        if (splitMaster[pointI] >= 0)
+        {
+            splitPoints[nSplitPoints++] = pointI;
+        }
+    }
+
+    return splitPoints;
+}
+
+
+//void Foam::polyRef::markIndex
+//(
+//    const label maxLevel,
+//    const label level,
+//    const label index,
+//    const label markValue,
+//    labelList& indexValues
+//) const
+//{
+//    if (level < maxLevel && indexValues[index] == -1)
+//    {
+//        // Mark
+//        indexValues[index] = markValue;
+//
+//        // Mark parent
+//        const splitCell8& split = history_.splitCells()[index];
+//
+//        if (split.parent_ >= 0)
+//        {
+//            markIndex
+//            (
+//              maxLevel, level+1, split.parent_, markValue, indexValues);
+//            )
+//        }
+//    }
+//}
+//
+//
+//// Get all cells which (down to level) originate from the same cell.
+//// level=0 returns cell only, level=1 returns the 8 cells this cell
+//// originates from, level=2 returns 64 cells etc.
+//// If the cell does not originate from refinement returns just itself.
+//void Foam::polyRef::markCellClusters
+//(
+//    const label maxLevel,
+//    labelList& cluster
+//) const
+//{
+//    cluster.setSize(mesh_.nCells());
+//    cluster = -1;
+//
+//    const DynamicList<splitCell8>& splitCells = history_.splitCells();
+//
+//    // Mark all splitCells down to level maxLevel with a cell originating from
+//    // it.
+//
+//    labelList indexLevel(splitCells.size(), -1);
+//
+//    forAll(visibleCells, cellI)
+//    {
+//        label index = visibleCells[cellI];
+//
+//        if (index >= 0)
+//        {
+//            markIndex(maxLevel, 0, index, cellI, indexLevel);
+//        }
+//    }
+//
+//    // Mark cells with splitCell
+//}
+
+
+Foam::labelList Foam::polyRef::consistentUnrefinement
+(
+    const labelList& pointsToUnrefine,
+    const bool maxSet
+) const
+{
+    if (debug)
+    {
+        Pout<< "polyRef::consistentUnrefinement :"
+            << " Determining 2:1 consistent unrefinement" << endl;
+    }
+
+    if (maxSet)
+    {
+        FatalErrorIn
+        (
+            "polyRef::consistentUnrefinement(const labelList&, const bool"
+        )   << "maxSet not implemented yet."
+            << abort(FatalError);
+    }
+
+    // Loop, modifying pointsToUnrefine, until no more changes to due to 2:1
+    // conflicts.
+    // maxSet = false : unselect points to refine
+    // maxSet = true: select points to refine
+
+    // Maintain boolList for pointsToUnrefine and cellsToUnrefine
+    PackedList<1> unrefinePoint(mesh_.nPoints(), 0);
+
+    forAll(pointsToUnrefine, i)
+    {
+        label pointI = pointsToUnrefine[i];
+
+        unrefinePoint.set(pointI, 1);
+    }
+
+
+    while (true)
+    {
+        // Construct cells to unrefine
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        PackedList<1> unrefineCell(mesh_.nCells(), 0);
+
+        forAll(unrefinePoint, pointI)
+        {
+            if (unrefinePoint.get(pointI) == 1)
+            {
+                const labelList& pCells = mesh_.pointCells()[pointI];
+
+                forAll(pCells, j)
+                {
+                    unrefineCell.set(pCells[j], 1);
+                }
+            }
+        }
+
+
+        label nChanged = 0;
+
+
+        // Check 2:1 consistency taking refinement into account
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        // Internal faces.
+        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+        {
+            label own = mesh_.faceOwner()[faceI];
+            label ownLevel = cellLevel_[own] - unrefineCell.get(own);
+
+            label nei = mesh_.faceNeighbour()[faceI];
+            label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
+
+            if (ownLevel < (neiLevel-1))
+            {
+                // Since was 2:1 this can only occur if own is marked for
+                // unrefinement.
+
+                if (maxSet)
+                {
+                    unrefineCell.set(nei, 1);
+                }
+                else
+                {
+                    if (unrefineCell.get(own) == 0)
+                    {
+                        FatalErrorIn("polyRef::consistentUnrefinement(..)")
+                            << "problem" << abort(FatalError);
+                    }
+
+                    unrefineCell.set(own, 0);
+                }
+                nChanged++;
+            }
+            else if (neiLevel < (ownLevel-1))
+            {
+                if (maxSet)
+                {
+                    unrefineCell.set(own, 1);
+                }
+                else
+                {
+                    if (unrefineCell.get(nei) == 0)
+                    {
+                        FatalErrorIn("polyRef::consistentUnrefinement(..)")
+                            << "problem" << abort(FatalError);
+                    }
+
+                    unrefineCell.set(nei, 0);
+                }
+                nChanged++;
+            }
+        }
+
+
+        // Coupled faces. Swap owner level to get neighbouring cell level.
+        labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+
+        forAll(neiLevel, i)
+        {
+            label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
+
+            neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
+        }
+
+        // Swap to neighbour
+        syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
+
+        forAll(neiLevel, i)
+        {
+            label faceI = i+mesh_.nInternalFaces();
+            label own = mesh_.faceOwner()[faceI];
+            label ownLevel = cellLevel_[own] - unrefineCell.get(own);
+
+            if (ownLevel < (neiLevel[i]-1))
+            {
+                if (!maxSet)
+                {
+                    if (unrefineCell.get(own) == 0)
+                    {
+                        FatalErrorIn("polyRef::consistentUnrefinement(..)")
+                            << "problem" << abort(FatalError);
+                    }
+
+                    unrefineCell.set(own, 0);
+                    nChanged++;
+                }
+            }
+            else if (neiLevel[i] < (ownLevel-1))
+            {
+                if (maxSet)
+                {
+                    if (unrefineCell.get(own) == 1)
+                    {
+                        FatalErrorIn("polyRef::consistentUnrefinement(..)")
+                            << "problem" << abort(FatalError);
+                    }
+
+                    unrefineCell.set(own, 1);
+                    nChanged++;
+                }
+            }
+        }
+
+        reduce(nChanged, sumOp<label>());
+
+        if (debug)
+        {
+            Pout<< "polyRef::consistentUnrefinement :"
+                << " Changed " << nChanged
+                << " refinement levels due to 2:1 conflicts."
+                << endl;
+        }
+
+        if (nChanged == 0)
+        {
+            break;
+        }
+
+
+        // Convert cellsToUnrefine back into points to unrefine
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        // Knock out any point whose cell neighbour cannot be unrefined.
+        forAll(unrefinePoint, pointI)
+        {
+            if (unrefinePoint.get(pointI) == 1)
+            {
+                const labelList& pCells = mesh_.pointCells()[pointI];
+
+                forAll(pCells, j)
+                {
+                    if (unrefineCell.get(pCells[j]) == 0)
+                    {
+                        unrefinePoint.set(pointI, 0);
+                        break;
+                    }
+                }
+            }
+        }
+    }
+
+
+    // Convert back to labelList.
+    label nSet = 0;
+
+    forAll(unrefinePoint, pointI)
+    {
+        if (unrefinePoint.get(pointI) == 1)
+        {
+            nSet++;
+        }
+    }
+
+    labelList newPointsToUnrefine(nSet);
+    nSet = 0;
+
+    forAll(unrefinePoint, pointI)
+    {
+        if (unrefinePoint.get(pointI) == 1)
+        {
+            newPointsToUnrefine[nSet++] = pointI;
+        }
+    }
+
+    return newPointsToUnrefine;
+}
+
+
+void Foam::polyRef::setUnrefinement
+(
+    const labelList& splitPointLabels,
+    directTopoChange& meshMod
+)
+{
+    if (!history_.active())
+    {
+        FatalErrorIn
+        (
+            "polyRef::setUnrefinement(const labelList&, polyTopoChange&)"
+        )   << "Only call if constructed with history capability"
+            << abort(FatalError);
+    }
+
+    if (debug)
+    {
+        Pout<< "polyRef::setUnrefinement :"
+            << " Checking initial mesh just to make sure" << endl;
+
+        checkMesh();
+
+        forAll(cellLevel_, cellI)
+        {
+            if (cellLevel_[cellI] < 0)
+            {
+                FatalErrorIn
+                (
+                    "polyRef::setUnrefinement(const labelList&, "
+                    "polyTopoChange&)"
+                )   << "Illegal cell level " << cellLevel_[cellI]
+                    << " for cell " << cellI
+                    << abort(FatalError);
+            }
+        }
+
+
+        // Write to sets.
+        pointSet pSet(mesh_, "splitPoints", labelHashSet(splitPointLabels));
+        pSet.write();
+
+        cellSet cSet
+        (
+            mesh_,
+            "splitPointCells",
+            splitPointLabels.size()
+        );
+
+        forAll(splitPointLabels, i)
+        {
+            const labelList& pCells = mesh_.pointCells()[splitPointLabels[i]];
+
+            forAll(pCells, j)
+            {
+                cSet.insert(pCells[j]);
+            }
+        }
+        cSet.write();
+
+        Pout<< "polyRef::setRefinement : Dumping " << pSet.size()
+            << " points and "
+            << cSet.size() << " cells for unrefinement to" << nl
+            << "    pointSet " << pSet.objectPath() << nl
+            << "    cellSet " << cSet.objectPath()
+            << endl;
+    }
+
+
+    labelList cellRegion;
+    labelList cellRegionMaster;
+    labelList facesToRemove;
+
+    {
+        labelHashSet splitFaces(12*splitPointLabels.size());
+
+        forAll(splitPointLabels, i)
+        {
+            const labelList& pFaces = mesh_.pointFaces()[splitPointLabels[i]];
+
+            forAll(pFaces, j)
+            {
+                splitFaces.insert(pFaces[j]);
+            }
+        }
+
+        // Check with faceRemover what faces will get removed. Note that this
+        // can be more (but never less) than splitFaces provided.
+        faceRemover_.compatibleRemoves
+        (
+            splitFaces.toc(),   // pierced faces
+            cellRegion,         // per cell -1 or region it is merged into
+            cellRegionMaster,   // per region the master cell
+            facesToRemove       // new faces to be removed.
+        );
+
+        if (facesToRemove.size() != splitFaces.size())
+        {
+            FatalErrorIn
+            (
+                "polyRef::setUnrefinement(const labelList&, polyTopoChange&)"
+            )   << "Ininitial set of split points to unrefine does not"
+                << " seem to be consistent or not mid points of refined cells"
+                << abort(FatalError);
+        }
+    }
+
+    // Redo the region master so it is consistent with our master.
+    // This will guarantee that the new cell (for which faceRemover uses
+    // the region master) is already compatible with our refinement structure.
+
+    forAll(splitPointLabels, i)
+    {
+        label pointI = splitPointLabels[i];
+
+        // Get original cell label
+
+        const labelList& pCells = mesh_.pointCells()[pointI];
+
+        // Check that the lowest numbered pCells is the master of the region
+        // (should be guaranteed by directRemoveFaces)
+        //if (debug)
+        {
+            label masterCellI = min(pCells);
+
+            forAll(pCells, j)
+            {
+                label cellI = pCells[j];
+
+                label region = cellRegion[cellI];
+
+                if (region == -1)
+                {
+                    FatalErrorIn("polyRef::setUnrefinement(..)")
+                        << "Ininitial set of split points to unrefine does not"
+                        << " seem to be consistent or not mid points"
+                        << " of refined cells" << nl
+                        << "cell:" << cellI << " on splitPoint " << pointI
+                        << " has no region to be merged into"
+                        << abort(FatalError);
+                }
+
+                if (masterCellI != cellRegionMaster[region])
+                {
+                    FatalErrorIn("polyRef::setUnrefinement(..)")
+                        << "cell:" << cellI << " on splitPoint:" << pointI
+                        << " in region " << region
+                        << " has master:" << cellRegionMaster[region]
+                        << " which is not the lowest numbered cell"
+                        << " among the pointCells:" << pCells
+                        << abort(FatalError);
+                }
+            }
+        }
+    }
+
+    // Insert all commands to combine cells. Never fails so don't have to
+    // test for success.
+    faceRemover_.setRefinement
+    (
+        facesToRemove,
+        cellRegion,
+        cellRegionMaster,
+        meshMod
+    );
+
+    // Remove the n cells that originated from merging around the split point
+    // and adapt cell levels (not that pointLevels stay the same since points
+    // either get removed or stay at the same position).
+    forAll(splitPointLabels, i)
+    {
+        label pointI = splitPointLabels[i];
+
+        const labelList& pCells = mesh_.pointCells()[pointI];
+
+        label masterCellI = min(pCells);
+
+        forAll(pCells, j)
+        {
+            cellLevel_[pCells[j]]--;
+        }
+
+        history_.combineCells(masterCellI, pCells);
+    }
+
+    // Mark files as changed
+    setInstance(mesh_.facesInstance());
+
+    // history_.updateMesh will take care of truncating.
+}
+
+
+// Write refinement to polyMesh directory.
+bool Foam::polyRef::write() const
+{
+    bool writeOk = cellLevel_.write() && pointLevel_.write();
+
+    if (history_.active())
+    {
+        writeOk = writeOk && history_.write();
+    }
+
+    return writeOk;
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/polyRef.H b/src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/polyRef.H
new file mode 100644
index 000000000..10d3106e1
--- /dev/null
+++ b/src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/polyRef.H
@@ -0,0 +1,558 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::polyRef
+
+Description
+    Refinement of (split) polyhedra using directTopoChange.
+
+    Each polyhedron is split by the following procedure:
+    1. Adding points at the edge centre, face centre and cell centre,
+    2. Adding cells n cells where n is the number of points of the cell,
+    3. Splitting each face into multiple faces going from:
+       existing corner point -> new edge centre point -> new face centre
+       point -> other new edge centre point (sharing the same corner point)
+    4. Adding internal faces going from:
+       new edge centre point -> new face centre point -> new cell centre
+       point -> other new face centre point (sharing the same edge)
+
+SourceFiles
+    polyRef.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+Notes
+    Generalisation of hexRef8 for polyhedral cells
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyRef_H
+#define polyRef_H
+
+#include "labelIOList.H"
+#include "face.H"
+#include "HashSet.H"
+#include "DynamicList.H"
+#include "primitivePatch.H"
+#include "removeFaces.H"
+#include "polyRefinementHistory.H"
+#include "PackedList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyMesh;
+class polyPatch;
+class directTopoChange;
+class mapPolyMesh;
+class mapDistributePolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class polyRef Declaration
+\*---------------------------------------------------------------------------*/
+
+class polyRef
+{
+    // Private data
+
+        //- Reference to underlying mesh
+        const polyMesh& mesh_;
+
+        //- Per cell the refinement level
+        labelIOList cellLevel_;
+
+        //- Per point the refinement level
+        labelIOList pointLevel_;
+
+        //- Typical edge length between unrefined points
+        const scalar level0Edge_;
+
+        //- Refinement history
+        polyRefinementHistory history_;
+
+        //- Face remover engine
+        removeFaces faceRemover_;
+
+        //- Level of saved points
+        Map<label> savedPointLevel_;
+
+        //- Level of saved cells
+        Map<label> savedCellLevel_;
+
+
+    // Private Member Functions
+
+        //- Reorder according to map
+        static void reorder
+        (
+            const labelList& map,
+            const label len,
+            const label null,
+            labelList& elems
+        );
+
+        //- Get patch and zone info
+        void getFaceInfo
+        (
+            const label faceI,
+            label& patchID,
+            label& zoneID,
+            label& zoneFlip
+        ) const;
+
+        //- Adds a face on top of existing faceI. Reverses if nessecary
+        label addFace
+        (
+            directTopoChange& meshMod,
+            const label faceI,
+            const face& newFace,
+            const label own,
+            const label nei
+        ) const;
+
+        //- Adds internal face from point. No checks on reversal
+        label addInternalFace
+        (
+            directTopoChange& meshMod,
+            const label meshFaceI,
+            const label meshPointI,
+            const face& newFace,
+            const label own,
+            const label nei
+        ) const;
+
+        //- Modifies existing faceI for either new owner/neighbour or new face
+        //  points. Reverses if nessecary
+        void modFace
+        (
+            directTopoChange& meshMod,
+            const label faceI,
+            const face& newFace,
+            const label own,
+            const label nei
+        ) const;
+
+        scalar getLevel0EdgeLength() const;
+
+        //- Get cell added to point of cellI (if any)
+        label getAnchorCell
+        (
+            const labelListList& cellAnchorPoints,
+            const labelListList& cellAddedCells,
+            const label cellI,
+            const label faceI,
+            const label pointI
+        ) const;
+
+        //- Get new owner and neighbour (in unspecified order) of pointI
+        //  on faceI
+        void getFaceNeighbours
+        (
+            const labelListList& cellAnchorPoints,
+            const labelListList& cellAddedCells,
+            const label faceI,
+            const label pointI,
+
+            label& own,
+            label& nei
+        ) const;
+
+
+        //- Get index of point with minimum pointlevel
+        label findMinLevel(const labelList& f) const;
+
+        //- Get index of point with maximum pointlevel
+        label findMaxLevel(const labelList& f) const;
+
+        //- Count number of vertices <= anchorLevel
+        label countAnchors(const labelList&, const label) const;
+
+        //- Find index of point with wantedLevel, starting from fp
+        label findLevel
+        (
+            const face& f,
+            const label startFp,
+            const bool searchForward,
+            const label wantedLevel
+        ) const;
+
+        ////- Print levels of list of points.
+        //void printLevels(Ostream&, const labelList&) const;
+
+        //- debug: check orientation of added internal face
+        static void checkInternalOrientation
+        (
+            directTopoChange& meshMod,
+            const label cellI,
+            const label faceI,
+            const point& ownPt,
+            const point& neiPt,
+            const face& newFace
+        );
+
+        //- debug: check orientation of new boundary face
+        static void checkBoundaryOrientation
+        (
+            directTopoChange& meshMod,
+            const label cellI,
+            const label faceI,
+            const point& ownPt,
+            const point& boundaryPt,
+            const face& newFace
+        );
+
+        //- If p0 and p1 are existing vertices check if edge is split and insert
+        //  splitPoint
+        void insertEdgeSplit
+        (
+            const labelList& edgeMidPoint,
+            const label p0,
+            const label p1,
+            dynamicLabelList& verts
+        ) const;
+
+        //- Store in maps correspondence from midpoint to anchors and faces
+        label storeMidPointInfo
+        (
+            const labelListList& cellAnchorPoints,
+            const labelListList& cellAddedCells,
+            const labelList& cellMidPoint,
+            const labelList& edgeMidPoint,
+            const label cellI,
+            const label faceI,
+            const bool faceOrder,
+            const label midPointI,
+            const label anchorPointI,
+            const label faceMidPointI,
+
+            Map<edge>& midPointToAnchors,
+            Map<edge>& midPointToFaceMids,
+            directTopoChange& meshMod
+        ) const;
+
+        //- Create all internal faces from an unsplit face
+        void createInternalFromSplitFace
+        (
+            const labelListList& cellAnchorPoints,
+            const labelListList& cellAddedCells,
+            const labelList& cellMidPoint,
+            const labelList& faceMidPoint,
+            const labelList& edgeMidPoint,
+            const label cellI,
+            const label faceI,
+
+            Map<edge>& midPointToAnchors,
+            Map<edge>& midPointToFaceMids,
+            directTopoChange& meshMod,
+            label& nFacesAdded
+        ) const;
+
+        //- Create all internal faces to split cellI into n cells where n is the
+        //  number of cell points
+        void createInternalFaces
+        (
+            const labelListList& cellAnchorPoints,
+            const labelListList& cellAddedCells,
+            const labelList& cellMidPoint,
+            const labelList& faceMidPoint,
+            const labelList& faceAnchorLevel,
+            const labelList& edgeMidPoint,
+            const label cellI,
+            directTopoChange& meshMod
+        ) const;
+
+        //- Store vertices from startFp upto face split point.
+        //  Used when splitting face into n faces where n is the number of
+        //  points in a face (or number of edges)
+        void walkFaceToMid
+        (
+            const labelList& edgeMidPoint,
+            const label cLevel,
+            const label faceI,
+            const label startFp,
+            dynamicLabelList& faceVerts
+        ) const;
+
+        //- Same as walkFaceToMid but now walk back
+        void walkFaceFromMid
+        (
+            const labelList& edgeMidPoint,
+            const label cLevel,
+            const label faceI,
+            const label startFp,
+            dynamicLabelList& faceVerts
+        ) const;
+
+        //- Updates refineCell so consistent 2:1 refinement. Returns local
+        //  number of cells changed
+        label faceConsistentRefinement
+        (
+            const bool maxSet,
+            PackedList<1>& refineCell
+        ) const;
+
+        //- Check wanted refinement for 2:1 consistency
+        void checkWantedRefinementLevels(const labelList&) const;
+
+
+        // Copy control
+
+            //- Disallow default bitwise copy construct
+            polyRef(const polyRef&);
+
+            //- Disallow default bitwise assignment
+            void operator=(const polyRef&);
+
+
+public:
+
+    //- Runtime type information
+    ClassName("polyRef");
+
+
+    // Constructors
+
+        //- Construct from mesh, read_if_present refinement data
+        //  (from write below)
+        polyRef(const polyMesh& mesh);
+
+        //- Construct from mesh and un/refinement data
+        polyRef
+        (
+            const polyMesh& mesh,
+            const labelList& cellLevel,
+            const labelList& pointLevel,
+            const polyRefinementHistory& history
+        );
+
+        //- Construct from mesh and refinement data.
+        polyRef
+        (
+            const polyMesh& mesh,
+            const labelList& cellLevel,
+            const labelList& pointLevel
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            const labelIOList& cellLevel() const
+            {
+                return cellLevel_;
+            }
+
+            const labelIOList& pointLevel() const
+            {
+                return pointLevel_;
+            }
+
+            const polyRefinementHistory& history() const
+            {
+                return history_;
+            }
+
+            //- Typical edge length between unrefined points
+            scalar level0EdgeLength() const
+            {
+                return level0Edge_;
+            }
+
+        // Refinement
+
+            // Get cell level such that the face has the most points that are
+            // <= level (at least three points)
+            label getAnchorLevel(const label faceI) const;
+
+            //- Helper: get points of a cell without using cellPoints addressing
+            labelList cellPoints(const label cellI) const;
+
+            //- Given valid mesh and current cell level and proposed
+            //  cells to refine calculate any clashes (due to 2:1) and return
+            //  ok list of cells to refine.
+            //  Either adds cells to refine to set (maxSet = true) or
+            //  removes cells to refine (maxSet = false)
+            labelList consistentRefinement
+            (
+                const labelList& cellsToRefine,
+                const bool maxSet
+            ) const;
+
+            //- Like consistentRefinement but slower:
+            //  - specify number of cells between consecutive refinement levels
+            //    (consistentRefinement equivalent to 1)
+            //  - specify max level difference between point-connected cells.
+            //    (-1 to disable). Note that with normal 2:1 limitation
+            //    (maxFaceDiff=1) there can be 8:1 size difference across point
+            //    connected cells so maxPointDiff allows you to make that less.
+            //  cellsToRefine : cells we're thinking about refining. It will
+            //                  extend this set. All refinement levels will be
+            //                  at least maxFaceDiff layers thick.
+            //  facesToCheck  : additional faces where to implement the
+            //                  maxFaceDiff thickness (usually only boundary
+            //                  faces)
+            labelList consistentSlowRefinement
+            (
+                const label maxFaceDiff,
+                const labelList& cellsToRefine,
+                const labelList& facesToCheck,
+                const label maxPointDiff,
+                const labelList& pointsToCheck
+            ) const;
+
+            //- Like consistentSlowRefinement but uses different meshWave
+            //  (proper distance instead of toplogical count). No point checks
+            //  yet
+            labelList consistentSlowRefinement2
+            (
+                const label maxFaceDiff,
+                const labelList& cellsToRefine,
+                const labelList& facesToCheck
+            ) const;
+
+            //- Insert refinement. All selected cells will be split into n cells
+            //  where n is the number of points per cell.
+            //  Returns per element in cells the n cells they were split into.
+            //  Guarantees that the 0th element is the original cell label.
+            //  Mapping:
+            //  -split cells: n new ones get added from original
+            //  -split faces: original gets modified; n - 1 new ones get added
+            //                from original
+            //  -added internal faces: added from original cell face (if
+            //   that was internal) or created out-of-nothing (so will not
+            //   get mapped!). Note: could make this inflate from point but
+            //   that will allocate interpolation.
+            //  -points added to split edge: added from edge start()
+            //  -midpoints added: added from cellPoints[0].
+            labelListList setRefinement
+            (
+                const labelList& cells,
+                directTopoChange&
+            );
+
+            //- Update local numbering for changed mesh
+            void updateMesh(const mapPolyMesh&);
+
+
+            // Restoring : is where other processes delete and reinsert data.
+            //             These callbacks allow this to restore the cellLevel
+            //             and pointLevel for reintroduced points.
+            //             Is not related to undoing my refinement
+
+                //- Signal points/face/cells for which to store data
+                void storeData
+                (
+                    const labelList& pointsToStore,
+                    const labelList& facesToStore,
+                    const labelList& cellsToStore
+                );
+
+                //- Update local numbering + undo
+                //  Data to restore given as new pointlabel + stored pointlabel
+                //  (i.e. what was in pointsToStore)
+                void updateMesh
+                (
+                    const mapPolyMesh&,
+                    const Map<label>& pointsToRestore,
+                    const Map<label>& facesToRestore,
+                    const Map<label>& cellsToRestore
+                );
+
+
+            //- Update local numbering for subsetted mesh.
+            //  Gets new-to-old maps. Not compatible with unrefinement.
+            void subset
+            (
+                const labelList& pointMap,
+                const labelList& faceMap,
+                const labelList& cellMap
+            );
+
+            //- Update local numbering for mesh redistribution
+            void distribute(const mapDistributePolyMesh&);
+
+            //- Debug: Check coupled mesh for correctness
+            void checkMesh() const;
+
+            //- Debug: Check 2:1 consistency across faces.
+            //  maxPointDiff==-1 : only check 2:1 across faces
+            //  maxPointDiff!=-1 : check point-connected cells.
+            void checkRefinementLevels
+            (
+                const label maxPointDiff,
+                const labelList& pointsToCheck
+            ) const;
+
+
+        // Unrefinement (undoing refinement, not arbitrary coarsening)
+
+            //- Return the points at the centre of top-level split cells
+            //  that can be unsplit.
+            labelList getSplitPoints() const;
+
+            //- Given proposed splitPoints to unrefine according to calculate
+            //  any clashes (due to 2:1) and return ok list of points to
+            //  unrefine. Either adds points to refine to set (maxSet = true)
+            //  or removes points to refine (maxSet = false)
+            labelList consistentUnrefinement
+            (
+                const labelList& pointsToUnrefine,
+                const bool maxSet
+            ) const;
+
+            //- Remove some refinement. Needs to be supplied output of
+            //  consistentUnrefinement. Only call if undoable set.
+            //  All n pointCells of a split point will be combined into
+            //  the lowest numbered cell of those n.
+            void setUnrefinement
+            (
+                const labelList& splitPointLabels,
+                directTopoChange&
+            );
+
+
+        // Write
+
+            // Set instance for mesh files
+            void setInstance(const fileName& inst);
+
+            //- Force writing refinement+history to polyMesh directory.
+            bool write() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/dynamicMesh/refinementData/polyRefinementHistory.C b/src/dynamicMesh/dynamicMesh/refinementData/polyRefinementHistory.C
new file mode 100644
index 000000000..a487cadf7
--- /dev/null
+++ b/src/dynamicMesh/dynamicMesh/refinementData/polyRefinementHistory.C
@@ -0,0 +1,1160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "DynamicList.H"
+#include "polyRefinementHistory.H"
+#include "ListOps.H"
+#include "mapPolyMesh.H"
+#include "mapDistributePolyMesh.H"
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(polyRefinementHistory, 0);
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::polyRefinementHistory::writeEntry
+(
+    const List<splitPolyCell>& splitCells,
+    const splitPolyCell& split
+)
+{
+    // Write me:
+    if (split.addedCellsPtr_.valid())
+    {
+        Pout<< "parent:" << split.parent_
+            << " subCells:" << split.addedCellsPtr_()
+            << endl;
+    }
+    else
+    {
+        Pout<< "parent:" << split.parent_
+            << " no subcells"
+            << endl;
+    }
+    if (split.parent_ >= 0)
+    {
+        Pout<< "parent data:" << endl;
+        // Write my parent
+        string oldPrefix = Pout.prefix();
+        Pout.prefix() = "  " + oldPrefix;
+        writeEntry(splitCells, splitCells[split.parent_]);
+        Pout.prefix() = oldPrefix;
+    }
+}
+
+
+void Foam::polyRefinementHistory::writeDebug
+(
+    const labelList& visibleCells,
+    const List<splitPolyCell>& splitCells
+)
+{
+    string oldPrefix = Pout.prefix();
+    Pout.prefix() = "";
+
+    forAll(visibleCells, cellI)
+    {
+        label index = visibleCells[cellI];
+
+        if (index >= 0)
+        {
+            Pout<< "Cell from refinement:" << cellI << " index:" << index
+                << endl;
+
+            string oldPrefix = Pout.prefix();
+            Pout.prefix() = "  " + oldPrefix;
+            writeEntry(splitCells, splitCells[index]);
+            Pout.prefix() = oldPrefix;
+        }
+        else
+        {
+            Pout<< "Unrefined cell:" << cellI << " index:" << index << endl;
+        }
+    }
+    Pout.prefix() = oldPrefix;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+//- Construct null
+Foam::polyRefinementHistory::splitPolyCell::splitPolyCell()
+:
+    parent_(-1),
+    addedCellsPtr_(NULL)
+{}
+
+
+//- Construct as child element of parent
+Foam::polyRefinementHistory::splitPolyCell::splitPolyCell(const label parent)
+:
+    parent_(parent),
+    addedCellsPtr_(NULL)
+{}
+
+
+//- Construct from Istream
+Foam::polyRefinementHistory::splitPolyCell::splitPolyCell(Istream& is)
+{
+    is >> *this;
+}
+
+
+//- Construct as (deep) copy.
+Foam::polyRefinementHistory::splitPolyCell::splitPolyCell(const splitPolyCell& sc)
+:
+    parent_(sc.parent_),
+    addedCellsPtr_
+    (
+        sc.addedCellsPtr_.valid()
+      ? new labelList(sc.addedCellsPtr_())
+      : NULL
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>
+(
+    Istream& is, polyRefinementHistory::splitPolyCell& sc
+)
+{
+    labelList addedCells;
+
+    is >> sc.parent_ >> addedCells;
+
+    if (addedCells.size() > 0)
+    {
+        sc.addedCellsPtr_.reset(new labelList(addedCells));
+    }
+    else
+    {
+        sc.addedCellsPtr_.reset(NULL);
+    }
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const polyRefinementHistory::splitPolyCell& sc
+)
+{
+    if (sc.addedCellsPtr_.valid())
+    {
+        return os
+            << sc.parent_
+            << token::SPACE
+            << sc.addedCellsPtr_();
+    }
+    else
+    {
+        return os << sc.parent_ << token::SPACE << labelList(0);
+    }
+}
+
+
+void Foam::polyRefinementHistory::checkIndices() const
+{
+    // Check indices.
+    forAll(visibleCells_, i)
+    {
+        if (visibleCells_[i] < 0 && visibleCells_[i] >= splitCells_.size())
+        {
+            FatalErrorIn("polyRefinementHistory::checkIndices() const")
+                << "Illegal entry " << visibleCells_[i]
+                << " in visibleCells at location" << i << nl
+                << "It points outside the range of splitCells : 0.."
+                << splitCells_.size()-1
+                << abort(FatalError);
+        }
+    }
+}
+
+
+Foam::label Foam::polyRefinementHistory::allocateSplitCell
+(
+    const label parent,
+    const label i,
+    const label nCells
+)
+{
+    label index = -1;
+
+    if (freeSplitCells_.size() > 0)
+    {
+        index = freeSplitCells_.remove();
+
+        splitCells_[index] = splitPolyCell(parent);
+    }
+    else
+    {
+        index = splitCells_.size();
+
+        splitCells_.append(splitPolyCell(parent));
+    }
+
+
+    // Update the parent field
+    if (parent >= 0)
+    {
+        splitPolyCell& parentSplit = splitCells_[parent];
+
+        if (!parentSplit.addedCellsPtr_.valid())
+        {
+            // Allocate storage on parent for the new nCells cells
+            parentSplit.addedCellsPtr_.reset(new labelList(nCells, -1));
+        }
+
+
+        // Store me on my parent
+        labelList& parentSplits = parentSplit.addedCellsPtr_();
+
+        parentSplits[i] = index;
+    }
+
+    return index;
+}
+
+
+void Foam::polyRefinementHistory::freeSplitCell(const label index)
+{
+    splitPolyCell& split = splitCells_[index];
+
+    // Make sure parent does not point to me anymore.
+    if (split.parent_ >= 0)
+    {
+        autoPtr<labelList>& subCellsPtr =
+            splitCells_[split.parent_].addedCellsPtr_;
+
+        if (subCellsPtr.valid())
+        {
+            labelList& subCells = subCellsPtr();
+
+            label myPos = findIndex(subCells, index);
+
+            if (myPos == -1)
+            {
+                FatalErrorIn("polyRefinementHistory::freeSplitCell")
+                    << "Problem: cannot find myself in"
+                    << " parents' children" << abort(FatalError);
+            }
+            else
+            {
+                subCells[myPos] = -1;
+            }
+        }
+    }
+
+    // Mark splitCell as free
+    split.parent_ = -2;
+
+    // Add to cache of free splitCells
+    freeSplitCells_.append(index);
+}
+
+
+// Mark entry in splitCells. Recursively mark its parent and subs.
+void Foam::polyRefinementHistory::markSplit
+(
+    const label index,
+    labelList& oldToNew,
+    DynamicList<splitPolyCell>& newSplitCells
+) const
+{
+    if (oldToNew[index] == -1)
+    {
+        // Not yet compacted.
+
+        const splitPolyCell& split = splitCells_[index];
+
+        oldToNew[index] = newSplitCells.size();
+        newSplitCells.append(split);
+
+        if (split.parent_ >= 0)
+        {
+            markSplit(split.parent_, oldToNew, newSplitCells);
+        }
+        if (split.addedCellsPtr_.valid())
+        {
+            const labelList& splits = split.addedCellsPtr_();
+
+            forAll(splits, i)
+            {
+                if (splits[i] >= 0)
+                {
+                    markSplit(splits[i], oldToNew, newSplitCells);
+                }
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::polyRefinementHistory::polyRefinementHistory(const IOobject& io)
+:
+    regIOobject(io)
+{
+    if
+    (
+        io.readOpt() == IOobject::MUST_READ
+     || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
+    {
+        readStream(typeName) >> *this;
+        close();
+    }
+
+    if (debug)
+    {
+        Pout<< "polyRefinementHistory::polyRefinementHistory :"
+            << " constructed history from IOobject :"
+            << " splitCells:" << splitCells_.size()
+            << " visibleCells:" << visibleCells_.size()
+            << endl;
+    }
+}
+
+
+//- Read or construct
+Foam::polyRefinementHistory::polyRefinementHistory
+(
+    const IOobject& io,
+    const List<splitPolyCell>& splitCells,
+    const labelList& visibleCells
+)
+:
+    regIOobject(io),
+    splitCells_(splitCells),
+    freeSplitCells_(0),
+    visibleCells_(visibleCells)
+{
+    if
+    (
+        io.readOpt() == IOobject::MUST_READ
+     || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
+    {
+        readStream(typeName) >> *this;
+        close();
+    }
+
+    // Check indices.
+    checkIndices();
+
+    if (debug)
+    {
+        Pout<< "polyRefinementHistory::polyRefinementHistory :"
+            << " constructed history from IOobject or components :"
+            << " splitCells:" << splitCells_.size()
+            << " visibleCells:" << visibleCells_.size()
+            << endl;
+    }
+}
+
+
+// Construct from initial number of cells (all visible)
+Foam::polyRefinementHistory::polyRefinementHistory
+(
+    const IOobject& io,
+    const label nCells
+)
+:
+    regIOobject(io),
+    freeSplitCells_(0)
+{
+    if
+    (
+        io.readOpt() == IOobject::MUST_READ
+     || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
+    {
+        readStream(typeName) >> *this;
+        close();
+    }
+    else
+    {
+        splitCells_.setSize(nCells, splitPolyCell());
+        visibleCells_ = identity(nCells);
+    }
+
+    // Check indices.
+    checkIndices();
+
+    if (debug)
+    {
+        Pout<< "polyRefinementHistory::polyRefinementHistory :"
+            << " constructed history from IOobject or initial size :"
+            << " splitCells:" << splitCells_.size()
+            << " visibleCells:" << visibleCells_.size()
+            << endl;
+    }
+}
+
+
+// Construct as copy
+Foam::polyRefinementHistory::polyRefinementHistory
+(
+    const IOobject& io,
+    const polyRefinementHistory& rh
+)
+:
+    regIOobject(io),
+    splitCells_(rh.splitCells()),
+    freeSplitCells_(rh.freeSplitCells()),
+    visibleCells_(rh.visibleCells())
+{
+    if (debug)
+    {
+        Pout<< "polyRefinementHistory::polyRefinementHistory : constructed initial"
+            << " history." << endl;
+    }
+}
+
+
+// Construct from Istream
+Foam::polyRefinementHistory::polyRefinementHistory(const IOobject& io, Istream& is)
+:
+    regIOobject(io),
+    splitCells_(is),
+    freeSplitCells_(0),
+    visibleCells_(is)
+{
+    // Check indices.
+    checkIndices();
+
+    if (debug)
+    {
+        Pout<< "polyRefinementHistory::polyRefinementHistory :"
+            << " constructed history from Istream"
+            << " splitCells:" << splitCells_.size()
+            << " visibleCells:" << visibleCells_.size()
+            << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::polyRefinementHistory::resize(const label size)
+{
+    label oldSize = visibleCells_.size();
+
+    if (debug)
+    {
+        Pout<< "polyRefinementHistory::resize from " << oldSize << " to " << size
+            << " cells" << endl;
+    }
+
+    visibleCells_.setSize(size);
+
+    // Set additional elements to -1.
+    for (label i = oldSize; i < visibleCells_.size(); ++i)
+    {
+        visibleCells_[i] = -1;
+    }
+}
+
+
+void Foam::polyRefinementHistory::updateMesh(const mapPolyMesh& map)
+{
+    if (active())
+    {
+        const labelList& reverseCellMap = map.reverseCellMap();
+
+        // Note that only the live cells need to be renumbered.
+
+        labelList newVisibleCells(map.cellMap().size(), -1);
+
+        forAll(visibleCells_, cellI)
+        {
+            if (visibleCells_[cellI] != -1)
+            {
+                label index = visibleCells_[cellI];
+
+                // Check
+                if (splitCells_[index].addedCellsPtr_.valid())
+                {
+                    FatalErrorIn
+                    (
+                        "polyRefinementHistory::updateMesh(const mapPolyMesh&)"
+                    )   << "Problem" << abort(FatalError);
+                }
+
+                label newCellI = reverseCellMap[cellI];
+
+                if (newCellI >= 0)
+                {
+                    newVisibleCells[newCellI] = index;
+                }
+            }
+        }
+
+        if (debug)
+        {
+            Pout<< "polyRefinementHistory::updateMesh : from "
+                << visibleCells_.size()
+                << " to " << newVisibleCells.size()
+                << " cells" << endl;
+        }
+
+        visibleCells_.transfer(newVisibleCells);
+    }
+}
+
+
+// Update numbering for subsetting
+void Foam::polyRefinementHistory::subset
+(
+    const labelList& pointMap,
+    const labelList& faceMap,
+    const labelList& cellMap
+)
+{
+    if (active())
+    {
+        labelList newVisibleCells(cellMap.size(), -1);
+
+        forAll(newVisibleCells, cellI)
+        {
+            label oldCellI = cellMap[cellI];
+
+            label index = visibleCells_[oldCellI];
+
+            // Check that cell is live (so its parent has no refinement)
+            if (index >= 0 && splitCells_[index].addedCellsPtr_.valid())
+            {
+                FatalErrorIn
+                (
+                    "polyRefinementHistory::subset"
+                    "(const labelList&, const labelList&, const labelList&)"
+                )   << "Problem" << abort(FatalError);
+            }
+
+            newVisibleCells[cellI] = index;
+        }
+
+        if (debug)
+        {
+            Pout<< "polyRefinementHistory::updateMesh : from "
+                << visibleCells_.size()
+                << " to " << newVisibleCells.size()
+                << " cells" << endl;
+        }
+
+        visibleCells_.transfer(newVisibleCells);
+    }
+}
+
+
+void Foam::polyRefinementHistory::countProc
+(
+    const label index,
+    const label newProcNo,
+    labelList& splitCellProc,
+    labelList& splitCellNum
+) const
+{
+    if (splitCellProc[index] != newProcNo)
+    {
+        // Different destination processor from other cells using this
+        // parent. Reset count.
+        splitCellProc[index] = newProcNo;
+        splitCellNum[index] = 1;
+    }
+    else
+    {
+        splitCellNum[index]++;
+
+        // Increment parent if whole splitCell moves to same processor
+        if (splitCellNum[index] == 8)
+        {
+            Pout<< "Moving " << splitCellNum[index]
+                << " cells originating from cell " << index
+                << " from processor " << Pstream::myProcNo()
+                << " to processor " << splitCellProc[index]
+                << endl;
+
+            label parent = splitCells_[index].parent_;
+
+            if (parent >= 0)
+            {
+                string oldPrefix = Pout.prefix();
+                Pout.prefix() = "  " + oldPrefix;
+
+                countProc(parent, newProcNo, splitCellProc, splitCellNum);
+
+                Pout.prefix() = oldPrefix;
+            }
+        }
+    }
+}
+
+
+void Foam::polyRefinementHistory::distribute(const mapDistributePolyMesh& map)
+{
+    if (!active())
+    {
+        FatalErrorIn
+        (
+            "polyRefinementHistory::distribute(const mapDistributePolyMesh&)"
+        )   << "Calling distribute on inactive history" << abort(FatalError);
+    }
+
+
+    if (!Pstream::parRun())
+    {
+        return;
+    }
+
+    // Remove unreferenced history.
+    compact();
+
+    Pout<< nl << "--BEFORE:" << endl;
+    writeDebug();
+    Pout<< "---------" << nl << endl;
+
+
+    // Distribution is only partially functional.
+    // If all n cells resulting from a single parent are sent across in one
+    // go it will also send across that part of the refinement history.
+    // If however e.g. first 1 and then the other n - 1 are sent across the
+    // history will not be reconstructed.
+
+    // Determine clusters. This is per every entry in splitCells_ (that is
+    // a parent of some refinement) a label giving the processor it goes to
+    // if all its children are going to the same processor.
+
+    // Per visible cell the processor it goes to.
+    labelList destination(visibleCells_.size());
+
+    const labelListList& subCellMap = map.cellMap().subMap();
+
+    forAll(subCellMap, procI)
+    {
+        const labelList& newToOld = subCellMap[procI];
+
+        forAll(newToOld, i)
+        {
+            label oldCellI = newToOld[i];
+
+            destination[oldCellI] = procI;
+        }
+    }
+
+//Pout<< "polyRefinementHistory::distribute :"
+//    << " destination:" << destination << endl;
+
+    // Per splitCell entry the processor it moves to
+    labelList splitCellProc(splitCells_.size(), -1);
+    // Per splitCell entry the number of live cells that move to that processor
+    labelList splitCellNum(splitCells_.size(), 0);
+
+    forAll(visibleCells_, cellI)
+    {
+        label index = visibleCells_[cellI];
+
+        if (index >= 0)
+        {
+            countProc
+            (
+                splitCells_[index].parent_,
+                destination[cellI],
+                splitCellProc,
+                splitCellNum
+            );
+        }
+    }
+
+    Pout<< "polyRefinementHistory::distribute :"
+        << " splitCellProc:" << splitCellProc << endl;
+
+    Pout<< "polyRefinementHistory::distribute :"
+        << " splitCellNum:" << splitCellNum << endl;
+
+
+    // Create subsetted refinement tree consisting of all parents that
+    // move in their whole to other processor.
+    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+    {
+        Pout<< "-- Subetting for processor " << procI << endl;
+
+        // From uncompacted to compacted splitCells.
+        labelList oldToNew(splitCells_.size(), -1);
+
+        // Compacted splitCells. Similar to subset routine below.
+        DynamicList<splitPolyCell> newSplitCells(splitCells_.size());
+
+        // Loop over all entries. Note: could recurse like countProc so only
+        // visit used entries but is probably not worth it.
+
+        forAll(splitCells_, index)
+        {
+//            Pout<< "oldCell:" << index
+//                << " proc:" << splitCellProc[index]
+//                << " nCells:" << splitCellNum[index]
+//                << endl;
+
+            if (splitCellProc[index] == procI && splitCellNum[index] == 8)
+            {
+                // Entry moves in its whole to procI
+                oldToNew[index] = newSplitCells.size();
+                newSplitCells.append(splitCells_[index]);
+
+                Pout<< "Added oldCell " << index
+                    << " info " << newSplitCells[newSplitCells.size()-1]
+                    << " at position " << newSplitCells.size()-1
+                    << endl;
+            }
+        }
+
+        // Add live cells that are subsetted.
+        forAll(visibleCells_, cellI)
+        {
+            label index = visibleCells_[cellI];
+
+            if (index >= 0 && destination[cellI] == procI)
+            {
+                label parent = splitCells_[index].parent_;
+
+                Pout<< "Adding refined cell " << cellI
+                    << " since moves to "
+                    << procI << " old parent:" << parent << endl;
+
+                // Create new splitCell with parent
+                oldToNew[index] = newSplitCells.size();
+                newSplitCells.append(splitPolyCell(parent));
+            }
+        }
+
+        //forAll(oldToNew, index)
+        //{
+        //    Pout<< "old:" << index << " new:" << oldToNew[index]
+        //        << endl;
+        //}
+
+        newSplitCells.shrink();
+
+        // Renumber contents of newSplitCells
+        forAll(newSplitCells, index)
+        {
+            splitPolyCell& split = newSplitCells[index];
+
+            if (split.parent_ >= 0)
+            {
+                split.parent_ = oldToNew[split.parent_];
+            }
+            if (split.addedCellsPtr_.valid())
+            {
+                labelList& splits = split.addedCellsPtr_();
+
+                forAll(splits, i)
+                {
+                    if (splits[i] >= 0)
+                    {
+                        splits[i] = oldToNew[splits[i]];
+                    }
+                }
+            }
+        }
+
+
+        const labelList& subMap = subCellMap[procI];
+
+        // New visible cells.
+        labelList newVisibleCells(subMap.size(), -1);
+
+        forAll(subMap, newCellI)
+        {
+            label oldCellI = subMap[newCellI];
+
+            label oldIndex = visibleCells_[oldCellI];
+
+            if (oldIndex >= 0)
+            {
+                newVisibleCells[newCellI] = oldToNew[oldIndex];
+            }
+        }
+
+        //Pout<< nl << "--Subset for domain:" << procI << endl;
+        //writeDebug(newVisibleCells, newSplitCells);
+        //Pout<< "---------" << nl << endl;
+
+
+        // Send to neighbours
+        OPstream toNbr(Pstream::blocking, procI);
+        toNbr << newSplitCells << newVisibleCells;
+    }
+
+
+    // Receive from neighbours and merge
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Remove all entries. Leave storage intact.
+    splitCells_.clear();
+
+    visibleCells_.setSize(map.mesh().nCells());
+    visibleCells_ = -1;
+
+    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+    {
+        IPstream fromNbr(Pstream::blocking, procI);
+        List<splitPolyCell> newSplitCells(fromNbr);
+        labelList newVisibleCells(fromNbr);
+
+        //Pout<< nl << "--Received from domain:" << procI << endl;
+        //writeDebug(newVisibleCells, newSplitCells);
+        //Pout<< "---------" << nl << endl;
+
+
+        // newSplitCells contain indices only into newSplitCells so
+        // renumbering can be done here.
+        label offset = splitCells_.size();
+
+        Pout<< "**Renumbering data from proc " << procI << " with offset "
+            << offset << endl;
+
+        forAll(newSplitCells, index)
+        {
+            splitPolyCell& split = newSplitCells[index];
+
+            if (split.parent_ >= 0)
+            {
+                split.parent_ += offset;
+            }
+            if (split.addedCellsPtr_.valid())
+            {
+                labelList& splits = split.addedCellsPtr_();
+
+                forAll(splits, i)
+                {
+                    if (splits[i] >= 0)
+                    {
+                        splits[i] += offset;
+                    }
+                }
+            }
+
+            splitCells_.append(split);
+        }
+
+
+        // Combine visibleCell.
+        const labelList& constructMap = map.cellMap().constructMap()[procI];
+
+        forAll(newVisibleCells, i)
+        {
+            visibleCells_[constructMap[i]] = newVisibleCells[i] + offset;
+        }
+    }
+    splitCells_.shrink();
+
+    Pout<< nl << "--AFTER:" << endl;
+    writeDebug();
+    Pout<< "---------" << nl << endl;
+}
+
+
+void Foam::polyRefinementHistory::compact()
+{
+    if (debug)
+    {
+        Pout<< "polyRefinementHistory::compact() Entering with:"
+            << " freeSplitCells_:" << freeSplitCells_.size()
+            << " splitCells_:" << splitCells_.size()
+            << " visibleCells_:" << visibleCells_.size()
+            << endl;
+
+        // Check all free splitCells are marked as such
+        forAll(freeSplitCells_, i)
+        {
+            label index = freeSplitCells_[i];
+
+            if (splitCells_[index].parent_ != -2)
+            {
+                FatalErrorIn("polyRefinementHistory::compact()")
+                    << "Problem index:" << index
+                    << abort(FatalError);
+            }
+        }
+
+        // Check none of the visible cells are marked as free
+        forAll(visibleCells_, cellI)
+        {
+            if
+            (
+                visibleCells_[cellI] >= 0
+             && splitCells_[visibleCells_[cellI]].parent_ == -2
+            )
+            {
+                FatalErrorIn("polyRefinementHistory::compact()")
+                    << "Problem : visible cell:" << cellI
+                    << " is marked as being free." << abort(FatalError);
+            }
+        }
+    }
+
+    DynamicList<splitPolyCell> newSplitCells(splitCells_.size());
+
+    // From uncompacted to compacted splitCells.
+    labelList oldToNew(splitCells_.size(), -1);
+
+    // Mark all used splitCell entries. These are either indexed by visibleCells
+    // or indexed from other splitCell entries.
+
+    // Mark from visibleCells
+    forAll(visibleCells_, cellI)
+    {
+        label index = visibleCells_[cellI];
+
+        if (index >= 0)
+        {
+            // Make sure we only mark visible indices if they either have a
+            // parent or subsplits.
+            if
+            (
+                splitCells_[index].parent_ != -1
+             || splitCells_[index].addedCellsPtr_.valid()
+            )
+            {
+                markSplit(index, oldToNew, newSplitCells);
+            }
+        }
+    }
+
+    // Mark from splitCells
+    forAll(splitCells_, index)
+    {
+        if (splitCells_[index].parent_ == -2)
+        {
+            // freed cell.
+        }
+        else if
+        (
+            splitCells_[index].parent_ == -1
+         && !splitCells_[index].addedCellsPtr_.valid()
+        )
+        {
+            // recombined cell. No need to keep since no parent and no subsplits
+            // Note that gets marked if reachable from other index!
+        }
+        else
+        {
+            // Is used element.
+            markSplit(index, oldToNew, newSplitCells);
+        }
+    }
+
+
+    // Now oldToNew is fully complete and compacted elements are in
+    // newSplitCells.
+    // Renumber contents of newSplitCells and visibleCells.
+    forAll(newSplitCells, index)
+    {
+        splitPolyCell& split = newSplitCells[index];
+
+        if (split.parent_ >= 0)
+        {
+            split.parent_ = oldToNew[split.parent_];
+        }
+        if (split.addedCellsPtr_.valid())
+        {
+            labelList& splits = split.addedCellsPtr_();
+
+            forAll(splits, i)
+            {
+                if (splits[i] >= 0)
+                {
+                    splits[i] = oldToNew[splits[i]];
+                }
+            }
+        }
+    }
+
+
+    if (debug)
+    {
+        Pout<< "polyRefinementHistory::compact : compacted splitCells from "
+            << splitCells_.size() << " to " << newSplitCells.size() << endl;
+    }
+
+    splitCells_.transfer(newSplitCells);
+    freeSplitCells_.clearStorage();
+
+
+    if (debug)
+    {
+        Pout<< "polyRefinementHistory::compact() NOW:"
+            << " freeSplitCells_:" << freeSplitCells_.size()
+            << " splitCells_:" << splitCells_.size()
+            << " newSplitCells:" << newSplitCells.size()
+            << " visibleCells_:" << visibleCells_.size()
+            << endl;
+    }
+
+
+    // Adapt indices in visibleCells_
+    forAll(visibleCells_, cellI)
+    {
+        label index = visibleCells_[cellI];
+
+        if (index >= 0)
+        {
+            // Note that oldToNew can be -1 so it resets newVisibleCells.
+            visibleCells_[cellI] = oldToNew[index];
+        }
+        else
+        {
+            // Keep -1 value.
+        }
+    }
+}
+
+
+void Foam::polyRefinementHistory::writeDebug() const
+{
+    writeDebug(visibleCells_, splitCells_);
+}
+
+
+void Foam::polyRefinementHistory::storeSplit
+(
+    const label cellI,
+    const labelList& addedCells
+)
+{
+    label parentIndex = -1;
+
+    if (visibleCells_[cellI] != -1)
+    {
+        // Was already live. The current live cell becomes the
+        // parent of the cells split off from it.
+
+        parentIndex = visibleCells_[cellI];
+
+        // It is no longer live (note that actually cellI gets alive
+        // again below since is addedCells[0])
+        visibleCells_[cellI] = -1;
+    }
+    else
+    {
+        // Create 0th level. -1 parent to denote this.
+        parentIndex = allocateSplitCell(-1, -1, addedCells.size());
+    }
+
+    // Create live entries for added cells that point to the
+    // cell they were created from (parentIndex)
+    forAll(addedCells, i)
+    {
+        label addedCellI = addedCells[i];
+
+        // Create entries for the split off cells. All of them
+        // are visible.
+        visibleCells_[addedCellI] =
+            allocateSplitCell(parentIndex, i, addedCells.size());
+
+    }
+}
+
+
+void Foam::polyRefinementHistory::combineCells
+(
+    const label masterCellI,
+    const labelList& combinedCells
+)
+{
+    // Save the parent structure
+    label parentIndex = splitCells_[visibleCells_[masterCellI]].parent_;
+
+    // Remove the information for the combined cells
+    forAll(combinedCells, i)
+    {
+        label cellI = combinedCells[i];
+
+        freeSplitCell(visibleCells_[cellI]);
+        visibleCells_[cellI] = -1;
+    }
+
+    splitPolyCell& parentSplit = splitCells_[parentIndex];
+    parentSplit.addedCellsPtr_.reset(NULL);
+    visibleCells_[masterCellI] = parentIndex;
+}
+
+
+bool Foam::polyRefinementHistory::readData(Istream& is)
+{
+    is >> *this;
+    return !is.bad();
+}
+
+
+bool Foam::polyRefinementHistory::writeData(Ostream& os) const
+{
+    os << *this;
+
+    return os.good();
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, polyRefinementHistory& rh)
+{
+    rh.freeSplitCells_.clear();
+
+    is >> rh.splitCells_ >> rh.visibleCells_;
+
+    // Check indices.
+    rh.checkIndices();
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const polyRefinementHistory& rh)
+{
+    const_cast<polyRefinementHistory&>(rh).compact();
+
+    return os   << "// splitCells" << nl
+                << rh.splitCells_ << nl
+                << "// visibleCells" << nl
+                << rh.visibleCells_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/dynamicMesh/refinementData/polyRefinementHistory.H b/src/dynamicMesh/dynamicMesh/refinementData/polyRefinementHistory.H
new file mode 100644
index 000000000..4d120b7e1
--- /dev/null
+++ b/src/dynamicMesh/dynamicMesh/refinementData/polyRefinementHistory.H
@@ -0,0 +1,378 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::polyRefinementHistory
+
+Description
+    All refinement history. Used in unrefinement.
+
+    - visibleCells: valid for the current mesh and contains per cell -1
+      (cell unrefined) or an index into splitCells_.
+    - splitCells: for every split contains the parent (also index into
+      splitCells) and optionally a subsplit as n indices into splitCells, where
+      n is the number of points of a cell. Note that the numbers in splitCells
+      are not cell labels, they are purely indices into splitCells.
+
+    E.g. 2 cells, cell 1 gets refined so end up with 9 cells:
+    @verbatim
+        // splitCells
+        9
+        (
+        -1 (1 2 3 4 5 6 7 8)
+        0 0()
+        0 0()
+        0 0()
+        0 0()
+        0 0()
+        0 0()
+        0 0()
+        0 0()
+        )
+
+        // visibleCells
+        9(-1 1 2 3 4 5 6 7 8)
+    @endverbatim
+
+
+    So cell0 (visibleCells=-1) is unrefined.
+    Cells 1-8 have all valid splitCells entries which are:
+      - parent:0
+      - subsplits:0()
+
+    The parent 0 refers back to the splitcell entries.
+
+
+SourceFiles
+    polyRefinementHistory.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved
+
+Notes
+    Generalisation of refinementHistory for polyhedral cells
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyRefinementHistory_H
+#define polyRefinementHistory_H
+
+#include "DynamicList.H"
+#include "dynamicLabelList.H"
+#include "labelList.H"
+#include "FixedList.H"
+#include "SLList.H"
+#include "autoPtr.H"
+#include "regIOobject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class mapPolyMesh;
+class mapDistributePolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                     Class polyRefinementHistory Declaration
+\*---------------------------------------------------------------------------*/
+
+class polyRefinementHistory
+:
+    public regIOobject
+{
+public:
+
+    class splitPolyCell
+    {
+    public:
+
+        // Index to original splitCell this cell was refined off from
+        // -1: top level cell
+        // -2: free splitCell (so should also be in freeSplitCells_)
+        label parent_;
+
+        //- cells this cell was refined into
+        autoPtr<labelList> addedCellsPtr_;
+
+        //- Construct null (parent = -1)
+        splitPolyCell();
+
+        //- Construct from parent
+        splitPolyCell(const label parent);
+
+        //- Construct from Istream
+        splitPolyCell(Istream& is);
+
+        //- Construct as deep copy
+        splitPolyCell(const splitPolyCell&);
+
+        //- Copy operator since autoPtr otherwise 'steals' storage.
+        void operator=(const splitPolyCell& s)
+        {
+            // Check for assignment to self
+            if (this == &s)
+            {
+                FatalErrorIn
+                (
+                    "splitPolyCell::operator=(const Foam::splitPolyCell&)"
+                )   << "Attempted assignment to self"
+                    << abort(FatalError);
+            }
+
+            parent_ = s.parent_;
+
+            addedCellsPtr_.reset
+            (
+                s.addedCellsPtr_.valid()
+              ? new labelList(s.addedCellsPtr_())
+              : NULL
+            );
+        }
+
+        bool operator==(const splitPolyCell& s) const
+        {
+            if (addedCellsPtr_.valid() != s.addedCellsPtr_.valid())
+            {
+                return false;
+            }
+            else if (parent_ != s.parent_)
+            {
+                return false;
+            }
+            else if (addedCellsPtr_.valid())
+            {
+                return addedCellsPtr_() == s.addedCellsPtr_();
+            }
+            else
+            {
+                return true;
+            }
+        }
+
+        bool operator!=(const splitPolyCell& s) const
+        {
+            return !operator==(s);
+        }
+
+        friend Istream& operator>>(Istream&, splitPolyCell&);
+        friend Ostream& operator<<(Ostream&, const splitPolyCell&);
+    };
+
+
+private:
+
+    TypeName("polyRefinementHistory");
+
+    // Private data
+
+        //- Storage for splitCells
+        DynamicList<splitPolyCell> splitCells_;
+
+        //- Unused indices in splitCells
+        dynamicLabelList freeSplitCells_;
+
+        //- Currently visible cells. Indices into splitCells.
+        labelList visibleCells_;
+
+
+    // Private member functions
+
+        //- Debug write
+        static void writeEntry
+        (
+            const List<splitPolyCell>&,
+            const splitPolyCell&
+        );
+
+        //- Debug write
+        static void writeDebug
+        (
+            const labelList&,
+            const List<splitPolyCell>&
+        );
+
+        //- Check consistency of structure, i.e. indices into splitCells_.
+        void checkIndices() const;
+
+        //- Allocate a splitCell. Return index in splitCells_.
+        label allocateSplitCell
+        (
+            const label parent,
+            const label i,
+            const label nCells
+        );
+
+        //- Free a splitCell.
+        void freeSplitCell(const label index);
+
+        //- Mark entry in splitCells. Recursively mark its parent and subs.
+        void markSplit
+        (
+            const label,
+            labelList& oldToNew,
+            DynamicList<splitPolyCell>&
+        ) const;
+
+        void countProc
+        (
+            const label index,
+            const label newProcNo,
+            labelList& splitCellProc,
+            labelList& splitCellNum
+        ) const;
+
+public:
+
+    // Constructors
+
+        //- Construct (read) given an IOobject
+        polyRefinementHistory(const IOobject&);
+
+        //- Construct (read) or construct null
+        polyRefinementHistory
+        (
+            const IOobject&,
+            const List<splitPolyCell>& splitCells,
+            const labelList& visibleCells
+        );
+
+        //- Construct (read) or construct from initial number of cells
+        //  (all visible)
+        polyRefinementHistory(const IOobject&, const label nCells);
+
+        //- Construct as copy
+        polyRefinementHistory(const IOobject&, const polyRefinementHistory&);
+
+        //- Construct from Istream
+        polyRefinementHistory(const IOobject&, Istream&);
+
+
+    // Member Functions
+
+        //- Per cell in the current mesh (i.e. visible) either -1 (unrefined)
+        //  or an index into splitCells.
+        const labelList& visibleCells() const
+        {
+            return visibleCells_;
+        }
+
+        //- Storage for splitPolyCells.
+        const DynamicList<splitPolyCell>& splitCells() const
+        {
+            return splitCells_;
+        }
+
+        //- Cache of unused indices in splitCells
+        const dynamicLabelList& freeSplitCells() const
+        {
+            return freeSplitCells_;
+        }
+
+        //- Is there unrefinement history. Note that this will fall over if
+        //  there are 0 cells in the mesh. But this gives problems with
+        //  lots of other programs anyway.
+        bool active() const
+        {
+            return visibleCells_.size() > 0;
+        }
+
+        //- Get parent of cell
+        label parentIndex(const label cellI) const
+        {
+            label index = visibleCells_[cellI];
+
+            if (index < 0)
+            {
+                FatalErrorIn("polyRefinementHistory::parentIndex(const label)")
+                    << "Cell " << cellI << " is not visible"
+                    << abort(FatalError);
+            }
+            return splitCells_[index].parent_;
+        }
+
+        //- Store splitting of cell into n (number of points for a cell)
+        void storeSplit
+        (
+            const label cellI,
+            const labelList& addedCells
+        );
+
+        //- Store combining n cells into master
+        void combineCells
+        (
+            const label masterCellI,
+            const labelList& combinedCells
+        );
+
+        //- Update numbering for mesh changes
+        void updateMesh(const mapPolyMesh&);
+
+        //- Update numbering for subsetting
+        void subset
+        (
+            const labelList& pointMap,
+            const labelList& faceMap,
+            const labelList& cellMap
+        );
+
+        //- Update local numbering for mesh redistribution.
+        //  Can only distribute clusters sent across in one go; cannot
+        //  handle parts recombined in multiple passes.
+        void distribute(const mapDistributePolyMesh&);
+
+        //- Compact splitCells_. Removes all freeSplitCells_ elements.
+        void compact();
+
+        //- Extend/shrink storage. additional visibleCells_ elements get
+        //  set to -1.
+        void resize(const label nCells);
+
+        //- Debug write
+        void writeDebug() const;
+
+
+        //- ReadData function required for regIOobject read operation
+        virtual bool readData(Istream&);
+
+        //- WriteData function required for regIOobject write operation
+        virtual bool writeData(Ostream&) const;
+
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, polyRefinementHistory&);
+        friend Ostream& operator<<(Ostream&, const polyRefinementHistory&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/dynamicMesh/refinementData/refinementHistory.H b/src/dynamicMesh/dynamicMesh/refinementData/refinementHistory.H
index 92e1cb528..63dda1427 100644
--- a/src/dynamicMesh/dynamicMesh/refinementData/refinementHistory.H
+++ b/src/dynamicMesh/dynamicMesh/refinementData/refinementHistory.H
@@ -28,7 +28,7 @@ Description
     All refinement history. Used in unrefinement.
 
     - visibleCells: valid for the current mesh and contains per cell -1
-    (cell unrefined) or an index into splitCells_.
+      (cell unrefined) or an index into splitCells_.
     - splitCells: for every split contains the parent (also index into
       splitCells) and optionally a subsplit as 8 indices into splitCells.
       Note that the numbers in splitCells are not cell labels, they are purely
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index a5cd573b8..05f2aa816 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -186,6 +186,7 @@ $(derivedFvPatchFields)/uniformDensityHydrostaticPressure/uniformDensityHydrosta
 $(derivedFvPatchFields)/pulseFixedValue/pulseFixedValueFvPatchFields.C
 $(derivedFvPatchFields)/waveTransmissiveInlet/waveTransmissiveInletFvPatchFields.C
 $(derivedFvPatchFields)/noSlipWall/noSlipWallFvPatchVectorField.C
+$(derivedFvPatchFields)/noSlipMovingWall/noSlipMovingWallFvPatchVectorField.C
 
 fvsPatchFields = fields/fvsPatchFields
 $(fvsPatchFields)/fvsPatchField/fvsPatchFields.C
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/noSlipMovingWall/noSlipMovingWallFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/noSlipMovingWall/noSlipMovingWallFvPatchVectorField.C
new file mode 100644
index 000000000..5feb27a15
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/noSlipMovingWall/noSlipMovingWallFvPatchVectorField.C
@@ -0,0 +1,204 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "noSlipMovingWallFvPatchVectorField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "fvcMeshPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+noSlipMovingWallFvPatchVectorField::noSlipMovingWallFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(p, iF)
+{}
+
+
+noSlipMovingWallFvPatchVectorField::noSlipMovingWallFvPatchVectorField
+(
+    const noSlipMovingWallFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchVectorField(ptf, p, iF, mapper)
+{}
+
+
+noSlipMovingWallFvPatchVectorField::noSlipMovingWallFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchVectorField(p, iF)
+{
+    fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
+}
+
+
+noSlipMovingWallFvPatchVectorField::noSlipMovingWallFvPatchVectorField
+(
+    const noSlipMovingWallFvPatchVectorField& pivpvf
+)
+:
+    fixedValueFvPatchVectorField(pivpvf)
+{}
+
+
+noSlipMovingWallFvPatchVectorField::noSlipMovingWallFvPatchVectorField
+(
+    const noSlipMovingWallFvPatchVectorField& pivpvf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(pivpvf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void noSlipMovingWallFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const fvMesh& mesh = dimensionedInternalField().mesh();
+
+    if (mesh.changing())
+    {
+        const fvPatch& p = patch();
+        const polyPatch& pp = p.patch();
+        const pointField& oldPoints = mesh.oldPoints();
+
+        vectorField oldFc(pp.size());
+
+        forAll(oldFc, i)
+        {
+            oldFc[i] = pp[i].centre(oldPoints);
+        }
+
+        // Get wall-parallel mesh motion velocity from geometry
+        vectorField Up =
+            (pp.faceCentres() - oldFc)/mesh.time().deltaT().value();
+
+        const volVectorField& U =
+            mesh.lookupObject<volVectorField>
+            (
+                dimensionedInternalField().name()
+            );
+
+        scalarField phip =
+            p.patchField<surfaceScalarField, scalar>(fvc::meshPhi(U));
+
+        vectorField n = p.nf();
+        const scalarField& magSf = p.magSf();
+        scalarField Un = phip/(magSf + VSMALL);
+
+        // Adjust for surface-normal mesh motion flux
+        vectorField::operator=(Up + n*(Un - (n & Up)));
+    }
+
+    fixedValueFvPatchVectorField::updateCoeffs();
+}
+
+
+tmp<vectorField>
+noSlipMovingWallFvPatchVectorField::gradientInternalCoeffs() const
+{
+    const vectorField nPatch = this->patch().nf();
+
+    // Calculate the diagonal part of the transformation tensor
+    vectorField implicitCoeffs(nPatch.size(), vector::zero);
+    contractLinear(implicitCoeffs, I - nPatch*nPatch);
+
+    return -this->patch().deltaCoeffs()*implicitCoeffs;
+}
+
+
+tmp<vectorField>
+noSlipMovingWallFvPatchVectorField::gradientBoundaryCoeffs() const
+{
+    // Get necessary field
+    const vectorField& UPatch = *this;
+    const vectorField nPatch = this->patch().nf();
+    const tensorField T = I - nPatch*nPatch;
+    const vectorField UPatchInternal = this->patchInternalField();
+
+    // Calculate the explicit contribution in the loop for performance
+    vectorField explicitCoeffs(nPatch.size(), vector::zero);
+
+    forAll(explicitCoeffs, faceI)
+    {
+        const tensor& curT = T[faceI];
+        tensor curDiagT(tensor::zero);
+        expandLinear(curDiagT, contractLinear(curT));
+
+        explicitCoeffs[faceI] =
+            (curT & UPatch[faceI])
+          - (
+                (curT - curDiagT)
+              & UPatchInternal[faceI]
+            );
+    }
+
+    return this->patch().deltaCoeffs()*explicitCoeffs;
+}
+
+
+void noSlipMovingWallFvPatchVectorField::write(Ostream& os) const
+{
+    fvPatchVectorField::write(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchVectorField,
+    noSlipMovingWallFvPatchVectorField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/noSlipMovingWall/noSlipMovingWallFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/noSlipMovingWall/noSlipMovingWallFvPatchVectorField.H
new file mode 100644
index 000000000..6c6927b39
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/noSlipMovingWall/noSlipMovingWallFvPatchVectorField.H
@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::noSlipMovingWallFvPatchVectorField
+
+Description
+    Foam::noSlipMovingWallFvPatchVectorField
+
+SourceFiles
+    noSlipMovingWallFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noSlipMovingWallFvPatchVectorField_H
+#define noSlipMovingWallFvPatchVectorField_H
+
+#include "fvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class noSlipMovingWallFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class noSlipMovingWallFvPatchVectorField
+:
+    public fixedValueFvPatchVectorField
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("noSlipMovingWall");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        noSlipMovingWallFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        noSlipMovingWallFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given noSlipMovingWallFvPatchVectorField
+        //  onto a new patch
+        noSlipMovingWallFvPatchVectorField
+        (
+            const noSlipMovingWallFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        noSlipMovingWallFvPatchVectorField
+        (
+            const noSlipMovingWallFvPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new noSlipMovingWallFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        noSlipMovingWallFvPatchVectorField
+        (
+            const noSlipMovingWallFvPatchVectorField&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone
+        (
+            const DimensionedField<vector, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new noSlipMovingWallFvPatchVectorField(*this, iF)
+            );
+        }
+
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Return the matrix diagonal coefficients corresponding to the
+        //  evaluation of the gradient of this patchField
+        virtual tmp<vectorField> gradientInternalCoeffs() const;
+
+        //- Return the matrix source coefficients corresponding to the
+        //  evaluation of the gradient of this patchField
+        virtual tmp<vectorField> gradientBoundaryCoeffs() const;
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //