#ifndef _VOLTAGE_H_
#define _VOLTAGE_H_

#include "Problem.h"

using namespace dealii;

/**
 * The boundary condition for the voltage. It is just a child of the Function class template.
 */
template<int dim> class Vbc: public Function<dim>{
  public:
    /**
     * Set the test case and the initial time.
     * 
     * @param testcase : The test case. We need different voltage patterns to get different
     * effects.
     * @param t : The initial time.
     * Because of the way the solution scheme is written. The initial data is not
     * compatible. For this reason, we will for a little time advance without turning
     * on the voltage. How long is controlled by this parameter.
     */
    Vbc( const TestCase testcase, const double t = 0. );
    /**
     * Does nothing
     */
    virtual ~Vbc();
    /**
     * The value. As usual:
     * @param p : The point at which the function is evaluated.
     * @param component : Bogus parameter inherited from Function
     */
    virtual double value( const Point<dim> &p, unsigned component = 0 ) const;
  protected:
    /**
     * The testcase
     */
    const TestCase tc;
};

/**
 * An elliptic equation that models the electric potential:
 * \f[
 *   - \nabla\cdot\left( [ \varepsilon(\phi) ] \nabla V \right) = q \chi_\Omega,
 *    \quad \text{ in } \Omega^\star,
 * \f]
 * with boundary conditions
 * \f[
 *   V = V_0, \quad \text{ on }\partial^\star\Omega^\star,
 *  \qquad
 *  \partial_{\mathbf{n}} V = 0, \text{ on }\partial\Omega^\star \cap \Gamma,
 * \f]
 * where
 * \f[
 *  \varepsilon(\phi) = \begin{cases}
 *     \varepsilon(\phi), & \Omega, \\
 *     \varepsilon^\star, & \Omega^\star \setminus \Omega,
 *   \end{cases}
 * \f]
 * with $\varepsilon^\star$ a constant that denotes the permittivity at the plates.
 */
template<int dim> class Voltage: 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 space where the 
     *    voltage 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.
     * @param testcase : The testcase. Used to construct a suitable boundary condition.
     */
    Voltage( const Triangulation<dim> &tria, Material_Parameters &params, const unsigned deg,
             const unsigned u_prec, const double ee, const double thres, const unsigned sweeps,
             const TestCase testcase );
    /**
     * Clear the data
     */
    virtual ~Voltage();
  protected:
    const material_function &epsilon; ///< Permitivity
    const double epsilon_plate; ///< Permittivity of the plates
    Vbc<dim> volt; ///< boundary conditions
    typename FunctionMap<dim>::type BCs; ///< boundary indicators
    TrilinosWrappers::PreconditionAMG prec; ///< Preconditioner (AMG)
    TrilinosWrappers::PreconditionAMG::AdditionalData prec_data; ///< Preconditioner Data
    /**
     * We group the iterators of this problem and the problems this one 
     * takes information from as follows:
     * <ul>
     *    <li> Iterator[0]: <b>this</b>. </li>
     *    <li> Iterator[1]:  Charge. </li>
     *    <li> Iterator[2]: Phase Field, from CahnHilliard. </li> 
     * </ul>
     */
    typedef IteratorGroup<3, typename hp::DoFHandler<dim>::active_cell_iterator> Iterator;
    /**
     * Scratch Data for assembly of the matrix and right hand side.
     *
     * All the std::vector's are used to store local values. This is scratch data in the 
     * sense of the Parallell module of the deal.II 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> charge, ///< The AsFunction object that will be interpreted as Charge
                      phase; ///< Phase field, part of CahnHilliard. 
      std::vector<double> loc_charge, loc_phase;
      ScratchData( const hp::FECollection<dim> &fe, const std::vector<AsFunction<dim> *> &data,
                   const hp::QCollection<dim> &q, const UpdateFlags u_flags ) :
          quad(q), flags( u_flags ), fe_val( fe, quad, flags ),
          charge( *( data[0] ) ), phase( *( data[1] ) ),
          loc_charge( quad.max_n_quadrature_points() ), loc_phase( 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 ),
          charge( scratch.charge ), phase( scratch.phase ),
          loc_charge( quad.max_n_quadrature_points() ), loc_phase( 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> Charge. </li>
     *  <li> Phase Field, from CahnHilliard. </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 );
   /**
    * Apply boundary conditions.
    */
    virtual void SetupDoFsSuffix();
    /**
     * Reinitialize 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 value of the solution
     */
    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:
    /**
     * The voltage
     */
    AsFunction<dim> get_voltage;
};

#endif //_VOLTAGE_H_
