From 2f8c7782f43acdc8cb4471a1530e484de407e59c Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@Germany>
Date: Wed, 2 Dec 2009 14:02:07 +0100
Subject: [PATCH] removed
 dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C.save

- was not touched for a very long time, seems to be left over from
  the pre-git era
---
 .../dynamicRefineFvMesh.C.save                | 1507 -----------------
 1 file changed, 1507 deletions(-)
 delete mode 100644 src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C.save

diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C.save b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C.save
deleted file mode 100644
index 00d2ebe0c6..0000000000
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C.save
+++ /dev/null
@@ -1,1507 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 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 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
-
-\*---------------------------------------------------------------------------*/
-
-#include "SortableList.H"
-#include "dynamicRefineFvMesh.H"
-#include "addToRunTimeSelectionTable.H"
-#include "volFields.H"
-#include "polyTopoChange.H"
-#include "surfaceFields.H"
-#include "fvCFD.H"
-#include "syncTools.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(dynamicRefineFvMesh, 0);
-
-addToRunTimeSelectionTable(dynamicFvMesh, dynamicRefineFvMesh, IOobject);
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-label dynamicRefineFvMesh::count(const PackedList<1>& lst, const unsigned int v)
-{
-    label n = 0;
-
-    forAll(lst, i)
-    {
-        if (lst.get(i) == v)
-        {
-            n++;
-        }
-    }
-    return n;
-}
-
-
-void dynamicRefineFvMesh::calculateProtectedCells
-(
-    PackedList<1>& unrefineableCell
-) const
-{
-    if (protectedCell_.size() == 0)
-    {
-        unrefineableCell.clear();
-        return;
-    }
-
-    const labelList& cellLevel = meshCutter_.cellLevel();
-
-    unrefineableCell = protectedCell_;
-
-    // Get neighbouring cell level
-    labelList neiLevel(nFaces()-nInternalFaces());
-
-    for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
-    {
-        neiLevel[faceI-nInternalFaces()] = cellLevel[faceOwner()[faceI]];
-    }
-    syncTools::swapBoundaryFaceList(*this, neiLevel, false);
-
-
-    while (true)
-    {
-        // Pick up faces on border of protected cells
-        boolList seedFace(nFaces(), false);
-
-        forAll(faceNeighbour(), faceI)
-        {
-            label own = faceOwner()[faceI];
-            bool ownProtected = (unrefineableCell.get(own) == 1);
-            label nei = faceNeighbour()[faceI];
-            bool neiProtected = (unrefineableCell.get(nei) == 1);
-
-            if (ownProtected && (cellLevel[nei] > cellLevel[own]))
-            {
-                seedFace[faceI] = true;
-            }
-            else if (neiProtected && (cellLevel[own] > cellLevel[nei]))
-            {
-                seedFace[faceI] = true;
-            }
-        }
-        for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
-        {
-            label own = faceOwner()[faceI];
-            bool ownProtected = (unrefineableCell.get(own) == 1);
-            if
-            (
-                ownProtected
-             && (neiLevel[faceI-nInternalFaces()] > cellLevel[own])
-            )
-            {
-                seedFace[faceI] = true;
-            }
-        }
-
-        syncTools::syncFaceList(*this, seedFace, orEqOp<bool>(), false);
-
-
-        // Extend unrefineableCell
-        bool hasExtended = false;
-
-        for (label faceI = 0; faceI < nInternalFaces(); faceI++)
-        {
-            if (seedFace[faceI])
-            {
-                label own = faceOwner()[faceI];
-                if (unrefineableCell.get(own) == 0)
-                {
-                    unrefineableCell.set(own, 1);
-                    hasExtended = true;
-                }
-
-                label nei = faceNeighbour()[faceI];
-                if (unrefineableCell.get(nei) == 0)
-                {
-                    unrefineableCell.set(nei, 1);
-                    hasExtended = true;
-                }
-            }
-        }
-        for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
-        {
-            if (seedFace[faceI])
-            {
-                label own = faceOwner()[faceI];
-                if (unrefineableCell.get(own) == 0)
-                {
-                    unrefineableCell.set(own, 1);
-                    hasExtended = true;
-                }
-            }
-        }
-
-        if (!returnReduce(hasExtended, orOp<bool>()))
-        {
-            break;
-        }
-    }
-}
-
-
-void dynamicRefineFvMesh::readDict()
-{
-    dictionary refineDict
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "dynamicMeshDict",
-                time().constant(),
-                *this,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            )
-        ).subDict(typeName + "Coeffs")
-    );
-
-    correctFluxes_ = List<Pair<word> >(refineDict.lookup("correctFluxes"));
-
-    if (refineDict.found("projectSurfaces"))
-    {
-        surfacesPtr_.reset
-        (
-            new triSurfaceMeshes
-            (
-                IOobject
-                (
-                    "abc",                              // dummy name
-                    time().constant(),                  // directory
-                    "triSurface",                       // instance
-                    time(),                             // registry
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                ),
-                refineDict.lookup("projectSurfaces")
-            )
-        );
-
-        wordList projectPatchNames(refineDict.lookup("projectPatches"));
-
-        projectSurfaces_.setSize(boundaryMesh().size());
-
-        forAll(projectPatchNames, i)
-        {
-            label patchI = boundaryMesh().findPatchID(projectPatchNames[i]);
-
-            if (patchI == -1)
-            {
-                FatalErrorIn("dynamicRefineFvMesh::readDict()")
-                    << "Illegal projectPatch " << projectPatchNames[i] << endl
-                    << "Valid patch names are " << boundaryMesh().names()
-                    << exit(FatalError);
-            }
-
-            // For now have patches always refer to all surfaces
-            projectSurfaces_[patchI] = identity(surfacesPtr_().size());
-        }
-
-        matchDistSqr_ = Foam::sqr
-        (
-            readScalar(refineDict.lookup("projectDistance"))
-        );
-    }
-}
-
-
-// Project new vertices onto surface
-void dynamicRefineFvMesh::project
-(
-    const mapPolyMesh& map,
-    pointField& points
-) const
-{
-    if (!surfacesPtr_.valid())
-    {
-        return;
-    }
-
-
-    const polyBoundaryMesh& patches = boundaryMesh();
-    const labelList& pointMap = map.pointMap();
-    const labelList& reversePointMap = map.reversePointMap();
-
-    pointIndexHit hitInfo;
-
-    forAll(patches, patchI)
-    {
-        if (projectSurfaces_[patchI].size() > 0)
-        {
-            const labelList& meshPoints = patches[patchI].meshPoints();
-
-            forAll(meshPoints, i)
-            {
-                label pointI = meshPoints[i];
-
-                label oldPointI = pointMap[pointI];
-
-                if (oldPointI < 0 || reversePointMap[oldPointI] != pointI)
-                {
-                    // New point. Project.
-                    surfacesPtr_().findNearest
-                    (
-                        projectSurfaces_[patchI],
-                        points[pointI],
-                        matchDistSqr_,
-                        hitInfo
-                    );
-
-                    if (hitInfo.hit())
-                    {
-                        points[pointI] = hitInfo.hitPoint();
-                    }
-                }
-            }
-        }
-    }
-}
-
-
-// Refines cells, maps fields and recalculates (an approximate) flux
-autoPtr<mapPolyMesh> 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())
-            {
-                FatalErrorIn("dynamicRefineFvMesh::refine")
-                    << "New internal face:" << faceI
-                    << " fc:" << faceCentres()[faceI]
-                    << " originates from boundary oldFace:" << oldFaceI
-                    << abort(FatalError);
-            }
-        }
-    }
-
-
-    // Update fields
-    updateMesh(map);
-
-    // Move mesh
-    /*
-    pointField newPoints;
-    if (map().hasMotionPoints())
-    {
-        newPoints = map().preMotionPoints();
-    }
-    else
-    {
-        newPoints = points();
-    }
-    project(map, newPoints);
-    movePoints(newPoints);
-    */
-
-    // Correct the flux for modified/added faces. All the faces which only
-    // have been renumbered will already have been handled by the mapping.
-    {
-        const labelList& faceMap = map().faceMap();
-        const labelList& reverseFaceMap = map().reverseFaceMap();
-
-        // Storage for any master faces. These will be the original faces
-        // on the coarse cell that get split into four (or rather the
-        // master face gets modified and three faces get added from the master)
-        labelHashSet masterFaces(4*cellsToRefine.size());
-
-        forAll(faceMap, faceI)
-        {
-            label oldFaceI = faceMap[faceI];
-
-            if (oldFaceI >= 0)
-            {
-                label masterFaceI = reverseFaceMap[oldFaceI];
-
-                if (masterFaceI < 0)
-                {
-                    FatalErrorIn
-                    (
-                        "dynamicRefineFvMesh::refine(const labelList&)"
-                    )   << "Problem: should not have removed faces"
-                        << " when refining."
-                        << nl << "face:" << faceI << abort(FatalError);
-                }
-                else if (masterFaceI != faceI)
-                {
-                    masterFaces.insert(masterFaceI);
-                }
-            }
-        }
-        if (debug)
-        {
-            Info<< "Found " << returnReduce(masterFaces.size(), sumOp<label>())
-                << " split faces " << endl;
-        }
-
-        forAll(correctFluxes_, i)
-        {
-            if (debug)
-            {
-                Info<< "Mapping flux " << correctFluxes_[i][0]
-                    << " using interpolated flux " << correctFluxes_[i][1]
-                    << endl;
-            }
-            surfaceScalarField& phi = const_cast<surfaceScalarField&>
-            (
-                lookupObject<surfaceScalarField>(correctFluxes_[i][0])
-            );
-            surfaceScalarField phiU =
-                fvc::interpolate
-                (
-                    lookupObject<volVectorField>(correctFluxes_[i][1])
-                )
-              & Sf();
-
-            // Recalculate new internal faces.
-            for (label faceI = 0; faceI < nInternalFaces(); faceI++)
-            {
-                label oldFaceI = faceMap[faceI];
-
-                if (oldFaceI == -1)
-                {
-                    // Inflated/appended
-                    phi[faceI] = phiU[faceI];
-                }
-                else if (reverseFaceMap[oldFaceI] != faceI)
-                {
-                    // face-from-masterface
-                    phi[faceI] = phiU[faceI];
-                }
-            }
-
-            // Recalculate new boundary faces.
-            forAll(phi.boundaryField(), patchI)
-            {
-                fvsPatchScalarField& patchPhi = phi.boundaryField()[patchI];
-                const fvsPatchScalarField& patchPhiU =
-                    phiU.boundaryField()[patchI];
-
-                label faceI = patchPhi.patch().patch().start();
-
-                forAll(patchPhi, i)
-                {
-                    label oldFaceI = faceMap[faceI];
-
-                    if (oldFaceI == -1)
-                    {
-                        // Inflated/appended
-                        patchPhi[i] = patchPhiU[i];
-                    }
-                    else if (reverseFaceMap[oldFaceI] != faceI)
-                    {
-                        // face-from-masterface
-                        patchPhi[i] = patchPhiU[i];
-                    }
-
-                    faceI++;
-                }
-            }
-
-            // Update master faces
-            forAllConstIter(labelHashSet, masterFaces, iter)
-            {
-                label faceI = iter.key();
-
-                if (isInternalFace(faceI))
-                {
-                    phi[faceI] = phiU[faceI];
-                }
-                else
-                {
-                    label patchI = boundaryMesh().whichPatch(faceI);
-                    label i = faceI - boundaryMesh()[patchI].start();
-
-                    const fvsPatchScalarField& patchPhiU =
-                        phiU.boundaryField()[patchI];
-
-                    fvsPatchScalarField& patchPhi =
-                        phi.boundaryField()[patchI];
-
-                    patchPhi[i] = patchPhiU[i];
-                }
-            }
-        }
-    }            
-
-
-
-    // Update numbering of cells/vertices.
-    meshCutter_.updateMesh(map);
-
-    // Update numbering of protectedCell_
-    PackedList<1> newProtectedCell(nCells(), 0);
-
-    forAll(newProtectedCell, cellI)
-    {
-        label oldCellI = map().cellMap()[cellI];
-        newProtectedCell.set(cellI, protectedCell_.get(oldCellI));
-    }
-    protectedCell_.transfer(newProtectedCell);
-
-    // Debug: Check refinement levels (across faces only)
-    meshCutter_.checkRefinementLevels(-1, labelList(0));
-
-    return map;
-}
-
-
-// Combines previously split cells, maps fields and recalculates
-// (an approximate) flux
-autoPtr<mapPolyMesh> dynamicRefineFvMesh::unrefine
-(
-    const labelList& splitPoints
-)
-{
-    polyTopoChange 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 mesh.
-    //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
-    autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
-
-    Info<< "Unrefined from "
-        << returnReduce(map().nOldCells(), sumOp<label>())
-        << " to " << globalData().nTotalCells() << " cells."
-        << endl;
-
-    // Update fields
-    updateMesh(map);
-
-
-    // Move mesh
-    /*
-    pointField newPoints;
-    if (map().hasMotionPoints())
-    {
-        newPoints = map().preMotionPoints();
-    }
-    else
-    {
-        newPoints = points();
-    }
-    project(map, newPoints);
-    movePoints(newPoints);
-    */
-
-    // Correct the flux for modified faces.
-    {
-        const labelList& reversePointMap = map().reversePointMap();
-        const labelList& reverseFaceMap = map().reverseFaceMap();
-
-        forAll(correctFluxes_, i)
-        {
-            if (debug)
-            {
-                Info<< "Mapping flux " << correctFluxes_[i][0]
-                    << " using interpolated flux " << correctFluxes_[i][1]
-                    << endl;
-            }
-            surfaceScalarField& phi = const_cast<surfaceScalarField&>
-            (
-                lookupObject<surfaceScalarField>(correctFluxes_[i][0])
-            );
-            surfaceScalarField phiU =
-                fvc::interpolate
-                (
-                    lookupObject<volVectorField>(correctFluxes_[i][1])
-                )
-              & Sf();
-
-            forAllConstIter(Map<label>, faceToSplitPoint, iter)
-            {
-                label oldFaceI = iter.key();
-                label oldPointI = iter();
-
-                if (reversePointMap[oldPointI] < 0)
-                {
-                    // midpoint was removed. See if face still exists.
-                    label faceI = reverseFaceMap[oldFaceI];
-
-                    if (faceI >= 0)
-                    {
-                        if (isInternalFace(faceI))
-                        {
-                            phi[faceI] = phiU[faceI];
-                        }
-                        else
-                        {
-                            label patchI = boundaryMesh().whichPatch(faceI);
-                            label i = faceI - boundaryMesh()[patchI].start();
-
-                            const fvsPatchScalarField& patchPhiU =
-                                phiU.boundaryField()[patchI];
-
-                            fvsPatchScalarField& patchPhi =
-                                phi.boundaryField()[patchI];
-
-                            patchPhi[i] = patchPhiU[i];
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-
-    // Update numbering of cells/vertices.
-    meshCutter_.updateMesh(map);
-
-    // Update numbering of protectedCell_
-    PackedList<1> newProtectedCell(nCells(), 0);
-
-    forAll(newProtectedCell, cellI)
-    {
-        label oldCellI = map().cellMap()[cellI];
-        if (oldCellI >= 0)
-        {
-            newProtectedCell.set(cellI, protectedCell_.get(oldCellI));
-        }
-    }
-    protectedCell_.transfer(newProtectedCell);
-
-    // Debug: Check refinement levels (across faces only)
-    meshCutter_.checkRefinementLevels(-1, labelList(0));
-
-    return map;
-}
-
-
-//
-// Select based on cell values
-// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-//
-
-labelList dynamicRefineFvMesh::selectRefineCells
-(
-    const label maxCells,
-    const label maxRefinement,
-    const volScalarField& cellStatus
-) const
-{
-    DynamicList<label> candidateCells(nCells()/20);
-
-    // Every refined cell causes 7 extra cells
-    label nTotToRefine = (maxCells - globalData().nTotalCells()) / 7;
-
-    // Pick up just enough cells. In cellLevel order so coarsest cells
-    // are most likely to be refined if we hit the limit.
-    const labelList& cellLevel = meshCutter_.cellLevel();
-
-    // Mark cells that cannot be refined since they would trigger refinement
-    // of protected cells (since 2:1 cascade)
-    PackedList<1> unrefineableCell;
-    calculateProtectedCells(unrefineableCell);
-
-    label maxLevel = gMax(cellLevel);
-
-    for (label level = 0; level < maxRefinement && level <= maxLevel; level++)
-    {
-        forAll(cellStatus, cellI)
-        {
-            if
-            (
-                cellStatus[cellI] == 1
-             && cellLevel[cellI] == level
-             && (
-                    unrefineableCell.size() == 0
-                 || unrefineableCell.get(cellI) == 0
-                )
-            )
-            {
-                // Emergency stop: even just local refinement causes
-                // surpassing of max.
-                if (candidateCells.size() >= nTotToRefine)
-                {
-                    break;
-                }
-                candidateCells.append(cellI);
-            }
-        }
-
-        label nCandidates = returnReduce(candidateCells.size(), sumOp<label>());
-
-        if (nCandidates >= nTotToRefine)
-        {
-            break;
-        }
-    }
-    candidateCells.shrink();
-
-    // Guarantee 2:1 refinement across processor patches.
-    labelList cellsToRefine
-    (
-        meshCutter_.consistentRefinement
-        (
-            candidateCells,
-            true               // Add to set to guarantee 2:1
-        )
-    );
-
-
-    // Debug
-    if (debug)
-    {
-        forAll(cellsToRefine, i)
-        {
-            label cellI = cellsToRefine[i];
-
-            if (protectedCell_.get(cellI) == 1)
-            {
-                FatalErrorIn("dynamicRefineFvMesh::selectRefineCells")
-                    << "Problem:" << cellI
-                    << " should be protected from refinement."
-                    << abort(FatalError);
-            }
-        }
-    }
-
-    return cellsToRefine;
-}
-
-
-labelList dynamicRefineFvMesh::selectUnrefinePoints
-(
-    const volScalarField& cellStatus
-) const
-{
-    // All points that can be unrefined
-    const labelList splitPoints(meshCutter_.getSplitPoints());
-
-    DynamicList<label> newSplitPoints(splitPoints.size());
-
-    forAll(splitPoints, i)
-    {
-        label pointI = splitPoints[i];
-
-        const labelList& pCells = pointCells()[pointI];
-
-        // Check that all pCells are marked for unrefinement
-        bool allUnrefine = true;
-
-        forAll(pCells, pCellI)
-        {
-            if (cellStatus[pCells[pCellI]] != -1)
-            {
-                allUnrefine = false;
-                break;
-            }
-        }
-
-        if (allUnrefine)
-        {
-            newSplitPoints.append(pointI);
-        }
-    }
-    newSplitPoints.shrink();
-
-    return meshCutter_.consistentUnrefinement
-    (
-        newSplitPoints,
-        false
-    );
-}
-
-
-//
-// Select based on point values
-// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-//
-
-labelList dynamicRefineFvMesh::selectRefineCells
-(
-    const label maxCells,
-    const label maxRefinement,
-    const PackedList<1>& refinePoint
-) const
-{
-    // Every refined cell causes 7 extra cells
-    label nTotToRefine = (maxCells - globalData().nTotalCells()) / 7;
-
-    // Mark cells that cannot be refined since they would trigger refinement
-    // of protected cells (since 2:1 cascade)
-    PackedList<1> unrefineableCell;
-    calculateProtectedCells(unrefineableCell);
-
-    // Pick up just enough cells. In cellLevel order so coarsest cells
-    // are most likely to be refined if we hit the limit.
-    const labelList& cellLevel = meshCutter_.cellLevel();
-
-    labelList candidates;
-    {
-        labelHashSet candidateSet(nCells()/20);
-
-        forAll(refinePoint, pointI)
-        {
-            if (refinePoint.get(pointI) == 1)
-            {
-                const labelList& pCells = pointCells()[pointI];
-
-                forAll(pCells, i)
-                {
-                    label cellI = pCells[i];
-
-                    if
-                    (
-                        cellLevel[cellI] < maxRefinement
-                     && (
-                            unrefineableCell.size() == 0
-                         || unrefineableCell.get(cellI) == 0
-                        )
-                    )
-                    {
-                        candidateSet.insert(cellI);
-
-                        // Emergency stop: even just local refinement causes
-                        // surpassing of max.
-                        if (candidateSet.size() >= nTotToRefine)
-                        {
-                            break;
-                        }
-                    }
-                }
-            }
-        }
-
-        candidates = candidateSet.toc();
-    }
-
-    if (returnReduce(candidates.size(), sumOp<label>()) > nTotToRefine)
-    {
-        // Truncate. Refine lower levels first until reach max.
-
-        // Sort according to ascending cell level.
-        {
-            SortableList<label> sortedLevel(renumber(cellLevel, candidates));
-
-            const labelList& newToOld = sortedLevel.indices();
-
-            labelList sortedCells(candidates.size());
-            forAll(sortedCells, i)
-            {
-                sortedCells[i] = candidates[newToOld[i]];
-            }
-            candidates.transfer(sortedCells);
-        }
-
-        // Pick up per cell level until we reach total number of cells.
-        label maxLevel = gMax(cellLevel);
-
-        label i = 0;
-
-        for
-        (
-            label level = 0;
-            level < maxRefinement && level <= maxLevel;
-            level++
-        )
-        {
-            for
-            (
-                ;
-                i < candidates.size() && cellLevel[candidates[i]] == level;
-                i++
-            );
-
-            if (returnReduce(i, sumOp<label>()) >= nTotToRefine)
-            {
-                break;
-            }
-        }
-
-        candidates.setSize(i);
-    }
-
-    // Guarantee 2:1 refinement across processor patches.
-    labelList cellsToRefine
-    (
-        meshCutter_.consistentRefinement
-        (
-            candidates,
-            true               // Add to set to guarantee 2:1
-        )
-    );
-
-
-    /// Check if all correct
-    if (debug)
-    {
-        forAll(cellsToRefine, i)
-        {
-            label cellI = cellsToRefine[i];
-
-            if (protectedCell_.get(cellI) == 1)
-            {
-                FatalErrorIn("dynamicRefineFvMesh::selectRefineCells")
-                    << "Problem:" << cellI
-                    << " should be protected from refinement."
-                    << abort(FatalError);
-            }
-        }
-    }
-
-    return cellsToRefine;
-}
-
-
-labelList dynamicRefineFvMesh::selectUnrefinePoints
-(
-    const PackedList<1>& markedPoint
-) const
-{
-    // All points that can be unrefined
-    const labelList splitPoints(meshCutter_.getSplitPoints());
-
-    DynamicList<label> newSplitPoints(splitPoints.size());
-
-    forAll(splitPoints, i)
-    {
-        label pointI = splitPoints[i];
-
-        //if (gamma[pointI] < minGamma || gamma[pointI] > maxGamma)
-        if (markedPoint.get(pointI) == 0)
-        {
-            newSplitPoints.append(pointI);
-        }
-    }
-    newSplitPoints.shrink();
-
-    return meshCutter_.consistentUnrefinement
-    (
-        newSplitPoints,
-        false
-    );
-}
-
-
-//
-// Select with marker array.
-//
-void dynamicRefineFvMesh::markRefinePoints
-(
-    const scalarField& gamma,
-    const scalar minGamma,
-    const scalar maxGamma,
-    PackedList<1>& refinePoint
-) const
-{
-    forAll(gamma, pointI)
-    {
-        if (gamma[pointI] > minGamma && gamma[pointI] < maxGamma)
-        {
-            refinePoint.set(pointI, 1);
-        }
-    }
-}
-
-
-void dynamicRefineFvMesh::extendMarkedPoints
-(
-    PackedList<1>& markedPoint
-) const
-{
-    // Mark faces using any marked point
-
-    PackedList<1> markedFace(nFaces(), 0);
-
-    forAll(faces(), faceI)
-    {
-        const face& f = faces()[faceI];
-
-        forAll(f, fp)
-        {
-            if (markedPoint.get(f[fp]) == 1)
-            {
-                markedFace.set(faceI, 1);
-                break;
-            }
-        }
-    }
-
-    // Mark all points on marked faces
-
-    forAll(markedFace, faceI)
-    {
-        if (markedFace.get(faceI) == 1)
-        {
-            const face& f = faces()[faceI];
-
-            forAll(f, fp)
-            {
-                markedPoint.set(f[fp], 1);
-            }
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
-:
-    dynamicFvMesh(io),
-    meshCutter_(*this),
-    surfacesPtr_(NULL),
-    projectSurfaces_(0),
-    matchDistSqr_(-GREAT),
-    nRefinementIterations_(0),
-    protectedCell_(nCells(), 0)
-{
-    // Read the optional projection data
-    readDict();
-
-    const labelList& cellLevel = meshCutter_.cellLevel();
-    const labelList& pointLevel = meshCutter_.pointLevel();
-
-    // Set cells that should not be refined.
-    // This is currently any cell which does not have 8 anchor points or
-    // uses any face which does not have 4 anchor points.
-    // Note: do not use cellPoint addressing
-
-    // Count number of points <= cellLevel
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelList nAnchors(nCells(), 0);
-
-    label nProctected = 0;
-
-    forAll(pointCells(), pointI)
-    {
-        const labelList& pCells = pointCells()[pointI];
-
-        forAll(pCells, i)
-        {
-            label cellI = pCells[i];
-
-            if (protectedCell_.get(cellI) == 0)
-            {
-                if (pointLevel[pointI] <= cellLevel[cellI])
-                {
-                    nAnchors[cellI]++;
-
-                    if (nAnchors[cellI] > 8)
-                    {
-                        protectedCell_.set(cellI, 1);
-                        nProctected++;
-                    }
-                }
-            }
-        }
-    }
-
-
-    // Count number of points <= faceLevel
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Bit tricky since proc face might be one more refined than the owner since
-    // the coupled one is refined.
-
-    {
-        labelList neiLevel(nFaces());
-
-        for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
-        {
-            neiLevel[faceI] = cellLevel[faceOwner()[faceI]];
-        }
-        syncTools::swapFaceList(*this, neiLevel, false);
-
-
-        boolList protectedFace(nFaces(), false);
-
-        forAll(faceOwner(), faceI)
-        {
-            label faceLevel = max
-            (
-                cellLevel[faceOwner()[faceI]],
-                neiLevel[faceI]
-            );
-
-            const face& f = faces()[faceI];
-
-            label nAnchors = 0;
-
-            forAll(f, fp)
-            {
-                if (pointLevel[f[fp]] <= faceLevel)
-                {
-                    nAnchors++;
-
-                    if (nAnchors > 4)
-                    {
-                        protectedFace[faceI] = true;
-                        break;
-                    }
-                }
-            }
-        }
-
-        syncTools::syncFaceList
-        (
-            *this,
-            protectedFace,
-            orEqOp<bool>(),
-            false
-        );
-
-        for (label faceI = 0; faceI < nInternalFaces(); faceI++)
-        {
-            if (protectedFace[faceI])
-            {
-                protectedCell_.set(faceOwner()[faceI], 1);
-                nProctected++;
-                protectedCell_.set(faceNeighbour()[faceI], 1);
-                nProctected++;
-            }
-        }
-        for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
-        {
-            if (protectedFace[faceI])
-            {
-                protectedCell_.set(faceOwner()[faceI], 1);
-                nProctected++;
-            }
-        }
-    }
-
-    if (returnReduce(nProctected, sumOp<label>()) == 0)
-    {
-        protectedCell_.clear();
-    }
-
-    Info<< "dynamicRefineFvMesh constructed with:" << nl
-        << "    cellLevel min:" << gMin(cellLevel)
-        << "  max:" << gMax(cellLevel) << endl;
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-dynamicRefineFvMesh::~dynamicRefineFvMesh()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool dynamicRefineFvMesh::update()
-{
-    // 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"));
-
-    if (refineInterval == 0)
-    {
-        return false;
-    }
-    else if (refineInterval < 0)
-    {
-        FatalErrorIn("dynamicRefineFvMesh::update()")
-            << "Illegal refineInterval " << refineInterval << nl
-            << "The refineInterval setting in the dynamicMeshDict should"
-            << " be >= 1." << nl
-            << exit(FatalError);
-    }
-
-
-
-    bool hasChanged = false;
-
-
-    // Note: cannot refine at time 0 since no V0 present since mesh not
-    //       moved yet.
-
-    if (time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0)
-    {
-        label maxCells = readLabel(refineDict.lookup("maxCells"));
-
-        if (maxCells <= 0)
-        {
-            FatalErrorIn("dynamicRefineFvMesh::update()")
-                << "Illegal maximum number of cells " << maxCells << nl
-                << "The maxCells setting in the dynamicMeshDict should"
-                << " be > 0." << nl
-                << exit(FatalError);
-        }
-
-        label maxRefinement = readLabel(refineDict.lookup("maxRefinement"));
-
-        if (maxRefinement <= 0)
-        {
-            FatalErrorIn("dynamicRefineFvMesh::update()")
-                << "Illegal maximum refinement level " << maxRefinement << nl
-                << "The maxCells setting in the dynamicMeshDict should"
-                << " be > 0." << nl
-                << exit(FatalError);
-        }
-
-        const volScalarField& gamma = lookupObject<volScalarField>
-        (
-            refineDict.lookup("field")
-        );
-
-
-        //// Based on cell indicator field (1=refine , 0, -1=unrefine)
-        //{
-        //    labelList cellsToRefine
-        //    (
-        //        selectRefineCells
-        //        (
-        //            maxCells,
-        //            maxRefinement,
-        //            gamma
-        //        )
-        //    );
-        //
-        //    label nCellsToRefine = returnReduce
-        //    (
-        //        cellsToRefine.size(), sumOp<label>()
-        //    );
-        //
-        //    Info<< "Selected " << nCellsToRefine
-        //        << " cells to refine." << endl;
-        //
-        //    if (nCellsToRefine > 0)
-        //    {
-        //        // Refine/update mesh and map fields
-        //        autoPtr<mapPolyMesh> map = refine(cellsToRefine);
-        //
-        //        hasChanged = true;
-        //    }
-        //}
-        //{
-        //    labelList pointsToUnrefine(selectUnrefinePoints(gamma));
-        //
-        //    label nSplitPoints = returnReduce
-        //    (
-        //        pointsToUnrefine.size(),
-        //        sumOp<label>()
-        //    );
-        //
-        //    Info<< "Selected " << nSplitPoints
-        //        << " points where to merge surrounding cells." << endl;
-        //
-        //    if (nSplitPoints > 0)
-        //    {
-        //        // Refine/update mesh
-        //        unrefine(pointsToUnrefine);
-        //
-        //        hasChanged = true;
-        //    }
-        //}
-
-
-
-        const scalar minLevel = readScalar(refineDict.lookup("minLevel"));
-        const scalar maxLevel = readScalar(refineDict.lookup("maxLevel"));
-        const label nBufferLayers =
-            readLabel(refineDict.lookup("nBufferLayers"));
-
-        // Points marked for refinement
-        PackedList<1> refinePoint(nPoints(), 0);
-
-        {
-            // Do naive interpolation to get point values
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-            scalarField pointGamma(nPoints(), 0.0);
-            labelList nPointGamma(nPoints(), 0);
-
-            forAll(pointCells(), pointI)
-            {
-                const labelList& pCells = pointCells()[pointI];
-
-                forAll(pCells, i)
-                {
-                    pointGamma[pointI] += gamma[pCells[i]];
-                    nPointGamma[pointI]++;
-                }
-            }
-
-            syncTools::syncPointList
-            (
-                *this,
-                pointGamma,
-                plusEqOp<scalar>(), // combine op
-                0.0,                // null value
-                false               // no separation
-            );
-            syncTools::syncPointList
-            (
-                *this,
-                nPointGamma,
-                plusEqOp<label>(),  // combine op
-                0,                  // null value
-                false               // no separation
-            );
-
-            forAll(pointGamma, pointI)
-            {
-                pointGamma[pointI] /= nPointGamma[pointI];
-            }
-
-            // Mark all points with gamma witin refine range
-            markRefinePoints
-            (
-                pointGamma,
-                minLevel,
-                maxLevel,
-                refinePoint
-            );
-        }
-
-
-
-        {
-            labelList cellsToRefine
-            (
-                selectRefineCells
-                (
-                    maxCells,
-                    maxRefinement,
-                    refinePoint
-                )
-            );
-
-            label nCellsToRefine = returnReduce
-            (
-                cellsToRefine.size(), sumOp<label>()
-            );
-
-            if (nCellsToRefine > 0)
-            {
-                // Refine/update mesh and map fields
-                autoPtr<mapPolyMesh> map = refine(cellsToRefine);
-
-                // Update the interpolated field such that newly created points
-                // don't get unrefined.
-                const labelList& pointMap = map().pointMap();
-                const labelList& reversePointMap = map().reversePointMap();
-
-                // Map refinePoint. Set new points to a refine level
-                {
-                    PackedList<1> newRefinePoint(pointMap.size());
-
-                    forAll(pointMap, pointI)
-                    {
-                        label oldPointI = pointMap[pointI];
-
-                        if (oldPointI < 0)
-                        {
-                            newRefinePoint.set(pointI, 1);
-                        }
-                        else if (reversePointMap[oldPointI] != pointI)
-                        {
-                            newRefinePoint.set(pointI, 1);
-                        }
-                        else
-                        {
-                            newRefinePoint.set(pointI, refinePoint.get(pointI));
-                        }
-                    }
-                    refinePoint.transfer(newRefinePoint);
-                }
-
-
-                hasChanged = true;
-            }
-        }
-
-        // Extend with a buffer layer to prevent neighbouring points being
-        // unrefined.
-        for (label i = 0; i < nBufferLayers; i++)
-        {
-            extendMarkedPoints(refinePoint);
-        }
-
-        {
-            // Select unrefineable points that are not marked in refinePoint
-            labelList pointsToUnrefine
-            (
-                selectUnrefinePoints
-                (
-                    refinePoint
-                )
-            );
-
-            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 refinementHistory.
-            const_cast<refinementHistory&>(meshCutter().history()).compact();
-        }
-        nRefinementIterations_++;
-    }
-
-    mesh.changing(hasChanged);
-
-    return hasChanged;
-}
-
-
-bool dynamicRefineFvMesh::writeObject
-(
-    IOstream::streamFormat fmt,
-    IOstream::versionNumber ver,
-    IOstream::compressionType cmp
-) const
-{
-    // Force refinement data to go to the current time directory.
-    const_cast<hexRef8&>(meshCutter_).setInstance(time().timeName());
-
-    bool writeOk =
-        dynamicFvMesh::writeObjects(fmt, ver, cmp)
-     && meshCutter_.write();
-
-    if (debug)
-    {
-        volScalarField scalarCellLevel
-        (
-            IOobject
-            (
-                "cellLevel",
-                time().timeName(),
-                *this,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            *this,
-            dimensionedScalar("level", dimless, 0)
-        );
-
-        const labelList& cellLevel = meshCutter_.cellLevel();
-
-        forAll(cellLevel, cellI)
-        {
-            scalarCellLevel[cellI] = cellLevel[cellI];
-        }
-
-        Info<< "dynamicRefineFvMesh : Writing cellLevel ..." << endl;
-
-        writeOk = writeOk && scalarCellLevel.write();
-    }
-
-    return writeOk;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //