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

    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
    geophaseTestInterfacePositioning 

Description
    Interface positioning tests.  

    // Comparison with: 
    Chen, X., & Zhang, X. (2019). A predicted-Newton’s method for solving the
    interface positioning equation in the MoF method on general polyhedrons.
    Journal of Computational Physics, 384, 60–76.
    https://doi.org/10.1016/j.jcp.2018.12.038 
    
Authors
    Tomislav Maric maric@mma.tu-darmstadt.de
    
\*---------------------------------------------------------------------------*/

#include "IterativePositioner.hpp"

#include "Make.hpp"
#include "gtest/gtest.h"
#include "testUtilities.hpp"
#include "WriteVtkPolyData.hpp"
#include "Triangulate.hpp"

#include <fstream>
#include <vector>
#include <algorithm>
#include <chrono>
#include <sstream>
#include <filesystem>

std::vector<double> linspace(double start, double end, unsigned int size)
{
    std::vector<double> result(size, start);
 
    double h = (end - start) / (size - 1); 
    for (unsigned int i = 0; i < size; ++i)
       result[i] = start + h*i; 
    result.back() = end;
 
    return result;
}

static const double epsilonR = 1e-09; // Reconstruction tolerance.

using namespace geophase; 

template<typename Range> 
void testCubicSplinePolynomial(Range const& alphaStars)
{
    for (auto alphaStar : alphaStars)
    {
        // Test hermite interpolation.
        cubicHermitePolynomial hermitePoly(
            0, 1,
            0,1-alphaStar,0,
            1,-alphaStar,0,
            1-alphaStar
        );
        hermitePoly.findRoot();

        ASSERT_LE(
            std::abs(hermitePoly.rootValue()),
            1e-12
        )
            << "root value = " << alphaStar << "\n"
            << "hermite root = " << hermitePoly.root() << "\n"
            << "hermite root value = " 
            << hermitePoly.rootValue() << "\n"; 

    }
}

// Test function that models a fill-level function shifted by 
// alpha along the y-axis:
// f(s0) = 1 - alpha, df(s0) = 0, f(s1) = -alpha, df(s1) = 0.
// cosFillLevel(s*) = 0 gives alpha=alpha* 
struct cosFillLevel  
{
    double alphaStar_ = 0; 

    cosFillLevel (double alphaStar)
        : 
            alphaStar_(alphaStar)
    {}

    double value(double x) const
    {
        return -alphaStar_ - 0.5*cos(M_PI*x) + 0.5;
    }

    double f(double x) const
    {
        return value(x);
    }

    double df(double x) const
    {
        return 0.5*M_PI*sin(M_PI*x);
    }

    double root() const
    {
        double r1 = 2.0 - 0.318309886183791*acos(1.0 - 2.0*alphaStar_);
        double r2 = 0.318309886183791*acos(1.0 - 2.0*alphaStar_);

        if ((r1 >=  0) && (r1 <= 1))
            return r1;

        if ((r2 >=  0) && (r2 <= 1))
            return r2;

        assert(false);

        return std::numeric_limits<double>::max();
    }
};

