meshRefinementRefine.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2026 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.
 
    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
\*---------------------------------------------------------------------------*/
 
#include "meshRefinement.H"
#include "trackedParticle.H"
#include "syncTools.H"
#include "Time.H"
#include "refinementSurfaces.H"
#include "refinementFeatures.H"
#include "refinementRegions.H"
#include "faceSet.H"
#include "decompositionMethod.H"
#include "fvMeshDistribute.H"
#include "polyTopoChange.H"
#include "polyDistributionMap.H"
#include "Cloud.H"
#include "OBJstream.H"
#include "cellSet.H"
#include "meshSearch.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
    //- To compare normals
    static bool less(const vector& x, const vector& y)
    {
        for (direction i = 0; i < vector::nComponents; i++)
        {
            if (x[i] < y[i])
            {
                return true;
            }
            else if (x[i] > y[i])
            {
                return false;
            }
        }
        // All components the same
        return false;
    }
 
 
    //- To compare normals
    class normalLess
    {
        const vectorList& values_;
 
    public:
 
        normalLess(const vectorList& values)
        :
            values_(values)
        {}
 
        bool operator()(const label a, const label b) const
        {
            return less(values_[a], values_[b]);
        }
    };
 
 
    //- Template specialisation for pTraits<labelList> so we can have fields
    template<>
    class pTraits<labelList>
    {
 
    public:
 
        //- Component type
        typedef labelList cmptType;
    };
 
    //- Template specialisation for pTraits<labelList> so we can have fields
    template<>
    class pTraits<vectorList>
    {
 
    public:
 
        //- Component type
        typedef vectorList cmptType;
    };
}
 
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
// Get faces (on the new mesh) that have in some way been affected by the
// mesh change. Picks up all faces but those that are between two
// unrefined faces. (Note that of an unchanged face the edge still might be
// split but that does not change any face centre or cell centre.
Foam::labelList Foam::meshRefinement::getChangedFaces
(
    const polyTopoChangeMap& map,
    const labelList& oldCellsToRefine
)
{
    const polyMesh& mesh = map.mesh();
 
    labelList changedFaces;
 
    // For reporting: number of masterFaces changed
    label nMasterChanged = 0;
    PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
 
    {
        // Mark any face on a cell which has been added or changed
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        // Note that refining a face changes the face centre (for a warped face)
        // which changes the cell centre. This again changes the cellcentre-
        // cellcentre edge across all faces of the cell.
        // Note: this does not happen for unwarped faces but unfortunately
        // we don't have this information.
 
        const labelList& faceOwner = mesh.faceOwner();
        const labelList& faceNeighbour = mesh.faceNeighbour();
        const cellList& cells = mesh.cells();
        const label nInternalFaces = mesh.nInternalFaces();
 
        // Mark refined cells on old mesh
        PackedBoolList oldRefineCell(map.nOldCells());
 
        forAll(oldCellsToRefine, i)
        {
            oldRefineCell.set(oldCellsToRefine[i], 1u);
        }
 
        // Mark refined faces
        PackedBoolList refinedInternalFace(nInternalFaces);
 
        // 1. Internal faces
 
        for (label facei = 0; facei < nInternalFaces; facei++)
        {
            const label oldOwn = map.cellMap()[faceOwner[facei]];
            const label oldNei = map.cellMap()[faceNeighbour[facei]];
 
            if
            (
                oldOwn >= 0
             && oldRefineCell.get(oldOwn) == 0u
             && oldNei >= 0
             && oldRefineCell.get(oldNei) == 0u
            )
            {
                // Unaffected face since both neighbours were not refined.
            }
            else
            {
                refinedInternalFace.set(facei, 1u);
            }
        }
 
 
        // 2. Boundary faces
 
        boolList refinedBoundaryFace(mesh.nFaces()-nInternalFaces, false);
 
        forAll(mesh.boundary(), patchi)
        {
            const polyPatch& pp = mesh.boundary()[patchi];
 
            label facei = pp.start();
 
            forAll(pp, i)
            {
                label oldOwn = map.cellMap()[faceOwner[facei]];
 
                if (oldOwn >= 0 && oldRefineCell.get(oldOwn) == 0u)
                {
                    // owner did exist and wasn't refined.
                }
                else
                {
                    refinedBoundaryFace[facei - nInternalFaces] = true;
                }
                facei++;
            }
        }
 
        // Synchronise refined face status
        syncTools::syncBoundaryFaceList
        (
            mesh,
            refinedBoundaryFace,
            orEqOp<bool>()
        );
 
 
        // 3. Mark all faces affected by refinement. Refined faces are in
        //    - refinedInternalFace
        //    - refinedBoundaryFace
        boolList changedFace(mesh.nFaces(), false);
 
        forAll(refinedInternalFace, facei)
        {
            if (refinedInternalFace.get(facei) == 1u)
            {
                const cell& ownFaces = cells[faceOwner[facei]];
                forAll(ownFaces, owni)
                {
                    changedFace[ownFaces[owni]] = true;
                }
                const cell& neiFaces = cells[faceNeighbour[facei]];
                forAll(neiFaces, neii)
                {
                    changedFace[neiFaces[neii]] = true;
                }
            }
        }
 
        forAll(refinedBoundaryFace, i)
        {
            if (refinedBoundaryFace[i])
            {
                const cell& ownFaces = cells[faceOwner[i+nInternalFaces]];
                forAll(ownFaces, owni)
                {
                    changedFace[ownFaces[owni]] = true;
                }
            }
        }
 
        syncTools::syncFaceList
        (
            mesh,
            changedFace,
            orEqOp<bool>()
        );
 
 
        // Now we have in changedFace marked all affected faces. Pack.
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
        changedFaces = findIndices(changedFace, true);
 
        // Count changed master faces.
        nMasterChanged = 0;
 
        forAll(changedFace, facei)
        {
            if (changedFace[facei] && isMasterFace[facei])
            {
                nMasterChanged++;
            }
        }
 
    }
 
    if (debug&meshRefinement::MESH)
    {
        Pout<< "getChangedFaces : Detected "
            << " local:" << changedFaces.size()
            << " global:" << returnReduce(nMasterChanged, sumOp<label>())
            << " changed faces out of " << mesh.globalData().nTotalFaces()
            << endl;
 
        faceSet changedFacesSet(mesh, "changedFaces", changedFaces);
        Pout<< "getChangedFaces : Writing " << changedFaces.size()
            << " changed faces to faceSet " << changedFacesSet.name()
            << endl;
        changedFacesSet.write();
    }
 
    return changedFaces;
}
 
 
// Mark cell for refinement (if not already marked). Return false if
// refinelimit hit. Keeps running count (in nRefine) of cells marked for
// refinement
bool Foam::meshRefinement::markForRefine
(
    const label markValue,
    const label nAllowRefine,
 
    label& cellValue,
    label& nRefine
)
{
    if (cellValue == -1)
    {
        cellValue = markValue;
        nRefine++;
    }
 
    return nRefine <= nAllowRefine;
}
 
 
void Foam::meshRefinement::markFeatureCellLevel
(
    const List<point>& insidePoints,
    labelList& maxFeatureLevel
) const
{
    // We want to refine all cells containing a feature edge.
    // - don't want to search over whole mesh
    // - don't want to build octree for whole mesh
    // - so use tracking to follow the feature edges
    //
    // 1. find non-manifold points on feature edges (i.e. start of feature edge
    //    or 'knot')
    // 2. seed particle starting at insidePoint going to this non-manifold
    //    point
    // 3. track particles to their non-manifold point
    //
    // 4. track particles across their connected feature edges, marking all
    //    visited cells with their level (through trackingData)
    // 5. do 4 until all edges have been visited.
 
    // Find all start cells of features.
    // Is done by tracking from insidePoint.
    lagrangian::Cloud<trackedParticle> startPointCloud
    (
        mesh_,
        "startPointCloud",
        IDLList<trackedParticle>()
    );
 
    // Features are identical on all processors. Number them so we know
    // what to seed. Do this on only the processor that
    // holds the insidePoint.
    label nLocateBoundaryHits = 0;
 
    const meshSearch& searchEngine = meshSearch::New(mesh());
 
    forAll(insidePoints, i)
    {
        const point& insidePoint = insidePoints[i];
 
        const label celli = searchEngine.findCell(insidePoint);
 
        if (celli != -1)
        {
            // I am the processor that holds the insidePoint
 
            forAll(features_, feati)
            {
                const edgeMesh& featureMesh = features_[feati];
                const label featureLevel = features_.levels()[feati][0];
                const labelListList& pointEdges = featureMesh.pointEdges();
 
                // Find regions on edgeMesh
                labelList edgeRegion;
                const label nRegions = featureMesh.regions(edgeRegion);
 
 
                PackedBoolList regionVisited(nRegions);
 
 
                // 1. Seed all 'knots' in edgeMesh
 
 
                forAll(pointEdges, pointi)
                {
                    if (pointEdges[pointi].size() != 2)
                    {
                        if (debug&meshRefinement::FEATURESEEDS)
                        {
                            Pout<< "Adding particle from point:" << pointi
                                << " coord:" << featureMesh.points()[pointi]
                                << " since number of emanating edges:"
                                << pointEdges[pointi].size()
                                << endl;
                        }
 
                        // Non-manifold point. Create particle.
                        startPointCloud.addParticle
                        (
                            new trackedParticle
                            (
                                searchEngine,
                                insidePoint,
                                celli,
                                nLocateBoundaryHits,
                                featureMesh.points()[pointi],   // endpos
                                featureLevel,                   // level
                                feati,                          // featureMesh
                                pointi,                         // end point
                                -1                              // feature edge
                            )
                        );
 
                        // Mark
                        if (pointEdges[pointi].size() > 0)
                        {
                            label e0 = pointEdges[pointi][0];
                            label regioni = edgeRegion[e0];
                            regionVisited[regioni] = 1u;
                        }
                    }
                }
 
 
                // 2. Any regions that have not been visited at all? These can
                //    only be circular regions!
                forAll(featureMesh.edges(), edgei)
                {
                    if (regionVisited.set(edgeRegion[edgei], 1u))
                    {
                        const edge& e = featureMesh.edges()[edgei];
                        const label pointi = e.start();
                        if (debug&meshRefinement::FEATURESEEDS)
                        {
                            Pout<< "Adding particle from point:" << pointi
                                << " coord:" << featureMesh.points()[pointi]
                                << " on circular region:" << edgeRegion[edgei]
                                << endl;
                        }
 
                        // Non-manifold point. Create particle.
                        startPointCloud.addParticle
                        (
                            new trackedParticle
                            (
                                searchEngine,
                                insidePoint,
                                celli,
                                nLocateBoundaryHits,
                                featureMesh.points()[pointi],   // endpos
                                featureLevel,                   // level
                                feati,                          // featureMesh
                                pointi,                         // end point
                                -1                              // feature edge
                            )
                        );
                    }
                }
            }
        }
    }
 
 
    // Largest refinement level of any feature passed through
    maxFeatureLevel = labelList(mesh_.nCells(), -1);
 
    // Whether edge has been visited.
    List<PackedBoolList> featureEdgeVisited(features_.size());
 
    forAll(features_, feati)
    {
        featureEdgeVisited[feati].setSize(features_[feati].edges().size());
        featureEdgeVisited[feati] = 0u;
    }
 
    // Database to pass into trackedParticle::move
    trackedParticle::trackingData td
    (
        startPointCloud,
        2.0*mesh_.bounds().mag(),
        maxFeatureLevel,
        featureEdgeVisited
    );
 
 
    // Track all particles to their end position (= starting feature point)
    // Note that the particle might have started on a different processor
    // so this will transfer across nicely until we can start tracking proper.
    if (debug&meshRefinement::FEATURESEEDS)
    {
        Pout<< "Tracking " << startPointCloud.size()
            << " particles over distance " << td.maxTrackLen_
            << " to find the starting cell" << endl;
    }
    startPointCloud.move(startPointCloud, td);
 
 
    // Reset levels
    maxFeatureLevel = -1;
    forAll(features_, feati)
    {
        featureEdgeVisited[feati] = 0u;
    }
 
 
    lagrangian::Cloud<trackedParticle> cloud
    (
        mesh_,
        "featureCloud",
        IDLList<trackedParticle>()
    );
 
    if (debug&meshRefinement::FEATURESEEDS)
    {
        Pout<< "Constructing cloud for cell marking" << endl;
    }
 
    forAllIter(lagrangian::Cloud<trackedParticle>, startPointCloud, iter)
    {
        const trackedParticle& startTp = iter();
 
        const label feati = startTp.i();
        const label pointi = startTp.j();
 
        const edgeMesh& featureMesh = features_[feati];
        const labelList& pEdges = featureMesh.pointEdges()[pointi];
 
        // Now shoot particles down all pEdges.
        forAll(pEdges, pEdgei)
        {
            const label edgei = pEdges[pEdgei];
 
            if (featureEdgeVisited[feati].set(edgei, 1u))
            {
                // Unvisited edge. Make the particle go to the other point
                // on the edge.
 
                const edge& e = featureMesh.edges()[edgei];
                label otherPointi = e.otherVertex(pointi);
 
                trackedParticle* tp(new trackedParticle(startTp));
                tp->start() = tp->position(mesh_);
                tp->end() = featureMesh.points()[otherPointi];
                tp->fraction() = 1;
                tp->j() = otherPointi;
                tp->k() = edgei;
 
                if (debug&meshRefinement::FEATURESEEDS)
                {
                    Pout<< "Adding particle for point:" << pointi
                        << " coord:" << tp->position(mesh_)
                        << " feature:" << feati
                        << " to track to:" << tp->end()
                        << endl;
                }
 
                cloud.addParticle(tp);
            }
        }
    }
 
    startPointCloud.clear();
 
 
    while (true)
    {
        // Track all particles to their end position.
        if (debug&meshRefinement::FEATURESEEDS)
        {
            Pout<< "Tracking " << cloud.size()
                << " particles over distance " << td.maxTrackLen_
                << " to mark cells" << endl;
        }
        cloud.move(cloud, td);
 
        // Make particle follow edge.
        forAllIter(lagrangian::Cloud<trackedParticle>, cloud, iter)
        {
            trackedParticle& tp = iter();
 
            const label feati = tp.i();
            const label pointi = tp.j();
 
            const edgeMesh& featureMesh = features_[feati];
            const labelList& pEdges = featureMesh.pointEdges()[pointi];
 
            // Particle now at pointi. Check connected edges to see which one
            // we have to visit now.
 
            bool keepParticle = false;
 
            forAll(pEdges, i)
            {
                const label edgei = pEdges[i];
 
                if (featureEdgeVisited[feati].set(edgei, 1u))
                {
                    // Unvisited edge. Make the particle go to the other point
                    // on the edge.
 
                    const edge& e = featureMesh.edges()[edgei];
                    const label otherPointi = e.otherVertex(pointi);
 
                    tp.start() = tp.position(mesh_);
                    tp.end() = featureMesh.points()[otherPointi];
                    tp.fraction() = 1;
                    tp.j() = otherPointi;
                    tp.k() = edgei;
                    keepParticle = true;
                    break;
                }
            }
 
            if (!keepParticle)
            {
                // Particle at 'knot' where another particle already has been
                // seeded. Delete particle.
                cloud.deleteParticle(tp);
            }
        }
 
 
        if (debug&meshRefinement::FEATURESEEDS)
        {
            Pout<< "Remaining particles " << cloud.size() << endl;
        }
 
        if (returnReduce(cloud.size(), sumOp<label>()) == 0)
        {
            break;
        }
    }
 
 
 
    // if (debug&meshRefinement::FEATURESEEDS)
    //{
    //    forAll(maxFeatureLevel, celli)
    //    {
    //        if (maxFeatureLevel[celli] != -1)
    //        {
    //            Pout<< "Feature went through cell:" << celli
    //                << " coord:" << mesh_.cellCentres()[celli]
    //                << " level:" << maxFeatureLevel[celli]
    //                << endl;
    //        }
    //    }
    //}
}
 
 
// Calculates list of cells to refine based on intersection with feature edge.
Foam::label Foam::meshRefinement::markFeatureRefinement
(
    const List<point>& insidePoints,
    const label nAllowRefine,
 
    labelList& refineCell,
    label& nRefine
) const
{
    // Largest refinement level of any feature passed through
    labelList maxFeatureLevel;
    markFeatureCellLevel(insidePoints, maxFeatureLevel);
 
    // See which cells to refine. maxFeatureLevel will hold highest level
    // of any feature edge that passed through.
 
    const labelList& cellLevel = meshCutter_.cellLevel();
 
    const label oldNRefine = nRefine;
 
    forAll(maxFeatureLevel, celli)
    {
        if (maxFeatureLevel[celli] > cellLevel[celli])
        {
            // Mark
            if
            (
               !markForRefine
                (
                    0,                      // surface (n/a)
                    nAllowRefine,
                    refineCell[celli],
                    nRefine
                )
            )
            {
                // Reached limit
                break;
            }
        }
    }
 
    if
    (
        returnReduce(nRefine, sumOp<label>())
      > returnReduce(nAllowRefine, sumOp<label>())
    )
    {
        Info<< "Reached refinement limit." << endl;
    }
 
    return returnReduce(nRefine-oldNRefine,  sumOp<label>());
}
 
 
// Mark cells for distance-to-feature based refinement.
Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
(
    const label nAllowRefine,
 
    labelList& refineCell,
    label& nRefine
) const
{
    const labelList& cellLevel = meshCutter_.cellLevel();
    const pointField& cellCentres = mesh_.cellCentres();
 
    // Detect if there are any distance features
    if (features_.maxDistance() <= 0.0)
    {
        return 0;
    }
 
    label oldNRefine = nRefine;
 
    // Collect cells to test
    pointField testCc(cellLevel.size()-nRefine);
    labelList testLevels(cellLevel.size()-nRefine);
    label testi = 0;
 
    forAll(cellLevel, celli)
    {
        if (refineCell[celli] == -1)
        {
            testCc[testi] = cellCentres[celli];
            testLevels[testi] = cellLevel[celli];
            testi++;
        }
    }
 
    // Do test to see whether cells are in/near features with higher level
    labelList maxLevel;
    features_.findHigherLevel(testCc, testLevels, maxLevel);
 
    // Mark for refinement. Note that we didn't store the original cellID so
    // now just reloop in same order.
    testi = 0;
    forAll(cellLevel, celli)
    {
        if (refineCell[celli] == -1)
        {
            if (maxLevel[testi] > testLevels[testi])
            {
                const bool reachedLimit = !markForRefine
                (
                    maxLevel[testi],    // mark with any positive value
                    nAllowRefine,
                    refineCell[celli],
                    nRefine
                );
 
                if (reachedLimit)
                {
                    if (debug)
                    {
                        Pout<< "Stopped refining internal cells"
                            << " since reaching my cell limit of "
                            << mesh_.nCells()+7*nRefine << endl;
                    }
                    break;
                }
            }
            testi++;
        }
    }
 
    if
    (
        returnReduce(nRefine, sumOp<label>())
      > returnReduce(nAllowRefine, sumOp<label>())
    )
    {
        Info<< "Reached refinement limit." << endl;
    }
 
    return returnReduce(nRefine-oldNRefine, sumOp<label>());
}
 
 
Foam::label Foam::meshRefinement::markInternalRefinement
(
    const label nAllowRefine,
 
    labelList& refineCell,
    label& nRefine
) const
{
    const labelList& cellLevel = meshCutter_.cellLevel();
    const pointField& cellCentres = mesh_.cellCentres();
 
    const label oldNRefine = nRefine;
 
    // Collect cells to test
    pointField testCc(cellLevel.size()-nRefine);
    labelList testLevels(cellLevel.size()-nRefine);
    label testi = 0;
 
    forAll(cellLevel, celli)
    {
        if (refineCell[celli] == -1)
        {
            testCc[testi] = cellCentres[celli];
            testLevels[testi] = cellLevel[celli];
            testi++;
        }
    }
 
    // Do test to see whether cells is inside/outside shell with higher level
    labelList maxLevel;
    shells_.findHigherLevel
    (
        testCc,
        testLevels,
        meshCutter_.level0EdgeLength(),
        maxLevel
    );
 
    // Mark for refinement. Note that we didn't store the original cellID so
    // now just reloop in same order.
    testi = 0;
    forAll(cellLevel, celli)
    {
        if (refineCell[celli] == -1)
        {
            if (maxLevel[testi] > testLevels[testi])
            {
                bool reachedLimit = !markForRefine
                (
                    maxLevel[testi],    // mark with any positive value
                    nAllowRefine,
                    refineCell[celli],
                    nRefine
                );
 
                if (reachedLimit)
                {
                    if (debug)
                    {
                        Pout<< "Stopped refining internal cells"
                            << " since reaching my cell limit of "
                            << mesh_.nCells()+7*nRefine << endl;
                    }
                    break;
                }
            }
            testi++;
        }
    }
 
    if
    (
        returnReduce(nRefine, sumOp<label>())
      > returnReduce(nAllowRefine, sumOp<label>())
    )
    {
        Info<< "Reached refinement limit." << endl;
    }
 
    return returnReduce(nRefine-oldNRefine, sumOp<label>());
}
 
 
Foam::labelList Foam::meshRefinement::getRefineCandidateFaces
(
    const labelList& refineCell
) const
{
    labelList testFaces(mesh_.nFaces());
 
    label nTest = 0;
 
    forAll(surfaceIndex_, facei)
    {
        if (surfaceIndex_[facei] != -1)
        {
            label own = mesh_.faceOwner()[facei];
 
            if (mesh_.isInternalFace(facei))
            {
                const label nei = mesh_.faceNeighbour()[facei];
 
                if (refineCell[own] == -1 || refineCell[nei] == -1)
                {
                    testFaces[nTest++] = facei;
                }
            }
            else
            {
                if (refineCell[own] == -1)
                {
                    testFaces[nTest++] = facei;
                }
            }
        }
    }
    testFaces.setSize(nTest);
 
    return testFaces;
}
 
 
Foam::label Foam::meshRefinement::markSurfaceRefinement
(
    const label nAllowRefine,
    const labelList& neiLevel,
    const pointField& neiCc,
 
    labelList& refineCell,
    label& nRefine
) const
{
    const labelList& cellLevel = meshCutter_.cellLevel();
    const pointField& cellCentres = mesh_.cellCentres();
 
    const label oldNRefine = nRefine;
 
    // Use cached surfaceIndex_ to detect if any intersection. If so
    // re-intersect to determine level wanted.
 
    // Collect candidate faces
    // ~~~~~~~~~~~~~~~~~~~~~~~
 
    labelList testFaces(getRefineCandidateFaces(refineCell));
 
    // Collect segments
    // ~~~~~~~~~~~~~~~~
 
    pointField start(testFaces.size());
    pointField end(testFaces.size());
    labelList minLevel(testFaces.size());
 
    forAll(testFaces, i)
    {
        const label facei = testFaces[i];
 
        const label own = mesh_.faceOwner()[facei];
 
        if (mesh_.isInternalFace(facei))
        {
            label nei = mesh_.faceNeighbour()[facei];
 
            start[i] = cellCentres[own];
            end[i] = cellCentres[nei];
            minLevel[i] = min(cellLevel[own], cellLevel[nei]);
        }
        else
        {
            const label bFacei = facei - mesh_.nInternalFaces();
 
            start[i] = cellCentres[own];
            end[i] = neiCc[bFacei];
            minLevel[i] = min(cellLevel[own], neiLevel[bFacei]);
        }
    }
 
    // Extend segments a bit
    {
        const vectorField smallVec(rootSmall*(end-start));
        start -= smallVec;
        end += smallVec;
    }
 
 
    // Do test for higher intersections
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    labelList surfaceHit;
    labelList surfaceMinLevel;
    surfaces_.findHigherIntersection
    (
        start,
        end,
        minLevel,
 
        surfaceHit,
        surfaceMinLevel
    );
 
 
    // Mark cells for refinement
    // ~~~~~~~~~~~~~~~~~~~~~~~~~
 
    forAll(testFaces, i)
    {
        const label facei = testFaces[i];
 
        const label surfi = surfaceHit[i];
 
        if (surfi != -1)
        {
            // Found intersection with surface with higher wanted
            // refinement. Check if the level field on that surface
            // specifies an even higher level. Note:this is weird. Should
            // do the check with the surfaceMinLevel whilst intersecting the
            // surfaces?
 
            const label own = mesh_.faceOwner()[facei];
 
            if (surfaceMinLevel[i] > cellLevel[own])
            {
                // Owner needs refining
                if
                (
                   !markForRefine
                    (
                        surfi,
                        nAllowRefine,
                        refineCell[own],
                        nRefine
                    )
                )
                {
                    break;
                }
            }
 
            if (mesh_.isInternalFace(facei))
            {
                const label nei = mesh_.faceNeighbour()[facei];
 
                if (surfaceMinLevel[i] > cellLevel[nei])
                {
                    // Neighbour needs refining
                    if
                    (
                       !markForRefine
                        (
                            surfi,
                            nAllowRefine,
                            refineCell[nei],
                            nRefine
                        )
                    )
                    {
                        break;
                    }
                }
            }
        }
    }
 
    if
    (
        returnReduce(nRefine, sumOp<label>())
      > returnReduce(nAllowRefine, sumOp<label>())
    )
    {
        Info<< "Reached refinement limit." << endl;
    }
 
    return returnReduce(nRefine-oldNRefine, sumOp<label>());
}
 
 
// Count number of matches of first argument in second argument
Foam::label Foam::meshRefinement::countMatches
(
    const List<point>& normals1,
    const List<point>& normals2,
    const scalar tol
) const
{
    label nMatches = 0;
 
    forAll(normals1, i)
    {
        const vector& n1 = normals1[i];
 
        forAll(normals2, j)
        {
            const vector& n2 = normals2[j];
 
            if (magSqr(n1-n2) < tol)
            {
                nMatches++;
                break;
            }
        }
    }
 
    return nMatches;
}
 
 
// Mark cells for surface curvature based refinement.
Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
(
    const scalar curvature,
    const label nAllowRefine,
    const labelList& neiLevel,
    const pointField& neiCc,
 
    labelList& refineCell,
    label& nRefine
) const
{
    const labelList& cellLevel = meshCutter_.cellLevel();
    const pointField& cellCentres = mesh_.cellCentres();
 
    const label oldNRefine = nRefine;
 
    // 1. local test: any cell on more than one surface gets refined
    // (if its current level is < max of the surface max level)
 
    // 2. 'global' test: any cell on only one surface with a neighbour
    // on a different surface gets refined (if its current level etc.)
 
 
    const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
 
 
    // Collect candidate faces (i.e. intersecting any surface and
    // owner/neighbour not yet refined.
    labelList testFaces(getRefineCandidateFaces(refineCell));
 
    // Collect segments
    pointField start(testFaces.size());
    pointField end(testFaces.size());
    labelList minLevel(testFaces.size());
 
    forAll(testFaces, i)
    {
        const label facei = testFaces[i];
        const label own = mesh_.faceOwner()[facei];
 
        if (mesh_.isInternalFace(facei))
        {
            const label nei = mesh_.faceNeighbour()[facei];
 
            start[i] = cellCentres[own];
            end[i] = cellCentres[nei];
            minLevel[i] = min(cellLevel[own], cellLevel[nei]);
        }
        else
        {
            const label bFacei = facei - mesh_.nInternalFaces();
 
            start[i] = cellCentres[own];
            end[i] = neiCc[bFacei];
 
            if (!isMasterFace[facei])
            {
                Swap(start[i], end[i]);
            }
 
            minLevel[i] = min(cellLevel[own], neiLevel[bFacei]);
        }
    }
 
    // Extend segments a bit
    {
        const vectorField smallVec(rootSmall*(end-start));
        start -= smallVec;
        end += smallVec;
    }
 
 
    // Test for all intersections (with surfaces of higher max level than
    // minLevel) and cache per cell the interesting inter
    labelListList cellSurfLevels(mesh_.nCells());
    List<vectorList> cellSurfNormals(mesh_.nCells());
 
    {
        // Per segment the normals of the surfaces
        List<vectorList> surfaceNormal;
 
        // Per segment the list of levels of the surfaces
        labelListList surfaceLevel;
 
        surfaces_.findAllHigherIntersections
        (
            start,
            end,
            minLevel,           // max level of surface has to be bigger
                                // than min level of neighbouring cells
 
            surfaces_.maxLevel(),
 
            surfaceNormal,
            surfaceLevel
        );
 
 
        // Sort the data according to surface location. This will guarantee
        // that on coupled faces both sides visit the intersections in
        // the same order so will decide the same
        labelList visitOrder;
        forAll(surfaceNormal, pointi)
        {
            vectorList& pNormals = surfaceNormal[pointi];
            labelList& pLevel = surfaceLevel[pointi];
 
            sortedOrder(pNormals, visitOrder, normalLess(pNormals));
 
            pNormals = List<point>(pNormals, visitOrder);
            pLevel = UIndirectList<label>(pLevel, visitOrder);
        }
 
        // Clear out unnecessary data
        start.clear();
        end.clear();
        minLevel.clear();
 
        // Convert face-wise data to cell.
        forAll(testFaces, i)
        {
            const label facei = testFaces[i];
            const label own = mesh_.faceOwner()[facei];
 
            const vectorList& fNormals = surfaceNormal[i];
            const labelList& fLevels = surfaceLevel[i];
 
            forAll(fNormals, hiti)
            {
                if (fLevels[hiti] > cellLevel[own])
                {
                    cellSurfLevels[own].append(fLevels[hiti]);
                    cellSurfNormals[own].append(fNormals[hiti]);
                }
 
                if (mesh_.isInternalFace(facei))
                {
                    label nei = mesh_.faceNeighbour()[facei];
                    if (fLevels[hiti] > cellLevel[nei])
                    {
                        cellSurfLevels[nei].append(fLevels[hiti]);
                        cellSurfNormals[nei].append(fNormals[hiti]);
                    }
                }
            }
        }
    }
 
 
 
    // Bit of statistics
    if (debug)
    {
        label nSet = 0;
        label nNormals = 9;
        forAll(cellSurfNormals, celli)
        {
            const vectorList& normals = cellSurfNormals[celli];
            if (normals.size())
            {
                nSet++;
                nNormals += normals.size();
            }
        }
        reduce(nSet, sumOp<label>());
        reduce(nNormals, sumOp<label>());
        Info<< "markSurfaceCurvatureRefinement :"
            << " cells:" << mesh_.globalData().nTotalCells()
            << " of which with normals:" << nSet
            << " ; total normals stored:" << nNormals
            << endl;
    }
 
 
 
    bool reachedLimit = false;
 
 
    // 1. Check normals per cell
    // ~~~~~~~~~~~~~~~~~~~~~~~~~
 
    for
    (
        label celli = 0;
        !reachedLimit && celli < cellSurfNormals.size();
        celli++
    )
    {
        const vectorList& normals = cellSurfNormals[celli];
        const labelList& levels = cellSurfLevels[celli];
 
        // n^2 comparison of all normals in a cell
        for (label i = 0; !reachedLimit && i < normals.size(); i++)
        {
            for (label j = i+1; !reachedLimit && j < normals.size(); j++)
            {
                if ((normals[i] & normals[j]) < curvature)
                {
                    const label maxLevel = max(levels[i], levels[j]);
 
                    if (cellLevel[celli] < maxLevel)
                    {
                        if
                        (
                            !markForRefine
                            (
                                maxLevel,
                                nAllowRefine,
                                refineCell[celli],
                                nRefine
                            )
                        )
                        {
                            if (debug)
                            {
                                Pout<< "Stopped refining since reaching my cell"
                                    << " limit of " << mesh_.nCells()+7*nRefine
                                    << endl;
                            }
                            reachedLimit = true;
                            break;
                        }
                    }
                }
            }
        }
    }
 
 
 
    // 2. Find out a measure of surface curvature
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Look at normals between neighbouring surfaces
    // Loop over all faces. Could only be checkFaces, except if they're coupled
 
    // Internal faces
    for
    (
        label facei = 0;
        !reachedLimit && facei < mesh_.nInternalFaces();
        facei++
    )
    {
        const label own = mesh_.faceOwner()[facei];
        const label nei = mesh_.faceNeighbour()[facei];
 
        const vectorList& ownNormals = cellSurfNormals[own];
        const labelList& ownLevels = cellSurfLevels[own];
        const vectorList& neiNormals = cellSurfNormals[nei];
        const labelList& neiLevels = cellSurfLevels[nei];
 
 
        // Special case: owner normals set is a subset of the neighbour
        // normals. Do not do curvature refinement since other cell's normals
        // do not add any information. Avoids weird corner extensions of
        // refinement regions.
        bool ownisSubset =
            countMatches(ownNormals, neiNormals)
         == ownNormals.size();
 
        bool neiisSubset =
            countMatches(neiNormals, ownNormals)
         == neiNormals.size();
 
 
        if (!ownisSubset && !neiisSubset)
        {
            // n^2 comparison of between ownNormals and neiNormals
            for (label i = 0; !reachedLimit &&  i < ownNormals.size(); i++)
            {
                for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
                {
                    // Have valid data on both sides. Check curvature.
                    if ((ownNormals[i] & neiNormals[j]) < curvature)
                    {
                        // See which side to refine.
                        if (cellLevel[own] < ownLevels[i])
                        {
                            if
                            (
                                !markForRefine
                                (
                                    ownLevels[i],
                                    nAllowRefine,
                                    refineCell[own],
                                    nRefine
                                )
                            )
                            {
                                if (debug)
                                {
                                    Pout<< "Stopped refining since reaching"
                                        << " my cell limit of "
                                        << mesh_.nCells()+7*nRefine << endl;
                                }
                                reachedLimit = true;
                                break;
                            }
                        }
                        if (cellLevel[nei] < neiLevels[j])
                        {
                            if
                            (
                                !markForRefine
                                (
                                    neiLevels[j],
                                    nAllowRefine,
                                    refineCell[nei],
                                    nRefine
                                )
                            )
                            {
                                if (debug)
                                {
                                    Pout<< "Stopped refining since reaching"
                                        << " my cell limit of "
                                        << mesh_.nCells()+7*nRefine << endl;
                                }
                                reachedLimit = true;
                                break;
                            }
                        }
                    }
                }
            }
        }
    }
 
 
    // Send over surface normal to neighbour cell.
    List<vectorList> neiSurfaceNormals;
    syncTools::swapBoundaryCellList(mesh_, cellSurfNormals, neiSurfaceNormals);
 
    // Boundary faces
    for
    (
        label facei = mesh_.nInternalFaces();
        !reachedLimit && facei < mesh_.nFaces();
        facei++
    )
    {
        label own = mesh_.faceOwner()[facei];
        label bFacei = facei - mesh_.nInternalFaces();
 
        const vectorList& ownNormals = cellSurfNormals[own];
        const labelList& ownLevels = cellSurfLevels[own];
        const vectorList& neiNormals = neiSurfaceNormals[bFacei];
 
        // Special case: owner normals set is a subset of the neighbour
        // normals. Do not do curvature refinement since other cell's normals
        // do not add any information. Avoids weird corner extensions of
        // refinement regions.
        bool ownisSubset =
            countMatches(ownNormals, neiNormals)
         == ownNormals.size();
 
        bool neiisSubset =
            countMatches(neiNormals, ownNormals)
         == neiNormals.size();
 
 
        if (!ownisSubset && !neiisSubset)
        {
            // n^2 comparison of between ownNormals and neiNormals
            for (label i = 0; !reachedLimit &&  i < ownNormals.size(); i++)
            {
                for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
                {
                    // Have valid data on both sides. Check curvature.
                    if ((ownNormals[i] & neiNormals[j]) < curvature)
                    {
                        if (cellLevel[own] < ownLevels[i])
                        {
                            if
                            (
                                !markForRefine
                                (
                                    ownLevels[i],
                                    nAllowRefine,
                                    refineCell[own],
                                    nRefine
                                )
                            )
                            {
                                if (debug)
                                {
                                    Pout<< "Stopped refining since reaching"
                                        << " my cell limit of "
                                        << mesh_.nCells()+7*nRefine
                                        << endl;
                                }
                                reachedLimit = true;
                                break;
                            }
                        }
                    }
                }
            }
        }
    }
 
 
    if
    (
        returnReduce(nRefine, sumOp<label>())
      > returnReduce(nAllowRefine, sumOp<label>())
    )
    {
        Info<< "Reached refinement limit." << endl;
    }
 
    return returnReduce(nRefine-oldNRefine, sumOp<label>());
}
 
 
bool Foam::meshRefinement::isGap
(
    const scalar planarCos,
    const vector& point0,
    const vector& normal0,
 
    const vector& point1,
    const vector& normal1
) const
{
    //- Hits differ and angles are oppositeish and
    //  hits have a normal distance
    const vector d = point1 - point0;
    const scalar magD = mag(d);
 
    if (magD > mergeDistance())
    {
        scalar cosAngle = (normal0 & normal1);
 
        vector avg = Zero;
        if (cosAngle < (-1+planarCos))
        {
            // Opposite normals
            avg = 0.5*(normal0-normal1);
        }
        else if (cosAngle > (1-planarCos))
        {
            avg = 0.5*(normal0+normal1);
        }
 
        if (avg != vector::zero)
        {
            avg /= mag(avg);
 
            // Check normal distance of intersection locations
            if (mag(avg&d) > mergeDistance())
            {
                return true;
            }
            else
            {
                return  false;
            }
        }
        else
        {
            return false;
        }
    }
    else
    {
        return false;
    }
}
 
 
// Mark small gaps
bool Foam::meshRefinement::isNormalGap
(
    const scalar planarCos,
    const vector& point0,
    const vector& normal0,
 
    const vector& point1,
    const vector& normal1
) const
{
    //- Hits differ and angles are oppositeish and
    //  hits have a normal distance
    vector d = point1 - point0;
    const scalar magD = mag(d);
 
    if (magD > mergeDistance())
    {
        scalar cosAngle = (normal0 & normal1);
 
        vector avg = Zero;
        if (cosAngle < (-1+planarCos))
        {
            // Opposite normals
            avg = 0.5*(normal0-normal1);
        }
        else if (cosAngle > (1-planarCos))
        {
            avg = 0.5*(normal0+normal1);
        }
 
        if (avg != vector::zero)
        {
            avg /= mag(avg);
            d /= magD;
 
            // Check average normal with respect to intersection locations
            if (mag(avg&d) > Foam::cos(degToRad(45.0)))
            {
                return true;
            }
            else
            {
                return false;
            }
        }
        else
        {
            return false;
        }
    }
    else
    {
        return false;
    }
}
 
 
bool Foam::meshRefinement::checkProximity
(
    const scalar planarCos,
    const label nAllowRefine,
 
    const label surfaceLevel,       // current intersection max level
    const vector& surfaceLocation,  // current intersection location
    const vector& surfaceNormal,    // current intersection normal
 
    const label celli,
 
    label& cellMaxLevel,        // cached max surface level for this cell
    vector& cellMaxLocation,    // cached surface normal for this cell
    vector& cellMaxNormal,      // cached surface normal for this cell
 
    labelList& refineCell,
    label& nRefine
) const
{
    const labelList& cellLevel = meshCutter_.cellLevel();
 
    // Test if surface applicable
    if (surfaceLevel > cellLevel[celli])
    {
        if (cellMaxLevel == -1)
        {
            // First visit of cell. Store
            cellMaxLevel = surfaceLevel;
            cellMaxLocation = surfaceLocation;
            cellMaxNormal = surfaceNormal;
        }
        else
        {
            // Second or more visit.
            // Check if
            //  - different location
            //  - opposite surface
 
            bool closeSurfaces = isNormalGap
            (
                planarCos,
                cellMaxLocation,
                cellMaxNormal,
                surfaceLocation,
                surfaceNormal
            );
 
            // Set normal to that of highest surface. Not really necessary
            // over here but we reuse cellMax info when doing coupled faces.
            if (surfaceLevel > cellMaxLevel)
            {
                cellMaxLevel = surfaceLevel;
                cellMaxLocation = surfaceLocation;
                cellMaxNormal = surfaceNormal;
            }
 
 
            if (closeSurfaces)
            {
                // Pout<< "Found gap:" << nl
                //    << "    location:" << surfaceLocation
                //    << "\tnormal  :" << surfaceNormal << nl
                //    << "\tnormal  :" << cellMaxNormal << nl
                //    << "\tcos     :" << (surfaceNormal&cellMaxNormal) << nl
                //    << endl;
 
                return markForRefine
                (
                    surfaceLevel,   // mark with any non-neg number.
                    nAllowRefine,
                    refineCell[celli],
                    nRefine
                );
            }
        }
    }
 
    // Did not reach refinement limit.
    return true;
}
 
 
Foam::label Foam::meshRefinement::markProximityRefinement
(
    const scalar planarCos,
    const label nAllowRefine,
    const labelList& neiLevel,
    const pointField& neiCc,
 
    labelList& refineCell,
    label& nRefine
) const
{
    const labelList& cellLevel = meshCutter_.cellLevel();
    const pointField& cellCentres = mesh_.cellCentres();
 
    const label oldNRefine = nRefine;
 
    // 1. local test: any cell on more than one surface gets refined
    // (if its current level is < max of the surface max level)
 
    // 2. 'global' test: any cell on only one surface with a neighbour
    // on a different surface gets refined (if its current level etc.)
 
 
    // Collect candidate faces (i.e. intersecting any surface and
    // owner/neighbour not yet refined.
    labelList testFaces(getRefineCandidateFaces(refineCell));
 
    // Collect segments
    pointField start(testFaces.size());
    pointField end(testFaces.size());
    labelList minLevel(testFaces.size());
 
    forAll(testFaces, i)
    {
        const label facei = testFaces[i];
        const label own = mesh_.faceOwner()[facei];
 
        if (mesh_.isInternalFace(facei))
        {
            const label nei = mesh_.faceNeighbour()[facei];
 
            start[i] = cellCentres[own];
            end[i] = cellCentres[nei];
            minLevel[i] = min(cellLevel[own], cellLevel[nei]);
        }
        else
        {
            const label bFacei = facei - mesh_.nInternalFaces();
 
            start[i] = cellCentres[own];
            end[i] = neiCc[bFacei];
            minLevel[i] = min(cellLevel[own], neiLevel[bFacei]);
        }
    }
 
    // Extend segments a bit
    {
        const vectorField smallVec(rootSmall*(end-start));
        start -= smallVec;
        end += smallVec;
    }
 
 
    // Test for all intersections (with surfaces of higher gap level than
    // minLevel) and cache per cell the max surface level and the local normal
    // on that surface.
    labelList cellMaxLevel(mesh_.nCells(), -1);
    vectorField cellMaxNormal(mesh_.nCells(), Zero);
    pointField cellMaxLocation(mesh_.nCells(), Zero);
 
    {
        // Per segment the normals of the surfaces
        List<vectorList> surfaceLocation;
        List<vectorList> surfaceNormal;
        // Per segment the list of levels of the surfaces
        labelListList surfaceLevel;
 
        surfaces_.findAllHigherIntersections
        (
            start,
            end,
            minLevel,           // gap level of surface has to be bigger
                                // than min level of neighbouring cells
 
            surfaces_.gapLevel(),
 
            surfaceLocation,
            surfaceNormal,
            surfaceLevel
        );
        // Clear out unnecessary data
        start.clear();
        end.clear();
        minLevel.clear();
 
        // OBJstream str
        //(
        //    mesh_.time().path()
        //  / "findAllHigherIntersections_"
        //  + mesh_.time().name()
        //  + ".obj"
        //);
        // OBJstream str2
        //(
        //    mesh_.time().path()
        //  / "findAllHigherIntersections2_"
        //  + mesh_.time().name()
        //  + ".obj"
        //);
 
        forAll(testFaces, i)
        {
            label facei = testFaces[i];
            label own = mesh_.faceOwner()[facei];
 
            const labelList& fLevels = surfaceLevel[i];
            const vectorList& fPoints = surfaceLocation[i];
            const vectorList& fNormals = surfaceNormal[i];
 
            forAll(fLevels, hiti)
            {
                checkProximity
                (
                    planarCos,
                    nAllowRefine,
 
                    fLevels[hiti],
                    fPoints[hiti],
                    fNormals[hiti],
 
                    own,
                    cellMaxLevel[own],
                    cellMaxLocation[own],
                    cellMaxNormal[own],
 
                    refineCell,
                    nRefine
                );
            }
 
            if (mesh_.isInternalFace(facei))
            {
                label nei = mesh_.faceNeighbour()[facei];
 
                forAll(fLevels, hiti)
                {
                    checkProximity
                    (
                        planarCos,
                        nAllowRefine,
 
                        fLevels[hiti],
                        fPoints[hiti],
                        fNormals[hiti],
 
                        nei,
                        cellMaxLevel[nei],
                        cellMaxLocation[nei],
                        cellMaxNormal[nei],
 
                        refineCell,
                        nRefine
                    );
                }
            }
        }
    }
 
    // 2. Find out a measure of surface curvature and region edges.
    // Send over surface region and surface normal to neighbour cell.
 
    labelList neiBndMaxLevel(mesh_.nFaces()-mesh_.nInternalFaces());
    pointField neiBndMaxLocation(mesh_.nFaces()-mesh_.nInternalFaces());
    vectorField neiBndMaxNormal(mesh_.nFaces()-mesh_.nInternalFaces());
 
    for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
    {
        const label bFacei = facei-mesh_.nInternalFaces();
        const label own = mesh_.faceOwner()[facei];
 
        neiBndMaxLevel[bFacei] = cellMaxLevel[own];
        neiBndMaxLocation[bFacei] = cellMaxLocation[own];
        neiBndMaxNormal[bFacei] = cellMaxNormal[own];
    }
    syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLevel);
    syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLocation);
    syncTools::swapBoundaryFaceList(mesh_, neiBndMaxNormal);
 
    // Loop over all faces. Could only be checkFaces.. except if they're coupled
 
    // Internal faces
    for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
    {
        const label own = mesh_.faceOwner()[facei];
        const label nei = mesh_.faceNeighbour()[facei];
 
        if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
        {
            // Have valid data on both sides. Check planarCos.
            if
            (
                isNormalGap
                (
                    planarCos,
                    cellMaxLocation[own],
                    cellMaxNormal[own],
                    cellMaxLocation[nei],
                    cellMaxNormal[nei]
                )
            )
            {
                // See which side to refine
                if (cellLevel[own] < cellMaxLevel[own])
                {
                    if
                    (
                        !markForRefine
                        (
                            cellMaxLevel[own],
                            nAllowRefine,
                            refineCell[own],
                            nRefine
                        )
                    )
                    {
                        if (debug)
                        {
                            Pout<< "Stopped refining since reaching my cell"
                                << " limit of " << mesh_.nCells()+7*nRefine
                                << endl;
                        }
                        break;
                    }
                }
 
                if (cellLevel[nei] < cellMaxLevel[nei])
                {
                    if
                    (
                        !markForRefine
                        (
                            cellMaxLevel[nei],
                            nAllowRefine,
                            refineCell[nei],
                            nRefine
                        )
                    )
                    {
                        if (debug)
                        {
                            Pout<< "Stopped refining since reaching my cell"
                                << " limit of " << mesh_.nCells()+7*nRefine
                                << endl;
                        }
                        break;
                    }
                }
            }
        }
    }
    // Boundary faces
    for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
    {
        const label own = mesh_.faceOwner()[facei];
        const label bFacei = facei - mesh_.nInternalFaces();
 
        if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFacei] != -1)
        {
            // Have valid data on both sides. Check planarCos.
            if
            (
                isNormalGap
                (
                    planarCos,
                    cellMaxLocation[own],
                    cellMaxNormal[own],
                    neiBndMaxLocation[bFacei],
                    neiBndMaxNormal[bFacei]
                )
            )
            {
                if
                (
                    !markForRefine
                    (
                        cellMaxLevel[own],
                        nAllowRefine,
                        refineCell[own],
                        nRefine
                    )
                )
                {
                    if (debug)
                    {
                        Pout<< "Stopped refining since reaching my cell"
                            << " limit of " << mesh_.nCells()+7*nRefine
                            << endl;
                    }
                    break;
                }
            }
        }
    }
 
    if
    (
        returnReduce(nRefine, sumOp<label>())
      > returnReduce(nAllowRefine, sumOp<label>())
    )
    {
        Info<< "Reached refinement limit." << endl;
    }
 
    return returnReduce(nRefine-oldNRefine, sumOp<label>());
}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
// Calculate list of cells to refine. Gets for any edge (start - end)
// whether it intersects the surface. Does more accurate test and checks
// the wanted level on the surface intersected.
// Does approximate precalculation of how many cells can be refined before
// hitting overall limit maxGlobalCells.
Foam::labelList Foam::meshRefinement::refineCandidates
(
    const List<point>& insidePoints,
    const scalar curvature,
    const scalar planarAngle,
 
    const bool featureRefinement,
    const bool featureDistanceRefinement,
    const bool internalRefinement,
    const bool surfaceRefinement,
    const bool curvatureRefinement,
    const bool gapRefinement,
    const label maxGlobalCells,
    const label maxLocalCells
) const
{
    const label totNCells = mesh_.globalData().nTotalCells();
 
    labelList cellsToRefine;
 
    if (totNCells >= maxGlobalCells)
    {
        Info<< "No cells marked for refinement since reached limit "
            << maxGlobalCells << '.' << endl;
    }
    else
    {
        // Every cell I refine adds 7 cells. Estimate the number of cells
        // I am allowed to refine.
        // Assume perfect distribution so can only refine as much the fraction
        // of the mesh I hold. This prediction step prevents us having to do
        // lots of reduces to keep count of the total number of cells selected
        // for refinement.
 
        // scalar fraction = scalar(mesh_.nCells())/totNCells;
        // label myMaxCells = label(maxGlobalCells*fraction);
        // label nAllowRefine = (myMaxCells - mesh_.nCells())/7;
        //
        // Pout<< "refineCandidates:" << nl
        //    << "    total cells:" << totNCells << nl
        //    << "    local cells:" << mesh_.nCells() << nl
        //    << "    local fraction:" << fraction << nl
        //    << "    local allowable cells:" << myMaxCells << nl
        //    << "    local allowable refinement:" << nAllowRefine << nl
        //    << endl;
 
        //- Disable refinement shortcut. nAllowRefine is per processor limit.
        const label nAllowRefine = labelMax / Pstream::nProcs();
 
        // Marked for refinement (>= 0) or not (-1). Actual value is the
        // index of the surface it intersects.
        labelList refineCell(mesh_.nCells(), -1);
        label nRefine = 0;
 
 
        // Swap neighbouring cell centres and cell level
        labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
        pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
        calcNeighbourData(neiLevel, neiCc);
 
 
 
        // Cells pierced by feature edges
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
        if (featureRefinement)
        {
            const label nFeatures = markFeatureRefinement
            (
                insidePoints,
                nAllowRefine,
 
                refineCell,
                nRefine
            );
 
            Info<< "Marked for refinement due to explicit features             "
                << ": " << nFeatures << " cells."  << endl;
        }
 
        // Inside distance-to-feature
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
        if (featureDistanceRefinement)
        {
            const label nCellsFeat = markInternalDistanceToFeatureRefinement
            (
                nAllowRefine,
 
                refineCell,
                nRefine
            );
            Info<< "Marked for refinement due to distance to explicit features "
                ": " << nCellsFeat << " cells."  << endl;
        }
 
        // Inside refinement shells
        // ~~~~~~~~~~~~~~~~~~~~~~~~
 
        if (internalRefinement)
        {
            const label nShell = markInternalRefinement
            (
                nAllowRefine,
 
                refineCell,
                nRefine
            );
            Info<< "Marked for refinement due to refinement shells             "
                << ": " << nShell << " cells."  << endl;
        }
 
        // Refinement based on intersection of surface
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
        if (surfaceRefinement)
        {
            const label nSurf = markSurfaceRefinement
            (
                nAllowRefine,
                neiLevel,
                neiCc,
 
                refineCell,
                nRefine
            );
            Info<< "Marked for refinement due to surface intersection          "
                << ": " << nSurf << " cells."  << endl;
        }
 
        // Refinement based on curvature of surface
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
        if
        (
            curvatureRefinement
         && (curvature >= -1 && curvature <= 1)
         && (surfaces_.minLevel() != surfaces_.maxLevel())
        )
        {
            const label nCurv = markSurfaceCurvatureRefinement
            (
                curvature,
                nAllowRefine,
                neiLevel,
                neiCc,
 
                refineCell,
                nRefine
            );
            Info<< "Marked for refinement due to curvature/regions             "
                << ": " << nCurv << " cells."  << endl;
        }
 
 
        const scalar planarCos = Foam::cos(planarAngle);
 
        if
        (
            gapRefinement
         && (planarCos >= -1 && planarCos <= 1)
         && (max(surfaces_.gapLevel()) > -1)
        )
        {
            Info<< "Specified gap level : " << max(surfaces_.gapLevel())
                << ", planar angle " << radToDeg(planarAngle) << endl;
 
            const label nGap = markProximityRefinement
            (
                planarCos,
                nAllowRefine,
                neiLevel,
                neiCc,
 
                refineCell,
                nRefine
            );
            Info<< "Marked for refinement due to close opposite surfaces       "
                << ": " << nGap << " cells."  << endl;
        }
 
 
 
        // Pack cells-to-refine
        // ~~~~~~~~~~~~~~~~~~~~
 
        cellsToRefine.setSize(nRefine);
        nRefine = 0;
 
        forAll(refineCell, celli)
        {
            if (refineCell[celli] != -1)
            {
                cellsToRefine[nRefine++] = celli;
            }
        }
    }
 
    return cellsToRefine;
}
 
 
Foam::autoPtr<Foam::polyTopoChangeMap> Foam::meshRefinement::refine
(
    const labelList& cellsToRefine
)
{
    // Mesh changing engine.
    polyTopoChange meshMod(mesh_);
 
    // Play refinement commands into mesh changer.
    meshCutter_.setRefinement(cellsToRefine, meshMod);
 
    // Create mesh
    // return map from old to new mesh.
    autoPtr<polyTopoChangeMap> map = meshMod.changeMesh(mesh_);
 
    // Update fields
    mesh_.topoChange(map);
 
    // Reset the instance for if in overwrite mode
    mesh_.setInstance(name());
 
    // Update intersection info
    topoChange(map, getChangedFaces(map, cellsToRefine));
 
    return map;
}
 
 
Foam::autoPtr<Foam::polyDistributionMap>
Foam::meshRefinement::refineAndBalance
(
    const string& msg,
    decompositionMethod& decomposer,
    fvMeshDistribute& distributor,
    const labelList& cellsToRefine,
    const scalar maxLoadUnbalance
)
{
    // Do all refinement
    refine(cellsToRefine);
 
    if (debug&meshRefinement::MESH)
    {
        Pout<< "Writing refined but unbalanced " << msg
            << " mesh to time " << name() << endl;
        write
        (
            debugType(debug),
            writeType(writeLevel() | WRITEMESH),
            mesh_.time().path()/name()
        );
        Pout<< "Dumped debug data in = "
            << mesh_.time().cpuTimeIncrement() << " s" << endl;
 
        // test all is still synced across proc patches
        checkData();
    }
 
    Info<< "Refined mesh in = "
        << mesh_.time().cpuTimeIncrement() << " s" << endl;
    printMeshInfo(debug, "After refinement " + msg);
 
 
    // Load balancing
    // ~~~~~~~~~~~~~~
 
    autoPtr<polyDistributionMap> distMap;
 
    if (Pstream::nProcs() > 1)
    {
        const scalar nIdealCells =
            mesh_.globalData().nTotalCells()
          / Pstream::nProcs();
 
        const scalar unbalance = returnReduce
        (
            mag(1.0-mesh_.nCells()/nIdealCells),
            maxOp<scalar>()
        );
 
        if (unbalance <= maxLoadUnbalance)
        {
            Info<< "Skipping balancing since max unbalance " << unbalance
                << " is less than allowable " << maxLoadUnbalance
                << endl;
        }
        else
        {
            scalarField cellWeights(mesh_.nCells(), 1);
 
            distMap = balance
            (
                false,  // keepZoneFaces
                false,  // keepBaffles
                cellWeights,
                decomposer,
                distributor
            );
 
            Info<< "Balanced mesh in = "
                << mesh_.time().cpuTimeIncrement() << " s" << endl;
 
            printMeshInfo(debug, "After balancing " + msg);
 
 
            if (debug&meshRefinement::MESH)
            {
                Pout<< "Writing balanced " << msg
                    << " mesh to time " << name() << endl;
                write
                (
                    debugType(debug),
                    writeType(writeLevel() | WRITEMESH),
                    mesh_.time().path()/name()
                );
                Pout<< "Dumped debug data in = "
                    << mesh_.time().cpuTimeIncrement() << " s" << endl;
 
                // test all is still synced across proc patches
                checkData();
            }
        }
    }
 
    return distMap;
}
 
 
// Do load balancing followed by refinement of consistent set of cells.
Foam::autoPtr<Foam::polyDistributionMap>
Foam::meshRefinement::balanceAndRefine
(
    const string& msg,
    decompositionMethod& decomposer,
    fvMeshDistribute& distributor,
    const labelList& initCellsToRefine,
    const scalar maxLoadUnbalance
)
{
    labelList cellsToRefine(initCellsToRefine);
 
    //{
    //    globalIndex globalCells(mesh_.nCells());
    //
    //    Info<< "** Distribution before balancing/refining:" << endl;
    //    for (label proci = 0; proci < Pstream::nProcs(); proci++)
    //    {
    //        Info<< "    " << proci << '\t'
    //            << globalCells.localSize(proci) << endl;
    //    }
    //    Info<< endl;
    //}
    //{
    //    globalIndex globalCells(cellsToRefine.size());
    //
    //    Info<< "** Cells to be refined:" << endl;
    //    for (label proci = 0; proci < Pstream::nProcs(); proci++)
    //    {
    //        Info<< "    " << proci << '\t'
    //            << globalCells.localSize(proci) << endl;
    //    }
    //    Info<< endl;
    //}
 
 
    // Load balancing
    // ~~~~~~~~~~~~~~
 
    autoPtr<polyDistributionMap> distMap;
 
    if (Pstream::nProcs() > 1)
    {
        // First check if we need to balance at all. Precalculate number of
        // cells after refinement and see what maximum difference is.
        const scalar nNewCells =
            scalar(mesh_.nCells() + 7*cellsToRefine.size());
        const scalar nIdealNewCells =
            returnReduce(nNewCells, sumOp<scalar>())/Pstream::nProcs();
        const scalar unbalance = returnReduce
        (
            mag(1.0-nNewCells/nIdealNewCells),
            maxOp<scalar>()
        );
 
        if (unbalance <= maxLoadUnbalance)
        {
            Info<< "Skipping balancing since max unbalance " << unbalance
                << " is less than allowable " << maxLoadUnbalance
                << endl;
        }
        else
        {
            scalarField cellWeights(mesh_.nCells(), 1);
            forAll(cellsToRefine, i)
            {
                cellWeights[cellsToRefine[i]] += 7;
            }
 
            distMap = balance
            (
                false,  // keepZoneFaces
                false,  // keepBaffles
                cellWeights,
                decomposer,
                distributor
            );
 
            // Update cells to refine
            distMap().distributeCellIndices(cellsToRefine);
 
            Info<< "Balanced mesh in = "
                << mesh_.time().cpuTimeIncrement() << " s" << endl;
        }
 
        //{
        //    globalIndex globalCells(mesh_.nCells());
        //
        //    Info<< "** Distribution after balancing:" << endl;
        //    for (label proci = 0; proci < Pstream::nProcs(); proci++)
        //    {
        //        Info<< "    " << proci << '\t'
        //            << globalCells.localSize(proci) << endl;
        //    }
        //    Info<< endl;
        //}
 
        printMeshInfo(debug, "After balancing " + msg);
 
        if (debug&meshRefinement::MESH)
        {
            Pout<< "Writing balanced " << msg
                << " mesh to time " << name() << endl;
            write
            (
                debugType(debug),
                writeType(writeLevel() | WRITEMESH),
                mesh_.time().path()/name()
            );
            Pout<< "Dumped debug data in = "
                << mesh_.time().cpuTimeIncrement() << " s" << endl;
 
            // test all is still synced across proc patches
            checkData();
        }
    }
 
 
    // Refinement
    // ~~~~~~~~~~
 
    refine(cellsToRefine);
 
    if (debug&meshRefinement::MESH)
    {
        Pout<< "Writing refined " << msg
            << " mesh to time " << name() << endl;
        write
        (
            debugType(debug),
            writeType(writeLevel() | WRITEMESH),
            mesh_.time().path()/name()
        );
        Pout<< "Dumped debug data in = "
            << mesh_.time().cpuTimeIncrement() << " s" << endl;
 
        // test all is still synced across proc patches
        checkData();
    }
 
    Info<< "Refined mesh in = "
        << mesh_.time().cpuTimeIncrement() << " s" << endl;
 
    //{
    //    globalIndex globalCells(mesh_.nCells());
    //
    //    Info<< "** After refinement distribution:" << endl;
    //    for (label proci = 0; proci < Pstream::nProcs(); proci++)
    //    {
    //        Info<< "    " << proci << '\t'
    //            << globalCells.localSize(proci) << endl;
    //    }
    //    Info<< endl;
    //}
 
    printMeshInfo(debug, "After refinement " + msg);
 
    return distMap;
}
 
 
// ************************************************************************* //