From 4a1650335905e2f821aa8b50fa15829a06dab3e3 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 16 May 2018 15:10:47 +0100
Subject: [PATCH] Enabled dynamic refinement in immersed boundary dynamic mesh

---
 .../immersedBoundaryDynamicMesh/Make/files    |   3 +
 .../immersedBoundaryDynamicMesh/Make/options  |   2 +
 ...undaryDynamicRefineSolidBodyMotionFvMesh.C | 631 +-----------------
 ...undaryDynamicRefineSolidBodyMotionFvMesh.H |  54 +-
 .../immersedBoundaryRefinement.C              | 214 ++++++
 .../immersedBoundaryRefinement.H              | 114 ++++
 6 files changed, 341 insertions(+), 677 deletions(-)
 create mode 100644 src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryRefinement/immersedBoundaryRefinement.C
 create mode 100644 src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryRefinement/immersedBoundaryRefinement.H

diff --git a/src/immersedBoundary/immersedBoundaryDynamicMesh/Make/files b/src/immersedBoundary/immersedBoundaryDynamicMesh/Make/files
index a090e5ac8..963e02966 100644
--- a/src/immersedBoundary/immersedBoundaryDynamicMesh/Make/files
+++ b/src/immersedBoundary/immersedBoundaryDynamicMesh/Make/files
@@ -2,6 +2,9 @@ movingImmersedBoundary/movingImmersedBoundary.C
 refineImmersedBoundaryMesh/refineImmersedBoundaryMesh.C
 
 immersedBoundarySolidBodyMotionFvMesh/immersedBoundarySolidBodyMotionFvMesh.C
+
+immersedBoundaryRefinement/immersedBoundaryRefinement.C
+
 immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.C
 
 LIB = $(FOAM_LIBBIN)/libimmersedBoundaryDynamicMesh
diff --git a/src/immersedBoundary/immersedBoundaryDynamicMesh/Make/options b/src/immersedBoundary/immersedBoundaryDynamicMesh/Make/options
index 3fd50b686..10959e1e8 100644
--- a/src/immersedBoundary/immersedBoundaryDynamicMesh/Make/options
+++ b/src/immersedBoundary/immersedBoundaryDynamicMesh/Make/options
@@ -6,6 +6,7 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/topoChangerFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/meshMotion/solidBodyMotion/lnInclude
 
 LIB_LIBS = \
@@ -16,4 +17,5 @@ LIB_LIBS = \
     -lsampling \
     -ldynamicMesh \
     -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
     -lsolidBodyMotion