template<typename Range> 
void testCubicFillLevelPolynomial(Range const& alphaStars)
{
    for (auto alphaStar : alphaStars)
    {
        cosFillLevel func(alphaStar);

        double fs2 = func.f(alphaStar); 
        double dfs2 = func.df(alphaStar);

        // Test cubic polynomial tri-point interpolation.
        cubicFillLevelPolynomial poly(
            0, 1, // Interval containing the root.
            0,func.f(0), // s0, fs0
            1,func.f(1), // s1, fs1
            alphaStar, func.f(alphaStar), // s2, fs2
            func.df(alphaStar), // dfs2
            1-alphaStar // Secant first guess.
        );

        // Test the value and and derivative of the polynomial 
        // at s2. 
        double pfs2 = poly.f(alphaStar); 
        double pdfs2 = poly.df(alphaStar);  

        double fError = std::abs(fs2 - pfs2); 
        double dfError = std::abs(dfs2 - pdfs2); 

        ASSERT_LE(fError, 1e-12)
            << "alphaStar = " << alphaStar << "\n"
            << "f(0) = " << func.f(0) << "\n"
            << "p(0) = " << poly.f(0) << "\n"
            << "f(1) = " << func.f(1) << "\n"
            << "p(1) = " << poly.f(1) << "\n"
            << "f(alphaStar) = " << func.f(alphaStar) << "\n"
            << "p(alphaStar) = " << poly.f(alphaStar) << "\n"
            << "df(alphaStar) = " << func.df(alphaStar) << "\n"
            << "dp(alphaStar) = " << poly.df(alphaStar) << "\n"
            << "p = " << poly << "\n";

        ASSERT_LE(dfError, 1e-12)
            << "alphaStar = " << alphaStar << "\n"
            << "f(0) = " << func.f(0) << "\n"
            << "p(0) = " << poly.f(0) << "\n"
            << "f(1) = " << func.f(1) << "\n"
            << "p(1) = " << poly.f(1) << "\n"
            << "f(alphaStar) = " << func.f(alphaStar) << "\n"
            << "p(alphaStar) = " << poly.f(alphaStar) << "\n"
            << "df(alphaStar) = " << func.df(alphaStar) << "\n"
            << "dp(alphaStar) = " << poly.df(alphaStar) << "\n"
            << "p = " << poly << "\n";

        poly.findRoot(); 
        double root = poly.root();
        ASSERT_LE(std::abs(poly.value(root)), 1e-12)
            << "alphaStar = " << alphaStar << "\n"
            << "f(0) = " << func.f(0) << "\n"
            << "p(0) = " << poly.f(0) << "\n"
            << "f(1) = " << func.f(1) << "\n"
            << "p(1) = " << poly.f(1) << "\n"
            << "f(alphaStar) = " << func.f(alphaStar) << "\n"
            << "p(alphaStar) = " << poly.f(alphaStar) << "\n"
            << "df(alphaStar) = " << func.df(alphaStar) << "\n"
            << "dp(alphaStar) = " << poly.df(alphaStar) << "\n"
            << "p = " << poly << "\n";
    }
}


template<typename Range> 
void testQuinticSplinePolynomial(Range const& alphaStars)
{
    for (auto alphaStar : alphaStars)
    {
        cosFillLevel func(alphaStar);

        double s0 = 0; 
        double fs0 = func.f(s0); 
        double dfs0 = func.df(s0);

        double s2 = alphaStar; 
        double fs2 = func.f(s2); 
        double dfs2 = func.df(s2);

        double s1 = 1; 
        double fs1 = func.f(s1); 
        double dfs1 = func.df(s1);

        // Test fifth-order polynomial interpolation. 
        quinticSplinePolynomial poly(
            s0, s1,
            s0, fs0, dfs0,
            s1, fs1, dfs1, 
            s2, fs2, dfs2, 
            s1 - alphaStar // Secant initial guess.
        ); 

        // Test the value and and derivative of the polynomial 
        // at s2. 
        double pfs2 = poly.f(s2); 
        double pdfs2 = poly.df(s2);  

        double fs2Error = std::abs(fs2 - pfs2); 
        double dfs2Error = std::abs(dfs2 - pdfs2); 

        ASSERT_LE(fs2Error, 1e-10)
            << "alphaStar = " << alphaStar << "\n"
            << "f(0) = " << func.f(0) << "\n"
            << "p(0) = " << poly.f(0) << "\n"
            << "f(1) = " << func.f(1) << "\n"
            << "p(1) = " << poly.f(1) << "\n"
            << "f(alphaStar) = " << func.f(alphaStar) << "\n"
            << "p(alphaStar) = " << poly.f(alphaStar) << "\n"
            << "df(alphaStar) = " << func.df(alphaStar) << "\n"
            << "dp(alphaStar) = " << poly.df(alphaStar) << "\n"
            << "p = " << poly << "\n";

        ASSERT_LE(dfs2Error, 1e-10)
            << "alphaStar = " << alphaStar << "\n"
            << "f(0) = " << func.f(0) << "\n"
            << "p(0) = " << poly.f(0) << "\n"
            << "f(1) = " << func.f(1) << "\n"
            << "p(1) = " << poly.f(1) << "\n"
            << "f(alphaStar) = " << func.f(alphaStar) << "\n"
            << "p(alphaStar) = " << poly.f(alphaStar) << "\n"
            << "df(alphaStar) = " << func.df(alphaStar) << "\n"
            << "dp(alphaStar) = " << poly.df(alphaStar) << "\n"
            << "p = " << poly << "\n";

        poly.findRoot();
        double root = poly.root(); 
        ASSERT_LE(std::abs(poly.value(root)), 1e-12)
            << "alphaStar = " << alphaStar << "\n"
            << "f(0) = " << func.f(0) << "\n"
            << "p(0) = " << poly.f(0) << "\n"
            << "f(1) = " << func.f(1) << "\n"
            << "p(1) = " << poly.f(1) << "\n"
            << "f(alphaStar) = " << func.f(alphaStar) << "\n"
            << "p(alphaStar) = " << poly.f(alphaStar) << "\n"
            << "df(alphaStar) = " << func.df(alphaStar) << "\n"
            << "dp(alphaStar) = " << poly.df(alphaStar) << "\n"
            << "p = " << poly << "\n";
    }
}


