/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 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 "cellCuts.H"
#include "polyMesh.H"
#include "Time.H"
#include "ListOps.H"
#include "cellLooper.H"
#include "refineCell.H"
#include "meshTools.H"
#include "geomCellLooper.H"
#include "OFstream.H"
#include "plane.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
defineTypeNameAndDebug(cellCuts, 0);
}


// * * * * * * * * * * * * * Private Static Functions  * * * * * * * * * * * //

Foam::label Foam::cellCuts::findPartIndex
(
    const labelList& elems,
    const label nElems,
    const label val
)
{
    for (label i = 0; i < nElems; i++)
    {
        if (elems[i] == val)
        {
            return i;
        }
    }
    return -1;
}


Foam::boolList Foam::cellCuts::expand
(
    const label size,
    const labelList& labels
)
{
    boolList result(size, false);

    forAll(labels, labelI)
    {
        result[labels[labelI]] = true;
    }
    return result;
}


Foam::scalarField Foam::cellCuts::expand
(
    const label size,
    const labelList& labels,
    const scalarField& weights
)
{
    scalarField result(size, -GREAT);

    forAll(labels, labelI)
    {
        result[labels[labelI]] = weights[labelI];
    }
    return result;
}


Foam::label Foam::cellCuts::firstUnique
(
    const labelList& lst,
    const Map<label>& map
)
{
    forAll(lst, i)
    {
        if (!map.found(lst[i]))
        {
            return i;
        }
    }
    return -1;
}


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

void Foam::cellCuts::writeUncutOBJ
(
    const fileName& dir,
    const label cellI
) const
{
    //- Cell edges
    OFstream cutsStream(dir / "cell_" + name(cellI) + ".obj");

    Pout<< "Writing cell for time " <<  mesh().time().timeName()
        << " to " << cutsStream.name() << nl;

    meshTools::writeOBJ
    (
        cutsStream,
        mesh().cells(),
        mesh().faces(),
        mesh().points(),
        labelList(1, cellI)
    );

    //- Loop cutting cell in two
    OFstream cutStream(dir / "cellCuts_" + name(cellI) + ".obj");

    Pout<< "Writing raw cuts on cell for time " <<  mesh().time().timeName()
        << " to " << cutStream.name() << nl;

    const labelList& cPoints = mesh().cellPoints()[cellI];

    forAll(cPoints, i)
    {
        label pointI = cPoints[i];
        if (pointIsCut_[pointI])
        {
            meshTools::writeOBJ(cutStream, mesh().points()[pointI]);
        }
    }

    const pointField& pts = mesh().points();

    const labelList& cEdges = mesh().cellEdges()[cellI];

    forAll(cEdges, i)
    {
        label edgeI = cEdges[i];

        if (edgeIsCut_[edgeI])
        {
            const edge& e = mesh().edges()[edgeI];

            const scalar w = edgeWeight_[edgeI];

            meshTools::writeOBJ(cutStream, w*pts[e[1]] + (1-w)*pts[e[0]]);
        }
    }
}


void Foam::cellCuts::writeOBJ
(
    const fileName& dir,
    const label cellI,
    const pointField& loopPoints,
    const labelList& anchors
) const
{
    //- Cell edges
    OFstream cutsStream(dir / "cell_" + name(cellI) + ".obj");

    Pout<< "Writing cell for time " <<  mesh().time().timeName()
        << " to " << cutsStream.name() << nl;

    meshTools::writeOBJ
    (
        cutsStream,
        mesh().cells(),
        mesh().faces(),
        mesh().points(),
        labelList(1, cellI)
    );


    //- Loop cutting cell in two
    OFstream loopStream(dir / "cellLoop_" + name(cellI) + ".obj");

    Pout<< "Writing loop for time " <<  mesh().time().timeName()
        << " to " << loopStream.name() << nl;

    label vertI = 0;

    writeOBJ(loopStream, loopPoints, vertI);


    //- Anchors for cell
    OFstream anchorStream(dir / "anchors_" + name(cellI) + ".obj");

    Pout<< "Writing anchors for time " <<  mesh().time().timeName()
        << " to " << anchorStream.name() << endl;

    forAll(anchors, i)
    {
        meshTools::writeOBJ(anchorStream, mesh().points()[anchors[i]]);
    }
}


Foam::label Foam::cellCuts::edgeEdgeToFace
(
    const label cellI,
    const label edgeA,
    const label edgeB
) const
{
    const labelList& cFaces = mesh().cells()[cellI];

    forAll(cFaces, cFaceI)
    {
        label faceI = cFaces[cFaceI];

        const labelList& fEdges = mesh().faceEdges()[faceI];

        if
        (
            findIndex(fEdges, edgeA) != -1
         && findIndex(fEdges, edgeB) != -1
        )
        {
           return faceI;
        }
    }

    // Coming here means the loop is illegal since the two edges
    // are not shared by a face. We just mark loop as invalid and continue.

    WarningIn
    (
        "Foam::cellCuts::edgeEdgeToFace"
        "(const label cellI, const label edgeA,"
        "const label edgeB) const"
    )   << "cellCuts : Cannot find face on cell "
        << cellI << " that has both edges " << edgeA << ' ' << edgeB << endl
        << "faces : " << cFaces << endl
        << "edgeA : " << mesh().edges()[edgeA] << endl
        << "edgeB : " << mesh().edges()[edgeB] << endl
        << "Marking the loop across this cell as invalid" << endl;

    return -1;
}


Foam::label Foam::cellCuts::edgeVertexToFace
(
    const label cellI,
    const label edgeI,
    const label vertI
) const
{
    const labelList& cFaces = mesh().cells()[cellI];

    forAll(cFaces, cFaceI)
    {
        label faceI = cFaces[cFaceI];

        const face& f = mesh().faces()[faceI];

        const labelList& fEdges = mesh().faceEdges()[faceI];

        if
        (
            findIndex(fEdges, edgeI) != -1
         && findIndex(f, vertI) != -1
        )
        {
           return faceI;
        }
    }

    WarningIn
    (
        "Foam::cellCuts::edgeVertexToFace"
        "(const label cellI, const label edgeI, "
        "const label vertI) const"
    )   << "cellCuts : Cannot find face on cell "
        << cellI << " that has both edge " << edgeI << " and vertex "
        << vertI << endl
        << "faces : " << cFaces << endl
        << "edge : " << mesh().edges()[edgeI] << endl
        << "Marking the loop across this cell as invalid" << endl;

    return -1;
}


Foam::label Foam::cellCuts::vertexVertexToFace
(
    const label cellI,
    const label vertA,
    const label vertB
) const
{
    const labelList& cFaces = mesh().cells()[cellI];

    forAll(cFaces, cFaceI)
    {
        label faceI = cFaces[cFaceI];

        const face& f = mesh().faces()[faceI];

        if
        (
            findIndex(f, vertA) != -1
         && findIndex(f, vertB) != -1
        )
        {
           return faceI;
        }
    }

    WarningIn("Foam::cellCuts::vertexVertexToFace")
        << "cellCuts : Cannot find face on cell "
        << cellI << " that has vertex " << vertA << " and vertex "
        << vertB << endl
        << "faces : " << cFaces << endl
        << "Marking the loop across this cell as invalid" << endl;

    return -1;
}


void Foam::cellCuts::calcFaceCuts() const
{
    if (faceCutsPtr_)
    {
        FatalErrorIn("cellCuts::calcFaceCuts()")
            << "faceCuts already calculated" << abort(FatalError);
    }

    const faceList& faces = mesh().faces();


    faceCutsPtr_ = new labelListList(mesh().nFaces());
    labelListList& faceCuts = *faceCutsPtr_;

    for (label faceI = 0; faceI < mesh().nFaces(); faceI++)
    {
        const face& f = faces[faceI];

        // Big enough storage (possibly all points and all edges cut). Shrink
        // later on.
        labelList& cuts = faceCuts[faceI];

        cuts.setSize(2*f.size());

        label cutI = 0;

        // Do point-edge-point walk over face and collect all cuts.
        // Problem is that we want to start from one of the endpoints of a
        // string of connected cuts; we don't want to start somewhere in the
        // middle.

        // Pass1: find first point cut not preceeded by a cut.
        label startFp = -1;

        forAll(f, fp)
        {
            label v0 = f[fp];

            if (pointIsCut_[v0])
            {
                label vMin1 = f[f.rcIndex(fp)];

                if
                (
                    !pointIsCut_[vMin1]
                 && !edgeIsCut_[findEdge(faceI, v0, vMin1)]
                )
                {
                    cuts[cutI++] = vertToEVert(v0);
                    startFp = f.fcIndex(fp);
                    break;
                }
            }
        }

        // Pass2: first edge cut not preceeded by point cut
        if (startFp == -1)
        {
            forAll(f, fp)
            {
                label fp1 = f.fcIndex(fp);

                label v0 = f[fp];
                label v1 = f[fp1];

                label edgeI = findEdge(faceI, v0, v1);

                if (edgeIsCut_[edgeI] && !pointIsCut_[v0])
                {
                    cuts[cutI++] = edgeToEVert(edgeI);
                    startFp = fp1;
                    break;
                }
            }
        }

        // Pass3: nothing found so far. Either face is not cut at all or all
        // vertices are cut. Start from 0.
        if (startFp == -1)
        {
            startFp = 0;
        }

        // Store all cuts starting from startFp;
        label fp = startFp;

        bool allVerticesCut = true;

        forAll(f, i)
        {
            label fp1 = f.fcIndex(fp);

            // Get the three items: current vertex, next vertex and edge
            // inbetween
            label v0 = f[fp];
            label v1 = f[fp1];
            label edgeI = findEdge(faceI, v0, v1);

            if (pointIsCut_[v0])
            {
                cuts[cutI++] = vertToEVert(v0);
            }
            else
            {
                // For check if all vertices have been cut (= illegal)
                allVerticesCut = false;
            }

            if (edgeIsCut_[edgeI])
            {
                cuts[cutI++] = edgeToEVert(edgeI);
            }

            fp = fp1;
        }


        if (allVerticesCut)
        {
            WarningIn("Foam::cellCuts::calcFaceCuts() const")
                << "Face " << faceI << " vertices " << f
                << " has all its vertices cut. Not cutting face." << endl;

            cutI = 0;
        }

        // Remove duplicate starting point
        if (cutI > 1 && cuts[cutI-1] == cuts[0])
        {
            cutI--;
        }
        cuts.setSize(cutI);
    }
}


Foam::label Foam::cellCuts::findEdge
(
    const label faceI,
    const label v0,
    const label v1
) const
{
    const edgeList& edges = mesh().edges();

    const labelList& fEdges = mesh().faceEdges()[faceI];

    forAll(fEdges, i)
    {
        const edge& e = edges[fEdges[i]];

        if
        (
            (e[0] == v0 && e[1] == v1)
         || (e[0] == v1 && e[1] == v0)
        )
        {
            return fEdges[i];
        }
    }

    return -1;
}


Foam::label Foam::cellCuts::loopFace
(
    const label cellI,
    const labelList& loop
) const
{
    const cell& cFaces = mesh().cells()[cellI];

    forAll(cFaces, cFaceI)
    {
        label faceI = cFaces[cFaceI];

        const labelList& fEdges = mesh().faceEdges()[faceI];
        const face& f = mesh().faces()[faceI];

        bool allOnFace = true;

        forAll(loop, i)
        {
            label cut = loop[i];

            if (isEdge(cut))
            {
                if (findIndex(fEdges, getEdge(cut)) == -1)
                {
                    // Edge not on face. Skip face.
                    allOnFace = false;
                    break;
                }
            }
            else
            {
                if (findIndex(f, getVertex(cut)) == -1)
                {
                    // Vertex not on face. Skip face.
                    allOnFace = false;
                    break;
                }
            }
        }

        if (allOnFace)
        {
            // Found face where all elements of loop are on the face.
            return faceI;
        }
    }
    return -1;
}


bool Foam::cellCuts::walkPoint
(
    const label cellI,
    const label startCut,

    const label exclude0,
    const label exclude1,

    const label otherCut,

    label& nVisited,
    labelList& visited
) const
{
    label vertI = getVertex(otherCut);

    const labelList& pFaces = mesh().pointFaces()[vertI];

    forAll(pFaces, pFaceI)
    {
        label otherFaceI = pFaces[pFaceI];

        if
        (
            otherFaceI != exclude0
         && otherFaceI != exclude1
         && meshTools::faceOnCell(mesh(), cellI, otherFaceI)
        )
        {
            label oldNVisited = nVisited;

            bool foundLoop =
                walkCell
                (
                    cellI,
                    startCut,
                    otherFaceI,
                    otherCut,
                    nVisited,
                    visited
                );

            if (foundLoop)
            {
                return true;
            }

            // No success. Restore state and continue
            nVisited = oldNVisited;
        }
    }
    return false;
}


bool Foam::cellCuts::crossEdge
(
    const label cellI,
    const label startCut,
    const label faceI,
    const label otherCut,

    label& nVisited,
    labelList& visited
) const
{
    // Cross edge to other face
    label edgeI = getEdge(otherCut);

    label otherFaceI = meshTools::otherFace(mesh(), cellI, faceI, edgeI);

    // Store old state
    label oldNVisited = nVisited;

    // Recurse to otherCut
    bool foundLoop =
        walkCell
        (
            cellI,
            startCut,
            otherFaceI,
            otherCut,
            nVisited,
            visited
        );

    if (foundLoop)
    {
        return true;
    }
    else
    {
        // No success. Restore state (i.e. backtrack)
        nVisited = oldNVisited;

        return false;
    }
}


bool Foam::cellCuts::addCut
(
    const label cellI,
    const label cut,
    label& nVisited,
    labelList& visited
) const
{
    // Check for duplicate cuts.
    if (findPartIndex(visited, nVisited, cut) != -1)
    {
        // Truncate (copy of) visited for ease of printing.
        labelList truncVisited(visited);
        truncVisited.setSize(nVisited);

        Pout<< "For cell " << cellI << " : trying to add duplicate cut " << cut;
        labelList cuts(1, cut);
        writeCuts(Pout, cuts, loopWeights(cuts));

        Pout<< " to path:";
        writeCuts(Pout, truncVisited, loopWeights(truncVisited));
        Pout<< endl;

        return false;
    }
    else
    {
        visited[nVisited++] = cut;

        return true;
    }
}


bool Foam::cellCuts::walkFace
(
    const label cellI,
    const label startCut,
    const label faceI,
    const label cut,

    label& lastCut,
    label& beforeLastCut,
    label& nVisited,
    labelList& visited
) const
{
    const labelList& fCuts = faceCuts()[faceI];

    if (fCuts.size() < 2)
    {
        return false;
    }

    // Easy case : two cuts.
    if (fCuts.size() == 2)
    {
        if (fCuts[0] == cut)
        {
            if (!addCut(cellI, cut, nVisited, visited))
            {
                return false;
            }

            beforeLastCut = cut;
            lastCut = fCuts[1];

            return true;
        }
        else
        {
            if (!addCut(cellI, cut, nVisited, visited))
            {
                return false;
            }

            beforeLastCut = cut;
            lastCut = fCuts[0];

            return true;
        }
    }

    // Harder case: more than 2 cuts on face.
    // Should be start or end of string of cuts. Store all but last cut.
    if (fCuts[0] == cut)
    {
        // Walk forward
        for (label i = 0; i < fCuts.size()-1; i++)
        {
            if (!addCut(cellI, fCuts[i], nVisited, visited))
            {
                return false;
            }
        }
        beforeLastCut = fCuts[fCuts.size()-2];
        lastCut = fCuts[fCuts.size()-1];
    }
    else if (fCuts[fCuts.size()-1] == cut)
    {
        for (label i = fCuts.size()-1; i >= 1; --i)
        {
            if (!addCut(cellI, fCuts[i], nVisited, visited))
            {
                return false;
            }
        }
        beforeLastCut = fCuts[1];
        lastCut = fCuts[0];
    }
    else
    {
        WarningIn("Foam::cellCuts::walkFace")
            << "In middle of cut. cell:" << cellI << " face:" << faceI
            << " cuts:" << fCuts << " current cut:" << cut << endl;

        return false;
    }

    return true;
}



bool Foam::cellCuts::walkCell
(
    const label cellI,
    const label startCut,
    const label faceI,
    const label cut,

    label& nVisited,
    labelList& visited
) const
{
    // Walk across face, storing cuts. Return last two cuts visited.
    label lastCut = -1;
    label beforeLastCut = -1;


    if (debug & 2)
    {
        Pout<< "For cell:" << cellI << " walked across face " << faceI
            << " from cut ";
        labelList cuts(1, cut);
        writeCuts(Pout, cuts, loopWeights(cuts));
        Pout<< endl;
    }

    bool validWalk = walkFace
    (
        cellI,
        startCut,
        faceI,
        cut,

        lastCut,
        beforeLastCut,
        nVisited,
        visited
    );

    if (!validWalk)
    {
        return false;
    }

    if (debug & 2)
    {
        Pout<< "    to last cut ";
        labelList cuts(1, lastCut);
        writeCuts(Pout, cuts, loopWeights(cuts));
        Pout<< endl;
    }

    // Check if starting point reached.
    if (lastCut == startCut)
    {
        if (nVisited >= 3)
        {
            if (debug & 2)
            {
                // Truncate (copy of) visited for ease of printing.
                labelList truncVisited(visited);
                truncVisited.setSize(nVisited);

                Pout<< "For cell " << cellI << " : found closed path:";
                writeCuts(Pout, truncVisited, loopWeights(truncVisited));
                Pout<< " closed by " << lastCut << endl;
            }

            return true;
        }
        else
        {
            // Loop but too short. Probably folded back on itself.
            return false;
        }
    }


    // Check last cut and one before that to determine type.

    // From beforeLastCut to lastCut:
    //  - from edge to edge
    //      (- always not along existing edge)
    //  - from edge to vertex
    //      - not along existing edge
    //      - along existing edge: not allowed?
    //  - from vertex to edge
    //      - not along existing edge
    //      - along existing edge. Not allowed. See above.
    //  - from vertex to vertex
    //      - not along existing edge
    //      - along existing edge

    if (isEdge(beforeLastCut))
    {
        if (isEdge(lastCut))
        {
            // beforeLastCut=edge, lastCut=edge.

            // Cross lastCut (=edge) into face not faceI
            return crossEdge
            (
                cellI,
                startCut,
                faceI,
                lastCut,
                nVisited,
                visited
            );
        }
        else
        {
            // beforeLastCut=edge, lastCut=vertex.

            // Go from lastCut to all connected faces (except faceI)
            return walkPoint
            (
                cellI,
                startCut,
                faceI,          // exclude0: don't cross back on current face
                -1,             // exclude1
                lastCut,
                nVisited,
                visited
            );
        }
    }
    else
    {
        if (isEdge(lastCut))
        {
            // beforeLastCut=vertex, lastCut=edge.
            return crossEdge
            (
                cellI,
                startCut,
                faceI,
                lastCut,
                nVisited,
                visited
            );
        }
        else
        {
            // beforeLastCut=vertex, lastCut=vertex. Check if along existing
            // edge.
            label edgeI =
                findEdge
                (
                    faceI,
                    getVertex(beforeLastCut),
                    getVertex(lastCut)
                );

            if (edgeI != -1)
            {
                // Cut along existing edge. So is in fact on two faces.
                // Get faces on both sides of the edge to make
                // sure we dont fold back on to those.

                label f0, f1;
                meshTools::getEdgeFaces(mesh(), cellI, edgeI, f0, f1);

                return walkPoint
                (
                    cellI,
                    startCut,
                    f0,
                    f1,
                    lastCut,
                    nVisited,
                    visited
                );

            }
            else
            {
                // Cross cut across face.
                return walkPoint
                (
                    cellI,
                    startCut,
                    faceI,      // face to exclude
                    -1,         // face to exclude
                    lastCut,
                    nVisited,
                    visited
                );
            }
        }
    }
}


