/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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
    geomTestCorrectVolume
    

Description
    Test the volume correction algorithms. 

    Used by the correction for volume conservation. 

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

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

#include "fvCFD.H"

#include "Geometry.H"
#include "WriteVtkPolydata.H"

#include "gtest/gtest.h" 
#include "testUtilities.H"

using namespace GeometricalTransport;

const double EPSILON = 1e-15; 

// Correct the cube volume by 0.5 in the normal direction.
TEST(cube, simpleTest)
{
    typedef MixedPolyhedron<arrayTriangle, pointVector, std::vector> mixedPolyhedron;

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

    write_vtk_polydata(cube, fullTestName() + ".cube.vtk");  

    double volumeDifference = 0.5;

    auto correctedCube = correct_volume<mixedPolyhedron>(cube, cube.size() - 1, volumeDifference); 

    write_vtk_polydata(correctedCube, fullTestName() + ".correctedCube.vtk"); 

    ASSERT_LE(mag(volume(correctedCube) - volume(cube) - volumeDifference), Traits::tolerance<scalar>::value());
}

// Correct cube volume incrementally in the direction of the cap normal.
TEST(cube, incrementalCubeTest)
{
    typedef MixedPolyhedron<arrayTriangle, pointVector, std::vector> mixedPolyhedron;

    pointVector basePolygon = {point(0,0,0), point(0,1,0), point(1,1,0), point(1,0,0)};  

    vector sweepDirection(0,0,1); 
    sweepDirection /= mag(sweepDirection); 

    auto cube = build<pointVectorVector>(basePolygon, sweepDirection);  

    write_vtk_polydata(cube, fullTestName() + ".cube.vtk");

    double cubeVolume = volume(cube);  

    ASSERT_LE(mag(cubeVolume - 1), EPSILON)
       << "cubeVolume = " << cubeVolume << "\n"; 

    unsigned int MAX_IT = 1e06;
    std::ofstream errorFile; 
    errorFile.open(fullTestName() + ".dat"); 
    bool errorWithinTolerance = false; 
    double volumeError = 0;

    for (unsigned int i = 0; i < MAX_IT; ++i)
    {
        auto volumeDifference = 
            -(1.0 / 3.0) * cubeVolume + i * ((2.0 * cubeVolume) / MAX_IT); 

        auto correctedCube = correct_volume<mixedPolyhedron>(cube, cube.size() - 1, volumeDifference); 
        auto correctedCubeTri = flux_triangulate(correctedCube); 
        auto correctedCubeVolume = volume(correctedCubeTri); 

        volumeError = mag(correctedCubeVolume - cubeVolume - volumeDifference);

        errorFile << volumeDifference / cubeVolume  << " " << volumeError << "\n"; 

        errorWithinTolerance = (volumeError <= EPSILON); 

        if (!errorWithinTolerance)
        {
            write_vtk_polydata(correctedCubeTri, fullTestName() + ".cubeTriangulation", i);
            write_vtk_polydata(correctedCube, fullTestName() + ".correctedCube", i);
        }

        ASSERT_TRUE(errorWithinTolerance)
            << "cubeVolume = " << cubeVolume << "\n"
            << "volumeDifference = " << volumeDifference << "\n"
            << "correctedCubeVolume " << correctedCubeVolume << "\n"
            << "volumeError = " << volumeError << "\n";
    }
}

// FIXME: Either quantify proerly and select a triangulation that makes sense, or delete.
// Change angle of the volume correction for the cap face of the box polyhedron. 
//TEST(correctCubeVolume, directedIncrementalAngleTest)
//{
    //typedef MixedPolyhedron<arrayTriangle, pointVector, std::vector> mixedPolyhedron;

    //auto cube = build<pointVectorVector>(point(0,0,0), point(1,1,1));
    //// First polygon of the cube is the polygon with the z-coordinate of 1. 
    //const int polygonLabel = 0; 

    //write_vtk_polydata(cube, "correctCubeVolume.directedIncrementalAngleTest.cube.vtk"); 

    //double volumeDifference = 0.1;

    //const unsigned int MAX_IT = 1e04;
    //std::ofstream errorFile; 
    //errorFile.open("correctCubeVolume.directedIncrementalAngleTest.dat"); 

    //for (unsigned int increment = 0; increment < MAX_IT; ++increment)
    //{
        //double theta = increment * (3.14157 / MAX_IT); 
        //vector direction (0, 0.1 + std::cos(theta), 0.1 + std::sin(theta)); 

        //auto correctedCube = correct_volume<mixedPolyhedron>(cube, polygonLabel, direction, volumeDifference); 
        //auto correctedCubeTri = flux_triangulate(correctedCube); 
        //auto correctedCubeVolume = volume(correctedCubeTri); 

        //auto error = mag(correctedCubeVolume - volume(cube) - volumeDifference);  

        //if (error > EPSILON)
        //{
            //write_vtk_polydata(correctedCube, fullTestName() + ".correctedCube", increment);
            //write_vtk_polydata(correctedCubeTri, fullTestName() + ".cubeTriangulation", increment);
        //}

        //ASSERT_LE(error, EPSILON); 

        //errorFile << theta << " " << error << "\n"; 
    //}
