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

    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
    geophaseTestIntersection 

Description
    Test different intersection algorithms with different models.    
   
Author
    Tomislav Maric maric@mma.tu-darmstadt.de

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

#include "gtest/gtest.h" 
#include "PolygonIntersection.hpp"
#include "Halfspace.hpp"
#include "Intersect.hpp"
#include "Equal.hpp"
#include "Area.hpp"
#include "Volume.hpp"
#include "Make.hpp"
#include "WriteVtkPolyData.hpp"

#include "testUtilities.hpp"
#include "Perturb.hpp"

using namespace geophase;

// ******************************************       BEGIN Polygon 

TEST(polygonHalfspaceIntersection, point)
{
    auto doTest = [](const double LENGTH)
    {
        vectorPolygon poly = {
            point(0.,0.,0.), 
            point(LENGTH,0.,0.),
            point(LENGTH,LENGTH,0.), 
            point(0.,LENGTH,0.)
        };

        halfspace h (point(LENGTH, LENGTH, 0.), 
                     vector(1.,1.,0.)); 

        vectorPolygonIntersection isect;

        intersect_tolerance(isect, poly, h);

        const auto& iPoly = isect.polygon();

        ASSERT_TRUE(iPoly.size() == 1);
        ASSERT_TRUE(iPoly.front() == h.position());
    };

    doTest(1.0); 
    doTest(0.5); 
    doTest(0.1); 
    doTest(1e-03); 
    doTest(1e-06); 
    doTest(1e-10); 
    // Minimal working length for the intersect_tolerance algorithm.
    doTest(5 * std::numeric_limits<double>::epsilon()); 
}

TEST(polygonHalfspaceIntersection, edgeOnly)
{
    auto doTest = [](const double LENGTH)
    {
        vectorPolygon poly = {
            point(0.,0.,0.), 
            point(LENGTH,0.,0.),
            point(LENGTH,LENGTH,0.), 
            point(0.,LENGTH,0.)
        };

        halfspace h (point(LENGTH, LENGTH, 0.), 
                     vector(1.,0.,0.)); 

        vectorPolygonIntersection isect;
        intersect_tolerance(isect, poly, h);
        const auto& isectPoly = isect.polygon();

        ASSERT_TRUE(isectPoly.size() == 2);
        ASSERT_TRUE(isectPoly.front() == poly[1]); 
        ASSERT_TRUE(isectPoly.back() == poly[2]); 
    };

    doTest(1.0); 
    doTest(0.5); 
    doTest(0.1); 
    doTest(1e-03); 
    doTest(1e-06); 
    doTest(1e-10); 
    // Minimal working length for the intersect_tolerance algorithm.
    doTest(5 * std::numeric_limits<double>::epsilon()); 
}

TEST(polygonHalfspaceIntersection, edgeFullPolygon)
{
    auto doTest = [](const double LENGTH)
    {
        vectorPolygon poly = {
            point(0.,0.,0.), 
            point(LENGTH,0.,0.),
            point(LENGTH,LENGTH,0.), 
            point(0.,LENGTH,0.)
        };

        halfspace h (point(LENGTH, LENGTH, 0.), 
                     vector(-1.,0.,0.)); 

        vectorPolygonIntersection isect;
        intersect_tolerance(isect, poly, h);
        const auto& isectPoly = isect.polygon();

        ASSERT_TRUE(isectPoly.size() == 4);
        ASSERT_TRUE(isectPoly == poly); 
    };

    doTest(1.0); 
    doTest(0.5); 
    doTest(0.1); 
    doTest(1e-03); 
    doTest(1e-06); 
    doTest(1e-10); 
    // Minimal working length for the intersect_tolerance algorithm.
    doTest(5 * std::numeric_limits<double>::epsilon()); 
}

TEST(polygonHalfspaceIntersection, diagonalCut)
{
    auto doTest = [](const double LENGTH)
    {
        vectorPolygon poly = {
            point(0.,0.,0.), 
            point(LENGTH,0.,0.),
            point(LENGTH,LENGTH,0.), 
            point(0.,LENGTH,0.)
        };

        halfspace h (point(0.5 * LENGTH, 0.5 * LENGTH, 0.), 
                     vector(1.,-1.,0.)); 

        vectorPolygonIntersection isect; 
        intersect_tolerance(isect, poly, h);
        const auto& isectPoly = isect.polygon();

        ASSERT_TRUE(isectPoly.size() == 3)
            << "Original polygon = " << poly << std::endl
            << "Intersection polygon = " << isectPoly << std::endl;

        ASSERT_TRUE(isectPoly[0] == poly[0])
            << "Original polygon = " << poly << std::endl
            << "Intersection polygon = " << isectPoly << std::endl;

        ASSERT_TRUE(isectPoly[1] == poly[1])
            << "Original polygon = " << poly << std::endl
            << "Intersection polygon = " << isectPoly << std::endl;

        ASSERT_TRUE(isectPoly[2] == poly[2])
            << "Original polygon = " << poly << std::endl
            << "Intersection polygon = " << isectPoly << std::endl;

    };

    doTest(1.0); 
    doTest(0.5); 
    doTest(0.1); 
    doTest(1e-03); 
    doTest(1e-06); 
    doTest(1e-10); 
    // Minimal working length for the intersect_tolerance algorithm.
    doTest(5 * std::numeric_limits<double>::epsilon()); 
}

TEST(polygonHalfspaceIntersection, halvedPolygon)
{
    auto doTest = [](const double LENGTH)
    {
        vectorPolygon poly = {
            point(0.,0.,0.), 
            point(LENGTH,0.,0.),
            point(LENGTH,LENGTH,0.), 
            point(0.,LENGTH,0.)
        };

        halfspace hspace (point(0.5 * LENGTH, 0.5 * LENGTH, 0.), 
                          vector(0.,-1.,0.)); 

        vectorPolygonIntersection isect; 
        intersect_tolerance(isect, poly, hspace);
        const auto& isectPoly = isect .polygon();
        ASSERT_TRUE(isectPoly.size() == 4)
            << "Original polygon = " << poly << std::endl
            << "Intersection polygon = " << isectPoly << std::endl;

        ASSERT_TRUE(isectPoly[0] == poly[0]) 
            << "Original polygon = " << poly << std::endl
            << "Intersection polygon = " << isectPoly << std::endl;

        ASSERT_TRUE(isectPoly[1] == poly[1])
            << "Original polygon = " << poly << std::endl
            << "Intersection polygon = " << isectPoly << std::endl;

        ASSERT_TRUE(isectPoly[2] == point(LENGTH, 0.5 * LENGTH, 0.))
            << "Original polygon = " << poly << std::endl
            << "Intersection polygon = " << isectPoly << std::endl;

        ASSERT_TRUE(isectPoly[3] == point(0., 0.5 * LENGTH, 0.))
            << "Original polygon = " << poly << std::endl
            << "Intersection polygon = " << isectPoly << std::endl;
    };

    doTest(1.0); 
    doTest(0.5); 
    doTest(0.1); 
    doTest(1e-03); 
    doTest(1e-06); 
    doTest(1e-10); 
    // Minimal working length for the intersect_tolerance algorithm.
    doTest(5 * std::numeric_limits<double>::epsilon()); 
}

// ******************************************       END Polygon 
//
// ******************************************       BEGIN Polyhedron Halfspace 

TEST(polyhedronHalfspaceIntersection, halvedCube)
{
    auto cube = make_unit_cube();
    halfspace hspace(vector(0.5,0.5,0.5), vector(1.,0.,0.)); 

    auto cubeIntersection = intersect_tolerance(cube, hspace); 
    const auto& poly = cubeIntersection.polyhedron(); 

    ASSERT_EQ(volume_by_surf_tri(poly), 0.5); 
}

TEST(polyhedronHalfspaceIntersection, diagonallyHalvedCube)
{
    auto cube = make_unit_cube();
    halfspace hspace(vector(0.5,0.5,0.5), vector(1.,1.,0.)); 

    auto cubeIntersection = intersect_tolerance(cube, hspace); 
    const auto& poly = cubeIntersection.polyhedron(); 

    ASSERT_EQ(volume_by_surf_tri(poly), 0.5); 
}

TEST(polyhedronHalfspaceIntersection, faceIntersected)
{
    auto cube = make_unit_cube();
    halfspace hspace(vector(1.0 - EPSILON, 0., 0.), vector(1.,0.,0.)); 

    auto cubeIntersection = intersect_tolerance(cube, hspace); 
    const auto& poly = cubeIntersection.polyhedron(); 

    ASSERT_TRUE(poly.empty());
}

TEST(polyhedronHalfspaceIntersection, edgeIntersected)
{
    auto cube = make_unit_cube();
    halfspace hspace(vector(1.0 - EPSILON, 0., 0.), vector(1.,-1.,0.)); 

    auto cubeIntersection = intersect_tolerance(cube, hspace); 
    const auto& poly = cubeIntersection.polyhedron(); 

    ASSERT_TRUE(poly.empty());
}

// ******************************************       END Polyhedron Halfspace 


// FIXME: Check and clean up the tests.
//// Global intersection test helping functions.
//template<typename Point, typename Halfspace>
//void assertFalseIntersection(
    //const Point& p1, 
    //const Point& p2, 
    //const Halfspace h1 
//)
//{
    //typedef LineIntersection<Point> lineIntersection; 

    //auto result = intersect<lineIntersection>(p1, p2, h1); 

    //ASSERT_FALSE(result.isValid()); 

    //result = intersect<lineIntersection>(p1, h1, p2); 

    //ASSERT_FALSE(result.isValid()); 

    //result = intersect<lineIntersection>(h1, p1, p2); 

    //ASSERT_FALSE(result.isValid()); 
//};

//template<typename Point, typename Halfspace>
//void assertTrueIntersection(
    //const Point& p1, 
    //const Point& p2, 
    //const Halfspace h1 
//)
//{
    //typedef LineIntersection<Point> lineIntersection; 

    //auto result = intersect<lineIntersection>(p1, p2, h1); 

    //ASSERT_TRUE(result.isValid()); 

    //result = intersect<lineIntersection>(p1, h1, p2); 

    //ASSERT_TRUE(result.isValid()); 

    //result = intersect<lineIntersection>(h1, p1, p2); 

    //ASSERT_TRUE(result.isValid()); 
//};

//template<typename PolygonSequence, typename Point> 
//PolygonSequence build_box_star(Point p1, Point p2, double sweepFactor = 1)
//{
    //PolygonSequence result; 

    //AABBox<point> aabbox(p1, p2); 

    //auto polyhedron = build<PolygonSequence>(aabbox);  

    //typedef typename PolygonSequence::value_type polygonType; 
    //for (const auto& polygon : polyhedron)
    //{
        //auto polygonNormal = normal_area_vec(polygon); 
        //polygonNormal /= mag(polygonNormal); 
        //polygonNormal *= sweepFactor;
        //auto polygonCentroid = centroid(polygon); 
        //auto triangulationPoint = polygonCentroid + polygonNormal; 
        //auto polygonTriangulation = point_triangulate<triangleVector>(polygon, triangulationPoint); 

        //for (const auto& triangle : polygonTriangulation)
        //{
            //polygonType insertedTriangle;  
            //std::copy(triangle.begin(), triangle.end(), std::back_inserter(insertedTriangle));
            //result.push_back(insertedTriangle); 
        //}
    //}
    //return result; 
//}

//TEST(doIntersectPolygonHalfspace, diagonalIntersection)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef Halfspace<point, point> halfspace;

    //polygon poly = {point(0,0,0), point(1,0,0), point(1,1,0), point(0,1,0)};

    //halfspace h (point(0.5,0.5,0), vector(1,1,0)); 

    //ASSERT_TRUE(do_intersect(poly, h)); 

//}

//TEST(doIntersectPolygonHalfspace, regularIntersection)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef Halfspace<point, point> halfspace;

    //polygon poly = {point(0,0,0), point(1,0,0), point(1,1,0), point(0,1,0)};

    //halfspace h (point(0.75,0.75,0), vector(1,1,0)); 

    //ASSERT_TRUE(do_intersect(poly, h)); 
//}

