From 8ceec6456a57f93ab76f2e25608be8c07b2a922f Mon Sep 17 00:00:00 2001
From: Mattijs Janssens <m.janssens@opencfd.co.uk>
Date: Thu, 15 Nov 2018 14:58:56 +0000
Subject: [PATCH] Mapping injected faces through interpolation

---
 README.txt                                    |  47 +
 src/dynamicFvMesh/Make/files                  |   1 +
 src/dynamicFvMesh/Make/options                |   3 +
 .../balanceParDict                            |  29 +
 .../dynamicMeshDict                           | 125 +++
 .../dynamicRefineBalancedFvMesh.C             | 843 ++++++++++++++++++
 .../dynamicRefineBalancedFvMesh.H             | 206 +++++
 .../dynamicRefineBalancedFvMeshTemplates.C    | 217 +++++
 .../dynamicRefineFvMesh/dynamicRefineFvMesh.C | 309 ++++---
 .../dynamicRefineFvMesh/dynamicRefineFvMesh.H |  44 +-
 .../dynamicRefineFvMeshTemplates.C            | 196 ++++
 .../polyTopoChange/hexRef8/hexRef8.C          |   2 +-
 src/finiteVolume/fvMesh/fvMesh.C              |  10 +
 .../constant/dynamicMeshDict                  |   1 +
 14 files changed, 1906 insertions(+), 127 deletions(-)
 create mode 100644 README.txt
 create mode 100644 src/dynamicFvMesh/dynamicRefineBalancedFvMesh/balanceParDict
 create mode 100644 src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicMeshDict
 create mode 100644 src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C
 create mode 100644 src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.H
 create mode 100644 src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMeshTemplates.C
 create mode 100644 src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMeshTemplates.C

diff --git a/README.txt b/README.txt
new file mode 100644
index 0000000000..273b4dad20
--- /dev/null
+++ b/README.txt
@@ -0,0 +1,47 @@
+Notes from merging Integration-TUD-addOns and plus.
+
+0)
+- Changed
+dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
+to created internal faces out-of-nothing.
+- Added the mapNewInternalFaces to dynamicRefineFvMesh
+- unset FOAM_SETNAN, FOAM_SIGFPE and run
+
+    testCases/testAMRandLoadBalancing/damBreakWithObstacle
+
+- However when writing V0 we noticed a nan in the V0 field.
+  (this doesn't get written anymore but it might indicate a bug)
+
+
+1) Tried moving additional mapping of surface fields to
+
+    virtual dynamicRefineFvMesh::mapFields(const mapPolyMesh&);
+
+so it would get called from the fvMesh::updateMesh. However
+this mapping (explicitly) gets done using unadapted addressing
+and it causes the cell addressing to be wrong:
+
+
+[1] --> FOAM FATAL ERROR:
+[1] index 826015320 out of range 0 ... 71079
+[1]
+[1]     From function void Foam::UList<T>::checkIndex(Foam::label) const [with T = Foam::cell; Foam::label = int]
+[1]     in file /home/preston2/mattijs/OpenFOAM/work/OpenFOAM-plus.integration-TUD/src/OpenFOAM/lnInclude/UListI.H at line 106.
+
+
+From dynamicRefineFvMesh::mapNewInternalFaces : (I think)
+
+                    cell faceOwner = this->cells()[owner[facei]];
+
+So this adaptation should be done in a separate pass
+
+
+2) Mapping new internal faces:
+- currently done by averaging the values of 'properly' mapped
+  faces from owner and neighbour.
+- this would not work if a cell gets split into 3x3x3
+- instead this should be done by some geometric interpolation
+  from point values?
+- have selection mechanism instead or use oriented flag on surfaceFields
+
+3) Current averaging is valid only for internal faces & only scalar/vector
diff --git a/src/dynamicFvMesh/Make/files b/src/dynamicFvMesh/Make/files
index b8432c07c9..bb04e8ac64 100644
--- a/src/dynamicFvMesh/Make/files
+++ b/src/dynamicFvMesh/Make/files
@@ -5,6 +5,7 @@ dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C
 dynamicMultiMotionSolverFvMesh/dynamicMultiMotionSolverFvMesh.C
 dynamicInkJetFvMesh/dynamicInkJetFvMesh.C
 dynamicRefineFvMesh/dynamicRefineFvMesh.C
+dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C
 dynamicMotionSolverListFvMesh/dynamicMotionSolverListFvMesh.C
 
 simplifiedDynamicFvMesh/simplifiedDynamicFvMeshes.C
diff --git a/src/dynamicFvMesh/Make/options b/src/dynamicFvMesh/Make/options
index 3540c13bfc..718fde7d5c 100644
--- a/src/dynamicFvMesh/Make/options
+++ b/src/dynamicFvMesh/Make/options
@@ -1,10 +1,13 @@
 EXE_INC = \
+    -g -DFULLDEBUG -O0 -g \
+    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
+    -ldecompositionMethods \
     -lmeshTools \
     -ldynamicMesh \
     -lfiniteVolume