//}

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

    pointVector basePolygon = {point(0,0,0), point(0,1,0), point(1,1,0), point(1,0,0)};  

    vector sweepDirection(1,1,1); 
    sweepDirection /= mag(sweepDirection); 

    auto prism = build<pointVectorVector>(basePolygon, sweepDirection);  

    // Last polygon of the prism is the polygon with the z-coordinate of 1. 
    const int polygonLabel = prism.size() - 1; 

    write_vtk_polydata(prism, "correctPrismVolume.incrementalNormalTest.prism.vtk"); 

    auto prismVolume = volume(prism); 
    std::ofstream errorFile; 
    errorFile.open("correctPrismVolume.incrementalNormalTest.dat"); 
    const unsigned int MAX_IT = 1e04;

    bool passed = false; 
    const double MAX_VOL_DIFF = 0.15; 

    for (unsigned int increment = 0; increment < MAX_IT; ++increment)
    {
        double volumeDifference = 
            (-MAX_VOL_DIFF * prismVolume) + (
                increment *
                ((MAX_VOL_DIFF * prismVolume) / MAX_IT)
            ); 

        auto correctedPrism = correct_volume<mixedPolyhedron>(prism, polygonLabel, volumeDifference); 
        auto correctedPrismTri = flux_triangulate(correctedPrism); 
        auto correctedPrismVol = volume(correctedPrismTri); 

        auto error = mag(correctedPrismVol - prismVolume - volumeDifference);  

        passed = error <= Traits::tolerance<scalar>::value(); 

        if (!passed)
        {
            write_vtk_polydata(correctedPrismTri, fullTestName(), increment);  
            errorFile << volumeDifference / prismVolume << " " << error << "\n"; 
        }

        ASSERT_TRUE(passed)
            << "prismVolume = " << prismVolume << "\n"
            << "volDifference = " << volumeDifference << "\n"
            << "correctedPrismVol = " << correctedPrismVol << "\n";
    }
}

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

    pointVector basePolygon = {point(0,0,0), point(0,1,0), point(1,1,0), point(1,0,0)};  

    vector sweepDirection (1,1,1); 
    sweepDirection /= mag(sweepDirection); 

    auto prism = build<pointVectorVector>(basePolygon, sweepDirection);  

    // Last polygon of the prism is the polygon with the z-coordinate of 1. 
    const int polygonLabel = prism.size() - 1; 

    write_vtk_polydata(prism, "correctPrismVolume.incrementalDirectedTest.prism.vtk"); 

    auto prismVolume = volume(prism); 
    std::ofstream errorFile; 
    errorFile.open("correctPrismVolume.incrementalDirectedTest.dat"); 
    const unsigned int MAX_IT = 1e04;

    for (unsigned int increment = 0; increment < MAX_IT; ++increment)
    {
        double volumeDifference = -prismVolume + increment * ((2.0 * prismVolume) / MAX_IT); 

        auto correctedPrism = correct_volume<mixedPolyhedron>(prism, polygonLabel, sweepDirection, volumeDifference); 

        auto error = mag(volume(correctedPrism) - volume(prism) - volumeDifference);  

        if (error <= EPSILON)
        {
            //write_vtk_polydata(correctedPrism, "correctPrismVolume.incrementalDirectedTest.correctedPrism", increment); 
            errorFile << volumeDifference / prismVolume << " " << error << "\n"; 
        }
    }
}