class positioningParams 
{
    std::vector<double> alphas_; 
    std::vector<double> thetas_; 
    std::vector<double> phis_; 

    public:

        positioningParams(
            double alpha0 = 0, 
            double alpha1 = 1, 
            int nAlphas = 20,  
            double theta0 = 0, 
            double theta1 = M_PI,
            int nThetas = 20,
            double phi0 = 0, 
            double phi1 = 2*M_PI, 
            int nPhis = 40 
        )
            : 
                alphas_(linspace(alpha0, alpha1, nAlphas)), 
                thetas_(linspace(theta0, theta1, nThetas)), 
                phis_(linspace(phi0, phi1, nPhis))
        {}

        positioningParams(const std::vector<double>& alphas, int nThetas, int nPhis)
            :
                alphas_(std::move(alphas)), 
                thetas_(linspace(EPSILON, M_PI-EPSILON, nThetas)), 
                phis_(linspace(EPSILON, 2*M_PI-EPSILON, nPhis))
        {}

        positioningParams(
            std::vector<double>&& alphas, 
            std::vector<double>&& thetas, 
            std::vector<double>&& phis
        )
            :
                alphas_(std::move(alphas)), 
                thetas_(std::move(thetas)), 
                phis_(std::move(phis))
        {}

        const auto& alphas() const
        {
            return alphas_; 
        }

        const auto& thetas() const
        {
            return thetas_; 
        }

        const auto& phis() const
        {
            return phis_; 
        }
};


