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

// Cells to refine; name of cell set
set c0;

// Type of coordinate system:
// - global : coordinate system same for every cell. Usually aligned with
//   x,y,z axis. Specify in globalCoeffs section below.
// - patchLocal : coordinate system different for every cell. Specify in
//   patchLocalCoeffs section below.
// - fieldBased : uses the list of field names from the directions list for
//   selecting the directions to cut. Meant to be used with geometricCut, but
//   can also be used with useHexTopology.
coordinateSystem global;
//coordinateSystem patchLocal;
//coordinateSystem fieldBased;

// .. and its coefficients. x,y in this case. (normal direction is calculated
// as tan1^tan2)
globalCoeffs
{
    tan1 (1 0 0);
    tan2 (0 1 0);
}

patchLocalCoeffs
{
    patch outside;  // Normal direction is facenormal of zero'th face of patch
    tan1 (1 0 0);
}

// List of directions to refine, if global or patchLocal
directions
(
    tan1
    tan2
    normal
);

// List of directions to refine, if "fieldBased". Keep in mind that these
// fields must be of type "vectorField", not "volVectorField".
//directions
//(
//    radialDirectionFieldName
//    angularDirectionFieldName
//    heightDirectionFieldName
//);

// Whether to use hex topology. This will
// - if patchLocal: all cells on selected patch should be hex
// - split all hexes in 2x2x2 through the middle of edges.
useHexTopology  true;

// Cut purely geometric (will cut hexes through vertices) or take topology
// into account. Incompatible with useHexTopology
geometricCut    false;

// Write meshes from intermediate steps
writeMesh       false;

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

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "directions.H"
#include "polyMesh.H"
#include "twoDPointCorrector.H"
#include "directionInfo.H"
#include "MeshWave.H"
#include "OFstream.H"
#include "meshTools.H"
#include "hexMatcher.H"
#include "Switch.H"
#include "globalMeshData.H"

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

namespace Foam
{
    template<>
    const char* Foam::NamedEnum
    <
        Foam::directions::directionType,
        3
    >::names[] =
    {
        "tan1",
        "tan2",
        "normal"
    };
}

const Foam::NamedEnum<Foam::directions::directionType, 3>
    Foam::directions::directionTypeNames_;


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

// For debugging
void Foam::directions::writeOBJ(Ostream& os, const point& pt)
{
    os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
}

// For debugging
void Foam::directions::writeOBJ
(
    Ostream& os,
    const point& pt0,
    const point& pt1,
    label& vertI
)
{
    writeOBJ(os, pt0);
    writeOBJ(os, pt1);

    os << "l " << vertI + 1 << ' ' << vertI + 2 << endl;

    vertI += 2;
}


// Dump to file.
void Foam::directions::writeOBJ
(
    const fileName& fName,
    const primitiveMesh& mesh,
    const vectorField& dirs
)
{
    Pout<< "Writing cell info to " << fName << " as vectors at the cellCentres"
        << endl << endl;

    OFstream xDirStream(fName);

    label vertI = 0;

    forAll(dirs, cellI)
    {
        const point& ctr = mesh.cellCentres()[cellI];

        // Calculate local length scale
        scalar minDist = GREAT;

        const labelList& nbrs = mesh.cellCells()[cellI];

        forAll(nbrs, nbrI)
        {
            minDist = min(minDist, mag(mesh.cellCentres()[nbrs[nbrI]] - ctr));
        }

        scalar scale = 0.5*minDist;

        writeOBJ(xDirStream, ctr, ctr + scale*dirs[cellI], vertI);
    }
}


void Foam::directions::check2D
(
    const twoDPointCorrector* correct2DPtr,
    const vector& vec
)
{
    if (correct2DPtr)
    {
        if (mag(correct2DPtr->planeNormal() & vec) > 1e-6)
        {
            FatalErrorInFunction
                << "is not normal to plane defined in dynamicMeshDict."
                << endl
                << "Either make case 3D or adjust vector."
                << exit(FatalError);
        }
    }
}


