// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_DORIE_TRANSPORT_OPERATOR_HH
#define DUNE_DORIE_TRANSPORT_OPERATOR_HH

#include <vector>

#include <dune/geometry/referenceelements.hh>

#include <dune/pdelab/common/quadraturerules.hh>
#include <dune/pdelab/localoperator/pattern.hh>
#include <dune/pdelab/localoperator/flags.hh>
#include <dune/pdelab/localoperator/idefault.hh>
#include <dune/pdelab/localoperator/defaultimp.hh>

#include <dune/dorie/common/typedefs.hh>

namespace Dune {
namespace Dorie {
namespace Operator {

/*-------------------------------------------------------------------------*//**
 * @brief      Spatial local operator for the transport equation in unsaturated
 *             media in a finite volume scheme.
 * @details    It solves the spatial part of the transport equation:
 *
 *             @f{eqnarray*}{
 *             \partial_t[\theta C_w] +
 *             \nabla\cdot [\textbf{j}_w C_w] -
 *             \nabla [\theta \mathsf{D}_{\text{eff}}\nabla C_w]&=0 &\qquad \text{in }
 *             \Omega\\
 *             C_w &= g &\qquad \text{on } \Gamma_D
 *             \subseteq\partial\Omega\\
 *             \nabla C_w \cdot \textbf{n} &= \textbf{j}_{\scriptscriptstyle
 *             C_w}& \qquad \text{on } \Gamma_N =\partial\Omega \backslash
 *             \Gamma_D
 *             @f}
 *
 * @author     Santiago Ospina De Los Ríos
 * @date       2018
 * @ingroup    LocalOperators
 * @ingroup    TransportModel
 *
 * @todo Use diffusion coefficient from parameterization
 *
 * @tparam     Boundary        Type of the class providing boundary conditions
 * @tparam     GFWaterFlux     Type of a grid function which provides the water
 *                             flux
 * @tparam     GFWaterContent  Type of a grid function which provides the water
 *                             content
 */
template<class Boundary, class GFWaterFlux, class GFWaterContent>
class TransportFVSpatialOperator
    : public Dune::PDELab::NumericalJacobianVolume<
                        TransportFVSpatialOperator<
                            Boundary,
                            GFWaterFlux,
                            GFWaterContent
                        >
                    > 
    , public Dune::PDELab::NumericalJacobianSkeleton<
                        TransportFVSpatialOperator<
                            Boundary,
                            GFWaterFlux,
                            GFWaterContent
                        >
                    >
    , public Dune::PDELab::NumericalJacobianBoundary<
                        TransportFVSpatialOperator<
                            Boundary,
                            GFWaterFlux,
                            GFWaterContent
                        >
                    >
    , public Dune::PDELab::FullSkeletonPattern
    , public Dune::PDELab::FullVolumePattern
    , public Dune::PDELab::LocalOperatorDefaultFlags
    , public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>
{
public:
    // Pattern assembly flags
    enum { doPatternVolume = true };
    enum { doPatternSkeleton = true };

    // Residual assembly flags
    enum { doAlphaBoundary = true };
    enum { doAlphaSkeleton  = true };

private:
    static_assert(std::is_same<
            typename GFWaterFlux::Traits::GridViewType,
            typename GFWaterContent::Traits::GridViewType>::value, 
        "TransportFVSpatialOperator: GFWaterFlux and"
        "GFWaterContent has to use the same grid view.");

    // Define range type for water flux and water content
    using WaterFlux = typename GFWaterFlux::Traits::RangeType;
    using WaterContent = typename GFWaterContent::Traits::RangeType;

public:

    /**
     * @brief      Constructor of TransportFVSpatialOperator
     *
     * @param      boundary          The boundary terms
     * @param[in]  diff_coeff        The diffusion coefficient
     */
    TransportFVSpatialOperator(
        std::shared_ptr<Boundary> boundary, 
        double diff_coeff
    )   : _boundary(boundary)
        , _time(0.)
        , _diff_coeff(diff_coeff)
    { }

    /*---------------------------------------------------------------------*//**
     * @brief      Skeleton integral depending on test and ansatz functions.
     *             Each face is only visited once since this method is symmetric
     *
     * @param[in]  ig      The intersection entity of the grid (inside + outside
     *                     entities)
     * @param[in]  lfsu_i  The inside ansatz local function space
     * @param[in]  x_i     The coefficients of the lfsu_i
     * @param[in]  lfsv_i  The inside test local function space
     * @param[in]  lfsu_o  The outside ansatz local function space
     * @param[in]  x_o     The coefficients of the lfsu_o
     * @param[in]  lfsv_o  The outside test local function space
     * @param      r_i     The view of the residual vector w.r.t lfsu_i
     * @param      r_o     The view of the residual vector w.r.t lfsu_o
     *
     * @tparam     IG      The type for ig
     * @tparam     LFSU    The type for lfsu_i and lfsu_o
     * @tparam     X       The type for x_i and x_o
     * @tparam     LFSV    The type for lfsv_i and lfsv_o
     * @tparam     R       The type for r_i and r_o
     */
    template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
    void alpha_skeleton(const IG& ig, 
        const LFSU& lfsu_i, const X& x_i, const LFSV& lfsv_i,
        const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o, 
        R& r_i, R& r_o) const
    {
        // Get entity dim
        // const int dim = IG::Entity::dimension;

        // Get local basis traits from local function space
        using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
                                         Traits::LocalBasisType::Traits;

        // Get range, range field and gradient for trial space
        using RangeU = typename LocalBasisTraitsU::RangeType;
        // using RangeFieldU = typename LocalBasisTraitsU::RangeFieldType;
        // using GradientU = Dune::FieldVector<RangeFieldU,dim>;

        // Assert that polynomial degree is always 0
        assert(lfsu_i.finiteElement().localBasis().order()==0);
        assert(lfsu_o.finiteElement().localBasis().order()==0);
        assert(lfsv_i.finiteElement().localBasis().order()==0);
        assert(lfsv_o.finiteElement().localBasis().order()==0);

        // Define entities for each frame of reference
        const auto& entity_f = ig.intersection();
        const auto& entity_i = ig.inside();
        const auto& entity_o = ig.outside();

        // Get geometries
        auto geo_f = entity_f.geometry();
        auto geo_i = entity_i.geometry();
        auto geo_o = entity_o.geometry();

        // Get geometry of intersection in local coord. of neighbor entities
        auto geo_in_i = entity_f.geometryInInside();

        // Get volume of entity
        const auto volume_f = geo_f.volume();

        // Get normal
        auto normal_f = entity_f.centerUnitOuterNormal();

        // Entity centers in references elements
        const auto& ref_el_i = referenceElement(geo_i);
        const auto& ref_el_o = referenceElement(geo_o);
        const auto& center_position_i  = ref_el_i.position(0,0);
        const auto& center_position_o = ref_el_o.position(0,0);
        
        // Get center position of the face w.r.t inside entity 
        auto face_position_i = geo_in_i.center();

        // Evaluate water content
        WaterContent water_content_i,water_content_o;
        _gf_water_content->evaluate(entity_i,center_position_i,water_content_i);
        _gf_water_content->evaluate(entity_o,center_position_o,water_content_o);

        // Eval flux field assume H(div) velocity field => may choose any side
        WaterFlux water_flux;
        _gf_water_flux->evaluate(entity_i, face_position_i, water_flux);
        auto water_flux_n = water_flux*normal_f;
        // auto water_vel_i = water_flux_n*water_content_i;
        // auto water_vel_o = water_flux_n*water_content_o;

        // Compute the effective hydrodynamic dispersion

        // bind parameterization and retrieve functions
        // _param->bind(entity_i);
        // const auto hydr_disp_func_i = _param->hydrodynamic_dispersion_f();
        // auto hydr_disp_i = hydr_disp_func_i(water_vel_i,water_content_i);

        // _param->bind(entity_o);
        // const auto hydr_disp_func_o = _param->hydrodynamic_dispersion_f();
        // auto hydr_disp_o = hydr_disp_func_o(water_vel_o,water_content_o);

        // Evaluate diffusion coeff. from both sides
        // GradientU diff_coeff_i, diff_coeff_o;
        // hydr_disp_i.mv(normal_f,diff_coeff_i);
        // hydr_disp_o.mv(normal_f,diff_coeff_o);

        auto diff_coeff_i = _diff_coeff * water_content_i; // FIXME
        auto diff_coeff_o = _diff_coeff * water_content_o; // FIXME
        
        // Take harmonic average
        auto diff_coeff_f = 2.0/(   1.0/(diff_coeff_i + 1E-30) 
                                  + 1.0/(diff_coeff_o + 1E-30) );

        // Inside/outside solute value
        RangeU u_i = x_i(lfsu_i,0);
        RangeU u_o = x_o(lfsu_o,0);
        
        // Upwinding
        const auto& u_upwind = (water_flux_n>=0) ? u_i : u_o;

        // Entity centers in global coordinates
        auto center_position_i_g = geo_i.global(center_position_i);
        auto center_position_o_g = geo_o.global(center_position_o);

        // Distance between the two entity centers
        center_position_i_g -= center_position_o_g;
        auto distance = center_position_i_g.two_norm();

        // Finite difference of u between the two entities
        RangeU dudn = (u_o-u_i)/distance;

        // Solute flux in normal direction w.r.t the intersection
        auto soulte_flux_n = u_upwind*water_flux_n - diff_coeff_f*dudn;

        // Symmetric contribution to residual on inside element
        r_i.accumulate(lfsv_i, 0,  soulte_flux_n*volume_f );
        r_o.accumulate(lfsv_o, 0, -soulte_flux_n*volume_f );
    }


    /*---------------------------------------------------------------------*//**
     * @brief      Boundary integral depending on test and ansatz functions.
     *
     * @param[in]  ig      The intersection entity of the grid (inside + outside
     *                     entities)
     * @param[in]  lfsu_i  The inside ansatz local function space
     * @param[in]  x_i     The coefficients of the lfsu_i
     * @param[in]  lfsv_i  The inside test local function space
     * @param      r_i     The view of the residual vector w.r.t lfsu_i
     *
     * @tparam     IG      The type for ig
     * @tparam     LFSU    The type for lfsu_i
     * @tparam     X       The type for x_i
     * @tparam     LFSV    The type for lfsv_i
     * @tparam     R       The type for r_i
     */
    template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
    void alpha_boundary (const IG& ig, const LFSU& lfsu_i, const X& x_i,
                                       const LFSV& lfsv_i, R& r_i) const
    {
        // Get entity dim
        // const int dim = IG::Entity::dimension;

        // Get local basis traits from local function space
        using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
                            Traits::LocalBasisType::Traits;

        // Get range, range field and gradient for trial space
        using RangeU = typename LocalBasisTraitsU::RangeType;
        // using RangeFieldU = typename LocalBasisTraitsU::RangeFieldType;
        // using GradientU = Dune::FieldVector<RangeFieldU,dim>;

        // Assert that polynomial degree is always 0
        assert(lfsu_i.finiteElement().localBasis().order()==0);
        assert(lfsv_i.finiteElement().localBasis().order()==0);

        // Get entities
        const auto& entity_f = ig.intersection(); 
        const auto& entity_i = ig.inside();

        // Get geometries
        auto geo_f = entity_f.geometry();
        const auto& ref_el_f = referenceElement(geo_f);
        const auto& center_position_f = ref_el_f.position(0,0);

        // Get normal 
        auto normal_f = entity_f.centerUnitOuterNormal();

        // Get boundary condition type
        auto bc = _boundary->bc(entity_f,center_position_f,_time);

        // Get geometry of intersection in local coord. of neighbor entity
        auto geo_in_i = ig.geometryInInside();
        auto face_position_i = geo_in_i.center();

        // Evaluate convective term (i.e water flux)
        WaterFlux water_flux;
        _gf_water_flux->evaluate(entity_i, face_position_i, water_flux);
        auto water_flux_n = water_flux*normal_f;
        // auto water_vel_i = water_flux_n*water_content_i;

        // bind parameterization and retrieve functions
        // _param->bind(entity_i);
        // const auto hydr_disp_func_i = _param->hydrodynamic_dispersion_f();

        // Face volume for integration
        auto volume_f = geo_f.volume();

        if (BoundaryCondition::isOutflow(bc)) 
        {
            // Check that water is going out
            if (Dune::FloatCmp::lt(water_flux_n,0.))
                return;

            // Inside unknown value
            RangeU u = x_i(lfsu_i,0);

            // Eval flux boundary condition
            auto o = _boundary->o(entity_f,center_position_f,_time);

            // convection term
            auto conv_num_flux = (u * water_flux_n + o) * volume_f;

            // integrate o
            r_i.accumulate(lfsu_i,0,conv_num_flux);
        }
        else if (BoundaryCondition::isDirichlet(bc))
        {
            // Only allow dirichlet values for influx cases
            if (Dune::FloatCmp::gt(water_flux_n,0.))
                return;

            // Get center of entity in global coordinates
            auto geo_i = entity_i.geometry();
            auto center_position_i_g = geo_i.center();

            // Face center in global coordinates
            auto center_position_f_g = geo_f.center();

            // Compute distance of these two points
            center_position_i_g -= center_position_f_g;
            auto distance = center_position_i_g.two_norm();

            // Get geometry of intersection in local coord. of neighbor entity
            auto geo_in_i = ig.geometryInInside();
            auto face_position_i = geo_in_i.center();

            // Evaluate water content
            WaterContent water_content_i;
            _gf_water_content->evaluate(entity_i, face_position_i, water_content_i);

            // Evaluate Dirichlet condition
            auto g = _boundary->g(entity_f,center_position_f,_time);

            // Compute the effective hydrodynamic dispersion
            // auto hydr_disp_i = hydr_disp_func_i(water_vel_i,water_content_i);

            // Calculate diffusion coefficient from inside
            // GradientU diff_coeff_i;
            // hydr_disp_i.mv(normal_f,diff_coeff_i);

            auto diff_coeff_i = _diff_coeff; // FIXME

            // Inside unknown value
            RangeU u = x_i(lfsu_i,0);

            // Solute flux in normal direction w.r.t the intersection
            auto soulte_flux_n = (g*water_flux_n)
                                -(diff_coeff_i*water_content_i*(g-u)/distance);

            // Contribution to residual from Dirichlet boundary
            r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f);

        }
        else if (BoundaryCondition::isNeumann(bc))
        {
            // Solute flux in normal direction w.r.t the intersection
            auto soulte_flux_n = _boundary->j(entity_f,center_position_f,_time);

            r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f );
        }
    }