void Foam::cellCuts::calcCellLoops(const labelList& cutCells)
{
    // Determine for every cut cell the loop (= face) it is cut by. Done by
    // starting from a cut edge or cut vertex and walking across faces, from
    // cut to cut, until starting cut hit.
    // If multiple loops are possible across a cell circumference takes the
    // first one found.

    // Calculate cuts per face.
    const labelListList& allFaceCuts = faceCuts();

    // Per cell the number of faces with valid cuts. Is used as quick
    // rejection to see if cell can be cut.
    labelList nCutFaces(mesh().nCells(), 0);

    forAll(allFaceCuts, faceI)
    {
        const labelList& fCuts = allFaceCuts[faceI];

        if (fCuts.size() == mesh().faces()[faceI].size())
        {
            // Too many cuts on face. WalkCell would get very upset so disable.
            nCutFaces[mesh().faceOwner()[faceI]] = labelMin;

            if (mesh().isInternalFace(faceI))
            {
                nCutFaces[mesh().faceNeighbour()[faceI]] = labelMin;
            }
        }
        else if (fCuts.size() >= 2)
        {
            // Could be valid cut. Update count for owner and neighbour.
            nCutFaces[mesh().faceOwner()[faceI]]++;

            if (mesh().isInternalFace(faceI))
            {
                nCutFaces[mesh().faceNeighbour()[faceI]]++;
            }
        }
    }


    // Stack of visited cuts (nVisited used as stack pointer)
    // Size big enough.
    labelList visited(mesh().nPoints());

    forAll(cutCells, i)
    {
        label cellI = cutCells[i];

        bool validLoop = false;

        // Quick rejection: has enough faces that are cut?
        if (nCutFaces[cellI] >= 3)
        {
            const labelList& cFaces = mesh().cells()[cellI];

            if (debug & 2)
            {
                Pout<< "cell:" << cellI << " cut faces:" << endl;
                forAll(cFaces, i)
                {
                    label faceI = cFaces[i];
                    const labelList& fCuts = allFaceCuts[faceI];

                    Pout<< "    face:" << faceI << " cuts:";
                    writeCuts(Pout, fCuts, loopWeights(fCuts));
                    Pout<< endl;
                }
            }

            label nVisited = 0;

            // Determine the first cut face to start walking from.
            forAll(cFaces, cFaceI)
            {
                label faceI = cFaces[cFaceI];

                const labelList& fCuts = allFaceCuts[faceI];

                // Take first or last cut of multiple on face.
                // Note that in calcFaceCuts
                // we have already made sure this is the start or end of a
                // string of cuts.
                if (fCuts.size() >= 2)
                {
                    // Try walking from start of fCuts.
                    nVisited = 0;

                    if (debug & 2)
                    {
                        Pout<< "cell:" << cellI
                            << " start walk at face:" << faceI
                            << " cut:";
                        labelList cuts(1, fCuts[0]);
                        writeCuts(Pout, cuts, loopWeights(cuts));
                        Pout<< endl;
                    }

                    validLoop =
                        walkCell
                        (
                            cellI,
                            fCuts[0],
                            faceI,
                            fCuts[0],

                            nVisited,
                            visited
                        );

                    if (validLoop)
                    {
                        break;
                    }

                    // No need to try and walk from end of cuts since
                    // all paths already tried by walkCell.
                }
            }

            if (validLoop)
            {
                // Copy nVisited elements out of visited (since visited is
                // never truncated for efficiency reasons)

                labelList& loop = cellLoops_[cellI];

                loop.setSize(nVisited);

                for (label i = 0; i < nVisited; i++)
                {
                    loop[i] = visited[i];
                }
            }
            else
            {
                // Invalid loop. Leave cellLoops_[cellI] zero size which
                // flags this.
                Pout<< "calcCellLoops(const labelList&) : did not find valid"
                    << " loop for cell " << cellI << endl;
                // Dump cell and cuts on cell.
                writeUncutOBJ(".", cellI);

                cellLoops_[cellI].setSize(0);
            }
        }
        else
        {
            //Pout<< "calcCellLoops(const labelList&) : did not find valid"
            //    << " loop for cell " << cellI << " since not enough cut faces"
            //    << endl;
            cellLoops_[cellI].setSize(0);
        }
    }
}


void Foam::cellCuts::walkEdges
(
    const label cellI,
    const label pointI,
    const label status,

    Map<label>& edgeStatus,
    Map<label>& pointStatus
) const
{
    if (pointStatus.insert(pointI, status))
    {
        // First visit to pointI

        const labelList& pEdges = mesh().pointEdges()[pointI];

        forAll(pEdges, pEdgeI)
        {
            label edgeI = pEdges[pEdgeI];

            if
            (
                meshTools::edgeOnCell(mesh(), cellI, edgeI)
             && edgeStatus.insert(edgeI, status)
            )
            {
                // First visit to edgeI so recurse.

                label v2 = mesh().edges()[edgeI].otherVertex(pointI);

                walkEdges(cellI, v2, status, edgeStatus, pointStatus);
            }
        }
    }
}


Foam::labelList Foam::cellCuts::nonAnchorPoints
(
    const labelList& cellPoints,
    const labelList& anchorPoints,
    const labelList& loop
) const
{
    labelList newElems(cellPoints.size());
    label newElemI = 0;

    forAll(cellPoints, i)
    {
        label pointI = cellPoints[i];

        if
        (
            findIndex(anchorPoints, pointI) == -1
         && findIndex(loop, vertToEVert(pointI)) == -1
        )
        {
            newElems[newElemI++] = pointI;
        }
    }

    newElems.setSize(newElemI);

    return newElems;
}


bool Foam::cellCuts::loopAnchorConsistent
(
    const label cellI,
    const pointField& loopPts,
    const labelList& anchorPoints
) const
{
    // Create identity face for ease of calculation of normal etc.
    face f(identity(loopPts.size()));

    vector normal = f.normal(loopPts);
    point ctr = f.centre(loopPts);


    // Get average position of anchor points.
    vector avg(vector::zero);

    forAll(anchorPoints, ptI)
    {
        avg += mesh().points()[anchorPoints[ptI]];
    }
    avg /= anchorPoints.size();


    if (((avg - ctr) & normal) > 0)
    {
        return true;
    }
    else
    {
        return false;
    }
}


bool Foam::cellCuts::calcAnchors
(
    const label cellI,
    const labelList& loop,
    const pointField& loopPts,

    labelList& anchorPoints
) const
{
    const edgeList& edges = mesh().edges();

    const labelList& cPoints = mesh().cellPoints()[cellI];
    const labelList& cEdges = mesh().cellEdges()[cellI];
    const cell& cFaces = mesh().cells()[cellI];

    // Points on loop

    // Status of point:
    // 0 - on loop
    // 1 - point set 1
    // 2 - point set 2
    Map<label> pointStatus(2*cPoints.size());
    Map<label> edgeStatus(2*cEdges.size());

    // Mark loop vertices
    forAll(loop, i)
    {
        label cut = loop[i];

        if (isEdge(cut))
        {
            edgeStatus.insert(getEdge(cut), 0);
        }
        else
        {
            pointStatus.insert(getVertex(cut), 0);
        }
    }
    // Since edges between two cut vertices have not been marked, mark them
    // explicitly
    forAll(cEdges, i)
    {
        label edgeI = cEdges[i];
        const edge& e = edges[edgeI];

        if (pointStatus.found(e[0]) && pointStatus.found(e[1]))
        {
            edgeStatus.insert(edgeI, 0);
        }
    }


    // Find uncut starting vertex
    label uncutIndex = firstUnique(cPoints, pointStatus);

    if (uncutIndex == -1)
    {
        WarningIn("Foam::cellCuts::calcAnchors")
            << "Invalid loop " << loop << " for cell " << cellI << endl
            << "Can not find point on cell which is not cut by loop."
            << endl;

        writeOBJ(".", cellI, loopPts, labelList(0));

        return false;
    }

    // Walk unset vertices and edges and mark with 1 in pointStatus, edgeStatus
    walkEdges(cellI, cPoints[uncutIndex], 1, edgeStatus, pointStatus);

    // Find new uncut starting vertex
    uncutIndex = firstUnique(cPoints, pointStatus);

    if (uncutIndex == -1)
    {
        // All vertices either in loop or in anchor. So split is along single
        // face.
        WarningIn("Foam::cellCuts::calcAnchors")
            << "Invalid loop " << loop << " for cell " << cellI << endl
            << "All vertices of cell are either in loop or in anchor set"
            << endl;

        writeOBJ(".", cellI, loopPts, labelList(0));

        return false;
    }

    // Walk unset vertices and edges and mark with 2. These are the
    // pointset 2.
    walkEdges(cellI, cPoints[uncutIndex], 2, edgeStatus, pointStatus);

    // Collect both sets in lists.
    DynamicList<label> connectedPoints(cPoints.size());
    DynamicList<label> otherPoints(cPoints.size());

    forAllConstIter(Map<label>, pointStatus, iter)
    {
        if (iter() == 1)
        {
            connectedPoints.append(iter.key());
        }
        else if (iter() == 2)
        {
            otherPoints.append(iter.key());
        }
    }
    connectedPoints.shrink();
    otherPoints.shrink();

    // Check that all points have been used.
    uncutIndex = firstUnique(cPoints, pointStatus);

    if (uncutIndex != -1)
    {
        WarningIn("Foam::cellCuts::calcAnchors")
            << "Invalid loop " << loop << " for cell " << cellI
            << " since it splits the cell into more than two cells" << endl;

        writeOBJ(".", cellI, loopPts, connectedPoints);

        return false;
    }


    // Check that both parts (connectedPoints, otherPoints) have enough faces.
    labelHashSet connectedFaces(2*cFaces.size());
    labelHashSet otherFaces(2*cFaces.size());

    forAllConstIter(Map<label>, pointStatus, iter)
    {
        label pointI = iter.key();

        const labelList& pFaces = mesh().pointFaces()[pointI];

        if (iter() == 1)
        {
            forAll(pFaces, pFaceI)
            {
                if (meshTools::faceOnCell(mesh(), cellI, pFaces[pFaceI]))
                {
                    connectedFaces.insert(pFaces[pFaceI]);
                }
            }
        }
        else if (iter() == 2)
        {
            forAll(pFaces, pFaceI)
            {
                if (meshTools::faceOnCell(mesh(), cellI, pFaces[pFaceI]))
                {
                    otherFaces.insert(pFaces[pFaceI]);
                }
            }
        }
    }

    if (connectedFaces.size() < 3)
    {
        WarningIn("Foam::cellCuts::calcAnchors")
            << "Invalid loop " << loop << " for cell " << cellI
            << " since would have too few faces on one side." << nl
            << "All faces:" << cFaces << endl;

        writeOBJ(".", cellI, loopPts, connectedPoints);

        return false;
    }

    if (otherFaces.size() < 3)
    {
        WarningIn("Foam::cellCuts::calcAnchors")
            << "Invalid loop " << loop << " for cell " << cellI
            << " since would have too few faces on one side." << nl
            << "All faces:" << cFaces << endl;

        writeOBJ(".", cellI, loopPts, otherPoints);

        return false;
    }


    // Check that faces are split into two regions and not more.
    // When walking across the face the only transition of pointStatus is
    // from set1 to loop to set2 or back. Not allowed is from set1 to loop to
    // set1.
    {
        forAll(cFaces, i)
        {
            label faceI = cFaces[i];

            const face& f = mesh().faces()[faceI];

            bool hasSet1 = false;
            bool hasSet2 = false;

            label prevStat = pointStatus[f[0]];

            forAll(f, fp)
            {
                label v0 = f[fp];
                label pStat = pointStatus[v0];

                if (pStat == prevStat)
                {
                }
                else if (pStat == 0)
                {
                    // Loop.
                }
                else if (pStat == 1)
                {
                    if (hasSet1)
                    {
                        // Second occurence of set1.
                        WarningIn("Foam::cellCuts::calcAnchors")
                            << "Invalid loop " << loop << " for cell " << cellI
                            << " since face " << f << " would be split into"
                            << " more than two faces" << endl;

                        writeOBJ(".", cellI, loopPts, otherPoints);

                        return false;
                    }

                    hasSet1 = true;
                }
                else if (pStat == 2)
                {
                    if (hasSet2)
                    {
                        // Second occurence of set1.
                        WarningIn("Foam::cellCuts::calcAnchors")
                            << "Invalid loop " << loop << " for cell " << cellI
                            << " since face " << f << " would be split into"
                            << " more than two faces" << endl;

                        writeOBJ(".", cellI, loopPts, otherPoints);

                        return false;
                    }

                    hasSet2 = true;
                }
                else
                {
                    FatalErrorIn("Foam::cellCuts::calcAnchors")
                        << abort(FatalError);
                }

                prevStat = pStat;


                label v1 = f.nextLabel(fp);
                label edgeI = findEdge(faceI, v0, v1);

                label eStat = edgeStatus[edgeI];

                if (eStat == prevStat)
                {
                }
                else if (eStat == 0)
                {
                    // Loop.
                }
                else if (eStat == 1)
                {
                    if (hasSet1)
                    {
                        // Second occurence of set1.
                        WarningIn("Foam::cellCuts::calcAnchors")
                            << "Invalid loop " << loop << " for cell " << cellI
                            << " since face " << f << " would be split into"
                            << " more than two faces" << endl;

                        writeOBJ(".", cellI, loopPts, otherPoints);

                        return false;
                    }

                    hasSet1 = true;
                }
                else if (eStat == 2)
                {
                    if (hasSet2)
                    {
                        // Second occurence of set1.
                        WarningIn("Foam::cellCuts::calcAnchors")
                            << "Invalid loop " << loop << " for cell " << cellI
                            << " since face " << f << " would be split into"
                            << " more than two faces" << endl;

                        writeOBJ(".", cellI, loopPts, otherPoints);

                        return false;
                    }

                    hasSet2 = true;
                }
                prevStat = eStat;
            }
        }
    }




    // Check which one of point sets to use.
    bool loopOk = loopAnchorConsistent(cellI, loopPts, connectedPoints);

    //if (debug)
    {
        // Additional check: are non-anchor points on other side?
        bool otherLoopOk = loopAnchorConsistent(cellI, loopPts, otherPoints);

        if (loopOk == otherLoopOk)
        {
            // Both sets of points are supposedly on the same side as the
            // loop normal. Oops.

            WarningIn("Foam::cellCuts::calcAnchors")
                << " For cell:" << cellI
                << " achorpoints and nonanchorpoints are geometrically"
                << " on same side!" << endl
                << "cellPoints:" << cPoints << endl
                << "loop:" << loop << endl
                << "anchors:" << connectedPoints << endl
                << "otherPoints:" << otherPoints << endl;

            writeOBJ(".", cellI, loopPts, connectedPoints);
        }
    }

    if (loopOk)
    {
        // connectedPoints on 'outside' of loop so these are anchor points
        anchorPoints.transfer(connectedPoints);
        connectedPoints.clear();
    }
    else
    {
        anchorPoints.transfer(otherPoints);
        otherPoints.clear();
    }
    return true;
}


Foam::pointField Foam::cellCuts::loopPoints
(
    const labelList& loop,
    const scalarField& loopWeights
) const
{
    pointField loopPts(loop.size());

    forAll(loop, fp)
    {
        loopPts[fp] = coord(loop[fp], loopWeights[fp]);
    }
    return loopPts;
}


Foam::scalarField Foam::cellCuts::loopWeights(const labelList& loop) const
{
    scalarField weights(loop.size());

    forAll(loop, fp)
    {
        label cut = loop[fp];

        if (isEdge(cut))
        {
            label edgeI = getEdge(cut);

            weights[fp] = edgeWeight_[edgeI];
        }
        else
        {
            weights[fp] = -GREAT;
        }
    }
    return weights;
}


bool Foam::cellCuts::validEdgeLoop
(
    const labelList& loop,
    const scalarField& loopWeights
) const
{
    forAll(loop, fp)
    {
        label cut = loop[fp];

        if (isEdge(cut))
        {
            label edgeI = getEdge(cut);

            // Check: cut compatible only if can be snapped to existing one.
            if (edgeIsCut_[edgeI])
            {
                scalar edgeLen = mesh().edges()[edgeI].mag(mesh().points());

                if
                (
                    mag(loopWeights[fp] - edgeWeight_[edgeI])
                  > geomCellLooper::snapTol()*edgeLen
                )
                {
                    // Positions differ too much->would create two cuts.
                    return false;
                }
            }
        }
    }
    return true;
}


Foam::label Foam::cellCuts::countFaceCuts
(
    const label faceI,
    const labelList& loop
) const
{
    // Includes cuts through vertices and through edges.
    // Assumes that if edge is cut both in edgeIsCut and in loop that the
    // position of the cut is the same.

    label nCuts = 0;

    // Count cut vertices
    const face& f = mesh().faces()[faceI];

    forAll(f, fp)
    {
        label vertI = f[fp];

        // Vertex already cut or mentioned in current loop.
        if
        (
            pointIsCut_[vertI]
         || (findIndex(loop, vertToEVert(vertI)) != -1)
        )
        {
            nCuts++;
        }
    }

    // Count cut edges.
    const labelList& fEdges = mesh().faceEdges()[faceI];

    forAll(fEdges, fEdgeI)
    {
        label edgeI = fEdges[fEdgeI];

        // Edge already cut or mentioned in current loop.
        if
        (
            edgeIsCut_[edgeI]
         || (findIndex(loop, edgeToEVert(edgeI)) != -1)
        )
        {
            nCuts++;
        }
    }

    return nCuts;
}


bool Foam::cellCuts::conservativeValidLoop
(
    const label cellI,
    const labelList& loop
) const
{

    if (loop.size() < 2)
    {
        return false;
    }

    forAll(loop, cutI)
    {
        if (isEdge(loop[cutI]))
        {
            label edgeI = getEdge(loop[cutI]);

            if (edgeIsCut_[edgeI])
            {
                // edge compatibility already checked.
            }
            else
            {
                // Quick rejection: vertices of edge itself cannot be cut.
                const edge& e = mesh().edges()[edgeI];

                if (pointIsCut_[e.start()] || pointIsCut_[e.end()])
                {
                    return false;
                }


                // Check faces using this edge
                const labelList& eFaces = mesh().edgeFaces()[edgeI];

                forAll(eFaces, eFaceI)
                {
                    label nCuts = countFaceCuts(eFaces[eFaceI], loop);

                    if (nCuts > 2)
                    {
                        return false;
                    }
                }
            }
        }
        else
        {
            // Vertex cut

            label vertI = getVertex(loop[cutI]);

            if (!pointIsCut_[vertI])
            {
                // New cut through vertex.

                // Check edges using vertex.
                const labelList& pEdges = mesh().pointEdges()[vertI];

                forAll(pEdges, pEdgeI)
                {
                    label edgeI = pEdges[pEdgeI];

                    if (edgeIsCut_[edgeI])
                    {
                        return false;
                    }
                }

                // Check faces using vertex.
                const labelList& pFaces = mesh().pointFaces()[vertI];

                forAll(pFaces, pFaceI)
                {
                    label nCuts = countFaceCuts(pFaces[pFaceI], loop);

                    if (nCuts > 2)
                    {
                        return false;
                    }
                }
            }
        }
    }

    // All ok.
    return true;
}


bool Foam::cellCuts::validLoop
(
    const label cellI,
    const labelList& loop,
    const scalarField& loopWeights,

    Map<edge>& newFaceSplitCut,
    labelList& anchorPoints
) const
{
    // Determine compatibility of loop with existing cut pattern. Does not use
    // derived cut-addressing (faceCuts), only pointIsCut, edgeIsCut.
    // Adds any cross-cuts found to newFaceSplitCut and sets cell points on
    // one side of the loop in anchorPoints.

    if (loop.size() < 2)
    {
        return false;
    }

    if (debug & 4)
    {
        // Allow as fallback the 'old' loop checking where only a single
        // cut per face is allowed.
        if (!conservativeValidLoop(cellI, loop))
        {
            Info << "Invalid conservative loop: " << loop << endl;
            return  false;
        }
    }

    forAll(loop, fp)
    {
        label cut = loop[fp];
        label nextCut = loop[(fp+1) % loop.size()];

        // Label (if any) of face cut (so cut not along existing edge)
        label meshFaceI = -1;

        if (isEdge(cut))
        {
            label edgeI = getEdge(cut);

            // Look one cut ahead to find if it is along existing edge.

            if (isEdge(nextCut))
            {
                // From edge to edge -> cross cut
                label nextEdgeI = getEdge(nextCut);

                // Find face and mark as to be split.
                meshFaceI = edgeEdgeToFace(cellI, edgeI, nextEdgeI);

                if (meshFaceI == -1)
                {
                    // Can't find face using both cut edges.
                    return false;
                }
            }
            else
            {
                // From edge to vertex -> cross cut only if no existing edge.

                label nextVertI = getVertex(nextCut);
                const edge& e = mesh().edges()[edgeI];

                if (e.start() != nextVertI && e.end() != nextVertI)
                {
                    // New edge. Find face and mark as to be split.
                    meshFaceI = edgeVertexToFace(cellI, edgeI, nextVertI);

                    if (meshFaceI == -1)
                    {
                        // Can't find face. Ilegal.
                        return false;
                    }
                }
            }
        }
        else
        {
            // Cut is vertex.
            label vertI = getVertex(cut);

            if (isEdge(nextCut))
            {
                // From vertex to edge -> cross cut only if no existing edge.
                label nextEdgeI = getEdge(nextCut);

                const edge& nextE = mesh().edges()[nextEdgeI];

                if (nextE.start() != vertI && nextE.end() != vertI)
                {
                    // Should be cross cut. Find face.
                    meshFaceI = edgeVertexToFace(cellI, nextEdgeI, vertI);

                    if (meshFaceI == -1)
                    {
                        return false;
                    }
                }
            }
            else
            {
                // From vertex to vertex -> cross cut only if no existing edge.
                label nextVertI = getVertex(nextCut);

                if (meshTools::findEdge(mesh(), vertI, nextVertI) == -1)
                {
                    // New edge. Find face.
                    meshFaceI = vertexVertexToFace(cellI, vertI, nextVertI);

                    if (meshFaceI == -1)
                    {
                        return false;
                    }
                }
            }
        }

        if (meshFaceI != -1)
        {
            // meshFaceI is cut across along cut-nextCut (so not along existing
            // edge). Check if this is compatible with existing pattern.
            edge cutEdge(cut, nextCut);

            Map<edge>::const_iterator iter = faceSplitCut_.find(meshFaceI);

            if (iter == faceSplitCut_.end())
            {
                // Face not yet cut so insert.
                newFaceSplitCut.insert(meshFaceI, cutEdge);
            }
            else
            {
                // Face already cut. Ok if same edge.
                if (iter() != cutEdge)
                {
                    return false;
                }
            }
        }
    }

    // Is there a face on which all cuts are?
    label faceContainingLoop = loopFace(cellI, loop);

    if (faceContainingLoop != -1)
    {
        WarningIn("Foam::cellCuts::validLoop")
            << "Found loop on cell " << cellI << " with all points"
            << " on face " << faceContainingLoop << endl;

        //writeOBJ(".", cellI, loopPoints(loop, loopWeights), labelList(0));

        return false;
    }

    // Calculate anchor points
    // Final success is determined by whether anchor points can be determined.
    return calcAnchors
    (
        cellI,
        loop,
        loopPoints(loop, loopWeights),
        anchorPoints
    );
}