template<typename Polyhedron, typename ParamSet, typename Positioner>
void test_positioning(
    Polyhedron const& testPoly, 
    ParamSet const& testParams, 
    Positioner& positioner
)
{
    std::ofstream errorFile(fullTestName() + ".csv");
    errorFile 
        << "#TEST_NAME,POSITIONER,POSITIONING_TOLERANCE,MAX_ITERATIONS\n" 
        << "#" << fullTestName() << "," << positioner.name() << ","
        << positioner.tolerance() << "," << positioner.maxIterations() << "\n";
    errorFile << "ALPHA_STAR,THETA,PHI,ITERATIONS,ROOT_VALUE,ROOT_FINDING_MICROSECONDS,"
        << "POLYNOMIAL_INTERPOLATION_MICROSECONDS,GEOMETRICAL_OPERATIONS_MICROSECONDS,"
        << "CPU_TIME_MICROSECONDS\n"; 

    unsigned long int totalCPUtime = 0;
    unsigned long int totalIterations = 0;
    unsigned long int totalTests = 0; 
#ifdef TIMING
    unsigned long int totalRootTime = 0;
    unsigned long int totalPolyTime = 0;
    unsigned long int totalGeomTime = 0;
#endif
    vector normal(0.,0.,0.);
    for (auto alphaStar : testParams.alphas())
    {
        for (auto theta : testParams.thetas())
        {
            for (auto phi : testParams.phis())
            {
                ++totalTests; 

                normal = vector(
                    sin(theta)*cos(phi), 
                    sin(theta)*sin(phi), 
                    cos(theta)
                );

                auto time0 = std::chrono::steady_clock::now();
                positioner.positionInterface(
                    alphaStar, 
                    normal, 
                    extreme_points(testPoly, normal)
                ); 
                auto time1 = std::chrono::steady_clock::now();
                auto cpuTimeDelta = 
                    std::chrono::duration_cast<std::chrono::microseconds>(time1-time0).count();

                EXPECT_LE(
                    positioner.rootValue(), 
                    positioner.tolerance()
                )
                    << "alphaStar = " << alphaStar << "\n"
                    << "normal = " << normal << "\n" 
                    << "theta = " << theta << "\n"
                    << "phi = " << phi << "\n";

                errorFile << std::setprecision(20) 
                    << alphaStar << "," 
                    << theta << "," 
                    << phi << "," 
                    << positioner.iterations() << "," 
                    << positioner.rootValue() << ","
#ifdef TIMING
                    << positioner.rootTime() << "," 
                    << positioner.polyTime() << "," 
                    << positioner.geomTime() << "," 
#endif
                    << cpuTimeDelta << "\n"; 

#ifdef TIMING
                totalRootTime += positioner.rootTime(); 
                totalPolyTime += positioner.polyTime(); 
                totalGeomTime += positioner.geomTime(); 
#endif
                totalCPUtime  += cpuTimeDelta;
                totalIterations += positioner.iterations(); 
            }
        }
    }

    const std::string fileName(
        testInfo()->test_case_name() + 
        std::string("_TOTAL_DATA.csv")
    );

    std::stringstream totalsStream; 
    totalsStream << testInfo()->name() << "," 
#ifdef TIMING
        << totalRootTime << ","
        << totalPolyTime << ","
        << totalGeomTime << "," 
        << double(totalRootTime) / totalIterations << ","
        << double(totalPolyTime) / totalIterations << ","
        << double(totalGeomTime) / totalIterations << "," 
#endif 
        << totalCPUtime << "," 
        << totalIterations << ","
        << double(totalCPUtime) / totalTests << ","
        << double(totalIterations) / totalTests << "\n";

    if (!std::filesystem::exists(fileName))
    {
        std::ofstream totalsFile(fileName);
        totalsFile << "TEST,"
#ifdef TIMING
            << "ROOT_FINDING_MICROSECONDS_TOTAL,"
            << "POLYNOMIAL_INTERPOLATION_MICROSECONDS_TOTAL,"
            << "GEOMETRICAL_INTERSECTION_MICROSECONDS_TOTAL,"
            << "ROOT_FINDING_MICROSECONDS_AVERAGE,"
            << "POLYNOMIAL_INTERPOLATION_MICROSECONDS_AVERAGE,"
            << "GEOMETRICAL_INTERSECTION_MICROSECONDS_AVERAGE,"
#endif
            << "CPU_TIME_MICROSECONDS_TOTAL,ITERATIONS_TOTAL,"
            << "AVERAGE_CPU_TIME_MICROSECONDS,AVERAGE_ITERATIONS\n"; 
        totalsFile << totalsStream.str();
    } else
    {
        std::ofstream totalsFile(fileName, std::ios::app);
        totalsFile << totalsStream.str(); 
    }
}