//TEST(doIntersectPolygonHalfspace, insideHalfspace)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef Halfspace<point, point> halfspace;

    //polygon poly = {point(0,0,0), point(1,0,0), point(1,1,0), point(0,1,0)};

    //halfspace h (point(0,0,0), vector(1,1,0)); 

    //ASSERT_TRUE(do_intersect(poly, h)); 
//}

//TEST(doIntersectPolygonHalfspace, pointIntersection)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef Halfspace<point, point> halfspace;

    //polygon poly = {point(0,0,0), point(1,0,0), point(1,1,0), point(0,1,0)};

    //halfspace h (point(1,1,0), vector(1,1,0)); 

    //ASSERT_TRUE(do_intersect(poly, h)); 
//}

//TEST(doIntersectPolyhedronHalfspace, diagonalIntersection)
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef Halfspace<point, point> halfspace;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;

    //auto poly = build<polyhedron>(aabbox(point(0,0,0), point(1,1,1))); 

    //halfspace h (point(0.5,0.5,0), vector(1,1,0)); 

    //ASSERT_TRUE(do_intersect(poly, h)); 
//}

//TEST(doIntersectPolyhedronHalfspace, regularIntersection)
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef Halfspace<point, point> halfspace;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;

    //auto poly = build<polyhedron>(aabbox(point(0,0,0), point(1,1,1))); 

    //halfspace h (point(0.75,0.75,0.75), vector(1,1,1)); 

    //ASSERT_TRUE(do_intersect(poly, h)); 
//}

//TEST(doIntersectPolyhedronHalfspace, insideHalfspace)
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef Halfspace<point, point> halfspace;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;

    //auto poly = build<polyhedron>(aabbox(point(0,0,0), point(1,1,1))); 

    //halfspace h (point(0,0,0), vector(1,1,1)); 

    //ASSERT_TRUE(do_intersect(poly, h)); 
//}

//TEST(doIntersectPolyhedronHalfspace, pointIntersection)
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef Halfspace<point, point> halfspace;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;

    //auto poly = build<polyhedron>(aabbox(point(0,0,0), point(1,1,1))); 

    //halfspace h (point(1,1,1), vector(1,1,1)); 

    //ASSERT_TRUE(do_intersect(poly, h)); 
//}

///// @brief Test axis aligned box intersection for two disjoint boxes
/////
///// @param doIntersectBoxes
///// @param nonOverlappingBoxes
//TEST(doIntersectBoxes, nonOverlappingBoxes)
//{
    //aabbox firstBox (point(0,0,0), point(1,1,1)); 
    //aabbox secondBox (point(3,3,3), point(4,4,4)); 

    //ASSERT_FALSE(do_intersect(firstBox, secondBox)); 
//}

///// @brief Test AABBox intersection between intersecting boxes
//TEST(doIntersectBoxes, intersectingBoxes)
//{
    //aabbox firstBox (point(0,0,0), point(1,1,1));
    //aabbox secondBox (point(0.5,0.5,0.5), point(1.5,1.5,1.5)); 
    //ASSERT_TRUE(do_intersect(firstBox, secondBox)); 
//}

///// @brief Test AABBox intersection between overlapping boxes.
//TEST(doIntersectBoxes, overlappingBoxes)
//{
    //aabbox firstBox (point(0,0,0), point(1,1,1)); 
    //aabbox secondBox (firstBox); 

    //ASSERT_TRUE(do_intersect(firstBox, secondBox)); 
//}

///// @brief Two AABBoxes that touch faces do not intersect.
//TEST(doIntersectBoxes, boxesTouchingFaces)
//{
    //aabbox firstBox (point(0,0,0), point(1,1,1));
    //aabbox secondBox (point(1,0,0), point(2,1,1)); 
    //ASSERT_FALSE(do_intersect(firstBox, secondBox)); 
//}

///// @brief Two AABBoxes that touch edges do not intersect.
//TEST(doIntersectBoxes, boxesTouchingEdges)
//{
    //aabbox firstBox (point(0,0,0), point(1,1,1));
    //aabbox secondBox (point(1,1,0), point(2,2,1)); 
    //ASSERT_FALSE(do_intersect(firstBox, secondBox)); 
//}

///// @brief Two AABBoxes that touch a point do not intersect.
//TEST(doIntersectBoxes, boxesTouchingPoint)
//{
    //aabbox firstBox (point(0,0,0), point(1,1,1));
    //aabbox secondBox (point(1,1,0), point(2,2,1)); 
    //ASSERT_FALSE(do_intersect(firstBox, secondBox)); 
//}

//TEST(twoPointLineIntersection, noIntersection)
//{
    //typedef Halfspace<point, vector> halfspace;

    //point p1(0,0,0), p2(1,0,0);  

    //// No intersection, orthogonal normal.
    //halfspace h1(point(2,0,0), vector(1,0,0));
    //assertFalseIntersection(p1, p2, h1); 

    //// Coplanar edge.
    //h1.setDirection(vector(0, 1, 0));
    //assertFalseIntersection(p1,p2,h1); 

    //// No intersection, but plane intersects line outside of the segment.
    //h1.setPosition(point(0,1,0));
    //h1.setDirection(vector(0.5, 1, 0));
    //assertFalseIntersection(p1,p2,h1); 
//}

//TEST(twoPointLineIntersection, samePointIntersection)
//{
    //typedef Halfspace<point, vector> halfspace;

    //point p1(1,0,0), p2(1,0,0);  

    //halfspace h1(point(1,0,0), vector(1,0,0));

    //assertFalseIntersection(p1, p2, h1); 
//}

//TEST(twoPointLineIntersection, pointIntersection)
//{
    //typedef Halfspace<point, vector> halfspace;

    //point p1(0,0,0), p2(1,0,0);  

    //halfspace h1(point(1,0,0), vector(1,0,0));

    //assertTrueIntersection(p1, p2, h1); 
//}

//TEST(lineSegmentIntersection, halfspaceNoIntersection)
//{
    //typedef Halfspace<point, vector> halfspace;
    //typedef LineSegment<point> lineSegment;

    //lineSegment l1 (point(0,0,0), point(1,0,0));
    //halfspace h1 (point(-1,0,0), vector(-1,0,0));

    //auto result = intersect<point>(l1, h1); 

    //auto resultReverse = intersect<point>(h1, l1); 

    //ASSERT_TRUE(result == resultReverse);
//}

//TEST(lineSegmentIntersection, halfspacePointIntersection)
//{
    //typedef Halfspace<point, vector> halfspace;
    //typedef LineSegment<point> lineSegment;

    //lineSegment l1 (point(0,0,0), point(1,0,0));
    //halfspace h1 (l1.firstPoint(), vector(-1,0,0));

    //auto result = intersect<point>(l1, h1); 

    //auto resultReverse = intersect<point>(h1, l1); 


    //ASSERT_TRUE(result == resultReverse);
    //ASSERT_TRUE(result == h1.position());
//}

//TEST(lineSegmentIntersection, halfspaceRegularIntersection)
//{
    //typedef Halfspace<point, vector> halfspace;
    //typedef LineSegment<point> lineSegment;

    //lineSegment l1 (point(0,0,0), point(1,0,0));

    //auto testResult = 0.5 * (l1.firstPoint() + l1.secondPoint());

    //halfspace h1 (testResult, vector(-1,0,0));

    //auto result = intersect<point>(l1, h1); 

    //auto resultReverse = intersect<point>(h1, l1); 

    //ASSERT_TRUE(result == resultReverse);
    //ASSERT_TRUE(result == testResult);  
//}

//TEST(polygonIntersection, halfspaceNoIntersection)
//{
    //const double LENGTH = 1.0; 

    //polygon poly = {
        //point(0.,0.,0.), 
        //point(LENGTH,0.,0.),
        //point(LENGTH,LENGTH,0.), 
        //point(0.,LENGTH,0.)
    //};

    //halfspace hspace (point(2, 2, 0), vector(1,1,0)); 

    //typedef PointSequenceIntersection<polygon> polygonIntersection;

    //auto intersection = intersect<polygonIntersection>(poly, h);
    //const auto& result = intersection.polygon();

    //write_vtk_polydata(result, fullTestFileName());  

    //ASSERT_TRUE(result.size() == 0); 

    //auto intersectionReverse = intersect<polygonIntersection>(h, poly);

    //const auto& resultReverse = intersectionReverse.polygon();
    //ASSERT_TRUE(std::equal(result.begin(), result.end(), resultReverse.begin()));
//}

//TEST(polygonIntersection, halfspacePointFullIntersection)
//{
    //auto testLengthScale = [](double LENGTH, std::size_t NEPSILON)
    //{
        //vectorPolygon poly = {
            //point{0.,0.,0.}, 
            //point{LENGTH, 0.,0.},
            //point{LENGTH, LENGTH,0.}, 
            //point{0.,LENGTH,0.}
        //};

        //halfspace hspace (vector(LENGTH, LENGTH, 0.), vector(-1.,-1.,0.)); 

        //auto test = [](const auto& result, const auto& testPolygon)
        //{
            //bool sameSizes = result.polygon_size() == testPolygon.size();
            //ASSERT_TRUE(sameSizes);
                ////<< result.polygon() << std::endl
                ////<< testPolygon << std::endl;

            //bool exactlyEqual = std::equal(result.polygon_begin(), 
                                           //result.polygon_end(), testPolygon.begin());
            //ASSERT_TRUE(exactlyEqual); 
        //};

        //auto result = intersect_tolerance(poly, hspace);
        //auto resultReverse = intersect_tolerance(hspace, poly);

        //test(result, poly); 
        //test(resultReverse, poly); 

        //auto hspacePos = hspace.position(); 
        //auto hspaceOri = hspace.orientation();
        //for (std::size_t testI = 0; testI < 1000; ++testI)
        //{
            //// Perturb halfspace position for NEPSILON * EPSILON.
            //auto hspacePosPert = tolerance_perturbed(hspacePos, NEPSILON * EPSILON); 
            //ASSERT_NE(hspacePos, hspacePosPert);
            //auto hspaceOriPert = tolerance_perturbed(hspaceOri, NEPSILON * EPSILON); 
            //ASSERT_NE(hspaceOri, hspaceOriPert);
            //hspace.setPosition(hspacePosPert); 
            //hspace.setOrientation(hspaceOriPert); 
            //auto result = intersect_tolerance(poly, hspace);
            //test(result, poly);

        //}

        ////auto hspaceOri = hspace.orientation();
        ////for (std::size_t testI = 0; testI < 1000; ++testI)
        ////{
            ////// Perturb halfspce orientation for NEPSILON * EPSILON.
            ////auto hspaceOriPert = tolerance_perturbed(hspaceOri, NEPSILON * EPSILON); 
            ////ASSERT_NE(hspaceOri, hspaceOriPert);
            ////hspace.setOrientation(hspaceOriPert); 
            ////result = intersect_tolerance(poly, hspace);
            ////test(result, poly);
        ////}
    //};

    //testLengthScale(0.125, 1); 
    //testLengthScale(0.5, 1); 
    //testLengthScale(1.0, 1); 
    //testLengthScale(2.0, 1); 
    //testLengthScale(4.0, 1); 
//}

// ******************************************       BEGIN Polygon perturbed

// TODO: Perturbed tests 
//TEST(polygonIntersection, halfspaceRegularIntersectionPerturbed)
//{
    //vectorPolygon poly = {
        //point(0.,0.,0.), 
        //point(1.,0.,0.),
        //point(1.,1.,0.), 
        //point(0.,1.,0.)
    //};

    //auto originalPoly = poly;  

    //halfspace h (point(0.75, 0.75, 0.), vector(-1.,-1.,0.)); 

    //for (unsigned int i = 0; i < MAX_IT; ++i)
    //{
        //// TODO: Perturb
        //perturb(poly, EPSILON); 

        //auto isect = intersect(poly, h);

        //const auto& isectPoly = isect.polygon();

        //ASSERT_TRUE(isectPoly.size() == 5);

        //// Tolearnce-equal check for a vector.
        ////ASSERT_TRUE(equal_under_tolerance(result[2], point(1,0.5,0), EPSILON));
        ////ASSERT_TRUE(equal_under_tolerance(result[3], point(0.5,1,0), EPSILON));

        ////auto intersectionReverse = intersect<polygonIntersection>(h, poly);
        ////const auto& resultReverse = intersectionReverse.polygon();
        ////ASSERT_TRUE(std::equal(result.begin(), result.end(), resultReverse.begin()));

        //poly = originalPoly; 
    //}
//}

//TEST(polygonIntersection, halfspaceEdgeFullIntersectionPerturbed) 
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef Halfspace<point, vector> halfspace;
    //typedef PointSequenceIntersection<polygon> polygonIntersection;

    //polygon inputPolygon = {
        //point(0,0,0), 
        //point(1,0,0),
        //point(1,1,0), 
        //point(0,1,0)
    //};

    //auto originalPolygon = inputPolygon; 
    //const double EPSILON = 1e-15; 
    //halfspace hspace(point(1,1,0), vector(-1,0,0)); 

    //for (unsigned int i = 0; i < MAX_IT; ++i)
    //{
        //perturb(inputPolygon, EPSILON); 
        ////Info << "perturbed : " << inputPolygon << endl; 
        
        //auto intersection = intersect<polygonIntersection>(inputPolygon, hspace, EPSILON); 
        //auto intersectionPolygon = intersection.polygon(); 
        //bool testFailed = intersectionPolygon.size() != 4; 
        //if (testFailed) {
            //write_vtk_polydata(intersectionPolygon, fullTestName(), i); 
            //Info << intersection.polygon();
        //}

        //inputPolygon = originalPolygon; 
        ////Info << "reset : " << inputPolygon << endl;

        //ASSERT_FALSE(testFailed)
            //<< "Polygon size : " << intersectionPolygon.size() << "\n" 
            //<< "Halfspace point : " << hspace.position() << "\n"
            //<< "EPSILON : " << EPSILON << "\n"; 
    //}
//}


//// A unit square is intersected with a plane cutting the square along a diagonal.  
//// Upper left triangle should be the result, consisting of three points and one 
//// half of the original input polygon area. 
//TEST(polygonIntersection, halfspaceDiagonalIntersectionPerturbed)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef Halfspace<point, vector> halfspace;

    //polygon inputPolygon = {
        //point(0,0,0), 
        //point(1,0,0),
        //point(1,1,0), 
        //point(0,1,0)
    //};

    //auto originalPolygon = inputPolygon; 

    //halfspace hspace(point(0.5, 0.5, 0), vector(1,1,0)); 

    //typedef PointSequenceIntersection<polygon> polygonIntersection;

    //for (unsigned int i = 0; i < MAX_IT; ++i)
    //{
        //perturb(inputPolygon, EPSILON); 

        ////Info << "perturbed input polygon : " << inputPolygon << endl;
        //auto intersection = intersect<polygonIntersection>(inputPolygon, hspace);
        //const auto& intersectionPolygon = intersection.polygon();

        ////Info << "intersection polygon : " << intersectionPolygon << endl;

        //if (intersectionPolygon.size() != 3) 
            //write_vtk_polydata(intersectionPolygon, fullTestFileName()); 

        //ASSERT_TRUE(intersectionPolygon.size() == 3)
            //<< "intersection polygon size : " 
            //<< intersectionPolygon.size() << "\n";

        ////auto intersectionReverse = intersect<polygonIntersection>(h, poly);
        ////const auto& resultReverse = intersectionReverse.polygon();
        ////ASSERT_TRUE(std::equal(result.begin(), result.end(), resultReverse.begin()));

        //inputPolygon = originalPolygon; 
    //}
//}


//TEST(polygonIntersection, halfspaceRegularIntersectionPerturbed)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef Halfspace<point, vector> halfspace;

    //polygon poly = {
        //point(0,0,0), 
        //point(1,0,0),
        //point(1,1,0), 
        //point(0,1,0)
    //};

    //halfspace h (point(0.75, 0.75, 0), vector(-1,-1,0)); 

    //typedef PointSequenceIntersection<polygon> polygonIntersection;

    //auto intersection = intersect<polygonIntersection>(poly, h);

    //const auto& result = intersection.polygon();

    //write_vtk_polydata(result, fullTestFileName()); 

    //ASSERT_TRUE(result.size() == 5);
    //ASSERT_TRUE(result[2] == point(1,0.5,0));
    //ASSERT_TRUE(result[3] == point(0.5,1,0));

    //auto intersectionReverse = intersect<polygonIntersection>(h, poly);
    //const auto& resultReverse = intersectionReverse.polygon();
    //ASSERT_TRUE(std::equal(result.begin(), result.end(), resultReverse.begin()));
//}

// ******************************************       END Polygon perturbed 
// ******************************************       BEGIN AABB 

//TEST(aabboxIntersection, AABBtrueIntersection)
//{
    //typedef AABBox<point> aabbox; 

    //aabbox b1(point(0,0,0), point(1,1,1)); 
    //aabbox b2(point(0.9,0.9,0.9), point(2,2,2)); 

    //ASSERT_TRUE(aabb_do_intersect(b1,b2));
//}

//TEST(aabboxIntersection, AABBpointIntersection)
//{
    //typedef AABBox<point> aabbox; 

    //aabbox b1(point(0,0,0), point(1,1,1)); 
    //aabbox b2(point(1,1,1), point(2,2,2)); 

    //ASSERT_FALSE(aabb_do_intersect(b1,b2))
       //<< "AABBs that share points should not intersect. \n";
//}

//TEST(aabboxIntersection, AABBedgeIntersection)
//{
    //typedef AABBox<point> aabbox; 

    //aabbox b1(point(0,0,0), point(1,1,1)); 
    //aabbox b2(point(1,1,0), point(2,2,1)); 

    //ASSERT_FALSE(aabb_do_intersect(b1,b2))
       //<< "AABBs that share edges should not intersect. \n"; 
//}

//TEST(aabboxIntersection, AABBfaceIntersection)
//{
    //typedef AABBox<point> aabbox; 

    //aabbox b1(point(0,0,0), point(1,1,1)); 
    //aabbox b2(point(1,0,0), point(2,1,1)); 

    //ASSERT_FALSE(aabb_do_intersect(b1,b2))
       //<< "AABBs that share faces should not intersect. \n"; 
//}

//TEST(aabboxIntersection, AABBnoIntersection)
//{
    //typedef AABBox<point> aabbox; 

    //aabbox b1(point(-10,-10,-10), point(2,3,4)); 
    //aabbox b2(point(10,15,20), point(17,24,35)); 

    //ASSERT_FALSE(aabb_do_intersect(b1,b2));
//}

// ******************************************       END AABB 
// ******************************************       BEGIN Polyhedron AABB 

//TEST(polyhedronIntersection, AABBpointIntersection) 
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;

    //point pMin = point(0,0,0); 
    //point pMax = point(1,1,1);

    //auto P = build<polyhedron>(aabbox(pMin, pMax));

    //write_vtk_polydata(P, fullTestName() + "P.vtk");

    //point qMin = point(1,1,1); 
    //point qMax = point(2,2,2);

    //auto Q = build<polyhedron>(aabbox(qMin, qMax));

    //write_vtk_polydata(Q, fullTestName() + "Q.vtk");

    //auto Pbox = build<aabbox>(P); 

    //ASSERT_TRUE(Pbox.minPoint() == pMin);  
    //ASSERT_TRUE(Pbox.maxPoint() == pMax); 

    //auto Qbox = build<aabbox>(Q);

    //ASSERT_TRUE(Qbox.minPoint() == qMin);  
    //ASSERT_TRUE(Qbox.maxPoint() == qMax); 

    //ASSERT_FALSE(aabb_do_intersect(P, Q))
       //<< "Polyhedra that share points should not intersect. \n"; 
//}

//TEST(polyhedronIntersection, AABBedgeIntersection) 
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;

    //point pMin = point(0,0,0); 
    //point pMax = point(1,1,1);

    //auto P = build<polyhedron>(aabbox(pMin, pMax));

    //write_vtk_polydata(P, fullTestName() + "P.vtk");

    //point qMin = point(1,1,0); 
    //point qMax = point(2,2,1);

    //auto Q = build<polyhedron>(aabbox(qMin, qMax));

    //write_vtk_polydata(Q, fullTestName() + "Q.vtk");

    //auto Pbox = build<aabbox>(P); 

    //ASSERT_TRUE(Pbox.minPoint() == pMin);  
    //ASSERT_TRUE(Pbox.maxPoint() == pMax); 

    //auto Qbox = build<aabbox>(Q);

    //ASSERT_TRUE(Qbox.minPoint() == qMin);  
    //ASSERT_TRUE(Qbox.maxPoint() == qMax); 

    //ASSERT_FALSE(aabb_do_intersect(P, Q))
       //<< "Polyhedra that share edges should not intersect. \n"; 
//}

//TEST(polyhedronIntersection, AABBfaceIntersection) 
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;

    //point pMin = point(0,0,0); 
    //point pMax = point(1,1,1);

    //auto P = build<polyhedron>(aabbox(pMin, pMax));

    //write_vtk_polydata(P, fullTestName() + "P.vtk");

    //point qMin = point(1,0,0); 
    //point qMax = point(2,1,1);

    //auto Q = build<polyhedron>(aabbox(qMin, qMax));

    //write_vtk_polydata(Q, fullTestName() + "Q.vtk");

    //auto Pbox = build<aabbox>(P); 

    //ASSERT_TRUE(Pbox.minPoint() == pMin);  
    //ASSERT_TRUE(Pbox.maxPoint() == pMax); 

    //auto Qbox = build<aabbox>(Q);

    //ASSERT_TRUE(Qbox.minPoint() == qMin);  
    //ASSERT_TRUE(Qbox.maxPoint() == qMax); 

    //ASSERT_FALSE(aabb_do_intersect(P, Q))
       //<< "Polyhedra that share faces should not intersect. \n"; 
//}


//TEST(polyhedronIntersection, AABBtrueIntersection) 
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;

    //auto P = build<polyhedron>(aabbox(point(0,0,0), point(1,1,1)));
    
    //auto Q = build<polyhedron>(aabbox(point(0.5,0.5,0.5), point(1.5,1.5,1.5)));  

    //ASSERT_TRUE(aabb_do_intersect(P,Q)); 
//}


//// Polyhedron/halfspace intersection calculations. 

//TEST(polyhedronHalfspaceIntersection, point) 
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef Halfspace<point, vector> halfspace;

    //auto P = build<polyhedron>(aabbox(point(0,0,0), point(1,1,1)));

    //halfspace h(point(1,1,1), vector(1,1,1)); 
    
    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection; 

    //auto intersection = intersect<polyhedronIntersection>(P, h); 

    //const auto& poly = intersection.polyhedron(); 

    //write_vtk_polydata(poly, fullTestFileName());
//}

// ******************************************       END Polyhedron AABB 
//TEST(polyhedronCube, edge) 
//{
    //auto cube = make<vectorPolyhedron>(point(0.,0.,0.), point(1.,1.,1.));

    //halfspace hspace(point(1.,1.,1.), vector(1.,1.,0.)); 

    //auto intersection = intersect_tolerance(cube, hspace); 
    //const auto& polyhedron = intersection.polyhedron(); 

    //ASSERT_TRUE(poly.size() == 0); 
//}

//TEST(polyhedronHalfspaceIntersection, face) 
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef Halfspace<point, vector> halfspace;

    //auto P = build<polyhedron>(aabbox(point(0,0,0), point(1,1,1)));

    //halfspace h(point(1,1,1), vector(1,0,0)); 

    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection; 
    
    //auto intersection = intersect<polyhedronIntersection>(P, h); 

    //auto poly = intersection.polyhedron(); 

    //write_vtk_polydata(poly, fullTestFileName()); 

    //ASSERT_TRUE(poly.size() == 0); 
//}

//TEST(polyhedronHalfspaceIntersection, faceIntersection) 
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef Halfspace<point, vector> halfspace;

    //auto P = build<polyhedron>(aabbox(point(0,0,0), point(1,1,1)));

    //halfspace h(point(1,1,1), vector(-1,0,0)); 

    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection; 
    
    //auto intersection = intersect<polyhedronIntersection>(P, h); 

    //auto poly = intersection.polyhedron(); 

    //write_vtk_polydata(poly, fullTestFileName()); 

    //ASSERT_TRUE(poly.size() == 6)
       //<< "poly.size() = " << poly.size() << "\n"; 

    //for (const auto & polygon : poly)
    //{
        //ASSERT_TRUE(polygon.size() == 4)
           //<< "polygon.size = " << polygon.size() << "\n"; 
    //}

    //ASSERT_LE(mag(volume(poly) - 1), EPSILON)
       //<< "polyhedron volume = " << volume(poly) << "\n";
