diff --git a/README.txt b/README.txt
new file mode 100644
index 0000000000000000000000000000000000000000..273b4dad2085f94fcdb458d476d3ca83ac62eee9
--- /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 b8432c07c905108747f9b94eb4cffefd26d1e0bd..bb04e8ac64705ac6dc159ca970a4bc5dc7320ca9 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 3540c13bfc8d77a85b5642874b6806561cfdacd8..718fde7d5c5de0e42db7e1a715730dd1bcc68f47 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 0000000000000000000000000000000000000000..08d1ef65388aa50771de1f3b78493a2d4bbdad0b
--- /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 0000000000000000000000000000000000000000..a8df1cb7efa27c0931063127ac6dd3d65dd2cefc
--- /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 0000000000000000000000000000000000000000..d5098ffe6805cbe0c7b58dd3c47c91423036ca70
--- /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 0000000000000000000000000000000000000000..f4944d5c850a8c0df73c88e02eeb1654bf4d7aa9
--- /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 0000000000000000000000000000000000000000..200e3f831d3356fb5316e8785ae55715f6c4add3
--- /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);
+        }
+    }
+}
+
+// ************************************************************************* //