TEST(CUBIC_SPLINE_POLYNOMIAL, ROOTS)
{
    testCubicSplinePolynomial(linspace(epsilonR, 1e-03, 1e04)); 
    testCubicSplinePolynomial(linspace(epsilonR, 1-epsilonR, 1e04)); 
    testCubicSplinePolynomial(linspace(1-1e-03, 1-epsilonR, 1e04)); 
}

TEST(CUBIC_FILL_LEVEL_POLYNOMIAL, ROOTS)
{
    testCubicFillLevelPolynomial(linspace(epsilonR, 1e-03, 1e04)); 
    testCubicFillLevelPolynomial(linspace(epsilonR, 1-epsilonR, 1e04)); 
    testCubicFillLevelPolynomial(linspace(1-1e-03, 1-epsilonR, 1e04)); 
}

TEST(QUINTIC_SPLINE_POLYNOMIAL, ROOTS)
{
    testQuinticSplinePolynomial(linspace(epsilonR, 1e-03, 1e04)); 
    testQuinticSplinePolynomial(linspace(epsilonR, 1-epsilonR, 1e04)); 
    testQuinticSplinePolynomial(linspace(1-1e-03, 1-epsilonR, 1e04)); 
}

TEST(INTERSECTION_VOLUME, CUBE)
{
    const auto testCube = make_unit_cube();

    auto tetrahedra =  centroid_vol_triangulate(testCube); 
    auto volumes = volumes_by_surf_tri(tetrahedra);
    halfspace hspace(vector{0.5,0.5,0.9},vector{0.,0.,1.});

    double intersectedVolume = intersect_volume(
        tetrahedra, 
        volumes,
        hspace
    );

    auto triangles = centroid_triangulate(make_unit_cube()); 
    ASSERT_LE(std::abs(intersectedVolume - 0.1), 2*EPSILON); 

    /* NEW Volume intersection underway. TM.
    // Intersect cube as a triangulated surface when it is completely inside.
    halfspace cube_inside_hspace(vector{0.,0.,0.},vector{0.,0.,1.});
    double cube_volume = intersect_volume_area(triangles, cube_inside_hspace);
    ASSERT_LE(std::abs(cube_volume - 1.0), 2*EPSILON)
        << "cube_volume = " << cube_volume << std::endl;

    std::cout << "***********************************" << std::endl;

    // Halve the cube as a triangulated surface. 
    halfspace halved_cube_hspace(vector{0.5,0.5,0.5},vector{0.,0.,1.});
    double halved_cube_vol = intersect_volume_area(triangles, halved_cube_hspace);
    ASSERT_LE(std::abs(halved_cube_vol - 0.5), 2*EPSILON)
        << "halved_cube_vol = " << halved_cube_vol << std::endl;

    std::cout << "***********************************" << std::endl;

    // Intersect the cube as a triangulated surface at 0.1. 
    halfspace cut_cube_hspace(vector{0.0,0.0,0.25},vector{0.,0.,-1.});
    double cut_cube_vol = intersect_volume_area(triangles, cut_cube_hspace);
    std::cout << cut_cube_vol << std::endl;
    ASSERT_LE(std::abs(cut_cube_vol - 0.25), 2*EPSILON); 
    */
}

// CONVERGENCE: Large parameter sets used for testing numerical convergence. 
class CONVERGENCE_TEST_POSITIONING : public ::testing::Test {

 protected:

  void SetUp() override {

    std::vector<double> leftAlphas { 
        1e-09, 1e-08, 1e-07, 1e-06, 
        1e-05, 1e-04, 1e-03, 1e-02,
    };
    std::vector<double> middleAlphas = linspace(1e-03,1-1e-03, 50); 

    std::vector<double> rightAlphas {
        1-1e-02, 1-1e-03, 1-1e-04, 
        1-1e-05, 1-1e-06, 1-1e-07, 
        1-1e-08, 1-1e-09, 
    };

    std::vector<double> testAlphas; 

    testAlphas.insert(
        testAlphas.end(), 
        leftAlphas.begin(), 
        leftAlphas.end()
    );

    testAlphas.insert(
        testAlphas.end(), 
        middleAlphas.begin(), 
        middleAlphas.end()
    );

    testAlphas.insert(
        testAlphas.end(), 
        rightAlphas.begin(), 
        rightAlphas.end()
    );

    params_ = positioningParams(
        testAlphas, // alpha values
        41, // Ntheta 
        21// Ntheta
    );

  }

