From 9e34a1e234dbe678749837c73e3c12e08cc3ebb7 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Fri, 12 Jan 2018 15:59:24 +0100
Subject: [PATCH] Updates to polyhedralRefinement.?

---
 .../polyhedralRefinement.C                    | 3832 +++++++++++++++--
 .../polyhedralRefinement.H                    |  142 +-
 2 files changed, 3604 insertions(+), 370 deletions(-)

diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C
index e15066faa..a2345579c 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C
@@ -21,12 +21,16 @@ License
     You should have received a copy of the GNU General Public License
     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Cell layer addition/removal mesh modifier
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+Notes
+    Generalisation of hexRef8 for polyhedral cells and refactorisation into mesh
+    modifier engine.
 
 \*---------------------------------------------------------------------------*/
 
-#include "layerAdditionRemoval.H"
+#include "polyhedralRefinement.H"
 #include "polyTopoChanger.H"
 #include "polyMesh.H"
 #include "foamTime.H"
@@ -38,128 +42,3196 @@ Description
 
 namespace Foam
 {
-    defineTypeNameAndDebug(layerAdditionRemoval, 0);
+    defineTypeNameAndDebug(polyhedralRefinement, 0);
     addToRunTimeSelectionTable
     (
         polyMeshModifier,
-        layerAdditionRemoval,
+        polyhedralRefinement,
         dictionary
     );
 }
 
 
-const Foam::scalar Foam::layerAdditionRemoval::addDelta_ = 0.3;
-const Foam::scalar Foam::layerAdditionRemoval::removeDelta_ = 0.1;
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::layerAdditionRemoval::checkDefinition()
+void Foam::polyhedralRefinement::setFaceInfo
+(
+    const label faceI,
+    label& patchID,
+    label& zoneID,
+    label& zoneFlip
+) const
 {
-    if (!faceZoneID_.active())
+    patchID = -1;
+
+    if (!mesh_.isInternalFace(faceI))
     {
-        FatalErrorIn
-        (
-            "void Foam::layerAdditionRemoval::checkDefinition()"
-        )   << "Master face zone named " << faceZoneID_.name()
-            << " cannot be found."
+        patchID = mesh_.boundaryMesh().whichPatch(faceI);
+    }
+
+    zoneID = mesh_.faceZones().whichZone(faceI);
+
+    zoneFlip = false;
+
+    if (zoneID > -1)
+    {
+        const faceZone& fZone = mesh_.faceZones()[zoneID];
+
+        zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
+    }
+}
+
+
+Foam::label Foam::polyhedralRefinement::getAnchorLevel
+(
+    const label faceI
+) const
+{
+    const face& f = mesh_.faces()[faceI];
+
+    if (f.size() <= 3)
+    {
+        return pointLevel_[f[findMinMaxLevel(f, maxOp<label>())]];
+    }
+    else
+    {
+        const label& ownLevel = cellLevel_[mesh_.faceOwner()[faceI]];
+
+        if (countAnchors(f, ownLevel) >= 3)
+        {
+            return ownLevel;
+        }
+        else if (countAnchors(f, ownLevel + 1) >= 3)
+        {
+            return ownLevel + 1;
+        }
+        else
+        {
+            return -1;
+        }
+    }
+}
+
+
+Foam::scalar Foam::polyhedralRefinement::calcLevel0EdgeLength() const
+{
+    if (cellLevel_.size() != mesh_.nCells())
+    {
+        FatalErrorIn("scalar polyhedralRefinement::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);
     }
 
-    const polyMesh& mesh = topoChanger().mesh();
+    // Determine minimum edge length per refinement level
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    if (debug > 1)
+    const scalar GREAT2 = sqr(GREAT);
+
+    const label nLevels = gMax(cellLevel_) + 1;
+
+    // Initialise typical edge length squared to dummy large value
+    scalarField typEdgeLenSqr(nLevels, GREAT2);
+
+
+    // 1. Look only at edges surrounded by cellLevel cells only
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Edge list containing cellLevel of connected cells.
+    // Note:
+    // - If not set = -1
+    // - Levels of connected cells = from 0 to nLevels - 1
+    // - If different levels = labelMax
+    labelList edgeLevel(mesh_.nEdges(), -1);
+
+    // Get cell edges
+    const labelListList& meshCellEdges = mesh_.cellEdges();
+
+    // Loop through cells
+    forAll(cellLevel_, cellI)
     {
-        fileName fvPath(mesh.time().path()/"VTK");
-        mkDir(fvPath);
+        // Get cell level and cell edges
+        const label& cLevel = cellLevel_[cellI];
+        const labelList& edges = meshCellEdges[cellI];
 
-        Info<< "Writing VTK files master face zone"
-            << "Layer addition/removal " << name()
-            << ", face zone " << faceZoneID_.name()
-            << "into " << fvPath
-            << endl;
+        // Loop through all edges
+        forAll(edges, i)
+        {
+            const label& edgeI = edges[i];
 
-        primitiveFacePatch::writeVTK
-        (
-            fvPath/fileName(faceZoneID_.name() + "FaceZone"),
-            mesh.faceZones()[faceZoneID_.index()]().localFaces(),
-            mesh.faceZones()[faceZoneID_.index()]().localPoints()
-        );
-
-        primitiveFacePatch::writeVTKNormals
-        (
-            fvPath/fileName(faceZoneID_.name() + "FaceZoneNormals"),
-            mesh.faceZones()[faceZoneID_.index()]().localFaces(),
-            mesh.faceZones()[faceZoneID_.index()]().localPoints()
-        );
+            if (edgeLevel[edgeI] == -1)
+            {
+                // Edge level not set, mark it with current cell level
+                edgeLevel[edgeI] = cLevel;
+            }
+            else if (edgeLevel[edgeI] == labelMax)
+            {
+                // Already marked as on different cellLevels
+            }
+            else if (edgeLevel[edgeI] != cLevel)
+            {
+                // Edge level with different cell levels
+                edgeLevel[edgeI] = labelMax;
+            }
+        }
     }
 
-    if
+    // Make sure that edges with different levels on different processors
+    // are also marked. Do the same test (edgeLevel != cLevel) on coupled
+    // edges.
+    syncTools::syncEdgeList
     (
-        minLayerThickness_ < VSMALL
-     || maxLayerThickness_ < minLayerThickness_
-    )
+        mesh_,
+        edgeLevel,
+        ifEqEqOp<labelMax>(),
+        labelMin,
+        false               // no separation
+    );
+
+    // Now use the edgeLevel with a valid value to determine the
+    // length per level.
+
+    // Get mesh edges and points
+    const edgeList& meshEdges = mesh_.edges();
+    const pointField& meshPoints = mesh_.points();
+
+    forAll(edgeLevel, edgeI)
     {
-        FatalErrorIn
-        (
-            "void Foam::layerAdditionRemoval::checkDefinition()"
-        )   << "Incorrect layer thickness definition."
-            << abort(FatalError);
+        // Get edge level
+        const label& eLevel = edgeLevel[edgeI];
+
+        // Edge has unique cell level (both cells have the same level)
+        if (eLevel > -1 && eLevel < labelMax)
+        {
+            // Get the edge and calculate the square of the length
+            const edge& e = meshEdges[edgeI];
+            const scalar edgeLenSqr = magSqr(e.vec(meshPoints));
+
+            typEdgeLenSqr[eLevel] = min(typEdgeLenSqr[eLevel], edgeLenSqr);
+        }
     }
 
-    // Check size of zones
-    label globalZoneSize =
-        returnReduce
-        (
-           mesh.faceZones()[faceZoneID_.index()].size(),
-            sumOp<label>()
-        );
+    // Get the minimum per level over all processors. Note: using minEqOp
+    // because if the cells are not cubic, we end up using the smallest edge
+    Pstream::listCombineGather(typEdgeLenSqr, minEqOp<scalar>());
+    Pstream::listCombineScatter(typEdgeLenSqr);
 
-    if (globalZoneSize == 0)
+    if (debug)
     {
-        FatalErrorIn
+        Pout<< "polyhedralRefinement::calcLevel0EdgeLength() :"
+            << " Initial edge lengths 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)
+    //    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Initialise maximum edge length to dummy large value
+    scalarField maxEdgeLenSqr(nLevels, -GREAT2);
+
+    // Loop through cells
+    forAll(cellLevel_, cellI)
+    {
+        // Get cell level and cell edges
+        const label& cLevel = cellLevel_[cellI];
+        const labelList& cEdges = meshCellEdges[cellI];
+
+        // Loop through edges
+        forAll(cEdges, i)
+        {
+            // Get the edge
+            const edge& e = meshEdges[cEdges[i]];
+            const scalar edgeLenSqr = magSqr(e.vec(meshPoints));
+
+            maxEdgeLenSqr[cLevel] = max(maxEdgeLenSqr[cLevel], edgeLenSqr);
+        }
+    }
+
+    // Get maximum per level over all processors
+    Pstream::listCombineGather(maxEdgeLenSqr, maxEqOp<scalar>());
+    Pstream::listCombineScatter(maxEdgeLenSqr);
+
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::calcLevel0EdgeLength() :"
+            << " Basic edge lengths squared per refinementlevel:"
+            << maxEdgeLenSqr << endl;
+    }
+
+
+    // 3. Combine the two sets of lengths
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    forAll(typEdgeLenSqr, levelI)
+    {
+        if
         (
-            "void Foam::layerAdditionRemoval::checkDefinition()"
-        )   << "Face extrusion zone contains no faces.  Please check your "
-            << "mesh definition."
-            << abort(FatalError);
+            equal(typEdgeLenSqr[levelI], GREAT2)
+         && maxEdgeLenSqr[levelI] >= 0
+        )
+        {
+            typEdgeLenSqr[levelI] = maxEdgeLenSqr[levelI];
+        }
     }
 
     if (debug)
     {
-        Pout<< "Cell layer addition/removal object " << name() << " :" << nl
-            << "    faceZoneID: " << faceZoneID_ << endl;
+        Pout<< "polyRef::getLevel0EdgeLength() :"
+            << " Final edge lengths squared per refinementlevel:"
+            << typEdgeLenSqr << endl;
+    }
+
+    // Find lowest length present across all levels
+    level0EdgeLength_ = -1;
+
+    forAll(typEdgeLenSqr, levelI)
+    {
+        const scalar& lenSqr = typEdgeLenSqr[levelI];
+
+        if (lenSqr < GREAT2)
+        {
+            // Note: use power instead of left shift operator to multiply the
+            // length with appropriate level. Easier to read. VV, 5/Jan/2018.
+            level0EdgeLength_ = Foam::sqrt(lenSqr)*(Foam::pow(2, levelI));
+
+            if (debug)
+            {
+                Pout<< "polyhedralRefinement::calcLevel0EdgeLength() :"
+                    << "Edge length: " << Foam::sqrt(lenSqr) << ", "
+                    << "with edge level: " << levelI << ", "
+                    << "has equivalent level0 lenght of:" << level0EdgeLength_
+                    << endl;
+            }
+
+            break;
+        }
+    }
+
+    if (level0EdgeLength_ == -1)
+    {
+        FatalErrorIn("polyhedralRefinement::calcLevel0EdgeLength()")
+            << "Problem in definition of typical edge length squared: "
+            << typEdgeLenSqr << abort(FatalError);
     }
 }
 
