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

    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/>.

Application
    geophaseTestTriangulation 

Description
    Triangulation tests.  
    
Authors
    Tomislav Maric maric@mma.tu-darmstadt.de

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

// Tested
#include "Triangulate.hpp"

// Testing utilities
#include "Make.hpp"
#include "gtest/gtest.h"
#include "testUtilities.hpp"
#include "WriteVtkPolyData.hpp"

using namespace geophase; 


// Helper functions for testing the triangulation. 
template<typename PolygonSequence>
vectorPolygon pointtri_polygon(PolygonSequence& triangulation)
{
    vectorPolygon result; 

    reserve(result, triangulation.size());  

    for (const auto& triangle : triangulation)
    {
        const auto P2it = std::next(triangle.begin()); 
        result.push_back(*P2it);
    }

    return result;
}

// Assumes the oriented triangulation was done starting with the first
// point of the polygon. Generally this must not be the case. TM.
template<typename PolygonSequence>
vectorPolygon orientedtri_polygon(PolygonSequence const& triangulation)
{
    vectorPolygon result; 

    // A triangulation only has a single triangle.
    if (triangulation.size() == 1)
    {
        const auto& tri = triangulation.front(); 
        auto pIt = tri.begin();
        result.push_back(*pIt);
        ++pIt; 
        result.push_back(*pIt);
        ++pIt; 
        result.push_back(*pIt);
        return result;
    }

    result.reserve(triangulation.size() + 2);  

    // The first two points of the first triangle 
    const auto& T0 = triangulation.front(); 
    result.push_back(*T0.begin()); 
    result.push_back(*std::next(T0.begin())); 

    auto triIt = std::next(triangulation.begin());
    const auto lastIt = std::prev(triangulation.end());

    for (;triIt != lastIt ; ++triIt)
    {
        auto triP1 = std::next(triIt->begin());
        result.push_back(*triP1);
    }

    const auto& TN = triangulation.back();
    auto TN1it = std::next(TN.begin());
    result.push_back(*TN1it); 
    result.push_back(TN.back());

    return result;
}

template<typename PolygonSequence>
vectorPolygon centroidtri_polygon(PolygonSequence const& triangulation)
{
    vectorPolygon result; 

    result.reserve(triangulation.size());  

    for (const auto& triangle : triangulation)
    {
        const auto TP1 = std::next(triangle.begin());
        result.push_back(*TP1);
    }

    return result;
}

template<typename Triangulation> 
bool triangulation_has_only_triangles(Triangulation const& triangulation)
{
    for (const auto& triangle : triangulation)
        if (triangle.size() != 3)
            return false; 
    return true;
}

template<typename Polyhedron, typename Polygon>
void test_polygon_triangulations(Polygon const& testPoly)
{
    auto centroidTri = centroid_triangulate<Polyhedron>(testPoly); 
    ASSERT_TRUE(triangulation_has_only_triangles(centroidTri));
    const auto centroidTriPoly = centroidtri_polygon(centroidTri); 
    ASSERT_TRUE(std::equal(centroidTriPoly.begin(), 
                           centroidTriPoly.end(), testPoly.begin()))
        << "Triangulation = " << centroidTri << std::endl
        << "Triangulation polygon = " << centroidTriPoly << std::endl
        << "Original polygon = " << testPoly << std::endl;
    
    auto orientedTri = oriented_triangulate<Polyhedron>(testPoly);   
    ASSERT_TRUE(triangulation_has_only_triangles(orientedTri));
    const auto orientedTriPoly = orientedtri_polygon(orientedTri); 
    ASSERT_TRUE(std::equal(orientedTriPoly.begin(), 
                           orientedTriPoly.end(), testPoly.begin()))
        << "Triangulation = " << orientedTri << std::endl
        << "Triangulation polygon = " << orientedTriPoly << std::endl
        << "Original polygon = " << testPoly << std::endl;

    ASSERT_TRUE(centroidTri.size() == testPoly.size());
    ASSERT_TRUE(orientedTri.size() == testPoly.size() - 2);
}

