/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

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

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

Application
    geomTransportTestWriteVtk 

Description
    Test the build algorithms. 

Authors
    Tomislav Maric maric@csi.tu-darmstadt.de, tomislav@sourceflux.de
    Mathematical Modeling and Analysis
    Center of Smart Interfaces, TU Darmstadt

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

//#include <AxisAlignedBoundingBox.hpp>

#include "tag.hpp"
#include "AxisAlignedBoundingBox.hpp"

//#include "fvCFD.H"
//#include "geomIOhelperFunctions.H"
//#include "Geometry.H"
//#include "testUtilities.H"
//
//using namespace GeometricalTransport;
//const double EPSILON = 1e-15;  
//const unsigned int MAX_IT = 1e06; 

#include "gtest/gtest.h" 

//TEST(visualization, polyDataFieldOutput)
//{
    //typedef pointVectorVector polyhedron; 

    //auto p1 = build<polyhedron>(point(0,0,0), point(1,1,1)); 
    //auto p2 = build<polyhedron>(point(1,0,0), point(2,1,1)); 
    //auto p3 = build<polyhedron>(point(2,0,0), point(3,1,1)); 
    
    //// Prepare the stream and the output file for flux polyhedra.
    //vtk_polydata_stream polyhedronStream(fullTestName() + ".vtk"); 

    //polyhedronStream.streamGeometry(p1); 
    //polyhedronStream.streamGeometry(p2); 
    //polyhedronStream.streamGeometry(p3); 

    //scalarField s(3,0); 
    //s[0] = 5; 
    //s[1] = 6; 
    //s[2] = 7; 

    //polyhedronStream.streamField(s, "values"); 
//}

//TEST(perturbedVector, halfspace)
//{
    //vector O(0,0,0); // Origin vector.

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

    //for (unsigned int i = 0; i < MAX_IT; ++i)
    //{
        //perturb(O, EPSILON); // Perturb the origin within epsilon-ball.

        //EXPECT_LE(mag(O[0]), EPSILON) 
            //<<  "O : " << O << "\n";

        //EXPECT_LE(mag(O[1]), EPSILON) 
            //<<  "O : " << O << "\n";

        //EXPECT_LE(mag(O[2]), EPSILON) 
            //<<  "O : " << O << "\n";

        //EXPECT_LE(distance(O,h), EPSILON) << 
           //"distance(O,h) : " <<  distance(O,h) << "\n"
           //"O : " << O << "\n";

        //O = vector(0,0,0); 
    //}
//}

//TEST(pointSequence, halfspace)
//{
    //halfspace plane(point(1,1,1), point(1,1,1)); 

    //auto poly = build<pointVector>(plane); 

    //write_vtk_polydata(plane, "halfspace.vtk"); 

    //// Polygon normal area vector and plane normal are collinear.
    //EXPECT_LE(mag(plane.direction() ^ normal_area_vec(poly)), EPSILON); 

    //// Polygon normal area vector and plane normal have the same direction.
    //EXPECT_GE(plane.direction() & normal_area_vec(poly), 0); 

    //// All polygon points lie on the plane.
    //for (const auto& point : poly)
    //{
        //auto pointToPlaneDistance = distance(point, plane);  
        //EXPECT_TRUE(
            //(pointToPlaneDistance > -EPSILON) && 
            //(pointToPlaneDistance < EPSILON)
        //) << "pointToPlaneDistance : " << pointToPlaneDistance << "\n";  
    //}
//}

//TEST(perturbedPolygon, centroid)
//{
    //halfspace plane(point(1,1,1), point(1,1,1)); 

    //auto poly = build<pointVector>(plane); 

    //write_vtk_polydata(plane, "halfspace.vtk"); 

    //// Polygon normal area vector and plane normal are collinear.
    //EXPECT_LE(mag(plane.direction() ^ normal_area_vec(poly)), EPSILON); 

    //// Polygon normal area vector and plane normal have the same direction.
    //EXPECT_GE(plane.direction() & normal_area_vec(poly), 0); 

    //// All polygon points lie on the plane.
    //for (const auto& point : poly)
    //{
        //auto pointToPlaneDistance = distance(point, plane);  
        //EXPECT_TRUE(
            //(pointToPlaneDistance > -EPSILON) && 
            //(pointToPlaneDistance < EPSILON)
        //) << "pointToPlaneDistance : " << pointToPlaneDistance << "\n";  
    //}
//}


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

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

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

    //write_vtk_polydata(mixed, "buildMixedPolyhedron.base.mixed.vtk");

    //EXPECT_TRUE(mixed.triangles_size() == 4);
    //EXPECT_TRUE(mixed.polygons_size() == 5);
//}

