#ifndef _CAHN_HILLIARD_H_
#define _CAHN_HILLIARD_H_

#include "Problem.h"

/**
 * The initial conditiion for the phase field.
 * It is just a child of the <code> Function </code> class template.
 */
template<int dim> class Phase0: public Function<dim>{
  public:
    /**
     * Constructor.
     * @param dd : The interface thickness.
     * @param _tc : The test case.
     */
    Phase0( const double dd, const TestCase _tc );
    /**
     * Destructor. does nothing
     */
    virtual ~Phase0();
    /**
     * The value. As usual:
     * @param p : The point at which the function is evaluated.
     * @param component : The component of the (vector valued) function.
     */
    virtual double value( const Point<dim> &p, const unsigned component = 0 ) const;
  protected:
    /**
     * The interface thickness
     */
    const double delta;
    /**
     * The testcase
     */
    const TestCase tc;
};

/**
 * A Cahn Hilliard equation with surface transport terms:
 * \f[
 *   \phi_t + \mathbf{u} \cdot \nabla \phi = \nabla\cdot( M(\phi) \nabla \mu ), \quad 
 *   \text{ in } \Omega
 * \f]
 * \f[
 * \mu = \gamma \left( \frac1\delta \mathcal{W}'(\phi) - \delta \Delta \phi \right)
 *  - \frac12\varepsilon'(\phi)|\nabla V|^2 + \frac12 \rho'(\phi)|\mathbf{u}|^2,  \quad 
 *   \text{ in } \Omega
 * \f]
 * with boundary conditions
 * \f[
 *   \alpha\left( \phi_t + \mathbf{u}_{\boldsymbol{\tau}} \partial_{\boldsymbol{\tau}} \phi \right)
 *   + \gamma \left( \Theta_{fs}'(\phi) + \delta \partial_{\mathbf{n}} \phi \right) = 0,
 *   \quad
 *   M(\phi) \partial_n \mu = 0, \quad \text{on }\Gamma,
 * \f]
 */
template<int dim> class CahnHilliard: public Problem<dim>{
  public:
    /**
     * Constructor. 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 phase field 
     *    and chemical potential are defined.
     * @param _tc : The testcase. Needed to know which initial condition to invoke.
     */
    CahnHilliard( const Triangulation<dim> &tria, Material_Parameters &params, const unsigned deg,
                  const TestCase _tc );
    /**
     * Clear the data.
     */
    virtual ~CahnHilliard();
    /**
     * Advance in time.
     */
    void advance();
    /**
     * compute the iterative error.
     */
    double it_error();
  protected:
    const material_function &mobility, ///<Mobility of the phases.
                            &permitivity, ///< Electric permittivity of the phases.
                            &density; ///< Density of the phases.
    const GLPotential &W; ///< The derivative of the Ginzburg Landau potential.
    const InterfaceEnergy &gamma_fs; ///< The derivative of the interface fee energy.
    const double gamma, ///< Suface tension.
                 delta, ///< Interface thickness.
                 alpha; ///< Surface transport parameter.
    const TestCase tc; ///< The test case.
    TrilinosWrappers::Vector oldsol, ///< The old solution
                             olditsol; ///< Solution at the previous iterative step in the inner loop.
    /**
     * What to do after the LA data is initialized
     */
    virtual void InitLADataSuffix();
    /**
     * We group the iterators of this problem and the problems this one takes information from
     * as follows:
     * <ul>
     *   <li> <code>Iterator[0]</code>: This.</li>
     *   <li> <code>Iterator[1]</code>: Voltage.</li>
     *   <li> <code>Iterator[2]</code>: Velocity.</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 <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 local FEValues object
      hp::QCollection<dim-1> f_quad; ///< The collection of face quadrature formulas
      UpdateFlags f_flags; ///< The collection of face update flags
      hp::FEFaceValues<dim> fe_face_val; ///< The local FEFaceValues object
      AsFunction<dim> voltage, ///< The AsFunction object that will be interpreted as Voltage
                      velocity; ///< The AsFunction object that will be interpreted as Velocity
      std::vector< Tensor<1,dim> > loc_grad_voltage, loc_vel, loc_face_vel, phase_grad, phase_face_grad;
      std::vector< Point<dim> > normals;
      std::vector<double> phase_val, phase_face_val;
      ScratchData( const hp::FECollection<dim> &fe, std::vector<AsFunction<dim> *> &data,
                   const hp::QCollection<dim> &q, const UpdateFlags u_flags,
                   const hp::QCollection<dim-1> &fq, const UpdateFlags u_f_flags ) :
          quad(q), flags( u_flags ), fe_val( fe, quad, flags ),
          f_quad(fq), f_flags( u_f_flags ), fe_face_val( fe, f_quad, f_flags ),
          voltage( *(data[0]) ), velocity( *(data[1]) ),
          loc_grad_voltage( quad.max_n_quadrature_points() ), loc_vel( quad.max_n_quadrature_points() ),
          loc_face_vel( f_quad.max_n_quadrature_points() ), phase_grad( quad.max_n_quadrature_points() ),
          phase_face_grad( f_quad.max_n_quadrature_points() ),
          normals( f_quad.max_n_quadrature_points() ),
          phase_val( quad.max_n_quadrature_points() ),
          phase_face_val( f_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 ),
          f_quad( scratch.f_quad ), f_flags( scratch.f_flags ),
          fe_face_val( scratch.fe_face_val.get_fe_collection(), f_quad, f_flags ),
          voltage( scratch.voltage ), velocity( scratch.velocity ),
          loc_grad_voltage( quad.max_n_quadrature_points() ), loc_vel( quad.max_n_quadrature_points() ),
          loc_face_vel( f_quad.max_n_quadrature_points() ), phase_grad( quad.max_n_quadrature_points() ),
          phase_face_grad( f_quad.max_n_quadrature_points() ),
          normals( f_quad.max_n_quadrature_points() ),
          phase_val( quad.max_n_quadrature_points() ),
          phase_face_val( f_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> Velocity. </li>
     * </ol>
     * @param time : The current time.
     * @param m_theaded : 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_theaded = 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 );
    /**
     * Does nothing. Needed for compatibility
     */
    virtual void ReinitPrec();
    /**
     * What to do before we solve the system
     */
    virtual void SolveSystemPreffix();
    /**
     * Solve the system to given accuracy
     */
    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_phase, ///< The phase field
                    get_potential, ///< The chemical potential
                    get_old_phase; ///< The old phase field
};

#endif