diff --git a/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.C b/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.C
index b0833651e..1bcb5c10e 100644
--- a/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.C
+++ b/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.C
@@ -55,444 +55,13 @@ addToRunTimeSelectionTable
 } // End namespace Foam
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::label Foam::immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::count
-(
-    const PackedBoolList& l,
-    const unsigned int val
-)
-{
-    label n = 0;
-    forAll (l, i)
-    {
-        if (l.get(i) == val)
-        {
-            n++;
-        }
-
-        // debug also serves to get-around Clang compiler trying to optimise
-        // out this forAll loop under O3 optimisation
-
-        if (debug)
-        {
-            Info<< "n=" << n << endl;
-        }
-    }
-
-    return n;
-}
-
-
-Foam::autoPtr<Foam::mapPolyMesh>
-Foam::immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::refine
-(
-    const labelList& cellsToRefine
-)
-{
-    // Mesh changing engine.
-    directTopoChange meshMod(*this);
-
-    // Play refinement commands into mesh changer.
-    meshCutter_.setRefinement(cellsToRefine, meshMod);
-
-    // Create mesh (with inflation), return map from old to new mesh.
-    //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
-    autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
-
-    Info<< "Refined from "
-        << returnReduce(map().nOldCells(), sumOp<label>())
-        << " to " << globalData().nTotalCells() << " cells." << endl;
-
-    if (debug)
-    {
-        // Check map.
-        for (label faceI = 0; faceI < nInternalFaces(); faceI++)
-        {
-            label oldFaceI = map().faceMap()[faceI];
-
-            if (oldFaceI >= nInternalFaces())
-            {
-                FatalErrorIn
-                (
-                    "immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::"
-                    "refine(const labelList&)"
-                )   << "New internal face:" << faceI
-                    << " fc:" << faceCentres()[faceI]
-                    << " originates from boundary oldFace:" << oldFaceI
-                    << abort(FatalError);
-            }
-        }
-    }
-
-
-    // Update fields
-    updateMesh(map);
-
-    // Update numbering of cells/vertices.
-    meshCutter_.updateMesh(map);
-
-    // Debug: Check refinement levels (across faces only)
-    meshCutter_.checkRefinementLevels(-1, labelList(0));
-
-    return map;
-}
-
-
-Foam::autoPtr<Foam::mapPolyMesh>
-Foam::immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::unrefine
-(
-    const labelList& splitPoints
-)
-{
-    directTopoChange meshMod(*this);
-
-    // Play refinement commands into mesh changer.
-    meshCutter_.setUnrefinement(splitPoints, meshMod);
-
-    // Save information on faces that will be combined
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Find the faceMidPoints on cells to be combined.
-    // for each face resulting of split of face into four store the
-    // midpoint
-    Map<label> faceToSplitPoint(3*splitPoints.size());
-
-    {
-        forAll (splitPoints, i)
-        {
-            label pointI = splitPoints[i];
-
-            const labelList& pEdges = pointEdges()[pointI];
-
-            forAll (pEdges, j)
-            {
-                label otherPointI = edges()[pEdges[j]].otherVertex(pointI);
-
-                const labelList& pFaces = pointFaces()[otherPointI];
-
-                forAll (pFaces, pFaceI)
-                {
-                    faceToSplitPoint.insert(pFaces[pFaceI], otherPointI);
-                }
-            }
-        }
-    }
-
-    // Change mesh and generate map.
-    autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
-
-    Info<< "Unrefined from "
-        << returnReduce(map().nOldCells(), sumOp<label>())
-        << " to " << globalData().nTotalCells() << " cells."
-        << endl;
-
-    // Update fields
-    updateMesh(map);
-
-    // Update numbering of cells/vertices.
-    meshCutter_.updateMesh(map);
-
-    // Debug: Check refinement levels (across faces only)
-    meshCutter_.checkRefinementLevels(-1, labelList(0));
-
-    return map;
-}
-
-
-Foam::labelList
-Foam::immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::
-selectRefineCells
-(
-    const label maxCells,
-    const label maxRefLevel,
-    const scalar refinementDistance,
-    PackedBoolList& candidateCell
-) const
-{
-    // Calculate distance to immersed boundary
-    // Note: minimum distance is inside (WET) and implies a live cell
-    scalarField cellDistance(nCells(), -GREAT);
-
-    const vector span(GREAT, GREAT, GREAT);
-
-    forAll (boundaryMesh(), patchI)
-    {
-        if (isA<immersedBoundaryPolyPatch>(boundaryMesh()[patchI]))
-        {
-            const immersedBoundaryPolyPatch& ibPatch =
-                refCast<const immersedBoundaryPolyPatch>
-                (
-                    boundaryMesh()[patchI]
-                );
-
-            // Create a distance function
-            triSurfaceDistance tsd
-            (
-                ibPatch.triSurfSearch(),
-                span,
-                ibPatch.internalFlow(),
-                false                     // Do not iterate
-            );
-
-            // Calculate distance
-            cellDistance = Foam::max
-            (
-                cellDistance,
-                tsd.distance(cellCentres())
-            );
-        }
-    }
-
-    Info<< "Cell distance: " << min(cellDistance) << " " << max(cellDistance) << endl;
-        
-    // Create a sorting criterion
-    forAll (cellDistance, cellI)
-    {
-        if
-        (
-             cellDistance[cellI] > -refinementDistance
-           && cellDistance[cellI] < 0
-        )
-        {
-            // Found a refinement cell
-            candidateCell.set(cellI, 1);
-        }
-    }
-
-    // Note: Assuming predominantly hex mesh, i.e. every refined cell causes 7
-    // extra cells. This might lead to slight over shoot of maximum number of
-    // cells.
-    label nTotToRefine = (maxCells - globalData().nTotalCells())/7;
-
-    const labelList& cellLevel = meshCutter_.cellLevel();
-
-    // Count current selection
-    label nCandidates = returnReduce(count(candidateCell, 1), sumOp<label>());
-
-    // Collect all cells
-    dynamicLabelList candidates(nCells());
-
-    if (nCandidates < nTotToRefine)
-    {
-        // We won't exceed the maximum number of cells. Collect all candidate
-        // cells that have refinement level smaller than max level
-        forAll (candidateCell, cellI)
-        {
-            if
-            (
-                cellLevel[cellI] < maxRefLevel
-             && candidateCell.get(cellI) == 1
-            )
-            {
-                candidates.append(cellI);
-            }
-        }
-    }
-    else
-    {
-        // We will exceed the maximum number of cells. Simply truncate the
-        // list. It is possible to prefer certain cells based on error field.
-        for (label level = 0; level < maxRefLevel; level++)
-        {
-            forAll (candidateCell, cellI)
-            {
-                if
-                (
-                    cellLevel[cellI] == level
-                 && candidateCell.get(cellI) == 1
-                )
-                {
-                    candidates.append(cellI);
-                }
-            }
-
-            if (returnReduce(candidates.size(), sumOp<label>()) > nTotToRefine)
-            {
-                break;
-            }
-        }
-    }
-
-    // Guarantee 2:1 refinement after refinement
-    labelList consistentSet
-    (
-        meshCutter_.consistentRefinement
-        (
-            candidates.shrink(),
-            true               // Add to set to guarantee 2:1
-        )
-    );
-
-    Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
-        << " cells for refinement out of " << globalData().nTotalCells()
-        << "." << endl;
-
-    return consistentSet;
-}
-
-
-Foam::labelList
-Foam::immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::selectUnrefinePoints
-(
-    const scalar coarsenDistance,
-    const PackedBoolList& markedCell
-) const
-{
-    // All points that can be unrefined
-    const labelList splitPointLabels = meshCutter_.getSplitPoints();
-
-    // Get split points from all points
-    pointField splitPoints(points(), splitPointLabels);
-    
-    // Calculate distance to immersed boundary
-    // Note: minimum distance is inside (WET) and implies a live cell
-    scalarField splitPointDistance(splitPoints.size(), -GREAT);
-
-    const vector span(GREAT, GREAT, GREAT);
-
-    forAll (boundaryMesh(), patchI)
-    {
-        if (isA<immersedBoundaryPolyPatch>(boundaryMesh()[patchI]))
-        {
-            const immersedBoundaryPolyPatch& ibPatch =
-                refCast<const immersedBoundaryPolyPatch>
-                (
-                    boundaryMesh()[patchI]
-                );
-
-            // Create a distance function
-            triSurfaceDistance tsd
-            (
-                ibPatch.triSurfSearch(),
-                span,
-                ibPatch.internalFlow(),
-                false                     // Do not iterate
-            );
-
-            // Calculate distance
-            splitPointDistance = Foam::max
-            (
-                splitPointDistance,
-                tsd.distance(splitPoints)
-            );
-        }
-    }
-
-    Info<< "Split point distance: " << min(splitPointDistance) << " " << max(splitPointDistance) << endl;
-        
-    dynamicLabelList newSplitPoints(splitPointLabels.size());
-
-    const labelListList& pc = pointCells();
-
-    // Only checking split points
-    forAll (splitPointLabels, pointI)
-    {
-        // Coarse live cells next to immersed boundary or dead cells
-        if
-        (
-            splitPointDistance[pointI] < -coarsenDistance
-         || splitPointDistance[pointI] > SMALL
-        )
-        {
-            // Check that all cells are not marked
-            const labelList& pCells = pc[splitPointLabels[pointI]];
-
-            bool hasMarked = false;
-
-            forAll (pCells, pCellI)
-            {
-                if (markedCell.get(pCells[pCellI]) == 1)
-                {
-                    hasMarked = true;
-                    break;
-                }
-            }
-
-            if (!hasMarked)
-            {
-                newSplitPoints.append(splitPointLabels[pointI]);
-            }
-        }
-    }
-
-    newSplitPoints.shrink();
-
-    // Guarantee 2:1 refinement after unrefinement
-    labelList consistentSet
-    (
-        meshCutter_.consistentUnrefinement
-        (
-            newSplitPoints,
-            false
-        )
-    );
-
-    Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
-        << " split points out of a possible "
-        << returnReduce(splitPointLabels.size(), sumOp<label>())
-        << "." << endl;
-
-    return consistentSet;
-}
-
-
-void Foam::immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::extendMarkedCells
-(
-    PackedBoolList& markedCell
-) const
-{
-    // Mark faces using any marked cell
-    boolList markedFace(nFaces(), false);
-
-    forAll (markedCell, cellI)
-    {
-        if (markedCell.get(cellI) == 1)
-        {
-            const cell& cFaces = cells()[cellI];
-
-            forAll (cFaces, i)
-            {
-                markedFace[cFaces[i]] = true;
-            }
-        }
-    }
-
-    syncTools::syncFaceList(*this, markedFace, orEqOp<bool>(), false);
-
-    const labelList& own = faceOwner();
-    const labelList& nei = faceNeighbour();
-    
-    // Update cells using any markedFace
-    for (label faceI = 0; faceI < nInternalFaces(); faceI++)
-    {
-        if (markedFace[faceI])
-        {
-            markedCell.set(own[faceI], 1);
-            markedCell.set(nei[faceI], 1);
-        }
-    }
-
-    for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
-    {
-        if (markedFace[faceI])
-        {
-            markedCell.set(own[faceI], 1);
-        }
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::
 immersedBoundaryDynamicRefineSolidBodyMotionFvMesh(const IOobject& io)
 :
-    dynamicFvMesh(io),
-    ibMotions_(),
-    meshCutter_(*this),
-    curTimeIndex_(-1),
-    nRefinementIterations_(0)
+    dynamicPolyRefinementFvMesh(io),
+    ibMotions_()
 {
     // Read motion function for all regions
     dictionary dynamicMeshCoeffs
@@ -510,7 +79,7 @@ immersedBoundaryDynamicRefineSolidBodyMotionFvMesh(const IOobject& io)
         ).subDict(typeName + "Coeffs")
     );
 
-    PtrList<entry> motionDicts(dynamicMeshCoeffs.lookup("motionFunctions"));
+    PtrList<entry> motionDicts(refinementDict().lookup("motionFunctions"));
 
     ibMotions_.setSize(motionDicts.size());
 
@@ -542,76 +111,16 @@ immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::
 bool immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::update()
 {
     // Handling multiple calls in a single time step
-    if (curTimeIndex_ == this->time().timeIndex())
+    if (!firstUpdate())
     {
         // This is not the first call to update, simply return false
         return false;
     }
-
-    // Re-read dictionary. Choosen since usually -small so trivial amount
-    // of time compared to actual refinement. Also very useful to be able
-    // to modify on-the-fly.
-    dictionary refineDict
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "dynamicMeshDict",
-                time().constant(),
-                *this,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            )
-        ).subDict(typeName + "Coeffs")
-    );
-
-    label refineInterval = readLabel(refineDict.lookup("refineInterval"));
-    const label maxRefLevel = readLabel(refineDict.lookup("maxRefLevel"));
-    const label maxCells = readLabel(refineDict.lookup("maxCells"));
-
-    if (refineInterval < 0)
-    {
-        WarningIn
-        (
-            "immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::update()"
-        )   << "Illegal refineInterval " << refineInterval << nl
-            << "The refineInterval setting in the dynamicMeshDict should"
-            << " be >= 1." << nl
-            << exit(FatalError);
-    }
-
-    if (maxCells <= 0)
-    {
-        WarningIn
-        (
-            "immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::update()"
-        )   << "Illegal maximum number of cells " << maxCells << nl
-            << "The maxCells setting in the dynamicMeshDict should"
-            << " be > 0." << nl
-            << exit(FatalError);
-    }
-
-    if (maxRefLevel <= 0)
-    {
-        WarningIn
-        (
-            "immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::update()"
-        )   << "Illegal maximum refinement level " << maxRefLevel << nl
-            << "The maxCells setting in the dynamicMeshDict should"
-            << " be > 0." << nl
-            << exit(FatalError);
-    }
-
-    // Move points based on new motion
-    if (curTimeIndex_ < this->time().timeIndex())
+    else
     {
         // Grab old volumes before moving the mesh
         // This MUST be followed by mesh motion.  HJ, 29/Dec/2017
         setV0();
-
-        curTimeIndex_ = this->time().timeIndex();
     }
 
     forAll (ibMotions_, ibI)