void Foam::cellCuts::setFromCellLoops()
{
    // 'Uncut' edges/vertices that are not used in loops.
    pointIsCut_ = false;

    edgeIsCut_ = false;

    faceSplitCut_.clear();

    forAll(cellLoops_, cellI)
    {
        const labelList& loop = cellLoops_[cellI];

        if (loop.size())
        {
            // Storage for cross-face cuts
            Map<edge> faceSplitCuts(loop.size());
            // Storage for points on one side of cell.
            labelList anchorPoints;

            if
            (
               !validLoop
                (
                    cellI,
                    loop,
                    loopWeights(loop),
                    faceSplitCuts,
                    anchorPoints
                )
            )
            {
                //writeOBJ(".", cellI, loopPoints(cellI), anchorPoints);

                //FatalErrorIn("cellCuts::setFromCellLoops()")
                WarningIn("cellCuts::setFromCellLoops")
                    << "Illegal loop " << loop
                    << " when recreating cut-addressing"
                    << " from existing cellLoops for cell " << cellI
                    << endl;
                    //<< abort(FatalError);

                cellLoops_[cellI].setSize(0);
                cellAnchorPoints_[cellI].setSize(0);
            }
            else
            {
                // Copy anchor points.
                cellAnchorPoints_[cellI].transfer(anchorPoints);

                // Copy faceSplitCuts into overall faceSplit info.
                forAllConstIter(Map<edge>, faceSplitCuts, iter)
                {
                    faceSplitCut_.insert(iter.key(), iter());
                }

                // Update edgeIsCut, pointIsCut information
                forAll(loop, cutI)
                {
                    label cut = loop[cutI];

                    if (isEdge(cut))
                    {
                        edgeIsCut_[getEdge(cut)] = true;
                    }
                    else
                    {
                        pointIsCut_[getVertex(cut)] = true;
                    }
                }
            }
        }
    }

    // Reset edge weights
    forAll(edgeIsCut_, edgeI)
    {
        if (!edgeIsCut_[edgeI])
        {
            // Weight not used. Set to illegal value to make any use fall over.
            edgeWeight_[edgeI] = -GREAT;
        }
    }
}


bool Foam::cellCuts::setFromCellLoop
(
    const label cellI,
    const labelList& loop,
    const scalarField& loopWeights
)
{
    // Update basic cut information from single cellLoop. Returns true if loop
    // was valid.

    // Dump loop for debugging.
    if (debug)
    {
        OFstream str("last_cell.obj");

        str<< "# edges of cell " << cellI << nl;

        meshTools::writeOBJ
        (
            str,
            mesh().cells(),
            mesh().faces(),
            mesh().points(),
            labelList(1, cellI)
        );


        OFstream loopStr("last_loop.obj");

        loopStr<< "# looppoints for cell " << cellI << nl;

        pointField pointsOfLoop = loopPoints(loop, loopWeights);

        forAll(pointsOfLoop, i)
        {
            meshTools::writeOBJ(loopStr, pointsOfLoop[i]);
        }

        str << 'l';

        forAll(pointsOfLoop, i)
        {
            str << ' ' << i + 1;
        }
        str << ' ' << 1 << nl;
    }

    bool okLoop = false;

    if (validEdgeLoop(loop, loopWeights))
    {
        // Storage for cuts across face
        Map<edge> faceSplitCuts(loop.size());
        // Storage for points on one side of cell
        labelList anchorPoints;

        okLoop =
            validLoop(cellI, loop, loopWeights, faceSplitCuts, anchorPoints);

        if (okLoop)
        {
            // Valid loop. Copy cellLoops and anchorPoints
            cellLoops_[cellI] = loop;
            cellAnchorPoints_[cellI].transfer(anchorPoints);

            // Copy split cuts
            forAllConstIter(Map<edge>, faceSplitCuts, iter)
            {
                faceSplitCut_.insert(iter.key(), iter());
            }


            // Update edgeIsCut, pointIsCut information
            forAll(loop, cutI)
            {
                label cut = loop[cutI];

                if (isEdge(cut))
                {
                    label edgeI = getEdge(cut);

                    edgeIsCut_[edgeI] = true;

                    edgeWeight_[edgeI] = loopWeights[cutI];
                }
                else
                {
                    label vertI = getVertex(cut);

                    pointIsCut_[vertI] = true;
                }
            }
        }
    }

    return okLoop;
}


void Foam::cellCuts::setFromCellLoops
(
    const labelList& cellLabels,
    const labelListList& cellLoops,
    const List<scalarField>& cellLoopWeights
)
{
    // 'Uncut' edges/vertices that are not used in loops.
    pointIsCut_ = false;

    edgeIsCut_ = false;

    forAll(cellLabels, cellLabelI)
    {
        label cellI = cellLabels[cellLabelI];

        const labelList& loop = cellLoops[cellLabelI];

        if (loop.size())
        {
            const scalarField& loopWeights = cellLoopWeights[cellLabelI];

            if (setFromCellLoop(cellI, loop, loopWeights))
            {
                // Valid loop. Call above will have upated all already.
            }
            else
            {
                // Clear cellLoops
                cellLoops_[cellI].setSize(0);
            }
        }
    }
}


void Foam::cellCuts::setFromCellCutter
(
    const cellLooper& cellCutter,
    const List<refineCell>& refCells
)
{
    // 'Uncut' edges/vertices that are not used in loops.
    pointIsCut_ = false;

    edgeIsCut_ = false;

    // storage for loop of cuts (cut vertices and/or cut edges)
    labelList cellLoop;
    scalarField cellLoopWeights;

    // For debugging purposes
    DynamicList<label> invalidCutCells(2);
    DynamicList<labelList> invalidCutLoops(2);
    DynamicList<scalarField> invalidCutLoopWeights(2);


    forAll(refCells, refCellI)
    {
        const refineCell& refCell = refCells[refCellI];

        label cellI = refCell.cellNo();

        const vector& refDir = refCell.direction();


        // Cut cell. Determines cellLoop and cellLoopWeights
        bool goodCut =
            cellCutter.cut
            (
                refDir,
                cellI,

                pointIsCut_,
                edgeIsCut_,
                edgeWeight_,

                cellLoop,
                cellLoopWeights
            );

        // Check whether edge refinement is on a per face basis compatible with
        // current pattern.
        if (goodCut)
        {
            if (setFromCellLoop(cellI, cellLoop, cellLoopWeights))
            {
                // Valid loop. Will have updated all info already.
            }
            else
            {
                cellLoops_[cellI].setSize(0);

                WarningIn("Foam::cellCuts::setFromCellCutter")
                    << "Found loop on cell " << cellI
                    << " that resulted in an unexpected bad cut."
                    << "    Suggestions:" << nl
                    << "      - Turn on the debug switch for 'cellCuts' to get"
                    << " geometry files that identify this cell." << nl
                    << "      - Also keep in mind to check the defined"
                    << " reference directions, as these are most likely the"
                    << " origin of the problem."
                    << nl << endl;

                // Discarded by validLoop
                if (debug)
                {
                    invalidCutCells.append(cellI);
                    invalidCutLoops.append(cellLoop);
                    invalidCutLoopWeights.append(cellLoopWeights);
                }
            }
        }
        else
        {
            // Clear cellLoops
            cellLoops_[cellI].setSize(0);
        }
    }

    if (debug && invalidCutCells.size())
    {
        invalidCutCells.shrink();
        invalidCutLoops.shrink();
        invalidCutLoopWeights.shrink();

        fileName cutsFile("invalidLoopCells.obj");

        Pout<< "cellCuts : writing inValidLoops cells to " << cutsFile << endl;

        OFstream cutsStream(cutsFile);

        meshTools::writeOBJ
        (
            cutsStream,
            mesh().cells(),
            mesh().faces(),
            mesh().points(),
            invalidCutCells
        );

        fileName loopsFile("invalidLoops.obj");

        Pout<< "cellCuts : writing inValidLoops loops to " << loopsFile << endl;

        OFstream loopsStream(loopsFile);

        label vertI = 0;

        forAll(invalidCutLoops, i)
        {
            writeOBJ
            (
                loopsStream,
                loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
                vertI
            );
        }
    }
}


void Foam::cellCuts::setFromCellCutter
(
    const cellLooper& cellCutter,
    const labelList& cellLabels,
    const List<plane>& cellCutPlanes
)
{
    // 'Uncut' edges/vertices that are not used in loops.
    pointIsCut_ = false;

    edgeIsCut_ = false;

    // storage for loop of cuts (cut vertices and/or cut edges)
    labelList cellLoop;
    scalarField cellLoopWeights;

    // For debugging purposes
    DynamicList<label> invalidCutCells(2);
    DynamicList<labelList> invalidCutLoops(2);
    DynamicList<scalarField> invalidCutLoopWeights(2);


    forAll(cellLabels, i)
    {
        label cellI = cellLabels[i];

        // Cut cell. Determines cellLoop and cellLoopWeights
        bool goodCut =
            cellCutter.cut
            (
                cellCutPlanes[i],
                cellI,

                pointIsCut_,
                edgeIsCut_,
                edgeWeight_,

                cellLoop,
                cellLoopWeights
            );

        // Check whether edge refinement is on a per face basis compatible with
        // current pattern.
        if (goodCut)
        {
            if (setFromCellLoop(cellI, cellLoop, cellLoopWeights))
            {
                // Valid loop. Will have updated all info already.
            }
            else
            {
                cellLoops_[cellI].setSize(0);

                // Discarded by validLoop
                if (debug)
                {
                    invalidCutCells.append(cellI);
                    invalidCutLoops.append(cellLoop);
                    invalidCutLoopWeights.append(cellLoopWeights);
                }
            }
        }
        else
        {
            // Clear cellLoops
            cellLoops_[cellI].setSize(0);
        }
    }

    if (debug && invalidCutCells.size())
    {
        invalidCutCells.shrink();
        invalidCutLoops.shrink();
        invalidCutLoopWeights.shrink();

        fileName cutsFile("invalidLoopCells.obj");

        Pout<< "cellCuts : writing inValidLoops cells to " << cutsFile << endl;

        OFstream cutsStream(cutsFile);

        meshTools::writeOBJ
        (
            cutsStream,
            mesh().cells(),
            mesh().faces(),
            mesh().points(),
            invalidCutCells
        );

        fileName loopsFile("invalidLoops.obj");

        Pout<< "cellCuts : writing inValidLoops loops to " << loopsFile << endl;

        OFstream loopsStream(loopsFile);

        label vertI = 0;

        forAll(invalidCutLoops, i)
        {
            writeOBJ
            (
                loopsStream,
                loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
                vertI
            );
        }
    }
}


void Foam::cellCuts::orientPlanesAndLoops()
{
    // Determine anchorPoints if not yet done by validLoop.
    forAll(cellLoops_, cellI)
    {
        const labelList& loop = cellLoops_[cellI];

        if (loop.size() && cellAnchorPoints_[cellI].empty())
        {
            // Leave anchor points empty if illegal loop.
            calcAnchors
            (
                cellI,
                loop,
                loopPoints(cellI),
                cellAnchorPoints_[cellI]
            );
        }
    }

    if (debug & 2)
    {
        Pout<< "cellAnchorPoints:" << endl;
    }
    forAll(cellAnchorPoints_, cellI)
    {
        if (cellLoops_[cellI].size())
        {
            if (cellAnchorPoints_[cellI].empty())
            {
                FatalErrorIn("orientPlanesAndLoops()")
                    << "No anchor points for cut cell " << cellI << endl
                    << "cellLoop:" << cellLoops_[cellI] << abort(FatalError);
            }

            if (debug & 2)
            {
                Pout<< "    cell:" << cellI << " anchored at "
                    << cellAnchorPoints_[cellI] << endl;
            }
        }
    }

    // Calculate number of valid cellLoops
    nLoops_ = 0;

    forAll(cellLoops_, cellI)
    {
        if (cellLoops_[cellI].size())
        {
            nLoops_++;
        }
    }
}


void Foam::cellCuts::calcLoopsAndAddressing(const labelList& cutCells)
{
    // Sanity check on weights
    forAll(edgeIsCut_, edgeI)
    {
        if (edgeIsCut_[edgeI])
        {
            scalar weight = edgeWeight_[edgeI];

            if (weight < 0 || weight > 1)
            {
                FatalErrorIn
                (
                    "cellCuts::calcLoopsAndAddressing(const labelList&)"
                )   << "Weight out of range [0,1]. Edge " << edgeI
                    << " verts:" << mesh().edges()[edgeI]
                    << " weight:" << weight << abort(FatalError);
            }
        }
        else
        {
            // Weight not used. Set to illegal value to make any use fall over.
            edgeWeight_[edgeI] = -GREAT;
        }
    }


    // Calculate faces that split cells in two
    calcCellLoops(cutCells);

    if (debug & 2)
    {
        Pout<< "-- cellLoops --" << endl;
        forAll(cellLoops_, cellI)
        {
            const labelList& loop = cellLoops_[cellI];

            if (loop.size())
            {
                Pout<< "cell:" << cellI << "  ";
                writeCuts(Pout, loop, loopWeights(loop));
                Pout<< endl;
            }
        }
    }

    // Redo basic cut information (pointIsCut, edgeIsCut, faceSplitCut)
    // using cellLoop only.
    setFromCellLoops();
}


void Foam::cellCuts::check() const
{
    // Check weights for unsnapped values
    forAll(edgeIsCut_, edgeI)
    {
        if (edgeIsCut_[edgeI])
        {
            if
            (
                edgeWeight_[edgeI] < geomCellLooper::snapTol()
             || edgeWeight_[edgeI] > (1 - geomCellLooper::snapTol())
            )
            {
                // Should have been snapped.
                //FatalErrorIn("cellCuts::check()")
                WarningIn("cellCuts::check()")
                    << "edge:" << edgeI << " vertices:"
                    << mesh().edges()[edgeI]
                    << " weight:" << edgeWeight_[edgeI] << " should have been"
                    << " snapped to one of its endpoints"
                    //<< abort(FatalError);
                    << endl;
            }
        }
        else
        {
            if (edgeWeight_[edgeI] > - 1)
            {
                FatalErrorIn("cellCuts::check()")
                    << "edge:" << edgeI << " vertices:"
                    << mesh().edges()[edgeI]
                    << " weight:" << edgeWeight_[edgeI] << " is not cut but"
                    << " its weight is not marked invalid"
                    << abort(FatalError);
            }
        }
    }

    // Check that all elements of cellloop are registered
    forAll(cellLoops_, cellI)
    {
        const labelList& loop = cellLoops_[cellI];

        forAll(loop, i)
        {
            label cut = loop[i];

            if
            (
                (isEdge(cut) && !edgeIsCut_[getEdge(cut)])
             || (!isEdge(cut) && !pointIsCut_[getVertex(cut)])
            )
            {
                labelList cuts(1, cut);
                writeCuts(Pout, cuts, loopWeights(cuts));

                FatalErrorIn("cellCuts::check()")
                    << "cell:" << cellI << " loop:"
                    << loop
                    << " cut:" << cut << " is not marked as cut"
                    << abort(FatalError);
            }
        }
    }

    // Check that no elements of cell loop are anchor point.
    forAll(cellLoops_, cellI)
    {
        const labelList& anchors = cellAnchorPoints_[cellI];

        const labelList& loop = cellLoops_[cellI];

        if (loop.size() && anchors.empty())
        {
            FatalErrorIn("cellCuts::check()")
                << "cell:" << cellI << " loop:" << loop
                << " has no anchor points"
                << abort(FatalError);
        }


        forAll(loop, i)
        {
            label cut = loop[i];

            if
            (
                !isEdge(cut)
             && findIndex(anchors, getVertex(cut)) != -1
            )
            {
                FatalErrorIn("cellCuts::check()")
                    << "cell:" << cellI << " loop:" << loop
                    << " anchor points:" << anchors
                    << " anchor:" << getVertex(cut) << " is part of loop"
                    << abort(FatalError);
            }
        }
    }


    // Check that cut faces have a neighbour that is cut.
    forAllConstIter(Map<edge>, faceSplitCut_, iter)
    {
        label faceI = iter.key();

        if (mesh().isInternalFace(faceI))
        {
            label own = mesh().faceOwner()[faceI];
            label nei = mesh().faceNeighbour()[faceI];

            if (cellLoops_[own].empty() && cellLoops_[nei].empty())
            {
                FatalErrorIn("cellCuts::check()")
                    << "Internal face:" << faceI << " cut by " << iter()
                    << " has owner:" << own
                    << " and neighbour:" << nei
                    << " that are both uncut"
                    << abort(FatalError);
            }
        }
        else
        {
            label own = mesh().faceOwner()[faceI];

            if (cellLoops_[own].empty())
            {
                FatalErrorIn("cellCuts::check()")
                    << "Boundary face:" << faceI << " cut by " << iter()
                    << " has owner:" << own
                    << " that is uncut"
                    << abort(FatalError);
            }
        }
    }
}


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

Foam::cellCuts::cellCuts
(
    const polyMesh& mesh,
    const labelList& cutCells,
    const labelList& meshVerts,
    const labelList& meshEdges,
    const scalarField& meshEdgeWeights
)
:
    edgeVertex(mesh),
    pointIsCut_(expand(mesh.nPoints(), meshVerts)),
    edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
    edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
    faceCutsPtr_(NULL),
    faceSplitCut_(cutCells.size()),
    cellLoops_(mesh.nCells()),
    nLoops_(-1),
    cellAnchorPoints_(mesh.nCells())
{
    if (debug)
    {
        Pout<< "cellCuts : constructor from cut verts and edges" << endl;
    }

    calcLoopsAndAddressing(cutCells);

    // Calculate planes and flip cellLoops if necessary
    orientPlanesAndLoops();

    if (debug)
    {
        check();
    }

    clearOut();

    if (debug)
    {
        Pout<< "cellCuts : leaving constructor from cut verts and edges"
            << endl;
    }
}


Foam::cellCuts::cellCuts
(
    const polyMesh& mesh,
    const labelList& meshVerts,
    const labelList& meshEdges,
    const scalarField& meshEdgeWeights
)
:
    edgeVertex(mesh),
    pointIsCut_(expand(mesh.nPoints(), meshVerts)),
    edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
    edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
    faceCutsPtr_(NULL),
    faceSplitCut_(mesh.nFaces()/10 + 1),
    cellLoops_(mesh.nCells()),
    nLoops_(-1),
    cellAnchorPoints_(mesh.nCells())
{
    // Construct from pattern of cuts. Finds out itself which cells are cut.
    // (can go wrong if e.g. all neighbours of cell are refined)

    if (debug)
    {
        Pout<< "cellCuts : constructor from cellLoops" << endl;
    }

    calcLoopsAndAddressing(identity(mesh.nCells()));

    // Calculate planes and flip cellLoops if necessary
    orientPlanesAndLoops();

    if (debug)
    {
        check();
    }

    clearOut();

    if (debug)
    {
        Pout<< "cellCuts : leaving constructor from cellLoops" << endl;
    }
}


Foam::cellCuts::cellCuts
(
    const polyMesh& mesh,
    const labelList& cellLabels,
    const labelListList& cellLoops,
    const List<scalarField>& cellEdgeWeights
)
:
    edgeVertex(mesh),
    pointIsCut_(mesh.nPoints(), false),
    edgeIsCut_(mesh.nEdges(), false),
    edgeWeight_(mesh.nEdges(), -GREAT),
    faceCutsPtr_(NULL),
    faceSplitCut_(cellLabels.size()),
    cellLoops_(mesh.nCells()),
    nLoops_(-1),
    cellAnchorPoints_(mesh.nCells())
{
    if (debug)
    {
        Pout<< "cellCuts : constructor from cellLoops" << endl;
    }

    // Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
    // Makes sure cuts are consistent
    setFromCellLoops(cellLabels, cellLoops, cellEdgeWeights);

    // Calculate planes and flip cellLoops if necessary
    orientPlanesAndLoops();

    if (debug)
    {
        check();
    }

    clearOut();

    if (debug)
    {
        Pout<< "cellCuts : leaving constructor from cellLoops" << endl;
    }
}


Foam::cellCuts::cellCuts
(
    const polyMesh& mesh,
    const cellLooper& cellCutter,
    const List<refineCell>& refCells
)
:
    edgeVertex(mesh),
    pointIsCut_(mesh.nPoints(), false),
    edgeIsCut_(mesh.nEdges(), false),
    edgeWeight_(mesh.nEdges(), -GREAT),
    faceCutsPtr_(NULL),
    faceSplitCut_(refCells.size()),
    cellLoops_(mesh.nCells()),
    nLoops_(-1),
    cellAnchorPoints_(mesh.nCells())
{
    if (debug)
    {
        Pout<< "cellCuts : constructor from cellCutter" << endl;
    }

    // Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
    // Makes sure cuts are consistent
    setFromCellCutter(cellCutter, refCells);

    // Calculate planes and flip cellLoops if necessary
    orientPlanesAndLoops();

    if (debug)
    {
        check();
    }

    clearOut();

    if (debug)
    {
        Pout<< "cellCuts : leaving constructor from cellCutter" << endl;
    }
}


Foam::cellCuts::cellCuts
(
    const polyMesh& mesh,
    const cellLooper& cellCutter,
    const labelList& cellLabels,
    const List<plane>& cutPlanes
)
:
    edgeVertex(mesh),
    pointIsCut_(mesh.nPoints(), false),
    edgeIsCut_(mesh.nEdges(), false),
    edgeWeight_(mesh.nEdges(), -GREAT),
    faceCutsPtr_(NULL),
    faceSplitCut_(cellLabels.size()),
    cellLoops_(mesh.nCells()),
    nLoops_(-1),
    cellAnchorPoints_(mesh.nCells())
{
    if (debug)
    {
        Pout<< "cellCuts : constructor from cellCutter with prescribed plane"
            << endl;
    }

    // Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
    // Makes sure cuts are consistent
    setFromCellCutter(cellCutter, cellLabels, cutPlanes);

    // Calculate planes and flip cellLoops if necessary
    orientPlanesAndLoops();

    if (debug)
    {
        check();
    }

    clearOut();

    if (debug)
    {
        Pout<< "cellCuts : leaving constructor from cellCutter with prescribed"
            << " plane" << endl;
    }
}