template<typename Polygon> 
void test_sequence_polygon_triangulations(Polygon const& polygon)
{
    test_polygon_triangulations<vectorPolyhedron>(polygon);
    test_polygon_triangulations<listPolyhedron>(polygon);
    test_polygon_triangulations<dequePolyhedron>(polygon);
}

TEST(centroid, triangle)
{
    triangle arrayTriangle {point{0.,0.,0.}, point{1.,0.,0.}, point{0.,1.,0.}}; 
    vectorPolygon polyTriangle {point{0.,0.,0.}, point{1.,0.,0.}, point{0.,1.,0.}}; 

    auto arrCentroid = centroid(arrayTriangle); 
    auto polyCentroid = centroid(polyTriangle); 

    ASSERT_DOUBLE_EQ(mag(arrCentroid - polyCentroid), 0.); 

    ASSERT_DOUBLE_EQ(arrCentroid[0], 1.0 / 3.0);
    ASSERT_DOUBLE_EQ(arrCentroid[1], 1.0 / 3.0);
    ASSERT_DOUBLE_EQ(polyCentroid[0], 1.0 / 3.0);
    ASSERT_DOUBLE_EQ(polyCentroid[1], 1.0 / 3.0);
}

TEST(centroid, square)
{
    auto arraySquare = make_unit_square<rectangle>();
    auto polySquare =  make_unit_square<vectorPolygon>();

    auto arrCentroid = centroid(arraySquare); 
    auto polyCentroid = centroid(polySquare); 

    ASSERT_DOUBLE_EQ(mag(arrCentroid - polyCentroid), 0.); 
    ASSERT_DOUBLE_EQ(arrCentroid[0], 0.5);
    ASSERT_DOUBLE_EQ(arrCentroid[1], 0.5);
    ASSERT_DOUBLE_EQ(polyCentroid[0], 0.5);
    ASSERT_DOUBLE_EQ(polyCentroid[1], 0.5);

}

TEST(centroid, cube)
{
    auto cube = make_unit_cube(); 

    auto cubeCentroid = centroid(cube);

    ASSERT_EQ(mag(cubeCentroid - point(0.5,0.5,0.5)), 0);
}

TEST(triangulate, triangle)
{
    const auto triangle = make_unit_triangle();  
    test_sequence_polygon_triangulations(triangle);
}

TEST(triangulate, square)
{
    const auto square = make_unit_square();
    test_sequence_polygon_triangulations(square);
}
    
TEST(triangulate, hexagon)
{
    const auto hexa = make_unit<hexagon>();
    test_sequence_polygon_triangulations(hexa);
}

template<typename Triangulation, typename Polyhedron>
void test_polyhedron_triangulations(Polyhedron const& testPolyhedron)
{
    auto centroidTri = centroid_triangulate<Polyhedron>(testPolyhedron); 
    ASSERT_TRUE(triangulation_has_only_triangles(centroidTri));
    write_vtk(centroidTri, fullTestName() + ".centroidTriangulation.vtk");

    decltype(testPolyhedron.size()) nCentroidTriangles {0};  
    if (testPolyhedron.size() > 4)
    {
        for (const auto& polygon : testPolyhedron)
            nCentroidTriangles += polygon.size(); 
        ASSERT_EQ(centroidTri.size(), nCentroidTriangles);
    }
    else
        ASSERT_EQ(centroidTri.size(), 4);

    auto orientedTri = oriented_triangulate<Polyhedron>(testPolyhedron); 
    ASSERT_TRUE(triangulation_has_only_triangles(orientedTri));
    write_vtk(orientedTri, fullTestName() + ".orientTriangulation.vtk");

    decltype(testPolyhedron.size()) nOrientedTriangles {0};  
    if (testPolyhedron.size() > 4)
    {
        for (const auto& polygon : testPolyhedron)
            nOrientedTriangles += polygon.size() - 2; 
        ASSERT_EQ(orientedTri.size(), nOrientedTriangles);
    }
    else
        ASSERT_EQ(orientedTri.size(), 4);
}

