/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  dev
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    object      snappyHexMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// Which of the steps to run
castellatedMesh true;
snap            true;
addLayers       false;


//Optional: single region surfaces get patch names according to
//          surface only. Multi-region surfaces get patch name
//          surface "_ "region. Default is true
//singleRegionName false;


//Optional: preserve all generated patches. Default is to remove
//          zero-sized patches.
//keepPatches true;


// Geometry. Definition of all surfaces. All surfaces are of class
// searchableSurface.
// Surfaces are used
// - to specify refinement for any mesh cell intersecting it
// - to specify refinement for any mesh cell inside/outside/near
// - to 'snap' the mesh boundary to the surface
geometry
{
    box1x1x1
    {
        type searchableBox;
        min (1.5 1 -0.5);
        max (3.5 2 0.5);
    }

    sphere
    {
        type triSurfaceMesh;
        file "sphere.stl"

        // Optional: non-default tolerance on intersections
        // tolerance 1e-5;

        // Optional: depth of octree. Decrease only in case of memory
        // limitations.
        // maxTreeDepth 10;

        // Optional: scaling factor, e.g. unit conversion
        // scale 1;

        // Optional: discard triangles with low quality when getting normal
        // minQuality -1;

        // Per region the patchname. If not provided will be <surface>_<region>.
        // Note: this name cannot be used to identity this region in any
        //       other part of this dictionary; it is only a name
        //       for the combination of surface+region (which is only used
        //       when creating patches)
        regions
        {
            secondSolid
            {
                name mySecondPatch;
            }
        }
    }

    sphere2
    {
        type searchableSphere;
        centre  (1.5 1.5 1.5);
        radius  1.03;
    }
};

// Settings for the castellatedMesh generation.
castellatedMeshControls
{

    // Refinement parameters
    // ~~~~~~~~~~~~~~~~~~~~~

    // If local number of cells is >= maxLocalCells on any processor
    // switches from from refinement followed by balancing
    // (current method) to (weighted) balancing before refinement.
    maxLocalCells 100000;

    // Overall cell limit (approximately). Refinement will stop immediately
    // upon reaching this number so a refinement level might not complete.
    // Note that this is the number of cells before removing the part which
    // is not 'visible' from the keepPoint. The final number of cells might
    // actually be a lot less.
    maxGlobalCells 2000000;

    // The surface refinement loop might spend lots of iterations refining just a
    // few cells. This setting will cause refinement to stop if <= minimumRefine
    // are selected for refinement. Note: it will at least do one iteration
    // (unless the number of cells to refine is 0)
    minRefinementCells 0;

    // Allow a certain level of imbalance during refining
    // (since balancing is quite expensive)
    // Expressed as fraction of perfect balance (= overall number of cells /
    // nProcs). 0=balance always.
    maxLoadUnbalance 0.10;

    // Number of buffer layers between different levels.
    // 1 means normal 2:1 refinement restriction, larger means slower
    // refinement.
    nCellsBetweenLevels 1;


    // Explicit feature edge refinement
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    // Specifies a level for any cell intersected by explicitly provided
    // edges.
    // This is a featureEdgeMesh, read from constant/geometry for now.
    // Specify 'levels' in the same way as the 'distance' mode in the
    // refinementRegions (see below). The old specification
    //      level   2;
    // is equivalent to
    //      levels  ((0 2));

    features
    (
        //{
        //    file "someLine.eMesh";
        //    // level 2;
        //    levels ((0.0 2) (1.0 3));
        //}
    );


    // Surface based refinement
    // ~~~~~~~~~~~~~~~~~~~~~~~~

    // Specifies two levels for every surface. The first is the minimum level,
    // every cell intersecting a surface gets refined up to the minimum level.
    // The second level is the maximum level. Cells that 'see' multiple
    // intersections where the intersections make an
    // angle > resolveFeatureAngle get refined up to the maximum level.

    refinementSurfaces
    {
        sphere
        {
            // Surface-wise min and max refinement level
            level (2 2);

            // Optional region-wise level specification
            regions
            {
                secondSolid
                {
                    level (3 3);
                }
            }

            // Optional specification of patch type (default is wall). No
            // constraint types (cyclic, symmetry) etc. are allowed.
            patchInfo
            {
                type patch;
                inGroups (meshedPatches);
            }


            //- Optional increment (on top of max level) in small gaps
            // gapLevelIncrement 2;

            //- Optional angle to detect small-large cell situation
            //  perpendicular to the surface. Is the angle of face w.r.t.
            //  the local surface normal. Use on flat(ish) surfaces only.
            //  Otherwise leave out or set to negative number.
            // perpendicularAngle 10;

            //- Optional faceZone and (for closed surface) cellZone with
            //  how to select the cells that are in the cellZone
            //  (inside / outside / specified insidePoint)
            //  The orientation of the faceZone is
            //  - if on cellZone(s) : point out of (maximum) cellZone
            //  - if freestanding   : oriented according to surface

            // faceZone sphere;
            // cellZone sphere;
            // mode        inside;  // outside/insidePoint

            //- Optional specification of what to do with faceZone faces:
            //      internal : keep them as internal faces (default)
            //      baffle   : create baffles from them. This gives more
            //                 freedom in mesh motion
            //      boundary : create free-standing boundary faces (baffles
            //                 but without the shared points)
            // faceType baffle;
        }
    }

    // Feature angle:
    // - used if min and max refinement level of a surface differ
    // - used if feature snapping (see snapControls below) is used
    resolveFeatureAngle 30;

    //- Optional increment (on top of max level) in small gaps
    // gapLevelIncrement 2;


    // Planar angle:
    // - used to determine if surface normals
    //   are roughly the same or opposite. Used
    //      - in proximity refinement
    //      - to decide when to merge free-standing baffles
    //        (if e.g. running in surfaceSimplify mode set this to 180 to
    //         merge all baffles)
    //      - in snapping to avoid snapping to nearest on 'wrong' side
    //        of thin gap
    //
    // If not specified same as resolveFeatureAngle
    planarAngle 30;


    // Region-wise refinement
    // ~~~~~~~~~~~~~~~~~~~~~~

    // Specifies refinement level for cells in relation to a surface. One of
    // three modes
    // - distance. 'levels' specifies per distance to the surface the
    //   wanted refinement level. The distances need to be specified in
    //   increasing order.
    // - inside. 'levels' is only one entry and only the level is used. All
    //   cells inside the surface get refined up to the level. The surface
    //   needs to be closed for this to be possible.
    // - outside. Same but cells outside.

    refinementRegions
    {
        box1x1x1
        {
            mode inside;
            levels ((1.0 4));
        }
        // sphere
        //{
        //    mode distance;
        //    levels ((1.0 5) (2.0 3));
        //}
    }

    // Mesh selection
    // ~~~~~~~~~~~~~~

    // After refinement patches get added for all refinementSurfaces and
    // all cells intersecting the surfaces get put into these patches. The
    // section reachable from the insidePoint is kept.
    // NOTE: This point should never be on a face, always inside a cell, even
    // after refinement.
    insidePoint (5 0.28 0.43);

    // Whether any faceZones (as specified in the refinementSurfaces)
    // are only on the boundary of corresponding cellZones or also allow
    // free-standing zone faces. Not used if there are no faceZones.
    allowFreeStandingZoneFaces true;


    // Optional: do not remove cells likely to give snapping problems
    // handleSnapProblems false;

    // Optional: switch off topological test for cells to-be-squashed
    //           and use geometric test instead
    // useTopologicalSnapDetection false;

    // Extend refinement regions to the span of the surface facet/triangles
    extendedRefinementSpan true;
}