Foam::cellCuts::cellCuts
(
    const polyMesh& mesh,
    const boolList& pointIsCut,
    const boolList& edgeIsCut,
    const scalarField& edgeWeight,
    const Map<edge>& faceSplitCut,
    const labelListList& cellLoops,
    const label nLoops,
    const labelListList& cellAnchorPoints
)
:
    edgeVertex(mesh),
    pointIsCut_(pointIsCut),
    edgeIsCut_(edgeIsCut),
    edgeWeight_(edgeWeight),
    faceCutsPtr_(NULL),
    faceSplitCut_(faceSplitCut),
    cellLoops_(cellLoops),
    nLoops_(nLoops),
    cellAnchorPoints_(cellAnchorPoints)
{
    if (debug)
    {
        Pout<< "cellCuts : constructor from components" << endl;
        Pout<< "cellCuts : leaving constructor from components" << endl;
    }
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::cellCuts::~cellCuts()
{
    clearOut();
}


void Foam::cellCuts::clearOut()
{
    deleteDemandDrivenData(faceCutsPtr_);
}


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

Foam::pointField Foam::cellCuts::loopPoints(const label cellI) const
{
    const labelList& loop = cellLoops_[cellI];

    pointField loopPts(loop.size());

    forAll(loop, fp)
    {
        label cut = loop[fp];

        if (isEdge(cut))
        {
            label edgeI = getEdge(cut);

            loopPts[fp] = coord(cut, edgeWeight_[edgeI]);
        }
        else
        {
            loopPts[fp] = coord(cut, -GREAT);
        }
    }
    return loopPts;
}


void Foam::cellCuts::flip(const label cellI)
{
    labelList& loop = cellLoops_[cellI];

    reverse(loop);

    // Reverse anchor point set.
    cellAnchorPoints_[cellI] =
        nonAnchorPoints
        (
            mesh().cellPoints()[cellI],
            cellAnchorPoints_[cellI],
            loop
        );
}


void Foam::cellCuts::flipLoopOnly(const label cellI)
{
    labelList& loop = cellLoops_[cellI];

    reverse(loop);
}


void Foam::cellCuts::writeOBJ
(
    Ostream& os,
    const pointField& loopPts,
    label& vertI
) const
{
    label startVertI = vertI;

    forAll(loopPts, fp)
    {
        const point& pt = loopPts[fp];

        os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;

        vertI++;
    }

    os  << 'f';
    forAll(loopPts, fp)
    {
        os  << ' ' << startVertI + fp + 1;
    }
    os  << endl;
}


void Foam::cellCuts::writeOBJ(Ostream& os) const
{
    label vertI = 0;

    forAll(cellLoops_, cellI)
    {
        writeOBJ(os, loopPoints(cellI), vertI);
    }
}


void Foam::cellCuts::writeCellOBJ(const fileName& dir, const label cellI) const
{
    const labelList& anchors = cellAnchorPoints_[cellI];

    writeOBJ(dir, cellI, loopPoints(cellI), anchors);
}


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

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 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::cellCuts

Description
    Description of cuts across cells.

    Description of cut is given as list of vertices and
    list of edges to be cut (and position on edge).

    Does some checking of correctness/non-overlapping of cuts. 2x2x2
    refinement has to be done in three passes since cuts can not overlap
    (would make addressing too complicated)

    Introduces concept of 'cut' which is either an existing vertex
    or a edge.

    Input can either be
    -# list of cut vertices and list of cut edges. Constructs cell
       circumference walks ('cellLoops').

    -# list of cell circumference walks. Will filter them so they
       don't overlap.

    -# cellWalker and list of cells to refine (refineCell). Constructs
       cellLoops and does B. cellWalker is class which can cut a single
       cell using a plane through the cell centre and in a certain normal
       direction

    CellCuts constructed from cellLoops (B, C) can have multiple cut-edges
    and/or cut-point as long as there is per face only one (or none) cut across
    a face, i.e. a face can only be split into two faces.

    The information available after construction:
      - pointIsCut, edgeIsCut.
      - faceSplitCut : the cross-cut of a face.
      - cellLoops : the loop which cuts across a cell
      - cellAnchorPoints : per cell the vertices on one side which are
        considered the anchor points.

    AnchorPoints: connected loops have to be oriented in the same way to
    be able to grow new internal faces out of the same bottom faces.
    (limitation of the mapping procedure). The loop is cellLoop is oriented
    such that the normal of it points towards the anchorPoints.
    This is the only routine which uses geometric information.


    Limitations:
    - cut description is very precise. Hard to get right.
    - do anchor points need to be unique per cell? Very inefficient.
    - is orientation guaranteed to be correct? Uses geometric info so can go
      wrong on highly distorted meshes.
    - is memory inefficient. Probably need to use Maps instead of
      labelLists.

SourceFiles
    cellCuts.C

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

#ifndef cellCuts_H
#define cellCuts_H

#include "edgeVertex.H"
#include "labelList.H"
#include "boolList.H"
#include "scalarField.H"
#include "pointField.H"
#include "DynamicList.H"
#include "typeInfo.H"

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

namespace Foam
{

// Forward declaration of classes
class polyMesh;
class cellLooper;
class refineCell;
class plane;

/*---------------------------------------------------------------------------*\
                           Class cellCuts Declaration
\*---------------------------------------------------------------------------*/

class cellCuts
:
    public edgeVertex
{
    // Private data

        // Per point/edge status

            //- Is mesh point cut
            boolList pointIsCut_;

            //- Is edge cut
            boolList edgeIsCut_;

            //- If edge is cut gives weight (0->start() to 1->end())
            scalarField edgeWeight_;


        // Cut addressing

            //- Cuts per existing face (includes those along edge of face)
            //  Cuts in no particular order.
            mutable labelListList* faceCutsPtr_;

            //- Per face : cut across edge (so not along existing edge)
            //  (can only be one per face)
            Map<edge> faceSplitCut_;


        // Cell-cut addressing

            //- Loop across cell circumference
            labelListList cellLoops_;

            //- Number of valid loops in cellLoops_
            label nLoops_;

            //- For each cut cell the points on the 'anchor' side of the cell.
            labelListList cellAnchorPoints_;


    // Private Static Functions

        //- Find value in first n elements of list.
        static label findPartIndex
        (
            const labelList&,
            const label n,
            const label val
        );

        //- Create boolList with all labels specified set to true
        //  (and rest to false)
        static boolList expand(const label size, const labelList& labels);

        //- Create scalarField with all specified labels set to corresponding
        //  value in scalarField.
        static scalarField expand
        (
            const label,
            const labelList&,
            const scalarField&
        );

        //- Returns -1 or index of first element of lst that cannot be found
        //  in map.
        static label firstUnique
        (
            const labelList& lst,
            const Map<label>&
        );

    // Private Member Functions

        //- Debugging: write cell's edges and any cut vertices and edges
        //  (so no cell loop determined yet)
        void writeUncutOBJ(const fileName&, const label cellI) const;

        //- Debugging: write cell's edges, loop and anchors to directory.
        void writeOBJ
        (
            const fileName& dir,
            const label cellI,
            const pointField& loopPoints,
            const labelList& anchors
        ) const;

        //- Find face on cell using the two edges.
        label edgeEdgeToFace
        (
            const label cellI,
            const label edgeA,
            const label edgeB
        ) const;


        //- Find face on cell using an edge and a vertex.
        label edgeVertexToFace
        (
            const label cellI,
            const label edgeI,
            const label vertI
        ) const;

        //- Find face using two vertices (guaranteed not to be along edge)
        label vertexVertexToFace
        (
            const label cellI,
            const label vertA,
            const label vertB
        ) const;


        // Cut addressing

            //- Calculate faceCuts in face vertex order.
            void calcFaceCuts() const;


        // Loop (cuts on cell circumference) calculation

            //- Find edge (or -1) on faceI using vertices v0,v1
            label findEdge
            (
                const label faceI,
                const label v0,
                const label v1
            ) const;

            //- Find face on which all cuts are (very rare) or -1.
            label loopFace(const label cellI, const labelList& loop) const;

            //- Cross otherCut into next faces (not exclude0, exclude1)
            bool walkPoint
            (
                const label cellI,
                const label startCut,

                const label exclude0,
                const label exclude1,

                const label otherCut,

                label& nVisited,
                labelList& visited
            ) const;

            //- Cross cut (which is edge on faceI) onto next face
            bool crossEdge
            (
                const label cellI,
                const label startCut,
                const label faceI,
                const label otherCut,

                label& nVisited,
                labelList& visited
            ) const;

            // wrapper around visited[nVisited++] = cut. Checks for duplicate
            // cuts.
            bool addCut
            (
                const label cellI,
                const label cut,
                label& nVisited,
                labelList& visited
            ) const;

            //- Walk across faceI following cuts, starting at cut. Stores cuts
            //  visited
            // Returns true if valid walk.
            bool walkFace
            (
                const label cellI,
                const label startCut,
                const label faceI,
                const label cut,

                label& lastCut,
                label& beforeLastCut,
                label& nVisited,
                labelList& visited
            ) const;

            //- Walk across cuts (cut edges or cut vertices) of cell. Stops when
            //  hit cut  already visited. Returns true when loop of 3 or more
            //  vertices found.
            bool walkCell
            (
                const label cellI,
                const label startCut,   // overall starting cut
                const label faceI,
                const label prevCut,    // current cut
                label& nVisited,
                labelList& visited
            ) const;

            //- Determine for every cut cell the face it is cut by.
            void calcCellLoops(const labelList& cutCells);


        // Cell anchoring

            //- Are there enough faces on anchor side of cellI?
            bool checkFaces
            (
                const label cellI,
                const labelList& anchorPoints
            ) const;

            //- Walk unset edges of single cell from starting point and
            //  marks visited edges and vertices with status.
            void walkEdges
            (
                const label cellI,
                const label pointI,
                const label status,

                Map<label>& edgeStatus,
                Map<label>& pointStatus
            ) const;

            //- Check anchor points on 'outside' of loop
            bool loopAnchorConsistent
            (
                const label cellI,
                const pointField& loopPts,
                const labelList& anchorPoints
            ) const;

            //- Determines set of anchor points given a loop. The loop should
            //  split the cell into two. Returns true if a valid set of anchor
            //  points determined, false otherwise.
            bool calcAnchors
            (
                const label cellI,
                const labelList& loop,
                const pointField& loopPts,

                labelList& anchorPoints
            ) const;

            //- Returns coordinates of points on loop with explicitly provided
            //  weights.
            pointField loopPoints
            (
                const labelList& loop,
                const scalarField& loopWeights
            ) const;

            //- Returns weights of loop. Inverse of loopPoints.
            scalarField loopWeights(const labelList& loop) const;

            //- Check if cut edges in loop are compatible with ones in
            //  edgeIsCut_
            bool validEdgeLoop
            (
                const labelList& loop,
                const scalarField& loopWeights
            ) const;

            //- Counts number of cuts on face.
            label countFaceCuts
            (
                const label faceI,
                const labelList& loop
            ) const;

            //- Determine compatibility of loop with existing cut pattern.
            //  Does not use cut-addressing (faceCuts_, cutCuts_)
            bool conservativeValidLoop
            (
                const label cellI,
                const labelList& loop
            ) const;

            //- Check if loop is compatible with existing cut pattern in
            //  pointIsCut, edgeIsCut, faceSplitCut.
            //  Calculates and returns for current cell the cut faces and the
            //  points on one side of the loop.
            bool validLoop
            (
                const label cellI,
                const labelList& loop,
                const scalarField& loopWeights,
                Map<edge>& newFaceSplitCut,
                labelList& anchorPoints
            ) const;

            //- Update basic cut information from cellLoops.
            //  Assumes cellLoops_ and edgeWeight_ already set and consistent.
            //  Does not use any other information.
            void setFromCellLoops();

            //- Update basic cut information for single cell from cellLoop.
            bool setFromCellLoop
            (
                const label cellI,
                const labelList& loop,
                const scalarField& loopWeights
            );

            //- Update basic cut information from cellLoops. Checks for
            //  consistency with existing cut pattern.
            void setFromCellLoops
            (
                const labelList& cellLabels,
                const labelListList& cellLoops,
                const List<scalarField>& cellLoopWeights
            );

            //- Cut cells and update basic cut information from cellLoops.
            //  Checks each loop for consistency with existing cut pattern.
            void setFromCellCutter
            (
                const cellLooper&,
                const List<refineCell>& refCells
            );

            //- Same as above but now cut with prescribed plane (instead of
            //  just normal).
            void setFromCellCutter
            (
                const cellLooper&,
                const labelList& cellLabels,
                const List<plane>&
            );

            //- Set orientation of loops
            void orientPlanesAndLoops();

            //- Top level driver: adressing calculation and loop detection
            //  (loops splitting cells).
            void calcLoopsAndAddressing(const labelList& cutCells);

            //- Check various consistencies.
            void check() const;


        //- Disallow default bitwise copy construct
        cellCuts(const cellCuts&);

        //- Disallow default bitwise assignment
        void operator=(const cellCuts&);


public:

    //- Runtime type information
    ClassName("cellCuts");


    // Constructors

        //- Construct from cells to cut and pattern of cuts
        cellCuts
        (
            const polyMesh& mesh,
            const labelList& cutCells,
            const labelList& meshVerts,
            const labelList& meshEdges,
            const scalarField& meshEdgeWeights
        );

        //- Construct from pattern of cuts. Detect cells to cut.
        cellCuts
        (
            const polyMesh& mesh,
            const labelList& meshVerts,
            const labelList& meshEdges,
            const scalarField& meshEdgeWeights
        );

        //- Construct from complete cellLoops through specified cells.
        //  Checks for consistency.
        //  Constructs cut-cut addressing and cellAnchorPoints.
        cellCuts
        (
            const polyMesh& mesh,
            const labelList& cellLabels,
            const labelListList& cellLoops,
            const List<scalarField>& cellEdgeWeights
        );

        //- Construct from list of cells to cut and direction to cut in
        //  (always through cell centre) and celllooper.
        cellCuts
        (
            const polyMesh& mesh,
            const cellLooper& cellCutter,
            const List<refineCell>& refCells
        );

        //- Construct from list of cells to cut and plane to cut with and
        //  celllooper. (constructor above always cuts through cell centre)
        cellCuts
        (
            const polyMesh& mesh,
            const cellLooper& cellCutter,
            const labelList& cellLabels,
            const List<plane>& cutPlanes
        );

        //- Construct from components
        cellCuts
        (
            const polyMesh& mesh,
            const boolList& pointIsCut,
            const boolList& edgeIsCut,
            const scalarField& edgeWeight,
            const Map<edge>& faceSplitCut,
            const labelListList& cellLoops,
            const label nLoops,
            const labelListList& cellAnchorPoints
        );


    //- Destructor
    ~cellCuts();

        //- Clear out demand driven storage
        void clearOut();


    // Member Functions

        // Access

            //- Is mesh point cut
            const boolList& pointIsCut() const
            {
                return pointIsCut_;
            }

            //- Is edge cut
            const boolList& edgeIsCut() const
            {
                return edgeIsCut_;
            }

            //- If edge is cut gives weight (ratio between start() and end())
            const scalarField& edgeWeight() const
            {
                return edgeWeight_;
            }

            //- Cuts per existing face (includes those along edge of face)
            //  Cuts in no particular order
            const labelListList& faceCuts() const
            {
                if (!faceCutsPtr_)
                {
                    calcFaceCuts();
                }
                return *faceCutsPtr_;
            }

            //- Gives for split face the two cuts that split the face into two.
            const Map<edge>& faceSplitCut() const
            {
                return faceSplitCut_;
            }

            //- For each cut cell the cut along the circumference.
            const labelListList& cellLoops() const
            {
                return cellLoops_;
            }

            //- Number of valid cell loops
            label nLoops() const
            {
                return nLoops_;
            }

            //- For each cut cell the points on the 'anchor' side of the cell.
            const labelListList& cellAnchorPoints() const
            {
                return cellAnchorPoints_;
            }


        // Other

            //- Returns coordinates of points on loop for given cell.
            //  Uses cellLoops_ and edgeWeight_
            pointField loopPoints(const label cellI) const;

            //- Invert anchor point selection.
            labelList nonAnchorPoints
            (
                const labelList& cellPoints,
                const labelList& anchorPoints,
                const labelList& loop
            ) const;

            //- Flip loop for cellI. Updates anchor points as well.
            void flip(const label cellI);

            //- Flip loop for cellI. Does not update anchors. Use with care
            //  (only if you're sure loop orientation is wrong)
            void flipLoopOnly(const label cellI);


        // Write

            //- debugging:Write list of cuts to stream in OBJ format
            void writeOBJ
            (
                Ostream& os,
                const pointField& loopPoints,
                label& vertI
            ) const;

            //- debugging:Write all of cuts to stream in OBJ format
            void writeOBJ(Ostream& os) const;

            //- debugging:Write edges of cell and loop
            void writeCellOBJ(const fileName& dir, const label cellI) const;

};


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

} // End namespace Foam

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

#endif

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

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 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 "cellCuts.H"
#include "polyMesh.H"
#include "Time.H"
#include "ListOps.H"
#include "cellLooper.H"
#include "refineCell.H"
#include "meshTools.H"
#include "geomCellLooper.H"
#include "OFstream.H"
#include "plane.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
defineTypeNameAndDebug(cellCuts, 0);
}


// * * * * * * * * * * * * * Private Static Functions  * * * * * * * * * * * //

Foam::label Foam::cellCuts::findPartIndex
(
    const labelList& elems,
    const label nElems,
    const label val
)
{
    for (label i = 0; i < nElems; i++)
    {
        if (elems[i] == val)
        {
            return i;
        }
    }
    return -1;
}


Foam::boolList Foam::cellCuts::expand
(
    const label size,
    const labelList& labels
)
{
    boolList result(size, false);

    forAll(labels, labelI)
    {
        result[labels[labelI]] = true;
    }
    return result;
}


Foam::scalarField Foam::cellCuts::expand
(
    const label size,
    const labelList& labels,
    const scalarField& weights
)
{
    scalarField result(size, -GREAT);

    forAll(labels, labelI)
    {
        result[labels[labelI]] = weights[labelI];
    }
    return result;
}


Foam::label Foam::cellCuts::firstUnique
(
    const labelList& lst,
    const Map<label>& map
)
{
    forAll(lst, i)
    {
        if (!map.found(lst[i]))
        {
            return i;
        }
    }
    return -1;
}


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

void Foam::cellCuts::writeUncutOBJ
(
    const fileName& dir,
    const label cellI
) const
{
    //- Cell edges
    OFstream cutsStream(dir / "cell_" + name(cellI) + ".obj");

    Pout<< "Writing cell for time " <<  mesh().time().timeName()
        << " to " << cutsStream.name() << nl;

    meshTools::writeOBJ
    (
        cutsStream,
        mesh().cells(),
        mesh().faces(),
        mesh().points(),
        labelList(1, cellI)
    );

    //- Loop cutting cell in two
    OFstream cutStream(dir / "cellCuts_" + name(cellI) + ".obj");

    Pout<< "Writing raw cuts on cell for time " <<  mesh().time().timeName()
        << " to " << cutStream.name() << nl;

    const labelList& cPoints = mesh().cellPoints()[cellI];

    forAll(cPoints, i)
    {
        label pointI = cPoints[i];
        if (pointIsCut_[pointI])
        {
            meshTools::writeOBJ(cutStream, mesh().points()[pointI]);
        }
    }

    const pointField& pts = mesh().points();

    const labelList& cEdges = mesh().cellEdges()[cellI];

    forAll(cEdges, i)
    {
        label edgeI = cEdges[i];

        if (edgeIsCut_[edgeI])
        {
            const edge& e = mesh().edges()[edgeI];

            const scalar w = edgeWeight_[edgeI];

            meshTools::writeOBJ(cutStream, w*pts[e[1]] + (1-w)*pts[e[0]]);
        }
    }
}


void Foam::cellCuts::writeOBJ
(
    const fileName& dir,
    const label cellI,
    const pointField& loopPoints,
    const labelList& anchors
) const
{
    //- Cell edges
    OFstream cutsStream(dir / "cell_" + name(cellI) + ".obj");

    Pout<< "Writing cell for time " <<  mesh().time().timeName()
        << " to " << cutsStream.name() << nl;

    meshTools::writeOBJ
    (
        cutsStream,
        mesh().cells(),
        mesh().faces(),
        mesh().points(),
        labelList(1, cellI)
    );


    //- Loop cutting cell in two
    OFstream loopStream(dir / "cellLoop_" + name(cellI) + ".obj");

    Pout<< "Writing loop for time " <<  mesh().time().timeName()
        << " to " << loopStream.name() << nl;

    label vertI = 0;

    writeOBJ(loopStream, loopPoints, vertI);


    //- Anchors for cell
    OFstream anchorStream(dir / "anchors_" + name(cellI) + ".obj");

    Pout<< "Writing anchors for time " <<  mesh().time().timeName()
        << " to " << anchorStream.name() << endl;

    forAll(anchors, i)
    {
        meshTools::writeOBJ(anchorStream, mesh().points()[anchors[i]]);
    }
}


Foam::label Foam::cellCuts::edgeEdgeToFace
(
    const label cellI,
    const label edgeA,
    const label edgeB
) const
{
    const labelList& cFaces = mesh().cells()[cellI];

    forAll(cFaces, cFaceI)
    {
        label faceI = cFaces[cFaceI];

        const labelList& fEdges = mesh().faceEdges()[faceI];

        if
        (
            findIndex(fEdges, edgeA) != -1
         && findIndex(fEdges, edgeB) != -1
        )
        {
           return faceI;
        }
    }

    // Coming here means the loop is illegal since the two edges
    // are not shared by a face. We just mark loop as invalid and continue.

    WarningInFunction
        << "cellCuts : Cannot find face on cell "
        << cellI << " that has both edges " << edgeA << ' ' << edgeB << endl
        << "faces : " << cFaces << endl
        << "edgeA : " << mesh().edges()[edgeA] << endl
        << "edgeB : " << mesh().edges()[edgeB] << endl
        << "Marking the loop across this cell as invalid" << endl;

    return -1;
}