-Foam::scalar Foam::layerAdditionRemoval::readOldThickness
+
+void Foam::polyhedralRefinement::setInstance(const fileName& inst)
+{
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::setInstance(const fileName& inst)"
+            << nl
+            << "Resetting file instance of refinement data to " << inst
+            << endl;
+    }
+
+    cellLevel_.instance() = inst;
+    pointLevel_.instance() = inst;
+    history_.instance() = inst;
+}
+
+
+void Foam::polyhedralRefinement::extendMarkedCells(boolList& refineCell) const
+{
+    // Mark all faces for all marked cells
+    const label nFaces = mesh_.nFaces();
+    boolList refineFace(nFaces, false);
+
+    // Get mesh cells
+    const cellList& meshCells = mesh_.cells();
+
+    // Loop through all cells
+    forAll (refineCell, cellI)
+    {
+        if (refineCell[cellI])
+        {
+            // This cell is a refinement candidate, get its faces
+            const cell& cFaces = meshCells[cellI];
+
+            forAll (cFaces, i)
+            {
+                refineFace[cFaces[i]] = true;
+            }
+        }
+    }
+
+    // Snyc the face list across processor boundaries
+    syncTools::syncFaceList(mesh_, refineFace, orEqOp<bool>(), false);
+
+    // Get necessary mesh data
+    const label nInternalFaces = mesh_.nInternalFaces();
+    const labelList& owner = mesh_.faceOwner();
+    const labelList& neighbour = mesh_.faceNeighbour();
+
+    // Internal faces
+    for (label faceI = 0; faceI < nInternalFaces; ++faceI)
+    {
+        if (refineFace[faceI])
+        {
+            // Face will be refined, mark both owner and neighbour
+            refineCell[owner[faceI]] = true;
+            refineCell[neighbour[faceI]] = true;
+        }
+    }
+
+    // Boundary faces
+    for (label faceI = nInternalFaces; faceI < nFaces; ++faceI)
+    {
+        if (refineFace[faceI])
+        {
+            // Face will be refined, mark owner
+            refineCell[owner[faceI]] = true;
+        }
+    }
+}
+
+
+void Foam::polyhedralRefinement::setPolyhedralRefinement
 (
-    const dictionary& dict
+    polyTopoChange& ref
 )
 {
-    if (dict.found("oldLayerThickness"))
+    // Note: assumes that cellsToRefine_ are set prior to the function call
+
+    // Do nothing if there are no cells to refine
+    if (cellsToRefine_.empty())
     {
-        return readScalar(dict.lookup("oldLayerThickness"));
+        Pout<< "polyehdralRefinement::setPolyhedralRefinement(...)" << nl
+            << "There are no cells selected for refinement. Returning... "
+            << endl;
+
+        return;
+    }
+
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)" << nl
+            << "Checking initial mesh before setting refinement" << endl;
+
+        checkMesh();
+    }
+
+    // New point and cell levels. Insert original lists into dynamic list for
+    // easy insertion. Note: dynamic lists shall be resized with multiplier 2 on
+    // the first insertion using operator() for non existing element
+    dynamicLabelList newCellLevel(cellLevel_);
+    dynamicLabelList newPointLevel(pointLevel_);
+
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)" << nl
+            << "Allocating " << cellsToRefine_.size() << " cell midpoints."
+            << endl;
+    }
+
+
+    // PART 1: Mark cells for refinement and add points at their cell centres
+
+    // Get necessary mesh data
+    const faceList& meshFaces = mesh_.faces();
+    const cellList& meshCells = mesh_.cells();
+    const vectorField& meshCellCentres = mesh_.cellCentres();
+
+    // Mid point for refined cell (points at cell centres):
+    // Not refined = -1
+    // Shall be refined > -1 (label of added mid point)
+    labelList cellMidPoint(mesh_.nCells(), -1);
+
+    // Loop through cells to refine
+    forAll(cellsToRefine_, i)
+    {
+        // Get cell idnex
+        const label& cellI = cellsToRefine_[i];
+
+        // Get anchor point: first point of a first face in a cell
+        const label& anchorPointI = meshFaces[meshCells[cellI][0]][0];
+
+        cellMidPoint[cellI] = ref.setAction
+        (
+            polyAddPoint
+            (
+                meshCellCentres[cellI], // Point to add (cell centre)
+                anchorPointI,           // Master point
+                -1,                     // Zone for point
+                true                    // Supports a cell
+            )
+        );
+
+        // Add a level for the new point into the dynamic list.
+        // Note: operator() resizes the list if necessary and returns the
+        // reference to newly created element
+        newPointLevel(cellMidPoint[cellI]) = cellLevel_[cellI] + 1;
+    }
+
+    // Write out split cells as a cell set for debug
+    if (debug)
+    {
+        // Note: cellSet is actually a hash table of labels
+        cellSet splitCells(mesh_, "splitCells", cellsToRefine_.size());
+
+        forAll(cellMidPoint, cellI)
+        {
+            if (cellMidPoint[cellI] > -1)
+            {
+                // Cell is marked for refinement, insert into cellSet
+                splitCells.insert(cellI);
+            }
+        }
+
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)" << nl
+            << "Writing " << splitCells.size()
+            << " cells to split to cellSet " << splitCells.objectPath()
+            << endl;
+
+        splitCells.write();
+    }
+
+
+    // PART 2: Mark edges for refinement and add points to edge centres
+
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)" << nl
+            << "Allocating edge midpoints."
+            << endl;
+    }
+
+    // Refined edges are defined by having both their point levels lower than
+    // the cell level, i.e. if any cell that gets split uses this edge, the edge
+    // needs to be split as well
+
+    // Get necessary mesh data
+    const labelListList& meshCellEdges = mesh_.cellEdges();
+    const edgeList& meshEdges = mesh_.edges();
+
+    // Mid points for refined edge:
+    // No need to split edge = -1
+    // Label of introduced mid point > -1
+    labelList edgeMidPoint(mesh_.nEdges(), -1);
+
+    // Note: Loop over all cells instead of all edges
+    forAll(cellMidPoint, cellI)
+    {
+        // Point is marked for refinement, proceed to look at the edges
+        if (cellMidPoint[cellI] > -1)
+        {
+            // Get edges of this cell
+            const labelList& cEdges = meshCellEdges[cellI];
+
+            forAll(cEdges, i)
+            {
+                // Get edge index and edge
+                const label& edgeI = cEdges[i];
+                const edge& e = meshEdges[edgeI];
+
+                if
+                (
+                    pointLevel_[e[0]] <= cellLevel_[cellI]
+                 && pointLevel_[e[1]] <= cellLevel_[cellI]
+                )
+                {
+                    // Point levels of both edge points are <= cell level, mark
+                    // the edge for splitting
+                    edgeMidPoint[edgeI] = 12345;
+                }
+            }
+        }
+    }
+
+    // Synchronize edgeMidPoint across coupled patches. Note: use 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();
+
+    // Memory management
+    {
+        // Phase 1: calculate midpoints and sync. This is necessary if we don't
+        // use binary format for writing and we slowly get differences.
+
+        // Allocate storage for edge points
+        pointField edgeMids(mesh_.nEdges(), point(-GREAT, -GREAT, -GREAT));
+
+        forAll(edgeMidPoint, edgeI)
+        {
+            if (edgeMidPoint[edgeI] > -1)
+            {
+                // Edge marked to be split. Get edge centre
+                edgeMids[edgeI] = meshEdges[edgeI].centre(meshPoints);
+            }
+        }
+
+        // Sync across processor boundaries
+        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] > -1)
+            {
+                // Edge marked to be split: replace edgeMidPoint with actual
+                // point label
+                const edge& e = meshEdges[edgeI];
+
+                edgeMidPoint[edgeI] = ref.setAction
+                (
+                    polyAddPoint
+                    (
+                        edgeMids[edgeI], // Point
+                        e[0],            // Master point
+                        -1,              // Zone for point
+                        true             // Supports a cell
+                    )
+                );
+
+                // Update level new point level: take maximum of the two points
+                // in the original edge and increment. Note: operator() resizes
+                // the list if necessary
+                newPointLevel(edgeMidPoint[edgeI]) =
+                    max
+                    (
+                        pointLevel_[e[0]],
+                        pointLevel_[e[1]]
+                    )
+                  + 1;
+            }
+        }
+    } // End memory management for syncing/adding edge points
+
+    // Write out edge mid points for split edges for debugging
+    if (debug)
+    {
+        OFstream str(mesh_.time().path()/"edgeMidPoint.obj");
+
+        forAll(edgeMidPoint, edgeI)
+        {
+            if (edgeMidPoint[edgeI] > -1)
+            {
+                // Get edge and write its cell centre
+                const edge& e = meshEdges[edgeI];
+                meshTools::writeOBJ(str, e.centre(meshPoints));
+            }
+        }
+
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)" <<
+            << "Writing centres of edges to split to file " << str.name()
+            << endl;
+    }
+
+
+    // PART 3: Calculate face level (after selected cells splitting)
+
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement" << nl
+            << "Allocating face midpoints."
+            << endl;
+    }
+
+    // Face anchor level. There are guaranteed at least 3 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);
+    }
+
+    // Mid points for refined face (points at face centres):
+    // Not refined = -1
+    // Shall be refined > -1 (label of added 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
+    // the face itself is 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] > -1 ? 1 : 0);
+
+        const label& nei = meshFaceNeighbour[faceI];
+        const label& neiLevel = cellLevel_[nei];
+        const label newNeiLevel = neiLevel + (cellMidPoint[nei] > -1 ? 1 : 0);
+
+        if
+        (
+            newOwnLevel > faceAnchorLevel[faceI]
+         || newNeiLevel > faceAnchorLevel[faceI]
+        )
+        {
+            // New level is higher than the face anchor level, mark for
+            // splitting
+            faceMidPoint[faceI] = 12345;
+        }
+    }
+
+    // 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)
+
+    // Memory management
+    {
+        // Create list for swapping boundary data
+        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] > -1 ? 1 : 0);
+
+            newNeiLevel[i] = newOwnLevel;
+        }
+
+        // Swap the list which now contains data from the other side
+        syncTools::swapBoundaryFaceList(mesh_, newNeiLevel, false);
+
+        forAll(newNeiLevel, i)
+        {
+            // Get face index
+            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] > -1 ? 1 : 0);
+
+            if
+            (
+                newOwnLevel > faceAnchorLevel[faceI]
+             || newNeiLevel[i] > faceAnchorLevel[faceI]
+            )
+            {
+                // New level is higher than the face anchor level, mark for
+                // splitting
+                faceMidPoint[faceI] = 12345;
+            }
+        }
+    } // End memory management for syncing owner/neighbour face levels
+
+    // Note: synronisation of faceMidPoints across coupled patches is not
+    // necessary since we have exchanged the neighbour data above using
+    // swapBoundaryFaceList, thus the faceMidPoint has to be the same on both
+    // sides. VV, 5/Jan/2018.
+//    // Synchronize faceMidPoint across coupled patches
+//    syncTools::syncFaceList
+//    (
+//        mesh_,
+//        faceMidPoint,
+//        maxEqOp<label>(),
+//        false
+//    );
+
+
+    // Introduce face points
+
+    // Get face centres
+    const vectorField& meshFaceCentres = mesh_.faceCentres();
+
+    // Memory management
+    {
+        // Phase 1: determine mid points and sync. Note: the same procedure has
+        // been used for syncing edge mid points
+
+        // Allocate storage for boundary face points
+        pointField bFaceMids
+        (
+            mesh_.nFaces() - mesh_.nInternalFaces(),
+            point(-GREAT, -GREAT, -GREAT)
+        );
+
+        // Loop through boundary face mids
+        forAll(bFaceMids, i)
+        {
+            // Get face index
+            const label faceI = i + mesh_.nInternalFaces();
+
+            if (faceMidPoint[faceI] > -1)
+            {
+                // This is a valid face mid, get the face centre
+                bFaceMids[i] = meshFaceCentres[faceI];
+            }
+        }
+
+        // Sync across coupled boundaries. Note: uses maximum of the components
+        // of the vector. This is completely arbitrary but it doesn't matter as
+        // long as we have same points on each side. VV, 5/Jan/2018.
+        syncTools::syncBoundaryFaceList
+        (
+            mesh_,
+            bFaceMids,
+            maxEqOp<vector>(),
+            true               // apply separation
+        );
+
+        // Loop through faces
+        forAll(faceMidPoint, faceI)
+        {
+            if (faceMidPoint[faceI] > -1)
+            {
+                // Face marked to be split. Add the point at face centre and
+                // replace faceMidPoint with actual point label
+
+                // Get mesh face
+                const face& f = meshFaces[faceI];
+
+                faceMidPoint[faceI] = ref.setAction
+                (
+                    polyAddPoint
+                    (
+                        (
+                            faceI < mesh_.nInternalFaces()
+                          ? meshFaceCentres[faceI]
+                          : bFaceMids[faceI - mesh_.nInternalFaces()]
+                        ),    // Point
+                        f[0], // Master point
+                        -1,   // Zone for point
+                        true  // Supports a cell
+                    )
+                );
+
+                // The point level of the added point will be one higher than
+                // the faceAnchorLevel (corner point). Note: operator() resizes
+                // the list if necessary
+                newPointLevel(faceMidPoint[faceI]) = faceAnchorLevel[faceI] + 1;
+            }
+        }
+    } // End memory management for syncing boundary data and adding face mids
+
+    // Write out split faces as a face set for debugging
+    if (debug)
+    {
+        // Create a faceSet with 3*cell sizes to prevent excessive resizing
+        faceSet splitFaces(mesh_, "splitFaces", 3*cellsToRefine_.size());
+
+        forAll(faceMidPoint, faceI)
+        {
+            if (faceMidPoint[faceI] > -1)
+            {
+                splitFaces.insert(faceI);
+            }
+        }
+
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)" << nl
+            << "Writing " splitFaces.size()
+            << " faces to split to faceSet " << splitFaces.objectPath()
+            << endl;
+
+        splitFaces.write();
+    }
+
+
+    // Now we have all the information we need to perform the refinement and we
+    // no longer need to refer to cellsToRefine_. The information is:
+    // - cellMidPoint >= 0 : cell needs to be split
+    // - faceMidPoint >= 0 : face needs to be split
+    // - edgeMidPoint >= 0 : edge needs to be split
+
+
+    // PART 4: Get corner and anchor points for all cells
+
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)" << nl
+            << "Finding cell anchorPoints" << 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] > -1)
+        {
+            // This cell shall be cut. Set capacity to 16 to prevent excessive
+            // resizing
+            cellAnchorPointsDynamic[cellI].setCapacity(16);
+        }
+    }
+
+    // Get necessary mesh data
+    const labelListList& meshPointCells = mesh_.pointCells();
+
+    // Loop through all points
+    forAll(pointLevel_, pointI)
+    {
+        // Get point cells
+        const labelList& pCells = meshPointCells[pointI];
+
+        // Loop through all cells sharing this point
+        forAll(pCells, pCellI)
+        {
+            // Get current cell index
+            const label& cellI = pCells[pCellI];
+
+            if
+            (
+                cellMidPoint[cellI] > -1
+             && pointLevel_[pointI] <= cellLevel_[cellI]
+            )
+            {
+                // This point cells is marked for refinement and its point level
+                // is smaller or equal to cell level, append the point
+                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)
+
+    // Get cell points for error output
+    const labelListList& meshCellPoints = mesh_.cellPoints();
+
+    forAll(cellMidPoint, cellI)
+    {
+        if (cellMidPoint[cellI] > -1)
+        {
+            // Cell selected for refinement
+            if (cellAnchorPointsDynamic[cellI].size() < 4)
+            {
+                // Cell has less than 4 anchor points. Issue an error and report
+                // cell points
+                const labelList& cPoints = meshCellPoints[cellI];
+
+                FatalErrorIn
+                (
+                    "polyhedralRefinement::setPolyhedralRefinement(...)"
+                )   << "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]);
+    }
+
+    // PART 5: Add the cells
+
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)" << nl
+            << " Adding cells."
+            << endl;
+    }
+
+    // We should have exactly n new cells per each split cell, 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] > -1)
+        {
+            // Get cell anchors
+            const labelList& cAnchors = cellAnchorPoints[cellI];
+
+            // Set the total number of added cells to number of anchors
+            labelList& cAdded = cellAddedCells[cellI];
+            cAdded.setSize(cAnchors.size());
+
+            // Original cell has index 0. Nothing to do except add the index
+            cAdded[0] = cellI;
+
+            // Update cell level
+            newCellLevel[cellI] = cellLevel_[cellI] + 1;
+
+            // Add other cells
+            for (label i = 1; i < cAdded.size(); ++i)
+            {
+                cAdded[i] = ref.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. Note: operator() resizes the list
+                // if necessary
+                newCellLevel(cAdded[i]) = cellLevel_[cellI] + 1;
+            }
+        }
+    }
+
+
+    // PART 6: Adding faces
+
+    // 6.1. Existing faces that get split (into n faces where n is the number of
+    //      points or edges)
+    // 6.2. Existing faces that do not get split but only edges get split
+    // 6.3. Existing faces that do not get split but get new owner/neighbour
+    // 6.4. New internal faces inside split cells.
+
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)" << nl
+            << " Marking faces to be handled"
+            << endl;
+    }
+
+    // Get all faces to split:
+    // a) All faces of a cell being split
+    // b) All faces that are being split
+    // c) Both faces of an edge that is being split
+    PackedBoolList facesToSplit(mesh_.nFaces(), 0);
+
+    // Get edge faces
+    const labelListList& meshEdgeFaces = mesh_.edgeFaces();
+
+    // a) All faces of a cell that is being split
+    forAll(cellMidPoint, cellI)
+    {
+        if (cellMidPoint[cellI] > -1)
+        {
+            const cell& cFaces = meshCells[cellI];
+
+            forAll(cFaces, i)
+            {
+                facesToSplit.set(cFaces[i], 1);
+            }
+        }
+    }
+
+    // b) All faces that are being split
+    forAll(faceMidPoint, faceI)
+    {
+        if (faceMidPoint[faceI] > -1)
+        {
+            facesToSplit.set(faceI, 1);
+        }
+    }
+
+    // c) Both faces of an edge that are being split
+    forAll(edgeMidPoint, edgeI)
+    {
+        if (edgeMidPoint[edgeI] > -1)
+        {
+            const labelList& eFaces = meshEdgeFaces[edgeI];
+
+            forAll(eFaces, i)
+            {
+                facesToSplit.set(eFaces[i], 1);
+            }
+        }
+    }
+
+    // Note: after splitting a certain face during parts 6.1. to 6.4.,
+    // facesToSplit for that face will be set back to 0, i.e. marked as finished
+
+
+    // PART 6.1. Add/modify faces for each face being split
+
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)" << endl
+            << "Splitting faces..." << endl;
+    }
+
+    forAll(faceMidPoint, faceI)
+    {
+        if (faceMidPoint[faceI] > -1 && facesToSplit.get(faceI) == 1)
+        {
+            // Face has not been split.
+            // Note: although facesToSplit can't be different than 1 here and
+            // the second check can be ommitted, it is left for clarity
+
+            // Get the face
+            const face& f = meshFaces[faceI];
+
+            // Flag to control whether the original faceI has been used
+            // Note: original face gets modified, other n - 1 faces are added,
+            // where n is the number of points/edges of a face
+            bool modifiedFace = false;
+
+            // Get anchor level for the face
+            const label& anchorLevel = faceAnchorLevel[faceI];
+
+            // New face always has four points/edges for arbitrary polygon
+            face newFace(4);
+
+            // Loop through all points of original face
+            forAll(f, fp)
+            {
+                const label& pointI = f[fp];
+
+                if (pointLevel_[pointI] <= anchorLevel)
+                {
+                    // This point is anchor, start collecting face
+
+                    // Create dynamic list (because of append) 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;
+                    }
+
+                    // Set new owner/neighbour indices based on split cells
+                    label own, nei;
+                    setNewFaceNeighbours
+                    (
+                        cellAnchorPoints,
+                        cellAddedCells,
+                        faceI,
+                        pointI,           // Anchor point
+
+                        own,
+                        nei
+                    );
+
+
+                    if (debug)
+                    {
+                        // Get mesh cell centres
+                        const vectorField& meshCellCentres =
+                            mesh_.cellCentres();
+
+                        if (mesh_.isInternalFace(faceI))
+                        {
+                            const label oldOwn = meshFaceOwner[faceI];
+                            const label oldNei = meshFaceNeighbour[faceI];
+
+                            checkInternalOrientation
+                            (
+                                ref,
+                                oldOwn,
+                                faceI,
+                                meshCellCentres[oldOwn],
+                                meshCellCentres[oldNei],
+                                newFace
+                            );
+                        }
+                        else
+                        {
+                            const label oldOwn = meshFaceOwner[faceI];
+
+                            checkBoundaryOrientation
+                            (
+                                ref,
+                                oldOwn,
+                                faceI,
+                                meshCellCentres[oldOwn],
+                                meshFaceCentres[faceI],
+                                newFace
+                            );
+                        }
+                    } // End debug
+
+
+                    if (!modifiedFace)
+                    {
+                        // Modify first face
+                        modifiedFace = true;
+                        modifyFace(ref, faceI, newFace, own, nei);
+                    }
+                    else
+                    {
+                        // Add additional faces
+                        addFace(ref, faceI, newFace, own, nei);
+                    }
+                } // End point anchor check
+            } // End for all points
+
+            // Mark face as handled
+            facesToSplit.set(faceI, 0);
+        }
+    }
+
+
+    // PART 6.2. Modify faces that do not get split but have edges that are
+    // being split
+
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)" << nl
+            << "Modifying faces with split edges..."
+            << endl;
+    }
+
+    // Get face edges
+    const labelListList& meshFaceEdges = mesh_.faceEdges();
+
+    forAll(edgeMidPoint, edgeI)
+    {
+        if (edgeMidPoint[edgeI] > -1)
+        {
+            // This is an edge that is going to be split, get edge faces
+            const labelList& eFaces = meshEdgeFaces[edgeI];
+
+            // Loop through all faces of an edge
+            forAll(eFaces, i)
+            {
+                // Get face index
+                const label faceI = eFaces[i];
+
+                // Check whether this is not a face that has been marked for
+                // splitting and that the face has not been handled yet. The
+                // second check is necessary since we go through edge faces
+                // instead of just faces
+                if (faceMidPoint[faceI] < 0 && facesToSplit.get(faceI) == 1)
+                {
+                    // This is unsplit face that has not been handled
+
+                    // Get face and face edges
+                    const face& f = meshFaces[faceI];
+                    const labelList& fEdges = meshFaceEdges[faceI];
+
+                    // Create a dynamic list containing new face vertices
+                    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] > -1)
+                        {
+                            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 = findMinMaxLevel(f, minOp<label>());
+
+                    label own, nei;
+                    setNewFaceNeighbours
+                    (
+                        cellAnchorPoints,
+                        cellAddedCells,
+                        faceI,
+                        f[anchorFp],      // Anchor point
+
+                        own,
+                        nei
+                    );
+
+
+                    if (debug)
+                    {
+                        // Get mesh cell centres
+                        const vectorField& meshCellCentres =
+                            mesh_.cellCentres();
+
+                        if (mesh_.isInternalFace(faceI))
+                        {
+                            const label oldOwn = meshFaceOwner[faceI];
+                            const label oldNei = meshFaceNeighbour[faceI];
+
+                            checkInternalOrientation
+                            (
+                                ref,
+                                oldOwn,
+                                faceI,
+                                meshCellCentres[oldOwn],
+                                meshCellCentres[oldNei],
+                                newFace
+                            );
+                        }
+                        else
+                        {
+                            const label oldOwn = meshFaceOwner[faceI];
+
+                            checkBoundaryOrientation
+                            (
+                                ref,
+                                oldOwn,
+                                faceI,
+                                meshCellCentres[oldOwn],
+                                meshFaceCentres[faceI],
+                                newFace
+                            );
+                        }
+                    } // End debug
+
+                    // Modify the face
+                    modifyFace(ref, faceI, newFace, own, nei);
+
+                    // Mark face as handled
+                    facesToSplit.set(faceI, 0);
+                }// End if unsplit, unhandled face
+            } // End for all edge faces
+        } // End if edge has been cut
+    } // End for all edges
+
+
+    // PART 6.3.: Modify faces that do not get split but whose owner/neighbour
+    // change due to splitting
+
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)" << nl
+            << " Changing owner/neighbour for otherwise unaffected faces..."
+            << endl;
+    }
+
+    forAll(facesToSplit, faceI)
+    {
+        // All remaining unnaffected faces are the ones whose owner/neighbour
+        // changed
+        if (facesToSplit.get(faceI) == 1)
+        {
+            // Get the face
+            const face& f = meshFaces[faceI];
+
+            // The point with the lowest level should be an anchor point of the
+            // neighbouring cells
+            label anchorFp = findMinMaxLevel(f, maxOp<label>());
+
+            label own, nei;
+            setNewFaceNeighbours
+            (
+                cellAnchorPoints,
+                cellAddedCells,
+                faceI,
+                f[anchorFp],      // Anchor point
+
+                own,
+                nei
+            );
+
+            // Modify the face, changing owner and neighbour
+            modifyFace(ref, faceI, f, own, nei);
+
+            // Mark face as handled
+            facesToSplit.set(faceI, 0);
+        }
+    }
+
+
+    // PART 6.4. Addd new internal faces inside split cells
+
+    // 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<< "polyhedralRefinement::setPolyhedralRefinement(...)" << nl
+            << " Adding new internal faces for split cells..."
+            << endl;
+    }
+
+    // Loop through cells
+    forAll(cellMidPoint, cellI)
+    {
+        if (cellMidPoint[cellI] > -1)
+        {
+            // Cell has been split, create internal faces
+            createInternalFaces
+            (
+                cellAnchorPoints,
+                cellAddedCells,
+                cellMidPoint,
+                faceMidPoint,
+                faceAnchorLevel,
+                edgeMidPoint,
+                cellI,
+                ref
+            );
+        }
+    }
+
+    // Debug: check minimum point index of added points, needs to be equal to
+    // number of points in the original mesh
+    if (debug)
+    {
+        label minPointI = labelMax;
+        label maxPointI = labelMin;
+
+        forAll(cellMidPoint, cellI)
+        {
+            if (cellMidPoint[cellI] > -1)
+            {
+                minPointI = min(minPointI, cellMidPoint[cellI]);
+                maxPointI = max(maxPointI, cellMidPoint[cellI]);
+            }
+        }
+        forAll(faceMidPoint, faceI)
+        {
+            if (faceMidPoint[faceI] > -1)
+            {
+                minPointI = min(minPointI, faceMidPoint[faceI]);
+                maxPointI = max(maxPointI, faceMidPoint[faceI]);
+            }
+        }
+        forAll(edgeMidPoint, edgeI)
+        {
+            if (edgeMidPoint[edgeI] > -1)
+            {
+                minPointI = min(minPointI, edgeMidPoint[edgeI]);
+                maxPointI = max(maxPointI, edgeMidPoint[edgeI]);
+            }
+        }
+
+        if (minPointI != labelMax && minPointI != mesh_.nPoints())
+        {
+            FatalErrorIn("polyhedralRefinement::setPolyhedralRefinement(...)")
+                << "Added point labels not consecutive to existing mesh points."
+                << nl
+                << "mesh_.nPoints():" << mesh_.nPoints()
+                << " minPointI: " << minPointI
+                << " maxPointI: " << maxPointI
+                << abort(FatalError);
+        }
+    }
+
+    // Update cell/point levels
+    pointLevel_.transfer(newPointLevel.shrink());
+    newPointLevel.clear();
+    cellLevel_.transfer(newCellLevel.shrink());
+    newCellLevel.clear();
+
+    // Mark files as changed
+    setInstance(mesh_.facesInstance());
+
+
+    // Update polyehdral refinement history (live split cells tree)
+    if (history_.active())
+    {
+        if (debug)
+        {
+            Pout<< "polyhedralRefinement::setPolyhedralRefinement(...)"
+                << " 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, store the split
+                history_.storeSplit(cellI, addedCells);
+            }
+        }
+    }
+}
+
+
+void Foam::polyhedralRefinement::setPolyhedralUnrefinement
+(
+    polyTopoChange& ref
+)
+{
+    // It is an error to attempt to unrefine if the history is inactive
+    if (!history_.active())
+    {
+        FatalErrorIn
+        (
+            "polyhedralRefinement::setPolyhedralUnrefinement"
+            << "(polyTopoChange& ref)"
+        )   << "Attempt to unrefine polyhedral cells without storing"
+            << " the refinement history. This is not allowed."
+            << abort(FatalError);
+    }
+
+    // Get point cells necessary for debug and face removal
+    const labelListList& meshPointCells = mesh_.pointCells();
+
+    if (debug)
+    {
+        Pout<< "polyhedralRefinement::setPolyhedralUnrefinement"
+            << "(polyTopoChange& ref)"
+            << nl
+            << "Checking initial mesh before setting unrefinement." << endl;
+
+        checkMesh();
+
+        forAll(cellLevel_, cellI)
+        {
+            if (cellLevel_[cellI] < 0)
+            {
+                FatalErrorIn
+                (
+                    "polyhedralRefinement::setPolyhedralUnrefinement"
+                    "(polyTopoChange& ref)"
+                )   << "Illegal cell level " << cellLevel_[cellI]
+                    << " for cell " << cellI
+                    << abort(FatalError);
+            }
+        }
+
+        // Write split points into a point set
+        pointSet pSet(mesh_, "splitPoints", labelHashSet(splitPointLabels_));
+        pSet.write();
+
+        // Write split point cells into a cell set
+        cellSet cSet
+        (
+            mesh_,
+            "splitPointCells",
+            splitPointLabels_.size()
+        );
+
+        forAll(splitPointLabels, i)
+        {
+            // Get point cells and insert them into cell set
+            const labelList& pCells = meshPointCells[splitPointLabels_[i]];
+
+            forAll(pCells, j)
+            {
+                cSet.insert(pCells[j]);
+            }
+        }
+        cSet.write();
+
+        Pout<< "polyhedralRefinement::setPolyhedralUnrefinement"
+            << "(polyTopoChange& ref)"
+            << nl
+            << "Writing " << pSet.size()
+            << " points and "
+            << cSet.size() << " cells for unrefinement to" << nl
+            << "pointSet " << pSet.objectPath() << nl
+            << "cellSet " << cSet.objectPath()
+            << endl;
+    }
+
+
+    // Create lists needed by face remover
+    labelList cellRegion;
+    labelList cellRegionMaster;
+    labelList facesToRemove;
+
+    // Memory management
+    {
+        // Collect split faces in the hast set, guess size to prevent excessive
+        // resizing
+        labelHashSet splitFaces(12*splitPointLabels_.size());
+
+        // Get point faces
+        const labelListList& meshPointFaces = mesh_.pointFaces();
+
+        forAll(splitPointsToUnrefine_, i)
+        {
+            // Loop through all faces of this point insert face index
+            const labelList& pFaces = meshPointFaces[splitPointsToUnrefine_[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,         // Region merged into (-1 for no region)
+            cellRegionMaster,   // Master cell for region
+            facesToRemove       // List of faces to be removed
+        );
+
+        if (facesToRemove.size() != splitFaces.size())
+        {
+            FatalErrorIn
+            (
+                "polyhedralRefinement::setPolyhedralUnrefinement"
+                "(polyTopoChange& ref)"
+            )   << "Either the initial set of split points to unrefine does not"
+                << " seem to be consistent or there are no mid points of"
+                << " refined cells."
+                << abort(FatalError);
+        }
+    }
+
+    // Insert all commands to combine cells
+    faceRemover_.setRefinement
+    (
+        facesToRemove,
+        cellRegion,
+        cellRegionMaster,
+        meshMod
+    );
+
+    // Remove the n cells that originated from merging around the split point
+    // and adapt cell levels. Note that pointLevels stay the same since points
+    // either get removed or stay at the same position.
+    forAll(splitPointsToUnrefine_, i)
+    {
+        // Get point index and point cells
+        const label& pointI = splitPointsToUnrefine_[i];
+        const labelList& pCells = meshPointCells[pointI];
+
+        // Get master cell index by searching for minimum of cell labels
+        const label masterCellI = min(pCells);
+
+        // Loop through cells and decrement their cell level
+        forAll(pCells, j)
+        {
+            --cellLevel_[pCells[j]];
+        }
+
+        // Combine cells in the history
+        history_.combineCells(masterCellI, pCells);
+    }
+
+    // Mark files as changed
+    setInstance(mesh_.facesInstance());
+}
+
+
+Foam::label Foam::polyhedralRefinement::addFace
+(
+    polyTopoChange& ref,
+    const label faceI,
+    const face& newFace,
+    const label own,
+    const label nei
+) const
+{
+    // Set face information
+    label patchID, zoneID, zoneFlip;
+    setFaceInfo(faceI, patchID, zoneID, zoneFlip);
+
+    // Set new face index to -1
+    label newFaceI = -1;
+
+    if ((nei == -1) || (own < nei))
+    {
+        // Ordering is ok, add the face
+        newFaceI = ref.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
     {
-        return -1.0;
+        // Ordering is flipped, reverse face and flip 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;
+}
+
+
+Foam::label Foam::polyhedralRefinement::addInternalFace
+(
+    polyTopoChange& ref,
+    const label meshFaceI,
+    const label meshPointI,
+    const face& newFace,
+    const label own,
+    const label nei
+) const
+{
+    // Check whether this is an internal face
+    if (mesh_.isInternalFace(meshFaceI))
+    {
+        return ref.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
+    {
+        // This is not an internal face. Add face out of nothing
+        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
+            )
+        );
     }
 }
 
 
-void Foam::layerAdditionRemoval::clearAddressing() const
+void Foam::polyhedralRefinement::modifyFace
+(
+    polyTopoChange& ref,
+    const label faceI,
+    const face& newFace,
+    const label own,
+    const label nei
+) const
 {
-    // Layer removal data
-    deleteDemandDrivenData(pointsPairingPtr_);
-    deleteDemandDrivenData(facesPairingPtr_);
+    // Set face info
+    label patchID, zoneID, zoneFlip;
+    setFaceInfo(faceI, patchID, zoneID, zoneFlip);
+
+    // Get owner/neighbour addressing and mesh faces
+    const labelList& owner = mesh_.faceOwner();
+    const labelList& neighbour = mesh_.faceNeighbour();
+
+    const faceList& meshFaces = mesh_.faces();
+
+    if
+    (
+        (own != owner[faceI])
+     || (
+            mesh_.isInternalFace(faceI)
+         && (nei != neighbour[faceI])
+        )
+     || (newFace != meshFaces[faceI])
+    )
+    {
+        // Either:
+        // 1. Owner index does not correspond to mesh owner,
+        // 2. Neighbour index does not correspond to mesh neighbour,
+        // 3. New face does not correspond to mesh face
+        // So we need to modify this face
+        if ((nei == -1) || (own < nei))
+        {
+            // Ordering is ok, add the face
+            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
+        {
+            // Ordering is flipped, reverse face and flip 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
+                )
+            );
+        }
+    }
+}
+
+
+void Foam::polyhedralRefinement::createInternalFaces
+(
+    const labelListList& cellAnchorPoints,
+    const labelListList& cellAddedCells,
+    const labelList& cellMidPoint,
+    const labelList& faceMidPoint,
+    const labelList& faceAnchorLevel,
+    const labelList& edgeMidPoint,
+    const label cellI,
+
+    polyTopoChange& ref
+) const
+{
+    // Find in every face the cellLevel + 1 points (from edge subdivision)
+    // and the anchor points
+
+    // Get current cell and its level
+    const cell& meshCells = mesh_.cells()[cellI];
+    const label& cLevel = cellLevel_[cellI];
+
+    // Get mesh faces and face edges
+    const faceList& meshFaces = mesh_.faces();
+    const labelListList& meshFaceEdges = mesh_.faceEdges();
+
+    // Get mesh cell points
+    const labelListList& meshCellPoints = mesh_.cellPoints();
+
+    // Map from edge mid to anchor points
+    Map<edge> midPointToAnchors(24);
+    // Map from edge mid to face mids
+    Map<edge> midPointToFaceMids(24);
+
+    // Running count of number of internal faces added so far
+    label nFacesAdded = 0;
+
+    // Loop through faces of the cell
+    forAll(meshCells, i)
+    {
+        // Get face index
+        const label& faceI = meshCells[i];
+
+        // Get current face and its edges
+        const face& f = meshFaces[faceI];
+        const labelList& fEdges = meshFaceEdges[faceI];
+
+        // We are on the cellI side of face f. The face will have 1 or n
+        // cLevel points (where n is the number of points/edges of a face)
+        // and lots of higher numbered ones
+
+        // Index of face mid point
+        label faceMidPointI = -1;
+
+        // Get number of anchors for the face
+        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;
+
+            // Loop through face points
+            forAll(f, fp)
+            {
+                if (pointLevel_[f[fp]] <= cLevel)
+                {
+                    // Point level is smaller than cLevel + 1, this is the
+                    // anchor point
+                    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);
+
+            // Set face mid point index
+            faceMidPointI = f[faceMid];
+        }
+        else
+        {
+            // There is no face middle yet but the face will be split. Set face
+            // mid point index
+            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)
+        {
+            // Get point index
+            const label& point0 = f[fp0];
+
+            if (pointLevel_[point0] <= cLevel)
+            {
+                // This is anchor point
+
+                // Walk forward to cLevel + 1 or edgeMidPoint of this level
+                label edgeMidPointI = -1;
+
+                const label fp1 = f.fcIndex(fp0);
+
+                if (pointLevel_[f[fp1]] <= cLevel)
+                {
+                    // Another anchor: edge will be split
+                    const label& edgeI = fEdges[fp0];
+
+                    edgeMidPointI = edgeMidPoint[edgeI];
+
+                    // Sanity check
+                    if (edgeMidPointI == -1)
+                    {
+                        const labelList& cPoints = meshCellPoints[cellI];
+
+                        FatalErrorIn
+                        (
+                            "polyhedralRefinement::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,
+                    ref
+                );
+
+                if (newFaceI != -1)
+                {
+                    ++nFacesAdded;
+                }
+
+
+                // Now walk backward
+
+                label fpMin1 = f.rcIndex(fp0);
+
+                if (pointLevel_[f[fpMin1]] <= cLevel)
+                {
+                    // Another anchor: edge will be split
+                    const label& edgeI = fEdges[fpMin1];
+
+                    edgeMidPointI = edgeMidPoint[edgeI];
+
+                    // Sanity check
+                    if (edgeMidPointI == -1)
+                    {
+                        const labelList& cPoints = meshCellPoints[cellI];
+
+                        FatalErrorIn
+                        (
+                            "polyhedralRefinement::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 in face to clevel + 1
+                    const 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
+}
+
+
+Foam::label Foam::polyhedralRefinement::getAnchorCell
+(
+    const labelListList& cellAnchorPoints,
+    const labelListList& cellAddedCells,
+    const label cellI,
+    const label faceI,
+    const label pointI
+) const
+{
+    // Check whether pointI is an anchor of cellI. If it is not, check whether
+    // any other point on the face is an anchor for the cell
+
+    if (cellAnchorPoints[cellI].size() > 0)
+    {
+        const label index = findIndex(cellAnchorPoints[cellI], pointI);
+
+        if (index != -1)
+        {
+            return cellAddedCells[cellI][index];
+        }
+
+
+        // pointI is not an anchor for the cell. Maybe we already refined the
+        // face so check all the face vertices
+        const face& f = mesh_.faces()[faceI];
+
+        forAll(f, fp)
+        {
+            const label index = findIndex(cellAnchorPoints[cellI], f[fp]);
+
+            if (index != -1)
+            {
+                return cellAddedCells[cellI][index];
+            }
+        }
+
+        // Problem: the point does not seem to be an anchor for cell
+
+        // Pick up points of the cell
+        const labelList& cPoints = mesh_.cellPoints()[cellI];
+
+        Perr<< "cell: " << cellI << ", points: " << endl;
+        forAll(cPoints, i)
+        {
+            const label pointI = cPoints[i];
+
+            Perr<< "    " << pointI << " coord: " << mesh_.points()[pointI]
+                << nl;
+        }
+
+        Perr<< "cell: " << cellI << " anchorPoints: " << cellAnchorPoints[cellI]
+            << endl;
+
+        FatalErrorIn("polyhedralRefinement::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;
+    }
+}
+
+
+void Foam::polyhedralRefinement::setNewFaceNeighbours
+(
+    const labelListList& cellAnchorPoints,
+    const labelListList& cellAddedCells,
+    const label faceI,
+    const label pointI,
+
+    label& own,
+    label& nei
+) const
+{
+    // Get anchor cell for this anchor point on owner side
+    own = getAnchorCell
+    (
+        cellAnchorPoints,
+        cellAddedCells,
+        mesh_.faceOwner()[faceI],
+        faceI,
+        pointI
+    );
+
+    if (mesh_.isInternalFace(faceI))
+    {
+        // Get nachor cell for this anchor point on neighbour side
+        nei = getAnchorCell
+        (
+            cellAnchorPoints,
+            cellAddedCells,
+            mesh_.faceNeighbour()[faceI],
+            faceI,
+            pointI
+        );
+    }
+    else
+    {
+        // Boundary face: set neighbour to -1
+        nei = -1;
+    }
+}
+
+
+void Foam::polyhedralRefinement::walkFaceToMid
+(
+    const labelList& edgeMidPoint,
+    const label cLevel,
+    const label faceI,
+    const label startFp,
+    dynamicLabelList& faceVerts
+) const
+{
+    // Get the face and its edges
+    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]] > -1)
+        {
+            // Edge is split, append its mid point
+            faceVerts.append(edgeMidPoint[fEdges[fp]]);
+        }
+
+        fp = f.fcIndex(fp);
+
+        if (pointLevel_[f[fp]] <= cLevel)
+        {
+            // Found next anchor. Already appended the split point just above
+            return;
+        }
+        else if (pointLevel_[f[fp]] == cLevel + 1)
+        {
+            // Mid level, append and return
+            faceVerts.append(f[fp]);
+
+            return;
+        }
+        else if (pointLevel_[f[fp]] == cLevel + 2)
+        {
+            // Store and continue to cLevel + 1
+            faceVerts.append(f[fp]);
+        }
+    }
+}
+
+
+void Foam::polyRef::walkFaceFromMid
+(
+    const labelList& edgeMidPoint,
+    const label cLevel,
+    const label faceI,
+    const label startFp,
+    dynamicLabelList& faceVerts
+) const
+{
+    // Get the face and its edges
+    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 point, break
+            break;
+        }
+        else if (pointLevel_[f[fp]] == cLevel + 1)
+        {
+            // Mid level, append and break
+            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]] > -1)
+        {
+            // Edge is split, append its mid point
+            faceVerts.append(edgeMidPoint[fEdges[fp]]);
+        }
+
+        fp = f.fcIndex(fp);
+
+        if (fp == startFp)
+        {
+            break;
+        }
+        faceVerts.append(f[fp]);
+    }
+}
+
+
+Foam::label Foam::polyhedralRefinement::countAnchors
+(
+    const labelList& f,
+    const label anchorLevel
+) const
+{
+    label nAnchors = 0;
+
+    forAll(f, fp)
+    {
+        if (pointLevel_[f[fp]] <= anchorLevel)
+        {
+            ++nAnchors;
+        }
+    }
+    return nAnchors;
+}
+
+
+Foam::label Foam::polyhedralRefinement::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("polyhedralRefinement::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("polyhedralRefinement::findLevel(...)")
+        << "face:" << f
+        << " level:" << IndirectList<label>(pointLevel_, f)()
+        << " startFp:" << startFp
+        << " wantedLevel:" << wantedLevel
+        << abort(FatalError);
+
+    return -1;
+}
+
+
+Foam::label Foam::polyhedralRefinement::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
+{
+    // A single internal face is added per edge inbetween anchor points,
+    // i.e. one face per midPoint between anchor points. The information is
+    // stored on the midPoint and if we have enough information (finished
+    // collecting two anchors and two face mid points), we add the face.
+    // Note that this member function 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
+
+
+    // 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 cell mid and face mids
+
+        dynamicLabelList newFaceVerts(4);
+        if (faceOrder == (mesh_.faceOwner()[faceI] == cellI))
+        {
+            newFaceVerts.append(faceMidPointI);
+
+            // Check and 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;
+    }
+}
+
+
+void Foam::polyhedralRefinement::insertEdgeSplit
+(
+    const labelList& edgeMidPoint,
+    const label p0,
+    const label p1,
+    dynamicLabelList& verts
+) const
+{
+    // Get number of points
+    const label nPoints = mesh_.nPoints();
+
+    if (p0 < nPoints && p1 < nPoints)
+    {
+        label edgeI = meshTools::findEdge(mesh_, p0, p1);
+
+        if (edgeI != -1 && edgeMidPoint[edgeI] != -1)
+        {
+            verts.append(edgeMidPoint[edgeI]);
+        }
+    }
+}
+
+
+void Foam::polyhedralRefinement::checkInternalOrientation
+(
+    polyTopoChange& ref,
+    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>(ref.points(), newFace)()
+    );
+
+    const vector n(compactFace.normal(compactPoints));
+    const vector dir(neiPt - ownPt);
+
+    // Check orientation error
+    if ((dir & n) < 0)
+    {
+        FatalErrorIn
+        (
+            "polyhedralRefinement::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
+        (
+            "polyhedralRefinement::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::polyhedralRefinement::checkBoundaryOrientation
+(
+    polyTopoChange& ref,
+    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
+        (
+            "polyhedralRefinement::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
+        (
+            "polyhedralRefinement::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;
+    }
+}
+
+
+Foam::label Foam::polyhedralRefinement::faceConsistentRefinement
+(
+    boolList& cellsToRefine
+) const
+{
+    // Count number of cells that will be added
+    label nAddCells = 0;
+
+    // Get necessary mesh data
+    const label nFaces = mesh_.nFaces();
+    const label nInternalFaces = mesh_.nInternalFaces();
+
+    const labelList& owner = mesh_.faceOwner();
+    const labelList& neighbour = mesh_.faceNeighbour();
+
+    // Loop through internal faces and check consistency
+    for (label faceI = 0; faceI < nInternalFaces; ++faceI)
+    {
+        // Get owner and neighbour labels
+        const label& own = owner[faceI];
+        const label& nei = neighbour[faceI];
+
+        // Get owner and neighbour cell levels
+        // Note: If the cell is marked for refinement, the level is current
+        // level + 1, otherwise it is equal to the current level
+        const label ownLevel =
+            cellsToRefine[own] ? cellLevel_[own] + 1 : cellLevel_[own];
+        const label neiLevel =
+            cellsToRefine[nei] ? cellLevel_[nei] + 1 : cellLevel_[nei];
+
+        if (ownLevel > (neiLevel + 1))
+        {
+            // Owner level is higher than neighbour level + 1, neighbour must be
+            // marked for refinement
+            cellsToRefine[nei] = true;
+            ++nAddCells;
+        }
+        else if (neiLevel > (ownLevel + 1))
+        {
+            // Neighbour level is higher than owner level + 1, owner must be
+            // marked for refinement
+            cellsToRefine[own] = true;
+            ++nAddCells;
+        }
+    }
+
+    // Create owner level for boundary faces to prepare for swapping on coupled
+    // boundaries
+    labelList ownLevel(nFaces - nInternalFaces);
+    forAll (ownLevel, i)
+    {
+        // Get owner of the face and update owner cell levels
+        const label& own = owner[i + nInternalFaces];
+        ownLevel[i] =
+            cellsToRefine[own] ? cellLevel_[own] + 1 : cellLevel_[own];
+    }
+
+    // Swap boundary face lists (coupled boundary update)
+    syncTools::swapBoundaryFaceList(mesh_, ownLevel, false);
+
+    // Note: now the ownLevel list actually contains the neighbouring level
+    // (from the other side), use alias (reference) for clarity from now on
+    const labelList& neiLevel = ownLevel;
+
+    // Loop through boundary faces
+    forAll (neiLevel, i)
+    {
+        // Get owner of the face and owner level
+        const label& own = owner[i + nInternalFaces];
+        const label curOwnLevel =
+            cellsToRefine[own] ? cellLevel_[own] + 1 : cellLevel_[own];
+
+        if (neiLevel[i] > (curOwnLevel + 1))
+        {
+            // Neighbour level is higher than owner level + 1, owner must be
+            // marked for refinement
+            cellsToRefine[own] = true;
+            ++nAddCells;
+        }
+        // Note: other possibility (that owner level is higher than neighbour
+        // level + 1) is taken into account on the other side automatically
+    }
+
+    // Return number of added cells
+    return nAddCells;
+}
+
+
+Foam::label Foam::polyhedralRefinement::pointConsistentRefinement
+(
+    boolList& cellsToRefine
+) const
+{
+    // Count number of cells that will be added
+    label nAddCells = 0;
+
+    // Collect all points from cells to refine. Assume that 10% of mesh poitns
+    // are going to be affected to prevent excessive resizing
+    labelHashSet pointsToConsider(mesh_.nPoints()/10);
+
+    // Get cell points
+    const labelListList& meshCellPoints = mesh_.cellPoints();
+
+    // Collect all points
+    forAll (meshCellPoints, cellI)
+    {
+        // Get current cell points and insert them into hash set
+        const labelList& curPoints = meshCellPoints[cellI];
+        forAll (curPoints, pointI)
+        {
+            pointsToConsider.insert(curPoints[pointI]);
+        }
+    }
+
+    // Maximum surrounding cell refinement level for each point
+    labelList maxRefLevel(mesh_.nPoints(), 0);
+
+    // Get point cells
+    const labelListList& meshPointCells = mesh_.pointCells();
+
+    // Loop through all points and collect maximum point level for each point
+    forAllConstIter (labelHashSet, pointsToConsider, iter)
+    {
+        // Get point index
+        const label& pointI = iter.key();
+
+        // Get the cells for this point
+        const labelList& curCells = pointCells[pointI];
+
+        // Find maximum refinement level for this points
+        forAll (curCells, cellI)
+        {
+            // Get cell index and cell level
+            const label& curCellI = curCells[cellI];
+            const label curLevel =
+                cellsToRefine[curCellI]
+              ? cellLevel_[curCellI] + 1
+              : cellLevel_[curCellI];
+
+            // Set the maximum point refinement level
+            if (curLevel > maxRefLevel[pointI])
+            {
+                maxRefLevel[pointI] = curLevel;
+            }
+        }
+    }
+
+    // Sync maximum refinement level across coupled boundaries
+    syncTools::syncPointList
+    (
+        mesh_,
+        maxRefLevel,
+        maxEqOp<label>(),
+        0,   // Null value
+        true // Apply separation for parallel cyclics
+    );
+
+    // Now that the levels are synced, go through considered points and add
+    // cells to refine
+    forAllConstIter (labelHashSet, pointsToConsider, iter)
+    {
+        // Get point index
+        const label& pointI = iter.key();
+
+        // Get the cells for this point
+        const labelList& curCells = pointCells[pointI];
+
+        // Loop through these point cells and set cells for refinement which
+        // would end up having refinement level smaller than maximum level - 1
+        forAll (curCells, cellI)
+        {
+            // Get cell index, reference to refinement flag and cell level
+            const label& curCellI = curCells[cellI];
+            bool& willBeRefined = cellsToRefine[curCellI];
+            const label curLevel =
+                willBeRefined
+              ? cellLevel_[curCellI] + 1
+              : cellLevel_[curCellI];
+
+            if (curLevel < maxRefLevel[pointI] - 1)
+            {
+                if (!willBeRefined)
+                {
+                    // Cell has not been marked for refinement, mark the cell for
+                    // refinement and increment the counter
+                    willBeRefined = true;
+                    ++nAddCells;
+                }
+                else
+                {
+                    FatalErrorIn
+                    (
+                        "label polyhedralRefinement::pointConsistentRefinement"
+                        "(boolList cellsToRefine) const"
+                    )   << "Cell is marked for refinement, but the 4:1 point"
+                        << " consistency cannot be ensured." << nl
+                        << "Something went wrong before this step."
+                        << endl;
+                }
+            }
+        }
+    }
+
+    return nAddCells;
+}
+
+
+Foam::label Foam::polyhedralRefinement::faceConsistentUnrefinement
+(
+    boolList& cellsToUnrefine
+) const
+{
+    // Count number of removed cells from unrefinement
+    label nRemCells = 0;
+
+    // Get necessary mesh data
+    const label nFaces = mesh_.nFaces();
+    const label nInternalFaces = mesh_.nInternalFaces();
+
+    const labelList& owner = mesh_.faceOwner();
+    const labelList& neighbour = mesh_.faceNeighbour();
+
+    // Loop through internal faces and check consistency
+    for (label faceI = 0; faceI < nInternalFaces; ++faceI)
+    {
+        // Get owner and neighbour labels
+        const label& own = owner[faceI];
+        const label& nei = neighbour[faceI];
+
+        // Get owner and neighbour cell levels
+        // Note: If the cell is marked for unrefinement, the level is current
+        // level - 1, otherwise it is equal to the current level
+        const label ownLevel =
+            cellsToUnrefine[own] ? cellLevel_[own] - 1 : cellLevel_[own];
+        const label neiLevel =
+            cellsToUnrefine[nei] ? cellLevel_[nei] - 1 : cellLevel_[nei];
+
+        if (ownLevel < (neiLevel -1))
+        {
+            // Owner level is smaller than neighbour level - 1, we must not
+            // unrefine owner
+
+            // Check whether the cell has not been marked for unrefinement
+            if (!cellsToUnrefine[own])
+            {
+                FatalErrorIn
+                (
+                    "label polyhedralRefinement::faceConsistentUnrefinement"
+                    "(...)"
+                )   << "Cell not marked for unrefinement, indicating a"
+                    << " previous unnoticed problem with unrefinement."
+                    << abort(FatalError);
+            }
+
+            cellsToUnrefine[own] = false;
+            ++nRemCells;
+        }
+        else if (neiLevel < (ownLevel - 1))
+        {
+            // Neighbour level is smaller than owner level - 1, we must not
+            // unrefine neighbour
+
+            // Check whether the cell has not been marked for unrefinement
+            if (!cellsToUnrefine[nei])
+            {
+                FatalErrorIn
+                (
+                    "label polyhedralRefinement::faceConsistentUnrefinement"
+                    "(...)"
+                )   << "Cell not marked for unrefinement, indicating a"
+                    << " previous unnoticed problem with unrefinement."
+                    << abort(FatalError);
+            }
+
+            cellsToUnrefine[nei] = false;
+            ++nRemCells;
+        }
+    }
+
+    // Create owner level for boundary faces to prepare for swapping on coupled
+    // boundaries
+    labelList ownLevel(nFaces - nInternalFaces);
+    forAll (ownLevel, i)
+    {
+        // Get owner of the face and update owner cell levels
+        const label& own = owner[i + nInternalFaces];
+        ownLevel[i] =
+            cellsToUnrefine[own] ? cellLevel_[own] - 1 : cellLevel_[own];
+    }
+
+    // Swap boundary face lists (coupled boundary update)
+    syncTools::swapBoundaryFaceList(mesh_, ownLevel, false);
+
+    // Note: now the ownLevel list actually contains the neighbouring level
+    // (from the other side), use alias (reference) for clarity from now on
+    const labelList& neiLevel = ownLevel;
+
+    // Loop through boundary faces
+    forAll (neiLevel, i)
+    {
+        // Get owner of the face and owner level
+        const label& own = owner[i + nInternalFaces];
+        const label curOwnLevel =
+            cellsToUnrefine[own] ? cellLevel_[own] - 1 : cellLevel_[own];
+
+        if (curOwnLevel < (neiLevel[i] - 1))
+        {
+            // Owner level is smaller than neighbour level - 1, we must not
+            // unrefine owner
+
+            // Check whether the cell has not been marked for unrefinement
+            if (!cellsToUnrefine[own])
+            {
+                FatalErrorIn
+                (
+                    "label polyhedralRefinement::faceConsistentUnrefinement"
+                    "(...)"
+                )   << "Cell not marked for unrefinement, indicating a"
+                    << " previous unnoticed problem with unrefinement."
+                    << abort(FatalError);
+            }
+
+            cellsToUnrefine[own] = false;
+            ++nRemCells;
+        }
+
+        // Note: other possibility (that neighbour level is smaller than owner
+        // level - 1) is taken into account on the other side automatically
+    }
+
+    // Return number of local cells removed from unrefinement
+    return nRemCells;
+}
+
+
+Foam::label Foam::polyhedralRefinement::pointConsistentUnrefinement
+(
+    const boolList& splitPointsToUnrefine,
+    boolList& cellsToUnrefine
+) const
+{
+    // Count number of cells removed from unrefinement
+    label nRemCells = 0;
+
+    // Get necessary mesh data
+    const label nPoints = mesh_.nPoints();
+    const labelListList& meshPointCells = mesh_.pointCells();
+
+    // Minimum cell refinement level for each point. Note: initialise with
+    // labelMax
+    labelList minRefLevel(nPoints, labelMax);
+
+    // Loop through all points and collect minimum point level for each point
+    forAll (splitPointsToUnrefine, pointI)
+    {
+        // Check whether the point is marked for unrefinement
+        if (splitPointsToUnrefine[pointI])
+        {
+            // Get cells for this point
+            const labelList& curCells = meshPointCells[pointI];
+
+            // Find minimum refinement level for this points
+            forAll (curCells, cellI)
+            {
+                // Get current cell and its new level
+                const label curCellI = curCells[cellI];
+                const label curLevel =
+                    cellsToUnrefine[curCellI]
+                  ? cellLevel_[curCellI] - 1
+                  : cellLevel_[curCellI];
+
+                // If the current level is smaller than current minimum
+                // refinement level, update the minimum refinement level
+                if (curLevel < minRefLevel[pointI])
+                {
+                    minRefLevel[pointI] = curLevel;
+                }
+            }
+        }
+    }
+
+    // Sync minimum refinement level across coupled boundaries
+    syncTools::syncPointList
+    (
+        mesh_,
+        minRefLevel,
+        minEqOp<label>(),
+        0,   // Null value
+        true // Apply separation for parallel cyclics
+    );
+
+    // Now that the levels are synced, go through split point candidates and add
+    // cells to unrefine
+    forAll (splitPointsToUnrefine, pointI)
+    {
+        // Check whether the point is marked for unrefinement
+        if (splitPointsToUnrefine[pointI])
+        {
+            // Get the cells for this point
+            const labelList& curCells = meshPointCells[pointI];
+
+            // Loop through these point cells and set cells for unrefinement which
+            // would end up having refinement level greater than level + 1
+            forAll (curCells, cellI)
+            {
+                // Get current cell index and level
+                const label curCellI = curCells[cellI];
+                const label curLevel =
+                    cellsToUnrefine[curCellI]
+                  ? cellLevel_[curCellI] - 1
+                  : cellLevel_[curCellI];
+
+                if (curLevel > minRefLevel[pointI] + 1)
+                {
+                    // Check whether the cell has not been marked for unrefinement
+                    if (!cellsToUnrefine[curCellI])
+                    {
+                        // Set the cell for unrefinement and increment the
+                        // counter
+                        cellsToUnrefine[curCellI] = true;
+                        ++nRemCells;
+                    }
+                    else
+                    {
+                        FatalErrorIn
+                        (
+                            "label"
+                            "polyhedralRefinement::pointConsistentUnrefinement"
+                            "\n("
+                            "\n    const boolList& splitPointsToUnrefine,"
+                            "\n    boolList& cellsToUnrefine"
+                            "\n) const"
+                        )   << "Cell is marked for unrefinement, but the 4:1 point"
+                            << " consistency cannot be ensured." << nl
+                            << "Something went wrong before this step."
+                            << endl;
+                    }
+                }
+            }
+        }
+    }
+
+    return nRemCells;
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-Foam::layerAdditionRemoval::layerAdditionRemoval
+Foam::polyhedralRefinement::polyhedralRefinement
 (
     const word& name,
     const label index,
@@ -186,7 +3258,7 @@ Foam::layerAdditionRemoval::layerAdditionRemoval
 
 
 // Construct from dictionary
-Foam::layerAdditionRemoval::layerAdditionRemoval
+Foam::polyhedralRefinement::polyhedralRefinement
 (
     const word& name,
     const dictionary& dict,
@@ -211,7 +3283,7 @@ Foam::layerAdditionRemoval::layerAdditionRemoval
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::layerAdditionRemoval::~layerAdditionRemoval()
+Foam::polyhedralRefinement::~polyhedralRefinement()
 {
     clearAddressing();
 }
@@ -219,258 +3291,410 @@ Foam::layerAdditionRemoval::~layerAdditionRemoval()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void  Foam::layerAdditionRemoval::setRemoval()
+void Foam::polyhedralRefinement::setCellsToRefine
+(
+    const labelList& refinementCellCandidates
+)
 {
-    triggerRemoval_ = topoChanger().morphIndex();
+    if (debug)
+    {
+        Info<< "polyhedralRefinement::setCellsToRefine"
+            << "(const labelList& refinementCellCandidates)" << nl
+            << "Setting cells to refine" << endl;
+    }
+
+    if (refinementCellCandidates.empty())
+    {
+        // Set cells to refine to empty list and return
+        cellsToRefine_.clear();
+        return;
+    }
+
+    // Create a mark-up field for cells to refine
+    boolList refineCell(mesh_.nCells(), false);
+
+    // Mark initial refinement candidates for refinement
+    forAll (refinementCellCandidates, i)
+    {
+        refineCell[refinementCellCandidates[i]] = true;
+    }
+
+    // Extend cells using a specified number of buffer layers
+    for (label i = 0; i < nBufferLayers_; ++i)
+    {
+        extendMarkedCells(refineCell);
+    }
+
+    // Make sure that the refinement is face consistent (2:1 consistency) and
+    // point consistent (4:1 consistency) if necessary
+    do
+    {
+        // Number of additional cells to refine due to consistency in this
+        // iteration
+        label nAddCells = 0;
+
+        if (pointBasedRefinement_)
+        {
+            // Check for 4:1 point based consistent refinement. Updates
+            // cellsToRefine and returns number of cells added in this iteration
+            nAddCells += pointConsistentRefinement(cellsToRefine);
+        }
+
+        // Check for 2:1 face based consistent refinement. Updates cellsToRefine
+        // and returns number of cells added in this iteration
+        nAddCells += faceConsistentRefinement(cellsToRefine);
+
+        // Global reduction
+        reduce(nAddCells, sumOp<label>());
+
+        if (debug)
+        {
+            Info<< "Added " << nAddCells
+                << " cells to obtain consistent refinement."
+                << endl;
+        }
+
+    } while (nAddCells > 0);
+
+    // Collect all cells to refine in a dynamic list
+    dynamicLabelList cellsToRefineDynamic(mesh_.nCells());
+
+    forAll (refineCell, cellI)
+    {
+        if (refineCell[cellI])
+        {
+            // Cell marked for refinement, append it
+            cellsToRefineDynamic.append(cellI);
+        }
+    }
+
+    // Transfer the contents into the data member (ordinary list)
+    cellsToRefine_.transfer(cellsToRefineDynamic);
 }
 
 
-void  Foam::layerAdditionRemoval::setAddition()
+void Foam::polyhedralRefinement::setSplitPointsToUnrefine
+(
+    const labelList& unrefinementPointCandidates
+)
 {
-    triggerAddition_ = topoChanger().morphIndex();
+    if (debug)
+    {
+        Info<< "polyhedralRefinement::setSplitPointsToUnrefine"
+            << "(const labelList& unrefinementPointCandidates)" << nl
+            << "Setting cells to refine" << endl;
+    }
+
+    if (unrefinementPointCandidates.empty())
+    {
+        // Set split points to unrefine to empty list and return
+        splitPointsToUnrefine_.clear();
+        return;
+    }
+
+    // Get necessary mesh data
+    const label nPoints = mesh_.nPoints();
+    const labelListList& meshCellPoints = mesh_.cellPoints();
+
+    // PART 1: Mark all split points in the mesh (points that can be unrefined)
+    boolList splitPointsMarkup(nPoints, false);
+
+    // Field indicating whether the point has been split:
+    // -1 = undetermined yet
+    // -2 = certainly not a split point
+    // > -1 = label of master cell
+    labelList splitMaster(nPoints, -1);
+
+    // Level of split point
+    labelList splitMasterLevel(nPoints, 0);
+
+    // Get visible cells from history: -1 if the cell is unrefined and otherwise
+    // index into history::splitCells_
+    const labelList& visibleCells = history_.visibleCells();
+
+    // Unmark all points with different master cell
+    forAll (visibleCells, cellI)
+    {
+        // Get cell points of this cell
+        const labelList& cPoints = meshCellPoints[cellI];
+
+        if (visibleCells[cellI] != -1 && history_.parentIndex(cellI) > -1)
+        {
+            // Get parent index
+            const label parentIndex = history_.parentIndex(cellI);
+
+            // Check whether the master is the same
+            forAll (cPoints, i)
+            {
+                // Get point index and master cell index
+                const label& pointI = cPoints[i];
+                const label& masterCellI = splitMaster[pointI];
+
+                if (masterCellI == -1)
+                {
+                    // This is the first time we visisted the 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 definitely not split point,
+                    // see below
+                }
+                else if
+                (
+                    (masterCellI != parentIndex)
+                 || (splitMasterLevel[pointI] != cellLevel_[cellI] - 1)
+                )
+                {
+                    // Point has different masters so it can't be a split point
+                    // as it is on two different refinement patterns
+                    splitMaster[pointI] = -2;
+                }
+            } // End for all cell points
+        }
+        else
+        {
+            // The cell is either not visible or it is unrefined cell, all
+            // points can't be split points
+            forAll (cPoints, i)
+            {
+                splitMaster[cPoints[i]] = -2;
+            }
+        }
+    } // End for all visible cells
+
+    // Points at boundaries can't be split points, unmark them
+
+    // Get necessary mesh data for this step
+    const label nFaces = mesh_.nFaces();
+    const label nInternalFaces = mesh_.nInternalFaces();
+
+    const faceList& meshFaces = mesh_.faces();
+
+    // Loop through all boundary faces
+    for (label faceI < nInternalFaces; faceI < nFaces; ++faceI)
+    {
+        // Get current boundary face and mark all its points as ordinary (not
+        // split) points
+        const face& f = meshFaces[faceI];
+        forAll (f, fpI)
+        {
+            splitMaster[f[fpI]] = -2;
+        }
+    }
+
+    // Finally, mark all split points
+    forAll (splitMaster, point)
+    {
+        if (splitMaster[pointI] > -1)
+        {
+            splitPointsMarkup[pointI] = true;
+        }
+    }
+
+
+    // PART 2: Mark all unrefinement point candidates that are split points at
+    // the same time (basically the intersection of split points and candidates)
+
+    // Create markup field of split points to unrefine
+    // True: this is a split point which should be unrefined
+    // False: this is either not a split point or it shouldn't be unrefined
+    boolList splitPointsToUnrefine(nPoints, false);
+
+    // Loop through all unrefinement candidates
+    forAll (unrefinementPointCandidates, i)
+    {
+        // Get point index
+        const label& pointI = unrefinementPointCandidates[i];
+
+        if (splitPointsMarkup[pointI])
+        {
+            // This is a split point, mark it for unrefinement
+            splitPointsToUnrefine[pointI] = true;
+        }
+    }
+
+
+    // PART 3: Make sure that we skip unrefining around split points that
+    // possibly have cells around that will be refined. This might happen if
+    // someone uses an inconsistent refinement/unrefinement selection procedure
+
+    // Get cell points
+    const labelListList& meshCellPoints = mesh_.cellPoints();
+
+    // Loop through cells to refine
+    forAll (cellsToRefine_, i)
+    {
+        // Get cell index and the list of cell points
+        const label& cellI = cellsToRefine_[i];
+        const labelList& cPoints = meshCellPoints[cellI];
+
+        // Loop through cell points and make sure they are not marked for
+        // unrefinement
+        forAll (cPoints, j)
+        {
+            splitPointsToUnrefine[cPoints[j]] = false;
+        }
+    }
+
+
+    // PART 4: Ensure face consistent (2:1 constraint) and possibly point
+    // consistent (4:1 constraint) unrefinement
+
+    // Get necessary mesh data
+    const label nCells = mesh_.nCells();
+    const labelListList& meshPointCells = mesh_.pointCells();
+
+    do
+    {
+        // First, create cells to unrefine (all cells sharing point to unrefine)
+        boolList cellsToUnrefine(nCells, false);
+
+        // Loop through all split points to unrefine
+        forAll (splitPointsToUnrefine, pointI)
+        {
+            if (splitPointsToUnrefine[pointI])
+            {
+                // This split point is marked for unrefinement, collect all of
+                // its cells
+                const labelList& pCells = meshPointCells[pointI];
+                forAll (pCells, i)
+                {
+                    cellsToUnrefine[pCells[i]] = true;
+                }
+            }
+        }
+
+        // Count number of removed cells from unrefinement (cells that will
+        // not be unrefined)
+        label nRemCells = 0;
+
+        // Note: No need to check for point consistent unrefinement if a point
+        // consistent refinement process has been ensured. Skip it.
+//        if (pointBasedRefinement_)
+//        {
+//            // Check for 4:1 point based consistent unrefinement. Updates
+//            // cellsToUnrefine and returns number of removed cells from
+//            // unrefinement in this iteration
+//            nRemCells +=
+//                pointConsistentUnrefinement
+//                (
+//                    splitPointsToUnrefine,
+//                    cellsToUnrefine
+//                );
+//        }
+
+        // Check for 2:1 face based consistent unrefinement. Updates
+        // cellsToUnrefine and returns number of removed cells from unrefinement
+        // in this iteration
+        nRemCells += faceConsistentUnrefinement(cellsToUnrefine);
+
+        // Global reduction
+        reduce(nRemCells, sumOp<label>());
+
+        if (debug)
+        {
+            Info<< "Removed " << nRemCells
+                << " cells to obtain consistent unrefinement."
+                << endl;
+        }
+
+        // If we have removed at least one cell from unrefinement, we need to
+        // protect its split points as well from unrefinement
+        if (nRemCells > 0)
+        {
+            // Get point cells
+            const labelListList& meshPointCells = mesh_.pointCells();
+
+            // Loop through all split points to unrefine
+            forAll (splitPointsToUnrefine, pointI)
+            {
+                if (splitPointsToUnrefine[pointI])
+                {
+                    // This is a split point for unrefinement, get the cells
+                    const labelList& pCells = meshPointCells[pointI];
+
+                    // Loop through all point cells
+                    forAll (pCells, i)
+                    {
+                        if (!cellsToUnrefine[own])
+                        {
+                            // Cell must not be refined, remove point from
+                            // unrefinement as well
+                            splitPointsToUnrefine[pointI] = false;
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+
+    } while (nRemCells > 0);
+
+    // Collect all split points to unrefine in a dynamic list
+    dynamicLabelList splitPointsToUnrefineDynamic(nPoints);
+
+    forAll (splitPointsToUnrefine, pointI)
+    {
+        if (splitPointsToUnrefine[pointI])
+        {
+            // Split point marked for unrefinement, append it
+            splitPointsToUnrefineDynamic.append(pointI);
+        }
+    }
+
+    // Transfer the contents into the data member (ordinary list)
+    splitPointsToUnrefine_.transfer(splitPointsToUnrefineDynamic);
+
+    // Make the refinement history more compact
+    history_.compact();
 }
 
 
-bool Foam::layerAdditionRemoval::changeTopology() const
+bool Foam::polyhedralRefinement::changeTopology() const
 {
-    // Protect from multiple calculation in the same time-step
-    if (triggerRemoval_ > -1 || triggerAddition_ > -1)
+    // If either cellsToRefine_ or splitPointsToUnrefine_ have at least one
+    // entry, the topology is changing
+    if (!cellsToRefine_.empty() || !splitPointsToUnrefine_.empty())
     {
         return true;
     }
 
-    // Go through all the cells in the master layer and calculate
-    // approximate layer thickness as the ratio of the cell volume and
-    // face area in the face zone.
-    // Layer addition:
-    //     When the max thickness exceeds the threshold, trigger refinement.
-    // Layer removal:
-    //     When the min thickness falls below the threshold, trigger removal.
-
-    const faceZone& fz = topoChanger().mesh().faceZones()[faceZoneID_.index()];
-    const labelList& mc = fz.masterCells();
-
-    if (debug)
-    {
-        // Check master cell addressing
-        if (min(mc) < 0)
-        {
-            const polyMesh& mesh = topoChanger().mesh();
-
-            fileName fvPath(mesh.time().path()/"VTK");
-            mkDir(fvPath);
-
-            Info<< "Writing VTK files master face zone"
-                << "Layer addition/removal " << name()
-                << ", face zone " << faceZoneID_.name()
-                << "into " << fvPath
-                << endl;
-
-            primitiveFacePatch::writeVTK
-            (
-                fvPath/fileName(faceZoneID_.name() + "FaceZone"),
-                mesh.faceZones()[faceZoneID_.index()]().localFaces(),
-                mesh.faceZones()[faceZoneID_.index()]().localPoints()
-            );
-
-            primitiveFacePatch::writeVTKNormals
-            (
-                fvPath/fileName(faceZoneID_.name() + "FaceZoneNormals"),
-                mesh.faceZones()[faceZoneID_.index()]().localFaces(),
-                mesh.faceZones()[faceZoneID_.index()]().localPoints()
-            );
-
-            FatalErrorIn("bool layerAdditionRemoval::changeTopology() const")
-                << "Error in master cell addressing for face zone "
-                    << faceZoneID_.name()
-                    << abort(FatalError);
-        }
-    }
-
-    const scalarField& V = topoChanger().mesh().cellVolumes();
-    const vectorField& S = topoChanger().mesh().faceAreas();
-
-    if (min(V) < -VSMALL)
-    {
-        FatalErrorIn("bool layerAdditionRemoval::changeTopology() const")
-            << "negative cell volume. Error in mesh motion before "
-            << "topological change."
-            << abort(FatalError);
-    }
-
-    scalar minDelta = GREAT;
-    scalar maxDelta = 0;
-    scalar avgDelta = 0;
-    scalar nAvg = 0;
-
-    if (fz.size())
-    {
-        forAll (fz, faceI)
-        {
-            scalar curDelta = V[mc[faceI]]/mag(S[fz[faceI]]);
-            avgDelta += curDelta;
-            minDelta = min(minDelta, curDelta);
-            maxDelta = max(maxDelta, curDelta);
-        }
-
-        nAvg += fz.size();
-    }
-
-    // If the patch is empty on a processor in a parallel simulation,
-    // original values will be preserved.  HJ, 7/Mar/2011
-
-    reduce(minDelta, minOp<scalar>());
-    reduce(maxDelta, maxOp<scalar>());
-    reduce(avgDelta, sumOp<scalar>());
-    reduce(nAvg, sumOp<scalar>());
-
-    avgDelta /= nAvg;
-
-    if (debug)
-    {
-        Info<< "bool layerAdditionRemoval::changeTopology() const "
-            << " for object " << name() << " : " << nl
-            << "Layer thickness: min: " << minDelta
-            << " max: " << maxDelta << " avg: " << avgDelta
-            << " old thickness: " << oldLayerThickness_ << nl
-            << "Removal threshold: " << minLayerThickness_
-            << " addition threshold: " << maxLayerThickness_ << endl;
-    }
-
-    bool topologicalChange = false;
-
-    // If the thickness is decreasing and crosses the min thickness,
-    // trigger removal
-    if (oldLayerThickness_ < 0)
-    {
-        if (debug)
-        {
-            Info<< "First step. No addition/removal" << endl;
-        }
-
-        // No topological changes allowed before first mesh motion
-        // HJ, 8/Oct/2002
-        oldLayerThickness_ = avgDelta;
-
-        topologicalChange = false;
-    }
-    else if (avgDelta < oldLayerThickness_)
-    {
-        // Layers moving towards removal
-        if (minDelta < minLayerThickness_)
-        {
-            // Check layer pairing
-            if (setLayerPairing())
-            {
-                // A mesh layer detected.  Check that collapse is valid
-                if (validCollapse())
-                {
-                    // At this point, info about moving the old mesh
-                    // in a way to collapse the cells in the removed
-                    // layer is available.  Not sure what to do with
-                    // it.  HJ, 3/Nov/2003
-
-                    if (debug)
-                    {
-                        Info<< "bool layerAdditionRemoval::changeTopology() "
-                            << " const for object " << name() << " : "
-                            << "Triggering layer removal" << endl;
-                    }
-
-                    triggerRemoval_ = topoChanger().morphIndex();
-
-                    // Old thickness looses meaning.
-                    // Set it up to indicate layer removal
-                    oldLayerThickness_ = GREAT;
-
-                    topologicalChange = true;
-                }
-                else
-                {
-                    // No removal, clear addressing
-                    clearAddressing();
-                }
-            }
-        }
-        else
-        {
-            oldLayerThickness_ = avgDelta;
-        }
-    }
-    else
-    {
-        // Layers moving towards addition
-        if (maxDelta > maxLayerThickness_)
-        {
-            if (debug)
-            {
-                Info<< "bool layerAdditionRemoval::changeTopology() const "
-                    << " for object " << name() << " : "
-                    << "Triggering layer addition" << endl;
-            }
-
-            triggerAddition_ = topoChanger().morphIndex();
-
-            // Old thickness looses meaning.
-            // Set it up to indicate layer removal
-            oldLayerThickness_ = 0;
-
-            topologicalChange = true;
-        }
-        else
-        {
-            oldLayerThickness_ = avgDelta;
-        }
-    }
-
-    return topologicalChange;
+    // Otherwise, there is no topo change (at least on this processor in this
+    // time step): return false
+    return false;
 }
 
 
-void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
+void Foam::polyhedralRefinement::setRefinement(polyTopoChange& ref) const
 {
-    // Insert the layer addition/removal instructions
-    // into the topological change
-
-    if (triggerRemoval_ == topoChanger().morphIndex())
+    if (!cellsToRefine_.empty())
     {
-        removeCellLayer(ref);
-
-        // Clear addressing.  This also resets the addition/removal data
-        if (debug)
-        {
-            Info<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
-                << " for object " << name() << " : "
-                << "Clearing addressing after layer removal. " << endl;
-        }
-
-        triggerRemoval_ = -1;
-        clearAddressing();
+        // There are some cells to refine, insert the refinement instruction and
+        // clear the list of cells
+        setPolyhedralRefinement(ref);
+        cellsToRefine_.clear();
     }
 
-    if (triggerAddition_ == topoChanger().morphIndex())
+    if (!splitPointsToUnrefine_.empty())
     {
-        addCellLayer(ref);
-
-        // Clear addressing.  This also resets the addition/removal data
-        if (debug)
-        {
-            Info<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
-                << " for object " << name() << " : "
-                << "Clearing addressing after layer addition. " << endl;
-        }
-
-        triggerAddition_ = -1;
-        clearAddressing();
+        // There are some split points to unrefine around, inser the
+        // unrefinement instruction and clear the list of split points
+        setPolyhedralUnrefinement(ref);
+        splitPointsToUnrefine_.clear();
     }
 }
 
 
-void Foam::layerAdditionRemoval::updateMesh(const mapPolyMesh&)
+void Foam::polyhedralRefinement::updateMesh(const mapPolyMesh&)
 {
     if (debug)
     {
-        Info<< "layerAdditionRemoval::updateMesh(const mapPolyMesh&) "
+        Info<< "polyhedralRefinement::updateMesh(const mapPolyMesh&) "
             << " for object " << name() << " : "
             << "Clearing addressing on external request. ";
 
@@ -491,7 +3715,7 @@ void Foam::layerAdditionRemoval::updateMesh(const mapPolyMesh&)
 }
 
 
-void Foam::layerAdditionRemoval::setMinLayerThickness(const scalar t) const
+void Foam::polyhedralRefinement::setMinLayerThickness(const scalar t) const
 {
     if
     (
@@ -501,7 +3725,7 @@ void Foam::layerAdditionRemoval::setMinLayerThickness(const scalar t) const
     {
         FatalErrorIn
         (
-            "void layerAdditionRemoval::setMinLayerThickness("
+            "void polyhedralRefinement::setMinLayerThickness("
             "const scalar t) const"
         )   << "Incorrect layer thickness definition."
             << abort(FatalError);
@@ -511,13 +3735,13 @@ void Foam::layerAdditionRemoval::setMinLayerThickness(const scalar t) const
 }
 
 
-void Foam::layerAdditionRemoval::setMaxLayerThickness(const scalar t) const
+void Foam::polyhedralRefinement::setMaxLayerThickness(const scalar t) const
 {
     if (t < minLayerThickness_)
     {
         FatalErrorIn
         (
-            "void layerAdditionRemoval::setMaxLayerThickness("
+            "void polyhedralRefinement::setMaxLayerThickness("
             "const scalar t) const"
         )   << "Incorrect layer thickness definition."
             << abort(FatalError);
@@ -527,7 +3751,7 @@ void Foam::layerAdditionRemoval::setMaxLayerThickness(const scalar t) const
 }
 
 
-void Foam::layerAdditionRemoval::write(Ostream& os) const
+void Foam::polyhedralRefinement::write(Ostream& os) const
 {
     os  << nl << type() << nl
         << name()<< nl
@@ -538,7 +3762,7 @@ void Foam::layerAdditionRemoval::write(Ostream& os) const
 }
 
 
-void Foam::layerAdditionRemoval::writeDict(Ostream& os) const
+void Foam::polyhedralRefinement::writeDict(Ostream& os) const
 {
     os  << nl << name() << nl << token::BEGIN_BLOCK << nl
         << "    type " << type()
diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H
index e777eeccf..61443fdcc 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H
@@ -55,7 +55,7 @@ Notes
 #include "polyMeshModifier.H"
 #include "labelIOList.H"
 #include "polyRefinementHistory.H"
-#include "removeFaces.H"
+#include "removePolyFaces.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -72,26 +72,41 @@ class polyhedralRefinement
 {
     // Private data
 
-        //- Cell refinement level
-        labelIOList cellLevel_;
+        //- Reference to polyMesh for easy access in helper functions
+        polyMesh& mesh_;
 
-        //- Point refinement level
-        labelIOList pointLevel_;
 
-        //- Typical edge length between unrefined points
-        const scalar level0EdgeLength_;
+        // Refinement control and handling
 
-        //- Refinement history
-        polyRefinementHistory history_;
+            //- List of cells to refine in this time step
+            labelList cellsToRefine_;
 
-        //- Face remover engine
-        removeFaces faceRemover_;
+            //- List of split point labels to unrefine in this time step
+            labelList splitPointsToUnrefine_;
 
-        //- List of cells to refine in this time step
-        labelList cellsToRefine_;
+            //- Cell refinement level
+            labelIOList cellLevel_;
 
-        //- List of split point labels to unrefine in this time step
-        labelList splitPointLabels_;
+            //- Point refinement level
+            labelIOList pointLevel_;
+
+            //- Typical edge length between unrefined points
+            const scalar level0EdgeLength_;
+
+            //- Refinement history
+            polyRefinementHistory history_;
+
+            //- Face remover engine
+            removePolyFaces faceRemover_;
+
+            //- Switch whether to use point based consistency on refinement
+            Switch pointBasedRefinement_;
+
+            //- Number of buffer layers for refinement
+            label nBufferLayers_;
+
+            //- Number of refinement/unrefinement iterations done so far
+            label nRefinementIterations_;
 
 
     // Private Member Functions
@@ -112,13 +127,17 @@ class polyhedralRefinement
             //  points smaller than the level
             label getAnchorLevel(const label faceI) const;
 
-            //- Get points of a cell by going through its faces (without using
-            //  cellPoints addressing)
-            Xfer<labelList> cellPoints(const label cellI) const;
-
             //- Calculate level0EdgeLength_ (constructor helper)
             void calcLevel0EdgeLength() const;
 
+            //- Set file instance for cellLevel_, pointLevel_ and history_
+            void setInstance(const fileName& inst);
+
+            //- Extend cells to refine. Given a markup field of refinement
+            //  candidates (true for cells to refine), marks the neigbhours of
+            //  marked cells as well
+            void extendMarkedCells(boolList& refineCell) const;
+
 
         // Global topology modification functions (operate on whole polyMesh)
         // and directly use local topology modification functions below
@@ -136,13 +155,13 @@ class polyhedralRefinement
             //   that will allocate interpolation.
             //  -points added to split edge: added from edge start()
             //  -midpoints added: added from cellPoints[0].
-            void setPolyhedralRefinement() const;
+            void setPolyhedralRefinement(polyTopoChange& ref);
 
             //- Remove some of the previously performed refinement. Uses
             //  splitPointLabels_ to determine the unrefinement.
             //  All n pointCells of a split point will be combined into
             //  the lowest numbered cell of those n.
-            void setUnrefinement() const;
+            void setPolyhedralUnrefinement(polyTopoChange& ref);
 
 
         // Local topology modification functions (operate on cells/faces)
@@ -150,7 +169,7 @@ class polyhedralRefinement
             //- Adds a face on top of existing faceI. Reverses if nessecary
             label addFace
             (
-                polyTopoChange& meshMod,
+                polyTopoChange& ref,
                 const label faceI,
                 const face& newFace,
                 const label own,
@@ -160,7 +179,7 @@ class polyhedralRefinement
             //- Adds internal face from point. No checks on reversal
             label addInternalFace
             (
-                polyTopoChange& meshMod,
+                polyTopoChange& ref,
                 const label meshFaceI,
                 const label meshPointI,
                 const face& newFace,
@@ -172,7 +191,7 @@ class polyhedralRefinement
             //  points. Reverses if nessecary
             void modifyFace
             (
-                polyTopoChange& meshMod,
+                polyTopoChange& ref,
                 const label faceI,
                 const face& newFace,
                 const label own,
@@ -190,7 +209,7 @@ class polyhedralRefinement
                 const labelList& faceAnchorLevel,
                 const labelList& edgeMidPoint,
                 const label cellI,
-                polyTopoChange& meshMod
+                polyTopoChange& ref
             ) const;
 
 
@@ -243,7 +262,7 @@ class polyhedralRefinement
 
             //- Get index of point with minimum or maximum point level
             template<class BinaryOp>
-            label findMinMaxLevel(const labelList& f) const;
+            label findMinMaxLevel(const labelList& f, BinaryOp& op) const;
 
             //- Count number of vertices <= anchorLevel for a given face
             label countAnchors
@@ -278,7 +297,7 @@ class polyhedralRefinement
 
                 Map<edge>& midPointToAnchors,
                 Map<edge>& midPointToFaceMids,
-                polyTopoChange& meshMod
+                polyTopoChange& ref
             ) const;
 
             //- If p0 and p1 are existing vertices check if edge is split and insert
@@ -297,7 +316,7 @@ class polyhedralRefinement
             //- Check orientation of added internal face
             void checkInternalOrientation
             (
-                polyTopoChange& meshMod,
+                polyTopoChange& ref,
                 const label cellI,
                 const label faceI,
                 const point& ownPt,
@@ -308,7 +327,7 @@ class polyhedralRefinement
             //- Check orientation of a new boundary face
             void checkBoundaryOrientation
             (
-                polyTopoChange& meshMod,
+                polyTopoChange& ref,
                 const label cellI,
                 const label faceI,
                 const point& ownPt,
@@ -319,45 +338,24 @@ class polyhedralRefinement
 
         // Refinement/unrefinement consistency checks
 
-            //- Given valid mesh, current cell level and proposed
-            //  cells to refine calculate any clashes (due to 2:1) and return
-            //  ok list of cells to refine.
-            Xfer<labelList> consistentRefinement
-            (
-                const labelList& cellsToRefine,
-                const bool pointBasedRefinement
-            ) const;
-
-            //- Given proposed splitPoints to unrefine calculate
-            //  any clashes (due to 2:1) and return ok list of points to
-            //  unrefine.
-            Xfer<labelList> consistentUnrefinement
-            (
-                const labelList& pointsToUnrefine,
-                const bool pointBasedUnrefinement
-            ) const;
-
-            //- Updates refineCell such that a face consistent 2:1 refinement
+            //- Updates cellsToRefine such that a face consistent 2:1 refinement
             //  is obtained. Returns local number of cells changed
-            label faceConsistentRefinement(PackedBoolList& refineCell) const;
+            label faceConsistentRefinement(boolList& cellsToRefine) const;
 
-            //- Updates refineCell such that a point consistent 4:1 refinement
+            //- Updates cellsToRefine such that a point consistent 4:1 refinement
             //  is obtained. Returns local number of cells changed
-            label pointConsistentRefinement(PackedBoolList& refineCell) const;
+            label pointConsistentRefinement(boolList& cellsToRefine) const;
 
-            //- Updates unrefineCell such that a face consistent 2:1
-            //  unrefinement is obtained. Returns local number of cells changed
-            label faceConsistentUnrefinement
-            (
-                PackedBoolList& unrefineCell
-            ) const;
+            //- Updates cellsToUnrefine such that a face consistent 2:1
+            //  unrefinement is obtained. Returns local number of points changed
+            label faceConsistentUnrefinement(boolList& cellsToUnrefine) const;
 
-            //- Updates unrefineCell such that a point consistent 4:1
-            //  unrefinement is obtained. Returns local number of cells changed
+            //- Updates cellsToUnrefine such that a point consistent 4:1
+            //  unrefinement is obtained. Returns local number of points changed
             label pointConsistentUnrefinement
             (
-                const PackedBoolList& unrefinePoints,
-                PackedBoolList& unrefineCell
+                const boolList& splitPointsToUnrefine,
+                boolList& cellsToUnrefine
             ) const;
 
 
@@ -421,14 +419,26 @@ public:
 
         // Edit
 
-            //- Set cells to refine
-            void setCellsToRefine(const labelList& cellsToRefine);
+            //- Set cells to refine given a list of refinement
+            //  candidates. Refinement candidates are extended within the
+            //  function due to possible 4:1 conflicts and specified number of
+            //  buffer layers.
+            //  Note: must be called BEFORE setSplitPointsToUnrefine
+            void setCellsToRefine(const labelList& refinementCellCandidates);
 
-            //- Update split points to unrefine
-            void updateSplitPointsToUnrefine();
+            //- Set split points to unrefine given a list of all mesh points
+            //  that are candidates for unrefinement. Split points are
+            //  determined as a subset of unrefinement candidates, avoiding
+            //  splitting points of cells that are going to be refined at the
+            //  same time and ensuring consistent unrefinement.
+            //  Note: must be called AFTER setCellsToRefine
+            void setSplitPointsToUnrefine
+            (
+                const labelList& unrefinementPointCandidates
+            );
 
 
-        // polyMeshModifier interface
+        // Inherited interface from polyMeshModifier
 
             //- Check for topology change
             virtual bool changeTopology() const;