// Settings for the snapping.
snapControls
{
    // Number of patch smoothing iterations before finding correspondence
    // to surface
    nSmoothPatch 3;

    // Maximum relative distance for points to be attracted by surface.
    // True distance is this factor times local maximum edge length.
    // Note: changed(corrected) w.r.t 17x! (17x used 2* tolerance)
    tolerance 2.0;

    // Number of mesh displacement relaxation iterations.
    nSolveIter 30;

    // Maximum number of snapping relaxation iterations. Should stop
    // before upon reaching a correct mesh.
    nRelaxIter 5;

    // Feature snapping

        // Number of feature edge snapping iterations.
        // Leave out altogether to disable.
        nFeatureSnapIter 10;

        // Detect (geometric only) features by sampling the surface
        // (default=false).
        implicitFeatureSnap false;

        // Use castellatedMeshControls::features (default = true)
        explicitFeatureSnap true;

        // Detect features between multiple surfaces
        // (only for explicitFeatureSnap, default = false)
        multiRegionFeatureSnap false;


    // wip: disable snapping to opposite near surfaces (revert to 22x behaviour)
    // detectNearSurfacesSnap false;
}

// Settings for the layer addition.
addLayersControls
{
    // Are the thickness parameters below relative to the undistorted
    // size of the refined cell outside layer (true) or absolute sizes (false).
    relativeSizes true;

    // Layer thickness specification. This can be specified in one of following
    // ways:
    // - expansionRatio and finalLayerThickness (cell nearest internal mesh)
    // - expansionRatio and firstLayerThickness (cell on surface)
    // - overall thickness and firstLayerThickness
    // - overall thickness and finalLayerThickness
    // - overall thickness and expansionRatio
    //
    // Note: the mode thus selected is global, i.e. one cannot override the
    //       mode on a per-patch basis (only the values can be overridden)

        // Expansion factor for layer mesh
        expansionRatio 1.0;

        // Wanted thickness of the layer furthest away from the wall.
        // If relativeSizes this is relative to undistorted size of cell
        // outside layer.
        finalLayerThickness 0.3;

        // Wanted thickness of the layer next to the wall.
        // If relativeSizes this is relative to undistorted size of cell
        // outside layer.
        // firstLayerThickness 0.3;

        // Wanted overall thickness of layers.
        // If relativeSizes this is relative to undistorted size of cell
        // outside layer.
        // thickness 0.5


    // Minimum overall thickness of total layers. If for any reason layer
    // cannot be above minThickness do not add layer.
    // If relativeSizes this is relative to undistorted size of cell
    // outside layer..
    minThickness 0.25;


    // Per final patch (so not geometry!) the layer information
    // Note: This behaviour changed after 21x. Any non-mentioned patches
    //       now slide unless:
    //          - nSurfaceLayers is explicitly mentioned to be 0.
    //          - angle to nearest surface < slipFeatureAngle (see below)
    layers
    {
        sphere_firstSolid
        {
            nSurfaceLayers 1;

        }
        maxY
        {
            nSurfaceLayers 1;
            // Per patch layer data
            expansionRatio      1.3;
            finalLayerThickness 0.3;
            minThickness        0.1;
        }

        // Disable any mesh shrinking and layer addition on any point of
        // a patch by setting nSurfaceLayers to 0
        frozenPatches
        {
            nSurfaceLayers 0;
        }
    }

    // If points get not extruded do nGrow layers of connected faces that are
    // also not grown. This helps convergence of the layer addition process
    // close to features.
    // Note: changed(corrected) w.r.t 17x! (didn't do anything in 17x)
    nGrow 0;

    // Advanced settings


    // Static analysis of starting mesh

        // When not to extrude surface. 0 is flat surface, 90 is when two faces
        // are perpendicular
        featureAngle 130;

        // Stop layer growth on highly warped cells
        maxFaceThicknessRatio 0.5;


    // Patch displacement

        // Number of smoothing iterations of surface normals
        nSmoothSurfaceNormals 1;

        // Smooth layer thickness over surface patches
        nSmoothThickness 10;



    // Medial axis analysis

        // Angle used to pick up medial axis points
        // Note: changed(corrected) w.r.t 17x! 90 degrees corresponds to 130
        // in 17x.
        minMedialAxisAngle 90;

        // Reduce layer growth where ratio thickness to medial
        // distance is large
        maxThicknessToMedialRatio 0.3;

        // Number of smoothing iterations of interior mesh movement direction
        nSmoothNormals 3;

        // Optional: limit the number of steps walking away from the surface.
        // Default is unlimited.
        // nMedialAxisIter 10;

        // Optional: smooth displacement after medial axis determination.
        // default is 0.
        // nSmoothDisplacement 90;

        // (wip)Optional: do not extrude a point if none of the surrounding points is
        // not extruded. Default is false.
        // detectExtrusionIsland true;


    // Mesh shrinking

        // Optional: at non-patched sides allow mesh to slip if extrusion
        // direction makes angle larger than slipFeatureAngle. Default is
        // 0.5*featureAngle.
        slipFeatureAngle 30;

        // Maximum number of snapping relaxation iterations. Should stop
        // before upon reaching a correct mesh.
        nRelaxIter 5;

        // Create buffer region for new layer terminations
        nBufferCellsNoExtrude 0;

        // Overall max number of layer addition iterations. The mesher will
        // exit if it reaches this number of iterations; possibly with an
        // illegal mesh.
        nLayerIter 50;

        // Max number of iterations after which relaxed meshQuality controls
        // get used. Up to nRelaxedIter it uses the settings in
        // meshQualityControls,
        // after nRelaxedIter it uses the values in
        // meshQualityControls::relaxed.
        nRelaxedIter 20;

        // Additional reporting: if there are just a few faces where there
        // are mesh errors (after adding the layers) print their face centres.
        // This helps in tracking down problematic mesh areas.
        // additionalReporting true;
}