Foam::label Foam::cellCuts::edgeVertexToFace
(
    const label cellI,
    const label edgeI,
    const label vertI
) const
{
    const labelList& cFaces = mesh().cells()[cellI];

    forAll(cFaces, cFaceI)
    {
        label faceI = cFaces[cFaceI];

        const face& f = mesh().faces()[faceI];

        const labelList& fEdges = mesh().faceEdges()[faceI];

        if
        (
            findIndex(fEdges, edgeI) != -1
         && findIndex(f, vertI) != -1
        )
        {
           return faceI;
        }
    }

    WarningInFunction
        << "cellCuts : Cannot find face on cell "
        << cellI << " that has both edge " << edgeI << " and vertex "
        << vertI << endl
        << "faces : " << cFaces << endl
        << "edge : " << mesh().edges()[edgeI] << endl
        << "Marking the loop across this cell as invalid" << endl;

    return -1;
}


Foam::label Foam::cellCuts::vertexVertexToFace
(
    const label cellI,
    const label vertA,
    const label vertB
) const
{
    const labelList& cFaces = mesh().cells()[cellI];

    forAll(cFaces, cFaceI)
    {
        label faceI = cFaces[cFaceI];

        const face& f = mesh().faces()[faceI];

        if
        (
            findIndex(f, vertA) != -1
         && findIndex(f, vertB) != -1
        )
        {
           return faceI;
        }
    }

    WarningInFunction
        << "cellCuts : Cannot find face on cell "
        << cellI << " that has vertex " << vertA << " and vertex "
        << vertB << endl
        << "faces : " << cFaces << endl
        << "Marking the loop across this cell as invalid" << endl;

    return -1;
}


void Foam::cellCuts::calcFaceCuts() const
{
    if (faceCutsPtr_)
    {
        FatalErrorInFunction
            << "faceCuts already calculated" << abort(FatalError);
    }

    const faceList& faces = mesh().faces();


    faceCutsPtr_ = new labelListList(mesh().nFaces());
    labelListList& faceCuts = *faceCutsPtr_;

    for (label faceI = 0; faceI < mesh().nFaces(); faceI++)
    {
        const face& f = faces[faceI];

        // Big enough storage (possibly all points and all edges cut). Shrink
        // later on.
        labelList& cuts = faceCuts[faceI];

        cuts.setSize(2*f.size());

        label cutI = 0;

        // Do point-edge-point walk over face and collect all cuts.
        // Problem is that we want to start from one of the endpoints of a
        // string of connected cuts; we don't want to start somewhere in the
        // middle.

        // Pass1: find first point cut not preceeded by a cut.
        label startFp = -1;

        forAll(f, fp)
        {
            label v0 = f[fp];

            if (pointIsCut_[v0])
            {
                label vMin1 = f[f.rcIndex(fp)];

                if
                (
                    !pointIsCut_[vMin1]
                 && !edgeIsCut_[findEdge(faceI, v0, vMin1)]
                )
                {
                    cuts[cutI++] = vertToEVert(v0);
                    startFp = f.fcIndex(fp);
                    break;
                }
            }
        }

        // Pass2: first edge cut not preceeded by point cut
        if (startFp == -1)
        {
            forAll(f, fp)
            {
                label fp1 = f.fcIndex(fp);

                label v0 = f[fp];
                label v1 = f[fp1];

                label edgeI = findEdge(faceI, v0, v1);

                if (edgeIsCut_[edgeI] && !pointIsCut_[v0])
                {
                    cuts[cutI++] = edgeToEVert(edgeI);
                    startFp = fp1;
                    break;
                }
            }
        }

        // Pass3: nothing found so far. Either face is not cut at all or all
        // vertices are cut. Start from 0.
        if (startFp == -1)
        {
            startFp = 0;
        }

        // Store all cuts starting from startFp;
        label fp = startFp;

        bool allVerticesCut = true;

        forAll(f, i)
        {
            label fp1 = f.fcIndex(fp);

            // Get the three items: current vertex, next vertex and edge
            // inbetween
            label v0 = f[fp];
            label v1 = f[fp1];
            label edgeI = findEdge(faceI, v0, v1);

            if (pointIsCut_[v0])
            {
                cuts[cutI++] = vertToEVert(v0);
            }
            else
            {
                // For check if all vertices have been cut (= illegal)
                allVerticesCut = false;
            }

            if (edgeIsCut_[edgeI])
            {
                cuts[cutI++] = edgeToEVert(edgeI);
            }

            fp = fp1;
        }


        if (allVerticesCut)
        {
            WarningInFunction
                << "Face " << faceI << " vertices " << f
                << " has all its vertices cut. Not cutting face." << endl;

            cutI = 0;
        }

        // Remove duplicate starting point
        if (cutI > 1 && cuts[cutI-1] == cuts[0])
        {
            cutI--;
        }
        cuts.setSize(cutI);
    }
}


Foam::label Foam::cellCuts::findEdge
(
    const label faceI,
    const label v0,
    const label v1
) const
{
    const edgeList& edges = mesh().edges();

    const labelList& fEdges = mesh().faceEdges()[faceI];

    forAll(fEdges, i)
    {
        const edge& e = edges[fEdges[i]];

        if
        (
            (e[0] == v0 && e[1] == v1)
         || (e[0] == v1 && e[1] == v0)
        )
        {
            return fEdges[i];
        }
    }

    return -1;
}


Foam::label Foam::cellCuts::loopFace
(
    const label cellI,
    const labelList& loop
) const
{
    const cell& cFaces = mesh().cells()[cellI];

    forAll(cFaces, cFaceI)
    {
        label faceI = cFaces[cFaceI];

        const labelList& fEdges = mesh().faceEdges()[faceI];
        const face& f = mesh().faces()[faceI];

        bool allOnFace = true;

        forAll(loop, i)
        {
            label cut = loop[i];

            if (isEdge(cut))
            {
                if (findIndex(fEdges, getEdge(cut)) == -1)
                {
                    // Edge not on face. Skip face.
                    allOnFace = false;
                    break;
                }
            }
            else
            {
                if (findIndex(f, getVertex(cut)) == -1)
                {
                    // Vertex not on face. Skip face.
                    allOnFace = false;
                    break;
                }
            }
        }

        if (allOnFace)
        {
            // Found face where all elements of loop are on the face.
            return faceI;
        }
    }
    return -1;
}


bool Foam::cellCuts::walkPoint
(
    const label cellI,
    const label startCut,

    const label exclude0,
    const label exclude1,

    const label otherCut,

    label& nVisited,
    labelList& visited
) const
{
    label vertI = getVertex(otherCut);

    const labelList& pFaces = mesh().pointFaces()[vertI];

    forAll(pFaces, pFaceI)
    {
        label otherFaceI = pFaces[pFaceI];

        if
        (
            otherFaceI != exclude0
         && otherFaceI != exclude1
         && meshTools::faceOnCell(mesh(), cellI, otherFaceI)
        )
        {
            label oldNVisited = nVisited;

            bool foundLoop =
                walkCell
                (
                    cellI,
                    startCut,
                    otherFaceI,
                    otherCut,
                    nVisited,
                    visited
                );

            if (foundLoop)
            {
                return true;
            }

            // No success. Restore state and continue
            nVisited = oldNVisited;
        }
    }
    return false;
}


bool Foam::cellCuts::crossEdge
(
    const label cellI,
    const label startCut,
    const label faceI,
    const label otherCut,

    label& nVisited,
    labelList& visited
) const
{
    // Cross edge to other face
    label edgeI = getEdge(otherCut);

    label otherFaceI = meshTools::otherFace(mesh(), cellI, faceI, edgeI);

    // Store old state
    label oldNVisited = nVisited;

    // Recurse to otherCut
    bool foundLoop =
        walkCell
        (
            cellI,
            startCut,
            otherFaceI,
            otherCut,
            nVisited,
            visited
        );

    if (foundLoop)
    {
        return true;
    }
    else
    {
        // No success. Restore state (i.e. backtrack)
        nVisited = oldNVisited;

        return false;
    }
}


bool Foam::cellCuts::addCut
(
    const label cellI,
    const label cut,
    label& nVisited,
    labelList& visited
) const
{
    // Check for duplicate cuts.
    if (findPartIndex(visited, nVisited, cut) != -1)
    {
        // Truncate (copy of) visited for ease of printing.
        labelList truncVisited(visited);
        truncVisited.setSize(nVisited);

        Pout<< "For cell " << cellI << " : trying to add duplicate cut " << cut;
        labelList cuts(1, cut);
        writeCuts(Pout, cuts, loopWeights(cuts));

        Pout<< " to path:";
        writeCuts(Pout, truncVisited, loopWeights(truncVisited));
        Pout<< endl;

        return false;
    }
    else
    {
        visited[nVisited++] = cut;

        return true;
    }
}


bool Foam::cellCuts::walkFace
(
    const label cellI,
    const label startCut,
    const label faceI,
    const label cut,

    label& lastCut,
    label& beforeLastCut,
    label& nVisited,
    labelList& visited
) const
{
    const labelList& fCuts = faceCuts()[faceI];

    if (fCuts.size() < 2)
    {
        return false;
    }

    // Easy case : two cuts.
    if (fCuts.size() == 2)
    {
        if (fCuts[0] == cut)
        {
            if (!addCut(cellI, cut, nVisited, visited))
            {
                return false;
            }

            beforeLastCut = cut;
            lastCut = fCuts[1];

            return true;
        }
        else
        {
            if (!addCut(cellI, cut, nVisited, visited))
            {
                return false;
            }

            beforeLastCut = cut;
            lastCut = fCuts[0];

            return true;
        }
    }

    // Harder case: more than 2 cuts on face.
    // Should be start or end of string of cuts. Store all but last cut.
    if (fCuts[0] == cut)
    {
        // Walk forward
        for (label i = 0; i < fCuts.size()-1; i++)
        {
            if (!addCut(cellI, fCuts[i], nVisited, visited))
            {
                return false;
            }
        }
        beforeLastCut = fCuts[fCuts.size()-2];
        lastCut = fCuts[fCuts.size()-1];
    }
    else if (fCuts[fCuts.size()-1] == cut)
    {
        for (label i = fCuts.size()-1; i >= 1; --i)
        {
            if (!addCut(cellI, fCuts[i], nVisited, visited))
            {
                return false;
            }
        }
        beforeLastCut = fCuts[1];
        lastCut = fCuts[0];
    }
    else
    {
        WarningInFunction
            << "In middle of cut. cell:" << cellI << " face:" << faceI
            << " cuts:" << fCuts << " current cut:" << cut << endl;

        return false;
    }

    return true;
}



bool Foam::cellCuts::walkCell
(
    const label cellI,
    const label startCut,
    const label faceI,
    const label cut,

    label& nVisited,
    labelList& visited
) const
{
    // Walk across face, storing cuts. Return last two cuts visited.
    label lastCut = -1;
    label beforeLastCut = -1;


    if (debug & 2)
    {
        Pout<< "For cell:" << cellI << " walked across face " << faceI
            << " from cut ";
        labelList cuts(1, cut);
        writeCuts(Pout, cuts, loopWeights(cuts));
        Pout<< endl;
    }

    bool validWalk = walkFace
    (
        cellI,
        startCut,
        faceI,
        cut,

        lastCut,
        beforeLastCut,
        nVisited,
        visited
    );

    if (!validWalk)
    {
        return false;
    }

    if (debug & 2)
    {
        Pout<< "    to last cut ";
        labelList cuts(1, lastCut);
        writeCuts(Pout, cuts, loopWeights(cuts));
        Pout<< endl;
    }

    // Check if starting point reached.
    if (lastCut == startCut)
    {
        if (nVisited >= 3)
        {
            if (debug & 2)
            {
                // Truncate (copy of) visited for ease of printing.
                labelList truncVisited(visited);
                truncVisited.setSize(nVisited);

                Pout<< "For cell " << cellI << " : found closed path:";
                writeCuts(Pout, truncVisited, loopWeights(truncVisited));
                Pout<< " closed by " << lastCut << endl;
            }

            return true;
        }
        else
        {
            // Loop but too short. Probably folded back on itself.
            return false;
        }
    }


    // Check last cut and one before that to determine type.

    // From beforeLastCut to lastCut:
    //  - from edge to edge
    //      (- always not along existing edge)
    //  - from edge to vertex
    //      - not along existing edge
    //      - along existing edge: not allowed?
    //  - from vertex to edge
    //      - not along existing edge
    //      - along existing edge. Not allowed. See above.
    //  - from vertex to vertex
    //      - not along existing edge
    //      - along existing edge

    if (isEdge(beforeLastCut))
    {
        if (isEdge(lastCut))
        {
            // beforeLastCut=edge, lastCut=edge.

            // Cross lastCut (=edge) into face not faceI
            return crossEdge
            (
                cellI,
                startCut,
                faceI,
                lastCut,
                nVisited,
                visited
            );
        }
        else
        {
            // beforeLastCut=edge, lastCut=vertex.

            // Go from lastCut to all connected faces (except faceI)
            return walkPoint
            (
                cellI,
                startCut,
                faceI,          // exclude0: don't cross back on current face
                -1,             // exclude1
                lastCut,
                nVisited,
                visited
            );
        }
    }
    else
    {
        if (isEdge(lastCut))
        {
            // beforeLastCut=vertex, lastCut=edge.
            return crossEdge
            (
                cellI,
                startCut,
                faceI,
                lastCut,
                nVisited,
                visited
            );
        }
        else
        {
            // beforeLastCut=vertex, lastCut=vertex. Check if along existing
            // edge.
            label edgeI =
                findEdge
                (
                    faceI,
                    getVertex(beforeLastCut),
                    getVertex(lastCut)
                );

            if (edgeI != -1)
            {
                // Cut along existing edge. So is in fact on two faces.
                // Get faces on both sides of the edge to make
                // sure we dont fold back on to those.

                label f0, f1;
                meshTools::getEdgeFaces(mesh(), cellI, edgeI, f0, f1);

                return walkPoint
                (
                    cellI,
                    startCut,
                    f0,
                    f1,
                    lastCut,
                    nVisited,
                    visited
                );

            }
            else
            {
                // Cross cut across face.
                return walkPoint
                (
                    cellI,
                    startCut,
                    faceI,      // face to exclude
                    -1,         // face to exclude
                    lastCut,
                    nVisited,
                    visited
                );
            }
        }
    }
}


void Foam::cellCuts::calcCellLoops(const labelList& cutCells)
{
    // Determine for every cut cell the loop (= face) it is cut by. Done by
    // starting from a cut edge or cut vertex and walking across faces, from
    // cut to cut, until starting cut hit.
    // If multiple loops are possible across a cell circumference takes the
    // first one found.

    // Calculate cuts per face.
    const labelListList& allFaceCuts = faceCuts();

    // Per cell the number of faces with valid cuts. Is used as quick
    // rejection to see if cell can be cut.
    labelList nCutFaces(mesh().nCells(), 0);

    forAll(allFaceCuts, faceI)
    {
        const labelList& fCuts = allFaceCuts[faceI];

        if (fCuts.size() == mesh().faces()[faceI].size())
        {
            // Too many cuts on face. WalkCell would get very upset so disable.
            nCutFaces[mesh().faceOwner()[faceI]] = labelMin;

            if (mesh().isInternalFace(faceI))
            {
                nCutFaces[mesh().faceNeighbour()[faceI]] = labelMin;
            }
        }
        else if (fCuts.size() >= 2)
        {
            // Could be valid cut. Update count for owner and neighbour.
            nCutFaces[mesh().faceOwner()[faceI]]++;

            if (mesh().isInternalFace(faceI))
            {
                nCutFaces[mesh().faceNeighbour()[faceI]]++;
            }
        }
    }


    // Stack of visited cuts (nVisited used as stack pointer)
    // Size big enough.
    labelList visited(mesh().nPoints());

    forAll(cutCells, i)
    {
        label cellI = cutCells[i];

        bool validLoop = false;

        // Quick rejection: has enough faces that are cut?
        if (nCutFaces[cellI] >= 3)
        {
            const labelList& cFaces = mesh().cells()[cellI];

            if (debug & 2)
            {
                Pout<< "cell:" << cellI << " cut faces:" << endl;
                forAll(cFaces, i)
                {
                    label faceI = cFaces[i];
                    const labelList& fCuts = allFaceCuts[faceI];

                    Pout<< "    face:" << faceI << " cuts:";
                    writeCuts(Pout, fCuts, loopWeights(fCuts));
                    Pout<< endl;
                }
            }

            label nVisited = 0;

            // Determine the first cut face to start walking from.
            forAll(cFaces, cFaceI)
            {
                label faceI = cFaces[cFaceI];

                const labelList& fCuts = allFaceCuts[faceI];

                // Take first or last cut of multiple on face.
                // Note that in calcFaceCuts
                // we have already made sure this is the start or end of a
                // string of cuts.
                if (fCuts.size() >= 2)
                {
                    // Try walking from start of fCuts.
                    nVisited = 0;

                    if (debug & 2)
                    {
                        Pout<< "cell:" << cellI
                            << " start walk at face:" << faceI
                            << " cut:";
                        labelList cuts(1, fCuts[0]);
                        writeCuts(Pout, cuts, loopWeights(cuts));
                        Pout<< endl;
                    }

                    validLoop =
                        walkCell
                        (
                            cellI,
                            fCuts[0],
                            faceI,
                            fCuts[0],

                            nVisited,
                            visited
                        );

                    if (validLoop)
                    {
                        break;
                    }

                    // No need to try and walk from end of cuts since
                    // all paths already tried by walkCell.
                }
            }

            if (validLoop)
            {
                // Copy nVisited elements out of visited (since visited is
                // never truncated for efficiency reasons)

                labelList& loop = cellLoops_[cellI];

                loop.setSize(nVisited);

                for (label i = 0; i < nVisited; i++)
                {
                    loop[i] = visited[i];
                }
            }
            else
            {
                // Invalid loop. Leave cellLoops_[cellI] zero size which
                // flags this.
                Pout<< "calcCellLoops(const labelList&) : did not find valid"
                    << " loop for cell " << cellI << endl;
                // Dump cell and cuts on cell.
                writeUncutOBJ(".", cellI);

                cellLoops_[cellI].setSize(0);
            }
        }
        else
        {
            //Pout<< "calcCellLoops(const labelList&) : did not find valid"
            //    << " loop for cell " << cellI << " since not enough cut faces"
            //    << endl;
            cellLoops_[cellI].setSize(0);
        }
    }
}


void Foam::cellCuts::walkEdges
(
    const label cellI,
    const label pointI,
    const label status,

    Map<label>& edgeStatus,
    Map<label>& pointStatus
) const
{
    if (pointStatus.insert(pointI, status))
    {
        // First visit to pointI

        const labelList& pEdges = mesh().pointEdges()[pointI];

        forAll(pEdges, pEdgeI)
        {
            label edgeI = pEdges[pEdgeI];

            if
            (
                meshTools::edgeOnCell(mesh(), cellI, edgeI)
             && edgeStatus.insert(edgeI, status)
            )
            {
                // First visit to edgeI so recurse.

                label v2 = mesh().edges()[edgeI].otherVertex(pointI);

                walkEdges(cellI, v2, status, edgeStatus, pointStatus);
            }
        }
    }
}


Foam::labelList Foam::cellCuts::nonAnchorPoints
(
    const labelList& cellPoints,
    const labelList& anchorPoints,
    const labelList& loop
) const
{
    labelList newElems(cellPoints.size());
    label newElemI = 0;

    forAll(cellPoints, i)
    {
        label pointI = cellPoints[i];

        if
        (
            findIndex(anchorPoints, pointI) == -1
         && findIndex(loop, vertToEVert(pointI)) == -1
        )
        {
            newElems[newElemI++] = pointI;
        }
    }

    newElems.setSize(newElemI);

    return newElems;
}


bool Foam::cellCuts::loopAnchorConsistent
(
    const label cellI,
    const pointField& loopPts,
    const labelList& anchorPoints
) const
{
    // Create identity face for ease of calculation of normal etc.
    face f(identity(loopPts.size()));

    vector normal = f.normal(loopPts);
    point ctr = f.centre(loopPts);


    // Get average position of anchor points.
    vector avg(vector::zero);

    forAll(anchorPoints, ptI)
    {
        avg += mesh().points()[anchorPoints[ptI]];
    }
    avg /= anchorPoints.size();


    if (((avg - ctr) & normal) > 0)
    {
        return true;
    }
    else
    {
        return false;
    }
}


