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

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

Description
    Volume calculation. 

Author
    Tomislav Maric maric@mma.tu-darmstadt.de

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

#ifndef Volume_HPP
#define Volume_HPP

#include "Area.hpp"

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

namespace geophase { 

    template <int N = 0, typename AABBox, typename RT=double>
    enable_if_type_tagged_t<AABBox, axis_box_tag, RT>
    volume(AABBox const& box)
    {
        RT result = box.maxPoint()[N] - box.minPoint()[N]; 

        if constexpr (N < box.maxPoint().size() - 1)
            result *= volume<N+1, AABBox, RT>(box); 

        return result;
    }

    template <typename Point, typename RT=double>
    enable_if_type_tagged_t <Point, vector_tag, RT>
    volume(
        Point const& P0, 
        Point const& P1, 
        Point const& P2, 
        Point const& P3
    )
    {
        return (1.0 / 6.0) * 
            dot(P1 - P0, cross(P2-P0,P3-P0)); 
    }

    template <typename Polyhedron, typename RT=double>
    std::enable_if_t<tag_base_of<Polyhedron, sequence_polyhedron_tag>(), RT>
    volume_by_surf_tri(Polyhedron const& polyhedron)
    {
        RT result(0); 

        for (const auto& polygon : polyhedron)
            result += dot(area_normal_vector(polygon), centroid(polygon));

        return (1.0 / 3.0) * result; 
    }


    template <typename Polyhedron, typename RT=double>
    std::enable_if_t<tag_base_of<Polyhedron, sequence_polyhedron_tag>(), RT>
    volume_by_vol_tri(Polyhedron const& polyhedron)
    {
        RT result(0); 

        const auto polyhCenter = centroid(polyhedron); 

        for (const auto& polygon : polyhedron)
        {
            const auto polygCenter = centroid(polygon); 

            auto it = polygon.begin(); 
            auto itnext = std::next(polygon.begin()); 

            const auto polyhPolyg = polygCenter - polyhCenter;  

            for(; itnext != polygon.end(); ++it, ++itnext)
            {
                result += dot(
                    polyhPolyg, 
                    cross(
                        *it - polyhCenter, 
                        *itnext - polyhCenter
                    )
                );
            }
            result += dot(
                polyhPolyg, 
                cross(
                    *it - polyhCenter, 
                    polygon.front() - polyhCenter
                )
            );
        }

        result *= 1.0 / 6.0;

        return result;
    }

    template <typename Polyhedra, typename RT=double>
    std::enable_if_t<
        type_tagged<Polyhedra, sequence_tag>() &&
        tag_base_of<typename Polyhedra::value_type, sequence_polyhedron_tag>(), 
        RT 
    >
    volume_by_vol_tri(Polyhedra const& polyhedra)
    {
        RT volume {0}; 

        for (const auto& polyhedron : polyhedra)
            volume += volume_by_vol_tri(polyhedron);

        return volume; 
    }

    template <typename Polyhedra, typename RT=double>
    std::enable_if_t<
        type_tagged<Polyhedra, sequence_tag>() &&
        tag_base_of<typename Polyhedra::value_type, sequence_polyhedron_tag>(), 
        RT 
    >
    volume_by_surf_tri(Polyhedra const& polyhedra)
    {
        RT volume {0}; 

        for (const auto& polyhedron : polyhedra)
            volume += volume_by_surf_tri(polyhedron);

        return volume; 
    }

    template <typename Polyhedra, typename RT=double>
    std::enable_if_t<
        type_tagged<Polyhedra, sequence_tag>() &&
        tag_base_of<typename Polyhedra::value_type, sequence_polyhedron_tag>(), 
        std::vector<RT> 
    >
    volumes_by_vol_tri(Polyhedra const& polyhedra)
    {
        std::vector<RT> result {}; 
        result.reserve(polyhedra.size()); 

        for (const auto& polyhedron : polyhedra)
            result.push_back(volume_by_vol_tri(polyhedron));

        return result; 
    }

    template <typename Polyhedra, typename RT=double>
    std::enable_if_t<
        type_tagged<Polyhedra, sequence_tag>() &&
        tag_base_of<typename Polyhedra::value_type, sequence_polyhedron_tag>(), 
        std::vector<RT> 
    >
    volumes_by_surf_tri(Polyhedra const& polyhedra)
    {
        std::vector<RT> result {}; 
        result.reserve(polyhedra.size()); 

        for (const auto& polyhedron : polyhedra)
            result.push_back(volume_by_surf_tri(polyhedron));

        return result; 
    }

} // End namespace geophase 

#endif

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