//}


//// A unit box polyhedron intersected with a plane positioned at (1,1,1)
//// oriented at (-1, 0, 0). A result should be the the polyhedron itself. 
//// The polyhedron points are individually randomly perturbed.
//TEST(polyhedronHalfspaceIntersection, faceIntersectionPerturbed) 
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef Halfspace<point, vector> halfspace;
    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection; 

    //auto inputPolyhedron = build<polyhedron>(aabbox(point(0,0,0), point(1,1,1)));

    //write_vtk_polydata(inputPolyhedron, fullTestName() + ".inputPolyhedron.vtk"); 

    //auto originalPolyhedron = inputPolyhedron; 
    //halfspace hspace(point(1,1,1), vector(-1,0,0)); 

    //write_vtk_polydata(hspace, fullTestName() + ".halfspace.vtk"); 

    //for (unsigned int i = 0; i < MAX_IT; ++i)
    //{
        //perturb(inputPolyhedron, EPSILON); 

        //auto intersection = intersect<polyhedronIntersection>(inputPolyhedron, hspace, EPSILON); 

        //auto intersectionPolyhedron = intersection.polyhedron(); 

        //if 
        //( 
            //(intersectionPolyhedron.size() != 6) || 
            //(mag(volume(intersectionPolyhedron) - 1) >  (6*EPSILON))
        //)
        //{
            //write_vtk_polydata(intersectionPolyhedron, fullTestName() + ".vtk"); 
        //}


        //ASSERT_TRUE(intersectionPolyhedron.size() == 6) 
            //<< "Polyhedron size : " << intersectionPolyhedron.size() << "\n";

        //for (const auto & polygon : intersectionPolyhedron)
        //{
            //if (polygon.size() != 4)
            //{
                //write_vtk_polydata(polygon, fullTestName() + ".polygon.vtk"); 
                //write_vtk_polydata(intersectionPolyhedron, fullTestName() + ".vtk"); 
            //}
            //ASSERT_TRUE(polygon.size() == 4)
               //<< "Polygon size : " << polygon.size() << "\n"; 
        //}
        
        //// Volume inflated on all 6 sides with epsilon, where side edges have
        //// unit lengths, result in a 6 epsilon maximal volume change introduced
        //// by the perturbation. TM.  
        //ASSERT_LE(mag(volume(intersectionPolyhedron) - 1),  (6 * EPSILON))
            //<< "Volume : " << volume(intersectionPolyhedron) << "\n"
            //<< "Polyhedron size: " << intersectionPolyhedron.size() << "\n";

        //inputPolyhedron = originalPolyhedron; 
    //}
//}

//TEST(polyhedronIntersection, cubeSpatialDiagonalIntersection) 
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef Halfspace<point, vector> halfspace;

    //auto P= build<polyhedron>(aabbox(point(0,0,0), point(1,1,1)));

    //halfspace h(point(0.5,0.5,0.5), vector(1,1,1)); 

    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection; 
    
    //auto intersection = intersect<polyhedronIntersection>(P, h); 

    //auto poly = intersection.polyhedron(); 

    //write_vtk_polydata(poly, fullTestFileName()); 

    //ASSERT_TRUE(poly.size() == 7); 

    //auto capPolygon = poly.back();  

    //ASSERT_TRUE(capPolygon.size() == 6); 
//}

//TEST(polyhedronIntersection, cubeHalvedIntersection) 
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef Halfspace<point, vector> halfspace;

    //auto P = build<polyhedron>(aabbox(point(0,0,0), point(1,1,1)));

    //halfspace h(point(0.5,0.5,0.5), vector(0,0,1)); 
    
    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection; 

    //auto intersection = intersect<polyhedronIntersection>(P, h); 

    //const auto& poly = intersection.polyhedron(); 

    //write_vtk_polydata(poly, fullTestFileName()); 

    //ASSERT_TRUE(poly.size() == 6); 

    //for (const auto& polygon : poly)
    //{
        //ASSERT_TRUE(polygon.size() == 4);
    //}
//}

//TEST(polyhedronIntersection, cubeHalvedDiagonally) 
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef Halfspace<point, vector> halfspace;

    //auto P = build<polyhedron>(aabbox(point(0,0,0), point(1,1,1)));

    //halfspace h(point(0.5,0.5,0.5), vector(0,-1,1)); 
    
    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection; 

    //auto intersection = intersect<polyhedronIntersection>(P, h); 

    //const auto& poly = intersection.polyhedron(); 

    //write_vtk_polydata(poly, fullTestFileName()); 

    //ASSERT_TRUE(poly.size() == 5); 
//}

//TEST(isInsideTest, triangleInsideHalfspace)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef Halfspace<point, vector> halfspace;
    
    //halfspace h(point(2,2,2), vector(-1,-1,-1)); 

    //auto t = build<triangle>(point(0,0,0), point(1,0,0), point(0,1,0));

    //// Fully inside.
    //ASSERT_TRUE(is_inside(h, t)); 

    //// Touching halfspace with a point is considered intersection.
    //h = halfspace(point(1,0,0), vector(-1,0,0)); 
    //ASSERT_FALSE(is_inside(h, t)); 

    //// Touching halfspace with an edge is considered as intersection.
    //h = halfspace(point(0,0,0), vector(0,1,0)); 
    //ASSERT_FALSE(is_inside(h, t)); 
//}

//TEST(isInsideTest, tetrahedronInsideHalfspace)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef Halfspace<point, vector> halfspace;
    
    //halfspace h(point(2,2,2), vector(-1,-1,-1)); 

    //auto t = build<tetrahedron>(point(0,0,0), point(1,0,0), point(0,1,0), point(0,0,1));

    //// Fully inside. 
    //ASSERT_TRUE(is_inside(h, t)); 

    //// Touching halfspace with a point is an intersection. TM.
    //h = halfspace(point(1,0,0),vector(-1,0,0));
    //ASSERT_FALSE(is_inside(h, t)); 
    
    //// Touching halfspace with an edge is an intersection. TM.
    //h = halfspace(point(0,0,0),vector(1,1,0));
    //ASSERT_FALSE(is_inside(h, t)); 

    //// Touching halfspace with a face is an intersection. TM.
    //h = halfspace(point(0,0,0),vector(0,0,1));
    //ASSERT_FALSE(is_inside(h, t)); 
//}

/// @brief For each point of each square polygon, check if the point is inside the other polygon. 
/// Different configurations are tested for overlapping squares and the number of inside points
/// is checked against a known expected value. 
/// @param polygonIsInsideTest
/// @param overlappingSquares
//TEST(polygonIsInsideTest, overlappingSquares)
//{
    //const pointVector basePolygon = {point(0,0,0), point(1,0,0),point(1,1,0),point(0,1,0)}; 
    ///// @brief Return points of basePolygon and capProjection that reside in both of them. 
    /////
    ///// @param basePolygon
    ///// @param capPolygon
    ///// @param tolerance
    /////
    ///// @return Inside points. 
    //auto calcInsidePoints = [] (
        //const auto& basePolygon, 
        //const auto& capPolygon,
        //double tolerance = EPSILON 
    //) -> auto
    //{
        //std::vector<point> insidePoints; 

        //auto baseIt = basePolygon.begin(); 
        //auto capIt = capPolygon.begin(); 

        //auto baseNormal = normal_area_vec(basePolygon); 
        //baseNormal /= mag(baseNormal); 
        
        //for(; baseIt != basePolygon.end(); ++baseIt, ++capIt)
        //{
            //if (is_inside(capPolygon, *baseIt, baseNormal, tolerance))
                //insidePoints.push_back(*baseIt); 

            //if (is_inside(basePolygon, *capIt, baseNormal, tolerance))
                //insidePoints.push_back(*capIt); 
        //}

        //return insidePoints; 
    //};
    
    //auto insidePoints = calcInsidePoints(basePolygon, basePolygon);
    //ASSERT_TRUE(insidePoints.size() == basePolygon.size() * 2)
        //<< "Number of inside points = " << insidePoints.size() << "\n";

    //auto basePolygonPerturbed(basePolygon); 
    //perturb(basePolygonPerturbed, EPSILON); 
    //insidePoints = calcInsidePoints(basePolygonPerturbed, basePolygon, EPSILON);
    //ASSERT_TRUE(insidePoints.size() == basePolygon.size() * 2)
        //<< "Number of inside points = " << insidePoints.size() << "\n";

    //auto quarterOverlap = displace(basePolygon, vector(0.5, 0.5, 0));  
    //insidePoints = calcInsidePoints(basePolygon, quarterOverlap);

    //ASSERT_TRUE(insidePoints.size() == 2)
        //<< "Number of inside points = " << insidePoints.size() << "\n";

    //auto edgeOverlap = displace(basePolygon, vector(1, 0, 0));  
    //insidePoints = calcInsidePoints(basePolygon, edgeOverlap);

    //ASSERT_TRUE(insidePoints.size() == 4)
        //<< "Number of inside points = " << insidePoints.size() << "\n";

    //auto pointOverlap = displace(basePolygon, vector(1, 1, 0));  
    //insidePoints = calcInsidePoints(basePolygon, pointOverlap);

    //ASSERT_TRUE(insidePoints.size() == 2)
        //<< "Number of inside points = " << insidePoints.size() << "\n";
//}


//TEST(triangleHalfspaceIntersection, triangleInsideHalfspace)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef Halfspace<point, vector> halfspace;
    //typedef TriangleIntersection<triangle, pointVector> triangleHalfspaceIntersection;

    //auto inputTriangle = build<triangle>(point(0,0,0), point(1,0,0), point(0,1,0));

    //halfspace inputHalfspace(point(2,2,2), vector(-1,-1,-1)); 

    //auto intersectionResult = intersect<triangleHalfspaceIntersection>(inputTriangle, inputHalfspace); 

    //const auto& resultPolygon = intersectionResult.polygon(); 

    //write_vtk_polydata(resultPolygon, fullTestFileName()); 

    //ASSERT_TRUE(resultPolygon.size() == 3); 
    //ASSERT_TRUE(std::equal(resultPolygon.begin(), resultPolygon.end(), inputTriangle.begin()));
//} 

//TEST(triangleHalfspaceIntersection, pointIntersection)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef Halfspace<point, vector> halfspace;
    //typedef TriangleIntersection<triangle, pointVector> triangleHalfspaceIntersection;

    //auto inputTriangle = build<triangle>(point(0,0,0), point(1,0,0), point(0,1,0));

    //halfspace inputHalfspace(point(1,0,0), vector(-1,0,0)); 

    //auto intersectionResult = intersect<triangleHalfspaceIntersection>(inputTriangle, inputHalfspace); 

    //const auto& resultPolygon = intersectionResult.polygon(); 

    //write_vtk_polydata(resultPolygon, fullTestFileName()); 

    //ASSERT_TRUE(resultPolygon.size() == 3); 
    //ASSERT_TRUE(std::equal(resultPolygon.begin(), resultPolygon.end(), inputTriangle.begin()));

    //const auto intersectionPoints = intersectionResult.intersectionPoints(); 

    //ASSERT_TRUE(intersectionPoints.size() == 1);
    //ASSERT_TRUE(intersectionPoints[0] == point(1,0,0));
//} 


//TEST(triangleHalfspaceIntersection, edgeIntersection)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef Halfspace<point, vector> halfspace;
    //typedef TriangleIntersection<triangle, pointVector> triangleHalfspaceIntersection;

    //auto inputTriangle = build<triangle>(point(0,0,0), point(1,0,0), point(0,1,0));

    //halfspace inputHalfspace(point(0,1,0), vector(1,0,0)); 

    //auto intersectionResult = intersect<triangleHalfspaceIntersection>(inputTriangle, inputHalfspace); 

    //const auto& resultPolygon = intersectionResult.polygon(); 

    //write_vtk_polydata(resultPolygon, fullTestFileName()); 

    //ASSERT_TRUE(resultPolygon.size() == 3); 
    //ASSERT_TRUE(std::equal(resultPolygon.begin(), resultPolygon.end(), inputTriangle.begin()));

    //const auto intersectionPoints = intersectionResult.intersectionPoints(); 

    //ASSERT_TRUE(intersectionPoints.size() == 2);
    //ASSERT_TRUE(intersectionPoints[0] == point(0,0,0));
    //ASSERT_TRUE(intersectionPoints[1] == point(0,1,0));
