/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | cfMesh: A library for mesh generation
   \\    /   O peration     |
    \\  /    A nd           | www.cfmesh.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2014-2017 Creative Fields, Ltd.
-------------------------------------------------------------------------------
Author
     Franjo Juretic (franjo.juretic@c-fields.com)

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 "demandDrivenData.H"
#include "meshUntangler.H"

//#define DEBUGSmooth

# ifdef DEBUGSmooth
#include "helperFunctions.H"
# endif

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

void Foam::Module::meshUntangler::cutRegion::findNewFaces()
{
    #ifdef DEBUGSmooth
    Info<< "Finding new faces " << endl;
    #endif

    bool finished;
    do
    {
        finished = true;

        const DynList<DynList<label, 8>, 64>& fcs = *facesPtr_;
        DynList<edge, 128>& cEdges = *cEdgesPtr_;

        cFacesPtr_ = new DynList<DynList<label, 8>, 64>();
        DynList<DynList<label, 8>, 64>& cFaces = *cFacesPtr_;

        DynList<label, 8> faceInPlane;

        DynList<label, 64> pointUsage;
        pointUsage.setSize(cPtsPtr_->size());

        forAll(fcs, fI)
        {
            const DynList<label, 8>& f = fcs[fI];

            # ifdef DEBUGSmooth
            Info<< "Creating new face from face " << fI
                << " consisting of edges " << f << endl;
            # endif

            pointUsage = 0;

            DynList<label, 8> newFace;

            forAll(f, eI)
            {
                # ifdef DEBUGSmooth
                const DynList<edge>& edges = *edgesPtr_;
                Info<< "Vertex types for face edge " << eI << " are "
                    << label(vertexTypes_[edges[f[eI]].start()]) << " and "
                    << label(vertexTypes_[edges[f[eI]].end()]) << endl;
                # endif

                const label edgeLabel = newEdgeLabel_[f[eI]];

                if (edgeLabel != -1)
                {
                    # ifdef DEBUGSmooth
                    Info<< "Orig edge " << eI << " " << edges[f[eI]]
                        << " is replaced with " << cEdges[edgeLabel] << endl;
                    # endif

                    const edge& e = cEdges[edgeLabel];
                    ++pointUsage[e[0]];
                    ++pointUsage[e[1]];
                    newFace.append(edgeLabel);
                }
            }

            if (newFace.size() > 1)
            {
                DynList<label, 4> newEdge;
                forAll(pointUsage, pI)
                    if (pointUsage[pI] == 1)
                        newEdge.append(pI);

                if (newEdge.size() == 2)
                {
                    # ifdef DEBUGSmooth
                    Info<< "Storing new edge " << newEdge << endl;
                    # endif

                    newFace.append(cEdges.size());
                    cEdges.append(edge(newEdge[0], newEdge[1]));
                }
                else if (newEdge.size() > 2)
                {
                    # ifdef DEBUGSmooth
                    Info<< "New edge " << newEdge << endl;
                    # endif

                    tieBreak(f);
                    if (!valid_) return;
                    finished = false;
                    break;

                    FatalErrorInFunction
                        << "Edge has more than two nodes!"
                        << abort(FatalError);
                }

                cFaces.append(newFace);
            }
        }

        if (!finished) continue;

        // find edges which form the faceInPlane
        DynList<label, 128> edgeUsage;
        edgeUsage.setSize(cEdges.size());
        edgeUsage = 0;
        forAll(cFaces, fI)
        {
            const DynList<label, 8>& f = cFaces[fI];

            forAll(f, eI)
                ++edgeUsage[f[eI]];
        }

        forAll(edgeUsage, eI)
            if (edgeUsage[eI] == 1)
                faceInPlane.append(eI);

        if (faceInPlane.size() > 2)
        {
            # ifdef DEBUGSmooth
            Info<< "Adding face in plane " << faceInPlane << endl;
            Info<< "Face in plane consists of edges " << endl;
            forAll(faceInPlane, eI)
                Info<< "Edge " << eI << " is "
                    << cEdges[faceInPlane[eI]] << endl;
            # endif

            cFaces.append(faceInPlane);
        }

        # ifdef DEBUGSmooth
        Info<< "cEdges " << cEdges << endl;
        Info<< "Number of faces before cutting " << fcs.size() << endl;
        Info<< "Found " << cFaces.size() << " new faces" << endl;
        forAll(fcs, fI)
        {
            Info<< "Old face " << fI << " contains edges " << fcs[fI] << endl;
        }

        forAll(cFaces, fI)
        {
            Info<< "New face " << fI << " contains edges " << cFaces[fI] << endl;
        }

        // test if the region is closed
        List<DynList<label, 4>> eFaces(cEdges.size());
        forAll(cFaces, fI)
        {
            const DynList<label, 8>& f = cFaces[fI];
            forAll(f, eI)
                eFaces[f[eI]].append(fI);
        }

        if (eFaces.size() > 5)
        forAll(eFaces, fI)
            if (eFaces[fI].size() != 2)
            {
                Info<< "eFaces " << eFaces << endl;
                Info<< "cEdges " << cEdges << endl;
                Info<< "cFaces " << cFaces << endl;

                FatalErrorInFunction
                    << "Cell is not topologically closed!" << abort(FatalError);
            }
        # endif

    } while (!finished);
}


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