diff --git a/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/balanceParDict b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/balanceParDict
new file mode 100644
index 0000000000..08d1ef6538
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/balanceParDict
@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.1                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      balanceParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 2;
+
+method          ptscotch;
+
+constraints
+{
+    refinementHistory
+    {
+        type    refinementHistory;
+    }
+}
+// ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicMeshDict b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicMeshDict
new file mode 100644
index 0000000000..a8df1cb7ef
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicMeshDict
@@ -0,0 +1,125 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   dynamicRefineBalancedFvMesh;
+
+refinementControls
+{
+    enableRefinementControl  true;
+    
+    fields // must be scalarFields
+    (
+        //alpha (min max refineLevel)
+        alpha (0.01 0.99 2) // refine cells where alpha in [0.01:0.99] with maximal 2 refinement layers
+    );
+    
+    interface // must be a scalarField (only one dictionary!)
+    (
+        alpha // refine interface (found based on snGrad of alpha > 0.1) 
+        {
+            innerRefLayers 2; // describes how many cell layers inside phase alpha are to be refined
+            outerRefLayers 5; // describes how many cell layers outside phase alpha are to be refined
+            
+            // optional settings:
+            maxRefineLevel 4; // max refinement layers; Default: maxRefinement from dynamicRefineFvMeshCoeffs is used
+            // to get slower than 2:1 refinement; add #nAddLayers between each refinement level at the interface
+            nAddLayers 1; //Default: 0 
+        }
+    );
+    
+    gradients // must be scalars
+    (
+        // arguments as in 'fields'
+        // min/max values are based on mag(fvc::grad(volScalarField)) * cellVolume
+        T    (0.01 10 1)
+    );
+    
+    curls // must be vectors
+    (
+        // arguments as in 'fields'
+        // min/max values are based on mag(fvc::curl(volVectorField))
+        U   (0.5 1 2)
+    );
+    
+    regions
+    (
+        boxToCell
+        {
+            minLevel 1;
+            
+            box (-1 0.001 0.002)(1 0.005 0.003);
+        }
+        
+    );
+}
+
+dynamicRefineFvMeshCoeffs
+{
+    // Extra entries for balancing
+    enableBalancing true;
+    allowableImbalance 0.15;
+
+    // How often to refine
+    refineInterval  10;
+    
+    // Field to be refinement on (set it to 'internalRefinementField' to use the
+    // refinementControls dictionary entries above)
+    field           internalRefinementField;
+    
+    // Refine field inbetween lower..upper
+    lowerRefineLevel 0.5; // do not change
+    upperRefineLevel 3.5; // maxRefinement+0.5
+    
+    // If value < unrefineLevel unrefine
+    unrefineLevel   -0.5; // do not change
+    
+    // Have slower than 2:1 refinement
+    nBufferLayers   4;
+    
+    // Refine cells only up to maxRefinement levels
+    maxRefinement   3;
+    
+    // Stop refinement if maxCells reached
+    maxCells        200000;
+    
+    // Flux field and corresponding velocity field. Fluxes on changed
+    // faces get recalculated by interpolating the velocity. Use 'none'
+    // on surfaceScalarFields that do not need to be reinterpolated.
+    correctFluxes
+    (
+        (phi Urel)
+        (phiAbs U)
+        (phiAbs_0 U_0)
+        (nHatf none)
+        (rho*phi none)
+        (ghf none)
+    );
+
+    // List of non-flux surface<Type>Fields to be mapped
+    // only for new internal faces (AMR refine)
+    mapSurfaceFields
+    (
+        Uf
+        Uf_0
+    );
+    
+    // Write the refinement level as a volScalarField
+    dumpLevel       true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C
new file mode 100644
index 0000000000..d5098ffe68
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C
@@ -0,0 +1,843 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 Tyler Voskuilen
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is a derivative work of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+
+\*---------------------------------------------------------------------------*/
+
+#include "dynamicRefineBalancedFvMesh.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceInterpolate.H"
+#include "volFields.H"
+#include "polyTopoChange.H"
+#include "surfaceFields.H"
+#include "syncTools.H"
+#include "pointFields.H"
+#include "fvCFD.H"
+#include "volPointInterpolation.H"
+#include "pointMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(dynamicRefineBalancedFvMesh, 0);
+    addToRunTimeSelectionTable(dynamicFvMesh, dynamicRefineBalancedFvMesh, IOobject);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+Foam::label Foam::dynamicRefineBalancedFvMesh::topParentID(label p)
+{
+    label nextP = meshCutter().history().splitCells()[p].parent_;
+    if( nextP < 0 )
+    {
+        return p;
+    }
+    else
+    {
+        return topParentID(nextP);
+    }
+}
+
+Foam::List<Foam::scalar> Foam::dynamicRefineBalancedFvMesh::readRefinementPoints()
+{
+    dictionary refineDict
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "dynamicMeshDict",
+                time().constant(),
+                *this,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).subDict("dynamicRefineFvMeshCoeffs")
+    );
+    
+    List<scalar> refData(4, scalar(0));
+    
+    refData[0] = readScalar(refineDict.lookup("unrefineLevel"));
+    refData[1] = readScalar(refineDict.lookup("lowerRefineLevel"));
+    refData[2] = readScalar(refineDict.lookup("refineInterval"));
+    refData[3] = readScalar(refineDict.lookup("maxRefinement"));
+        
+    return refData;
+}
+
+void Foam::dynamicRefineBalancedFvMesh::updateRefinementField()
+{
+    Info<< "Calculating internal refinement field" << endl;
+    
+    volScalarField& intRefFld = *internalRefinementFieldPtr_;
+    volScalarField& targetFld = *targetLevelPtr_;
+    volScalarField& currFld = *isLevelPtr_;
+    
+    // Set the internal refinement field to zero to start with
+    intRefFld = dimensionedScalar("zero",dimless,0.0);
+    
+    // Get the cell level field from dynamicRefineFvMesh
+    const labelList& cellLevel = meshCutter().cellLevel();
+    
+    // Read the points at which refinement and unrefinement occur from the
+    // dynamicMeshDict entries
+    List<scalar> refinePoints = readRefinementPoints();
+    
+    // init list for refine/unrefine field (-1=unrefine, 1=refine, 0=do nothing)
+    labelList markRefineCells (this->nCells(), 0);
+    
+    // init list for target refinement level per cell
+    labelList targetLevel (this->nCells(), 0);
+    
+    // First fields
+    List<word> fieldNames = fields_.toc();
+    Field<scalar> refFld(nCells(),0.0);
+    
+    forAll(fieldNames, i)
+    {
+        word fldName = fieldNames[i];
+        scalar minValue = fields_[fldName][0];
+        scalar maxValue = fields_[fldName][1];
+        label refineLevel = static_cast<label>(fields_[fldName][2]);
+        
+        const volScalarField& fld = this->lookupObject<volScalarField>(fldName);
+        
+        // Limit the value of refFld based on its max level
+        forAll(fld, cellI)
+        {
+            if ((fld[cellI] >= minValue) && (fld[cellI] <= maxValue))
+            {
+                // increase targetLevel up to refineLevel 
+                // BUT: do not decrease if cell already marked for higher refinement level by previous criterion
+                targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
+            }
+        }
+    }
+    
+    // Then gradients
+    List<word> gradFieldNames = gradFields_.toc();
+    
+    Field<scalar> cubeRtV = Foam::pow(this->V(),1.0/3.0);
+    
+    forAll(gradFieldNames, i)
+    {
+        word fldName = gradFieldNames[i];
+        scalar minValue = gradFields_[fldName][0];
+        scalar maxValue = gradFields_[fldName][1];
+        label refineLevel = static_cast<label>(gradFields_[fldName][2]);
+        
+        const volScalarField& fld = this->lookupObject<volScalarField>(fldName);
+
+        refFld = mag(fvc::grad(fld)) * cubeRtV;
+        
+        // Limit the value of refFld based on its max level
+        forAll(refFld, cellI)
+        {
+            if ((refFld[cellI] >= minValue) && (refFld[cellI] <= maxValue))
+            {
+                targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
+            }
+        }
+    }
+    
+    // Then curls
+    List<word> curlFieldNames = curlFields_.toc();
+    
+    
+    forAll(curlFieldNames, i)
+    {
+        word fldName = curlFieldNames[i];
+        scalar minValue = curlFields_[fldName][0];
+        scalar maxValue = curlFields_[fldName][1];
+        label refineLevel = static_cast<label>(curlFields_[fldName][2]);
+        
+        const volVectorField& fld = this->lookupObject<volVectorField>(fldName);
+        
+        refFld = mag(fvc::curl(fld));
+        
+        // Limit the value of refFld based on its max level
+        forAll(refFld, cellI)
+        {
+            if ((refFld[cellI] >= minValue) && (refFld[cellI] <= maxValue))
+            {
+                targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
+            }
+        }
+    }
+    
+    // At the interface, assumed to be always the maximum refinement level
+    List<word> interfaceRefineField = interface_.toc();
+    forAll(interfaceRefineField, i)
+    {
+        word fldName = interfaceRefineField[i];
+        // read region of maximum refinement levels inside and outside of interface indicator field 
+        // (does not need to be alpha, can also be a concentration field)
+        
+        scalar innerRefLayers = readScalar(interface_[fldName].lookup("innerRefLayers"));
+        scalar outerRefLayers = readScalar(interface_[fldName].lookup("outerRefLayers"));
+        
+        scalar nAddLayers(0);
+        if (interface_[fldName].found("nAddLayers"))
+        {
+            nAddLayers = readScalar(interface_[fldName].lookup("nAddLayers"));
+        }
+        
+        label refineLevel = refinePoints[3];
+        
+        if (interface_[fldName].found("maxRefineLevel"))
+        {
+            refineLevel = readScalar(interface_[fldName].lookup("maxRefineLevel"));
+
+            // to avoid user input mistakes, limit the value with the maximum allowed
+            refineLevel = min(refinePoints[3], refineLevel);
+        }
+        
+        const volScalarField& fld = this->lookupObject<volScalarField>(fldName);
+        
+        volScalarField isInterface(intRefFld * 0.0);
+        isInterface = dimensionedScalar("isInterface",dimless,0.0);
+        
+        surfaceScalarField deltaAlpha = mag(fvc::snGrad(fld) / this->deltaCoeffs());
+        
+        const labelUList& owner = this->owner();
+        const labelUList& neighbour = this->neighbour();
+        
+        forAll(deltaAlpha, faceI)
+        {
+        //TODO: add a min max value of the specified field to be marked for refinement
+        // NOW: hard-coded method snGradAlpha: checks mag(alpha_N - alpha_P) > 0.1 ? 1:0
+            if (deltaAlpha[faceI] > 0.1) // currently a fixed prescribed value; should be read from díctionary
+            {
+                label own = owner[faceI];
+                label nei = neighbour[faceI];
+                
+                // set isInterface field to one
+                isInterface[own] = 1.0;
+                isInterface[nei] = 1.0;
+            }
+        }
+        
+        // assumed fld=0.5*(fldMax+fldMin) defines the interface
+        dimensionedScalar fldInterfaceValue(0.5*(max(fld)+min(fld)));
+        
+        //-DD: old implementation based on face interpolation
+        //     which results in slower transport in diagonal direction
+        // add inner refinement layers
+        //for(label i=0; i < innerRefLayers; i++)
+        //{
+        //    isInterface += neg(- fvc::average(fvc::interpolate(isInterface)) * pos(fld - fldInterfaceValue)); 
+        //    isInterface = neg(- isInterface);
+        //}
+        //
+        // add outer refinement layers
+        //for(label i=0; i < outerRefLayers; i++)
+        //{
+        //    isInterface += neg(- fvc::average(fvc::interpolate(isInterface)) * pos(fldInterfaceValue - fld));
+        //    isInterface = neg(- isInterface);
+        //}
+        //
+        //forAll(isInterface, cellI)
+        //{
+        //    if (isInterface[cellI] > 0.5)
+        //    {
+        //        targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
+        //    }
+        //}
+
+        //-DD: new version using volPointInterpolation (direction independent buffer layer)
+        const volPointInterpolation& pInterp = volPointInterpolation::New(*this);
+        const fvMesh& mesh = fld.mesh();
+     
+        // add inner refinement layers
+        for(label i=0; i < innerRefLayers; i++)
+        {
+            volScalarField markInner(isInterface*pos(fld - fldInterfaceValue));
+            pointScalarField markLayerP(pInterp.interpolate(markInner));
+
+            forAll(mesh.C(), cellI)
+            {
+                scalar sum = 0.;
+                label nPoints = 0;
+
+                forAll(mesh.cellPoints()[cellI], pointI)
+                {
+                    sum += markLayerP[mesh.cellPoints()[cellI][pointI]];
+                    nPoints++;
+                }
+                if (nPoints > 0)
+                {
+                    sum /= nPoints;
+                }
+                isInterface[cellI] += sum;
+            }
+        }
+        isInterface = pos(isInterface - SMALL);
+
+        // add outer refinement layers
+        for(label i=0; i < outerRefLayers; i++)
+        {
+            volScalarField markOuter(isInterface*pos(fldInterfaceValue - fld));
+            pointScalarField markLayerP(pInterp.interpolate(markOuter));
+
+            forAll(mesh.C(), cellI)
+            {
+                scalar sum = 0.;
+                label nPoints = 0;
+
+                forAll(mesh.cellPoints()[cellI], pointI)
+                {
+                    sum += markLayerP[mesh.cellPoints()[cellI][pointI]];
+                    nPoints++;
+                }
+                if (nPoints > 0)
+                {
+                    sum /= nPoints;
+                }
+                isInterface[cellI] += sum;
+            }
+        }
+        isInterface = pos(isInterface - SMALL);
+
+        forAll(isInterface, cellI)
+        {
+            if (isInterface[cellI] > 0.5)
+            {
+                targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
+            }
+        }
+        
+        //-DD: old implementation based on face interpolation
+        //     which results in slower transport in diagonal direction
+        // expand additional layers if specified:
+        if (nAddLayers > 0)
+        {
+           // loop over additional layers
+           for(label i=1; i < refineLevel; i++)
+           {
+               // #nAddLayers cells per additional level
+               for(label j=0; j < nAddLayers; j++)
+               {
+                   isInterface += neg(- fvc::average(fvc::interpolate(isInterface)));
+                   isInterface = neg(- isInterface);
+               }
+               
+               forAll(isInterface, cellI)
+               {
+                   if (isInterface[cellI] == 1.0)
+                   {
+                       targetLevel[cellI] = max(targetLevel[cellI], (refineLevel - i));
+                   }
+               }
+           }
+        }
+    }
+    
+    // The set refinement physical regions (force the mesh to stay refined
+    // near key features)
+    forAll(refinedRegions_, regionI)
+    {
+        const entry& region = refinedRegions_[regionI];
+        
+        autoPtr<topoSetSource> source = 
+            topoSetSource::New(region.keyword(), *this, region.dict());
+        
+        cellSet selectedCellSet
+        (
+            *this,
+            "cellSet",
+            nCells()/10+1 //estimate
+        );
+        
+        source->applyToSet
+        (
+            topoSetSource::NEW,
+            selectedCellSet
+        );
+        
+        const labelList cells = selectedCellSet.toc();
+        
+        label minLevel = readLabel(region.dict().lookup("minLevel"));
+        
+        forAll(cells, i)
+        {
+            const label& cellI = cells[i];
+            
+            targetLevel[cellI] = max(targetLevel[cellI], minLevel);
+        }
+    }
+    
+    //-DD: add buffer layer based on targetLevel field to prevent 2-to-1 refinement
+    
+    volScalarField blockedLevel = targetFld * 0.;
+    
+    for (label currLayer=refinePoints[3]; currLayer>=1; currLayer--)
+    {
+        forAll (targetLevel, cellI)
+        {
+            if (targetLevel[cellI] >= currLayer)
+            {
+                blockedLevel[cellI] = targetLevel[cellI];
+            }
+        }
+        
+        for (label i=0; i<nBufferLayers_; i++)
+        {
+            blockedLevel = max(blockedLevel, pos(fvc::average(fvc::interpolate(blockedLevel)) - SMALL)*(currLayer-1));
+        }
+
+        labelList blockLev(blockedLevel.internalField().size(), 0);
+        forAll (blockLev, cellI)
+        {
+            blockLev[cellI] = blockedLevel[cellI];
+        }
+        targetLevel = max(targetLevel, blockLev);
+    }
+    
+    // mark cells to be refined/unrefined based on above criteria:
+    // simply check, if targetLevel lower or higher cellLevel
+    forAll (intRefFld.internalField(), cellI)
+    {
+        intRefFld.primitiveFieldRef()[cellI] = targetLevel[cellI] - cellLevel[cellI];
+        targetFld.primitiveFieldRef()[cellI] = targetLevel[cellI];
+        currFld.primitiveFieldRef()[cellI] = cellLevel[cellI];
+    }
+
+    intRefFld.correctBoundaryConditions();
+    
+    Info<<"Min,max refinement field = " << Foam::min(intRefFld).value() << ", "
+        << Foam::max(intRefFld).value() << endl;
+}
+
+void Foam::dynamicRefineBalancedFvMesh::readRefinementDict()
+{
+    dictionary dynamicMeshDict
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "dynamicMeshDict",
+                time().constant(),
+                *this,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        )
+    );
+
+    const dictionary refineDict
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "dynamicMeshDict",
+                time().constant(),
+                *this,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).subDict("dynamicRefineFvMeshCoeffs")
+    );
+
+    {
+        wordList surfFlds;
+        if (refineDict.readIfPresent("mapSurfaceFields", surfFlds))
+        {
+            // Rework into hashtable.
+            mapSurfaceFields_.resize(surfFlds.size());
+            forAll(surfFlds, i)
+            {
+                mapSurfaceFields_.insert(surfFlds[i], surfFlds[i]);
+            }
+        }
+    }
+
+    nBufferLayers_ = readLabel(refineDict.lookup("nBufferLayers"));
+
+    if( dynamicMeshDict.isDict("refinementControls") )
+    {
+        dictionary refineControlDict = 
+            dynamicMeshDict.subDict("refinementControls");
+        
+        enableRefinementControl_ = 
+            Switch(refineControlDict.lookup("enableRefinementControl"));
+        
+        if( enableRefinementControl_ )
+        {
+            // Overwrite field name entry in dynamicRefineFvMeshCoeffs?
+            // For now you just have to be smart and enter 
+            // 'internalRefinementField' for the field name manually
+            
+            // Read HashTable of field-refinement scalars
+            if( refineControlDict.found("fields") )
+            {
+                fields_ = HashTable< List<scalar> >
+                (
+                    refineControlDict.lookup("fields")
+                );
+            }
+            
+            // Read HashTable of gradient-refinement scalars
+            if( refineControlDict.found("gradients") )
+            {
+                gradFields_ = HashTable< List<scalar> >
+                (
+                    refineControlDict.lookup("gradients")
+                );
+            }
+            
+            // Read HashTable of curl-refinement vectors
+            if( refineControlDict.found("curls") )
+            {
+                curlFields_ = HashTable< List<scalar> >
+                (
+                    refineControlDict.lookup("curls")
+                );
+            }
+            
+            // Read interface refinement data
+            if( refineControlDict.found("interface") )
+            {
+                interface_ = HashTable< dictionary >
+                (
+                    refineControlDict.lookup("interface")
+                );
+            }
+            
+            // Read refinement regions
+            if( refineControlDict.found("regions") )
+            {
+                refinedRegions_ = PtrList<entry>
+                (
+                    refineControlDict.lookup("regions")
+                );
+            }
+        }
+    }
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+            
+Foam::dynamicRefineBalancedFvMesh::dynamicRefineBalancedFvMesh
+(
+    const IOobject& io
+)
+:
+    dynamicRefineFvMesh(io),
+    internalRefinementFieldPtr_(NULL),
+    targetLevelPtr_(NULL),
+    isLevelPtr_(NULL),
+    fields_(),
+    gradFields_(),
+    curlFields_(),
+    interface_(),
+    refinedRegions_(),
+    enableRefinementControl_(false),
+    nBufferLayers_(0),
+    rebalance_(false)
+{
+    readRefinementDict();
+    
+    if( enableRefinementControl_ )
+    {
+        internalRefinementFieldPtr_ = new volScalarField
+        (
+            IOobject
+            (
+                "internalRefinementField",
+                time().timeName(),
+                *this,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            *this,
+            dimensionedScalar("zero", dimless, 0.0)
+        );
+
+        targetLevelPtr_ = new volScalarField
+        (
+            IOobject
+            (
+                "targetLevel",
+                time().timeName(),
+                *this,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            *this,
+            dimensionedScalar("zero", dimless, 0.0)
+        );
+
+        isLevelPtr_ = new volScalarField
+        (
+            IOobject
+            (
+                "isLevel",
+                time().timeName(),
+                *this,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            *this,
+            dimensionedScalar("zero", dimless, 0.0)
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::dynamicRefineBalancedFvMesh::~dynamicRefineBalancedFvMesh()
+{
+
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::dynamicRefineBalancedFvMesh::mapFields(const mapPolyMesh& mpm)
+{
+DebugVar(mpm.nOldCells());
+    dynamicRefineFvMesh::mapFields(mpm);
+
+    // Correct surface fields on introduced internal faces. These get
+    // created out-of-nothing so get an interpolated value.
+    mapNewInternalFaces<scalar>(mpm.faceMap());
+    mapNewInternalFaces<vector>(mpm.faceMap());
+    mapNewInternalFaces<sphericalTensor>(mpm.faceMap());
+    mapNewInternalFaces<symmTensor>(mpm.faceMap());
+    mapNewInternalFaces<tensor>(mpm.faceMap());
+}
+
+
+bool Foam::dynamicRefineBalancedFvMesh::update()
+{
+    //Part 0 - Update internally calculated refinement field
+    readRefinementDict();
+    
+    List<scalar> refinePoints = readRefinementPoints();
+    label refineInterval = refinePoints[2];
+    
+    if( enableRefinementControl_ && time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0 )
+    {
+        updateRefinementField();
+    }
+
+    //Part 1 - Call normal update from dynamicRefineFvMesh
+    bool hasChanged = dynamicRefineFvMesh::update();
+
+    if( Pstream::parRun() && hasChanged)
+    {
+        //Correct values on all coupled patches
+        correctBoundaries<scalar>();
+        correctBoundaries<vector>();
+        correctBoundaries<sphericalTensor>();
+        correctBoundaries<symmTensor>();
+        correctBoundaries<tensor>();      
+    }
+
+    dictionary decomposeParDict
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "decomposeParDict",
+                time().system(),
+                *this,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        )
+    );
+    
+    rebalance_ = false; 
+
+    // Part 2 - Load Balancing
+    {    
+        dictionary refineDict
+        (
+            IOdictionary
+            (
+                IOobject
+                (
+                    "dynamicMeshDict",
+                    time().constant(),
+                    *this,
+                    IOobject::MUST_READ_IF_MODIFIED,
+                    IOobject::NO_WRITE,
+                    false
+                )
+            ).subDict("dynamicRefineFvMeshCoeffs")
+        );
+    
+        Switch enableBalancing = refineDict.lookup("enableBalancing");
+    
+        if ( Pstream::parRun() && hasChanged )
+        {
+            const scalar allowableImbalance = 
+                readScalar(refineDict.lookup("allowableImbalance"));
+            
+            //First determine current level of imbalance - do this for all
+            // parallel runs with a changing mesh, even if balancing is disabled
+            label nGlobalCells = globalData().nTotalCells();
+            scalar idealNCells = scalar(nGlobalCells)/scalar(Pstream::nProcs());
+            scalar localImbalance = mag(scalar(nCells()) - idealNCells);
+            Foam::reduce(localImbalance, maxOp<scalar>());
+            scalar maxImbalance = localImbalance/idealNCells;
+        
+            Info<<"Maximum imbalance = " << 100*maxImbalance << " %" << endl;
+        
+            //If imbalanced, construct weighted coarse graph (level 0) with node
+            // weights equal to their number of subcells. This partitioning works
+            // as long as the number of level 0 cells is several times greater than
+            // the number of processors.
+            if( maxImbalance > allowableImbalance && enableBalancing)
+            {
+                Info << "\n**Solver hold for redistribution at time = "  << time().timeName() << " s" << endl;
+                
+                rebalance_ = true;   
+
+                const labelIOList& cellLevel = meshCutter().cellLevel();
+                Map<label> coarseIDmap(100);
+                labelList uniqueIndex(nCells(),0);
+            
+                label nCoarse = 0;
+
+                forAll(cells(), cellI)
+                {
+                    if( cellLevel[cellI] > 0 )
+                    {
+                        uniqueIndex[cellI] = nCells() + topParentID
+                        (
+                            meshCutter().history().parentIndex(cellI)
+                        );
+                    }
+                    else
+                    {
+                        uniqueIndex[cellI] = cellI;
+                    }
+                
+                    if( coarseIDmap.insert(uniqueIndex[cellI], nCoarse) )
+                    {
+                        ++nCoarse;
+                    }
+                }
+            
+                // Convert to local sequential indexing and calculate coarse
+                // points and weights
+                labelList localIndex(nCells(),0);
+                pointField coarsePoints(nCoarse,vector::zero);
+                scalarField coarseWeights(nCoarse,0.0);
+            
+                forAll(uniqueIndex, cellI)
+                {
+                    localIndex[cellI] = coarseIDmap[uniqueIndex[cellI]];
+                
+                    // If 2D refinement (quadtree) is ever implemented, this '3'
+                    // should be set in general as the number of refinement
+                    // dimensions.
+                    label w = (1 << (3*cellLevel[cellI]));
+                
+                    coarseWeights[localIndex[cellI]] += 1.0;
+                    coarsePoints[localIndex[cellI]] += C()[cellI]/w;
+                }
+            
+                // Set up decomposer - a separate dictionary is used here so
+                // you can use a simple partitioning for decomposePar and
+                // ptscotch for the rebalancing (or any chosen algorithms)
+                autoPtr<decompositionMethod> decomposer
+                (
+                    decompositionMethod::New
+                    (
+                        IOdictionary
+                        (
+                            IOobject
+                            (
+                                "balanceParDict",
+                                time().system(),
+                                *this,
+                                IOobject::MUST_READ_IF_MODIFIED,
+                                IOobject::NO_WRITE
+                            )
+                        )
+                    )
+                );
+            
+                labelList finalDecomp = decomposer().decompose
+                (
+                    *this, 
+                    localIndex,
+                    coarsePoints,
+                    coarseWeights
+                );
+
+                scalar tolDim = globalMeshData::matchTol_ * bounds().mag();
+            
+                Info<< "Distributing the mesh ..." << endl;
+                fvMeshDistribute distributor(*this, tolDim);
+            
+                Info<< "Mapping the fields ..." << endl;
+                autoPtr<mapDistributePolyMesh> map =
+                    distributor.distribute(finalDecomp);
+                  
+                Info<< "Distribute the map ..." << endl;
+                meshCutter_.distribute(map);
+
+
+                Info << "Successfully distributed mesh" << endl;
+
+                scalarList procLoadNew (Pstream::nProcs(), 0.0);
+                procLoadNew[Pstream::myProcNo()] = this->nCells();
+
+                reduce(procLoadNew, sumOp<List<scalar> >());
+
+                scalar overallLoadNew = sum(procLoadNew);
+                scalar averageLoadNew = overallLoadNew/double(Pstream::nProcs());
+
+                Info << "New distribution: " << procLoadNew << endl;
+                Info << "Max deviation: " << max(Foam::mag(procLoadNew-averageLoadNew)/averageLoadNew)*100.0 << " %" << endl;
+            }
+        }
+    }
+
+    if( Pstream::parRun() && rebalance_)
+    {
+        //Correct values on all coupled patches
+        correctBoundaries<scalar>();
+        correctBoundaries<vector>();
+        correctBoundaries<sphericalTensor>();
+        correctBoundaries<symmTensor>();
+        correctBoundaries<tensor>();      
+    }
+
+    return hasChanged;
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.H b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.H
new file mode 100644
index 0000000000..f4944d5c85
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.H
@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 Tyler Voskuilen
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is a derivative work of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+
+Class
+    Foam::dynamicRefineBalancedFvMesh
+
+SourceFiles
+    dynamicRefineBalancedFvMesh.C
+    dynamicRefineBalancedFvMeshTemplates.C
+
+Authors
+    T.G. Voskuilen ( https://github.com/tgvoskuilen/meshBalancing )
+    Daniel Deising <deising@mma.tu-darmstadt.de>
+    Daniel Rettenmaier <rettenmaier@gsc.tu-darmstadt.de>
+    All rights reserved.
+
+Description
+    A fvMesh with built-in multi-criterion refinement
+    and run-time load balancing.
+
+    You may refer to this software as :
+    //- full bibliographic data to be provided
+
+    This code has been developed by :
+        Daniel Rettenmaier (main developer).
+    
+    Method Development and Intellectual Property :
+        T.G. Voskuilen (Purdue University)
+        Timothée Pourpoint <timothee@purdue.edu> (Purdue University
+        Daniel Rettenmaier <rettenmaier@gsc.tu-darmstadt.de>
+        Daniel Deising <deising@mma.tu-darmstadt.de>
+        Holger Marschall <marschall@csi.tu-darmstadt.de>
+        Dieter Bothe <bothe@csi.tu-darmstadt.de>
+        Cameron Tropea <ctropea@sla.tu-darmstadt.de>
+
+        School of Aeronautics and Astronautics
+        Purdue University
+        Mathematical Modeling and Analysis
+        Institute for Fluid Mechanics and Aerodynamics
+        Center of Smart Interfaces
+        Technische Universitaet Darmstadt
+       
+    If you use this software for your scientific work or your publications,
+    please don't forget to acknowledge explicitly the use of it.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dynamicRefineBalancedFvMesh_H
+#define dynamicRefineBalancedFvMesh_H
+
+
+#include "dynamicFvMesh.H"
+#include "dynamicRefineFvMesh.H"
+#include "hexRef8.H"
+#include "PackedBoolList.H"
+#include "Switch.H"
+#include "decompositionMethod.H"
+#include "fvMeshDistribute.H"
+#include "mapDistributePolyMesh.H"
+#include "HashTable.H"
+#include "topoSetSource.H"
+#include "cellSet.H"
+#include "PtrDictionary.H"
+#include "dictionaryEntry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class dynamicRefineBalancedFvMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class dynamicRefineBalancedFvMesh
+:
+    public dynamicRefineFvMesh
+{
+   
+private:
+
+        //- (non-flux)SurfaceFields to map
+        HashTable<word> mapSurfaceFields_;
+
+        //-
+        volScalarField* internalRefinementFieldPtr_;
+        volScalarField* targetLevelPtr_;
+        volScalarField* isLevelPtr_;
+        
+        //-
+        HashTable< List<scalar> > fields_;
+        
+        //-
+        HashTable< List<scalar> > gradFields_;
+        
+        //-
+        HashTable< List<scalar> > curlFields_;
+        
+        //-
+        HashTable< dictionary > interface_;
+        
+        //-
+        PtrList<entry> refinedRegions_; 
+        
+        //-
+        Switch enableRefinementControl_;
+        
+        //-
+        void updateRefinementField();
+        
+        //-
+        void readRefinementDict();
+        
+        //-
+        List<scalar> readRefinementPoints();
+        
+        //- Disallow default bitwise copy construct
+        dynamicRefineBalancedFvMesh(const dynamicRefineBalancedFvMesh&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const dynamicRefineBalancedFvMesh&);
+
+        //-
+        label topParentID(label p);
+
+        //-
+        label nBufferLayers_;
+
+        //-
+        bool rebalance_;
+
+        //- Map non-flux surface<Type>Fields for new internal faces
+        //  (from cell splitting)
+        template <class T>
+        void mapNewInternalFaces(const labelList& faceMap);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("dynamicRefineBalancedFvMesh");
+
+
+    // Constructors
+
+        //- Construct from IOobject
+        explicit dynamicRefineBalancedFvMesh(const IOobject& io);
+
+
+    //- Destructor
+    virtual ~dynamicRefineBalancedFvMesh();
+
+
+    // Member Functions
+
+        //- Update the mesh for both mesh motion and topology change
+        virtual bool update();
+
+        //- Map all fields in time using given map
+        virtual void mapFields(const mapPolyMesh& mpm);
+
+        //- Template to update all volField boundaries
+        template<class Type> void correctBoundaries();
+
+        //-
+        bool rebalance()
+        {
+            return rebalance_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "dynamicRefineBalancedFvMeshTemplates.C"
+#endif
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMeshTemplates.C b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMeshTemplates.C
new file mode 100644
index 0000000000..200e3f831d
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMeshTemplates.C
@@ -0,0 +1,217 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 Tyler Voskuilen
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "volMesh.H"
+#include "fvPatchField.H"
+#include "surfaceFields.H"
+
+template <class T>
+void Foam::dynamicRefineBalancedFvMesh::mapNewInternalFaces
+(
+    const labelList& faceMap
+)
+{
+    typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
+    HashTable<GeoField*> sFlds(this->objectRegistry::lookupClass<GeoField>());
+
+    const labelUList& owner = this->faceOwner();
+    const labelUList& neighbour = this->faceNeighbour();
+    const dimensionedScalar deltaN = 1e-8 / pow(average(this->V()), 1.0 / 3.0);
+
+    forAllIter(typename HashTable<GeoField*>, sFlds, iter)
+    {
+        GeoField& sFld = *iter();
+        if (mapSurfaceFields_.found(iter.key()))
+        {   
+            if (debug)
+            {
+                InfoInFunction << iter.key()<< endl;
+            }
+
+
+            // Create flat version of field
+            Field<T> tsFld(this->nFaces(), pTraits<T>::zero);
+            {
+                forAll(sFld.internalField(), iFace)
+                {
+                    tsFld[iFace] = sFld.internalField()[iFace];
+                }
+                forAll(sFld.boundaryField(), iPatch)
+                {
+                    const fvsPatchField<T>& pfld = sFld.boundaryField()[iPatch];
+
+                    label start = pfld.patch().start();
+                    forAll(pfld, faceI)
+                    {
+                        tsFld[faceI+start] = pfld[faceI];
+                    }
+                }
+            }
+
+
+            // Loop over all faces
+            for (label facei = 0; facei < nInternalFaces(); facei++)
+            {
+                label oldFacei = faceMap[facei];
+
+                //- map surface field on newly generated faces
+                if (oldFacei == -1)
+                {
+                    //- find owner and neighbour cell
+                    const cell& ownFaces = this->cells()[owner[facei]];
+                    const cell& neiFaces = this->cells()[neighbour[facei]];
+
+                    //- loop over all owner/neighbour cell faces
+                    //- and find already mapped ones (master-faces):
+                    T tmpValue = pTraits<T>::zero;
+                    scalar magFld = 0;
+                    label counter = 0;
+
+                    //- simple averaging of all neighbour master-faces
+                    forAll(ownFaces, iFace)
+                    {
+                        if (faceMap[ownFaces[iFace]] != -1)
+                        {
+                            tmpValue += tsFld[ownFaces[iFace]];
+                            magFld += mag(tsFld[ownFaces[iFace]]);
+                            counter++;
+                        }
+                    }
+
+                    forAll(neiFaces, iFace)
+                    {
+                        if (faceMap[neiFaces[iFace]] != -1)
+                        {
+                            tmpValue = tmpValue + tsFld[neiFaces[iFace]];
+                            magFld += mag(tsFld[neiFaces[iFace]]);
+                            counter++;
+                        }
+                    }
+
+                    if(counter > 0)
+                    {
+                        if 
+                        (   
+                            GeometricField<T, fvsPatchField, surfaceMesh>::typeName 
+                                == "surfaceScalarField"
+                        )
+                        {
+                            tmpValue /= counter;
+                        }                            
+                        else if 
+                        (    
+                            GeometricField<T, fvsPatchField, surfaceMesh>::typeName 
+                                == "surfaceVectorField"
+                        )
+                        {
+                            magFld /= counter;
+                            tmpValue *= magFld/(mag(tmpValue)+deltaN.value());
+                        }
+                        else
+                        {
+                            FatalErrorInFunction
+                                << "mapping implementation only valid for"
+                                << " scalar and vector fields! \n Field "
+                                << sFld.name() << " is of type: "
+                                << GeometricField<T, fvsPatchField, surfaceMesh>::typeName 
+                                << abort(FatalError);
+                        }
+                    }
+
+                    sFld[facei] = tmpValue;
+                }
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::dynamicRefineBalancedFvMesh::correctBoundaries()
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
+
+    HashTable<GeoField*> flds(this->objectRegistry::lookupClass<GeoField>());
+
+    forAllIter(typename HashTable<GeoField*>, flds, iter)
+    {
+        GeoField& fld = *iter();
+        
+        //mimic "evaluate" but only for coupled patches (processor or cyclic)
+        // and only for blocking or nonBlocking comms (no scheduled comms)
+        if
+        (
+            Pstream::defaultCommsType == Pstream::commsTypes::blocking
+         || Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking
+        )
+        {
+            label nReq = Pstream::nRequests();
+
+            forAll(fld.boundaryField(), patchi)
+            {
+                if(fld.boundaryField()[patchi].coupled())
+                {
+                    fld.boundaryFieldRef()[patchi].initEvaluate
+                    (
+                        Pstream::defaultCommsType
+                    );
+                }
+            }
+
+            // Block for any outstanding requests
+            if
+            (
+                Pstream::parRun()
+             && Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking
+            )
+            {
+                Pstream::waitRequests(nReq);
+            }
+
+            forAll(fld.boundaryField(), patchi)
+            {
+                if(fld.boundaryField()[patchi].coupled())
+                {
+                    fld.boundaryFieldRef()[patchi].evaluate
+                    (
+                        Pstream::defaultCommsType
+                    );
+                }
+            }
+        }
+        else
+        {
+            //Scheduled patch updates not supported
+            FatalErrorIn
+            (
+                "dynamicRefineBalancedFvMeshTemplates::correctBoundaries"
+            )   << "Unsuported communications type "
+                << Pstream::commsTypeNames[Pstream::defaultCommsType]
+                << exit(FatalError);
+        }
+    }
+}
+
+// ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
index 466b7d4e5e..91aad27d1e 100644
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
@@ -46,7 +46,7 @@ namespace Foam
 
 Foam::label Foam::dynamicRefineFvMesh::count
 (
-    const bitSet& l,
+    const PackedBoolList& l,
     const unsigned int val
 )
 {
@@ -58,9 +58,9 @@ Foam::label Foam::dynamicRefineFvMesh::count
             n++;
         }
 
-        // debug also serves to get-around Clang compiler trying to optimise
+        // debug also serves to get-around Clang compiler trying to optimsie
         // out this forAll loop under O3 optimisation
-        if (debug)
+        if (debug & 128)
         {
             Info<< "n=" << n << endl;
         }
@@ -72,7 +72,7 @@ Foam::label Foam::dynamicRefineFvMesh::count
 
 void Foam::dynamicRefineFvMesh::calculateProtectedCells
 (
-    bitSet& unrefineableCell
+    PackedBoolList& unrefineableCell
 ) const
 {
     if (protectedCell_.empty())
@@ -86,7 +86,7 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
     unrefineableCell = protectedCell_;
 
     // Get neighbouring cell level
-    labelList neiLevel(nBoundaryFaces());
+    labelList neiLevel(nFaces()-nInternalFaces());
 
     for (label facei = nInternalFaces(); facei < nFaces(); facei++)
     {
@@ -103,10 +103,9 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
         forAll(faceNeighbour(), facei)
         {
             label own = faceOwner()[facei];
+            bool ownProtected = unrefineableCell.get(own);
             label nei = faceNeighbour()[facei];
-
-            bool ownProtected = unrefineableCell.test(own);
-            bool neiProtected = unrefineableCell.test(nei);
+            bool neiProtected = unrefineableCell.get(nei);
 
             if (ownProtected && (cellLevel[nei] > cellLevel[own]))
             {
@@ -120,7 +119,7 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
         for (label facei = nInternalFaces(); facei < nFaces(); facei++)
         {
             label own = faceOwner()[facei];
-            bool ownProtected = unrefineableCell.test(own);
+            bool ownProtected = unrefineableCell.get(own);
             if
             (
                 ownProtected
@@ -142,14 +141,16 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
             if (seedFace[facei])
             {
                 label own = faceOwner()[facei];
-                if (unrefineableCell.set(own))
+                if (unrefineableCell.get(own) == 0)
                 {
+                    unrefineableCell.set(own, 1);
                     hasExtended = true;
                 }
 
                 label nei = faceNeighbour()[facei];
-                if (unrefineableCell.set(nei))
+                if (unrefineableCell.get(nei) == 0)
                 {
+                    unrefineableCell.set(nei, 1);
                     hasExtended = true;
                 }
             }
@@ -159,8 +160,9 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
             if (seedFace[facei])
             {
                 label own = faceOwner()[facei];
-                if (unrefineableCell.set(own))
+                if (unrefineableCell.get(own) == 0)
                 {
+                    unrefineableCell.set(own, 1);
                     hasExtended = true;
                 }
             }
@@ -205,63 +207,27 @@ void Foam::dynamicRefineFvMesh::readDict()
 }
 
 
-// Refines cells, maps fields and recalculates (an approximate) flux
-Foam::autoPtr<Foam::mapPolyMesh>
-Foam::dynamicRefineFvMesh::refine
-(
-    const labelList& cellsToRefine
-)
+void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
 {
-    // Mesh changing engine.
-    polyTopoChange meshMod(*this);
-
-    // Play refinement commands into mesh changer.
-    meshCutter_.setRefinement(cellsToRefine, meshMod);
-
-    // Create mesh (with inflation), return map from old to new mesh.
-    autoPtr<mapPolyMesh> mapPtr = meshMod.changeMesh(*this, false);
-    mapPolyMesh& map = *mapPtr;
-
-    Info<< "Refined from "
-        << returnReduce(map.nOldCells(), sumOp<label>())
-        << " to " << globalData().nTotalCells() << " cells." << endl;
-
-    if (debug)
-    {
-        // Check map.
-        for (label facei = 0; facei < nInternalFaces(); facei++)
-        {
-            label oldFacei = map.faceMap()[facei];
-
-            if (oldFacei >= nInternalFaces())
-            {
-                FatalErrorInFunction
-                    << "New internal face:" << facei
-                    << " fc:" << faceCentres()[facei]
-                    << " originates from boundary oldFace:" << oldFacei
-                    << abort(FatalError);
-            }
-        }
-    }
-
-    //    // Remove the stored tet base points
-    //    tetBasePtIsPtr_.clear();
-    //    // Remove the cell tree
-    //    cellTreePtr_.clear();
-
-    // Update fields
-    updateMesh(map);
-
+    dynamicFvMesh::mapFields(mpm);
     // Correct the flux for modified/added faces. All the faces which only
     // have been renumbered will already have been handled by the mapping.
     {
-        const labelList& faceMap = map.faceMap();
-        const labelList& reverseFaceMap = map.reverseFaceMap();
+        const labelList& faceMap = mpm.faceMap();
+        const labelList& reverseFaceMap = mpm.reverseFaceMap();
 
         // Storage for any master faces. These will be the original faces
         // on the coarse cell that get split into four (or rather the
         // master face gets modified and three faces get added from the master)
-        labelHashSet masterFaces(4*cellsToRefine.size());
+        // Estimate number of faces created
+        labelHashSet masterFaces
+        (
+            max
+            (
+                mag(nFaces()-mpm.nOldFaces())/4,
+                nFaces()/100
+            )
+        );
 
         forAll(faceMap, facei)
         {
@@ -393,8 +359,10 @@ Foam::dynamicRefineFvMesh::refine
             }
 
             // Update master faces
-            for (const label facei : masterFaces)
+            forAllConstIter(labelHashSet, masterFaces, iter)
             {
+                label facei = iter.key();
+
                 if (isInternalFace(facei))
                 {
                     phi[facei] = phiU[facei];
@@ -415,6 +383,86 @@ Foam::dynamicRefineFvMesh::refine
         }
     }
 
+    // Correct the flux for injected faces - these are the faces which have
+    // no correspondence to the old mesh (i.e. added without a masterFace, edge
+    // or point). An example is the internal faces from hexRef8.
+    {
+        const labelList& faceMap = mpm.faceMap();
+
+        mapNewInternalFaces<scalar>(this->Sf(), this->magSf(), faceMap);
+        mapNewInternalFaces<vector>(this->Sf(), this->magSf(), faceMap);
+
+        // No oriented fields of more complex type
+        mapNewInternalFaces<sphericalTensor>(faceMap);
+        mapNewInternalFaces<symmTensor>(faceMap);
+        mapNewInternalFaces<tensor>(faceMap);
+    }
+}
+
+
+// Refines cells, maps fields and recalculates (an approximate) flux
+Foam::autoPtr<Foam::mapPolyMesh>
+Foam::dynamicRefineFvMesh::refine
+(
+    const labelList& cellsToRefine
+)
+{
+    // Mesh changing engine.
+    polyTopoChange meshMod(*this);
+
+    // Play refinement commands into mesh changer.
+    meshCutter_.setRefinement(cellsToRefine, meshMod);
+
+    // Create mesh (with inflation), return map from old to new mesh.
+    //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
+
+    Info<< "Refined from "
+        << returnReduce(map().nOldCells(), sumOp<label>())
+        << " to " << globalData().nTotalCells() << " cells." << endl;
+
+    if (debug)
+    {
+        // Check map.
+        for (label facei = 0; facei < nInternalFaces(); facei++)
+        {
+            label oldFacei = map().faceMap()[facei];
+
+            if (oldFacei >= nInternalFaces())
+            {
+                FatalErrorInFunction
+                    << "New internal face:" << facei
+                    << " fc:" << faceCentres()[facei]
+                    << " originates from boundary oldFace:" << oldFacei
+                    << abort(FatalError);
+            }
+        }
+    }
+
+    //    // Remove the stored tet base points
+    //    tetBasePtIsPtr_.clear();
+    //    // Remove the cell tree
+    //    cellTreePtr_.clear();
+
+    // Update fields
+    updateMesh(map);
+
+
+    // Move mesh
+    /*
+    pointField newPoints;
+    if (map().hasMotionPoints())
+    {
+        newPoints = map().preMotionPoints();
+    }
+    else
+    {
+        newPoints = points();
+    }
+    movePoints(newPoints);
+    */
+
+
 
     // Update numbering of cells/vertices.
     meshCutter_.updateMesh(map);
@@ -422,12 +470,12 @@ Foam::dynamicRefineFvMesh::refine
     // Update numbering of protectedCell_
     if (protectedCell_.size())
     {
-        bitSet newProtectedCell(nCells());
+        PackedBoolList newProtectedCell(nCells());
 
         forAll(newProtectedCell, celli)
         {
-            const label oldCelli = map.cellMap()[celli];
-            newProtectedCell.set(celli, protectedCell_.test(oldCelli));
+            label oldCelli = map().cellMap()[celli];
+            newProtectedCell.set(celli, protectedCell_.get(oldCelli));
         }
         protectedCell_.transfer(newProtectedCell);
     }
@@ -435,7 +483,7 @@ Foam::dynamicRefineFvMesh::refine
     // Debug: Check refinement levels (across faces only)
     meshCutter_.checkRefinementLevels(-1, labelList(0));
 
-    return mapPtr;
+    return map;
 }
 
 
@@ -482,21 +530,36 @@ Foam::dynamicRefineFvMesh::unrefine
 
 
     // Change mesh and generate map.
-    autoPtr<mapPolyMesh> mapPtr = meshMod.changeMesh(*this, false);
-    mapPolyMesh& map = *mapPtr;
+    //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
 
     Info<< "Unrefined from "
-        << returnReduce(map.nOldCells(), sumOp<label>())
+        << returnReduce(map().nOldCells(), sumOp<label>())
         << " to " << globalData().nTotalCells() << " cells."
         << endl;
 
     // Update fields
     updateMesh(map);
 
+
+    // Move mesh
+    /*
+    pointField newPoints;
+    if (map().hasMotionPoints())
+    {
+        newPoints = map().preMotionPoints();
+    }
+    else
+    {
+        newPoints = points();
+    }
+    movePoints(newPoints);
+    */
+
     // Correct the flux for modified faces.
     {
-        const labelList& reversePointMap = map.reversePointMap();
-        const labelList& reverseFaceMap = map.reverseFaceMap();
+        const labelList& reversePointMap = map().reversePointMap();
+        const labelList& reverseFaceMap = map().reverseFaceMap();
 
         HashTable<surfaceScalarField*> fluxes
         (
@@ -585,14 +648,14 @@ Foam::dynamicRefineFvMesh::unrefine
     // Update numbering of protectedCell_
     if (protectedCell_.size())
     {
-        bitSet newProtectedCell(nCells());
+        PackedBoolList newProtectedCell(nCells());
 
         forAll(newProtectedCell, celli)
         {
-            label oldCelli = map.cellMap()[celli];
+            label oldCelli = map().cellMap()[celli];
             if (oldCelli >= 0)
             {
-                newProtectedCell.set(celli, protectedCell_.test(oldCelli));
+                newProtectedCell.set(celli, protectedCell_.get(oldCelli));
             }
         }
         protectedCell_.transfer(newProtectedCell);
@@ -601,7 +664,7 @@ Foam::dynamicRefineFvMesh::unrefine
     // Debug: Check refinement levels (across faces only)
     meshCutter_.checkRefinementLevels(-1, labelList(0));
 
-    return mapPtr;
+    return map;
 }
 
 
@@ -688,7 +751,7 @@ void Foam::dynamicRefineFvMesh::selectRefineCandidates
     const scalar lowerRefineLevel,
     const scalar upperRefineLevel,
     const scalarField& vFld,
-    bitSet& candidateCell
+    PackedBoolList& candidateCell
 ) const
 {
     // Get error per cell. Is -1 (not to be refined) to >0 (to be refined,
@@ -711,7 +774,7 @@ void Foam::dynamicRefineFvMesh::selectRefineCandidates
     {
         if (cellError[celli] > 0)
         {
-            candidateCell.set(celli);
+            candidateCell.set(celli, 1);
         }
     }
 }
@@ -721,7 +784,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
 (
     const label maxCells,
     const label maxRefinement,
-    const bitSet& candidateCell
+    const PackedBoolList& candidateCell
 ) const
 {
     // Every refined cell causes 7 extra cells
@@ -731,7 +794,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
 
     // Mark cells that cannot be refined since they would trigger refinement
     // of protected cells (since 2:1 cascade)
-    bitSet unrefineableCell;
+    PackedBoolList unrefineableCell;
     calculateProtectedCells(unrefineableCell);
 
     // Count current selection
@@ -748,8 +811,11 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
             if
             (
                 cellLevel[celli] < maxRefinement
-             && candidateCell.test(celli)
-             && !unrefineableCell.test(celli)
+             && candidateCell.get(celli)
+             && (
+                    unrefineableCell.empty()
+                 || !unrefineableCell.get(celli)
+                )
             )
             {
                 candidates.append(celli);
@@ -766,8 +832,11 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
                 if
                 (
                     cellLevel[celli] == level
-                 && candidateCell.test(celli)
-                 && !unrefineableCell.test(celli)
+                 && candidateCell.get(celli)
+                 && (
+                        unrefineableCell.empty()
+                     || !unrefineableCell.get(celli)
+                    )
                 )
                 {
                     candidates.append(celli);
@@ -802,7 +871,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
 Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
 (
     const scalar unrefineLevel,
-    const bitSet& markedCell,
+    const PackedBoolList& markedCell,
     const scalarField& pFld
 ) const
 {
@@ -815,7 +884,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
     // If we have any protected cells make sure they also are not being
     // unrefined
 
-    bitSet protectedPoint(nPoints());
+    PackedBoolList protectedPoint(nPoints());
 
     if (protectedCell_.size())
     {
@@ -828,9 +897,9 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
             {
                 label cellI = pCells[pCellI];
 
-                if (protectedCell_.test(cellI))
+                if (protectedCell_[cellI])
                 {
-                    protectedPoint.set(pointI);
+                    protectedPoint[pointI] = true;
                     break;
                 }
             }
@@ -870,7 +939,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
 
             forAll(pCells, pCelli)
             {
-                if (markedCell.test(pCells[pCelli]))
+                if (markedCell.get(pCells[pCelli]))
                 {
                     hasMarked = true;
                     break;
@@ -907,7 +976,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
 
 void Foam::dynamicRefineFvMesh::extendMarkedCells
 (
-    bitSet& markedCell
+    PackedBoolList& markedCell
 ) const
 {
     // Mark faces using any marked cell
@@ -915,7 +984,7 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
 
     forAll(markedCell, celli)
     {
-        if (markedCell.test(celli))
+        if (markedCell.get(celli))
         {
             const cell& cFaces = cells()[celli];
 
@@ -933,15 +1002,15 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
     {
         if (markedFace[facei])
         {
-            markedCell.set(faceOwner()[facei]);
-            markedCell.set(faceNeighbour()[facei]);
+            markedCell.set(faceOwner()[facei], 1);
+            markedCell.set(faceNeighbour()[facei], 1);
         }
     }
     for (label facei = nInternalFaces(); facei < nFaces(); facei++)
     {
         if (markedFace[facei])
         {
-            markedCell.set(faceOwner()[facei]);
+            markedCell.set(faceOwner()[facei], 1);
         }
     }
 }
@@ -949,7 +1018,7 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
 
 void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
 (
-    bitSet& protectedCell,
+    PackedBoolList& protectedCell,
     label& nProtected
 ) const
 {
@@ -971,7 +1040,7 @@ void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
                 // Check if cell has already 8 anchor points -> protect cell
                 if (nAnchorPoints[celli] == 8)
                 {
-                    if (protectedCell.set(celli))
+                    if (protectedCell.set(celli, true))
                     {
                         nProtected++;
                     }
@@ -988,9 +1057,9 @@ void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
 
     forAll(protectedCell, celli)
     {
-        if (!protectedCell.test(celli) && nAnchorPoints[celli] != 8)
+        if (!protectedCell[celli] && nAnchorPoints[celli] != 8)
         {
-            protectedCell.set(celli);
+            protectedCell.set(celli, true);
             nProtected++;
         }
     }
@@ -1005,7 +1074,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
     meshCutter_(*this),
     dumpLevel_(false),
     nRefinementIterations_(0),
-    protectedCell_(nCells(), false)
+    protectedCell_(nCells(), 0)
 {
     // Read static part of dictionary
     readDict();
@@ -1034,7 +1103,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
         {
             label celli = pCells[i];
 
-            if (!protectedCell_.test(celli))
+            if (!protectedCell_.get(celli))
             {
                 if (pointLevel[pointi] <= cellLevel[celli])
                 {
@@ -1042,7 +1111,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
 
                     if (nAnchors[celli] > 8)
                     {
-                        protectedCell_.set(celli);
+                        protectedCell_.set(celli, 1);
                         nProtected++;
                     }
                 }
@@ -1105,9 +1174,9 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
         {
             if (protectedFace[facei])
             {
-                protectedCell_.set(faceOwner()[facei]);
+                protectedCell_.set(faceOwner()[facei], 1);
                 nProtected++;
-                protectedCell_.set(faceNeighbour()[facei]);
+                protectedCell_.set(faceNeighbour()[facei], 1);
                 nProtected++;
             }
         }
@@ -1115,7 +1184,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
         {
             if (protectedFace[facei])
             {
-                protectedCell_.set(faceOwner()[facei]);
+                protectedCell_.set(faceOwner()[facei], 1);
                 nProtected++;
             }
         }
@@ -1127,7 +1196,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
 
             if (cFaces.size() < 6)
             {
-                if (protectedCell_.set(celli))
+                if (protectedCell_.set(celli, 1))
                 {
                     nProtected++;
                 }
@@ -1138,7 +1207,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
                 {
                     if (faces()[cFaces[cFacei]].size() < 4)
                     {
-                        if (protectedCell_.set(celli))
+                        if (protectedCell_.set(celli, 1))
                         {
                             nProtected++;
                         }
@@ -1158,6 +1227,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
     }
     else
     {
+
         cellSet protectedCells(*this, "protectedCells", nProtected);
         forAll(protectedCell_, celli)
         {
@@ -1188,9 +1258,9 @@ Foam::dynamicRefineFvMesh::~dynamicRefineFvMesh()
 
 bool Foam::dynamicRefineFvMesh::update()
 {
-    // Re-read dictionary. Usually small so takes trivial amount of time
-    // compared to actual refinement. Also very useful to be able to modify
-    // on-the-fly.
+    // Re-read dictionary. Choosen since usually -small so trivial amount
+    // of time compared to actual refinement. Also very useful to be able
+    // to modify on-the-fly.
     dictionary refineDict
     (
         IOdictionary
@@ -1264,18 +1334,16 @@ bool Foam::dynamicRefineFvMesh::update()
             refineDict.get<scalar>("lowerRefineLevel");
         const scalar upperRefineLevel =
             refineDict.get<scalar>("upperRefineLevel");
-
         const scalar unrefineLevel = refineDict.lookupOrDefault<scalar>
         (
             "unrefineLevel",
             GREAT
         );
-
         const label nBufferLayers =
             refineDict.get<label>("nBufferLayers");
 
         // Cells marked for refinement or otherwise protected from unrefinement.
-        bitSet refineCell(nCells());
+        PackedBoolList refineCell(nCells());
 
         // Determine candidates for refinement (looking at field only)
         selectRefineCandidates
@@ -1316,7 +1384,7 @@ bool Foam::dynamicRefineFvMesh::update()
                     const labelList& cellMap = map().cellMap();
                     const labelList& reverseCellMap = map().reverseCellMap();
 
-                    bitSet newRefineCell(cellMap.size());
+                    PackedBoolList newRefineCell(cellMap.size());
 
                     forAll(cellMap, celli)
                     {
@@ -1324,15 +1392,15 @@ bool Foam::dynamicRefineFvMesh::update()
 
                         if (oldCelli < 0)
                         {
-                            newRefineCell.set(celli);
+                            newRefineCell.set(celli, 1);
                         }
                         else if (reverseCellMap[oldCelli] != celli)
                         {
-                            newRefineCell.set(celli);
+                            newRefineCell.set(celli, 1);
                         }
                         else
                         {
-                            newRefineCell.set(celli, refineCell.test(oldCelli));
+                            newRefineCell.set(celli, refineCell.get(oldCelli));
                         }
                     }
                     refineCell.transfer(newRefineCell);
@@ -1380,7 +1448,7 @@ bool Foam::dynamicRefineFvMesh::update()
 
         if ((nRefinementIterations_ % 10) == 0)
         {
-            // Compact refinement history occasionally (how often?).
+            // Compact refinement history occassionally (how often?).
             // Unrefinement causes holes in the refinementHistory.
             const_cast<refinementHistory&>(meshCutter().history()).compact();
         }
@@ -1403,8 +1471,7 @@ bool Foam::dynamicRefineFvMesh::writeObject
 (
     IOstream::streamFormat fmt,
     IOstream::versionNumber ver,
-    IOstream::compressionType cmp,
-    const bool valid
+    IOstream::compressionType cmp
 ) const
 {
     // Force refinement data to go to the current time directory.
@@ -1412,8 +1479,8 @@ bool Foam::dynamicRefineFvMesh::writeObject
 
     bool writeOk =
     (
-        dynamicFvMesh::writeObject(fmt, ver, cmp, valid)
-     && meshCutter_.write(valid)
+        dynamicFvMesh::writeObject(fmt, ver, cmp)
+     && meshCutter_.write()
     );
 
     if (dumpLevel_)
@@ -1430,7 +1497,7 @@ bool Foam::dynamicRefineFvMesh::writeObject
                 false
             ),
             *this,
-            dimensionedScalar(dimless, Zero)
+            dimensionedScalar("level", dimless, 0)
         );
 
         const labelList& cellLevel = meshCutter_.cellLevel();
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H
index babbceabea..2ccd324486 100644
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,6 +56,7 @@ Description
             (rho*phi none)   //none)
             (ghf none)  //NaN)   //none)
         );
+
         // Write the refinement level as a volScalarField
         dumpLevel       true;
 
@@ -117,10 +118,10 @@ protected:
 
 
         //- Refine cells. Update mesh and fields.
-        autoPtr<mapPolyMesh> refine(const labelList&);
+        virtual autoPtr<mapPolyMesh> refine(const labelList&);
 
         //- Unrefine cells. Gets passed in centre points of cells to combine.
-        autoPtr<mapPolyMesh> unrefine(const labelList&);
+        virtual autoPtr<mapPolyMesh> unrefine(const labelList&);
 
 
         // Selection of cells to un/refine
@@ -185,6 +186,32 @@ protected:
                 label& nProtected
             ) const;
 
+            //- Map single non-flux surface<Type>Field
+            //  for new internal faces (e.g. AMR refine). This currently
+            //  interpolates values from surrounding faces (faces on
+            //  neighbouring cells) that do have a value.
+            template <class T>
+            void mapNewInternalFaces
+            (
+                const labelList& faceMap,
+                GeometricField<T, fvsPatchField, surfaceMesh>&
+            );
+
+            //- Correct surface fields for new faces
+            template <class T>
+            void mapNewInternalFaces(const labelList& faceMap);
+
+            //- Correct surface fields for new faces. Converts any oriented
+            //  fields into non-oriented (i.e. phi -> Uf) before mapping
+            template <class T>
+            void mapNewInternalFaces
+            (
+                const surfaceVectorField& Sf,
+                const surfaceScalarField& magSf,
+                const labelList& faceMap
+            );
+
+
 private:
 
         //- No copy construct
@@ -232,6 +259,9 @@ public:
         //- Update the mesh for both mesh motion and topology change
         virtual bool update();
 
+        //- Map all fields in time using given map.
+        virtual void mapFields(const mapPolyMesh& mpm);
+
 
     // Writing
 
@@ -243,16 +273,20 @@ public:
             IOstream::compressionType cmp,
             const bool valid
         ) const;
-
 };
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+    #include "dynamicRefineFvMeshTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMeshTemplates.C b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMeshTemplates.C
new file mode 100644
index 0000000000..6bf655a7ec
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMeshTemplates.C
@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class T>
+void Foam::dynamicRefineFvMesh::mapNewInternalFaces
+(
+    const labelList& faceMap,
+    GeometricField<T, fvsPatchField, surfaceMesh>& sFld
+)
+{
+    typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
+
+    //- Make flat field for ease of looping
+    Field<T> tsFld(this->nFaces(), pTraits<T>::zero);
+    SubField<T>(tsFld, this->nInternalFaces()) = sFld.internalField();
+
+    const typename GeoField::Boundary& bFld = sFld.boundaryField();
+    forAll(bFld, patchi)
+    {
+        label facei = this->boundaryMesh()[patchi].start();
+        for (const T& val : bFld[patchi])
+        {
+            tsFld[facei++] = val;
+        }
+    }
+
+    const labelUList& owner = this->faceOwner();
+    const labelUList& neighbour = this->faceNeighbour();
+    const cellList& cells = this->cells();
+
+    for (label facei = 0; facei < nInternalFaces(); facei++)
+    {
+        label oldFacei = faceMap[facei];
+
+        // Map surface field on newly generated faces by obtaining the
+        // hull of the outside faces
+        if (oldFacei == -1)
+        {
+            // Loop over all owner/neighbour cell faces
+            // and find already mapped ones (master-faces):
+            T tmpValue = pTraits<T>::zero;
+            label counter = 0;
+
+            const cell& ownFaces = cells[owner[facei]];
+            for (auto ownFacei : ownFaces)
+            {
+                if (faceMap[ownFacei] != -1)
+                {
+                    tmpValue += tsFld[ownFacei];
+                    counter++;
+                }
+            }
+
+            const cell& neiFaces = cells[neighbour[facei]];
+            for (auto neiFacei : neiFaces)
+            {
+                if (faceMap[neiFacei] != -1)
+                {
+                    tmpValue += tsFld[neiFacei];
+                    counter++;
+                }
+            }
+
+            if (counter > 0)
+            {
+                sFld[facei] = tmpValue/counter;
+            }
+        }
+    }
+}
+
+
+template<class T>
+void Foam::dynamicRefineFvMesh::mapNewInternalFaces
+(
+    const labelList& faceMap
+)
+{
+    typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
+    HashTable<GeoField*> sFlds(this->objectRegistry::lookupClass<GeoField>());
+
+    forAllIter(typename HashTable<GeoField*>, sFlds, iter)
+    {
+        //if (mapSurfaceFields_.found(iter.key()))
+        {
+            if (debug)
+            {
+                Info<< "dynamicRefineFvMesh::mapNewInternalFaces():"
+                    << " Mapping new internal faces by interpolation on "
+                    << iter.key()<< endl;
+            }
+
+            GeoField& sFld = *iter();
+
+            if (sFld.oriented()())
+            {
+                WarningInFunction << "Ignoring mapping oriented field "
+                    << sFld.name() << " since of type " << sFld.type()
+                    << endl;
+            }
+            else
+            {
+                mapNewInternalFaces(faceMap, sFld);
+            }
+        }
+    }
+}
+
+template<class T>
+void Foam::dynamicRefineFvMesh::mapNewInternalFaces
+(
+    const surfaceVectorField& Sf,
+    const surfaceScalarField& magSf,
+    const labelList& faceMap
+)
+{
+    typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
+    HashTable<GeoField*> sFlds(this->objectRegistry::lookupClass<GeoField>());
+
+    forAllIter(typename HashTable<GeoField*>, sFlds, iter)
+    {
+        //if (mapSurfaceFields_.found(iter.key()))
+        {
+            if (debug)
+            {
+                Info<< "dynamicRefineFvMesh::mapNewInternalFaces():"
+                    << " Mapping new internal faces by interpolation on "
+                    << iter.key()<< endl;
+            }
+
+            GeoField& sFld = *iter();
+
+            if (sFld.oriented()())
+            {
+                if (debug)
+                {
+                    Info<< "dynamicRefineFvMesh::mapNewInternalFaces(): "
+                        << "Converting oriented field " << iter.key()
+                        << " to intensive field and mapping" << endl;
+                }
+
+                // Assume any oriented field is face area weighted (i.e. a flux)
+                // Convert to intensive (& oriented) before mapping. Untested.
+
+                typedef GeometricField
+                <
+                    typename outerProduct<vector, T>::type,
+                    fvsPatchField,
+                    surfaceMesh
+                > NormalGeoField;
+
+                // Convert to intensive and non oriented
+                NormalGeoField fFld(sFld*Sf/Foam::sqr(magSf));
+
+                // Interpolate
+                mapNewInternalFaces(faceMap, fFld);
+
+                // Convert back to extensive and oriented
+                sFld = (fFld & Sf);
+            }
+            else
+            {
+                mapNewInternalFaces(faceMap, sFld);
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
index a66bd7c486..f535a71cd4 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
@@ -210,7 +210,7 @@ Foam::label Foam::hexRef8::addInternalFace
                 nei,                        // neighbour
                 -1,                         // master point
                 -1,                         // master edge
-                meshFacei,                  // master face for addition
+                -1,                         // master face for addition
                 false,                      // flux flip
                 -1,                         // patch for face
                 -1,                         // zone for face
diff --git a/src/finiteVolume/fvMesh/fvMesh.C b/src/finiteVolume/fvMesh/fvMesh.C
index a0a032935b..41dc37f784 100644
--- a/src/finiteVolume/fvMesh/fvMesh.C
+++ b/src/finiteVolume/fvMesh/fvMesh.C
@@ -571,6 +571,13 @@ const Foam::lduAddressing& Foam::fvMesh::lduAddr() const
 {
     if (!lduPtr_)
     {
+        if (debug)
+        {
+            InfoInFunction
+                << " calculating fvMeshLduAddressing from nFaces:"
+                << nFaces() << endl;
+        }
+
         lduPtr_ = new fvMeshLduAddressing(*this);
     }
 
@@ -821,6 +828,9 @@ void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm)
     // Update polyMesh. This needs to keep volume existent!
     polyMesh::updateMesh(mpm);
 
+    // Our slice of the addressing is no longer valid
+    deleteDemandDrivenData(lduPtr_);
+
     if (VPtr_)
     {
         // Grab old time volumes if the time has been incremented
diff --git a/tutorials/multiphase/interFoam/laminar/damBreakWithObstacle/constant/dynamicMeshDict b/tutorials/multiphase/interFoam/laminar/damBreakWithObstacle/constant/dynamicMeshDict
index 88613b395f..88a8b82bb8 100644
--- a/tutorials/multiphase/interFoam/laminar/damBreakWithObstacle/constant/dynamicMeshDict
+++ b/tutorials/multiphase/interFoam/laminar/damBreakWithObstacle/constant/dynamicMeshDict
@@ -49,6 +49,7 @@ correctFluxes
     (rhoPhi none)
     (alphaPhi0.water none)
     (ghf none)
+    (alphaPhiUn none)
 );
 
 // Write the refinement level as a volScalarField