template<typename Polyhedron>
void test_sequence_polyhedron_triangulations(Polyhedron const& polyhedron)
{
    test_polyhedron_triangulations<vectorPolyhedron>(polyhedron);
    test_polyhedron_triangulations<listPolyhedron>(polyhedron);
    test_polyhedron_triangulations<dequePolyhedron>(polyhedron);
}

TEST(triangulate, tetrahedron)
{
    const auto tetrahedron = make_unit_tetrahedron();
    test_sequence_polyhedron_triangulations(tetrahedron);
}

TEST(triangulate, cube)
{
    const auto cube = make_unit_cube();
    test_sequence_polyhedron_triangulations(cube);
}

TEST(triangulate, dodecahedron)
{
    const auto dodecahedron = make_unit_dodecahedron();
    test_sequence_polyhedron_triangulations(dodecahedron);
}

TEST(triangulate, nonPlanarDodecahedron)
{
    const auto dodecahedron = make_unit_dodecahedron(0.2);
    test_sequence_polyhedron_triangulations(dodecahedron);
}

TEST(triangulate, endoDodecahedron)
{
    const auto endoDodecahedron = make_unit_endo_dodecahedron();
    test_sequence_polyhedron_triangulations(endoDodecahedron);
}

TEST(triangulate, nonPlanarEndoDodecahedron)
{
    const auto endoDodecahedron = make_unit_endo_dodecahedron(0.2);
    test_sequence_polyhedron_triangulations(endoDodecahedron);
}

TEST(volTriangulate, tetrahedron)
{
    const auto tetrahedron = make_unit_tetrahedron();

    const auto triangulation = centroid_vol_triangulate(tetrahedron); 

    static_assert(std::is_same_v<tag_t<decltype(triangulation)>,void>);

    write_vtk(triangulation, fullTestName() + ".centroidTriangulation.vtk");

    ASSERT_TRUE(triangulation.size() == 1);

    for (const auto& tetra : triangulation)
    {
        ASSERT_TRUE(tetra.size() == 4); 
        for (const auto& tri : tetra)
            ASSERT_TRUE(tri.size() == 3); 
    }
}

TEST(volTriangulate, cube)
{
    const auto cube = make_unit_cube();

    auto triangulation = centroid_vol_triangulate(cube); 

    write_vtk(triangulation, fullTestName() + ".centroidTriangulation.vtk");

    ASSERT_TRUE(triangulation.size() == 24);

    for (const auto& tetra : triangulation)
    {
        ASSERT_TRUE(tetra.size() == 4); 
        for (const auto& tri : tetra)
            ASSERT_TRUE(tri.size() == 3); 
    }
}

TEST(volTriangulate, dodecahedron)
{
    const auto dodecahedron = make_unit_dodecahedron();

    auto triangulation = centroid_vol_triangulate(dodecahedron); 

    write_vtk(triangulation, fullTestName() + ".centroidTriangulation.vtk");

    ASSERT_TRUE(triangulation.size() == 12 * 5);

    for (const auto& tetra : triangulation)
    {
        ASSERT_TRUE(tetra.size() == 4); 
        for (const auto& tri : tetra)
            ASSERT_TRUE(tri.size() == 3); 
    }
}

TEST(volTriangulate, endoDodecahedron)
{
    const auto endoDodecahedron = make_unit_endo_dodecahedron();

    auto triangulation = centroid_vol_triangulate(endoDodecahedron); 

    write_vtk(triangulation, fullTestName() + ".centroidTriangulation.vtk");

    ASSERT_TRUE(triangulation.size() == 12 * 5);

    for (const auto& tetra : triangulation)
    {
        ASSERT_TRUE(tetra.size() == 4); 
        for (const auto& tri : tetra)
            ASSERT_TRUE(tri.size() == 3); 
    }
}