//TEST(aabbox, expand)
//{
    //aabbox b; 

    //pointVector points = {point(0,0,0),point(1,0,0),point(1,1,1)};

    //b.expand(points); 

    //EXPECT_TRUE(b.minPoint() == point(0,0,0)) 
        //<< b.minPoint() << "\n"; 

    //EXPECT_TRUE(b.maxPoint() == point(1,1,1))
        //<< b.maxPoint() << "\n"; 
//}

//TEST(openFoamBuild, meshEntities)
//{
    //extern int mainArgc; 
    //extern char** mainArgv; 

    //int argc = mainArgc; 
    //char** argv = mainArgv; 

    //#include "setRootCase.H"
    //#include "createTime.H"
    //#include "createMesh.H"

    //const faceList& meshFaces = mesh.faces(); 

    //const surfaceVectorField& Sf = mesh.Sf(); 
    //forAll(meshFaces, I)
    //{
        //// Build a polygon from a face and test the normal.
        //auto facePolygon = build<pointVector>(meshFaces[I], mesh); 
        //EXPECT_LE(mag(normal_area_vec(facePolygon) - Sf[I]), EPSILON);  

        //// Build a polygon from a flipped face and test the normal.
        //face flipFace(meshFaces[I]); 
        //flipFace.flip(); 
        //auto flipFacePolygon = build<pointVector>(flipFace, mesh); 
        //EXPECT_LE(mag(normal_area_vec(flipFacePolygon) + Sf[I]), EPSILON);  

        //// Build a polyhedron from a polygon and an OpenFOAM vector/pointField. 
        //vectorField displacements(facePolygon.size(), Sf[I] / mag(Sf[I])); 

        //// Test the sweep build. Use both orientations for displacements.
        //auto sweptFacePolyhedron = build<pointVectorVector>(facePolygon, displacements);

        //auto volumeError = mag(volume(sweptFacePolyhedron) - mag(Sf[I])); 
        //if (volumeError > EPSILON)
            //write_vtk_polydata(sweptFacePolyhedron, fullTestName() + ".badSweptFacePolyhedron.vtk"); 
        //EXPECT_LE(volumeError, EPSILON); 

        //displacements *= -1; 
        //sweptFacePolyhedron = build<pointVectorVector>(facePolygon, displacements);

        //volumeError = mag(volume(sweptFacePolyhedron) - mag(Sf[I])); 
        //if (volumeError > EPSILON)
            //write_vtk_polydata(sweptFacePolyhedron, fullTestName() + ".badSweptFacePolyhedron.vtk"); 

        //EXPECT_LE(volumeError, EPSILON); 

    //}

    ////Build and test a halfspace collection from all mesh cells.
    //const auto& C = mesh.C();  
    //const auto& V = mesh.V();  
    //forAll(C, cellI)
    //{
        //auto cellHalfspaces = build<halfspaceVector>(cellI, mesh); 

        //// Test 
        //for (const auto& hspace : cellHalfspaces)
        //{
            //auto dist = distance(hspace, C[cellI]); 
            //if (dist < 0)
            //{
                //Info << "Distance = " << dist << endl;
                //Info << "Cell = " <<  cellI << endl;
                //write_vtk_polydata(build<pointVectorVector>(cellI, mesh), fullTestName() + ".cell.vtk");
                //write_vtk_polydata(hspace, fullTestName() + ".hspace.vtk"); 
                //Info << "Center = " << C[cellI] << endl;
                //Info << "Hspace position = " << hspace.position() << endl;
                //Info << "Hspace direction = " << hspace.direction() << endl;
            //}
            //EXPECT_TRUE(dist > 0);
        //}

        //auto polyhedron = build<pointVectorVector>(cellI, mesh); 

        //// Test polyhedron volume.
        //auto polyVolume = volume(polyhedron); 

        //auto volDiff = mag(polyVolume - V[cellI]); 

        //bool differentVolumes = volDiff > EPSILON; 

        //EXPECT_FALSE(differentVolumes);
        //if (differentVolumes)
        //{
            //Info << "Relative volume error = " << volDiff / polyVolume  << endl;
            //write_vtk_polydata(
                //barycentric_triangulate<tetrahedronVector>(polyhedron), 
                //"failedVolumeTri", cellI
            //);
            //write_vtk_polydata(polyhedron, "failedVolumePoly", cellI);
        //}

        //// Test polygon normal orientation.
        //for (const auto& polygon : polyhedron)
        //{
            //auto n = normal_area_vec(polygon); 
            //auto dist = (C[cellI] - polygon.front()) & n; 
            //if (dist > 0)
                //write_vtk_polydata(polyhedron, "geomTestBuild.misorientedPolyhedron", cellI);
            //EXPECT_TRUE(dist < 0);
        //}
    //}

//}

int mainArgc; 
char** mainArgv; 

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

int main(int argc, char *argv[])
{
    ::testing::InitGoogleTest(&argc, argv);
   
    mainArgc = argc; 
    mainArgv = argv;

    return RUN_ALL_TESTS();
}


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