/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "polyMeshGenChecks.H"
#include "polyMeshGenAddressing.H"
#include "cell.H"
#include "Map.H"

# ifdef USE_OMP
#include <omp.h>
# endif

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

bool Foam::Module::polyMeshGenChecks::checkPoints
(
    const polyMeshGen& mesh,
    const bool report,
    labelHashSet* setPtr
)
{
    label nFaceErrors = 0;
    label nCellErrors = 0;

    const VRWGraph& pf = mesh.addressingData().pointFaces();

    forAll(pf, pointI)
    {
        if (pf.sizeOfRow(pointI) == 0)
        {
            WarningInFunction
                << "Point " << pointI << " not used by any faces." << endl;

            if (setPtr)
            {
                setPtr->insert(pointI);
            }

            ++nFaceErrors;
        }
    }

    const VRWGraph& pc = mesh.addressingData().pointCells();

    forAll(pc, pointI)
    {
        if (pc.sizeOfRow(pointI) == 0)
        {
            WarningInFunction
                << "Point " << pointI << " not used by any cells." << endl;

            if (setPtr)
            {
                setPtr->insert(pointI);
            }

            ++nCellErrors;
        }
    }

    reduce(nFaceErrors, sumOp<label>());
    reduce(nCellErrors, sumOp<label>());

    if (nFaceErrors > 0 || nCellErrors > 0)
    {
        WarningInFunction
            << "Error in point usage detected: " << nFaceErrors
            << " unused points found in the mesh.  This mesh is invalid."
            << endl;

        return true;
    }
    else
    {
        if (report)
        {
            Info<< "Point usage check OK.\n" << endl;
        }

        return false;
    }
}


bool Foam::Module::polyMeshGenChecks::checkUpperTriangular
(
    const polyMeshGen& mesh,
    const bool report,
    labelHashSet* setPtr
)
{
    // Check whether internal faces are ordered in the upper triangular order
    const labelList& own = mesh.owner();
    const labelList& nei = mesh.neighbour();

    const cellListPMG& cells = mesh.cells();
    const VRWGraph& cc = mesh.addressingData().cellCells();

    const label internal = mesh.nInternalFaces();

    labelList checkInternalFaces(internal, -1);

    label nChecks = 0;

    bool error = false;

    // Loop through faceCells once more and make sure that for internal cell
    // the first label is smaller
    for (label faceI = 0; faceI < internal; ++faceI)
    {
        if (own[faceI] >= nei[faceI])
        {
            if (report)
            {
                Pout<< FUNCTION_NAME << endl
                    << "face " << faceI
                    << " has the owner label greater than neighbour:" << endl
                    << own[faceI] << tab << nei[faceI] << endl;
            }

            if (setPtr)
            {
                setPtr->insert(faceI);
            }

            error = true;
        }
    }

    // Loop through all cells. For each cell, find the face that is internal and
    // add it to the check list (upper triangular order).
    // Once the list is completed, check it against the faceCell list
    forAll(cells, cellI)
    {
        const labelList& curFaces = cells[cellI];

        // Using the fact that cell neighbour always appear
        // in the increasing order
        boolList usedNbr(cc.sizeOfRow(cellI), false);

        for (label nSweeps = 0; nSweeps < usedNbr.size(); ++nSweeps)
        {
            // Find the lowest neighbour which is still valid
            label nextNei = -1;
            label minNei = cells.size();

            forAllRow(cc, cellI, nbrI)
            {
                const label neiI = cc(cellI, nbrI);
                if ((neiI > cellI) && !usedNbr[nbrI] && (neiI < minNei))
                {
                    nextNei = nbrI;
                    minNei = neiI;
                }
            }

            if (nextNei > -1)
            {
                // Mark this neighbour as used
                usedNbr[nextNei] = true;

                forAll(curFaces, faceI)
                {
                    if (curFaces[faceI] < internal)
                    {
                        if (nei[curFaces[faceI]] == cc(cellI, nextNei))
                        {
                            checkInternalFaces[nChecks] = curFaces[faceI];
                            ++nChecks;

                            break;
                        }
                    }
                }
            }
        }
    }

    // Check list created. If everything is OK, the face label is equal to index
    forAll(checkInternalFaces, faceI)
    {
        if (checkInternalFaces[faceI] != faceI)
        {
            error = true;

            Pout<< FUNCTION_NAME << endl
                << "face " << faceI << " out of position. Markup label: "
                << checkInternalFaces[faceI] << ". All subsequent faces will "
                << "also be out of position. Please check the mesh manually."
                << endl;

            if (setPtr)
            {
                setPtr->insert(faceI);
            }

            break;
        }
    }

    reduce(error, orOp<bool>());

    if (error)
    {
        WarningInFunction
            << "Error in face ordering: faces not in upper triangular order!"
            << endl;

        return true;
    }
    else
    {
        if (report)
        {
            Info<< "Upper triangular ordering OK.\n" << endl;
        }

        return false;
    }
}