//} 


//TEST(triangleHalfspaceIntersection, diagonalIntersection)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef Halfspace<point, vector> halfspace;
    //typedef TriangleIntersection<triangle, pointVector> triangleHalfspaceIntersection;

    //auto inputTriangle = build<triangle>(point(0,0,0), point(1,0,0), point(0,1,0));

    //write_vtk_polydata(inputTriangle, fullTestName() + ".inputTriangle.vtk"); 

    //halfspace inputHalfspace(point(0,0,0), vector(-1,1,0)); 

    //auto intersectionResult = intersect<triangleHalfspaceIntersection>(inputTriangle, inputHalfspace); 

    //const auto& resultPolygon = intersectionResult.polygon(); 

    //write_vtk_polydata(resultPolygon, fullTestName() + "resultPolygon.vtk"); 

    //ASSERT_TRUE(resultPolygon.size() == 3); 

    //const auto intersectionPoints = intersectionResult.intersectionPoints(); 

    //auto inputTriangleNormal = normal_area_vec(inputTriangle); 
    //inputTriangleNormal /= mag(inputTriangleNormal); 

    //auto resultPolygonNormal = normal_area_vec(resultPolygon); 
    //resultPolygonNormal /= mag(resultPolygonNormal); 

    //ASSERT_LE(mag(inputTriangleNormal - resultPolygonNormal), Traits::tolerance<scalar>::value());

    //ASSERT_TRUE(intersectionPoints.size() == 2);
    //ASSERT_TRUE(intersectionPoints[0] == point(0,0,0))
        //<< "intersectionPoints[0]" << intersectionPoints[0] << "\n";
    //ASSERT_TRUE(intersectionPoints[1] == point(0.5,0.5,0))
        //<< "intersectionPoints[1]" << intersectionPoints[1] << "\n";
//} 

//TEST(triangleHalfspaceIntersection, regularIntersection)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef Halfspace<point, vector> halfspace;
    //typedef TriangleIntersection<triangle, pointVector> triangleHalfspaceIntersection;

    //auto inputTriangle = build<triangle>(point(0,0,0), point(1,0,0), point(0,1,0));

    //halfspace inputHalfspace(point(0,0.5,0), vector(0,-1,0)); 

    //auto intersectionResult = intersect<triangleHalfspaceIntersection>(inputTriangle, inputHalfspace); 

    //const auto& resultPolygon = intersectionResult.polygon(); 

    //write_vtk_polydata(resultPolygon, fullTestFileName()); 

    //ASSERT_TRUE(resultPolygon.size() == 4); 
    //ASSERT_TRUE(resultPolygon[0] == point(0,0.5,0));
    //ASSERT_TRUE(resultPolygon[1] == point(0,0,0));
    //ASSERT_TRUE(resultPolygon[2] == point(1,0,0));
    //ASSERT_TRUE(resultPolygon[3] == point(0.5,0.5,0));

    //const auto intersectionPoints = intersectionResult.intersectionPoints(); 

    //ASSERT_TRUE(intersectionPoints.size() == 2);
    //ASSERT_TRUE(intersectionPoints[0] == point(0,0.5,0));
    //ASSERT_TRUE(intersectionPoints[1] == point(0.5,0.5,0));
//} 

//TEST(tetrahedronHalfspaceIntersection, tetrahedronInsideHalfspace)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef Halfspace<point, vector> halfspace;
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection;
    
    //halfspace hspace (point(2,2,2), vector(-1,-1,-1)); 

    //auto tet = build<tetrahedron>(point(0,0,0), point(1,0,0), point(0,1,0), point(0,0,1));
    
    //auto result = intersect<polyhedronIntersection>(tet, hspace);  

    //const auto& resultPolyhedron = result.polyhedron(); 

    //write_vtk_polydata(resultPolyhedron, fullTestFileName()); 

    //ASSERT_TRUE(resultPolyhedron.size() == 4);

    //// Test to see that the result polyhedron is exactly the same as the
    //// input tetrahedron.
    //for (char I = 0; I < 4; ++I)
    //{
        //ASSERT_TRUE(resultPolyhedron[I][0] == tet[I][0]);
        //ASSERT_TRUE(resultPolyhedron[I][1] == tet[I][1]);
        //ASSERT_TRUE(resultPolyhedron[I][2] == tet[I][2]);
    //}
//}

//TEST(tetrahedronHalfspaceIntersection, pointIntersection)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef Halfspace<point, vector> halfspace;
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection;
    
    //halfspace hspace (point(1,0,0), vector(-1,1,0)); 

    //auto tet = build<tetrahedron>(point(0,0,0), point(1,0,0), point(0,1,0), point(0,0,1));

    //auto result = intersect<polyhedronIntersection>(tet, hspace);  

    //const auto& resultPolyhedron = result.polyhedron(); 

    //write_vtk_polydata(resultPolyhedron, fullTestFileName()); 

    //ASSERT_TRUE(resultPolyhedron.size() == 4);

    //// Test to see that the result polyhedron is exactly the same as the
    //// input tetrahedron.
    //for (char I = 0; I < 4; ++I)
    //{
        //ASSERT_TRUE(resultPolyhedron[I][0] == tet[I][0]);
        //ASSERT_TRUE(resultPolyhedron[I][1] == tet[I][1]);
        //ASSERT_TRUE(resultPolyhedron[I][2] == tet[I][2]);
    //}
//}

//TEST(tetrahedronHalfspaceIntersection, pointRegularIntersection)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef Halfspace<point, vector> halfspace;
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection;
    
    //halfspace hspace (point(0,0,0), vector(0,1,0)); 

    //auto tet = build<tetrahedron>(point(0,0,0), point(0.5,-0.5,0), point(0.5,0.5,0), point(0.5,-0.25,1));

    //write_vtk_polydata(tet, fullTestName() + ".tet.vtk");

    //auto result = intersect<polyhedronIntersection>(tet, hspace);  

    //const auto& resultPolyhedron = result.polyhedron(); 

    //write_vtk_polydata(resultPolyhedron, fullTestFileName()); 

    //ASSERT_TRUE(resultPolyhedron.size() == 4);

    //halfspace hspacePoly (point(0.5,-0.25,1), vector(0,1,0)); 

    //auto resultPoly = intersect<polyhedronIntersection>(tet, hspacePoly);  

    //const auto& resultPolyPolyhedron = resultPoly.polyhedron(); 

    //write_vtk_polydata(resultPolyPolyhedron, fullTestName() + ".polyhedronResult.vtk"); 

    //ASSERT_TRUE(resultPolyPolyhedron.size() == 5);
//}

//TEST(tetrahedronHalfspaceIntersection, edgeIntersection)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef Halfspace<point, vector> halfspace;
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection;
    
    //halfspace hspace (point(0,0,0), vector(0,1,0.5)); 

    //auto tet = build<tetrahedron>(point(0,0,0), point(1,0,0), point(0,1,0), point(0,0,1));

    //auto result = intersect<polyhedronIntersection>(tet, hspace);  

    //const auto& resultPolyhedron = result.polyhedron(); 

    //write_vtk_polydata(resultPolyhedron, fullTestFileName()); 

    //ASSERT_TRUE(resultPolyhedron.size() == 4);

    //// Test to see that the result polyhedron is exactly the same as the
    //// input tetrahedron.
    //for (char I = 0; I < 4; ++I)
    //{
        //ASSERT_TRUE(resultPolyhedron[I][0] == tet[I][0]);
        //ASSERT_TRUE(resultPolyhedron[I][1] == tet[I][1]);
        //ASSERT_TRUE(resultPolyhedron[I][2] == tet[I][2]);
    //}
//}

//TEST(tetrahedronHalfspaceIntersection, faceIntersection)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef Halfspace<point, vector> halfspace;
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection;
    
    //halfspace hspace (point(0,0,0), vector(0,0,1)); 

    //write_vtk_polydata(hspace, fullTestName() + ".halfspace.vtk"); 

    //auto tet = build<tetrahedron>(point(0,0,0), point(1,0,0), point(0,1,0), point(0,0,1));

    //write_vtk_polydata(tet, fullTestName() + ".tetrahedron.vtk"); 

    //auto result = intersect<polyhedronIntersection>(tet, hspace);  

    //const auto& resultPolyhedron = result.polyhedron(); 

    //write_vtk_polydata(resultPolyhedron, fullTestName() + ".polyhedron.vtk"); 

    //ASSERT_TRUE(resultPolyhedron.size() == 4)
        //<< "result polyhedron size: " << resultPolyhedron.size() << "\n";

    //// Test triangle sizes. 
    //for (char I = 0; I < 4; ++I)
    //{
        //ASSERT_TRUE(resultPolyhedron[I].size() == 3); 
    //}

    //// Test volume.
    //ASSERT_TRUE(mag(volume(tet) - volume(resultPolyhedron)) < Traits::tolerance<scalar>::value());
//}


//TEST(tetrahedronHalfspaceIntersection, tetrahedronFaceIntersectionPerturbed)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef Halfspace<point, vector> halfspace;
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection;
    
    //auto tet = build<tetrahedron>(point(0,0,0), point(1,0,0), point(0,1,0), point(0,0,1));
    //auto originalTet = tet; 

    //point hspacePoint(0,0,0); 
    //auto originalPoint = hspacePoint; 

    //for (unsigned int i = 0; i < MAX_IT; ++i)
    //{
        //perturb(hspacePoint, EPSILON); 

        //halfspace hspace (hspacePoint, vector(1,0,0)); 

        //perturb(tet, EPSILON);

        //auto result = intersect<polyhedronIntersection>(tet, hspace, EPSILON);  

        //const auto& resultPolyhedron = result.polyhedron(); 

        //write_vtk_polydata(resultPolyhedron, fullTestFileName()); 

        //ASSERT_TRUE(resultPolyhedron.size() == 4)
            //<< "polyhedron size : " << resultPolyhedron.size() << "\n"
            //<< "hspace point: " 
            //<< hspacePoint[0] << " " 
            //<< hspacePoint[1] << " " 
            //<< hspacePoint[2] << " " 
            //<< "\n";

        //for (const auto& polygon : resultPolyhedron)
        //{
            //ASSERT_TRUE(polygon.size() == 3);
        //}

        //// Test volume.
        //auto volumeError = volume(resultPolyhedron) - volume(tet); 

        //ASSERT_LE(mag(volumeError), 6*EPSILON)
           //<< "relative volume error: " << volumeError << "\n"
           //<< "hspace point: " 
           //<< hspacePoint[0] << " " 
           //<< hspacePoint[1] << " " 
           //<< hspacePoint[2] << " " 
           //<< "\n";

        //hspacePoint = originalPoint;
        //tet = originalTet; 
    //}
//}

//TEST(tetrahedronHalfspaceIntersection, incrementalIntersection)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef Halfspace<point, vector> halfspace;
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection;
    
    //halfspace hspace (point(0,0,0), vector(1,0,0)); 

    //double a = 1;
    //double b = 1;
    //double h = 1;

    //auto tet = build<tetrahedron>(point(0,0,0), point(a,0,0), point(0,b,0), point(0,0,h));

    //auto tetVolume = volume(tet); 

    //double exactVolume = 1. / 6. * (a * b * h ); 

    //ASSERT_TRUE(mag(volume(tet) - exactVolume) < Traits::tolerance<scalar>::value());

    //const unsigned int intersectionCount = 1000; 

    //std::ofstream dataFile;
    //dataFile.open(fullTestName() + ".dat"); 

    //for (unsigned int I = 0; I <= intersectionCount; ++I)
    //{
        //double hTrunc = I * 1. / intersectionCount;  

        //halfspace hspace (point(0,0,hTrunc), vector(0,0, -1)); 

        //auto result = intersect<polyhedronIntersection>(tet, hspace);  

        //const auto& resultPolyhedron = result.polyhedron(); 

        //write_vtk_polydata(
            //resultPolyhedron, 
            //fullTestName(), 
            //I
        //);

        //// Test the volume calculation.
        //auto resultVolume = volume(resultPolyhedron); 

        //ASSERT_TRUE(resultVolume >= 0); 

        //auto exactResultVolume = exactVolume - (pow(h - hTrunc, 3) * a * b)  / (6 * h * h);

        //auto volumeError = mag(resultVolume - exactResultVolume); 

        //// Output the volume values into a file.
        //dataFile << h << " " << resultVolume / tetVolume << "\n"; 
        
        //ASSERT_TRUE(volumeError < Traits::tolerance<scalar>::value());
    //}