//TEST(triangulation, mixedPolyhedron)
//{
    //typedef MixedPolyhedron<arrayTriangle, pointVector, std::vector> mixedPolyhedron;

    //auto cube = build<pointVectorVector>(point(0,0,0), point(1,1,1));

    //auto mixedPoly = build<mixedPolyhedron>(cube, point(0.5,0.5,1.5), 0);

    //write_vtk_polydata(mixedPoly, fullTestName() + ".mixedPolyhedron.vtk"); 

    //auto mixedPolyTri =  barycentric_triangulate<tetrahedronVector>(mixedPoly); 
    //write_vtk_polydata(mixedPolyTri, fullTestName() + ".mixedPolyTri.vtk");  

    //ASSERT_TRUE(mixedPolyTri.size() == 24);

    //auto mixedPolyFluxTri = flux_triangulate<tetrahedronVector>(mixedPoly); 
    //write_vtk_polydata(mixedPolyFluxTri, fullTestName() + ".mixedPolyFluxTri.vtk");  

    //ASSERT_LE(mag(volume(mixedPolyFluxTri) - volume(mixedPolyTri)), EPSILON); 
//}

//TEST(triangulation, displacementSweep)
//{
    //pointVector basePolygon = {point(0,0,0), point(1,0,0), point(1,1,0), point(0,1,0)};

    //vector displacement (1,1,1); 
    
    //auto sweptPolyhedron = build<pointVectorVector>(basePolygon, displacement);  

    //write_vtk_polydata(sweptPolyhedron, fullTestName() + ".sweptPlyhedron.vtk"); 

    //ASSERT_LE(mag(volume(sweptPolyhedron) - 1), EPSILON); 
//}

//TEST(prismTriangulation, uniformDisplacementListSweep)
//{
    //pointVector basePolygon = {point(0,0,0), point(0,1,0), point(1,1,0), point(1,0,0)};
    
    //pointVector displacements = {point(1,1,1), point(1,1,1), point(1,1,1), point(1,1,1)};
     
    //auto sweptPolyhedron = build<pointVectorVector>(basePolygon, displacements);  
    //auto barycentricTriangulation = barycentric_triangulate<tetrahedronVector>(sweptPolyhedron);
    //auto fluxTriangulation = flux_triangulate(sweptPolyhedron);

    //write_vtk_polydata(sweptPolyhedron, fullTestName() + ".sweptPolyhedron.vtk"); 
    //write_vtk_polydata(fluxTriangulation, fullTestName() + ".fluxTriangulation.vtk");

    //auto baryTriVolume = volume(barycentricTriangulation); 
    //auto fluxTriVolume = volume(fluxTriangulation); 
    //auto fluxTriVolumeError = mag(fluxTriVolume - 1);  

    //ASSERT_LE(mag(baryTriVolume - fluxTriVolume), EPSILON)
        //<< "barycentric triangulation volume = " << baryTriVolume << "\n"
        //<< "flux triangulation volume = " << fluxTriVolume << "\n";

    //ASSERT_LE(fluxTriVolumeError, EPSILON)
        //<< "flux triangulation volume error = " 
        //<< fluxTriVolumeError << "\n";
//}

//TEST(prismTriangulation, uniformDisplacementListHalfSweep)
//{
    //pointVector basePolygon = {point(0,0,0), point(0,1,0), point(1,1,0), point(1,0,0)};
    
    //pointVector displacements = {point(0.5,0.5,0.5), point(0.5,0.5,0.5), point(0.5,0.5,0.5), point(0.5,0.5,0.5)};
     
    //auto sweptPolyhedron = build<pointVectorVector>(basePolygon, displacements);  
    //auto barycentricTriangulation = barycentric_triangulate<tetrahedronVector>(sweptPolyhedron);
    //auto fluxTriangulation = flux_triangulate(sweptPolyhedron);

    //write_vtk_polydata(sweptPolyhedron, fullTestName() + ".sweptPolyhedron.vtk"); 
    //write_vtk_polydata(fluxTriangulation, fullTestName() + ".fluxTriangulation.vtk");

    //auto baryTriVolume = volume(barycentricTriangulation); 
    //auto fluxTriVolume = volume(fluxTriangulation); 
    //auto fluxTriVolumeError = mag(fluxTriVolume - 0.5);  

    //ASSERT_LE(mag(baryTriVolume - fluxTriVolume), EPSILON)
        //<< "barycentric triangulation volume = " << baryTriVolume << "\n"
        //<< "flux triangulation volume = " << fluxTriVolume << "\n";

    //ASSERT_LE(fluxTriVolumeError, EPSILON)
        //<< "flux triangulation volume error = " 
        //<< fluxTriVolumeError << "\n";
