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

    geophase: A C++ template library for geometrical operations on non-convex
    polyhedra with non-planar faces
     
    Copyright (C) 2019 Tomislav Maric, TU Darmstadt 

    This program 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.

    This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.

Description
    Triangulation of flux polyhedra. Take into account the topological structure
    of the flux polyhedron. Make it possible to triangulate non-convex objects 
    that arise from face-coplanar point velocities.

Author
    Tomislav Maric maric@mma.tu-darmstadt.de

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

#ifndef FluxTriangulate_H 
#define FluxTriangulate_H 

#include "Triangulate.H" 
#include "DistanceSqr.H"
#include <algorithm>

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

namespace Foam { namespace GeometricalTransport { 

    // Assumes the following structure for fluxPolyhedron 
    // [basePolygon, sidePolygon0, sidePolygon1, ... sidePolygonS, capPolygon]
    // 0,1,..S indices are oriented like the basePolygon. TM.
    template
    <
        typename Triangulation = tetrahedronVector, 
        typename Polyhedron
    >
    typename std::enable_if
    <
        Traits::tag_enabled<Triangulation, tetrahedron_sequence_tag>::value &&
        Traits::tag_enabled<Polyhedron, polygon_collection_tag>::value,
        Triangulation
    >::type 
    flux_triangulate(Polyhedron const & fluxPolyhedron)
    {
        Triangulation result;

        typedef typename Triangulation::value_type tetrahedron; 
        typedef typename Polyhedron::value_type::value_type point;  

        // Triangulate the sides of the fluxPolyhedron and compute the volume centroid.  
        triangleVector sidesTriangulation;
        auto sideIter = std::next(fluxPolyhedron.begin()); 
        auto sideEnd = std::prev(fluxPolyhedron.end()); 
        for (; sideIter != sideEnd; ++sideIter)
        {
            auto sideTriangles = edge_triangulate<triangleVector>(*sideIter); 

            for (const auto& triangle : sideTriangles)
                sidesTriangulation.push_back(triangle); 
        }

        // Triangulate the base and cap polygon of the flux polyhedron. 
        auto capTriangulation = 
            barycentric_triangulate<triangleVector>(fluxPolyhedron.back()); 

        auto capProjection(fluxPolyhedron.back()); 
        const auto& basePolygon = fluxPolyhedron.front(); 
        // Compute the normal of the base polygon.
        point baseNormal = normal_area_vec(fluxPolyhedron.front());  
        baseNormal /= mag(baseNormal); 
        // Project the cap polygon onto to the base polygon plane.
        for (auto & projectedPoint : capProjection)
            projectedPoint -= 
                (projectedPoint - basePolygon.front()) & baseNormal * baseNormal; 

        std::vector<point> intersectionPoints; 
        auto baseIt = basePolygon.begin(); 
        auto capIt = capProjection.begin(); 
        
        for(; baseIt != basePolygon.end(); ++baseIt, ++capIt)
        {
            if (is_inside(capProjection, *baseIt, baseNormal))
                intersectionPoints.push_back(*baseIt); 

            if (is_inside(basePolygon, *capIt, baseNormal))
                intersectionPoints.push_back(*capIt); 
        }
       
        point volumeCentroid = std::accumulate(
            intersectionPoints.begin(), 
            intersectionPoints.end(), 
            point(0,0,0)
        );

        if (intersectionPoints.size() > 0)
            volumeCentroid /= intersectionPoints.size(); 
        else 
            volumeCentroid = centroid(capProjection); 

        // Triangulate the volume from the sides triangulation.
        for (const auto& triangle : sidesTriangulation)
        {
            result.push_back(
                build<tetrahedron>(
                    volumeCentroid, 
                    triangle[0], 
                    triangle[1],
                    triangle[2]
                )
            ); 
        }

        // Trianngulate volume from the cap and base triangulation.
        for (const auto& triangle : capTriangulation)
        {
            result.push_back(
                build<tetrahedron>(
                    volumeCentroid, 
                    triangle[0], 
                    triangle[1], 
                    triangle[2]
                )
            ); 
        }

        return result; 
    }

    template
    <
        typename Triangulation = tetrahedronVector, 
        typename MixedPolyhedron
    >
    typename std::enable_if
    <
        Traits::tag_enabled<Triangulation, tetrahedron_sequence_tag>::value &&
        Traits::tag_enabled<MixedPolyhedron, mixed_polyhedron_tag>::value,
        Triangulation
    >::type 
    flux_triangulate(MixedPolyhedron const & fluxMixedPolyhedron)
    {
        Triangulation result;
        
        typedef typename Triangulation::value_type tetrahedron; 
        typedef typename MixedPolyhedron::PointType point;  

        // Compute the normal of the base polygon.
        point baseNormal = normal_area_vec(fluxMixedPolyhedron.polygons_front());
        baseNormal /= mag(baseNormal); 

        // Project the cap points to the base polygon. 
        const auto& basePolygon = fluxMixedPolyhedron.polygons_front(); 
        typename MixedPolyhedron::PolygonType capProjection;
        // Built the cap polygon from the cap triangulation. 
        for (auto& triangle : fluxMixedPolyhedron.triangles())
        {
            capProjection.push_back(triangle[1]);
        }
        // Project the cap polygon to the base polygon plane. 
        for (auto & projectedPoint : capProjection)
            projectedPoint -= 
                (projectedPoint - basePolygon.front()) & baseNormal * baseNormal; 

        // Pick up the intersection points: inside both polygons.
        std::vector<point> intersectionPoints; 
        auto baseIt = basePolygon.begin(); 
        auto capIt = capProjection.begin(); 
        
        for(; baseIt != basePolygon.end(); ++baseIt, ++capIt)
        {
            if (is_inside(capProjection, *baseIt, baseNormal))
                intersectionPoints.push_back(*baseIt); 

            if (is_inside(basePolygon, *capIt, baseNormal))
                intersectionPoints.push_back(*capIt); 
        }
       
        // Set volume centroid as the average intersection points. 
        point volumeCentroid = std::accumulate(
                intersectionPoints.begin(), 
                intersectionPoints.end(), 
                point(0,0,0)
            );

        if (intersectionPoints.size() > 0)
            volumeCentroid /= intersectionPoints.size(); 
        else 
            volumeCentroid = centroid(capProjection); 

        //for (const auto& polygon : fluxMixedPolyhedron.polygons())
        for (auto it = std::next(fluxMixedPolyhedron.polygons_begin()); 
             it != fluxMixedPolyhedron.polygons_end(); 
             ++it)
        {
            auto triangles = edge_triangulate<triangleVector>(*it); 

            for (const auto& triangle : triangles)
            {
                result.push_back(
                    build<tetrahedron>(
                        volumeCentroid, 
                        triangle[0], 
                        triangle[1],
                        triangle[2]
                    )
                ); 
            }
        }

        for (const auto& triangle : fluxMixedPolyhedron.triangles())
        {
            result.push_back(
                build<tetrahedron>(
                    volumeCentroid, 
                    triangle[0], 
                    triangle[1],
                    triangle[2]
                )
            ); 
        }
        return result; 
    }

    template <typename Triangulation> 
    struct flux_triangulation 
    {
        template<typename Geometry>
        static Triangulation apply(Geometry const & geo)
        {
            return flux_triangulate<Triangulation>(geo);
        }

        template<typename Geometry>
        Triangulation operator()(Geometry const & geo) const
        {
            return apply(geo); 
        }
    };

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

}} // End namespace Foam::GeometricalTransport

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

#endif

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