/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "triSurfaceImportSurfaceAsSubset.H"
#include "meshOctree.H"
#include "meshOctreeCreator.H"
#include "helperFunctions.H"
#include "demandDrivenData.H"

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

void Foam::Module::triSurfaceImportSurfaceAsSubset::createOctree
(
    const triSurf& /*surf*/,
    meshOctree& /*octree*/
)
{}


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

Foam::Module::triSurfaceImportSurfaceAsSubset::triSurfaceImportSurfaceAsSubset
(
    triSurf& surface
)
:
    surf_(surface),
    octreePtr_(nullptr)
{}


Foam::Module::triSurfaceImportSurfaceAsSubset::
~triSurfaceImportSurfaceAsSubset()
{
    deleteDemandDrivenData(octreePtr_);
}


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

void Foam::Module::triSurfaceImportSurfaceAsSubset::addSurfaceAsSubset
(
    const triSurf& importSurf,
    const word& subsetName,
    const scalar angleTol
)
{
    if (!octreePtr_)
    {
        octreePtr_ = new meshOctree(surf_);
        meshOctreeCreator(*octreePtr_).createOctreeWithRefinedBoundary
        (
            direction(20),
            15
        );
    }

    const pointField& points = surf_.points();
    const vectorField& fNornals = surf_.facetNormals();
    const vectorField& fCentres = surf_.facetCentres();

    labelList nearestTriangle(importSurf.size(), -1);

    // check which triangles in the surface fit best to the centres of the
    // triangles in the import surface
    const pointField& importSurfPoints = importSurf.points();
    const vectorField& importFaceCentres = importSurf.facetCentres();
    const vectorField& importFaceNormals = importSurf.facetNormals();
    # ifdef USE_OMP
    # pragma omp parallel for schedule(dynamic, 40)
    # endif
    forAll(nearestTriangle, triI)
    {
        point np;
        scalar dSq;
        label nt, patch;

        octreePtr_->findNearestSurfacePoint
        (
            np,
            dSq,
            nt,
            patch,
            importFaceCentres[triI]
        );

        // find the longest edge distance
        scalar maxEdgeDSq(0.);
        const labelledTri& tri = importSurf[triI];
        forAll(tri, pI)
        {
            const point& s = importSurfPoints[tri[pI]];
            const point& e = importSurfPoints[tri[(pI + 1)%3]];

            maxEdgeDSq = max(maxEdgeDSq, magSqr(e - s));
        }

        // check if the triangle has been found
        if ((nt < 0) || (dSq > 0.09*maxEdgeDSq))
        {
            Warning << "Could not find a matching triangle " << endl;
            Warning << "It seems that your surface meshes do not overlap" << endl;
            continue;
        }

        vector nTri = importFaceNormals[triI];
        const scalar magSqrTri = magSqr(nTri);

        // skip sliver triangles
        if (magSqrTri < VSMALL)
            continue;

        vector normal = fNornals[nt];
        const scalar dSqNormal = magSqr(normal);

        // skip sliver triangles
        if (dSqNormal < VSMALL)
            continue;

        if (((nTri & normal)/(magSqrTri*dSqNormal)) > angleTol)
            nearestTriangle[triI] = nt;
    }

    meshOctree otherSurfOctree(importSurf);
    meshOctreeCreator(otherSurfOctree).createOctreeWithRefinedBoundary(20, 15);

    // search for nearest facets in the import surface
    DynList<label> containedTriangles;
    # ifdef USE_OMP
    # pragma omp parallel for schedule(dynamic, 40) private(containedTriangles)
    # endif
    forAll(surf_, triI)
    {
        // find the bounding box and the ize of the triangle
        boundBox bb(fCentres[triI], fCentres[triI]);

        scalar maxEdgeDSq(0.);
        const labelledTri& tri = surf_[triI];
        forAll(tri, pI)
        {
            // bounding box of the surface triangle
            bb.min() = min(bb.min(), points[tri[pI]]);
            bb.max() = max(bb.max(), points[tri[pI]]);

            const point& s = points[tri[pI]];
            const point& e = points[tri[(pI + 1)%3]];

            maxEdgeDSq = max(maxEdgeDSq, magSqr(e-s));
        }

        // find the nearest triangle in the surface which shall be imported
        otherSurfOctree.findTrianglesInBox(bb, containedTriangles);

        label nt(-1);
        scalar dSq(VGREAT);
        forAll(containedTriangles, ctI)
        {
            const point p =
                help::nearestPointOnTheTriangle
                (
                    containedTriangles[ctI],
                    importSurf,
                    fCentres[triI]
                );

            const scalar distSq = magSqr(p - fCentres[triI]);

            if (distSq < dSq)
            {
                nt = containedTriangles[ctI];
                dSq = distSq;
            }
        }

        // check if the triangle has been found
        if ((nt < 0) || (dSq > 0.09*maxEdgeDSq))
            continue;

        // skip firther checkes f it has found the same triangle
        if (nearestTriangle[nt] == triI)
            continue;

        vector nTri = fNornals[triI];
        const scalar magSqrTri = magSqr(nTri);

        // skip sliver triangles
        if (magSqrTri < VSMALL)
            continue;

        vector normal = importFaceNormals[nt];
        const scalar dSqNormal = magSqr(normal);

        // skip sliver triangles
        if (dSqNormal < VSMALL)
            continue;

        if (((nTri & normal)/(magSqrTri*dSqNormal)) > angleTol)
            nearestTriangle[nt] = triI;
    }

    // create a facet subset in the surface mesh and add the facets into it
    const label subsetId = surf_.addFacetSubset(subsetName);

    forAll(nearestTriangle, triI)
    {
        if (nearestTriangle[triI] < 0)
            continue;

        surf_.addFacetToSubset(subsetId, nearestTriangle[triI]);
    }
}


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