// A non convex prism where one edge pair is copplanar with the basePolygon, and the other one
// is rotating with increasing spherical angle theta. A parameterized volume correction is tested. 
TEST(correctVolume, coplanarEdgePairTrianglePrismParameterised)
{
    typedef MixedPolyhedron<arrayTriangle, pointVector, std::vector> mixedPolyhedron;

    double Ymax = 1; 
    double Xmax = 1; 
    double edgeLength = 0.5;  // Prevent self-intersection with 0.5.

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

    const unsigned int N = 100; 

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

    auto testRelativeDifference = [
        Ymax, 
        Xmax, 
        edgeLength, 
        basePolygon, 
        N,
        &outputFile
    ](double relVolDiff)
    {
        unsigned int counter = 0; 
        for (unsigned int I = 0; I < N+1; ++I)
        {
            double s = I * (edgeLength / N);

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

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

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

                auto correctVolume = (0.5 * mag(a ^ b) * Ymax); 
                auto volumeDifference = relVolDiff * correctVolume; 
                 
                auto sweptPolyhedron = build<pointVectorVector>(basePolygon, displacements);  

                auto sweptPolyhedronTri = flux_triangulate(sweptPolyhedron); 

                // Correct the swept polyhedron in the direction of the cap polygon normal. 
                auto correctedPolyhedron = correct_volume<mixedPolyhedron>(
                        sweptPolyhedron, 
                        sweptPolyhedron.size() - 1, 
                        volumeDifference
                ); 

                auto correctedVolume = volume<flux_triangulation>(correctedPolyhedron); 
                auto volumeError = mag(correctVolume + volumeDifference - correctedVolume);  

                bool passedTest = (volumeError <= EPSILON);

                // Uncomment to visualize allways. CAREFUL: Generates a lot of files.
                if (! passedTest)
                {
                    write_vtk_polydata(sweptPolyhedronTri, fullTestName() + ".sweptPolyhedronTriangulation", counter); 
                    outputFile <<  s << " " << theta << " " 
                        <<  volumeDifference / correctVolume << " " 
                        << volumeError / correctVolume << "\n";
                }

                EXPECT_TRUE(passedTest)
                    << "theta = " << theta << "\n"
                    << "s = " << s << "\n"
                    << "displacementX = " << displacementX << "\n"
                    << "displacementX = " << displacementZ << "\n"
                    << "correct volume = " << correctVolume << "\n" 
                    << "normalized flux triangulation volume error = " 
                    << volumeError / correctVolume << "\n";

                ++counter; 
            }
        }
    };

    // TODO: Implement the variation.
    testRelativeDifference(0.2); 
}

