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

Application
    surfaceHookUp

Description
    Find close open edges and stitches the surface along them

Usage
    - surfaceHookUp hookDistance [OPTION]

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

#include "argList.H"
#include "Time.H"

#include "triSurfaceMesh.H"
#include "indexedOctree.H"
#include "treeBoundBox.H"
#include "PackedBoolList.H"
#include "unitConversion.H"
#include "searchableSurfaces.H"

using namespace Foam;

// Split faceI along edgeI at position newPointI
void greenRefine
(
    const triSurface& surf,
    const label faceI,
    const label edgeI,
    const label newPointI,
    DynamicList<labelledTri>& newFaces
)
{
    const labelledTri& f = surf.localFaces()[faceI];
    const edge& e = surf.edges()[edgeI];

    // Find index of edge in face.

    label fp0 = findIndex(f, e[0]);
    label fp1 = f.fcIndex(fp0);
    label fp2 = f.fcIndex(fp1);

    if (f[fp1] == e[1])
    {
        // Edge oriented like face
        newFaces.append
        (
            labelledTri
            (
                f[fp0],
                newPointI,
                f[fp2],
                f.region()
            )
        );
        newFaces.append
        (
            labelledTri
            (
                newPointI,
                f[fp1],
                f[fp2],
                f.region()
            )
        );
    }
    else
    {
        newFaces.append
        (
            labelledTri
            (
                f[fp2],
                newPointI,
                f[fp1],
                f.region()
            )
        );
        newFaces.append
        (
            labelledTri
            (
                newPointI,
                f[fp0],
                f[fp1],
                f.region()
            )
        );
    }
}


//scalar checkEdgeAngle
//(
//    const triSurface& surf,
//    const label edgeIndex,
//    const label pointIndex,
//    const scalar& angle
//)
//{
//    const edge& e = surf.edges()[edgeIndex];

//    vector eVec = e.vec(surf.localPoints());
//    eVec /= mag(eVec) + SMALL;

//    const labelList& pEdges = surf.pointEdges()[pointIndex];
//
//    forAll(pEdges, eI)
//    {
//        const edge& nearE = surf.edges()[pEdges[eI]];

//        vector nearEVec = nearE.vec(surf.localPoints());
//        nearEVec /= mag(nearEVec) + SMALL;

//        const scalar dot = eVec & nearEVec;
//        const scalar minCos = degToRad(angle);

//        if (mag(dot) > minCos)
//        {
//            return false;
//        }
//    }

//    return true;
//}


void createBoundaryEdgeTrees
(
    const PtrList<triSurfaceMesh>& surfs,
    PtrList<indexedOctree<treeDataEdge> >& bEdgeTrees,
    labelListList& treeBoundaryEdges
)
{
    forAll(surfs, surfI)
    {
        const triSurface& surf = surfs[surfI];

        // Boundary edges
        treeBoundaryEdges[surfI] =
            labelList
            (
                identity(surf.nEdges() - surf.nInternalEdges())
              + surf.nInternalEdges()
            );

        Random rndGen(17301893);

        // Slightly extended bb. Slightly off-centred just so on symmetric
        // geometry there are less face/edge aligned items.
        treeBoundBox bb
        (
            treeBoundBox(UList<point>(surf.localPoints())).extend(rndGen, 1e-4)
        );

        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);

        bEdgeTrees.set
        (
            surfI,
            new indexedOctree<treeDataEdge>
            (
                treeDataEdge
                (
                    false,                      // cachebb
                    surf.edges(),               // edges
                    surf.localPoints(),         // points
                    treeBoundaryEdges[surfI]    // selected edges
                ),
                bb,     // bb
                8,      // maxLevel
                10,     // leafsize
                3.0     // duplicity
            )
       );
    }
}


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