// Get direction on all cells
Foam::vectorField Foam::directions::propagateDirection
(
    const polyMesh& mesh,
    const bool useTopo,
    const polyPatch& pp,
    const vectorField& ppField,
    const vector& defaultDir
)
{
    // Seed all faces on patch
    labelList changedFaces(pp.size());
    List<directionInfo> changedFacesInfo(pp.size());

    if (useTopo)
    {
        forAll(pp, patchFaceI)
        {
            label meshFaceI = pp.start() + patchFaceI;

            label cellI = mesh.faceOwner()[meshFaceI];

            if (!hexMatcher().isA(mesh, cellI))
            {
                FatalErrorInFunction
                    << "useHexTopology specified but cell " << cellI
                    << " on face " << patchFaceI << " of patch " << pp.name()
                    << " is not a hex" << exit(FatalError);
            }

            const vector& cutDir = ppField[patchFaceI];

            // Get edge(bundle) on cell most in direction of cutdir
            label edgeI = meshTools::cutDirToEdge(mesh, cellI, cutDir);

            // Convert edge into index on face
            label faceIndex =
                directionInfo::edgeToFaceIndex
                (
                    mesh,
                    cellI,
                    meshFaceI,
                    edgeI
                );

            // Set initial face and direction
            changedFaces[patchFaceI] = meshFaceI;
            changedFacesInfo[patchFaceI] =
                directionInfo
                (
                    faceIndex,
                    cutDir
                );
        }
    }
    else
    {
        forAll(pp, patchFaceI)
        {
            changedFaces[patchFaceI] = pp.start() + patchFaceI;
            changedFacesInfo[patchFaceI] =
                directionInfo
                (
                    -2,         // Geometric information only
                    ppField[patchFaceI]
                );
        }
    }

    MeshWave<directionInfo> directionCalc
    (
        mesh,
        changedFaces,
        changedFacesInfo,
        mesh.globalData().nTotalCells()+1
    );

    const List<directionInfo>& cellInfo = directionCalc.allCellInfo();

    vectorField dirField(cellInfo.size());

    label nUnset = 0;
    label nGeom = 0;
    label nTopo = 0;

    forAll(cellInfo, cellI)
    {
        label index = cellInfo[cellI].index();

        if (index == -3)
        {
            // Never visited
            WarningInFunction
                << "Cell " << cellI << " never visited to determine "
                << "local coordinate system" << endl
                << "Using direction " << defaultDir << " instead" << endl;

            dirField[cellI] = defaultDir;

            nUnset++;
        }
        else if (index == -2)
        {
            // Geometric direction
            dirField[cellI] = cellInfo[cellI].n();

            nGeom++;
        }
        else if (index == -1)
        {
            FatalErrorInFunction
                << "Illegal index " << index << endl
                << "Value is only allowed on faces" << abort(FatalError);
        }
        else
        {
            // Topological edge cut. Convert into average cut direction.
            dirField[cellI] = meshTools::edgeToCutDir(mesh, cellI, index);

            nTopo++;
        }
    }

    Pout<< "Calculated local coords for " << defaultDir
        << endl
        << "    Geometric cut cells   : " << nGeom << endl
        << "    Topological cut cells : " << nTopo << endl
        << "    Unset cells           : " << nUnset << endl
        << endl;

    return dirField;
}


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