// Generic mesh quality settings. At any undoable phase these determine
// where to undo.
meshQualityControls
{
    // Specify mesh quality constraints in separate dictionary so can
    // be reused (e.g. checkMesh -meshQuality)
    #include "meshQualityDict"


    // Optional : some meshing phases allow usage of relaxed rules.
    // See e.g. addLayersControls::nRelaxedIter.
    relaxed
    {
        // Maximum non-orthogonality allowed. Set to 180 to disable.
        maxNonOrtho 75;
    }
}

// Advanced

//// Debug flags
//debugFlags
//(
//    mesh            // write intermediate meshes
//    intersections   // write current mesh intersections as .obj files
//    featureSeeds    // write information about explicit feature edge
//                    // refinement
//    attraction      // write attraction as .obj files
//    layerInfo       // write information about layers
//);
//
//// Write flags
//writeFlags
//(
//    scalarLevels    // write volScalarField with cellLevel for postprocessing
//    layerSets       // write cellSets, faceSets of faces in layer
//    layerFields     // write volScalarField for layer coverage
//);

// Merge tolerance. Is fraction of overall bounding box of initial mesh.
// Note: the write tolerance needs to be higher than this.
mergeTolerance 1e-6;

// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2022 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 "mapDistributePolyMesh.H"
#include "Cloud.H"
#include "OBJstream.H"
#include "cellSet.H"
#include "treeDataCell.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 mapPolyMesh& 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.boundaryMesh(), patchi)
        {
            const polyPatch& pp = mesh.boundaryMesh()[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.
    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.

    forAll(insidePoints, i)
    {
        const point& insidePoint = insidePoints[i];

        const label celli = mesh_.cellTree().findInside(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
                            (
                                mesh_,
                                insidePoint,
                                celli,
                                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
                            (
                                mesh_,
                                insidePoint,
                                celli,
                                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,
        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.
    scalar maxTrackLen = 2.0*mesh_.bounds().mag();

    if (debug&meshRefinement::FEATURESEEDS)
    {
        Pout<< "Tracking " << startPointCloud.size()
            << " particles over distance " << maxTrackLen
            << " to find the starting cell" << endl;
    }
    startPointCloud.move(startPointCloud, td, maxTrackLen);


    // Reset levels
    maxFeatureLevel = -1;
    forAll(features_, feati)
    {
        featureEdgeVisited[feati] = 0u;
    }


    Cloud<trackedParticle> cloud
    (
        mesh_,
        "featureCloud",
        IDLList<trackedParticle>()
    );

    if (debug&meshRefinement::FEATURESEEDS)
    {
        Pout<< "Constructing cloud for cell marking" << endl;
    }

    forAllIter(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();
                tp->end() = featureMesh.points()[otherPointi];
                tp->j() = otherPointi;
                tp->k() = edgei;

                if (debug&meshRefinement::FEATURESEEDS)
                {
                    Pout<< "Adding particle for point:" << pointi
                        << " coord:" << tp->position()
                        << " 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 " << maxTrackLen
                << " to mark cells" << endl;
        }
        cloud.move(cloud, td, maxTrackLen);

        // Make particle follow edge.
        forAllIter(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();
                    tp.end() = featureMesh.points()[otherPointi];
                    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]
    //                << " leve:" << 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
                ///    << "    location:" << cellMaxLocation
                //    << "\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();

        //// Extract per cell information on the surface with the highest max
        // OBJstream str
        //(
        //    mesh_.time().path()
        //  / "findAllHigherIntersections_"
        //  + mesh_.time().timeName()
        //  + ".obj"
        //);
        //// All intersections
        // OBJstream str2
        //(
        //    mesh_.time().path()
        //  / "findAllHigherIntersections2_"
        //  + mesh_.time().timeName()
        //  + ".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;
        ////label nAllowRefine = (maxLocalCells - 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(degToRad(planarAngle));

        if
        (
            gapRefinement
         && (planarCos >= -1 && planarCos <= 1)
         && (max(surfaces_.gapLevel()) > -1)
        )
        {
            Info<< "Specified gap level : " << max(surfaces_.gapLevel())
                << ", planar angle " << 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::mapPolyMesh> Foam::meshRefinement::refine
(
    const labelList& cellsToRefine
)
{
    // Mesh changing engine.
    polyTopoChange meshMod(mesh_);

    // Play refinement commands into mesh changer.
    meshCutter_.setRefinement(cellsToRefine, meshMod);

    // Create mesh (no inflation), return map from old to new mesh.
    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false);

    // Update fields
    mesh_.updateMesh(map);

    // Optionally inflate mesh
    if (map().hasMotionPoints())
    {
        mesh_.movePoints(map().preMotionPoints());
    }
    else
    {
        // Delete mesh volumes.
        mesh_.clearOut();
    }

    // Reset the instance for if in overwrite mode
    mesh_.setInstance(timeName());

    // Update intersection info
    updateMesh(map, getChangedFaces(map, cellsToRefine));

    return map;
}


Foam::autoPtr<Foam::mapDistributePolyMesh>
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 " << timeName() << endl;
        write
        (
            debugType(debug),
            writeType(writeLevel() | WRITEMESH),
            mesh_.time().path()/timeName()
        );
        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<mapDistributePolyMesh> 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 " << timeName() << endl;
                write
                (
                    debugType(debug),
                    writeType(writeLevel() | WRITEMESH),
                    mesh_.time().path()/timeName()
                );
                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::mapDistributePolyMesh>
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<mapDistributePolyMesh> 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 " << timeName() << endl;
            write
            (
                debugType(debug),
                writeType(writeLevel() | WRITEMESH),
                mesh_.time().path()/timeName()
            );
            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 " << timeName() << endl;
        write
        (
            debugType(debug),
            writeType(writeLevel() | WRITEMESH),
            mesh_.time().path()/timeName()
        );
        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;
}


// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2022 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 "refinementFeatures.H"
#include "Time.H"
#include "Tuple2.H"
#include "DynamicField.H"
#include "featureEdgeMesh.H"

// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

void Foam::refinementFeatures::read
(
    const objectRegistry& io,
    const PtrList<dictionary>& featDicts
)
{
    forAll(featDicts, feati)
    {
        const dictionary& dict = featDicts[feati];

        fileName featFileName(dict.lookup("file"));


        // Try reading extendedEdgeMesh first

        typeIOobject<extendedFeatureEdgeMesh> extFeatObj
        (
            featFileName,                       // name
            io.time().constant(),               // instance
            "extendedFeatureEdgeMesh",          // local
            io.time(),                          // registry
            IOobject::MUST_READ,
            IOobject::NO_WRITE,
            false
        );

        const fileName fName(extFeatObj.filePath());

        if (!fName.empty() && extendedEdgeMesh::canRead(fName))
        {
            autoPtr<extendedEdgeMesh> eMeshPtr = extendedEdgeMesh::New
            (
                fName
            );

            Info<< "Read extendedFeatureEdgeMesh " << extFeatObj.name()
                << nl << incrIndent;
            eMeshPtr().writeStats(Info);
            Info<< decrIndent << endl;

            set(feati, new extendedFeatureEdgeMesh(extFeatObj, eMeshPtr()));
        }
        else
        {
            // Try reading edgeMesh

            typeIOobject<featureEdgeMesh> featObj
            (
                featFileName,
                io.time().constant(),
                searchableSurface::geometryDir(io.time()),
                io.time(),
                IOobject::MUST_READ,
                IOobject::NO_WRITE,
                false
            );

            const fileName fName(featObj.filePath());

            if (fName.empty())
            {
                FatalIOErrorInFunction
                (
                    dict
                )   << "Could not open "
                    << featObj.objectPath()
                    << exit(FatalIOError);
            }


            // Read as edgeMesh
            autoPtr<edgeMesh> eMeshPtr = edgeMesh::New(fName);
            const edgeMesh& eMesh = eMeshPtr();

            Info<< "Read edgeMesh " << featObj.name() << nl
                << incrIndent;
            eMesh.writeStats(Info);
            Info<< decrIndent << endl;


            // Analyse for feature points. These are all classified as mixed
            // points for lack of anything better
            const labelListList& pointEdges = eMesh.pointEdges();

            labelList oldToNew(eMesh.points().size(), -1);
            DynamicField<point> newPoints(eMesh.points().size());
            forAll(pointEdges, pointi)
            {
                if (pointEdges[pointi].size() > 2)
                {
                    oldToNew[pointi] = newPoints.size();
                    newPoints.append(eMesh.points()[pointi]);
                }
                // else if (pointEdges[pointi].size() == 2)
                // MEJ: do something based on a feature angle?
            }
            label nFeatures = newPoints.size();
            forAll(oldToNew, pointi)
            {
                if (oldToNew[pointi] == -1)
                {
                    oldToNew[pointi] = newPoints.size();
                    newPoints.append(eMesh.points()[pointi]);
                }
            }


            const edgeList& edges = eMesh.edges();
            edgeList newEdges(edges.size());
            forAll(edges, edgei)
            {
                const edge& e = edges[edgei];
                newEdges[edgei] = edge
                (
                    oldToNew[e[0]],
                    oldToNew[e[1]]
                );
            }

            // Construct an extendedEdgeMesh with
            // - all points on more than 2 edges : mixed feature points
            // - all edges : external edges

            extendedEdgeMesh eeMesh
            (
                newPoints,          // pts
                newEdges,           // eds
                0,                  // (point) concaveStart
                0,                  // (point) mixedStart
                nFeatures,          // (point) nonFeatureStart
                edges.size(),       // (edge) internalStart
                edges.size(),       // (edge) flatStart
                edges.size(),       // (edge) openStart
                edges.size(),       // (edge) multipleStart
                vectorField(0),     // normals
                List<extendedEdgeMesh::sideVolumeType>(0),// normalVolumeTypes
                vectorField(0),     // edgeDirections
                labelListList(0),   // normalDirections
                labelListList(0),   // edgeNormals
                labelListList(0),   // featurePointNormals
                labelListList(0),   // featurePointEdges
                identity(newEdges.size())   // regionEdges
            );

            // Info<< "Constructed extendedFeatureEdgeMesh " << featObj.name()
            //    << nl << incrIndent;
            // eeMesh.writeStats(Info);
            // Info<< decrIndent << endl;

            set(feati, new extendedFeatureEdgeMesh(featObj, eeMesh));
        }

        const extendedEdgeMesh& eMesh = operator[](feati);

        if (dict.found("levels"))
        {
            List<Tuple2<scalar, label>> distLevels(dict["levels"]);

            if (dict.size() < 1)
            {
                FatalErrorInFunction
                    << " : levels should be at least size 1" << endl
                    << "levels : "  << dict["levels"]
                    << exit(FatalError);
            }

            distances_[feati].setSize(distLevels.size());
            levels_[feati].setSize(distLevels.size());

            forAll(distLevels, j)
            {
                distances_[feati][j] = distLevels[j].first();
                levels_[feati][j] = distLevels[j].second();

                // Check in incremental order
                if (j > 0)
                {
                    if
                    (
                        (distances_[feati][j] <= distances_[feati][j-1])
                     || (levels_[feati][j] > levels_[feati][j-1])
                    )
                    {
                        FatalErrorInFunction
                            << " : Refinement should be specified in order"
                            << " of increasing distance"
                            << " (and decreasing refinement level)." << endl
                            << "Distance:" << distances_[feati][j]
                            << " refinementLevel:" << levels_[feati][j]
                            << exit(FatalError);
                    }
                }
            }
        }
        else
        {
            // Look up 'level' for single level
            levels_[feati] = labelList(1, dict.lookup<label>("level"));
            distances_[feati] = scalarField(1, 0.0);
        }

        Info<< "Refinement level according to distance to "
            << featFileName << " (" << eMesh.points().size() << " points, "
            << eMesh.edges().size() << " edges)." << endl;
        forAll(levels_[feati], j)
        {
            Info<< "    level " << levels_[feati][j]
                << " for all cells within " << distances_[feati][j]
                << " metre." << endl;
        }
    }
}


void Foam::refinementFeatures::buildTrees(const label feati)
{
    const extendedEdgeMesh& eMesh = operator[](feati);
    const pointField& points = eMesh.points();
    const edgeList& edges = eMesh.edges();

    // Calculate bb of all points
    treeBoundBox bb(points);

    // Slightly extended bb. Slightly off-centred just so on symmetric
    // geometry there are less face/edge aligned items.
    bb = bb.extend(1e-4);

    edgeTrees_.set
    (
        feati,
        new indexedOctree<treeDataEdge>
        (
            treeDataEdge
            (
                false,                  // do not cache bb
                edges,
                points,
                identity(edges.size())
            ),
            bb,     // overall search domain
            8,      // maxLevel
            10,     // leafsize
            3.0     // duplicity
        )
    );


    labelList featurePoints(identity(eMesh.nonFeatureStart()));

    pointTrees_.set
    (
        feati,
        new indexedOctree<treeDataPoint>
        (
            treeDataPoint(points, featurePoints),
            bb,     // overall search domain
            8,      // maxLevel
            10,     // leafsize
            3.0     // duplicity
        )
    );
}


// Find maximum level of a feature edge.
void Foam::refinementFeatures::findHigherLevel
(
    const pointField& pt,
    const label feati,
    labelList& maxLevel
) const
{
    const labelList& levels = levels_[feati];

    const scalarField& distances = distances_[feati];

    // Collect all those points that have a current maxLevel less than
    // (any of) the feature edge. Also collect the furthest distance allowable
    // to any feature edge with a higher level.

    pointField candidates(pt.size());
    labelList candidateMap(pt.size());
    scalarField candidateDistSqr(pt.size());
    label candidatei = 0;

    forAll(maxLevel, pointi)
    {
        forAllReverse(levels, leveli)
        {
            if (levels[leveli] > maxLevel[pointi])
            {
                candidates[candidatei] = pt[pointi];
                candidateMap[candidatei] = pointi;
                candidateDistSqr[candidatei] = sqr(distances[leveli]);
                candidatei++;
                break;
            }
        }
    }
    candidates.setSize(candidatei);
    candidateMap.setSize(candidatei);
    candidateDistSqr.setSize(candidatei);

    // Do the expensive nearest test only for the candidate points.
    const indexedOctree<treeDataEdge>& tree = edgeTrees_[feati];

    List<pointIndexHit> nearInfo(candidates.size());
    forAll(candidates, candidatei)
    {
        nearInfo[candidatei] = tree.findNearest
        (
            candidates[candidatei],
            candidateDistSqr[candidatei]
        );
    }

    // Update maxLevel
    forAll(nearInfo, candidatei)
    {
        if (nearInfo[candidatei].hit())
        {
            // Check which level it actually is in.
            label minDistI = findLower
            (
                distances,
                mag(nearInfo[candidatei].hitPoint()-candidates[candidatei])
            );

            label pointi = candidateMap[candidatei];

            // pt is in between feature[minDistI] and feature[minDistI+1]
            maxLevel[pointi] = levels[minDistI+1];
        }
    }
}


// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

const Foam::PtrList<Foam::indexedOctree<Foam::treeDataEdge>>&
Foam::refinementFeatures::regionEdgeTrees() const
{
    if (!regionEdgeTreesPtr_.valid())
    {
        regionEdgeTreesPtr_.reset
        (
            new PtrList<indexedOctree<treeDataEdge>>(size())
        );
        PtrList<indexedOctree<treeDataEdge>>& trees = regionEdgeTreesPtr_();

        forAll(*this, feati)
        {
            const extendedEdgeMesh& eMesh = operator[](feati);
            const pointField& points = eMesh.points();
            const edgeList& edges = eMesh.edges();

            // Calculate bb of all points
            treeBoundBox bb(points);

            // Slightly extended bb. Slightly off-centred just so on symmetric
            // geometry there are less face/edge aligned items.
            bb = bb.extend(1e-4);

            trees.set
            (
                feati,
                new indexedOctree<treeDataEdge>
                (
                    treeDataEdge
                    (
                        false,                  // do not cache bb
                        edges,
                        points,
                        eMesh.regionEdges()
                    ),
                    bb,     // overall search domain
                    8,      // maxLevel
                    10,     // leafsize
                    3.0     // duplicity
                )
            );
        }
    }
    return regionEdgeTreesPtr_();
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::refinementFeatures::refinementFeatures
(
    const objectRegistry& io,
    const PtrList<dictionary>& featDicts
)
:
    PtrList<extendedFeatureEdgeMesh>(featDicts.size()),
    distances_(featDicts.size()),
    levels_(featDicts.size()),
    edgeTrees_(featDicts.size()),
    pointTrees_(featDicts.size())
{
    // Read features
    read(io, featDicts);

    // Search engines
    forAll(*this, i)
    {
        buildTrees(i);
    }
}


//Foam::refinementFeatures::refinementFeatures
//(
//    const objectRegistry& io,
//    const PtrList<dictionary>& featDicts,
//    const scalar minCos
//)
//:
//    PtrList<extendedFeatureEdgeMesh>(featDicts.size()),
//    distances_(featDicts.size()),
//    levels_(featDicts.size()),
//    edgeTrees_(featDicts.size()),
//    pointTrees_(featDicts.size())
//{
//    // Read features
//    read(io, featDicts);
//
//    // Search engines
//    forAll(*this, i)
//    {
//        const edgeMesh& eMesh = operator[](i);
//        const pointField& points = eMesh.points();
//        const edgeList& edges = eMesh.edges();
//        const labelListList& pointEdges = eMesh.pointEdges();
//
//        DynamicList<label> featurePoints;
//        forAll(pointEdges, pointi)
//        {
//            const labelList& pEdges = pointEdges[pointi];
//            if (pEdges.size() > 2)
//            {
//                featurePoints.append(pointi);
//            }
//            else if (pEdges.size() == 2)
//            {
//                // Check the angle
//                const edge& e0 = edges[pEdges[0]];
//                const edge& e1 = edges[pEdges[1]];
//
//                const point& p = points[pointi];
//                const point& p0 = points[e0.otherVertex(pointi)];
//                const point& p1 = points[e1.otherVertex(pointi)];
//
//                vector v0 = p-p0;
//                scalar v0Mag = mag(v0);
//
//                vector v1 = p1-p;
//                scalar v1Mag = mag(v1);
//
//                if
//                (
//                    v0Mag > small
//                 && v1Mag > small
//                 && ((v0/v0Mag & v1/v1Mag) < minCos)
//                )
//                {
//                    featurePoints.append(pointi);
//                }
//            }
//        }
//
//        Info<< "Detected " << featurePoints.size()
//            << " featurePoints out of " << points.size()
//            << " points on feature " << i   // eMesh.name()
//            << " when using feature cos " << minCos << endl;
//
//        buildTrees(i, featurePoints);
//    }
//}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::refinementFeatures::findNearestEdge
(
    const pointField& samples,
    const scalarField& nearestDistSqr,
    labelList& nearFeature,
    List<pointIndexHit>& nearInfo,
    vectorField& nearNormal
) const
{
    nearFeature.setSize(samples.size());
    nearFeature = -1;
    nearInfo.setSize(samples.size());
    nearInfo = pointIndexHit();
    nearNormal.setSize(samples.size());
    nearNormal = Zero;

    forAll(edgeTrees_, feati)
    {
        const indexedOctree<treeDataEdge>& tree = edgeTrees_[feati];

        if (tree.shapes().size() > 0)
        {
            forAll(samples, samplei)
            {
                const point& sample = samples[samplei];

                scalar distSqr;
                if (nearInfo[samplei].hit())
                {
                    distSqr = magSqr(nearInfo[samplei].hitPoint()-sample);
                }
                else
                {
                    distSqr = nearestDistSqr[samplei];
                }

                pointIndexHit info = tree.findNearest(sample, distSqr);

                if (info.hit())
                {
                    nearFeature[samplei] = feati;
                    nearInfo[samplei] = pointIndexHit
                    (
                        info.hit(),
                        info.hitPoint(),
                        tree.shapes().edgeLabels()[info.index()]
                    );

                    const treeDataEdge& td = tree.shapes();
                    const edge& e = td.edges()[nearInfo[samplei].index()];
                    nearNormal[samplei] =  e.vec(td.points());
                    nearNormal[samplei] /= mag(nearNormal[samplei])+vSmall;
                }
            }
        }
    }
}


void Foam::refinementFeatures::findNearestRegionEdge
(
    const pointField& samples,
    const scalarField& nearestDistSqr,
    labelList& nearFeature,
    List<pointIndexHit>& nearInfo,
    vectorField& nearNormal
) const
{
    nearFeature.setSize(samples.size());
    nearFeature = -1;
    nearInfo.setSize(samples.size());
    nearInfo = pointIndexHit();
    nearNormal.setSize(samples.size());
    nearNormal = Zero;


    const PtrList<indexedOctree<treeDataEdge>>& regionTrees =
        regionEdgeTrees();

    forAll(regionTrees, feati)
    {
        const indexedOctree<treeDataEdge>& regionTree = regionTrees[feati];

        forAll(samples, samplei)
        {
            const point& sample = samples[samplei];

            scalar distSqr;
            if (nearInfo[samplei].hit())
            {
                distSqr = magSqr(nearInfo[samplei].hitPoint()-sample);
            }
            else
            {
                distSqr = nearestDistSqr[samplei];
            }

            // Find anything closer than current best
            pointIndexHit info = regionTree.findNearest(sample, distSqr);

            if (info.hit())
            {
                const treeDataEdge& td = regionTree.shapes();

                nearFeature[samplei] = feati;
                nearInfo[samplei] = pointIndexHit
                (
                    info.hit(),
                    info.hitPoint(),
                    regionTree.shapes().edgeLabels()[info.index()]
                );

                const edge& e = td.edges()[nearInfo[samplei].index()];
                nearNormal[samplei] =  e.vec(td.points());
                nearNormal[samplei] /= mag(nearNormal[samplei])+vSmall;
            }
        }
    }
}


//void Foam::refinementFeatures::findNearestPoint
//(
//    const pointField& samples,
//    const scalarField& nearestDistSqr,
//    labelList& nearFeature,
//    labelList& nearIndex
//) const
//{
//    nearFeature.setSize(samples.size());
//    nearFeature = -1;
//    nearIndex.setSize(samples.size());
//    nearIndex = -1;
//
//    forAll(pointTrees_, feati)
//    {
//        const indexedOctree<treeDataPoint>& tree = pointTrees_[feati];
//
//        if (tree.shapes().pointLabels().size() > 0)
//        {
//            forAll(samples, samplei)
//            {
//                const point& sample = samples[samplei];
//
//                scalar distSqr;
//                if (nearFeature[samplei] != -1)
//                {
//                    label nearFeatI = nearFeature[samplei];
//                    const indexedOctree<treeDataPoint>& nearTree =
//                        pointTrees_[nearFeatI];
//                    label featPointi =
//                        nearTree.shapes().pointLabels()[nearIndex[samplei]];
//                    const point& featPt =
//                        operator[](nearFeatI).points()[featPointi];
//                    distSqr = magSqr(featPt-sample);
//                }
//                else
//                {
//                    distSqr = nearestDistSqr[samplei];
//                }
//
//                pointIndexHit info = tree.findNearest(sample, distSqr);
//
//                if (info.hit())
//                {
//                    nearFeature[samplei] = feati;
//                    nearIndex[samplei] = info.index();
//                }
//            }
//        }
//    }
//}


void Foam::refinementFeatures::findNearestPoint
(
    const pointField& samples,
    const scalarField& nearestDistSqr,
    labelList& nearFeature,
    List<pointIndexHit>& nearInfo
) const
{
    nearFeature.setSize(samples.size());
    nearFeature = -1;
    nearInfo.setSize(samples.size());
    nearInfo = pointIndexHit();

    forAll(pointTrees_, feati)
    {
        const indexedOctree<treeDataPoint>& tree = pointTrees_[feati];

        if (tree.shapes().pointLabels().size() > 0)
        {
            forAll(samples, samplei)
            {
                const point& sample = samples[samplei];

                scalar distSqr;
                if (nearFeature[samplei] != -1)
                {
                    distSqr = magSqr(nearInfo[samplei].hitPoint()-sample);
                }
                else
                {
                    distSqr = nearestDistSqr[samplei];
                }

                pointIndexHit info = tree.findNearest(sample, distSqr);

                if (info.hit())
                {
                    nearFeature[samplei] = feati;
                    nearInfo[samplei] = pointIndexHit
                    (
                        info.hit(),
                        info.hitPoint(),
                        tree.shapes().pointLabels()[info.index()]
                    );
                }
            }
        }
    }
}


void Foam::refinementFeatures::findHigherLevel
(
    const pointField& pt,
    const labelList& ptLevel,
    labelList& maxLevel
) const
{
    // Maximum level of any feature edge. Start off with level of point.
    maxLevel = ptLevel;

    forAll(*this, feati)
    {
        findHigherLevel(pt, feati, maxLevel);
    }
}


Foam::scalar Foam::refinementFeatures::maxDistance() const
{
    scalar overallMax = -great;
    forAll(distances_, feati)
    {
        overallMax = max(overallMax, max(distances_[feati]));
    }
    return overallMax;
}


// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2022 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/>.

Class
    Foam::refinementFeatures

Description
    Encapsulates queries for features.

SourceFiles
    refinementFeatures.C

\*---------------------------------------------------------------------------*/

#ifndef refinementFeatures_H
#define refinementFeatures_H

#include "extendedFeatureEdgeMesh.H"
#include "indexedOctree.H"
#include "treeDataEdge.H"
#include "treeDataPoint.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
                     Class refinementFeatures Declaration
\*---------------------------------------------------------------------------*/

class refinementFeatures
:
    public PtrList<extendedFeatureEdgeMesh>
{
    // Private Data

        //- Per feature the list of ranges
        List<scalarField> distances_;

        //- Per feature per distance the refinement level
        labelListList levels_;

        //- Edge
        PtrList<indexedOctree<treeDataEdge>> edgeTrees_;

        //- Features points
        PtrList<indexedOctree<treeDataPoint>> pointTrees_;

        //- Region edge trees (demand driven)
        mutable autoPtr<PtrList<indexedOctree<treeDataEdge>>>
            regionEdgeTreesPtr_;


    // Private Member Functions

        //- Read set of feature edge meshes
        void read(const objectRegistry&, const PtrList<dictionary>&);

        //- Build edge tree and feature point tree
        void buildTrees(const label);

        //- Find feature level higher than ptLevel
        void findHigherLevel
        (
            const pointField& pt,
            const label feati,
            labelList& maxLevel
        ) const;


protected:

        const PtrList<indexedOctree<treeDataEdge>>& edgeTrees() const
        {
            return edgeTrees_;
        }

        const PtrList<indexedOctree<treeDataPoint>>& pointTrees() const
        {
            return pointTrees_;
        }


        const PtrList<indexedOctree<treeDataEdge>>& regionEdgeTrees() const;

public:

    // Constructors

        //- Construct from description
        refinementFeatures
        (
            const objectRegistry& io,
            const PtrList<dictionary>& featDicts
        );


    // Member Functions

        // Access

            //- Per featureEdgeMesh the list of level
            const labelListList& levels() const
            {
                return levels_;
            }

            //- Per featureEdgeMesh the list of ranges
            const List<scalarField>& distances() const
            {
                return distances_;
            }


        // Query

            //- Highest distance of all features
            scalar maxDistance() const;

            //- Find nearest point on nearest feature edge. Sets:
            //
            //    - nearFeature: index of feature mesh
            //    - nearInfo   : location on feature edge and edge index
            //        (note: not feature edge index but index into edges()
            //        directly)
            //    - nearNormal : local feature edge normal
            void findNearestEdge
            (
                const pointField& samples,
                const scalarField& nearestDistSqr,
                labelList& nearFeature,
                List<pointIndexHit>& nearInfo,
                vectorField& nearNormal
            ) const;

            //- Find nearest point on nearest region edge. Sets:
            //
            //    - nearFeature: index of feature mesh
            //    - nearInfo   : location on feature edge and edge index
            //        (note: not feature edge index but index into edges()
            //        directly)
            //    - nearNormal : local feature edge normal
            void findNearestRegionEdge
            (
                const pointField& samples,
                const scalarField& nearestDistSqr,
                labelList& nearFeature,
                List<pointIndexHit>& nearInfo,
                vectorField& nearNormal
            ) const;

            //- Find nearest feature point. Sets:
            //
            //    - nearFeature: index of feature mesh
            //    - nearInfo   : location on feature point and point index.
            //        (note: not index into shapes().pointLabels() but index
            //        into points() directly)
            void findNearestPoint
            (
                const pointField& samples,
                const scalarField& nearestDistSqr,
                labelList& nearFeature,
                List<pointIndexHit>& nearInfo
            ) const;

            //- Find feature level higher than ptLevel
            void findHigherLevel
            (
                const pointField& pt,
                const labelList& ptLevel,
                labelList& maxLevel
            ) const;

};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

diff --git a/etc/caseDicts/annotated/snappyHexMeshDict b/etc/caseDicts/annotated/snappyHexMeshDict
index 0407691..cb98025 100644
--- a/etc/caseDicts/annotated/snappyHexMeshDict
+++ b/etc/caseDicts/annotated/snappyHexMeshDict
@@ -50,9 +50,18 @@ geometry
         type triSurfaceMesh;
         file "sphere.stl"
 
-        // tolerance   1e-5;   // optional:non-default tolerance on intersections
-        // maxTreeDepth 10;    // optional:depth of octree. Decrease only in case
-                              // of memory limitations.
+        // Optional: non-default tolerance on intersections
+        // tolerance 1e-5;
+
+        // Optional: depth of octree. Decrease only in case of memory
+        // limitations.
+        // maxTreeDepth 10;
+
+        // Optional: scaling factor, e.g. unit conversion
+        // scale 1;
+
+        // Optional: discard triangles with low quality when getting normal
+        // minQuality -1;
 
         // Per region the patchname. If not provided will be <surface>_<region>.
         // Note: this name cannot be used to identity this region in any
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
index fea1f26..d65cb97 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -707,7 +707,7 @@ Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
     const labelList& cellLevel = meshCutter_.cellLevel();
     const pointField& cellCentres = mesh_.cellCentres();
 
-    // Detect if there are any distance shells
+    // Detect if there are any distance features
     if (features_.maxDistance() <= 0.0)
     {
         return 0;
@@ -730,7 +730,7 @@ Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
         }
     }
 
-    // Do test to see whether cells is inside/outside shell with higher level
+    // Do test to see whether cells are in/near features with higher level
     labelList maxLevel;
     features_.findHigherLevel(testCc, testLevels, maxLevel);
 
@@ -2086,7 +2086,7 @@ Foam::labelList Foam::meshRefinement::refineCandidates
 
 
 
-        // Cells pierced by feature lines
+        // Cells pierced by feature edges
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
         if (featureRefinement)
@@ -2104,12 +2104,12 @@ Foam::labelList Foam::meshRefinement::refineCandidates
                 << ": " << nFeatures << " cells."  << endl;
         }
 
-        // Inside distance-to-feature shells
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // Inside distance-to-feature
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
         if (featureDistanceRefinement)
         {
-            const label nShell = markInternalDistanceToFeatureRefinement
+            const label nCellsFeat = markInternalDistanceToFeatureRefinement
             (
                 nAllowRefine,
 
@@ -2117,7 +2117,7 @@ Foam::labelList Foam::meshRefinement::refineCandidates
                 nRefine
             );
             Info<< "Marked for refinement due to distance to explicit features "
-                ": " << nShell << " cells."  << endl;
+                ": " << nCellsFeat << " cells."  << endl;
         }
 
         // Inside refinement shells
diff --git a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
index 2acdbf1..60289d9 100644
--- a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
+++ b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -295,7 +295,7 @@ void Foam::refinementFeatures::buildTrees(const label feati)
 }
 
 
-// Find maximum level of a shell.
+// Find maximum level of a feature edge.
 void Foam::refinementFeatures::findHigherLevel
 (
     const pointField& pt,
@@ -308,8 +308,8 @@ void Foam::refinementFeatures::findHigherLevel
     const scalarField& distances = distances_[feati];
 
     // Collect all those points that have a current maxLevel less than
-    // (any of) the shell. Also collect the furthest distance allowable
-    // to any shell with a higher level.
+    // (any of) the feature edge. Also collect the furthest distance allowable
+    // to any feature edge with a higher level.
 
     pointField candidates(pt.size());
     labelList candidateMap(pt.size());
@@ -361,7 +361,7 @@ void Foam::refinementFeatures::findHigherLevel
 
             label pointi = candidateMap[candidatei];
 
-            // pt is in between shell[minDistI] and shell[minDistI+1]
+            // pt is in between feature[minDistI] and feature[minDistI+1]
             maxLevel[pointi] = levels[minDistI+1];
         }
     }
@@ -747,7 +747,7 @@ void Foam::refinementFeatures::findHigherLevel
     labelList& maxLevel
 ) const
 {
-    // Maximum level of any shell. Start off with level of point.
+    // Maximum level of any feature edge. Start off with level of point.
     maxLevel = ptLevel;
 
     forAll(*this, feati)
diff --git a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H
index a7c8286..1343004 100644
--- a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H
+++ b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,10 +55,10 @@ class refinementFeatures
 {
     // Private Data
 
-        //- Per shell the list of ranges
+        //- Per feature the list of ranges
         List<scalarField> distances_;
 
-        //- Per shell per distance the refinement level
+        //- Per feature per distance the refinement level
         labelListList levels_;
 
         //- Edge
@@ -80,7 +80,7 @@ class refinementFeatures
         //- Build edge tree and feature point tree
         void buildTrees(const label);
 
-        //- Find shell level higher than ptLevel
+        //- Find feature level higher than ptLevel
         void findHigherLevel
         (
             const pointField& pt,
@@ -184,7 +184,7 @@ public:
                 List<pointIndexHit>& nearInfo
             ) const;
 
-            //- Find shell level higher than ptLevel
+            //- Find feature level higher than ptLevel
             void findHigherLevel
             (
                 const pointField& pt,