int main(int argc, char *argv[])
{
    argList::addNote
    (
        "hook surfaces to other surfaces by moving and retriangulating their"
        "boundary edges to match other surface boundary edges"
    );
    argList::noParallel();
    argList::validArgs.append("hookTolerance");

#   include "addDictOption.H"

#   include "setRootCase.H"
#   include "createTime.H"

    const word dictName("surfaceHookUpDict");
#   include "setSystemRunTimeDictionaryIO.H"

    Info<< "Reading " << dictName << nl << endl;

    const IOdictionary dict(dictIO);

    const scalar dist(args.argRead<scalar>(1));
    const scalar matchTolerance(SMALL);

    Info<< "Hooking distance = " << dist << endl;

    searchableSurfaces surfs
    (
        IOobject
        (
            "surfacesToHook",
            runTime.constant(),
            "triSurface",
            runTime
        ),
        dict
    );

    Info<< nl << "Reading surfaces: " << endl;
    forAll(surfs, surfI)
    {
        Info<< incrIndent;
        Info<< nl << indent << "Surface     = " << surfs.names()[surfI] << endl;

        const wordList& regions = surfs[surfI].regions();
        forAll(regions, surfRegionI)
        {
            Info<< incrIndent;
            Info<< indent << "Regions = " << regions[surfRegionI] << endl;
            Info<< decrIndent;
        }
        Info<< decrIndent;
    }

    PtrList<indexedOctree<treeDataEdge> > bEdgeTrees(surfs.size());
    labelListList treeBoundaryEdges(surfs.size());

    List<DynamicList<labelledTri> > newFaces(surfs.size());
    List<DynamicList<point> > newPoints(surfs.size());
    List<PackedBoolList> visitedFace(surfs.size());

    PtrList<triSurfaceMesh> newSurfaces(surfs.size());
    forAll(surfs, surfI)
    {
        const triSurfaceMesh& surf =
            refCast<const triSurfaceMesh>(surfs[surfI]);

        newSurfaces.set
        (
            surfI,
            new triSurfaceMesh
            (
                IOobject
                (
                    "hookedSurface_" + surfs.names()[surfI],
                    runTime.constant(),
                    "triSurface",
                    runTime
                ),
                surf
            )
        );
    }

    label nChanged = 0;
    label nIters = 0;

    do
    {
        Info<< nl << "Iteration = " << nIters++ << endl;
        nChanged = 0;

        createBoundaryEdgeTrees(newSurfaces, bEdgeTrees, treeBoundaryEdges);

        forAll(newSurfaces, surfI)
        {
            const triSurface& newSurf = newSurfaces[surfI];

            newFaces[surfI] = newSurf.localFaces();
            newPoints[surfI] = newSurf.localPoints();
            visitedFace[surfI] = PackedBoolList(newSurf.size(), false);
        }

        forAll(newSurfaces, surfI)
        {
            const triSurface& surf = newSurfaces[surfI];

            List<pointIndexHit> bPointsTobEdges(surf.boundaryPoints().size());
            labelList bPointsHitTree(surf.boundaryPoints().size(), -1);

            const labelListList& pointEdges = surf.pointEdges();

            forAll(bPointsTobEdges, bPointI)
            {
                pointIndexHit& nearestHit = bPointsTobEdges[bPointI];

                const label pointI = surf.boundaryPoints()[bPointI];
                const point& samplePt = surf.localPoints()[pointI];

                const labelList& pEdges = pointEdges[pointI];

                // Add edges connected to the edge to the shapeMask
                DynamicList<label> shapeMask;
                shapeMask.append(pEdges);

                forAll(bEdgeTrees, treeI)
                {
                    const indexedOctree<treeDataEdge>& bEdgeTree =
                        bEdgeTrees[treeI];

                    pointIndexHit currentHit =
                        bEdgeTree.findNearest
                        (
                            samplePt,
                            sqr(dist),
                            treeDataEdge::findNearestOpSubset
                            (
                                bEdgeTree,
                                shapeMask
                            )
                        );

                    if
                    (
                        currentHit.hit()
                     &&
                        (
                            !nearestHit.hit()
                         ||
                            (
                                magSqr(currentHit.hitPoint() - samplePt)
                              < magSqr(nearestHit.hitPoint() - samplePt)
                            )
                        )
                    )
                    {
                        nearestHit = currentHit;
                        bPointsHitTree[bPointI] = treeI;
                    }
                }

                scalar dist2 = magSqr(nearestHit.rawPoint() - samplePt);

                if (nearestHit.hit())
                {
    //                bool rejectEdge =
    //                    checkEdgeAngle
    //                    (
    //                        surf,
    //                        nearestHit.index(),
    //                        pointI,
    //                        30
    //                    );

                    if (dist2 > Foam::sqr(dist))
                    {
                        nearestHit.setMiss();
                    }
                }
            }

            forAll(bPointsTobEdges, bPointI)
            {
                const pointIndexHit& eHit = bPointsTobEdges[bPointI];

                if (eHit.hit())
                {
                    const label hitSurfI = bPointsHitTree[bPointI];
                    const triSurface& hitSurf = newSurfaces[hitSurfI];

                    const label eIndex =
                        treeBoundaryEdges[hitSurfI][eHit.index()];
                    const edge& e = hitSurf.edges()[eIndex];

                    const label pointI = surf.boundaryPoints()[bPointI];

                    const labelList& eFaces = hitSurf.edgeFaces()[eIndex];

                    if (eFaces.size() != 1)
                    {
                        WarningIn(args.executable())
                            << "Edge is attached to " << eFaces.size()
                            << " faces." << endl;

                        continue;
                    }

                    const label faceI = eFaces[0];

                    if (visitedFace[hitSurfI][faceI])
                    {
                        continue;
                    }

                    DynamicList<labelledTri> newFacesFromSplit(2);

                    const point& pt = surf.localPoints()[pointI];

                    if
                    (
                        (
                            magSqr(pt - hitSurf.localPoints()[e.start()])
                          < matchTolerance
                        )
                     || (
                            magSqr(pt - hitSurf.localPoints()[e.end()])
                          < matchTolerance
                        )
                    )
                    {
                        continue;
                    }

                    nChanged++;

                    // Keep the points in the same place and move the edge
                    newPoints[hitSurfI].append(newPoints[surfI][pointI]);

                    // Move the points to the edges
                    //newPoints[pointI] = eHit.hitPoint();
                    //newPoints.append(eHit.hitPoint());

                    visitedFace[hitSurfI][faceI] = true;

                    // Split the other face.
                    greenRefine
                    (
                        hitSurf,
                        faceI,
                        eIndex,
                        newPoints[hitSurfI].size() - 1,
                        newFacesFromSplit
                    );

                    forAll(newFacesFromSplit, newFaceI)
                    {
                        if (newFaceI == 0)
                        {
                            newFaces[hitSurfI][faceI] = newFacesFromSplit[0];
                        }
                        else
                        {
                            newFaces[hitSurfI].append
                            (
                                newFacesFromSplit[newFaceI]
                            );
                        }
                    }
                }
            }
        }

        Info<< "    Number of edges split = " << nChanged << endl;

        forAll(newSurfaces, surfI)
        {
            newSurfaces.set
            (
                surfI,
                new triSurfaceMesh
                (
                    IOobject
                    (
                        "hookedSurface_" + surfs.names()[surfI],
                        runTime.constant(),
                        "triSurface",
                        runTime
                    ),
                    triSurface
                    (
                        newFaces[surfI],
                        newSurfaces[surfI].patches(),
                        pointField(newPoints[surfI])
                    )
                )
            );
        }

    } while (nChanged > 0);

    Info<< endl;

    forAll(newSurfaces, surfI)
    {
        const triSurfaceMesh& newSurf = newSurfaces[surfI];

        Info<< "Writing hooked surface " << newSurf.searchableSurface::name()
            << endl;

        newSurf.searchableSurface::write();
    }

    Info<< "\nEnd\n" << endl;

    return 0;
}


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