  positioningParams params_; 
};

// UNIT: Small tests used in development for fast / integration testing.
class UNIT_TEST_POSITIONING : public ::testing::Test {
    
 protected:

  void SetUp() override {

    params_ = positioningParams({
        1e-09, 1e-08, 1e-07, 1e-06, 
        1e-05, 1e-04, 1e-03, 1e-02,
        0.1,0.25,0.5,0.75,0.9,
        1-1e-02, 1-1e-03, 1-1e-04, 
        1-1e-05, 1-1e-06, 1-1e-07, 
        1-1e-08, 1-1e-09, 
        }, // alpha values 
        10, // Ntheta
        20 // Nphi
    );
  }

  positioningParams params_; 
};

// DEBUG: Smoke tests, very small datasets used for debugging. 
class DEBUG_TEST_POSITIONING : public ::testing::Test {

 protected:

  void SetUp() override {

    params_ = positioningParams({
        1e-09, 1e-07,  1e-05, 1e-03, 
        0.1,0.5,0.75,0.9,
        1-1e-03, 1-1e-05, 1-1e-07, 1-1e-09, 
        }, // alpha values 
        5, // Ntheta
        10 // Nphi
    );
  }

  positioningParams params_; 
};

// **************** BEGIN BISECTION **************** //
// Bisection is too slow for practicle use, so only tet and 
// cube unit tests are performed.

template<typename Polyhedron, typename Parameter>
void test_bisection_positioner(
    Polyhedron const& polyhedron, 
    Parameter const& params
)
{
    bisectionPositioner bisection(
        centroid_vol_triangulate(polyhedron), 
        volume_by_surf_tri(polyhedron)
    );

    test_positioning(polyhedron, params, bisection);
}

TEST_F(UNIT_TEST_POSITIONING, TETRAHEDRON_BISECTION)
{
    const auto testTetrahedron = make_unit_tetrahedron();
    test_bisection_positioner(testTetrahedron, params_);
}

TEST_F(UNIT_TEST_POSITIONING, CUBE_BISECTION)
{
    const auto testCube = make_unit_cube();
    test_bisection_positioner(testCube, params_);
}

// **************** END BISECTION **************** //

// **************** BEGIN NEWTON CUBIC SPLINE **************** //

template<typename Polyhedron, typename Parameter>
void test_newton_positioner(Polyhedron const& polyhedron, Parameter const& params)
{
    write_vtk(polyhedron, fullTestName() + ".vtk");

    newtonCubicSplinePositioner newtonCubicSpline(
        centroid_vol_triangulate(polyhedron), 
        volume_by_surf_tri(polyhedron)
    );

    test_positioning(polyhedron, params, newtonCubicSpline);
}

// **************** UNIT TESTS **************** //