//}

//TEST(tetPolyStarHalfspaceIntersection, fillLevelFunction)
//{
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef Halfspace<point, vector> halfspace;
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef PolygonSequenceIntersection<polyhedron> polyhedronIntersection;
    //typedef TriangulationIntersection<tetrahedronVector, polyhedron> triangulationIntersection;
    
    //halfspace hspace (point(0,0,0), vector(1,0,0)); 

    //double a = 1;
    //double b = 1;
    //double h = 1;

    //auto tet = build<tetrahedron>(point(0,0,0), point(a,0,0), point(0,b,0), point(0,0,h));
    //auto poly = build<polyhedron>(aabbox(point(0,0,0),point(a,b,h))); 
    //auto starPoly = build_box_star<polyhedron>(point(1.0/3.0,1.0/3.0,1.0/3.0),point(2.0/3.0,2.0/3.0,2.0/3.0),1.0/3.0);
    //auto star = oriented_triangulate<tetrahedronVector>(starPoly); 

    //auto tetVolume = volume(tet); 
    //auto polyVolume = volume(poly); 
    //auto starVolume = volume(star); 

    //const unsigned int intersectionCount = 1000; 

    //std::ofstream dataFile;
    //dataFile.open(fullTestName() + ".dat"); 

    //for (unsigned int I = 0; I <= intersectionCount; ++I)
    //{
        //double hTrunc = I * 2.0 / intersectionCount;  

        //halfspace hspace (point(hTrunc,0,0), vector(1,1,0)); 

        //auto tetResult = intersect<polyhedronIntersection>(tet, hspace);  
        //auto polyResult = intersect<polyhedronIntersection>(poly, hspace);  
        //auto starResult = intersect<triangulationIntersection>(star, hspace);  

        //const auto& tetResultPoly = tetResult.polyhedron(); 
        //const auto& polyResultPoly = polyResult.polyhedron(); 

        //write_vtk_polydata(tetResultPoly, fullTestName() + "-tet", I);
        //write_vtk_polydata(polyResultPoly, fullTestName() + "-poly", I);
        //write_vtk_polydata(starResult, fullTestName() + "-star", I);

        //auto tetResultVol = volume(tetResultPoly); 
        //auto polyResultVol = volume(polyResultPoly); 
        //auto starResultVol = volume(starResult); 

        //// Output the volume values into a file.
        //dataFile << hTrunc << " " << tetResultVol / tetVolume 
            //<< " " << polyResultVol / polyVolume 
            //<< " " << starResultVol / starVolume  << "\n"; 
    //}
//}

// FIXME: enable once the triangulation intersection tests are working again. TM.
//TEST(triangulationHalfspaceIntersection, triangulationInsideHalfspace)
//{
    //typedef AABBox<point> aabb_box; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef Halfspace<point, vector> halfspace;
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef TetrahedronSequence<tetrahedron, std::vector> triangulation;
    //typedef TriangulationIntersection<triangulation, polyhedron> triangulationIntersection;

    //polyhedron poly = build<polyhedron>(aabb_box(point(0,0,0), point(1,1,1)));

    //auto tri = barycentric_triangulate<triangulation>(poly); 

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

    //halfspace hspace (point(2,2,2), vector(-1,-1,-1)); 

    //auto intersection = intersect<triangulationIntersection>(tri, hspace); 

    //write_vtk_polydata(intersection, fullTestFileName()); 

    //ASSERT_TRUE(intersection.polyhedra_size() == 0);
    //ASSERT_TRUE(intersection.tetrahedra_size() == tri.size());
    //ASSERT_TRUE(mag(volume(intersection) - 1) <= Traits::tolerance<scalar>::value());  
//}

//TEST(triangulationHalfspaceIntersection, faceIntersection)
//{
    //typedef AABBox<point> aabb_box; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef Halfspace<point, vector> halfspace;
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef TetrahedronSequence<tetrahedron, std::vector> triangulation;
    //typedef TriangulationIntersection<triangulation, polyhedron> triangulationIntersection;

    //polyhedron poly = build<polyhedron>(aabb_box(point(0,0,0), point(1,1,1)));

    //write_vtk_polydata(poly, fullTestName() + ".polyhedron.vtk"); 

    //auto tri = barycentric_triangulate<triangulation>(poly); 

    //write_vtk_polydata(tri, fullTestName() + ".triangulation.vtk"); 

    //auto intersectionError = [&](const halfspace& hspace)
    //{
        //auto intersection = intersect<triangulationIntersection>(tri, hspace, EPSILON); 
        //return mag(volume(intersection) - 1); 
    //};

    //halfspace h1 (point(0,0,0), vector(0,0,1)); 

    //double error = intersectionError(h1); 
    //ASSERT_LE(error, Traits::tolerance<scalar>::value())
       //<< "volume error = " << error << "\n"; 

    //h1.setDirection(vector(0,1,0));
    //error = intersectionError(h1); 
    //ASSERT_LE(error, Traits::tolerance<scalar>::value())
       //<< "volume error = " << error << "\n"; 

    //h1.setDirection(vector(1,0,0));
    //error = intersectionError(h1); 
    //ASSERT_LE(error, Traits::tolerance<scalar>::value())
       //<< "volume error = " << error << "\n"; 

    //h1.setPosition(vector(1,1,1));
    //h1.setDirection(vector(0,0,-1));
    //error = intersectionError(h1); 
    //ASSERT_LE(error, Traits::tolerance<scalar>::value())
       //<< "volume error = " << error << "\n"; 

    //h1.setDirection(vector(0,-1,0));
    //error = intersectionError(h1); 
    //ASSERT_LE(error, Traits::tolerance<scalar>::value())
       //<< "volume error = " << error << "\n"; 

    //h1.setDirection(vector(-1,0,0));
    //error = intersectionError(h1); 
    //ASSERT_LE(error, Traits::tolerance<scalar>::value())
       //<< "volume error = " << error << "\n"; 
    //error = intersectionError(h1); 
    //ASSERT_LE(error, Traits::tolerance<scalar>::value())
       //<< "volume error = " << error << "\n"; 
//}


//TEST(triangulationHalfspaceIntersection, successiveRegularIntersection)
//{
    //typedef AABBox<point> aabbox; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef TetrahedronSequence<tetrahedron, std::vector> triangulation;
    //typedef TriangulationIntersection<triangulation, polyhedron> triangulationIntersection;

    //polygon basePolygon = {point(0,0,0), point(1,0,0), point(1,1,0), point(0,1,0)};

    //vector sweepVector (1,1,1); 

    //polyhedron poly = build<polyhedron>(basePolygon, sweepVector);  

    //write_vtk_polydata(poly, fullTestName() + ".sweptPoly.vtk"); 

    //auto tri = barycentric_triangulate<triangulation>(poly); 

    //write_vtk_polydata(tri, fullTestName() + ".sweptTri.vtk"); 

    //polyhedron cutter = build<polyhedron>(aabbox(point(1,1,0.5), point(1.5,1.5,1))); 

    //write_vtk_polydata(cutter, fullTestName() + ".cutter.vtk"); 

    //auto intersection = intersect<triangulationIntersection>(tri, cutter); 

    //write_vtk_polydata(intersection, fullTestName() + ".result.vtk"); 

    //ASSERT_LE(mag(volume(intersection) - pow(0.5,3)), tolerance<scalar>::value())
       //<< "Resulting volume : " << volume(intersection) << "\n"
       //<< "Expected volume : " << pow(0.5,3) << "\n"; 
//}

//TEST(triangulationHalfspaceIntersection, incrementalIntersection)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef Halfspace<point, vector> halfspace;
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef TetrahedronSequence<tetrahedron, std::vector> triangulation;
    //typedef TriangulationIntersection<triangulation, polyhedron> triangulationIntersection;

    //double a = 1;
    //double b = 1;
    //double h = 1;

    //triangulation tri = {
        //build<tetrahedron>(point(0,0,0), point(a,0,0), point(0,b,0), point(0,0,h)),
        //build<tetrahedron>(point(a,0,0), point(a+1,0,0), point(a,b,0), point(a,0,h))
    //}; 

    //double exactVolume = tri.size() * (1. / 6. * (a * b * h )); 

    //const unsigned int intersectionCount = 100; 

    //for (unsigned int I = 0; I <= intersectionCount; ++I)
    //{
        //double hTrunc = I * 1. / intersectionCount;  

        //halfspace hspace (point(0,0,hTrunc), vector(0,0, -1)); 

        //// Intersect the triangulation with the halfspace.
        //auto result = intersect<triangulationIntersection>(tri, hspace); 

        //write_vtk_polydata(
            //result, 
            //fullTestName(), 
            //I
        //);

        //// Test the volume calculation.
        //auto resultVolume = volume(result); 

        //ASSERT_TRUE(resultVolume >= 0); 

        //auto exactResultVolume = exactVolume -  tri.size() * (pow(h - hTrunc, 3) * a * b)  / (6 * h * h);

        //auto volumeError = mag(resultVolume - exactResultVolume); 
        
        //ASSERT_TRUE(volumeError < Traits::tolerance<scalar>::value())
            //<< "volumeError = " << volumeError << "\n"; 

        //// Intersect the triangulation intersection with the same halfspace.
        //auto repeatedIntersectionResult = intersect<triangulationIntersection>(result, hspace);

        //write_vtk_polydata(
            //repeatedIntersectionResult, 
            //fullTestName() + ".repeated", 
            //I
        //);

        //auto repeatedIntersectionVolume = volume(repeatedIntersectionResult); 

        //auto repeatedIntersectionVolError = mag(repeatedIntersectionVolume - exactResultVolume); 

        //ASSERT_TRUE(repeatedIntersectionVolError < Traits::tolerance<scalar>::value()) 
            //<< "repeatedIntersectionVolError = " << repeatedIntersectionVolError << "\n";

    //}
//}