@@ -619,135 +128,7 @@ bool immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::update()
         ibMotions_[ibI].movePoints();
     }
 
-    bool hasChanged = false;
-
-    if
-    (
-        time().timeIndex() > 0
-     && refineInterval > 0
-     && maxCells > 1
-     && maxRefLevel > 0
-     && time().timeIndex() % refineInterval == 0)
-    {
-        const scalar coarseningDistance =
-            readScalar(refineDict.lookup("coarseningDistance"));
-
-        const scalar refinementDistance =
-            readScalar(refineDict.lookup("refinementDistance"));
-
-        const label nBufferLayers =
-            readLabel(refineDict.lookup("nBufferLayers"));
-
-        // Note:
-        // Positive cell distance implies a dead cell: coarsen if possible
-        // Negative cell distance implies live cell
-        // For negative distance greater than -refinementDistance, refine
-        // For neegative distance smaller than -coarseningDistance, coarsen
-        
-        // Cells marked for refinement or otherwise protected from unrefinement.
-        PackedBoolList refineCell(nCells());
-
-        if (globalData().nTotalCells() < maxCells)
-        {
-            // Determine candidates for refinement
-            labelList cellsToRefine =
-                selectRefineCells
-                (
-                    maxCells,
-                    maxRefLevel,
-                    refinementDistance,
-                    refineCell
-                );
-
-            label nCellsToRefine = returnReduce
-            (
-                cellsToRefine.size(),
-                sumOp<label>()
-            );
-
-            if (nCellsToRefine > 0)
-            {
-                // Refine/update mesh and map fields
-                autoPtr<mapPolyMesh> map = refine(cellsToRefine);
-
-                // Update refineCell. Note that some of the marked ones have
-                // not been refined due to constraints.
-                {
-                    const labelList& cellMap = map().cellMap();
-                    const labelList& reverseCellMap = map().reverseCellMap();
-
-                    PackedBoolList newRefineCell(cellMap.size());
-
-                    forAll (cellMap, cellI)
-                    {
-                        label oldCellI = cellMap[cellI];
-
-                        if (oldCellI < 0)
-                        {
-                            newRefineCell.set(cellI, 1);
-                        }
-                        else if (reverseCellMap[oldCellI] != cellI)
-                        {
-                            newRefineCell.set(cellI, 1);
-                        }
-                        else
-                        {
-                            newRefineCell.set(cellI, refineCell.get(oldCellI));
-                        }
-                    }
-
-                    refineCell.transfer(newRefineCell);
-                }
-
-                // Extend with a buffer layer to prevent neighbouring points
-                // being unrefined.
-                for (label i = 0; i < nBufferLayers; i++)
-                {
-                    extendMarkedCells(refineCell);
-                }
-
-                hasChanged = true;
-            }
-        }
-
-        // Coarsening
-        {
-            // Select unrefineable points that are not marked in refineCell
-            labelList pointsToUnrefine =
-                selectUnrefinePoints
-                (
-                    coarseningDistance,
-                    refineCell
-                );
-
-            label nSplitPoints = returnReduce
-            (
-                pointsToUnrefine.size(),
-                sumOp<label>()
-            );
-
-            if (nSplitPoints > 0)
-            {
-                // Refine/update mesh
-                unrefine(pointsToUnrefine);
-
-                hasChanged = true;
-            }
-        }
-
-
-        if ((nRefinementIterations_ % 10) == 0)
-        {
-            // Compact refinement history occassionally (how often?).
-            // Unrefinement causes holes in the polyRefinementHistory.
-            const_cast<polyRefinementHistory&>
-            (
-                meshCutter().history()
-            ).compact();
-        }
-
-        nRefinementIterations_++;
-    }
+    bool hasChanged = dynamicPolyRefinementFvMesh::update();
 
     changing(hasChanged);
 