bool Foam::cellCuts::calcAnchors
(
    const label cellI,
    const labelList& loop,
    const pointField& loopPts,

    labelList& anchorPoints
) const
{
    const edgeList& edges = mesh().edges();

    const labelList& cPoints = mesh().cellPoints()[cellI];
    const labelList& cEdges = mesh().cellEdges()[cellI];
    const cell& cFaces = mesh().cells()[cellI];

    // Points on loop

    // Status of point:
    // 0 - on loop
    // 1 - point set 1
    // 2 - point set 2
    Map<label> pointStatus(2*cPoints.size());
    Map<label> edgeStatus(2*cEdges.size());

    // Mark loop vertices
    forAll(loop, i)
    {
        label cut = loop[i];

        if (isEdge(cut))
        {
            edgeStatus.insert(getEdge(cut), 0);
        }
        else
        {
            pointStatus.insert(getVertex(cut), 0);
        }
    }
    // Since edges between two cut vertices have not been marked, mark them
    // explicitly
    forAll(cEdges, i)
    {
        label edgeI = cEdges[i];
        const edge& e = edges[edgeI];

        if (pointStatus.found(e[0]) && pointStatus.found(e[1]))
        {
            edgeStatus.insert(edgeI, 0);
        }
    }


    // Find uncut starting vertex
    label uncutIndex = firstUnique(cPoints, pointStatus);

    if (uncutIndex == -1)
    {
        WarningInFunction
            << "Invalid loop " << loop << " for cell " << cellI << endl
            << "Can not find point on cell which is not cut by loop."
            << endl;

        writeOBJ(".", cellI, loopPts, labelList(0));

        return false;
    }

    // Walk unset vertices and edges and mark with 1 in pointStatus, edgeStatus
    walkEdges(cellI, cPoints[uncutIndex], 1, edgeStatus, pointStatus);

    // Find new uncut starting vertex
    uncutIndex = firstUnique(cPoints, pointStatus);

    if (uncutIndex == -1)
    {
        // All vertices either in loop or in anchor. So split is along single
        // face.
        WarningInFunction
            << "Invalid loop " << loop << " for cell " << cellI << endl
            << "All vertices of cell are either in loop or in anchor set"
            << endl;

        writeOBJ(".", cellI, loopPts, labelList(0));

        return false;
    }

    // Walk unset vertices and edges and mark with 2. These are the
    // pointset 2.
    walkEdges(cellI, cPoints[uncutIndex], 2, edgeStatus, pointStatus);

    // Collect both sets in lists.
    DynamicList<label> connectedPoints(cPoints.size());
    DynamicList<label> otherPoints(cPoints.size());

    forAllConstIter(Map<label>, pointStatus, iter)
    {
        if (iter() == 1)
        {
            connectedPoints.append(iter.key());
        }
        else if (iter() == 2)
        {
            otherPoints.append(iter.key());
        }
    }
    connectedPoints.shrink();
    otherPoints.shrink();

    // Check that all points have been used.
    uncutIndex = firstUnique(cPoints, pointStatus);

    if (uncutIndex != -1)
    {
        WarningInFunction
            << "Invalid loop " << loop << " for cell " << cellI
            << " since it splits the cell into more than two cells" << endl;

        writeOBJ(".", cellI, loopPts, connectedPoints);

        return false;
    }


    // Check that both parts (connectedPoints, otherPoints) have enough faces.
    labelHashSet connectedFaces(2*cFaces.size());
    labelHashSet otherFaces(2*cFaces.size());

    forAllConstIter(Map<label>, pointStatus, iter)
    {
        label pointI = iter.key();

        const labelList& pFaces = mesh().pointFaces()[pointI];

        if (iter() == 1)
        {
            forAll(pFaces, pFaceI)
            {
                if (meshTools::faceOnCell(mesh(), cellI, pFaces[pFaceI]))
                {
                    connectedFaces.insert(pFaces[pFaceI]);
                }
            }
        }
        else if (iter() == 2)
        {
            forAll(pFaces, pFaceI)
            {
                if (meshTools::faceOnCell(mesh(), cellI, pFaces[pFaceI]))
                {
                    otherFaces.insert(pFaces[pFaceI]);
                }
            }
        }
    }

    if (connectedFaces.size() < 3)
    {
        WarningInFunction
            << "Invalid loop " << loop << " for cell " << cellI
            << " since would have too few faces on one side." << nl
            << "All faces:" << cFaces << endl;

        writeOBJ(".", cellI, loopPts, connectedPoints);

        return false;
    }

    if (otherFaces.size() < 3)
    {
        WarningInFunction
            << "Invalid loop " << loop << " for cell " << cellI
            << " since would have too few faces on one side." << nl
            << "All faces:" << cFaces << endl;

        writeOBJ(".", cellI, loopPts, otherPoints);

        return false;
    }


    // Check that faces are split into two regions and not more.
    // When walking across the face the only transition of pointStatus is
    // from set1 to loop to set2 or back. Not allowed is from set1 to loop to
    // set1.
    {
        forAll(cFaces, i)
        {
            label faceI = cFaces[i];

            const face& f = mesh().faces()[faceI];

            bool hasSet1 = false;
            bool hasSet2 = false;

            label prevStat = pointStatus[f[0]];

            forAll(f, fp)
            {
                label v0 = f[fp];
                label pStat = pointStatus[v0];

                if (pStat == prevStat)
                {
                }
                else if (pStat == 0)
                {
                    // Loop.
                }
                else if (pStat == 1)
                {
                    if (hasSet1)
                    {
                        // Second occurence of set1.
                        WarningInFunction
                            << "Invalid loop " << loop << " for cell " << cellI
                            << " since face " << f << " would be split into"
                            << " more than two faces" << endl;

                        writeOBJ(".", cellI, loopPts, otherPoints);

                        return false;
                    }

                    hasSet1 = true;
                }
                else if (pStat == 2)
                {
                    if (hasSet2)
                    {
                        // Second occurence of set1.
                        WarningInFunction
                            << "Invalid loop " << loop << " for cell " << cellI
                            << " since face " << f << " would be split into"
                            << " more than two faces" << endl;

                        writeOBJ(".", cellI, loopPts, otherPoints);

                        return false;
                    }

                    hasSet2 = true;
                }
                else
                {
                    FatalErrorInFunction
                        << abort(FatalError);
                }

                prevStat = pStat;


                label v1 = f.nextLabel(fp);
                label edgeI = findEdge(faceI, v0, v1);

                label eStat = edgeStatus[edgeI];

                if (eStat == prevStat)
                {
                }
                else if (eStat == 0)
                {
                    // Loop.
                }
                else if (eStat == 1)
                {
                    if (hasSet1)
                    {
                        // Second occurence of set1.
                        WarningInFunction
                            << "Invalid loop " << loop << " for cell " << cellI
                            << " since face " << f << " would be split into"
                            << " more than two faces" << endl;

                        writeOBJ(".", cellI, loopPts, otherPoints);

                        return false;
                    }

                    hasSet1 = true;
                }
                else if (eStat == 2)
                {
                    if (hasSet2)
                    {
                        // Second occurence of set1.
                        WarningInFunction
                            << "Invalid loop " << loop << " for cell " << cellI
                            << " since face " << f << " would be split into"
                            << " more than two faces" << endl;

                        writeOBJ(".", cellI, loopPts, otherPoints);

                        return false;
                    }

                    hasSet2 = true;
                }
                prevStat = eStat;
            }
        }
    }




    // Check which one of point sets to use.
    bool loopOk = loopAnchorConsistent(cellI, loopPts, connectedPoints);

    //if (debug)
    {
        // Additional check: are non-anchor points on other side?
        bool otherLoopOk = loopAnchorConsistent(cellI, loopPts, otherPoints);

        if (loopOk == otherLoopOk)
        {
            // Both sets of points are supposedly on the same side as the
            // loop normal. Oops.

            WarningInFunction
                << " For cell:" << cellI
                << " achorpoints and nonanchorpoints are geometrically"
                << " on same side!" << endl
                << "cellPoints:" << cPoints << endl
                << "loop:" << loop << endl
                << "anchors:" << connectedPoints << endl
                << "otherPoints:" << otherPoints << endl;

            writeOBJ(".", cellI, loopPts, connectedPoints);
        }
    }

    if (loopOk)
    {
        // connectedPoints on 'outside' of loop so these are anchor points
        anchorPoints.transfer(connectedPoints);
        connectedPoints.clear();
    }
    else
    {
        anchorPoints.transfer(otherPoints);
        otherPoints.clear();
    }
    return true;
}


Foam::pointField Foam::cellCuts::loopPoints
(
    const labelList& loop,
    const scalarField& loopWeights
) const
{
    pointField loopPts(loop.size());

    forAll(loop, fp)
    {
        loopPts[fp] = coord(loop[fp], loopWeights[fp]);
    }
    return loopPts;
}


Foam::scalarField Foam::cellCuts::loopWeights(const labelList& loop) const
{
    scalarField weights(loop.size());

    forAll(loop, fp)
    {
        label cut = loop[fp];

        if (isEdge(cut))
        {
            label edgeI = getEdge(cut);

            weights[fp] = edgeWeight_[edgeI];
        }
        else
        {
            weights[fp] = -GREAT;
        }
    }
    return weights;
}


bool Foam::cellCuts::validEdgeLoop
(
    const labelList& loop,
    const scalarField& loopWeights
) const
{
    forAll(loop, fp)
    {
        label cut = loop[fp];

        if (isEdge(cut))
        {
            label edgeI = getEdge(cut);

            // Check: cut compatible only if can be snapped to existing one.
            if (edgeIsCut_[edgeI])
            {
                scalar edgeLen = mesh().edges()[edgeI].mag(mesh().points());

                if
                (
                    mag(loopWeights[fp] - edgeWeight_[edgeI])
                  > geomCellLooper::snapTol()*edgeLen
                )
                {
                    // Positions differ too much->would create two cuts.
                    return false;
                }
            }
        }
    }
    return true;
}


Foam::label Foam::cellCuts::countFaceCuts
(
    const label faceI,
    const labelList& loop
) const
{
    // Includes cuts through vertices and through edges.
    // Assumes that if edge is cut both in edgeIsCut and in loop that the
    // position of the cut is the same.

    label nCuts = 0;

    // Count cut vertices
    const face& f = mesh().faces()[faceI];

    forAll(f, fp)
    {
        label vertI = f[fp];

        // Vertex already cut or mentioned in current loop.
        if
        (
            pointIsCut_[vertI]
         || (findIndex(loop, vertToEVert(vertI)) != -1)
        )
        {
            nCuts++;
        }
    }

    // Count cut edges.
    const labelList& fEdges = mesh().faceEdges()[faceI];

    forAll(fEdges, fEdgeI)
    {
        label edgeI = fEdges[fEdgeI];

        // Edge already cut or mentioned in current loop.
        if
        (
            edgeIsCut_[edgeI]
         || (findIndex(loop, edgeToEVert(edgeI)) != -1)
        )
        {
            nCuts++;
        }
    }

    return nCuts;
}


bool Foam::cellCuts::conservativeValidLoop
(
    const label cellI,
    const labelList& loop
) const
{

    if (loop.size() < 2)
    {
        return false;
    }

    forAll(loop, cutI)
    {
        if (isEdge(loop[cutI]))
        {
            label edgeI = getEdge(loop[cutI]);

            if (edgeIsCut_[edgeI])
            {
                // edge compatibility already checked.
            }
            else
            {
                // Quick rejection: vertices of edge itself cannot be cut.
                const edge& e = mesh().edges()[edgeI];

                if (pointIsCut_[e.start()] || pointIsCut_[e.end()])
                {
                    return false;
                }


                // Check faces using this edge
                const labelList& eFaces = mesh().edgeFaces()[edgeI];

                forAll(eFaces, eFaceI)
                {
                    label nCuts = countFaceCuts(eFaces[eFaceI], loop);

                    if (nCuts > 2)
                    {
                        return false;
                    }
                }
            }
        }
        else
        {
            // Vertex cut

            label vertI = getVertex(loop[cutI]);

            if (!pointIsCut_[vertI])
            {
                // New cut through vertex.

                // Check edges using vertex.
                const labelList& pEdges = mesh().pointEdges()[vertI];

                forAll(pEdges, pEdgeI)
                {
                    label edgeI = pEdges[pEdgeI];

                    if (edgeIsCut_[edgeI])
                    {
                        return false;
                    }
                }

                // Check faces using vertex.
                const labelList& pFaces = mesh().pointFaces()[vertI];

                forAll(pFaces, pFaceI)
                {
                    label nCuts = countFaceCuts(pFaces[pFaceI], loop);

                    if (nCuts > 2)
                    {
                        return false;
                    }
                }
            }
        }
    }

    // All ok.
    return true;
}


bool Foam::cellCuts::validLoop
(
    const label cellI,
    const labelList& loop,
    const scalarField& loopWeights,

    Map<edge>& newFaceSplitCut,
    labelList& anchorPoints
) const
{
    // Determine compatibility of loop with existing cut pattern. Does not use
    // derived cut-addressing (faceCuts), only pointIsCut, edgeIsCut.
    // Adds any cross-cuts found to newFaceSplitCut and sets cell points on
    // one side of the loop in anchorPoints.

    if (loop.size() < 2)
    {
        return false;
    }

    if (debug & 4)
    {
        // Allow as fallback the 'old' loop checking where only a single
        // cut per face is allowed.
        if (!conservativeValidLoop(cellI, loop))
        {
            Info << "Invalid conservative loop: " << loop << endl;
            return  false;
        }
    }

    forAll(loop, fp)
    {
        label cut = loop[fp];
        label nextCut = loop[(fp+1) % loop.size()];

        // Label (if any) of face cut (so cut not along existing edge)
        label meshFaceI = -1;

        if (isEdge(cut))
        {
            label edgeI = getEdge(cut);

            // Look one cut ahead to find if it is along existing edge.

            if (isEdge(nextCut))
            {
                // From edge to edge -> cross cut
                label nextEdgeI = getEdge(nextCut);

                // Find face and mark as to be split.
                meshFaceI = edgeEdgeToFace(cellI, edgeI, nextEdgeI);

                if (meshFaceI == -1)
                {
                    // Can't find face using both cut edges.
                    return false;
                }
            }
            else
            {
                // From edge to vertex -> cross cut only if no existing edge.

                label nextVertI = getVertex(nextCut);
                const edge& e = mesh().edges()[edgeI];

                if (e.start() != nextVertI && e.end() != nextVertI)
                {
                    // New edge. Find face and mark as to be split.
                    meshFaceI = edgeVertexToFace(cellI, edgeI, nextVertI);

                    if (meshFaceI == -1)
                    {
                        // Can't find face. Ilegal.
                        return false;
                    }
                }
            }
        }
        else
        {
            // Cut is vertex.
            label vertI = getVertex(cut);

            if (isEdge(nextCut))
            {
                // From vertex to edge -> cross cut only if no existing edge.
                label nextEdgeI = getEdge(nextCut);

                const edge& nextE = mesh().edges()[nextEdgeI];

                if (nextE.start() != vertI && nextE.end() != vertI)
                {
                    // Should be cross cut. Find face.
                    meshFaceI = edgeVertexToFace(cellI, nextEdgeI, vertI);

                    if (meshFaceI == -1)
                    {
                        return false;
                    }
                }
            }
            else
            {
                // From vertex to vertex -> cross cut only if no existing edge.
                label nextVertI = getVertex(nextCut);

                if (meshTools::findEdge(mesh(), vertI, nextVertI) == -1)
                {
                    // New edge. Find face.
                    meshFaceI = vertexVertexToFace(cellI, vertI, nextVertI);

                    if (meshFaceI == -1)
                    {
                        return false;
                    }
                }
            }
        }

        if (meshFaceI != -1)
        {
            // meshFaceI is cut across along cut-nextCut (so not along existing
            // edge). Check if this is compatible with existing pattern.
            edge cutEdge(cut, nextCut);

            Map<edge>::const_iterator iter = faceSplitCut_.find(meshFaceI);

            if (iter == faceSplitCut_.end())
            {
                // Face not yet cut so insert.
                newFaceSplitCut.insert(meshFaceI, cutEdge);
            }
            else
            {
                // Face already cut. Ok if same edge.
                if (iter() != cutEdge)
                {
                    return false;
                }
            }
        }
    }

    // Is there a face on which all cuts are?
    label faceContainingLoop = loopFace(cellI, loop);

    if (faceContainingLoop != -1)
    {
        WarningInFunction
            << "Found loop on cell " << cellI << " with all points"
            << " on face " << faceContainingLoop << endl;

        //writeOBJ(".", cellI, loopPoints(loop, loopWeights), labelList(0));

        return false;
    }

    // Calculate anchor points
    // Final success is determined by whether anchor points can be determined.
    return calcAnchors
    (
        cellI,
        loop,
        loopPoints(loop, loopWeights),
        anchorPoints
    );
}


void Foam::cellCuts::setFromCellLoops()
{
    // 'Uncut' edges/vertices that are not used in loops.
    pointIsCut_ = false;

    edgeIsCut_ = false;

    faceSplitCut_.clear();

    forAll(cellLoops_, cellI)
    {
        const labelList& loop = cellLoops_[cellI];

        if (loop.size())
        {
            // Storage for cross-face cuts
            Map<edge> faceSplitCuts(loop.size());
            // Storage for points on one side of cell.
            labelList anchorPoints;

            if
            (
               !validLoop
                (
                    cellI,
                    loop,
                    loopWeights(loop),
                    faceSplitCuts,
                    anchorPoints
                )
            )
            {
                //writeOBJ(".", cellI, loopPoints(cellI), anchorPoints);

                WarningInFunction
                    << "Illegal loop " << loop
                    << " when recreating cut-addressing"
                    << " from existing cellLoops for cell " << cellI
                    << endl;

                cellLoops_[cellI].setSize(0);
                cellAnchorPoints_[cellI].setSize(0);
            }
            else
            {
                // Copy anchor points.
                cellAnchorPoints_[cellI].transfer(anchorPoints);

                // Copy faceSplitCuts into overall faceSplit info.
                forAllConstIter(Map<edge>, faceSplitCuts, iter)
                {
                    faceSplitCut_.insert(iter.key(), iter());
                }

                // Update edgeIsCut, pointIsCut information
                forAll(loop, cutI)
                {
                    label cut = loop[cutI];

                    if (isEdge(cut))
                    {
                        edgeIsCut_[getEdge(cut)] = true;
                    }
                    else
                    {
                        pointIsCut_[getVertex(cut)] = true;
                    }
                }
            }
        }
    }

    // Reset edge weights
    forAll(edgeIsCut_, edgeI)
    {
        if (!edgeIsCut_[edgeI])
        {
            // Weight not used. Set to illegal value to make any use fall over.
            edgeWeight_[edgeI] = -GREAT;
        }
    }
}


bool Foam::cellCuts::setFromCellLoop
(
    const label cellI,
    const labelList& loop,
    const scalarField& loopWeights
)
{
    // Update basic cut information from single cellLoop. Returns true if loop
    // was valid.

    // Dump loop for debugging.
    if (debug)
    {
        OFstream str("last_cell.obj");

        str<< "# edges of cell " << cellI << nl;

        meshTools::writeOBJ
        (
            str,
            mesh().cells(),
            mesh().faces(),
            mesh().points(),
            labelList(1, cellI)
        );


        OFstream loopStr("last_loop.obj");

        loopStr<< "# looppoints for cell " << cellI << nl;

        pointField pointsOfLoop = loopPoints(loop, loopWeights);

        forAll(pointsOfLoop, i)
        {
            meshTools::writeOBJ(loopStr, pointsOfLoop[i]);
        }

        str << 'l';

        forAll(pointsOfLoop, i)
        {
            str << ' ' << i + 1;
        }
        str << ' ' << 1 << nl;
    }

    bool okLoop = false;

    if (validEdgeLoop(loop, loopWeights))
    {
        // Storage for cuts across face
        Map<edge> faceSplitCuts(loop.size());
        // Storage for points on one side of cell
        labelList anchorPoints;

        okLoop =
            validLoop(cellI, loop, loopWeights, faceSplitCuts, anchorPoints);

        if (okLoop)
        {
            // Valid loop. Copy cellLoops and anchorPoints
            cellLoops_[cellI] = loop;
            cellAnchorPoints_[cellI].transfer(anchorPoints);

            // Copy split cuts
            forAllConstIter(Map<edge>, faceSplitCuts, iter)
            {
                faceSplitCut_.insert(iter.key(), iter());
            }


            // Update edgeIsCut, pointIsCut information
            forAll(loop, cutI)
            {
                label cut = loop[cutI];

                if (isEdge(cut))
                {
                    label edgeI = getEdge(cut);

                    edgeIsCut_[edgeI] = true;

                    edgeWeight_[edgeI] = loopWeights[cutI];
                }
                else
                {
                    label vertI = getVertex(cut);

                    pointIsCut_[vertI] = true;
                }
            }
        }
    }

    return okLoop;
}


void Foam::cellCuts::setFromCellLoops
(
    const labelList& cellLabels,
    const labelListList& cellLoops,
    const List<scalarField>& cellLoopWeights
)
{
    // 'Uncut' edges/vertices that are not used in loops.
    pointIsCut_ = false;

    edgeIsCut_ = false;

    forAll(cellLabels, cellLabelI)
    {
        label cellI = cellLabels[cellLabelI];

        const labelList& loop = cellLoops[cellLabelI];

        if (loop.size())
        {
            const scalarField& loopWeights = cellLoopWeights[cellLabelI];

            if (setFromCellLoop(cellI, loop, loopWeights))
            {
                // Valid loop. Call above will have upated all already.
            }
            else
            {
                // Clear cellLoops
                cellLoops_[cellI].setSize(0);
            }
        }
    }
}