bool Foam::Module::polyMeshGenChecks::checkCellsZipUp
(
    const polyMeshGen& mesh,
    const bool report,
    labelHashSet* setPtr
)
{
    label nOpenCells = 0;

    const faceListPMG& faces = mesh.faces();
    const cellListPMG& cells = mesh.cells();

    # ifdef USE_OMP
    # pragma omp parallel for schedule(guided) reduction(+ : nOpenCells)
    # endif
    forAll(cells, cellI)
    {
        const labelList& c = cells[cellI];

        DynList<edge> cellEdges;
        DynList<label> edgeUsage;

        forAll(c, faceI)
        {
            const face& f = faces[c[faceI]];

            forAll(f, eI)
            {
                const edge e = f.faceEdge(eI);

                const label pos = cellEdges.find(e);

                if (pos < 0)
                {
                    cellEdges.append(e);
                    edgeUsage.append(1);
                }
                else
                {
                    ++edgeUsage[pos];
                }
            }
        }

        DynList<edge> singleEdges;

        forAll(edgeUsage, edgeI)
        {
            if (edgeUsage[edgeI] == 1)
            {
                singleEdges.append(cellEdges[edgeI]);
            }
            else if (edgeUsage[edgeI] != 2)
            {
                WarningInFunction
                    << "edge " << cellEdges[edgeI] << " in cell " << cellI
                    << " used " << edgeUsage[edgeI] << " times. " << endl
                    << "Should be 1 or 2 - serious error in mesh structure"
                    << endl;

                if (setPtr)
                {
                    # ifdef USE_OMP
                    # pragma omp critical
                    # endif
                    setPtr->insert(cellI);
                }
            }
        }

        if (singleEdges.size() > 0)
        {
            if (report)
            {
                Pout<< FUNCTION_NAME << endl
                    << "Cell " << cellI << " has got " << singleEdges.size()
                    << " unmatched edges: " << singleEdges << endl;
            }

            if (setPtr)
            {
                # ifdef USE_OMP
                # pragma omp critical
                # endif
                setPtr->insert(cellI);
            }

            ++nOpenCells;
        }
    }

    reduce(nOpenCells, sumOp<label>());

    if (nOpenCells > 0)
    {
        WarningInFunction
            << nOpenCells
            << " open cells found.  Please use the mesh zip - up tool. "
            << endl;

        return true;
    }
    else
    {
        if (report)
        {
            Info<< "Topological cell zip - up check OK.\n" << endl;
        }

        return false;
    }
}


bool Foam::Module::polyMeshGenChecks::checkFaceVertices
(
    const polyMeshGen& mesh,
    const bool report,
    labelHashSet* setPtr
)
{
    // Check that all vertex labels are valid
    const faceListPMG& faces = mesh.faces();

    label nErrorFaces = 0;
    const label nPoints = mesh.points().size();

    forAll(faces, fI)
    {
        const face& curFace = faces[fI];

        if (min(curFace) < 0 || max(curFace) > nPoints)
        {
            WarningInFunction
                << "Face " << fI << " contains vertex labels out of range: "
                << curFace << " Max point index = " << nPoints - 1 << endl;

            if (setPtr)
            {
                setPtr->insert(fI);
            }

            ++nErrorFaces;
        }

        // Uniqueness of vertices
        labelHashSet facePoints(2*curFace.size());

        forAll(curFace, fp)
        {
            bool inserted = facePoints.insert(curFace[fp]);

            if (!inserted)
            {
                WarningInFunction
                    << "Face " << fI << " contains duplicate vertex labels: "
                    << curFace << endl;

                if (setPtr)
                {
                    setPtr->insert(fI);
                }

                ++nErrorFaces;
            }
        }
    }

    reduce(nErrorFaces, sumOp<label>());

    if (nErrorFaces > 0)
    {
        SeriousErrorInFunction
            << nErrorFaces << " faces with invalid vertex labels found"
            << endl;

        return true;
    }
    else
    {
        if (report)
        {
            Info<< "Face vertices OK.\n" << endl;
        }

        return false;
    }
}


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