diff --git a/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.H b/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.H
index 9fb7ccd1d..3b0d16969 100644
--- a/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.H
+++ b/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.H
@@ -45,11 +45,8 @@ Notes
 #ifndef immersedBoundaryDynamicRefineSolidBodyMotionFvMesh_H
 #define immersedBoundaryDynamicRefineSolidBodyMotionFvMesh_H
 
-#include "dynamicFvMesh.H"
-#include "polyRef.H"
-#include "PackedBoolList.H"
+#include "dynamicPolyRefinementFvMesh.H"
 #include "movingImmersedBoundary.H"
-#include "triSurfaceDistance.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -62,7 +59,7 @@ namespace Foam
 
 class immersedBoundaryDynamicRefineSolidBodyMotionFvMesh
 :
-    public dynamicFvMesh
+    public dynamicPolyRefinementFvMesh
 {
 
     // Private data
@@ -70,15 +67,6 @@ class immersedBoundaryDynamicRefineSolidBodyMotionFvMesh
         //- Immersed boundary motion control function
         PtrList<movingImmersedBoundary> ibMotions_;
 
-        //- Mesh cutting engine
-        polyRef meshCutter_;
-
-        //- Current mesh motion time index
-        label curTimeIndex_;
-
-        //- Number of refinement/unrefinement steps done so far.
-        label nRefinementIterations_;
-
 
     // Private Member Functions
 
@@ -95,38 +83,6 @@ class immersedBoundaryDynamicRefineSolidBodyMotionFvMesh
         );
 
 
