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

    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
    Generic builders for different geometries from OpenFOAM mesh elements. 

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

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

#ifndef BuildFoam_H 
#define BuildFoam_H 

#include <vector>
#include <list>
#include <numeric>

#include "tag.H"
#include "pointMesh.H"
#include "AABBox.H"
#include "Build.H"
#include "surfaceFields.H"

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

namespace Foam { namespace GeometricalTransport { 

    template <typename AABBox> 
    typename std::enable_if
    <
        Traits::tag_enabled<AABBox, aabbox_tag>::value,
        AABBox 
    >::type 
    build(const cell& cell, polyMesh const & mesh)
    {
        const auto& meshFaces = mesh.faces(); 
        const auto& meshPoints = mesh.points(); 
        
        // Not creating a local copy of cell points for efficiency reasons. TM.
        std::pair<point,point> extremeCellPoints = std::make_pair(
            point(GREAT, GREAT, GREAT),
            point(-GREAT, -GREAT, -GREAT)
        );

        forAll(cell, faceI)
        {
            const face& meshFace = meshFaces[cell[faceI]]; 
            forAll(meshFace, pointI)
            {
                less_assign(meshPoints[meshFace[pointI]], extremeCellPoints.first); 
                greater_assign(meshPoints[meshFace[pointI]], extremeCellPoints.second);
            }
        }

        return AABBox(extremeCellPoints.first, extremeCellPoints.second);
    }

    template <typename PointSequence> 
    typename std::enable_if
    <
        Traits::tag_enabled<PointSequence, point_sequence_tag>::value,
        PointSequence 
    >::type 
    build(face const & face, polyMesh const & mesh)
    {
        PointSequence polygon;
        const auto& points = mesh.points(); 

        for (auto pointLabel : face)
        {
            polygon.push_back(points[pointLabel]); 
        }

        return polygon;
    }

    template <typename PolygonSequence>
    typename std::enable_if
    <
        Traits::tag_enabled<PolygonSequence, polygon_sequence_tag>::value, 
        PolygonSequence
    >::type 
    build(const face& meshFace, const polyMesh& mesh, pointField const& displacements)
    {
        auto basePolygon = 
            build<typename PolygonSequence::value_type>(meshFace, mesh); 

        return build<PolygonSequence>(basePolygon, displacements);
    }

    template <typename PointSequence> 
    typename std::enable_if
    <
        Traits::tag_enabled<PointSequence, point_sequence_tag>::value,
        PointSequence 
    >::type 
    build(label const & faceLabel, const polyMesh & mesh)
    {
        const auto& faces = mesh.faces(); 
        const auto& face = faces[faceLabel];

        return build<PointSequence>(face, mesh);
    }

    template <typename PolygonSequence> 
    typename std::enable_if
    <
        Traits::tag_enabled<PolygonSequence, polygon_sequence_tag>::value,
        PolygonSequence
    >::type 
    build(label cellLabel, const polyMesh & mesh)
    {
        typedef typename PolygonSequence::value_type PolygonType; 
        PolygonSequence polyhedron; 
        const auto& own = mesh.faceOwner(); 
        const auto& cell = mesh.cells()[cellLabel]; 

        for (auto faceLabel : cell)
        {
            auto polygon = build<PolygonType>(faceLabel, mesh);

            // If the face normal looks into the cell. 
            if (own[faceLabel] != cellLabel)
            {
                std::reverse(polygon.begin(), polygon.end()); 
            }

            // Append the poylygon to the PolygonCollection. 
            polyhedron.push_back(polygon); 
        }

        return polyhedron;
    }

    template <typename Halfspace>  
    typename std::enable_if
    <
        Traits::tag_enabled<Halfspace, halfspace_tag>::value,
        Halfspace 
    >::type 
    build(const face& meshFace, const pointField& meshPoints)
    {
        const point& P0 = meshPoints[meshFace[0]];  
        const point& P1 = meshPoints[meshFace[1]];  
        const point& P2 = meshPoints[meshFace[2]];  

        auto n = (P1 - P0) ^ (P2 - P0);

        return Halfspace(P0, n); 
    }
   
    template <typename Halfspace>  
    typename std::enable_if
    <
        Traits::tag_enabled<Halfspace, halfspace_tag>::value,
        Halfspace 
    >::type 
    build(const label& faceLabel, const fvMesh& mesh)
    {
        const auto& Sf = mesh.Sf(); 
        const auto& Cf = mesh.Cf(); 

        return Halfspace(Cf[faceLabel], Sf[faceLabel]); 
    }

    // Generate a sequence of outward oriented halfspaces from a mesh cell.
    template <typename HalfspaceCollection> 
    typename std::enable_if
    <
        Traits::tag_enabled<HalfspaceCollection, halfspace_collection_tag>::value,
        HalfspaceCollection 
    >::type 
    build(label cellLabel, const fvMesh& mesh)
    {
        HalfspaceCollection halfspaces;

        typedef typename HalfspaceCollection::value_type Halfspace; 

        const auto& Sf = mesh.Sf(); 
        const auto& Cf = mesh.Cf(); 

        const auto& own = mesh.faceOwner(); 
        const auto& cell = mesh.cells()[cellLabel]; 
        const auto& faces = mesh.faces(); 
        const auto& points = mesh.points(); 

        forAll(cell, faceI)
        {
            auto faceLabel = cell[faceI]; 

            if (mesh.isInternalFace(faceLabel))
            {
                halfspaces.push_back(Halfspace(Cf[cell[faceI]], Sf[cell[faceI]]));            
                // Flip the face normal if the cell owns the face.
                // This is to ensure that the set of halfspaces built from a cell
                // defines a convex set, since it is used as an intersection tool.
                // Otherwise a flip() of each halfspace is required. TM.
                if (cellLabel == own[cell[faceI]])
                    halfspaces.back().flip(); 
            }
            else
            {
                halfspaces.push_back(Halfspace(-1*faces[faceI].normal(points), points[faces[faceI][0]]));            
            }
        }

        return halfspaces; 
    }


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

}} // End namespace Foam::GeometricalTransport

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

#endif

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