Foam::directions::directions
(
    const polyMesh& mesh,
    const dictionary& dict,
    const twoDPointCorrector* correct2DPtr
)
:
    List<vectorField>(wordList(dict.lookup("directions")).size())
{
    const wordList wantedDirs(dict.lookup("directions"));
    const word coordSystem(dict.lookup("coordinateSystem"));

    bool wantNormal = false;
    bool wantTan1 = false;
    bool wantTan2 = false;
    label nDirs = 0;

    if (coordSystem != "fieldBased")
    {
        forAll(wantedDirs, i)
        {
            directionType wantedDir = directionTypeNames_[wantedDirs[i]];

            if (wantedDir == NORMAL)
            {
                wantNormal = true;
            }
            else if (wantedDir == TAN1)
            {
                wantTan1 = true;
            }
            else if (wantedDir == TAN2)
            {
                wantTan2 = true;
            }
        }
    }


    if (coordSystem == "global")
    {
        const dictionary& globalDict = dict.subDict("globalCoeffs");

        vector tan1(globalDict.lookup("tan1"));
        check2D(correct2DPtr, tan1);

        vector tan2(globalDict.lookup("tan2"));
        check2D(correct2DPtr, tan2);

        vector normal = tan1 ^ tan2;
        normal /= mag(normal);

        Pout<< "Global Coordinate system:" << endl
            << "     normal : " << normal << endl
            << "     tan1   : " << tan1 << endl
            << "     tan2   : " << tan2
            << endl << endl;

        if (wantNormal)
        {
            operator[](nDirs++) = vectorField(1, normal);
        }
        if (wantTan1)
        {
            operator[](nDirs++) = vectorField(1, tan1);
        }
        if (wantTan2)
        {
            operator[](nDirs++) = vectorField(1, tan2);
        }
    }
    else if (coordSystem == "patchLocal")
    {
        const dictionary& patchDict = dict.subDict("patchLocalCoeffs");

        const word patchName(patchDict.lookup("patch"));

        const label patchI = mesh.boundaryMesh().findPatchID(patchName);

        if (patchI == -1)
        {
            FatalErrorInFunction
                << "Cannot find patch "
                << patchName
                << exit(FatalError);
        }

        // Take zeroth face on patch
        const polyPatch& pp = mesh.boundaryMesh()[patchI];

        vector tan1(patchDict.lookup("tan1"));

        const vector& n0 = pp.faceNormals()[0];

        if (correct2DPtr)
        {
            tan1 = correct2DPtr->planeNormal() ^ n0;

            WarningInFunction
                << "Discarding user specified tan1 since 2D case." << endl
                << "Recalculated tan1 from face normal and planeNormal as "
                << tan1 << endl << endl;
        }

        Switch useTopo(dict.lookup("useHexTopology"));

        vectorField normalDirs;
        vectorField tan1Dirs;

        if (wantNormal || wantTan2)
        {
            normalDirs =
                propagateDirection
                (
                    mesh,
                    useTopo,
                    pp,
                    pp.faceNormals(),
                    n0
                );

            if (wantNormal)
            {
                this->operator[](nDirs++) = normalDirs;
            }
        }

        if (wantTan1 || wantTan2)
        {
            tan1Dirs =
                propagateDirection
                (
                    mesh,
                    useTopo,
                    pp,
                    vectorField(pp.size(), tan1),
                    tan1
                );


            if (wantTan1)
            {
                this->operator[](nDirs++) = tan1Dirs;
            }
        }
        if (wantTan2)
        {
            tmp<vectorField> tan2Dirs = normalDirs ^ tan1Dirs;

            this->operator[](nDirs++) = tan2Dirs;
        }
    }
    else if (coordSystem == "fieldBased")
    {
        forAll(wantedDirs, i)
        {
            operator[](nDirs++) =
                vectorIOField
                (
                    IOobject
                    (
                        mesh.instance()/wantedDirs[i],
                        mesh,
                        IOobject::MUST_READ,
                        IOobject::NO_WRITE
                    )
                );
        }
    }
    else
    {
        FatalErrorInFunction
            << "Unknown coordinate system "
            << coordSystem << endl
            << "Known types are global, patchLocal and fieldBased"
            << exit(FatalError);
    }
}


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

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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::directions

Description
    Set of directions for each cell in the mesh. Either uniform and size=1
    or one set of directions per cell.

    Used in splitting cells.
    Either all cells have similar refinement direction ('global') or
    direction is dependent on local cell geometry, or loads selected fields
    by name ('fieldBased'). Controlled by dictionary.

SourceFiles
    directions.C

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

#ifndef directions_H
#define directions_H

#include "List.H"
#include "vectorField.H"
#include "NamedEnum.H"
#include "point.H"

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

namespace Foam
{

// Forward declaration of classes
class polyMesh;
class twoDPointCorrector;
class primitiveMesh;
class polyPatch;
class dictionary;

/*---------------------------------------------------------------------------*\
                           Class directions Declaration
\*---------------------------------------------------------------------------*/

class directions
:
    public List<vectorField>
{
public:
    // Data types

        //- Enumeration listing the possible coordinate directions.
        enum directionType
        {
            TAN1,
            TAN2,
            NORMAL
        };

private:

        static const NamedEnum<directionType, 3> directionTypeNames_;


    // Private Member Functions


        //- For debugging. Write point coordinate.
        static void writeOBJ(Ostream& os, const point& pt);

        //- For debugging. Write edge between two points.
        static void writeOBJ
        (
            Ostream& os,
            const point& pt0,
            const point& pt1,
            label& vertI
        );

        //- For debugging. Write hedgehog display of vectorField as obj file.
        static void writeOBJ
        (
            const fileName& fName,
            const primitiveMesh& mesh,
            const vectorField& dirs
        );

        //- Check if vec has no component in 2D normal direction. Exits if
        //  so.
        static void check2D
        (
            const twoDPointCorrector* correct2DPtr,
            const vector& vec
        );

        //- Get coordinate direction for all cells in mesh by propagating from
        //  vector on patch.
        static vectorField propagateDirection
        (
            const polyMesh& mesh,
            const bool useTopo,
            const polyPatch& pp,
            const vectorField& ppField,
            const vector& defaultDir
        );

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

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


public:

    // Constructors

        //- Construct from mesh and dictionary and optional 2D corrector.
        directions
        (
            const polyMesh& mesh,
            const dictionary& dict,
            const twoDPointCorrector* correct2DPtr = NULL
        );

};


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

} // 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." << nl
                    << "    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);
}


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