-        //- Count set/unset elements in packedlist.
-        static label count(const PackedBoolList&, const unsigned int);
-
-        //- Refine cells. Update mesh and fields.
-        autoPtr<mapPolyMesh> refine(const labelList&);
-
-        //- Unrefine cells. Gets passed in centre points of cells to combine.
-        autoPtr<mapPolyMesh> unrefine(const labelList&);
-
-
-        // Selection of cells to refine and coarsen
-
-            //- Subset candidate cells for refinement
-            labelList selectRefineCells
-            (
-                const label maxCells,
-                const label maxRefLevel,
-                const scalar refinementDistance,
-                PackedBoolList& candidateCell
-            ) const;
-
-            //- Select points that can be unrefined
-            labelList selectUnrefinePoints
-            (
-                const scalar coarsenDistance,
-                const PackedBoolList& markedCell
-            ) const;
-
-            //- Extend markedCell with cell-face-cell
-            void extendMarkedCells(PackedBoolList& markedCell) const;
-
-
 public:
 
     //- Runtime type information
@@ -148,12 +104,6 @@ public:
 
     // Member Functions
 
-        //- Direct access to the refinement engine
-        const polyRef& meshCutter() const
-        {
-            return meshCutter_;
-        }
-
         //- Update the mesh for both mesh motion and topology change
         virtual bool update();
 };