TEST(UNIT_VOL_FRACTION, ENDO_DODECAHEDRON)
{
    const auto dodecahedron = make_unit_endo_dodecahedron();
    auto tetrahedra = centroid_vol_triangulate(dodecahedron); 
    double surfVolume = volume_by_surf_tri(dodecahedron);
    auto tetrasVolumes = volumes_by_surf_tri(tetrahedra); 
    double tetrasVolume = volume_by_surf_tri(tetrahedra); 
    double volVolume = volume_by_vol_tri(dodecahedron);

    vector normal (0., 0., 1.);
    auto extremePoints = extreme_points(dodecahedron, normal);
    double s10 = dot(extremePoints.second - extremePoints.first, normal); 

    std::ofstream alphaFile(fullTestName() + ".csv");
    alphaFile << "#TEST_NAME " << fullTestName() << "\n" 
        << "S,ALPHA_STAR,ALPHA_VOLUME,SURF_TRI_VOLUME,TETRAS_VOLUME,VOL_TRI_VOLUME\n"; 

    unsigned int NIT = 100; 
    for (unsigned int iter = 0; iter <= NIT; ++iter)
    {
        double s = iter * s10 / NIT; 
        halfspace hspace(extremePoints.first + s*normal, normal); 
        auto volArea = intersect_volume_area(tetrahedra, tetrasVolumes, hspace); 
        alphaFile << std::setprecision(20) << s << "," 
            << volArea.volume_ / tetrasVolume << "," << volArea.volume_ << ","
            << surfVolume << "," << tetrasVolume << "," << volVolume << "\n";

        double surfAlpha = volArea.volume_ / surfVolume; 
        double volAlpha = volArea.volume_ / volVolume; 
        double tetAlpha = volArea.volume_ / tetrasVolume;

        ASSERT_LE(std::abs(surfAlpha - volAlpha), 2e-15);
        ASSERT_LE(std::abs(volAlpha - tetAlpha), 2e-15);

        // Visualize the triangulation / halfspace intersection.
        vtkPolyDataOStream vtks(vtk_file_name(fullTestName(), iter)); 
        for (const auto& tetrahedron : tetrahedra)
        {
            auto intersection = intersect_tolerance(tetrahedron, hspace); 
            vtks << intersection.polyhedron();
        }
    }
}

TEST_F(UNIT_TEST_POSITIONING, TETRAHEDRON_NEWTON_CUBIC_SPLINE)
{
    const auto testTetrahedron = make_unit_tetrahedron();
    test_newton_positioner(testTetrahedron, params_);
}

TEST_F(UNIT_TEST_POSITIONING, CUBE_NEWTON_CUBIC_SPLINE)
{
    const auto testCube = make_unit_cube();
    test_newton_positioner(testCube, params_);
}

TEST_F(UNIT_TEST_POSITIONING, DODECAHEDRON_NEWTON_CUBIC_SPLINE)
{
    const auto testDodecahedron = make_unit_dodecahedron();
    test_newton_positioner(testDodecahedron, params_);
}

TEST_F(UNIT_TEST_POSITIONING, NON_PLANAR_DODECAHEDRON_NEWTON_CUBIC_SPLINE)
{
    const auto testDodecahedron = make_unit_dodecahedron(0.2);
    test_newton_positioner(testDodecahedron, params_);
}

TEST_F(UNIT_TEST_POSITIONING, ENDO_DODECAHEDRON_NEWTON_CUBIC_SPLINE)
{
    const auto testDodecahedron = make_unit_endo_dodecahedron();
    test_newton_positioner(testDodecahedron, params_);
}

// **************** CONVERGENCE TESTS **************** //


TEST_F(CONVERGENCE_TEST_POSITIONING, TETRAHEDRON_NEWTON_CUBIC_SPLINE)
{
    const auto testTetrahedron = make_unit_tetrahedron();
    test_newton_positioner(testTetrahedron, params_);
}

TEST_F(CONVERGENCE_TEST_POSITIONING, CUBE_NEWTON_CUBIC_SPLINE)
{
    const auto testCube = make_unit_cube();
    test_newton_positioner(testCube, params_);
}

TEST_F(CONVERGENCE_TEST_POSITIONING, DODECAHEDRON_NEWTON_CUBIC_SPLINE)
{
    const auto testDodecahedron = make_unit_dodecahedron();
    test_newton_positioner(testDodecahedron, params_);
}

TEST_F(CONVERGENCE_TEST_POSITIONING, NON_PLANAR_DODECAHEDRON_NEWTON_CUBIC_SPLINE)
{
    const auto testDodecahedron = make_unit_dodecahedron(0.2);
    test_newton_positioner(testDodecahedron, params_);
}