void Foam::cellCuts::setFromCellCutter
(
    const cellLooper& cellCutter,
    const List<refineCell>& refCells
)
{
    // 'Uncut' edges/vertices that are not used in loops.
    pointIsCut_ = false;

    edgeIsCut_ = false;

    // storage for loop of cuts (cut vertices and/or cut edges)
    labelList cellLoop;
    scalarField cellLoopWeights;

    // For debugging purposes
    DynamicList<label> invalidCutCells(2);
    DynamicList<labelList> invalidCutLoops(2);
    DynamicList<scalarField> invalidCutLoopWeights(2);


    forAll(refCells, refCellI)
    {
        const refineCell& refCell = refCells[refCellI];

        label cellI = refCell.cellNo();

        const vector& refDir = refCell.direction();


        // Cut cell. Determines cellLoop and cellLoopWeights
        bool goodCut =
            cellCutter.cut
            (
                refDir,
                cellI,

                pointIsCut_,
                edgeIsCut_,
                edgeWeight_,

                cellLoop,
                cellLoopWeights
            );

        // Check whether edge refinement is on a per face basis compatible with
        // current pattern.
        if (goodCut)
        {
            if (setFromCellLoop(cellI, cellLoop, cellLoopWeights))
            {
                // Valid loop. Will have updated all info already.
            }
            else
            {
                cellLoops_[cellI].setSize(0);

                WarningInFunction
                    << "Found loop on cell " << cellI
                    << " that resulted in an unexpected bad cut."
                    << "    Suggestions:" << nl
                    << "      - Turn on the debug switch for 'cellCuts' to get"
                    << " geometry files that identify this cell." << nl
                    << "      - Also keep in mind to check the defined"
                    << " reference directions, as these are most likely the"
                    << " origin of the problem."
                    << nl << endl;

                // Discarded by validLoop
                if (debug)
                {
                    invalidCutCells.append(cellI);
                    invalidCutLoops.append(cellLoop);
                    invalidCutLoopWeights.append(cellLoopWeights);
                }
            }
        }
        else
        {
            // Clear cellLoops
            cellLoops_[cellI].setSize(0);
        }
    }

    if (debug && invalidCutCells.size())
    {
        invalidCutCells.shrink();
        invalidCutLoops.shrink();
        invalidCutLoopWeights.shrink();

        fileName cutsFile("invalidLoopCells.obj");

        Pout<< "cellCuts : writing inValidLoops cells to " << cutsFile << endl;

        OFstream cutsStream(cutsFile);

        meshTools::writeOBJ
        (
            cutsStream,
            mesh().cells(),
            mesh().faces(),
            mesh().points(),
            invalidCutCells
        );

        fileName loopsFile("invalidLoops.obj");

        Pout<< "cellCuts : writing inValidLoops loops to " << loopsFile << endl;

        OFstream loopsStream(loopsFile);

        label vertI = 0;

        forAll(invalidCutLoops, i)
        {
            writeOBJ
            (
                loopsStream,
                loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
                vertI
            );
        }
    }
}


void Foam::cellCuts::setFromCellCutter
(
    const cellLooper& cellCutter,
    const labelList& cellLabels,
    const List<plane>& cellCutPlanes
)
{
    // 'Uncut' edges/vertices that are not used in loops.
    pointIsCut_ = false;

    edgeIsCut_ = false;

    // storage for loop of cuts (cut vertices and/or cut edges)
    labelList cellLoop;
    scalarField cellLoopWeights;

    // For debugging purposes
    DynamicList<label> invalidCutCells(2);
    DynamicList<labelList> invalidCutLoops(2);
    DynamicList<scalarField> invalidCutLoopWeights(2);


    forAll(cellLabels, i)
    {
        label cellI = cellLabels[i];

        // Cut cell. Determines cellLoop and cellLoopWeights
        bool goodCut =
            cellCutter.cut
            (
                cellCutPlanes[i],
                cellI,

                pointIsCut_,
                edgeIsCut_,
                edgeWeight_,

                cellLoop,
                cellLoopWeights
            );

        // Check whether edge refinement is on a per face basis compatible with
        // current pattern.
        if (goodCut)
        {
            if (setFromCellLoop(cellI, cellLoop, cellLoopWeights))
            {
                // Valid loop. Will have updated all info already.
            }
            else
            {
                cellLoops_[cellI].setSize(0);

                // Discarded by validLoop
                if (debug)
                {
                    invalidCutCells.append(cellI);
                    invalidCutLoops.append(cellLoop);
                    invalidCutLoopWeights.append(cellLoopWeights);
                }
            }
        }
        else
        {
            // Clear cellLoops
            cellLoops_[cellI].setSize(0);
        }
    }

    if (debug && invalidCutCells.size())
    {
        invalidCutCells.shrink();
        invalidCutLoops.shrink();
        invalidCutLoopWeights.shrink();

        fileName cutsFile("invalidLoopCells.obj");

        Pout<< "cellCuts : writing inValidLoops cells to " << cutsFile << endl;

        OFstream cutsStream(cutsFile);

        meshTools::writeOBJ
        (
            cutsStream,
            mesh().cells(),
            mesh().faces(),
            mesh().points(),
            invalidCutCells
        );

        fileName loopsFile("invalidLoops.obj");

        Pout<< "cellCuts : writing inValidLoops loops to " << loopsFile << endl;

        OFstream loopsStream(loopsFile);

        label vertI = 0;

        forAll(invalidCutLoops, i)
        {
            writeOBJ
            (
                loopsStream,
                loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
                vertI
            );
        }
    }
}


void Foam::cellCuts::orientPlanesAndLoops()
{
    // Determine anchorPoints if not yet done by validLoop.
    forAll(cellLoops_, cellI)
    {
        const labelList& loop = cellLoops_[cellI];

        if (loop.size() && cellAnchorPoints_[cellI].empty())
        {
            // Leave anchor points empty if illegal loop.
            calcAnchors
            (
                cellI,
                loop,
                loopPoints(cellI),
                cellAnchorPoints_[cellI]
            );
        }
    }

    if (debug & 2)
    {
        Pout<< "cellAnchorPoints:" << endl;
    }
    forAll(cellAnchorPoints_, cellI)
    {
        if (cellLoops_[cellI].size())
        {
            if (cellAnchorPoints_[cellI].empty())
            {
                FatalErrorInFunction
                    << "No anchor points for cut cell " << cellI << endl
                    << "cellLoop:" << cellLoops_[cellI] << abort(FatalError);
            }

            if (debug & 2)
            {
                Pout<< "    cell:" << cellI << " anchored at "
                    << cellAnchorPoints_[cellI] << endl;
            }
        }
    }

    // Calculate number of valid cellLoops
    nLoops_ = 0;

    forAll(cellLoops_, cellI)
    {
        if (cellLoops_[cellI].size())
        {
            nLoops_++;
        }
    }
}


void Foam::cellCuts::calcLoopsAndAddressing(const labelList& cutCells)
{
    // Sanity check on weights
    forAll(edgeIsCut_, edgeI)
    {
        if (edgeIsCut_[edgeI])
        {
            scalar weight = edgeWeight_[edgeI];

            if (weight < 0 || weight > 1)
            {
                FatalErrorInFunction
                    << "Weight out of range [0,1]. Edge " << edgeI
                    << " verts:" << mesh().edges()[edgeI]
                    << " weight:" << weight << abort(FatalError);
            }
        }
        else
        {
            // Weight not used. Set to illegal value to make any use fall over.
            edgeWeight_[edgeI] = -GREAT;
        }
    }


    // Calculate faces that split cells in two
    calcCellLoops(cutCells);

    if (debug & 2)
    {
        Pout<< "-- cellLoops --" << endl;
        forAll(cellLoops_, cellI)
        {
            const labelList& loop = cellLoops_[cellI];

            if (loop.size())
            {
                Pout<< "cell:" << cellI << "  ";
                writeCuts(Pout, loop, loopWeights(loop));
                Pout<< endl;
            }
        }
    }

    // Redo basic cut information (pointIsCut, edgeIsCut, faceSplitCut)
    // using cellLoop only.
    setFromCellLoops();
}


void Foam::cellCuts::check() const
{
    // Check weights for unsnapped values
    forAll(edgeIsCut_, edgeI)
    {
        if (edgeIsCut_[edgeI])
        {
            if
            (
                edgeWeight_[edgeI] < geomCellLooper::snapTol()
             || edgeWeight_[edgeI] > (1 - geomCellLooper::snapTol())
            )
            {
                // Should have been snapped.
                WarningInFunction
                    << "edge:" << edgeI << " vertices:"
                    << mesh().edges()[edgeI]
                    << " weight:" << edgeWeight_[edgeI] << " should have been"
                    << " snapped to one of its endpoints"
                    << endl;
            }
        }
        else
        {
            if (edgeWeight_[edgeI] > - 1)
            {
                FatalErrorInFunction
                    << "edge:" << edgeI << " vertices:"
                    << mesh().edges()[edgeI]
                    << " weight:" << edgeWeight_[edgeI] << " is not cut but"
                    << " its weight is not marked invalid"
                    << abort(FatalError);
            }
        }
    }

    // Check that all elements of cellloop are registered
    forAll(cellLoops_, cellI)
    {
        const labelList& loop = cellLoops_[cellI];

        forAll(loop, i)
        {
            label cut = loop[i];

            if
            (
                (isEdge(cut) && !edgeIsCut_[getEdge(cut)])
             || (!isEdge(cut) && !pointIsCut_[getVertex(cut)])
            )
            {
                labelList cuts(1, cut);
                writeCuts(Pout, cuts, loopWeights(cuts));

                FatalErrorInFunction
                    << "cell:" << cellI << " loop:"
                    << loop
                    << " cut:" << cut << " is not marked as cut"
                    << abort(FatalError);
            }
        }
    }

    // Check that no elements of cell loop are anchor point.
    forAll(cellLoops_, cellI)
    {
        const labelList& anchors = cellAnchorPoints_[cellI];

        const labelList& loop = cellLoops_[cellI];

        if (loop.size() && anchors.empty())
        {
            FatalErrorInFunction
                << "cell:" << cellI << " loop:" << loop
                << " has no anchor points"
                << abort(FatalError);
        }


        forAll(loop, i)
        {
            label cut = loop[i];

            if
            (
                !isEdge(cut)
             && findIndex(anchors, getVertex(cut)) != -1
            )
            {
                FatalErrorInFunction
                    << "cell:" << cellI << " loop:" << loop
                    << " anchor points:" << anchors
                    << " anchor:" << getVertex(cut) << " is part of loop"
                    << abort(FatalError);
            }
        }
    }


    // Check that cut faces have a neighbour that is cut.
    forAllConstIter(Map<edge>, faceSplitCut_, iter)
    {
        label faceI = iter.key();

        if (mesh().isInternalFace(faceI))
        {
            label own = mesh().faceOwner()[faceI];
            label nei = mesh().faceNeighbour()[faceI];

            if (cellLoops_[own].empty() && cellLoops_[nei].empty())
            {
                FatalErrorInFunction
                    << "Internal face:" << faceI << " cut by " << iter()
                    << " has owner:" << own
                    << " and neighbour:" << nei
                    << " that are both uncut"
                    << abort(FatalError);
            }
        }
        else
        {
            label own = mesh().faceOwner()[faceI];

            if (cellLoops_[own].empty())
            {
                FatalErrorInFunction
                    << "Boundary face:" << faceI << " cut by " << iter()
                    << " has owner:" << own
                    << " that is uncut"
                    << abort(FatalError);
            }
        }
    }
}


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

Foam::cellCuts::cellCuts
(
    const polyMesh& mesh,
    const labelList& cutCells,
    const labelList& meshVerts,
    const labelList& meshEdges,
    const scalarField& meshEdgeWeights
)
:
    edgeVertex(mesh),
    pointIsCut_(expand(mesh.nPoints(), meshVerts)),
    edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
    edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
    faceCutsPtr_(NULL),
    faceSplitCut_(cutCells.size()),
    cellLoops_(mesh.nCells()),
    nLoops_(-1),
    cellAnchorPoints_(mesh.nCells())
{
    if (debug)
    {
        Pout<< "cellCuts : constructor from cut verts and edges" << endl;
    }

    calcLoopsAndAddressing(cutCells);

    // Calculate planes and flip cellLoops if necessary
    orientPlanesAndLoops();

    if (debug)
    {
        check();
    }

    clearOut();

    if (debug)
    {
        Pout<< "cellCuts : leaving constructor from cut verts and edges"
            << endl;
    }
}


Foam::cellCuts::cellCuts
(
    const polyMesh& mesh,
    const labelList& meshVerts,
    const labelList& meshEdges,
    const scalarField& meshEdgeWeights
)
:
    edgeVertex(mesh),
    pointIsCut_(expand(mesh.nPoints(), meshVerts)),
    edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
    edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
    faceCutsPtr_(NULL),
    faceSplitCut_(mesh.nFaces()/10 + 1),
    cellLoops_(mesh.nCells()),
    nLoops_(-1),
    cellAnchorPoints_(mesh.nCells())
{
    // Construct from pattern of cuts. Finds out itself which cells are cut.
    // (can go wrong if e.g. all neighbours of cell are refined)

    if (debug)
    {
        Pout<< "cellCuts : constructor from cellLoops" << endl;
    }

    calcLoopsAndAddressing(identity(mesh.nCells()));

    // Calculate planes and flip cellLoops if necessary
    orientPlanesAndLoops();

    if (debug)
    {
        check();
    }

    clearOut();

    if (debug)
    {
        Pout<< "cellCuts : leaving constructor from cellLoops" << endl;
    }
}


Foam::cellCuts::cellCuts
(
    const polyMesh& mesh,
    const labelList& cellLabels,
    const labelListList& cellLoops,
    const List<scalarField>& cellEdgeWeights
)
:
    edgeVertex(mesh),
    pointIsCut_(mesh.nPoints(), false),
    edgeIsCut_(mesh.nEdges(), false),
    edgeWeight_(mesh.nEdges(), -GREAT),
    faceCutsPtr_(NULL),
    faceSplitCut_(cellLabels.size()),
    cellLoops_(mesh.nCells()),
    nLoops_(-1),
    cellAnchorPoints_(mesh.nCells())
{
    if (debug)
    {
        Pout<< "cellCuts : constructor from cellLoops" << endl;
    }

    // Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
    // Makes sure cuts are consistent
    setFromCellLoops(cellLabels, cellLoops, cellEdgeWeights);

    // Calculate planes and flip cellLoops if necessary
    orientPlanesAndLoops();

    if (debug)
    {
        check();
    }

    clearOut();

    if (debug)
    {
        Pout<< "cellCuts : leaving constructor from cellLoops" << endl;
    }
}


Foam::cellCuts::cellCuts
(
    const polyMesh& mesh,
    const cellLooper& cellCutter,
    const List<refineCell>& refCells
)
:
    edgeVertex(mesh),
    pointIsCut_(mesh.nPoints(), false),
    edgeIsCut_(mesh.nEdges(), false),
    edgeWeight_(mesh.nEdges(), -GREAT),
    faceCutsPtr_(NULL),
    faceSplitCut_(refCells.size()),
    cellLoops_(mesh.nCells()),
    nLoops_(-1),
    cellAnchorPoints_(mesh.nCells())
{
    if (debug)
    {
        Pout<< "cellCuts : constructor from cellCutter" << endl;
    }

    // Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
    // Makes sure cuts are consistent
    setFromCellCutter(cellCutter, refCells);

    // Calculate planes and flip cellLoops if necessary
    orientPlanesAndLoops();

    if (debug)
    {
        check();
    }

    clearOut();

    if (debug)
    {
        Pout<< "cellCuts : leaving constructor from cellCutter" << endl;
    }
}


Foam::cellCuts::cellCuts
(
    const polyMesh& mesh,
    const cellLooper& cellCutter,
    const labelList& cellLabels,
    const List<plane>& cutPlanes
)
:
    edgeVertex(mesh),
    pointIsCut_(mesh.nPoints(), false),
    edgeIsCut_(mesh.nEdges(), false),
    edgeWeight_(mesh.nEdges(), -GREAT),
    faceCutsPtr_(NULL),
    faceSplitCut_(cellLabels.size()),
    cellLoops_(mesh.nCells()),
    nLoops_(-1),
    cellAnchorPoints_(mesh.nCells())
{
    if (debug)
    {
        Pout<< "cellCuts : constructor from cellCutter with prescribed plane"
            << endl;
    }

    // Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
    // Makes sure cuts are consistent
    setFromCellCutter(cellCutter, cellLabels, cutPlanes);

    // Calculate planes and flip cellLoops if necessary
    orientPlanesAndLoops();

    if (debug)
    {
        check();
    }

    clearOut();

    if (debug)
    {
        Pout<< "cellCuts : leaving constructor from cellCutter with prescribed"
            << " plane" << endl;
    }
}


Foam::cellCuts::cellCuts
(
    const polyMesh& mesh,
    const boolList& pointIsCut,
    const boolList& edgeIsCut,
    const scalarField& edgeWeight,
    const Map<edge>& faceSplitCut,
    const labelListList& cellLoops,
    const label nLoops,
    const labelListList& cellAnchorPoints
)
:
    edgeVertex(mesh),
    pointIsCut_(pointIsCut),
    edgeIsCut_(edgeIsCut),
    edgeWeight_(edgeWeight),
    faceCutsPtr_(NULL),
    faceSplitCut_(faceSplitCut),
    cellLoops_(cellLoops),
    nLoops_(nLoops),
    cellAnchorPoints_(cellAnchorPoints)
{
    if (debug)
    {
        Pout<< "cellCuts : constructor from components" << endl;
        Pout<< "cellCuts : leaving constructor from components" << endl;
    }
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::cellCuts::~cellCuts()
{
    clearOut();
}


void Foam::cellCuts::clearOut()
{
    deleteDemandDrivenData(faceCutsPtr_);
}


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

Foam::pointField Foam::cellCuts::loopPoints(const label cellI) const
{
    const labelList& loop = cellLoops_[cellI];

    pointField loopPts(loop.size());

    forAll(loop, fp)
    {
        label cut = loop[fp];

        if (isEdge(cut))
        {
            label edgeI = getEdge(cut);

            loopPts[fp] = coord(cut, edgeWeight_[edgeI]);
        }
        else
        {
            loopPts[fp] = coord(cut, -GREAT);
        }
    }
    return loopPts;
}


void Foam::cellCuts::flip(const label cellI)
{
    labelList& loop = cellLoops_[cellI];

    reverse(loop);

    // Reverse anchor point set.
    cellAnchorPoints_[cellI] =
        nonAnchorPoints
        (
            mesh().cellPoints()[cellI],
            cellAnchorPoints_[cellI],
            loop
        );
}


void Foam::cellCuts::flipLoopOnly(const label cellI)
{
    labelList& loop = cellLoops_[cellI];

    reverse(loop);
}


void Foam::cellCuts::writeOBJ
(
    Ostream& os,
    const pointField& loopPts,
    label& vertI
) const
{
    label startVertI = vertI;

    forAll(loopPts, fp)
    {
        const point& pt = loopPts[fp];

        os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;

        vertI++;
    }

    os  << 'f';
    forAll(loopPts, fp)
    {
        os  << ' ' << startVertI + fp + 1;
    }
    os  << endl;
}


void Foam::cellCuts::writeOBJ(Ostream& os) const
{
    label vertI = 0;

    forAll(cellLoops_, cellI)
    {
        writeOBJ(os, loopPoints(cellI), vertI);
    }
}


void Foam::cellCuts::writeCellOBJ(const fileName& dir, const label cellI) const
{
    const labelList& anchors = cellAnchorPoints_[cellI];

    writeOBJ(dir, cellI, loopPoints(cellI), anchors);
}


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

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 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::cellCuts

Description
    Description of cuts across cells.

    Description of cut is given as list of vertices and
    list of edges to be cut (and position on edge).

    Does some checking of correctness/non-overlapping of cuts. 2x2x2
    refinement has to be done in three passes since cuts can not overlap
    (would make addressing too complicated)

    Introduces concept of 'cut' which is either an existing vertex
    or a edge.

    Input can either be
    -# list of cut vertices and list of cut edges. Constructs cell
       circumference walks ('cellLoops').

    -# list of cell circumference walks. Will filter them so they
       don't overlap.

    -# cellWalker and list of cells to refine (refineCell). Constructs
       cellLoops and does B. cellWalker is class which can cut a single
       cell using a plane through the cell centre and in a certain normal
       direction

    CellCuts constructed from cellLoops (B, C) can have multiple cut-edges
    and/or cut-point as long as there is per face only one (or none) cut across
    a face, i.e. a face can only be split into two faces.

    The information available after construction:
      - pointIsCut, edgeIsCut.
      - faceSplitCut : the cross-cut of a face.
      - cellLoops : the loop which cuts across a cell
      - cellAnchorPoints : per cell the vertices on one side which are
        considered the anchor points.

    AnchorPoints: connected loops have to be oriented in the same way to
    be able to grow new internal faces out of the same bottom faces.
    (limitation of the mapping procedure). The loop is cellLoop is oriented
    such that the normal of it points towards the anchorPoints.
    This is the only routine which uses geometric information.


    Limitations:
    - cut description is very precise. Hard to get right.
    - do anchor points need to be unique per cell? Very inefficient.
    - is orientation guaranteed to be correct? Uses geometric info so can go
      wrong on highly distorted meshes.
    - is memory inefficient. Probably need to use Maps instead of
      labelLists.

SourceFiles
    cellCuts.C

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

#ifndef cellCuts_H
#define cellCuts_H

#include "edgeVertex.H"
#include "labelList.H"
#include "boolList.H"
#include "scalarField.H"
#include "pointField.H"
#include "DynamicList.H"
#include "typeInfo.H"

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

namespace Foam
{

// Forward declaration of classes
class polyMesh;
class cellLooper;
class refineCell;
class plane;

/*---------------------------------------------------------------------------*\
                           Class cellCuts Declaration
\*---------------------------------------------------------------------------*/

class cellCuts
:
    public edgeVertex
{
    // Private data

        // Per point/edge status

            //- Is mesh point cut
            boolList pointIsCut_;

            //- Is edge cut
            boolList edgeIsCut_;

            //- If edge is cut gives weight (0->start() to 1->end())
            scalarField edgeWeight_;


        // Cut addressing

            //- Cuts per existing face (includes those along edge of face)
            //  Cuts in no particular order.
            mutable labelListList* faceCutsPtr_;

            //- Per face : cut across edge (so not along existing edge)
            //  (can only be one per face)
            Map<edge> faceSplitCut_;


        // Cell-cut addressing

            //- Loop across cell circumference
            labelListList cellLoops_;

            //- Number of valid loops in cellLoops_
            label nLoops_;

            //- For each cut cell the points on the 'anchor' side of the cell.
            labelListList cellAnchorPoints_;


    // Private Static Functions

        //- Find value in first n elements of list.
        static label findPartIndex
        (
            const labelList&,
            const label n,
            const label val
        );

        //- Create boolList with all labels specified set to true
        //  (and rest to false)
        static boolList expand(const label size, const labelList& labels);

        //- Create scalarField with all specified labels set to corresponding
        //  value in scalarField.
        static scalarField expand
        (
            const label,
            const labelList&,
            const scalarField&
        );

        //- Returns -1 or index of first element of lst that cannot be found
        //  in map.
        static label firstUnique
        (
            const labelList& lst,
            const Map<label>&
        );

    // Private Member Functions

        //- Debugging: write cell's edges and any cut vertices and edges
        //  (so no cell loop determined yet)
        void writeUncutOBJ(const fileName&, const label cellI) const;

        //- Debugging: write cell's edges, loop and anchors to directory.
        void writeOBJ
        (
            const fileName& dir,
            const label cellI,
            const pointField& loopPoints,
            const labelList& anchors
        ) const;

        //- Find face on cell using the two edges.
        label edgeEdgeToFace
        (
            const label cellI,
            const label edgeA,
            const label edgeB
        ) const;


        //- Find face on cell using an edge and a vertex.
        label edgeVertexToFace
        (
            const label cellI,
            const label edgeI,
            const label vertI
        ) const;

        //- Find face using two vertices (guaranteed not to be along edge)
        label vertexVertexToFace
        (
            const label cellI,
            const label vertA,
            const label vertB
        ) const;


        // Cut addressing

            //- Calculate faceCuts in face vertex order.
            void calcFaceCuts() const;


        // Loop (cuts on cell circumference) calculation

            //- Find edge (or -1) on faceI using vertices v0,v1
            label findEdge
            (
                const label faceI,
                const label v0,
                const label v1
            ) const;

            //- Find face on which all cuts are (very rare) or -1.
            label loopFace(const label cellI, const labelList& loop) const;

            //- Cross otherCut into next faces (not exclude0, exclude1)
            bool walkPoint
            (
                const label cellI,
                const label startCut,

                const label exclude0,
                const label exclude1,

                const label otherCut,

                label& nVisited,
                labelList& visited
            ) const;

            //- Cross cut (which is edge on faceI) onto next face
            bool crossEdge
            (
                const label cellI,
                const label startCut,
                const label faceI,
                const label otherCut,

                label& nVisited,
                labelList& visited
            ) const;

            // wrapper around visited[nVisited++] = cut. Checks for duplicate
            // cuts.
            bool addCut
            (
                const label cellI,
                const label cut,
                label& nVisited,
                labelList& visited
            ) const;

            //- Walk across faceI following cuts, starting at cut. Stores cuts
            //  visited
            // Returns true if valid walk.
            bool walkFace
            (
                const label cellI,
                const label startCut,
                const label faceI,
                const label cut,

                label& lastCut,
                label& beforeLastCut,
                label& nVisited,
                labelList& visited
            ) const;

            //- Walk across cuts (cut edges or cut vertices) of cell. Stops when
            //  hit cut  already visited. Returns true when loop of 3 or more
            //  vertices found.
            bool walkCell
            (
                const label cellI,
                const label startCut,   // overall starting cut
                const label faceI,
                const label prevCut,    // current cut
                label& nVisited,
                labelList& visited
            ) const;

            //- Determine for every cut cell the face it is cut by.
            void calcCellLoops(const labelList& cutCells);


        // Cell anchoring

            //- Are there enough faces on anchor side of cellI?
            bool checkFaces
            (
                const label cellI,
                const labelList& anchorPoints
            ) const;

            //- Walk unset edges of single cell from starting point and
            //  marks visited edges and vertices with status.
            void walkEdges
            (
                const label cellI,
                const label pointI,
                const label status,

                Map<label>& edgeStatus,
                Map<label>& pointStatus
            ) const;

            //- Check anchor points on 'outside' of loop
            bool loopAnchorConsistent
            (
                const label cellI,
                const pointField& loopPts,
                const labelList& anchorPoints
            ) const;

            //- Determines set of anchor points given a loop. The loop should
            //  split the cell into two. Returns true if a valid set of anchor
            //  points determined, false otherwise.
            bool calcAnchors
            (
                const label cellI,
                const labelList& loop,
                const pointField& loopPts,

                labelList& anchorPoints
            ) const;

            //- Returns coordinates of points on loop with explicitly provided
            //  weights.
            pointField loopPoints
            (
                const labelList& loop,
                const scalarField& loopWeights
            ) const;

            //- Returns weights of loop. Inverse of loopPoints.
            scalarField loopWeights(const labelList& loop) const;

            //- Check if cut edges in loop are compatible with ones in
            //  edgeIsCut_
            bool validEdgeLoop
            (
                const labelList& loop,
                const scalarField& loopWeights
            ) const;

            //- Counts number of cuts on face.
            label countFaceCuts
            (
                const label faceI,
                const labelList& loop
            ) const;

            //- Determine compatibility of loop with existing cut pattern.
            //  Does not use cut-addressing (faceCuts_, cutCuts_)
            bool conservativeValidLoop
            (
                const label cellI,
                const labelList& loop
            ) const;

            //- Check if loop is compatible with existing cut pattern in
            //  pointIsCut, edgeIsCut, faceSplitCut.
            //  Calculates and returns for current cell the cut faces and the
            //  points on one side of the loop.
            bool validLoop
            (
                const label cellI,
                const labelList& loop,
                const scalarField& loopWeights,
                Map<edge>& newFaceSplitCut,
                labelList& anchorPoints
            ) const;

            //- Update basic cut information from cellLoops.
            //  Assumes cellLoops_ and edgeWeight_ already set and consistent.
            //  Does not use any other information.
            void setFromCellLoops();

            //- Update basic cut information for single cell from cellLoop.
            bool setFromCellLoop
            (
                const label cellI,
                const labelList& loop,
                const scalarField& loopWeights
            );

            //- Update basic cut information from cellLoops. Checks for
            //  consistency with existing cut pattern.
            void setFromCellLoops
            (
                const labelList& cellLabels,
                const labelListList& cellLoops,
                const List<scalarField>& cellLoopWeights
            );

            //- Cut cells and update basic cut information from cellLoops.
            //  Checks each loop for consistency with existing cut pattern.
            void setFromCellCutter
            (
                const cellLooper&,
                const List<refineCell>& refCells
            );

            //- Same as above but now cut with prescribed plane (instead of
            //  just normal).
            void setFromCellCutter
            (
                const cellLooper&,
                const labelList& cellLabels,
                const List<plane>&
            );

            //- Set orientation of loops
            void orientPlanesAndLoops();

            //- Top level driver: adressing calculation and loop detection
            //  (loops splitting cells).
            void calcLoopsAndAddressing(const labelList& cutCells);

            //- Check various consistencies.
            void check() const;


        //- Disallow default bitwise copy construct
        cellCuts(const cellCuts&);

        //- Disallow default bitwise assignment
        void operator=(const cellCuts&);


public:

    //- Runtime type information
    ClassName("cellCuts");


    // Constructors

        //- Construct from cells to cut and pattern of cuts
        cellCuts
        (
            const polyMesh& mesh,
            const labelList& cutCells,
            const labelList& meshVerts,
            const labelList& meshEdges,
            const scalarField& meshEdgeWeights
        );

        //- Construct from pattern of cuts. Detect cells to cut.
        cellCuts
        (
            const polyMesh& mesh,
            const labelList& meshVerts,
            const labelList& meshEdges,
            const scalarField& meshEdgeWeights
        );

        //- Construct from complete cellLoops through specified cells.
        //  Checks for consistency.
        //  Constructs cut-cut addressing and cellAnchorPoints.
        cellCuts
        (
            const polyMesh& mesh,
            const labelList& cellLabels,
            const labelListList& cellLoops,
            const List<scalarField>& cellEdgeWeights
        );

        //- Construct from list of cells to cut and direction to cut in
        //  (always through cell centre) and celllooper.
        cellCuts
        (
            const polyMesh& mesh,
            const cellLooper& cellCutter,
            const List<refineCell>& refCells
        );

        //- Construct from list of cells to cut and plane to cut with and
        //  celllooper. (constructor above always cuts through cell centre)
        cellCuts
        (
            const polyMesh& mesh,
            const cellLooper& cellCutter,
            const labelList& cellLabels,
            const List<plane>& cutPlanes
        );

        //- Construct from components
        cellCuts
        (
            const polyMesh& mesh,
            const boolList& pointIsCut,
            const boolList& edgeIsCut,
            const scalarField& edgeWeight,
            const Map<edge>& faceSplitCut,
            const labelListList& cellLoops,
            const label nLoops,
            const labelListList& cellAnchorPoints
        );


    //- Destructor
    ~cellCuts();

        //- Clear out demand driven storage
        void clearOut();


    // Member Functions

        // Access

            //- Is mesh point cut
            const boolList& pointIsCut() const
            {
                return pointIsCut_;
            }

            //- Is edge cut
            const boolList& edgeIsCut() const
            {
                return edgeIsCut_;
            }

            //- If edge is cut gives weight (ratio between start() and end())
            const scalarField& edgeWeight() const
            {
                return edgeWeight_;
            }

            //- Cuts per existing face (includes those along edge of face)
            //  Cuts in no particular order
            const labelListList& faceCuts() const
            {
                if (!faceCutsPtr_)
                {
                    calcFaceCuts();
                }
                return *faceCutsPtr_;
            }

            //- Gives for split face the two cuts that split the face into two.
            const Map<edge>& faceSplitCut() const
            {
                return faceSplitCut_;
            }

            //- For each cut cell the cut along the circumference.
            const labelListList& cellLoops() const
            {
                return cellLoops_;
            }

            //- Number of valid cell loops
            label nLoops() const
            {
                return nLoops_;
            }

            //- For each cut cell the points on the 'anchor' side of the cell.
            const labelListList& cellAnchorPoints() const
            {
                return cellAnchorPoints_;
            }


        // Other

            //- Returns coordinates of points on loop for given cell.
            //  Uses cellLoops_ and edgeWeight_
            pointField loopPoints(const label cellI) const;

            //- Invert anchor point selection.
            labelList nonAnchorPoints
            (
                const labelList& cellPoints,
                const labelList& anchorPoints,
                const labelList& loop
            ) const;

            //- Flip loop for cellI. Updates anchor points as well.
            void flip(const label cellI);

            //- Flip loop for cellI. Does not update anchors. Use with care
            //  (only if you're sure loop orientation is wrong)
            void flipLoopOnly(const label cellI);


        // Write

            //- debugging:Write list of cuts to stream in OBJ format
            void writeOBJ
            (
                Ostream& os,
                const pointField& loopPoints,
                label& vertI
            ) const;

            //- debugging:Write all of cuts to stream in OBJ format
            void writeOBJ(Ostream& os) const;

            //- debugging:Write edges of cell and loop
            void writeCellOBJ(const fileName& dir, const label cellI) const;

};


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

} // End namespace Foam

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

#endif

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

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 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 "refinementIterator.H"
#include "polyMesh.H"
#include "Time.H"
#include "refineCell.H"
#include "undoableMeshCutter.H"
#include "polyTopoChange.H"
#include "mapPolyMesh.H"
#include "cellCuts.H"
#include "OFstream.H"
#include "meshTools.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
defineTypeNameAndDebug(refinementIterator, 0);
}


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

// Construct from components
Foam::refinementIterator::refinementIterator
(
    polyMesh& mesh,
    undoableMeshCutter& meshRefiner,
    const cellLooper& cellWalker,
    const bool writeMesh
)
:
    edgeVertex(mesh),
    mesh_(mesh),
    meshRefiner_(meshRefiner),
    cellWalker_(cellWalker),
    writeMesh_(writeMesh)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::refinementIterator::~refinementIterator()
{}


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

Foam::Map<Foam::label> Foam::refinementIterator::setRefinement
(
    const List<refineCell>& refCells
)
{
    Map<label> addedCells(2*refCells.size());

    Time& runTime = const_cast<Time&>(mesh_.time());

    label nRefCells = refCells.size();

    label oldRefCells = -1;

    // Operate on copy.
    List<refineCell> currentRefCells(refCells);

    bool stop = false;

    do
    {
        if (writeMesh_)
        {
            // Need different times to write meshes.
            runTime++;
        }

        polyTopoChange meshMod(mesh_);

        if (debug)
        {
            Pout<< "refinementIterator : refining "
                << currentRefCells.size() << " cells." << endl;
        }

        // Determine cut pattern.
        cellCuts cuts(mesh_, cellWalker_, currentRefCells);

        label nCuts = cuts.nLoops();
        reduce(nCuts, sumOp<label>());

        if (nCuts == 0)
        {
            if (debug)
            {
                Pout<< "refinementIterator : exiting iteration since no valid"
                    << " loops found for " << currentRefCells.size()
                    << " cells" << endl;


                fileName cutsFile("failedCuts_" + runTime.timeName() + ".obj");

                Pout<< "Writing cuts for time " <<  runTime.timeName()
                    << " to " << cutsFile << endl;

                OFstream cutsStream(cutsFile);


                labelList refCellsDebug(currentRefCells.size());
                forAll(currentRefCells, i)
                {
                    refCellsDebug[i] = currentRefCells[i].cellNo();
                }
                meshTools::writeOBJ
                (
                    cutsStream,
                    mesh().cells(),
                    mesh().faces(),
                    mesh().points(),
                    refCellsDebug
                );
            }

            break;
        }

        if (debug)
        {
            fileName cutsFile("cuts_" + runTime.timeName() + ".obj");

            Pout<< "Writing cuts for time " <<  runTime.timeName()
                << " to " << cutsFile << endl;

            OFstream cutsStream(cutsFile);
            cuts.writeOBJ(cutsStream);
        }


        // Insert mesh refinement into polyTopoChange.
        meshRefiner_.setRefinement(cuts, meshMod);


        //
        // Do all changes
        //

        autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh
        (
            mesh_,
            false
        );

        // Move mesh (since morphing does not do this)
        if (morphMap().hasMotionPoints())
        {
            mesh_.movePoints(morphMap().preMotionPoints());
        }

        // Update stored refinement pattern
        meshRefiner_.updateMesh(morphMap());

        // Write resulting mesh
        if (writeMesh_)
        {
            if (debug)
            {
                Pout<< "Writing refined polyMesh to time "
                    << runTime.timeName() << endl;
            }

            mesh_.write();
        }

        // Update currentRefCells for new cell numbers. Use helper function
        // in meshCutter class.
        updateLabels
        (
            morphMap->reverseCellMap(),
            currentRefCells
        );

        // Update addedCells for new cell numbers
        updateLabels
        (
            morphMap->reverseCellMap(),
            addedCells
        );

        // Get all added cells from cellCutter (already in new numbering
        // from meshRefiner.updateMesh call) and add to global list of added
        const Map<label>& addedNow = meshRefiner_.addedCells();

        forAllConstIter(Map<label>, addedNow, iter)
        {
            if (!addedCells.insert(iter.key(), iter()))
            {
                FatalErrorIn("refinementIterator")
                    << "Master cell " << iter.key()
                    << " already has been refined" << endl
                    << "Added cell:" << iter() << abort(FatalError);
            }
        }


        // Get failed refinement in new cell numbering and reconstruct input
        // to the meshRefiner. Is done by removing all refined cells from
        // current list of cells to refine.

        // Update refCells for new cell numbers.
        updateLabels
        (
            morphMap->reverseCellMap(),
            currentRefCells
        );

        // Pack refCells acc. to refined status
        nRefCells = 0;

        forAll(currentRefCells, refI)
        {
            const refineCell& refCell = currentRefCells[refI];

            if (!addedNow.found(refCell.cellNo()))
            {
                if (nRefCells != refI)
                {
                    currentRefCells[nRefCells++] =
                        refineCell
                        (
                            refCell.cellNo(),
                            refCell.direction()
                        );
                }
            }
        }

        oldRefCells = currentRefCells.size();

        currentRefCells.setSize(nRefCells);

        if (debug)
        {
            Pout<< endl;
        }

        // Stop only if all finished or all can't refine any further.
        stop = (nRefCells == 0) || (nRefCells == oldRefCells);
        reduce(stop, andOp<bool>());
    }
    while (!stop);


    if (nRefCells == oldRefCells)
    {
        WarningIn("refinementIterator")
            << "stopped refining."
            << "Did not manage to refine a single cell" << endl
            << "Wanted :" << oldRefCells << endl;
    }

    return addedCells;
}



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

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 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 "refinementIterator.H"
#include "polyMesh.H"
#include "Time.H"
#include "refineCell.H"
#include "undoableMeshCutter.H"
#include "polyTopoChange.H"
#include "mapPolyMesh.H"
#include "cellCuts.H"
#include "OFstream.H"
#include "meshTools.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
defineTypeNameAndDebug(refinementIterator, 0);
}


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

// Construct from components
Foam::refinementIterator::refinementIterator
(
    polyMesh& mesh,
    undoableMeshCutter& meshRefiner,
    const cellLooper& cellWalker,
    const bool writeMesh
)
:
    edgeVertex(mesh),
    mesh_(mesh),
    meshRefiner_(meshRefiner),
    cellWalker_(cellWalker),
    writeMesh_(writeMesh)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::refinementIterator::~refinementIterator()
{}


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

Foam::Map<Foam::label> Foam::refinementIterator::setRefinement
(
    const List<refineCell>& refCells
)
{
    Map<label> addedCells(2*refCells.size());

    Time& runTime = const_cast<Time&>(mesh_.time());

    label nRefCells = refCells.size();

    label oldRefCells = -1;

    // Operate on copy.
    List<refineCell> currentRefCells(refCells);

    bool stop = false;

    do
    {
        if (writeMesh_)
        {
            // Need different times to write meshes.
            runTime++;
        }

        polyTopoChange meshMod(mesh_);

        if (debug)
        {
            Pout<< "refinementIterator : refining "
                << currentRefCells.size() << " cells." << endl;
        }

        // Determine cut pattern.
        cellCuts cuts(mesh_, cellWalker_, currentRefCells);

        label nCuts = cuts.nLoops();
        reduce(nCuts, sumOp<label>());

        if (nCuts == 0)
        {
            if (debug)
            {
                Pout<< "refinementIterator : exiting iteration since no valid"
                    << " loops found for " << currentRefCells.size()
                    << " cells" << endl;


                fileName cutsFile("failedCuts_" + runTime.timeName() + ".obj");

                Pout<< "Writing cuts for time " <<  runTime.timeName()
                    << " to " << cutsFile << endl;

                OFstream cutsStream(cutsFile);


                labelList refCellsDebug(currentRefCells.size());
                forAll(currentRefCells, i)
                {
                    refCellsDebug[i] = currentRefCells[i].cellNo();
                }
                meshTools::writeOBJ
                (
                    cutsStream,
                    mesh().cells(),
                    mesh().faces(),
                    mesh().points(),
                    refCellsDebug
                );
            }

            break;
        }

        if (debug)
        {
            fileName cutsFile("cuts_" + runTime.timeName() + ".obj");

            Pout<< "Writing cuts for time " <<  runTime.timeName()
                << " to " << cutsFile << endl;

            OFstream cutsStream(cutsFile);
            cuts.writeOBJ(cutsStream);
        }


        // Insert mesh refinement into polyTopoChange.
        meshRefiner_.setRefinement(cuts, meshMod);


        //
        // Do all changes
        //

        autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh
        (
            mesh_,
            false
        );

        // Move mesh (since morphing does not do this)
        if (morphMap().hasMotionPoints())
        {
            mesh_.movePoints(morphMap().preMotionPoints());
        }

        // Update stored refinement pattern
        meshRefiner_.updateMesh(morphMap());

        // Write resulting mesh
        if (writeMesh_)
        {
            if (debug)
            {
                Pout<< "Writing refined polyMesh to time "
                    << runTime.timeName() << endl;
            }

            mesh_.write();
        }

        // Update currentRefCells for new cell numbers. Use helper function
        // in meshCutter class.
        updateLabels
        (
            morphMap->reverseCellMap(),
            currentRefCells
        );

        // Update addedCells for new cell numbers
        updateLabels
        (
            morphMap->reverseCellMap(),
            addedCells
        );

        // Get all added cells from cellCutter (already in new numbering
        // from meshRefiner.updateMesh call) and add to global list of added
        const Map<label>& addedNow = meshRefiner_.addedCells();

        forAllConstIter(Map<label>, addedNow, iter)
        {
            if (!addedCells.insert(iter.key(), iter()))
            {
                FatalErrorInFunction
                    << "Master cell " << iter.key()
                    << " already has been refined" << endl
                    << "Added cell:" << iter() << abort(FatalError);
            }
        }


        // Get failed refinement in new cell numbering and reconstruct input
        // to the meshRefiner. Is done by removing all refined cells from
        // current list of cells to refine.

        // Update refCells for new cell numbers.
        updateLabels
        (
            morphMap->reverseCellMap(),
            currentRefCells
        );

        // Pack refCells acc. to refined status
        nRefCells = 0;

        forAll(currentRefCells, refI)
        {
            const refineCell& refCell = currentRefCells[refI];

            if (!addedNow.found(refCell.cellNo()))
            {
                if (nRefCells != refI)
                {
                    currentRefCells[nRefCells++] =
                        refineCell
                        (
                            refCell.cellNo(),
                            refCell.direction()
                        );
                }
            }
        }

        oldRefCells = currentRefCells.size();

        currentRefCells.setSize(nRefCells);

        if (debug)
        {
            Pout<< endl;
        }

        // Stop only if all finished or all can't refine any further.
        stop = (nRefCells == 0) || (nRefCells == oldRefCells);
        reduce(stop, andOp<bool>());
    }
    while (!stop);


    if (nRefCells == oldRefCells)
    {
        WarningInFunction
            << "stopped refining."
            << "Did not manage to refine a single cell" << endl
            << "Wanted :" << oldRefCells << endl;
    }

    return addedCells;
}



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