    /*---------------------------------------------------------------------*//**
     * @brief      Sets the time.
     *
     * @param[in]  t           time of type RangeField
     *
     * @tparam     RangeField  type of the range field
     */
    template<class RangeField>
    void setTime (RangeField t)
    {
        _time = t;
        _gf_water_content->setTime(t);
        _gf_water_flux->setTime(t);
    }

    /*---------------------------------------------------------------------*//**
     * @brief      Sets the the water content grid function.
     *
     * @param[in]  gf_water_content  The water content grid function
     */
    void set_water_content(std::shared_ptr<GFWaterContent> gf_water_content)
    {
        _gf_water_content = gf_water_content;
    }

    /*---------------------------------------------------------------------*//**
     * @brief      Sets the the water flux grid function.
     *
     * @param[in]  gf_water_content  The water flux grid function
     */
    void set_water_flux(std::shared_ptr<GFWaterFlux> gf_water_flux)
    {
        _gf_water_flux = gf_water_flux;
    }

private:
    const std::shared_ptr<const Boundary> _boundary;
    std::shared_ptr<GFWaterFlux> _gf_water_flux;
    std::shared_ptr<GFWaterContent> _gf_water_content;
    double _time;
    double _diff_coeff;
};

/**
 * @brief      Temporal local operator for the transport equation in unsaturated
 *             media in a finite volume scheme.
 * @details    It solves the temporal part of the transport equation:
 *
 *             @f{eqnarray*}{
 *             \partial_t[\theta C_w] +
 *             \nabla\cdot [\textbf{j}_w C_w] -
 *             \nabla [\theta \mathsf{D}_{\text{eff}}\nabla C_w]&=0 &\qquad \text{in }
 *             \Omega\\
 *             C_w &= g &\qquad \text{on } \Gamma_D
 *             \subseteq\partial\Omega\\
 *             \nabla C_w \cdot \textbf{n} &= \textbf{j}_{\scriptscriptstyle
 *             C_w}& \qquad \text{on } \Gamma_N =\partial\Omega \backslash
 *             \Gamma_D
 *             @f}
 *
 * @author     Santiago Ospina De Los Ríos
 * @date       2018
 * @ingroup    LocalOperators
 * @ingroup    TransportModel
 *
 * @bug        The water content is not being updated by the set time to the
 *             right state of the grid function.
 *
 * @tparam     GFWaterContent  Type of a grid function which provides the water
 *                             content
 */
template<class GFWaterContent>
  class TransportFVTemporalOperator
  : public Dune::PDELab::NumericalJacobianVolume<
                            TransportFVTemporalOperator<GFWaterContent>>
  , public Dune::PDELab::NumericalJacobianApplyVolume<
                            TransportFVTemporalOperator<GFWaterContent>>
  , public Dune::PDELab::FullSkeletonPattern
  , public Dune::PDELab::FullVolumePattern
  , public Dune::PDELab::LocalOperatorDefaultFlags
  , public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>
{
public:

    // Pattern assembly flags
    enum { doPatternVolume = true };

    // Eesidual assembly flags
    enum { doAlphaVolume = true };

public:

    TransportFVTemporalOperator()
    : _time(0.)
    { }

    /*---------------------------------------------------------------------*//**
     * @brief      Volume integral depending on test and ansatz functions
     *
     * @param[in]  eg    THe entity of the grid
     * @param[in]  lfsu  The ansatz local function space
     * @param[in]  x     The coefficients of the lfsu
     * @param[in]  lfsv  The test local function space
     * @param      r     The view of the residual vector w.r.t lfsu
     *
     * @tparam     EG    The type for eg
     * @tparam     LFSU  The type for lfsu
     * @tparam     X     The type for x
     * @tparam     LFSV  The type for lfsv
     * @tparam     R     The type for r
     */
    template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
    void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, 
                                     const LFSV& lfsv, R& r) const
    {
        // Get local basis traits from local function space
        using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
                                                Traits::LocalBasisType::Traits;
        // Get range for trial space
        using RangeU = typename LocalBasisTraitsU::RangeType;

        RangeU u = x(lfsu,0);

        // entity geometry
        auto geo = eg.geometry();

        // inside cell center
        const auto& ref_el = referenceElement(geo);
        const auto& center_position = ref_el.position(0,0);

        typename GFWaterContent::Traits::RangeType water_content;

        _gf_water_content->evaluate(eg.entity(),center_position,water_content);

        // update residual
        r.accumulate(lfsv ,0 , water_content*u*geo.volume());
    }

    /*---------------------------------------------------------------------*//**
     * @brief      Sets the time.
     *
     * @param[in]  t           time of type RangeField
     *
     * @tparam     RangeField  type of the range field
     */
    template<class RF>
    void setTime (RF t)
    {
        _time = t;
        _gf_water_content->setTime(t);
    }

    /*---------------------------------------------------------------------*//**
     * @brief      Sets the the water content grid function.
     *
     * @param[in]  gf_water_content  The water content grid function
     */
    void set_water_content(std::shared_ptr<GFWaterContent> gf_water_content)
    {
        _gf_water_content = gf_water_content;
    }

private:
    std::shared_ptr<GFWaterContent> _gf_water_content;
    double _time;
};

} // namespace Operator
} // namespace Dorie
} // namespace Dune


#endif // DUNE_DORIE_TRANSPORT_OPERATOR_HH