template<typename PolygonSequence, typename Point> 
void checkSideTriangleOrientation (
    const PolygonSequence& polyhedron, 
    Point const& starPoint
) 
{
    for(const auto& polygon : polyhedron)
    {
        auto triangulation = edge_triangulate(polygon);   
        for (const auto& triangle : triangulation)
        {
            auto triCentroid = centroid(triangle); 
            auto triNormal = normal_area_vec(triangle); 
            triNormal /= mag(triNormal); 
            ASSERT_GE(mag((starPoint - triCentroid) & triNormal), 0); 
        }
    }
};

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

    pointVector basePolygon = {point(0,0,0), point(0,1,0), point(1,1,0), point(1,0,0)};
    pointVector displacements = {point(0,0,1), point(0,-0.5,0.1), point(-0.2,-0.5,1), point(0.2,-0.1,1)};
     
    auto sweptPolyhedron = build<pointVectorVector>(basePolygon, displacements); 

    // Check that the edge triangulation delivers triangles oriented outwards with 
    // respect to the volume star point. Choose star points that make sense:
    // those that you know are inside the swept polyhedron and can still "see" all 
    // the points without self-intersection.
    checkSideTriangleOrientation(sweptPolyhedron, point(1,1,0));
    checkSideTriangleOrientation(sweptPolyhedron, point(0,0,0));

    auto sweptPolyTri = flux_triangulate<tetrahedronVector>(sweptPolyhedron); 
    auto sweptPolyVol = volume(sweptPolyTri); 

    auto corrPolyhedron = correct_volume<mixedPolyhedron>(
            sweptPolyhedron, 
            sweptPolyhedron.size() - 1, 
            0.1 
    ); 

    auto corrPolyTri = flux_triangulate<tetrahedronVector>(corrPolyhedron); 
    auto corrPolyVol = volume(corrPolyTri);  

    write_vtk_polydata(basePolygon, fullTestName() + ".basePolygon.vtk"); 
    write_vtk_polydata(sweptPolyhedron, fullTestName() + ".sweptPolyhedron.vtk");
    write_vtk_polydata(sweptPolyTri, fullTestName() + ".sweptPolyhedronTriangulation.vtk");
    write_vtk_polydata(corrPolyhedron, fullTestName() + ".corrPolyhedron.vtk");
    write_vtk_polydata(corrPolyTri, fullTestName() + ".corrPolyhedronTriangulation.vtk");

    ASSERT_LE(mag(corrPolyVol - sweptPolyVol - 0.1), SMALL); 
}

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

    pointVector basePolygon = {point(0,0,0), point(0,1,0), point(1,1,0), point(1,0,0)};
    pointVector displacements = {point(0,0,1), point(-0.2,0,1), point(-0.2,-0.5,0.1), point(0.2,-0.5,1)};
     
    auto sweptPolyhedron = build<pointVectorVector>(basePolygon, displacements); 
    auto sweptPolyTri = flux_triangulate<tetrahedronVector>(sweptPolyhedron); 
    auto sweptPolyVol = volume(sweptPolyTri); 

    auto corrPolyhedron = correct_volume<mixedPolyhedron>(
            sweptPolyhedron, 
            sweptPolyhedron.size() - 1, 
            0.1 
    ); 
    auto corrPolyTri = flux_triangulate<tetrahedronVector>(corrPolyhedron); 

    auto corrPolyVol = volume(corrPolyTri);  

    write_vtk_polydata(basePolygon, fullTestName() + ".basePolygon.vtk"); 
    write_vtk_polydata(sweptPolyhedron, fullTestName() + ".sweptPolyhedron.vtk");
    write_vtk_polydata(sweptPolyTri, fullTestName() + ".sweptPolyhedronTriangulation.vtk");
    write_vtk_polydata(corrPolyhedron, fullTestName() + ".corrPolyhedron.vtk");
    write_vtk_polydata(corrPolyTri, fullTestName() + ".corrPolyhedronTriangulation.vtk");

    ASSERT_LE(mag(corrPolyVol - sweptPolyVol - 0.1), SMALL); 
}

TEST(correctVolume, nonPlanarNonConvexNeighboringSide)
{
    pointVector P1 = {point(0,0,0), point(0,1,0), point(1,1,0), point(1,0,0)};
    pointVector P2 = {point(1,0,0), point(1,1,0), point(2,1,0), point(2,0,0)};

    pointVector D = {point(-0.2,0.2,0.2), point(-0.4,0.2,0.2), point(-0.4,0.2,0.2), point(-0.2,0.2,0.2)};

    auto S1 = build<pointVectorVector>(P1, D); 
    auto S2 = build<pointVectorVector>(P2, D); 

    // Make sure that the normals are oriented outwards with respect to the
    // star point used to triangulate the volume. 
    checkSideTriangleOrientation(S1, point(1,0,0)); 
    checkSideTriangleOrientation(S2, point(2,0,0)); 

    // Check if the triangulation points of two "coplanar", adjacent sides are
    // equal up to EPSILON.
    const auto S1side = S1[3];
    const auto S2side = S2[1];

    write_vtk_polydata(S1side, fullTestName() + ".S1side.vtk");
    write_vtk_polydata(S2side, fullTestName() + ".S2side.vtk");

    auto S1sideTri = edge_triangulate(S1side);  
    auto S2sideTri = edge_triangulate(S2side);  

    write_vtk_polydata(S1sideTri, fullTestName() + ".S1sideTri.vtk");
    write_vtk_polydata(S2sideTri, fullTestName() + ".S2sideTri.vtk");

    ASSERT_LE(mag(S1sideTri[0][0] - S2sideTri[0][0]), EPSILON)
        << "first triangulation point = " << S1sideTri[0][0] << "\n" 
        << "second triangulation point = " << S2sideTri[0][0] << "\n";
}

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

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

    return RUN_ALL_TESTS();

    return 0;
}

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