//TEST(triangulationHalfspaceIntersection, sideIntersectionPerturbed)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef Halfspace<point, vector> halfspace;
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef TetrahedronSequence<tetrahedron, std::vector> triangulation;
    //typedef TriangulationIntersection<triangulation, polyhedron> triangulationIntersection;

    //double a = 1. / 3.;

    ////EPSILON=1e-10;

    //polygon basePolygon = {point(0,0,0), point(a,0,0), point(a,a,0), point(0,a,0)};

    //auto test_result = [](
        //const triangulationIntersection& result, 
        //const triangulation& sweptTriangulation
    //) -> bool  
    //{
        //double expectedVolume = volume(sweptTriangulation); 
        //double computedVolume = volume(result); 
        //double volumeError = computedVolume - expectedVolume;  

        //// In the worst case scenario, for random perturbation of the entire polyhedron, 
        //// 6 x length x width x epsilon will be the increase in the volume on all sides.
        //// TM.
        //bool volumeIsValid = mag(volumeError) <= 6 * EPSILON;

        //if (! volumeIsValid)
        //{
            //write_vtk_polydata(
                //result, 
                //fullTestName() + ".failedResult.vtk"
            //);

            //for (const auto& tetrahedron : result.tetrahedra())
            //{
                //if (tetrahedron.size() > 4)
                //{
                    //write_vtk_polydata(tetrahedron, fullTestName() + ".wrongTetrahedron.vtk");
                //}
            //}

            //for (const auto& polyhedron : result.polyhedra())
            //{
                //if (polyhedron.size() > 4)
                //{
                    //write_vtk_polydata(polyhedron, fullTestName() + ".wrongPolyhedron.vtk");
                //}
            //}
        //}

        //EXPECT_TRUE(volumeIsValid)
            //<< "expected volume : " << expectedVolume << "\n"  
            //<< "computed volume : " << computedVolume  << "\n"
            //<< "relative volume error : " << volumeError / expectedVolume << "\n"
            //<< "polyhedra count : first result : " << result.polyhedra_size() << "\n"   
            //<< "tetrahedra count : first result : " << result.tetrahedra_size() << "\n";  

        //return volumeIsValid;
    //};

    //point planePoint (1/3., 0 , 0); 
    //auto originalPlanePoint = planePoint;
    //vector sweepVector(0, 0, 0.1/.3);
    //polyhedron sweptPoly = build<polyhedron>(basePolygon, sweepVector);  
    //// Uncomment for additional data output. 
    ////write_vtk_polydata(sweptPoly, fullTestName() + ".sweptPolyhedron.vtk"); 
    //auto sweptTriangulation = barycentric_triangulate<triangulation>(sweptPoly);
    //auto originalSweptTriangulation = sweptTriangulation;

    //for (unsigned int i = 0; i < MAX_IT; ++i)
    //{
        //perturb(planePoint, EPSILON); 
        //halfspace h1(planePoint, vector(-1,0,0)); 

        //perturb(sweptTriangulation, EPSILON); 

        //// Uncomment for additional data output.  
        ////write_vtk_polydata( // WARNING: large amount of data.
            ////sweptTriangulation, 
            ////fullTestName() + 
            ////".sweptTriangulation.vtk"
        ////); 

        //auto result = intersect<triangulationIntersection>(sweptTriangulation, h1, EPSILON); 

        //bool successful = test_result(result, sweptTriangulation);
        //if (!successful)
        //{
            //break;
        //}

        //// Repeat the intersection for triangulation intersection. 
        //result = intersect<triangulationIntersection>(sweptTriangulation, h1, EPSILON); 

        //// Test the result.
        //successful = test_result(result, sweptTriangulation);
        //if (!successful)
        //{
            //break;
        //}

        //planePoint = originalPlanePoint;
        //sweptTriangulation = originalSweptTriangulation; 
    //}
//}

//TEST(triangulationPolyhedronIntersection, starIntersection)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef TriangulationIntersection<tetrahedronVector, polyhedron> triangulationIntersection;

    //// Build a box_star polyhedron.
    //point pb0 (-0.5, -0.5, -0.5); 
    //point pb1 (0.5, 0.5, 0.5); 
    //auto starPolyhedron = build_box_star<polyhedron>(pb0, pb1);

    //// Triangulate a box_star polyhedron.
    //auto starTriangulation = oriented_triangulate<tetrahedronVector>(starPolyhedron);

    //write_vtk_polydata(starTriangulation, fullTestName() + ".starTriangulation.vtk"); 
    //// Iterative box star cutting.  

    //point p0(-1.5, -1.5, -1.5); 
    //point p1(1.5, 1.5, 1.5); 

    //unsigned int N = 100;
    //vector delta (0.01, 0.01, 0.01); 

    //std::ofstream dataFile;
    //dataFile.open(fullTestName() + ".dat"); 

    //// Initialize pyramid side and height.
    //double a = 1; 
    //double ha = 1;

    //for (unsigned int I = 0; I <= N; ++I)
    //{
        //auto h = I * delta; 

        //auto t0 = p0 + h; 
        //auto t1 = p1 - h; 
        //auto starTriangulationCutter = build<polyhedron>(AABBox<point>(t0, t1));

        //// Write the cutter into VTK for each cut iteration.
        //write_vtk_polydata(
            //starTriangulationCutter, 
            //fullTestName(), 
            //I
        //);

        //// Compute the volume of the truncated pyramids + the volume of the base box of the star. 
        //double hx = h[0]; 
        //double pyramidsVolumeExact = 1 + ((2 * pow(a, 2) * ha)  -  (6 * a * pow(hx,3) / (3 * ha))); 

        //// Intersect the triangulation with the cutter.
        //auto intersection = intersect<triangulationIntersection>(starTriangulation, starTriangulationCutter); 

        //write_vtk_polydata(
            //intersection,
            //fullTestName() + ".intersection", 
            //I
        //);

        //double pyramidsVolume = volume(intersection); 

        //// Output the volume values into a file.
        //dataFile << hx << " " << pyramidsVolumeExact << " " << pyramidsVolume << "\n"; 

        //// Test the intersection volume against the exact volume.
        //ASSERT_LE(mag(pyramidsVolume - pyramidsVolumeExact), 5*tolerance<scalar>::value())
            //<< "pyramidsVolume : " << pyramidsVolume << "\n" 
            //<< "pyramidsVolumeExact : " << pyramidsVolumeExact << "\n" 
            //<< "diff : " << mag(pyramidsVolume - pyramidsVolumeExact) << "\n"; 
    //}
//}

// FIXME: Extract common parts into a lamda function and run the perturbation in
// a for loop for different polyhedra that have 0 intersection. TM.  Lambda
// arguments : point1, point2, tolerance, nTetra, nPolyhedra.  TM.
//TEST(triangulationPolyhedronIntersection, triangulationInsideOutsidePolyhedron)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef TetrahedronSequence<tetrahedron, std::vector> triangulation;
    //typedef TriangulationIntersection<triangulation, polyhedron> triangulationIntersection;

    //double tolerance = EPSILON; 

    //// FIXME: Generalize the point coordinates, go away from 0,1. TM.
    //auto polyInside = build<polyhedron>(point(0,0,0),point(1,1,1));
    //auto triInside = barycentric_triangulate<triangulation>(polyInside); 

    //auto polyOutside = build<polyhedron>(point(1,0,0),point(2,1,1));
    //auto triOutside = triInside;  

    //for (unsigned int i = 0; i < MAX_IT; ++i)
    //{
        //// FIXME: Using nonuniform perturbation introduces errors larger than 10 * EPSILON  
        //// for a perturbed polyhedron. The reason is the triangle area normal vector that
        //// is computed using an averaging operation. When the differences in point positions
        //// are EPSILON-sized, the averaging results with roundoff errors that is larger than
        //// 1e-15.  
        //// TODO: To avoid this, test the perturbed polygon normal calculation for a square.
        
        ////Uncomment this one to enable output of perturbation vectors. 
        ////perturb(polyInside, EPSILON, [](const point& p) {Info << p << endl;}); 
        ////perturb(polyInside, EPSILON);  

        //// Using uniform perturbation.
        //auto perturbationVec = perturbation_vector<vector>(tolerance); 
        //polyInside = displace(polyInside, perturbationVec); 
        
        //ASSERT_TRUE(is_inside(polyInside, triInside, EPSILON))
            //<< "perturbationVec : " << perturbationVec << "\n"; 

        ////Info << "Intersect triangulation with the polyhedron ..." << endl;
        //auto intersectionInside = intersect<triangulationIntersection>(triInside, polyInside, tolerance); 
        ////Info << "Done." << endl;

        ////Info << "Compute volume diff - polyhedron and triangualtion ..." << endl;
        //double polyInsideTriVolDiff = volume(polyInside) - volume(triInside); 
        ////Info << "Done." << endl;

        //ASSERT_LE(mag(polyInsideTriVolDiff), 10 * tolerance)
           //<< "polyInsideTriVolDiff : " 
           //<< polyInsideTriVolDiff << "\n"; 

        ////Info << "Compute volume diff - intersection and triangulation" << endl;
        ////double insideVolDiff = volume(intersectionInside) - volume(triInside); 
        ////Info << "Done." << endl;

        ////bool triangulationInsidePolyhedron = mag(insideVolDiff) <= 10 * tolerance;

        ////if (! triangulationInsidePolyhedron) 
        ////{
            ////write_vtk_polydata(polyInside, fullTestName() + ".polyhedronInside.vtk"); 
            ////write_vtk_polydata(triInside, fullTestName() + ".triangulationInside.vtk"); 
            ////write_vtk_polydata(intersectionInside, fullTestName() + ".intersectionInside.vtk"); 
        ////}

        //// Topological tests. The intersectionInside should hold 24 tetrahedra
        //// when barycentric decomposition is applied. TM.
        ////ASSERT_TRUE(intersectionInside.tetrahedra_size() == 24) 
           ////<< "Tetrahedra in the intersection : " 
           ////<< intersectionInside.tetrahedra_size() << "\n";
        //// The intersectionInside should hold no polyhedra. 
        ////ASSERT_TRUE(intersectionInside.polyhedra_size() == 0)
           ////<< "Polyhedra in the intersection : " 
           ////<< intersectionInside.polyhedra_size() << "\n";

        //// Volume test.
        ////ASSERT_TRUE(mag(insideVolDiff) <= 10 * tolerance) 
           ////<< "Volume difference between the result and the cutting polyhedron: " 
           ////<< insideVolDiff << "\n"; 


        //// Outside test

        //// Perturb the polyhedron that doesn't contain the triangulation.
        ////perturb(polyOutside, EPSILON); 

        ////Info << "Is inside assertion.." << endl;
        ////ASSERT_FALSE(is_inside(polyOutside, triOutside, EPSILON)); 
        ////Info << "Done." << endl;

        ////auto intersectionOutside = intersect<triangulationIntersection>(triOutside, polyOutside); 

        ////bool triangulationOutsidePolyhedron = 
            ////(intersectionOutside.tetrahedra_size() == 0) && 
            ////(intersectionOutside.polyhedra_size() == 0) &&
            ////(mag(volume(intersectionOutside)) <= tolerance);

        ////if (! triangulationOutsidePolyhedron) 
        ////{
            ////write_vtk_polydata(polyOutside, fullTestName() + ".polyhedronOutside.vtk"); 
            ////write_vtk_polydata(triOutside, fullTestName() + ".triangulationOutside.vtk"); 
            ////write_vtk_polydata(intersectionOutside, fullTestName() + ".intersectionOutside.vtk"); 
        ////}

        ////double polyOutsideTriVolDiff = volume(polyOutside) - volume(triOutside); 

        ////ASSERT_LE(mag(polyOutsideTriVolDiff), tolerance)
           ////<< "Volume difference between the polyhedron and the triangulation: " 
           ////<< polyOutsideTriVolDiff << "\n"; 

        ////// Topological tests.
        ////// The intersectionOutside should hold no tetrahedra. TM.
        ////ASSERT_TRUE(intersectionOutside.tetrahedra_size() == 0)
           ////<< "Tetrahedra in the intersection : " 
           ////<< intersectionOutside.tetrahedra_size() << "\n";
        ////// The intersectionOutside should hold no polyhedra. 
        ////ASSERT_TRUE(intersectionOutside.polyhedra_size() == 0)
           ////<< "Polyhedra in the intersection : " 
           ////<< intersectionOutside.polyhedra_size() << "\n";

        ////// Volume test.
        ////ASSERT_TRUE(mag(volume(intersectionOutside)) <= tolerance) 
           ////<< "Volume of the outside intersection should be zero: " 
           ////<< volume(intersectionOutside) << "\n"; 

    //}
//}

//TEST(triangulationPolyhedronIntersection, cornerIntersection)
//{
    //typedef AABBox<point> aabb_box; 
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef TetrahedronSequence<tetrahedron, std::vector> triangulation;
    //typedef TriangulationIntersection<triangulation, polyhedron> triangulationIntersection;

    //auto cutBox = aabb_box(point(0,0,0), point(0.75, 0.75, 1)); 
    //auto cutPoly = build<polyhedron>(cutBox);
    //auto triBox = aabb_box(point(0, 0, 0), point(1, 1, 1)); 
    //auto triPoly = build<polyhedron>(triBox);
    //auto tri = barycentric_triangulate<triangulation>(triPoly); 

    //auto intersection = intersect<triangulationIntersection>(tri, cutPoly); 

    //write_vtk_polydata(cutPoly, fullTestName() + ".cuttingPolyhedron.vtk"); 
    //write_vtk_polydata(tri, fullTestName() + ".triangulation.vtk"); 
    //write_vtk_polydata(intersection, fullTestName() + ".intersection.vtk"); 

    //double tolerance = Traits::tolerance<scalar>::value(); 
    //double boxCutPolyVolDiff = volume(cutBox) - volume(cutPoly); 

    //ASSERT_LE(mag(boxCutPolyVolDiff), tolerance)
       //<< "Volume difference between the box and the cutting polyhedron : " << boxCutPolyVolDiff << endl; 

    //double resultVolDiff = volume(intersection) - volume(cutPoly); 
    //ASSERT_TRUE(mag(resultVolDiff) <= tolerance) 
       //<< "Volume difference between the result and the cutting polyhedron: " << resultVolDiff << endl; 