diff --git a/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryRefinement/immersedBoundaryRefinement.C b/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryRefinement/immersedBoundaryRefinement.C
new file mode 100644
index 000000000..6bd81d862
--- /dev/null
+++ b/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryRefinement/immersedBoundaryRefinement.C
@@ -0,0 +1,214 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "immersedBoundaryRefinement.H"
+#include "fvMesh.H"
+#include "immersedBoundaryPolyPatch.H"
+#include "triSurfaceDistance.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(immersedBoundaryRefinement, 0);
+addToRunTimeSelectionTable
+(
+    refinementSelection,
+    immersedBoundaryRefinement,
+    dictionary
+);
+
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::immersedBoundaryRefinement::immersedBoundaryRefinement
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    refinementSelection(mesh, dict),
+    refinementDistance_(readScalar(coeffDict().lookup("refinementDistance"))),
+    coarseningDistance_(readScalar(coeffDict().lookup("coarseningDistance")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor* * * * * * * * * * * * * * * * //
+
+Foam::immersedBoundaryRefinement::~immersedBoundaryRefinement()
+{}
+
+
+// * * * * * * * * * * * * Public Member Functions * * * * * * * * * * * * * //
+
+Foam::Xfer<Foam::labelList>
+Foam::immersedBoundaryRefinement::refinementCellCandidates() const
+{
+    // Calculate distance to immersed boundary
+    // Note: minimum distance is inside (WET) and implies a live cell
+    scalarField cellDistance(mesh().nCells(), -GREAT);
+
+    const polyBoundaryMesh& bMesh = mesh().boundaryMesh();
+    
+    const vector span(GREAT, GREAT, GREAT);
+
+    forAll (bMesh, patchI)
+    {
+        if (isA<immersedBoundaryPolyPatch>(bMesh[patchI]))
+        {
+            const immersedBoundaryPolyPatch& ibPatch =
+                refCast<const immersedBoundaryPolyPatch>
+                (
+                    bMesh[patchI]
+                );
+
+            // Create a distance function
+            triSurfaceDistance tsd
+            (
+                ibPatch.triSurfSearch(),
+                span,
+                ibPatch.internalFlow(),
+                false                     // Do not iterate
+            );
+
+            // Calculate distance
+            cellDistance = Foam::max
+            (
+                cellDistance,
+                tsd.distance(mesh().cellCentres())
+            );
+        }
+    }
+
+    Info<< "Cell distance (min, max): (" << min(cellDistance)
+        << " " << max(cellDistance) << ")" << endl;
+        
+    // Create storage for collection of cells. Assume that almost all of the
+    // cells will be marked to prevent excessive resizing.
+    dynamicLabelList refinementCandidates(mesh().nCells());
+
+    // Create a sorting criterion
+    forAll (cellDistance, cellI)
+    {
+        if
+        (
+             cellDistance[cellI] > -refinementDistance_
+           && cellDistance[cellI] < 0
+        )
+        {
+            // Found a refinement cell
+            refinementCandidates.append(cellI);
+        }
+    }
+
+    // Print out some information
+    Info<< "Selection algorithm " << type() << " selected "
+        << returnReduce(refinementCandidates.size(), sumOp<label>())
+        << " cells as refinement candidates."
+        << endl;
+
+    // Return the list in the Xfer container to prevent copying
+    return refinementCandidates.xfer();
+}
+
+
+Foam::Xfer<Foam::labelList>
+Foam::immersedBoundaryRefinement::unrefinementPointCandidates() const
+{
+    // Mark all points as unrefinement candidates since only split points may be
+    // considered for actual unrefinement and since this refinement criterion
+    // will be usually used in combination with others. VV, 15/Mar/2018.
+
+    // All points are unrefinement candidates
+    dynamicLabelList unrefinementCandidates(mesh().nPoints());
+
+    // Calculate distance to immersed boundary
+    // Note: minimum distance is inside (WET) and implies a live cell
+    scalarField pointDistance(mesh().nPoints(), -GREAT);
+
+    const polyBoundaryMesh& bMesh = mesh().boundaryMesh();
+    
+    const vector span(GREAT, GREAT, GREAT);
+
+    forAll (bMesh, patchI)
+    {
+        if (isA<immersedBoundaryPolyPatch>(bMesh[patchI]))
+        {
+            const immersedBoundaryPolyPatch& ibPatch =
+                refCast<const immersedBoundaryPolyPatch>(bMesh[patchI]);
+
+            // Create a distance function
+            triSurfaceDistance tsd
+            (
+                ibPatch.triSurfSearch(),
+                span,
+                ibPatch.internalFlow(),
+                false                     // Do not iterate
+            );
+
+            // Calculate distance
+            pointDistance = Foam::max
+            (
+                pointDistance,
+                tsd.distance(mesh().points())
+            );
+        }
+    }
+
+    Info<< "Split point distance(min, max): (" << min(pointDistance)
+        << " " << max(pointDistance) << ")" << endl;
+
+    forAll (pointDistance, pointI)
+    {
+        // Coarse live cells next to immersed boundary or dead cells
+        if
+        (
+            pointDistance[pointI] < -coarseningDistance_
+         || pointDistance[pointI] > SMALL
+        )
+        {
+            unrefinementCandidates.append(pointI);
+        }
+    }
+
+    // Print out some information
+    Info<< "Selection algorithm " << type() << " selected "
+        << returnReduce(unrefinementCandidates.size(), sumOp<label>())
+        << " points as unrefinement candidates."
+        << endl;
+
+    // Return the list in the Xfer container to prevent copying
+    return unrefinementCandidates.xfer();
+}
+
+
+// ************************************************************************* //
diff --git a/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryRefinement/immersedBoundaryRefinement.H b/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryRefinement/immersedBoundaryRefinement.H
new file mode 100644
index 000000000..56332af8f
--- /dev/null
+++ b/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryRefinement/immersedBoundaryRefinement.H
@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::immersedBoundaryRefinement
+
+Description
+    Selection of refinement cells based on distance from immersed boundary
+SourceFiles
+    immersedBoundaryRefinement.C
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef immersedBoundaryRefinement_H
+#define immersedBoundaryRefinement_H
+
+#include "refinementSelection.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class immersedBoundaryRefinement Declaration
+\*---------------------------------------------------------------------------*/
+
+class immersedBoundaryRefinement
+:
+    public refinementSelection
+{
+    // Private data
+
+        //- Refinement distance to immersed boundary
+        scalar refinementDistance_;
+
+        //- Coarsening distance to immersed boundary
+        scalar coarseningDistance_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        immersedBoundaryRefinement(const immersedBoundaryRefinement&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const immersedBoundaryRefinement&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("immersedBoundaryRefinement");
+
+
+    // Constructors
+
+        //- Construct from components
+        immersedBoundaryRefinement
+        (
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~immersedBoundaryRefinement();
+
+
+    // Member Functions
+
+        // Selection of refinement/unrefinement candidates
+
+            //- Return transferable list of cells to refine
+            virtual Xfer<labelList> refinementCellCandidates() const;
+
+            //- Return transferable list of split points to unrefine
+            virtual Xfer<labelList> unrefinementPointCandidates() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //