#ifndef _CHARGE_H_
#define _CHARGE_H_

#include "Problem.h"

/**
 * A convection diffusion equation to model the transport of charges:
 * \f[
 *    q_t + \nabla\cdot (q \mathbf{u} ) = 
 *   \nabla \cdot \left[ K(\phi) \nabla 
 *          \left( \lambda  q + V \right)
 *          \right], \quad \text{ in } \Omega
 * \f]
 * with boundary condition
 * \f[
 *    K(\phi) \nabla\left( \lambda q + V \right) \cdot \mathbf{n} = 0,
 *    \quad \text{ on } \Gamma.
 * \f]
 */
template<int dim> class Charge: public Problem<dim>{
  public:
    /**
     * Copy the provided parameters to the members
     * @param tria : A reference to the mesh where the problem will live.
     * @param params : The collection of material parameters.
     * @param deg : The polynomial degree for the finite element spaces where the charge 
     *   is defined.
     * @param u_prec : How often we update the preconditioner.
     * @param ee : The accuracy to which we iteratively solve the problem.
     * @param thres : The aggregation threshold for the AMG preconditioner.
     * @param sweeps : The number of smoothing steps the AMG preconditioner must perform.
     */
    Charge( const Triangulation<dim> &tria, Material_Parameters &params, const unsigned deg,
            const unsigned u_prec, const double ee, const double thres, const unsigned sweeps );
    /**
     * Clear the data.
     */
    virtual ~Charge();
  protected:
    const material_function &conductivity; ///< Conductivity
    const double lambda; ///< Coupling parameter
    TrilinosWrappers::PreconditionAMG::AdditionalData prec_data; ///< Preconditioner (AMG)
    TrilinosWrappers::PreconditionAMG prec; ///< Preconditioner data.
    /**
     * We group the iterators of this problem and the problems this one takes information from
     * as follows:
     * <ul>
     *   <li> <code>Iterator[0]</code>: <b>this</b>.</li>
     *   <li> <code>Iterator[1]</code>: Voltage.</li>
     *   <li> <code>Iterator[2]</code>: Phase Field, from CahnHilliard.</li>
     *   <li> <code>Iterator[3]</code>: Velocity</li>
     * </ul>
     */
    typedef IteratorGroup<4,typename hp::DoFHandler<dim>::active_cell_iterator> Iterator;
    /**
     * Scratch Data for assembly of the matrix and right hand side.
     *
     * All the <code>std::vector</code>'s are used to store local values.
     * This is scratch data in the sense of the Parallell module of the <code>deal.II</code>
     * library, so that its only functionality is to encapsulate all the scratch data
     * and scratch vector and provide an explicit constructor and a copy constructor.
     */
    struct ScratchData{
      hp::QCollection<dim> quad; ///< The collection of quadrature formulas.
      UpdateFlags flags; ///< The update flags
      hp::FEValues<dim> fe_val; ///< The FEValues object
      AsFunction<dim> voltage, ///< The AsFunction object that will be interpreted as Voltage.
                      phase, ///< Phase field, part of CahnHilliard.
                      velocity; ///< The AsFunction object that will be interpreted as Velocity.
      std::vector<double> loc_charge, loc_phase, loc_div_vel;
      std::vector< Tensor<1, dim> > loc_vel, loc_grad_volt;
      ScratchData( const hp::FECollection<dim> &fe, std::vector<AsFunction<dim> *> &data,
                   const hp::QCollection<dim> &q, const UpdateFlags u_flags ) :
          quad(q), flags( u_flags ), fe_val( fe, quad, flags ),
          voltage( *( data[0] ) ), phase( *( data[1] ) ), velocity( *( data[2] ) ),
          loc_charge( quad.max_n_quadrature_points() ), loc_phase( quad.max_n_quadrature_points() ),
          loc_div_vel( quad.max_n_quadrature_points() ), loc_vel( quad.max_n_quadrature_points() ),
          loc_grad_volt( quad.max_n_quadrature_points() ) {}
      ScratchData( const ScratchData &scratch ) :
          quad( scratch.quad ), flags( scratch.flags ), fe_val( scratch.fe_val.get_fe_collection(), quad, flags ),
          voltage( scratch.voltage ), phase( scratch.phase ), velocity( scratch.velocity ),
          loc_charge( quad.max_n_quadrature_points() ), loc_phase( quad.max_n_quadrature_points() ),
          loc_div_vel( quad.max_n_quadrature_points() ), loc_vel( quad.max_n_quadrature_points() ),
          loc_grad_volt( quad.max_n_quadrature_points() )  {}
    };
    /**
     * Assemble the system matrix and right hand side for the current status of the problem
     * @param data : The views of the other problems this one needs information from to be
     *  able to assemble.
     * <ol>
       <li> Voltage. </li>
       <li> Phase field from CahnHilliard. </li>
       <li> Velocity. </li>
     * </ol>
     * @param time : The current time.
     * @param m_threaded : flag that indicates if the assembly should be done using multithreading.
     *
     */
    virtual void AssembleSystem( std::vector< AsFunction<dim> *> &data, const double time,
                                 const bool m_threaded = true );
    /**
     * Assemble the local problems.
     * @param Its : The cell iterators for the current problem and the data.
     * @param scratch : The scratch data.
     * @param data : The per task data.
     */
    void AssembleCell( const Iterator Its, ScratchData &scratch,
                       typename Problem<dim>::PerTaskData &data );
    /**
     * Reinit the preconditioner
     */
    virtual void ReinitPrec();
    /**
     * Solve the system.
     * @param control : Provides the maximal number of iterations and exit criterion.
     */
    virtual void DoSolve( SolverControl &control );
    /**
     * Set the initial data
     */
    virtual void SetInitialData();
    /**
     * Used to handle changes in the triangulation:
     * Copy to the solution transfer
     */
    virtual void PreRefinementPreffix();
    /**
     * Used to handle changes in the triangulation:
     * Copy from the solution transfer
     */
    virtual void PostRefinementSuffix( const std::vector<TrilinosWrappers::Vector> &sol_tmp );
  public:
    AsFunction<dim> get_charge; ///< The charge
};

#endif //_CHARGE_H_