//}

// Create a swept polyhedron from one of the base polyhedron polygons. 
// Make sure that the swept polyhedron is created so that it stays within the base polyhedron. 
// Loop MAX_IT times
    // **Perturb the base polyhedron by EPSILON.**
    // Check if the triangulation of the swept polyhedron is inside of the perturbed base polyhedron.
//TEST(triangulationPolyhedronIntersection, perturbedIsInsideTest)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;

    //// For the triangulation intersection operations, only elements within [-1e-15, 1e15] are 
    //// disregarded for accuracy reasons. TM.
    ////EPSILON = 1e-15; 

    //for (unsigned int i = 0; i < MAX_IT ; ++i)
    //{
        //auto basePolyhedron = build<polyhedron>(point(0,0,0), point(1,1,1));

        //const polygon& sweptPolygon = basePolyhedron.front(); 
        //vector sweepVector(0,0, -1);
        //auto sweptPolyhedron = build<polyhedron>(sweptPolygon, sweepVector);  

        //// Perturb the position of the basePolyhedron points within [-EPSILON, EPSILON] 
        
            //// Uncomment to output perturbed points for debugging.
        ////perturb(basePolyhedron, EPSILON,[](point& p) { Info << p << endl; }); 
        
            //// Uncomment to just perturb points without output.
        //perturb(basePolyhedron, EPSILON);

        //auto basePolyVol = volume(basePolyhedron); 
        //auto sweptPolyVol = volume(sweptPolyhedron); 
        //auto polyhedronVolDiff = mag(basePolyVol - sweptPolyVol); 

        //// A unit box, perturbed by EPSILON may have an error in volume  up to
        //// (1 + 2 * EPSILON)^2 * EPSILON per single faces. Six faces lead to 6
        //// * EPSILON. TM.
        //bool polyhedronVolDiffOk = polyhedronVolDiff <=  6 * EPSILON; 

        //if (!polyhedronVolDiffOk)
        //{
            //write_vtk_polydata(
                //sweptPolyhedron,
                //fullTestName() + ".sweptPolyhedron",
                //i
            //);
            //write_vtk_polydata(
                //basePolyhedron, 
                //fullTestName() + ".basePolyhedron",
                //i
            //);
        //}

        //ASSERT_TRUE(polyhedronVolDiffOk) 
           //<< "Base polyhedron volume : " << basePolyVol << "\n"
           //<< "Swept polyhedron volume : " << sweptPolyVol << "\n"
           //<< "Difference : " << polyhedronVolDiff << "\n" 
           //<< "Tolerance : " << EPSILON << endl; 
    //}
//}

/// @brief Test the intersection between a polyhedron stencil and a polyhedron triangulation. 
/// Sweep a polygon to create a polyhedron. Triangulate the swept polyhedron with different
/// triangulation strategies. Intersect the triangulation with the stencil of adjacent polyhedrons. 
/// Compute contributions from each individiual intersection. The total volume must be equal to 
/// the volume of the initial triangulation. 
/// 
/// Tested for : barycentric triangulation, flux triangulation.
///
/// @param triangulationPolyhedronIntersection
/// @param perturbedPolyhedronStencilTest
//TEST(triangulationPolyhedronIntersection, perturbedPolyhedronStencilTest)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef TetrahedronSequence<tetrahedron, std::vector> triangulation;
    //typedef TriangulationIntersection<triangulation, polyhedron> triangulationIntersection;

    //EPSILON=1e-15;
    //MAX_IT = 1e04;

    //auto basePolyhedron = build<polyhedron>(point(0,0,0), point(1,1,1));
    //auto originalBasePolyhedron = basePolyhedron;

    //const polygon& sweptPolygon = basePolyhedron.front(); 
    //vector sweepVector(0,0, -1);
    //auto sweptPolyhedron = build<polyhedron>(sweptPolygon, sweepVector);  

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

    //for (unsigned int i = 0; i < MAX_IT ; ++i)
    //{
        //perturb(basePolyhedron, EPSILON);

        //auto basePolyVol = volume(basePolyhedron); 
        //auto sweptPolyVol = volume(sweptPolyhedron); 
        //auto polyhedronVolDiff = mag(basePolyVol - sweptPolyVol); 

        //// FIXME: Check this! 
        //// A unit box, perturbed by EPSILON may have an error in volume  up to
        //// (1 + 2 * EPSILON)^2 * EPSILON per single faces. Six faces lead to 6
        //// * EPSILON. TM.
        //bool polyhedronVolDiffOk = polyhedronVolDiff <=  6 * EPSILON; 

        //if (!polyhedronVolDiffOk)
        //{
            //write_vtk_polydata(
                //sweptPolyhedron,
                //fullTestName() + ".sweptPolyhedron",
                //i
            //);
            //write_vtk_polydata(
                //basePolyhedron, 
                //fullTestName() + ".basePolyhedron",
                //i
            //);
        //}

        //ASSERT_TRUE(polyhedronVolDiffOk) 
           //<< "Base polyhedron volume : " << basePolyVol << "\n"
           //<< "Swept polyhedron volume : " << sweptPolyVol << "\n"
           //<< "Difference : " << polyhedronVolDiff << "\n" 
           //<< "Tolerance : " << EPSILON << endl; 

        //// Triangulate the flux polyhedron.
        //auto sweptTriangulation = barycentric_triangulate<triangulation>(sweptPolyhedron);

        //auto sweptTriVol = volume(sweptTriangulation);  
        //auto triVolDiff = mag(sweptTriVol - sweptPolyVol); 

        //ASSERT_LE(triVolDiff, EPSILON)
           //<< "Swept triangulation volume : " << sweptTriVol << "\n"
           //<< "Swept polyhedron volume : " << sweptPolyVol << endl;

        //// Build a "cell" stencil.
        //std::vector<polyhedron> stencil; 
        //stencil.push_back(basePolyhedron); 
        //stencil.push_back(displace(basePolyhedron, vector(1,0,0))); 
        //stencil.push_back(displace(basePolyhedron, vector(1,1,0))); 
        //stencil.push_back(displace(basePolyhedron, vector(0,1,0))); 

        //// Compute volume contributions.
        //double totalVolume = 0;  
        //for (const auto& polyhedron : stencil)
        //{
            //auto intersection = intersect<triangulationIntersection>(sweptTriangulation, polyhedron, EPSILON);
            //totalVolume += volume(intersection); 
        //}

        //double totalVolDiff = mag(totalVolume - sweptTriVol); 
        //bool volumesAreDifferent = totalVolDiff > 6*EPSILON; 

        //if (volumesAreDifferent)
        //{
            //write_vtk_polydata(sweptTriangulation, fullTestName() + ".triangulation", i); 
            //write_vtk_polydata(sweptPolyhedron, fullTestName() + ".polyhedron", i); 
            
            //// Visualize the intersection results.
            //unsigned int polyhedronCounter = 0;  
            //for (const auto& polyhedron : stencil)
            //{
                //auto intersection = intersect<triangulationIntersection>(sweptTriangulation, polyhedron, EPSILON);

                //Info << "Polyhedron counter : " << polyhedronCounter << "\n" 
                    //<< " Intersection volume: " << volume(intersection) << "\n";

                //write_vtk_polydata(polyhedron, fullTestName() + ".polyhedron", polyhedronCounter); 
                //write_vtk_polydata(intersection, fullTestName() + ".intersection", i, polyhedronCounter); 

                //// Uncomment to perform topological checks. Used for debugging the intersection algorithm. 
                //unsigned int polyhedronID = 0; 
                //for (const auto& polyhedron : intersection.polyhedra())
                //{
                    //// All polyhedra are tetrahedra with 4 sides.
                    //if (polyhedron.size() != 4)
                    //{
                        //write_vtk_polydata(polyhedron, fullTestName() + ".faultyTetrahedron", polyhedronID);
                    //}
                    //++polyhedronID;
                //}
                //++polyhedronCounter; 
            //}
        //}

        //ASSERT_FALSE(volumesAreDifferent)
           //<< "Total volume : " << totalVolume << "\n"
           //<< "Swept volume : " << sweptPolyVol << "\n"
           //<< "% volume difference : " << (totalVolDiff / sweptPolyVol) * 100 << "\n";

        //basePolyhedron = originalBasePolyhedron;
    //}
//}

//TEST(triangulationPolyhedronIntersection, pointIntersection)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef TetrahedronSequence<tetrahedron, std::vector> triangulation;
    //typedef TriangulationIntersection<triangulation, polyhedron> triangulationIntersection;

    //EPSILON=1e-10;
    //MAX_IT = 1e04;

    //auto basePolyhedron = build<polyhedron>(point(0,0,0), point(1,1,1));
    //auto originalBasePolyhedron = basePolyhedron;

    //double l = 1;  
    //double h = 0.1; 

    //polygon sweptPolygon = {point(0,0,0), point(0,l,0), point(l,l,0), point(l,0,0)}; 

    //vector sweepVector(h, h, h);
    //auto sweptPolyhedron = build<polyhedron>(sweptPolygon, sweepVector);  
    //auto exactVolume = (l * l) * h;  

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

    //auto intersectionIsValid = [&](const auto& doTriangulate, word triName)
    //{
        //auto sweptTri = doTriangulate(sweptPolyhedron); 
        //write_vtk_polydata(sweptTri, fullTestName() + "." + triName + "Triangulation.vtk"); 
        //auto sweptTriIntersect = intersect<triangulationIntersection>(sweptTri, basePolyhedron);
        //return (mag(volume(sweptTri) - exactVolume) <= EPSILON);
    //};

    //ASSERT_TRUE(intersectionIsValid(flux_triangulation<triangulation>(), "flux"));  
    //ASSERT_TRUE(intersectionIsValid(barycentric_triangulation<triangulation>(), "barycentric"));  
//}

//TEST(triangulationPolyhedronIntersection, triangulationInsidePolyhedron)
//{
    //typedef PointSequence<point, std::vector> polygon;
    //typedef PolygonSequence<polygon, std::vector> polyhedron;
    //typedef ArrayTriangle<point> triangle;
    //typedef ArrayTetrahedron<triangle> tetrahedron; 
    //typedef TetrahedronSequence<tetrahedron, std::vector> triangulation;
    //typedef TriangulationIntersection<triangulation, polyhedron> triangulationIntersection;

    //auto basePolyhedron = build<polyhedron>(point(0,0,0), point(1,1,1));
    //auto originalBasePolyhedron = basePolyhedron;

    //double l = 1;  
    //double h = 0.1; 

    //polygon sweptPolygon = {point(0,0,0), point(0,l,0), point(l,l,0), point(l,0,0)}; 

    //vector sweepVector(0, 0, h);

    //auto sweptPolyhedron = build<polyhedron>(sweptPolygon, sweepVector);  
    //auto exactVolume = (l * l) * h;  

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

    //auto intersectionIsValid = [&](const auto& doTriangulate, word triName)
    //{
        //auto sweptTri = doTriangulate(sweptPolyhedron); 

        //write_vtk_polydata(sweptTri, fullTestName() + "." + triName + "Triangulation.vtk"); 

        //auto sweptTriIntersect = intersect<triangulationIntersection>(sweptTri, basePolyhedron);

        //return (mag(volume(sweptTri) - exactVolume) <= EPSILON);
    //};

    //ASSERT_TRUE(intersectionIsValid(flux_triangulation<triangulation>(), "flux"));  
    //ASSERT_TRUE(intersectionIsValid(barycentric_triangulation<triangulation>(), "barycentric"));  

    //auto originalPolyhedron = basePolyhedron;

    //perturb(basePolyhedron, EPSILON); 

    //ASSERT_TRUE(intersectionIsValid(flux_triangulation<triangulation>(), "flux"));  
    //ASSERT_TRUE(intersectionIsValid(barycentric_triangulation<triangulation>(), "barycentric"));  

    //basePolyhedron = originalPolyhedron;
//}
    
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

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


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