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

    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 EdgeTriangulate_H 
#define EdgeTriangulate_H 

#include "Triangulate.H" 

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

namespace Foam { namespace GeometricalTransport { 

    template
    <
        typename Triangulation = triangleVector, 
        typename PointCollection
    >
    typename std::enable_if
    <
        Traits::tag_enabled<Triangulation, triangle_collection_tag>::value &&
        Traits::tag_enabled<PointCollection, point_collection_tag>::value,
        Triangulation
    >::type 
    edge_triangulate(PointCollection const& polygon)
    {
        Triangulation result; 

        const auto& P0 = polygon.front();  
        const auto& P1 = *std::next(polygon.begin()); 
        const auto& P2 = *std::next(polygon.begin(),2); 
        const auto& P3 = polygon.back();  

        const auto D02 = 0.5 * (P0 + P2); 
        const auto D13 = 0.5 * (P1 + P3); 

        const auto P0D02 = D02 - P0; 
        const auto P0D13 = D13 - P0; 

        auto E01 = P1 - P0; 
        E01 /= mag(E01); 

        // Testing : Arithmetic average.
        //auto sideCentroid = 0.5 * (D02 + D13);  

        // Minimal distance. 
        //const auto hD02 = P0D02 - (P0D02 & E01) * E01; 
        //const auto hD13 = P0D13 - (P0D13 & E01) * E01; 
        //auto sideCentroid(D13); 
        //if ((hD02 & hD02) < (hD13 & hD13))
           //sideCentroid = D02;  
           
        // Inversed distance weighted. 
        double wD02 = 1.0 / (mag(P0D02 - (P0D02 & E01) * E01) + 1e-15); 
        double wD13 = 1.0 / (mag(P0D13 - (P0D13 & E01) * E01) + 1e-15); 
        auto sideCentroid = ((D02 * wD02) + (D13 * wD13)) / (wD02 + wD13); 

        result = point_triangulate<Triangulation>(polygon, sideCentroid); 

        return result; 
    }

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

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

}} // End namespace Foam::GeometricalTransport

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

#endif

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