//}

//TEST(prismTriangulation, coplanarEdgePairTrianglePrism)
//{
    //double Ymax = 1; 
    //double Xmax = 1; 

    //pointVector basePolygon = {point(0,0,0), point(0,Ymax,0), point(Xmax,Ymax,0), point(Xmax,0,0)};

    //const int N = 100; 

    //OFstream outputFile (fullTestName() + ".dat"); 

    //for (unsigned int I = 0; I < N+1; ++I)
    //{
        //double s = I * (1.0 / N);
        //double correctVolume = 0.5 * (1 - s) * s * Ymax; 
        
        //pointVector displacements = {point(s,0,0), point(s,0,0), point(0,0,s), point(0,0,s)};
         
        //auto sweptPolyhedron = build<pointVectorVector>(basePolygon, displacements);  
        //auto fluxTriangulation = flux_triangulate(sweptPolyhedron);

        //// Uncomment to visualize. Generates a lot of files. 
        //write_vtk_polydata(fluxTriangulation, fullTestName() + ".fluxTriangulation", I);

        //auto fluxTriVolume = volume(fluxTriangulation); 
        //auto fluxTriVolumeError = mag(fluxTriVolume - correctVolume);  

        //ASSERT_LE(fluxTriVolumeError, EPSILON)
            //<< "normalized flux triangulation volume error = " 
            //<< fluxTriVolumeError / correctVolume << "\n"
            //<< "s " << s << "\n";

        //outputFile << s << " " << fluxTriVolumeError / correctVolume << "\n";
    //}
//}

// A non convex prism where one edge pair is coplanar with the basePolygon, and the other one
// is rotating with increasing spherical angle theta. 
//TEST(prismTriangulation, coplanarEdgePairTrianglePrismParameterised)
//{
    //double Ymax = 1; 
    //double Xmax = 1; 
    //double edgeLength = 0.5;  // Prevent self-intersection with 0.5.

    //pointVector basePolygon = {point(0,0,0), point(0,Ymax,0), point(Xmax,Ymax,0), point(Xmax,0,0)};

    //const int N = 50; 

    //OFstream outputFile (fullTestName() + ".dat"); 

    //unsigned int counter = 0; 
    //for (unsigned int I = 0; I < N+1; ++I)
    //{
        //double s = I * (edgeLength / N);

        //for (unsigned int J = 0; J < N+1; ++J)
        //{
            //double theta = (-0.5 * M_PI) + (J * (M_PI / N)); 

            //auto displacementX = s * Foam::sin(theta); 
            //auto displacementZ = s * Foam::cos(theta); 
            
            //pointVector displacements = {
                //point(0.5,0,0), 
                //point(0.5,0,0), 
                //point(displacementX, 0, displacementZ), 
                //point(displacementX, 0, displacementZ)
            //};

            //// Compute the sides of the triangle generated by the sweep. 
            //auto a = displacements[3];
            //auto b = (basePolygon[0] + displacements[0]) - basePolygon[3]; 

            //double correctVolume = 0.5 * mag(a ^ b) * Ymax;  
             
            //auto sweptPolyhedron = build<pointVectorVector>(basePolygon, displacements);  
            //auto fluxTriangulation = flux_triangulate(sweptPolyhedron);


            //auto fluxTriVolume = volume(fluxTriangulation); 
            //auto fluxTriVolumeError = mag(fluxTriVolume - correctVolume);  

            //outputFile << s << " " << fluxTriVolumeError / correctVolume << "\n";
            //// Uncomment control to visualize only in case of errors.
            //if (fluxTriVolumeError > EPSILON)
            //{
                //write_vtk_polydata(sweptPolyhedron, fullTestName() + ".sweptPolyhedron", counter); 
                //write_vtk_polydata(fluxTriangulation, fullTestName() + ".fluxTriangulation", counter);
            //}

            //EXPECT_LE(fluxTriVolumeError, EPSILON)
                //<< "theta = " << theta << "\n"
                //<< "s = " << s << "\n"
                //<< "displacementX = " << displacementX << "\n"
                //<< "displacementX = " << displacementZ << "\n"
                //<< "correct volume = " << correctVolume << "\n" 
                //<< "normalized flux triangulation volume error = " 
                //<< fluxTriVolumeError / correctVolume << "\n";
            //++counter; 
        //}
    //}