TEST_F(CONVERGENCE_TEST_POSITIONING, ENDO_DODECAHEDRON_NEWTON_CUBIC_SPLINE)
{
    const auto testDodecahedron = make_unit_endo_dodecahedron();
    test_newton_positioner(testDodecahedron, params_);
}

// **************** END NEWTON CUBIC SPLINE **************** //

// **************** BEGIN SUCCESSIVE CUBIC SPLINE **************** //

template<typename Polyhedron, typename Parameter>
void test_cubic_spline_positioner(
    Polyhedron const& polyhedron, 
    Parameter const& params
)
{
    cubicSplinePositioner cubicSpline(
        centroid_vol_triangulate(polyhedron), 
        volume_by_surf_tri(polyhedron)
    );

    test_positioning(polyhedron, params, cubicSpline);
}

// **************** UNIT TESTS **************** //

TEST_F(UNIT_TEST_POSITIONING, TETRAHEDRON_CONSECUTIVE_CUBIC_SPLINE)
{
    const auto testTetrahedron = make_unit_tetrahedron();
    test_cubic_spline_positioner(testTetrahedron, params_);
}

TEST_F(UNIT_TEST_POSITIONING, CUBE_CONSECUTIVE_CUBIC_SPLINE)
{
    const auto testCube = make_unit_cube();
    test_cubic_spline_positioner(testCube, params_);
}

TEST_F(UNIT_TEST_POSITIONING, DODECAHEDRON_CONSECUTIVE_CUBIC_SPLINE)
{
    const auto testDodecahedron = make_unit_dodecahedron();
    test_cubic_spline_positioner(testDodecahedron, params_);
}

TEST_F(UNIT_TEST_POSITIONING, NON_PLANAR_DODECAHEDRON_CONSECUTIVE_CUBIC_SPLINE)
{
    const auto testDodecahedron = make_unit_dodecahedron(0.2);
    test_cubic_spline_positioner(testDodecahedron, params_);
}

TEST_F(UNIT_TEST_POSITIONING, ENDO_DODECAHEDRON_CONSECUTIVE_CUBIC_SPLINE)
{
    const auto testDodecahedron = make_unit_endo_dodecahedron();
    test_cubic_spline_positioner(testDodecahedron, params_);
}

// **************** CONVERGENCE TESTS **************** //

TEST_F(CONVERGENCE_TEST_POSITIONING, TETRAHEDRON_CONSECUTIVE_CUBIC_SPLINE)
{
    const auto testTetrahedron = make_unit_tetrahedron();
    test_cubic_spline_positioner(testTetrahedron, params_);
}

TEST_F(CONVERGENCE_TEST_POSITIONING, CUBE_CONSECUTIVE_CUBIC_SPLINE)
{
    const auto testCube = make_unit_cube();
    test_cubic_spline_positioner(testCube, params_);
}

TEST_F(CONVERGENCE_TEST_POSITIONING, DODECAHEDRON_CONSECUTIVE_CUBIC_SPLINE)
{
    const auto testDodecahedron = make_unit_dodecahedron();
    test_cubic_spline_positioner(testDodecahedron, params_);
}

TEST_F(CONVERGENCE_TEST_POSITIONING, NON_PLANAR_DODECAHEDRON_CONSECUTIVE_CUBIC_SPLINE)
{
    const auto testDodecahedron = make_unit_dodecahedron(0.2);
    test_cubic_spline_positioner(testDodecahedron, params_);
}

TEST_F(CONVERGENCE_TEST_POSITIONING, ENDO_DODECAHEDRON_CONSECUTIVE_CUBIC_SPLINE)
{
    const auto testDodecahedron = make_unit_endo_dodecahedron();
    test_cubic_spline_positioner(testDodecahedron, params_);
}

// **************** END SUCCESSIVE CUBIC SPLINE **************** //

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

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

    return 0;
}


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