Skip to content
Snippets Groups Projects
snappyRefineDriver.C 115 KiB
Newer Older
  • Learn to ignore specific revisions
  • /*---------------------------------------------------------------------------*\
      =========                 |
      \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       \\    /   O peration     |
    
    OpenFOAM bot's avatar
    OpenFOAM bot committed
        \\  /    A nd           | www.openfoam.com
    
         \\/     M anipulation  |
    -------------------------------------------------------------------------------
    
    OpenFOAM bot's avatar
    OpenFOAM bot committed
        Copyright (C) 2011-2015 OpenFOAM Foundation
    
        Copyright (C) 2015-2025 OpenCFD Ltd.
    
    -------------------------------------------------------------------------------
    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/>.
    
    \*---------------------------------------------------------------------------*/
    
    mattijs's avatar
    mattijs committed
    #include "meshRefinement.H"
    
    #include "fvMesh.H"
    #include "Time.H"
    #include "cellSet.H"
    #include "syncTools.H"
    #include "refinementParameters.H"
    #include "refinementSurfaces.H"
    
    #include "refinementFeatures.H"
    
    #include "shellSurfaces.H"
    
    mattijs's avatar
    mattijs committed
    #include "mapDistributePolyMesh.H"
    
    #include "unitConversion.H"
    
    #include "snapParameters.H"
    
    #include "localPointRegion.H"
    
    #include "IOmanip.H"
    #include "labelVector.H"
    
    #include "profiling.H"
    
    #include "searchableSurfaces.H"
    
    #include "fvMeshSubset.H"
    #include "interpolationTable.H"
    
    #include "snappyVoxelMeshDriver.H"
    
    #include "regionSplit.H"
    #include "removeCells.H"
    
    #include "addPatchCellLayer.H"
    
    
    // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
    
    namespace Foam
    {
    
        defineTypeNameAndDebug(snappyRefineDriver, 0);
    
    } // End namespace Foam
    
    
    // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
    
    
    Foam::snappyRefineDriver::snappyRefineDriver
    
    (
        meshRefinement& meshRefiner,
        decompositionMethod& decomposer,
        fvMeshDistribute& distributor,
    
        const labelUList& globalToMasterPatch,
    
        const labelUList& globalToSlavePatch,
    
        coordSetWriter& setFormatter,
    
        refPtr<surfaceWriter>& surfFormatter,
    
    )
    :
        meshRefiner_(meshRefiner),
        decomposer_(decomposer),
        distributor_(distributor),
    
        globalToMasterPatch_(globalToMasterPatch),
    
        globalToSlavePatch_(globalToSlavePatch),
    
        setFormatter_(setFormatter),
    
        surfFormatter_(surfFormatter),
    
    {}
    
    
    // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
    
    
    Foam::label Foam::snappyRefineDriver::featureEdgeRefine
    
    (
        const refinementParameters& refineParams,
        const label maxIter,
        const label minRefine
    )
    {
    
        if (refineParams.minRefineCells() == -1)
        {
            // Special setting to be able to restart shm on meshes with inconsistent
            // cellLevel/pointLevel
            return 0;
        }
    
    
        addProfiling(edge, "snappyHexMesh::refine::edge");
    
        const fvMesh& mesh = meshRefiner_.mesh();
    
        label iter = 0;
    
    
        if (meshRefiner_.features().size() && maxIter > 0)
    
        {
            for (; iter < maxIter; iter++)
            {
                Info<< nl
                    << "Feature refinement iteration " << iter << nl
                    << "------------------------------" << nl
                    << endl;
    
                labelList candidateCells
                (
                    meshRefiner_.refineCandidates
                    (
    
                        refineParams.curvature(),
    
                        refineParams.planarAngle(),
    
    
                        true,               // featureRefinement
    
                        false,              // featureDistanceRefinement
    
                        false,              // internalRefinement
                        false,              // surfaceRefinement
                        false,              // curvatureRefinement
    
                        false,              // smallFeatureRefinement
    
                        false,              // gapRefinement
    
                        false,              // spreadGapSize
    
                        refineParams.maxGlobalCells(),
                        refineParams.maxLocalCells()
                    )
                );
                labelList cellsToRefine
                (
                    meshRefiner_.meshCutter().consistentRefinement
                    (
                        candidateCells,
                        true
                    )
                );
                Info<< "Determined cells to refine in = "
                    << mesh.time().cpuTimeIncrement() << " s" << endl;
    
    
    
    
                const label nCellsToRefine =
                    returnReduce(cellsToRefine.size(), sumOp<label>());
    
    
                Info<< "Selected for feature refinement : " << nCellsToRefine
                    << " cells (out of " << mesh.globalData().nTotalCells()
                    << ')' << endl;
    
                if (nCellsToRefine <= minRefine)
                {
                    Info<< "Stopping refining since too few cells selected."
                        << nl << endl;
                    break;
                }
    
    
                if (debug > 0)
                {
                    const_cast<Time&>(mesh.time())++;
                }
    
    
                if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
    
                {
                    meshRefiner_.balanceAndRefine
                    (
                        "feature refinement iteration " + name(iter),
                        decomposer_,
                        distributor_,
                        cellsToRefine,
    
                        refineParams.maxLoadUnbalance(),
                        refineParams.maxCellUnbalance()
    
                    );
                }
                else
                {
                    meshRefiner_.refineAndBalance
                    (
                        "feature refinement iteration " + name(iter),
                        decomposer_,
                        distributor_,
                        cellsToRefine,
    
                        refineParams.maxLoadUnbalance(),
                        refineParams.maxCellUnbalance()
    
    Foam::label Foam::snappyRefineDriver::smallFeatureRefine
    
    (
        const refinementParameters& refineParams,
        const label maxIter
    )
    {
    
        if (refineParams.minRefineCells() == -1)
        {
            // Special setting to be able to restart shm on meshes with inconsistent
            // cellLevel/pointLevel
            return 0;
        }
    
    
        addProfiling(feature, "snappyHexMesh::refine::smallFeature");
    
        const fvMesh& mesh = meshRefiner_.mesh();
    
        label iter = 0;
    
        // See if any surface has an extendedGapLevel
    
        const labelList surfaceMaxLevel(meshRefiner_.surfaces().maxGapLevel());
        const labelList shellMaxLevel(meshRefiner_.shells().maxGapLevel());
        const labelList curvMaxLevel(meshRefiner_.surfaces().maxCurvatureLevel());
    
        if
        (
            max(surfaceMaxLevel) == 0
         && max(shellMaxLevel) == 0
         && max(curvMaxLevel) == 0
        )
    
        {
            return iter;
        }
    
        for (; iter < maxIter; iter++)
        {
            Info<< nl
                << "Small surface feature refinement iteration " << iter << nl
                << "--------------------------------------------" << nl
                << endl;
    
    
            // Determine cells to refine
            // ~~~~~~~~~~~~~~~~~~~~~~~~~
    
            labelList candidateCells
            (
                meshRefiner_.refineCandidates
                (
                    refineParams.locationsInMesh(),
                    refineParams.curvature(),
                    refineParams.planarAngle(),
    
                    false,              // featureRefinement
                    false,              // featureDistanceRefinement
                    false,              // internalRefinement
                    false,              // surfaceRefinement
                    false,              // curvatureRefinement
                    true,               // smallFeatureRefinement
                    false,              // gapRefinement
                    false,              // bigGapRefinement
    
                    false,              // spreadGapSize
    
                    refineParams.maxGlobalCells(),
                    refineParams.maxLocalCells()
                )
            );
    
            labelList cellsToRefine
            (
                meshRefiner_.meshCutter().consistentRefinement
                (
                    candidateCells,
                    true
                )
            );
            Info<< "Determined cells to refine in = "
                << mesh.time().cpuTimeIncrement() << " s" << endl;
    
    
    
            const label nCellsToRefine =
                returnReduce(cellsToRefine.size(), sumOp<label>());
    
    
            Info<< "Selected for refinement : " << nCellsToRefine
                << " cells (out of " << mesh.globalData().nTotalCells()
                << ')' << endl;
    
            // Stop when no cells to refine or have done minimum necessary
            // iterations and not enough cells to refine.
            if (nCellsToRefine == 0)
            {
                Info<< "Stopping refining since too few cells selected."
                    << nl << endl;
                break;
            }
    
    
            if (debug)
            {
                const_cast<Time&>(mesh.time())++;
            }
    
    
    
            if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
    
            {
                meshRefiner_.balanceAndRefine
                (
                    "small feature refinement iteration " + name(iter),
                    decomposer_,
                    distributor_,
                    cellsToRefine,
    
                    refineParams.maxLoadUnbalance(),
                    refineParams.maxCellUnbalance()
    
                );
            }
            else
            {
                meshRefiner_.refineAndBalance
                (
                    "small feature refinement iteration " + name(iter),
                    decomposer_,
                    distributor_,
                    cellsToRefine,
    
                    refineParams.maxLoadUnbalance(),
                    refineParams.maxCellUnbalance()
    
    Foam::label Foam::snappyRefineDriver::surfaceOnlyRefine
    
    (
        const refinementParameters& refineParams,
    
        const label maxIter,
        const label leakBlockageIter
    
        if (refineParams.minRefineCells() == -1)
        {
            // Special setting to be able to restart shm on meshes with inconsistent
            // cellLevel/pointLevel
            return 0;
        }
    
    
        addProfiling(surface, "snappyHexMesh::refine::surface");
    
        const fvMesh& mesh = meshRefiner_.mesh();
    
        const refinementSurfaces& surfaces = meshRefiner_.surfaces();
    
    
        // Determine the maximum refinement level over all surfaces. This
    
    Andrew Heather's avatar
    Andrew Heather committed
        // determines the minimum number of surface refinement iterations.
    
        label overallMaxLevel = max(meshRefiner_.surfaces().maxLevel());
    
        label iter;
        for (iter = 0; iter < maxIter; iter++)
        {
            Info<< nl
                << "Surface refinement iteration " << iter << nl
                << "------------------------------" << nl
                << endl;
    
    
    
            // Do optional leak closing (by removing cells)
            if (iter >= leakBlockageIter)
            {
                // Block off intersections with unzoned surfaces with specified
                // leakLevel < iter
                const labelList unnamedSurfaces
                (
                    surfaceZonesInfo::getUnnamedSurfaces
                    (
                        surfaces.surfZones()
                    )
                );
    
                DynamicList<label> selectedSurfaces(unnamedSurfaces.size());
                for (const label surfi : unnamedSurfaces)
                {
                    const label regioni = surfaces.globalRegion(surfi, 0);
    
                    // Take shortcut: assume all cells on surface are refined to
                    // its refinement level at iteration iter. So just use the
                    // iteration to see if the surface is a candidate.
                    if (iter > surfaces.leakLevel()[regioni])
                    {
                        selectedSurfaces.append(surfi);
                    }
                }
    
                if
                (
                    selectedSurfaces.size()
                 && refineParams.locationsOutsideMesh().size()
                )
                {
                   meshRefiner_.blockLeakFaces
                    (
                        globalToMasterPatch_,
                        globalToSlavePatch_,
                        refineParams.locationsInMesh(),
                        refineParams.zonesInMesh(),
                        refineParams.locationsOutsideMesh(),
                        selectedSurfaces
                    );
                }
            }
    
    
    
            // Determine cells to refine
            // ~~~~~~~~~~~~~~~~~~~~~~~~~
            // Only look at surface intersections (minLevel and surface curvature),
            // do not do internal refinement (refinementShells)
    
            labelList candidateCells
            (
                meshRefiner_.refineCandidates
                (
    
                    refineParams.curvature(),
    
                    refineParams.planarAngle(),
    
    
                    false,              // featureRefinement
    
                    false,              // featureDistanceRefinement
    
                    false,              // internalRefinement
                    true,               // surfaceRefinement
                    true,               // curvatureRefinement
    
                    false,              // smallFeatureRefinement
    
                    false,              // gapRefinement
    
                    false,              // spreadGapSize
    
                    refineParams.maxGlobalCells(),
                    refineParams.maxLocalCells()
                )
            );
            labelList cellsToRefine
            (
                meshRefiner_.meshCutter().consistentRefinement
                (
                    candidateCells,
                    true
                )
            );
            Info<< "Determined cells to refine in = "
                << mesh.time().cpuTimeIncrement() << " s" << endl;
    
    
    
            const label nCellsToRefine =
                returnReduce(cellsToRefine.size(), sumOp<label>());
    
    
            Info<< "Selected for refinement : " << nCellsToRefine
                << " cells (out of " << mesh.globalData().nTotalCells()
                << ')' << endl;
    
    
    andy's avatar
    andy committed
            // Stop when no cells to refine or have done minimum necessary
    
            // iterations and not enough cells to refine.
            if
            (
                nCellsToRefine == 0
             || (
                    iter >= overallMaxLevel
                 && nCellsToRefine <= refineParams.minRefineCells()
                )
            )
            {
                Info<< "Stopping refining since too few cells selected."
                    << nl << endl;
                break;
            }
    
    
            if (debug)
            {
                const_cast<Time&>(mesh.time())++;
            }
    
    
            if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
    
            {
                meshRefiner_.balanceAndRefine
                (
                    "surface refinement iteration " + name(iter),
                    decomposer_,
                    distributor_,
                    cellsToRefine,
    
                    refineParams.maxLoadUnbalance(),
                    refineParams.maxCellUnbalance()
    
                );
            }
            else
            {
                meshRefiner_.refineAndBalance
                (
                    "surface refinement iteration " + name(iter),
                    decomposer_,
                    distributor_,
                    cellsToRefine,
    
                    refineParams.maxLoadUnbalance(),
                    refineParams.maxCellUnbalance()
    
    Foam::label Foam::snappyRefineDriver::gapOnlyRefine
    
    (
        const refinementParameters& refineParams,
        const label maxIter
    )
    {
    
        if (refineParams.minRefineCells() == -1)
        {
            // Special setting to be able to restart shm on meshes with inconsistent
            // cellLevel/pointLevel
            return 0;
        }
    
    
        const fvMesh& mesh = meshRefiner_.mesh();
    
        // Determine the maximum refinement level over all surfaces. This
    
    Andrew Heather's avatar
    Andrew Heather committed
        // determines the minimum number of surface refinement iterations.
    
    
        label maxIncrement = 0;
        const labelList& maxLevel = meshRefiner_.surfaces().maxLevel();
        const labelList& gapLevel = meshRefiner_.surfaces().gapLevel();
    
        forAll(maxLevel, i)
        {
            maxIncrement = max(maxIncrement, gapLevel[i]-maxLevel[i]);
        }
    
        label iter = 0;
    
        if (maxIncrement == 0)
        {
            return iter;
        }
    
        for (iter = 0; iter < maxIter; iter++)
        {
            Info<< nl
                << "Gap refinement iteration " << iter << nl
                << "--------------------------" << nl
                << endl;
    
    
            // Determine cells to refine
            // ~~~~~~~~~~~~~~~~~~~~~~~~~
            // Only look at surface intersections (minLevel and surface curvature),
            // do not do internal refinement (refinementShells)
    
            labelList candidateCells
            (
                meshRefiner_.refineCandidates
                (
    
                    refineParams.curvature(),
                    refineParams.planarAngle(),
    
                    false,              // featureRefinement
                    false,              // featureDistanceRefinement
                    false,              // internalRefinement
                    false,              // surfaceRefinement
                    false,              // curvatureRefinement
    
                    false,              // smallFeatureRefinement
    
                    true,               // gapRefinement
    
                    false,              // spreadGapSize
    
                    refineParams.maxGlobalCells(),
                    refineParams.maxLocalCells()
                )
            );
    
            if (debug&meshRefinement::MESH)
            {
    
                Pout<< "Writing current mesh to time "
                    << meshRefiner_.timeName() << endl;
                meshRefiner_.write
                (
                    meshRefinement::debugType(debug),
                    meshRefinement::writeType
                    (
                        meshRefinement::writeLevel()
                      | meshRefinement::WRITEMESH
                    ),
                    mesh.time().path()/meshRefiner_.timeName()
                );
                Pout<< "Dumped mesh in = "
                    << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
    
    
    
                Pout<< "Dumping " << candidateCells.size()
                    << " cells to cellSet candidateCellsFromGap." << endl;
                cellSet c(mesh, "candidateCellsFromGap", candidateCells);
                c.instance() = meshRefiner_.timeName();
                c.write();
            }
    
            // Grow by one layer to make sure we're covering the gap
            {
                boolList isCandidateCell(mesh.nCells(), false);
                forAll(candidateCells, i)
                {
                    isCandidateCell[candidateCells[i]] = true;
                }
    
                for (label i=0; i<1; i++)
                {
                    boolList newIsCandidateCell(isCandidateCell);
    
                    // Internal faces
    
                    for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
    
                        label own = mesh.faceOwner()[facei];
                        label nei = mesh.faceNeighbour()[facei];
    
    
                        if (isCandidateCell[own] != isCandidateCell[nei])
                        {
                            newIsCandidateCell[own] = true;
                            newIsCandidateCell[nei] = true;
                        }
                    }
    
                    // Get coupled boundary condition values
                    boolList neiIsCandidateCell;
                    syncTools::swapBoundaryCellList
                    (
                        mesh,
                        isCandidateCell,
                        neiIsCandidateCell
                    );
    
                    // Boundary faces
                    for
                    (
    
                        label facei = mesh.nInternalFaces();
                        facei < mesh.nFaces();
                        facei++
    
                        label own = mesh.faceOwner()[facei];
    
                        label bFacei = facei-mesh.nInternalFaces();
    
                        if (isCandidateCell[own] != neiIsCandidateCell[bFacei])
    
                        {
                            newIsCandidateCell[own] = true;
                        }
                    }
    
                    isCandidateCell.transfer(newIsCandidateCell);
                }
    
                label n = 0;
    
                forAll(isCandidateCell, celli)
    
                    if (isCandidateCell[celli])
    
                    {
                        n++;
                    }
                }
                candidateCells.setSize(n);
                n = 0;
    
                forAll(isCandidateCell, celli)
    
                    if (isCandidateCell[celli])
    
                        candidateCells[n++] = celli;
    
                    }
                }
            }
    
    
            if (debug&meshRefinement::MESH)
            {
                Pout<< "Dumping " << candidateCells.size()
                    << " cells to cellSet candidateCellsFromGapPlusBuffer." << endl;
                cellSet c(mesh, "candidateCellsFromGapPlusBuffer", candidateCells);
                c.instance() = meshRefiner_.timeName();
                c.write();
            }
    
    
            labelList cellsToRefine
            (
                meshRefiner_.meshCutter().consistentRefinement
                (
                    candidateCells,
                    true
                )
            );
            Info<< "Determined cells to refine in = "
                << mesh.time().cpuTimeIncrement() << " s" << endl;
    
    
    
            const label nCellsToRefine =
                returnReduce(cellsToRefine.size(), sumOp<label>());
    
    
            Info<< "Selected for refinement : " << nCellsToRefine
                << " cells (out of " << mesh.globalData().nTotalCells()
                << ')' << endl;
    
            // Stop when no cells to refine or have done minimum necessary
            // iterations and not enough cells to refine.
            if
            (
                nCellsToRefine == 0
             || (
                    iter >= maxIncrement
                 && nCellsToRefine <= refineParams.minRefineCells()
                )
            )
            {
                Info<< "Stopping refining since too few cells selected."
                    << nl << endl;
                break;
            }
    
    
            if (debug)
            {
                const_cast<Time&>(mesh.time())++;
            }
    
    
    
            if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
    
            {
                meshRefiner_.balanceAndRefine
                (
                    "gap refinement iteration " + name(iter),
                    decomposer_,
                    distributor_,
                    cellsToRefine,
    
                    refineParams.maxLoadUnbalance(),
                    refineParams.maxCellUnbalance()
    
                );
            }
            else
            {
                meshRefiner_.refineAndBalance
                (
                    "gap refinement iteration " + name(iter),
                    decomposer_,
                    distributor_,
                    cellsToRefine,
    
                    refineParams.maxLoadUnbalance(),
                    refineParams.maxCellUnbalance()
    
    Foam::label Foam::snappyRefineDriver::surfaceProximityBlock
    (
        const refinementParameters& refineParams,
        const label maxIter
    )
    {
        if (refineParams.minRefineCells() == -1)
        {
            // Special setting to be able to restart shm on meshes with inconsistent
            // cellLevel/pointLevel
            return 0;
        }
    
        fvMesh& mesh = meshRefiner_.mesh();
    
        if (min(meshRefiner_.surfaces().blockLevel()) == labelMax)
        {
            return 0;
        }
    
        label iter = 0;
    
        for (iter = 0; iter < maxIter; iter++)
        {
            Info<< nl
                << "Gap blocking iteration " << iter << nl
                << "------------------------" << nl
                << endl;
    
    
            // Determine cells to block
            // ~~~~~~~~~~~~~~~~~~~~~~~~
    
            meshRefiner_.removeGapCells
            (
                refineParams.planarAngle(),
                meshRefiner_.surfaces().blockLevel(),
                globalToMasterPatch_,
                refineParams.nFilterIter()
            );
    
            if (debug)
            {
                const_cast<Time&>(mesh.time())++;
    
                Pout<< "Writing gap blocking iteration "
                    << iter << " mesh to time " << meshRefiner_.timeName()
                    << endl;
                meshRefiner_.write
                (
                    meshRefinement::debugType(debug),
                    meshRefinement::writeType
                    (
                        meshRefinement::writeLevel()
                      | meshRefinement::WRITEMESH
                    ),
                    mesh.time().path()/meshRefiner_.timeName()
                );
    
    Foam::label Foam::snappyRefineDriver::bigGapOnlyRefine
    
    (
        const refinementParameters& refineParams,
    
        const bool spreadGapSize,
    
        if (refineParams.minRefineCells() == -1)
        {
            // Special setting to be able to restart shm on meshes with inconsistent
            // cellLevel/pointLevel
            return 0;
        }
    
    
        const fvMesh& mesh = meshRefiner_.mesh();
    
        label iter = 0;
    
        // See if any surface has an extendedGapLevel
        labelList surfaceMaxLevel(meshRefiner_.surfaces().maxGapLevel());
        labelList shellMaxLevel(meshRefiner_.shells().maxGapLevel());
    
        label overallMaxLevel(max(max(surfaceMaxLevel), max(shellMaxLevel)));
    
        if (overallMaxLevel == 0)
        {
            return iter;
        }
    
    
        for (; iter < maxIter; iter++)
        {
            Info<< nl
                << "Big gap refinement iteration " << iter << nl
                << "------------------------------" << nl
                << endl;
    
    
            // Determine cells to refine
            // ~~~~~~~~~~~~~~~~~~~~~~~~~
    
            labelList candidateCells
            (
                meshRefiner_.refineCandidates
                (
                    refineParams.locationsInMesh(),
                    refineParams.curvature(),
                    refineParams.planarAngle(),
    
                    false,              // featureRefinement
                    false,              // featureDistanceRefinement
                    false,              // internalRefinement
                    false,              // surfaceRefinement
                    false,              // curvatureRefinement
                    false,              // smallFeatureRefinement
                    false,              // gapRefinement
                    true,               // bigGapRefinement
    
                    spreadGapSize,      // spreadGapSize
    
                    refineParams.maxGlobalCells(),
                    refineParams.maxLocalCells()
                )
            );
    
    
            if (debug&meshRefinement::MESH)
            {
    
                Pout<< "Writing current mesh to time "
                    << meshRefiner_.timeName() << endl;
                meshRefiner_.write
                (
                    meshRefinement::debugType(debug),
                    meshRefinement::writeType
                    (
                        meshRefinement::writeLevel()
                      | meshRefinement::WRITEMESH
                    ),
                    mesh.time().path()/meshRefiner_.timeName()
                );
                Pout<< "Dumped mesh in = "
                    << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
    
    
                Pout<< "Dumping " << candidateCells.size()
                    << " cells to cellSet candidateCellsFromBigGap." << endl;
                cellSet c(mesh, "candidateCellsFromBigGap", candidateCells);
                c.instance() = meshRefiner_.timeName();
                c.write();
            }
    
            labelList cellsToRefine
            (
                meshRefiner_.meshCutter().consistentRefinement
                (
                    candidateCells,
                    true
                )
            );
            Info<< "Determined cells to refine in = "
                << mesh.time().cpuTimeIncrement() << " s" << endl;
    
    
    
            const label nCellsToRefine =
                returnReduce(cellsToRefine.size(), sumOp<label>());
    
    
            Info<< "Selected for refinement : " << nCellsToRefine
                << " cells (out of " << mesh.globalData().nTotalCells()
                << ')' << endl;
    
            // Stop when no cells to refine or have done minimum necessary
            // iterations and not enough cells to refine.
            if
            (
                nCellsToRefine == 0
             || (
                    iter >= overallMaxLevel
                 && nCellsToRefine <= refineParams.minRefineCells()
                )
            )
            {
                Info<< "Stopping refining since too few cells selected."
                    << nl << endl;
                break;
            }
    
    
            if (debug)
            {
                const_cast<Time&>(mesh.time())++;
            }
    
    
    
            if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
    
            {
                meshRefiner_.balanceAndRefine
                (
                    "big gap refinement iteration " + name(iter),
                    decomposer_,
                    distributor_,
                    cellsToRefine,
    
                    refineParams.maxLoadUnbalance(),
                    refineParams.maxCellUnbalance()
    
                );
            }
            else
            {
                meshRefiner_.refineAndBalance
                (
                    "big gap refinement iteration " + name(iter),
                    decomposer_,
                    distributor_,
                    cellsToRefine,
    
                    refineParams.maxLoadUnbalance(),
                    refineParams.maxCellUnbalance()