//}


// Diagonal mean point selected based on the edge height to triangulate a non-convex polygon. 
//TEST(triangulation, edgeTriangulation)
//{
    //double edgeLength = 1; 
    //const unsigned int N = 50; 

    //unsigned int counter = 0;
    //for (unsigned int I = 0; I < N + 1; ++I)
    //{
        //double s = I * 0.5 * (edgeLength / N);
        //for (unsigned int J = 0; J < N + 1; ++J)
        //{
            //double theta = (-0.5 * M_PI) + (J * (M_PI / N)); 

            //pointVector polygon = {
                //point(0,0,0), 
                //point(edgeLength,0,0), 
                //point(
                    //edgeLength + s * edgeLength * Foam::sin(theta), 
                    //s * edgeLength * Foam::cos(theta), 
                    //0 
                //), 
                //point(s * edgeLength,0,0)
            //};

            //auto a = polygon[2] - polygon[1];
            //auto b = polygon[3] - polygon[1]; 
            //double correctArea = 0.5 * mag(a ^ b);    

            //double polygonArea = area<edge_triangulation>(polygon); 
            //double areaError = mag(correctArea - polygonArea);  

            //ASSERT_LE(areaError, EPSILON)
                //<< "s = " << s << "\n" 
                //<< "theta = " << theta << "\n" 
                //<< "relative error = " << areaError / correctArea << "\n";  

            //write_vtk_polydata(
                    //edge_triangulate<triangleVector>(polygon), fullTestName() + ".triangulation", counter);  
            //++counter;
        //}
    //}
//}

// Triangulation point must be the same irrespective of the polygon orientation. 
//TEST(triangulation, edgeTriangulationOrientation)
//{
    //// A non-planar polygon
    //pointVector polygon = {point(0,0,0), point(1,0,0), point(1,0.1,1), point(0,-0.1,1)};
    //pointVector reversedPolygon = {point(0,0,0), point(1,0,0), point(1,0.1,1), point(0,-0.1,1)};

    //auto bTri = barycentric_triangulate<triangleVector>(polygon); 
    //write_vtk_polydata(bTri, fullTestName() + ".barycentricTriangulation.vtk"); 

    //auto bTriReversed = barycentric_triangulate<triangleVector>(reversedPolygon); 
    //write_vtk_polydata(bTriReversed, fullTestName() + ".barycentricTriangulationReversed.vtk"); 

    //auto eTri = edge_triangulate<triangleVector>(polygon); 
    //write_vtk_polydata(eTri, fullTestName() + ".edgeTriangulation.vtk"); 

    //auto eTriReversed = edge_triangulate<triangleVector>(reversedPolygon); 
    //write_vtk_polydata(eTriReversed, fullTestName() + ".edgeTriangulationReversed.vtk"); 

    //ASSERT_LE(mag(eTri[0][0] - eTriReversed[0][0]), SMALL); 
    //ASSERT_LE(mag(bTri[0][0] - bTriReversed[0][0]), SMALL); 
//}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

int main(int argc, char *argv[])
{
    ::testing::InitGoogleTest(&argc, argv);
    return RUN_ALL_TESTS();

    return 0;
}


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