#ifndef _NSE_H_
#define _NSE_H_

#include "Problem.h"

/**
 * The elliptic problem that determines the pressure. After time discretization, at each time
 * step, we solve
 * \f[
 * \Delta \xi^{k+1} = \frac{\varrho}{\Delta t} \nabla\cdot \mathbf{u}^{k+1},
 * \quad \text{ in } \Omega,
 * \f]
 * with boundary condition
 * \f[
 *    \partial_n \xi^{k+1} = 0, \quad \text{ on } \Gamma
 * \f]
 * After which, the pressure is determined by
 * \f[
 *    p^{k+1} = p^k + \xi^{k+1}.
 * \f]
 * For this reason, the "solution" in this problem is the variable $\xi$ and the pressure
 * is just a view, that is an instance of the AsFunction class template.
 */
template<int dim> class Pressure: public Problem<dim>{
  public:
    /**
     * Initialize the members:
     * @param tria : A reference to the mesh where this problem is defined.
     * @param params : The list of material parameters.
     * @param deg : The polynomial degree
     * @param ee : The accuracy to which we iteraritvely solve the problem.
     * @param thres : The aggregation threshold for the AMG preconditioner.
     * @param sweeps  : The number of smoothing steps the AMG preconditioner must perform.
     *
     * Notice that in this problem the matrix does not change and, therefore, the preconditioner
     * does not need to be updated.
     */
    Pressure( const Triangulation<dim> &tria, Material_Parameters &params, const unsigned deg,
              const double ee, const double thres, const unsigned sweeps );
    /**
     * Clear parameters
     */
    virtual ~Pressure();
  protected:
    const double rho_min; ///< Penalization parameter
    TrilinosWrappers::PreconditionAMG prec; ///< Preconditioner (AMG)
    TrilinosWrappers::PreconditionAMG::AdditionalData prec_data; ///< Preconditioner data
    TrilinosWrappers::Vector p, ///< pressure
                             pres_extr; ///< extrapolated pressure $=2p^k - p^{k-1}$.
    /**
    * Since the problem we are solving has a homogeneous Neumann boundary condition
    * we obtain a singular matrix. To remove the singularity we add to the ConstraintMatrix
    * one extra constraint, fixing the value of the solution at one degree of freedom.
    */
    virtual void SetupDoFsSuffix();
    /**
     * We just set the auxiliary vectors to the proper size
     */
    virtual void InitLADataSuffix();
    /**
     * 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]: Velocity </li>
     * </ul>
     */
    typedef IteratorGroup<2, 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 <code>FEValues</code> object
      AsFunction<dim> velocity; ///< The AsFunction object that is interpreted as Velocity
      std::vector<double> loc_div_vel;
      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 ), velocity( *( data[0] ) ),
          loc_div_vel( 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 ), velocity( scratch.velocity ),
          loc_div_vel( 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> 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 );
    /**
     * Since $p^{k+1}$ and $p^{\star,k+1}$ are actually views of the problem. 
     * After solving the system, we need to
     * update them.
     */
    virtual void SolveSystemSuffix();
    /**
     * 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_pressure, ///< pressure $p^{k+1}$
                    get_extrapolated_pressure; ///< extrapolated pressure $p^k-p^{k-1}$
};

/**
 * The vector valued convection diffusion problem that defines the evolution of the velocity
 * \f[
 *   \frac{ D(\rho(\phi){\mathbf{u}}) }{Dt} - \nabla\cdot\left( \eta(\phi) \mathbf{S}({\mathbf{u}}) \right)
 *   + \nabla p = \mu \nabla \phi - q \nabla \left( V + \lambda q \right)
 *   + \frac12 \rho'(\phi)\phi_t {\mathbf{u}}, \quad \text{on }\Omega,
 * \f]
 * where
 * \f[
 *   \frac{ D(\rho(\phi){\mathbf{u}}) }{Dt} = \sigma(\sigma \mathbf{u})_t + \frac12 \nabla\cdot(\rho \mathbf{u})\mathbf{u},
 *   \quad
 *   \mathbf{S}(\mathbf{u}) = \frac12 \left( \nabla \mathbf{u} + {\nabla \mathbf{u}}^T \right).
 * \f]
 * This system is supplemented with the boundary conditions
 * \f[
 *  {\mathbf{u}}\cdot{\mathbf{n}} = 0, \qquad
 *    \beta(\phi) {\mathbf{u}}_{\boldsymbol{\tau}} + \eta(\phi) \mathbf{S}({\mathbf{u}})_{{\mathbf{n}}{\boldsymbol{\tau}}}
 *        = \gamma\left( \Theta_{fs}'(\phi) + \delta \partial_{\mathbf{n}} \phi \right) \partial_{\boldsymbol{\tau}} \phi,
 *        \quad \text{on }\Gamma,
 * \f]
 */
template<int dim> class Velocity: 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 v_sweeps  : The number of smoothing steps the AMG preconditioner must perform.
     * @param K_size : The dimension of the Krylov subspace (used in GMRES).
     */
    Velocity( const Triangulation<dim> &tria, Material_Parameters &params, const unsigned deg, const unsigned u_prec,
              const double ee, const double thres, const unsigned v_sweeps, const unsigned K_size );
    /**
     * Clear the data
     */
    virtual ~Velocity();
    /**
     * Advance in time. CHNSE calls this function once the inner iteration loop has converged.
     */
    void advance();
    /**
     * Compute the inner iteration error. CahnHilliard calls this function to verify wether the
     * inner iteration has converged.
     */
    double it_error();
  protected:
    const material_function &density, ///< density
                            &viscosity, ///< viscosity
                            &slip_coeff; ///< slip coefficient
    InterfaceEnergy &gamma_fs; ///< interface free energy
    const double alpha, ///< surface transport parameter
                 lambda, ///< charge coupling parameter
                 gamma, ///< surface tension
                 delta; ///< interface thickness
    const unsigned Krylov_size; ///< dimension of the Krylov subspace
    TrilinosWrappers::PreconditionAMG prec; // Preconditioner (AMG)
    TrilinosWrappers::PreconditionAMG::AdditionalData prec_data; ///< Preconditioner data
    TrilinosWrappers::Vector oldsol, ///< the old solution
                             olditsol; ///< the old iterative solution
    /**
     * What to do after the Linear Algebra data has been initialized.
     * 
     * We set <code>oldsol</code> and <code>olditsol</code> to the correct size.
     */
    virtual void InitLADataSuffix();
    /**
     * What to do after the DoFs are set.
     * 
     * We add the impermeability conditions as constraints to the ConstraintMatrix object.
     * Notice that since this is applied not only at physical boundaries, but also at the
     * interface between fluid and plates, special care must be taken. In general this seems like
     * a difficult process, since the normals might not be well defined. What we do is assume
     * that the interface is horizontal, so that by using the <code>FESystem::system_to_component_index</code>
     * we know which local DOF corresponds to the vertical velocity.
     */
    virtual void SetupDoFsSuffix();
    /**
     * We group the iterators of this problem and the problems this one takes information from as follows:
     * <ul>
     *  <li> <code>Iterator</code>[0]: Phase field from CahnHilliard </li>
     *    <li>  <code>Iterator</code>[1]: Chemical Potential from CahnHilliard </li>
     *    <li>  <code>Iterator</code>[2]: Charge </li>
     *    <li>  <code>Iterator</code>[3]: Voltage </li>
     *    <li>  <code>Iterator</code>[4]: %Pressure extrapolation from Pressure </li>
     *    <li>  <code>Iterator</code>[5]:  Value of the phase field variable at previous time step from CahnHilliard</li>
     * </ul>
     */
    typedef IteratorGroup<7,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 Parallel 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
      hp::QCollection<dim-1> fquad; ///< The collection of face quadrature formulas
      UpdateFlags f_flags; ///< The face update flags
      hp::FEFaceValues<dim> fe_face_val; ///< The FEFaceValues object
      AsFunction<dim> phase, ///< Phase field, part of CahnHilliard
                      mu, ///< Chemical potential, part of CahnHilliard
                      q, ///< Charge
                      V, ///< Voltage
                      psharp, ///< Extrapolated pressure, part of Pressure
                      oldphase; ///< Phase field at previous time, part of CahnHilliard
      std::vector<double> loc_phase, loc_old_phase, face_phase, face_old_phase, loc_q,
                          loc_psharp, loc_divu, loc_mu;
      std::vector< Tensor<1,dim> > loc_grad_q, loc_grad_V, loc_grad_phase, face_grad_phase, loc_u, face_base,
                                   loc_face_u;
      Point<dim> normal;
      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 ),
          fquad( fq ), f_flags( u_f_flags ), fe_face_val( fe, fquad, f_flags ),
          phase( *( data[0] ) ), mu( *( data[1] ) ), q( *( data[2] ) ), V( *( data[3] ) ),
          psharp( *( data[4] ) ), oldphase( *( data[5] ) ),
          loc_phase( quad.max_n_quadrature_points() ), loc_old_phase( quad.max_n_quadrature_points() ),
          face_phase( fquad.max_n_quadrature_points() ), face_old_phase( fquad.max_n_quadrature_points() ),
          loc_q( quad.max_n_quadrature_points() ), loc_psharp( quad.max_n_quadrature_points() ),
          loc_divu( quad.max_n_quadrature_points() ), loc_mu( quad.max_n_quadrature_points() ),
          loc_grad_q( quad.max_n_quadrature_points() ), loc_grad_V( quad.max_n_quadrature_points() ),
          loc_grad_phase( quad.max_n_quadrature_points() ), face_grad_phase( fquad.max_n_quadrature_points() ),
          loc_u( quad.max_n_quadrature_points() ), face_base( fe.max_dofs_per_cell() ),
          loc_face_u( fquad.max_n_quadrature_points() ) {}
      ScratchData( const ScratchData &scratch ) :
          quad( scratch.quad ), flags( scratch.flags ), fe_val( scratch.fe_val.get_fe_collection(), quad, flags ),
          fquad( scratch.fquad ), f_flags( scratch.f_flags ),
          fe_face_val( scratch.fe_face_val.get_fe_collection(), fquad, f_flags ),
          phase( scratch.phase ), mu( scratch.mu ), q( scratch.q ), V( scratch.V ),
          psharp( scratch.psharp ), oldphase( scratch.oldphase ),
          loc_phase( quad.max_n_quadrature_points() ), loc_old_phase( quad.max_n_quadrature_points() ),
          face_phase( fquad.max_n_quadrature_points() ), face_old_phase( fquad.max_n_quadrature_points() ),
          loc_q( quad.max_n_quadrature_points() ), loc_psharp( quad.max_n_quadrature_points() ),
          loc_divu( quad.max_n_quadrature_points() ), loc_mu( quad.max_n_quadrature_points() ),
          loc_grad_q( quad.max_n_quadrature_points() ), loc_grad_V( quad.max_n_quadrature_points() ),
          loc_grad_phase( quad.max_n_quadrature_points() ), face_grad_phase( fquad.max_n_quadrature_points() ),
          loc_u( quad.max_n_quadrature_points() ), face_base( fe_val.get_fe_collection().max_dofs_per_cell() ),
          loc_face_u( fquad.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> Phase field from CahnHilliard </li>
     * <li> Chemical potential from CahnHilliard </li>
     * <li> Charge </li>
     * <li> Voltage </li>
     * <li> Extrapolated pressure from Pressure </li>
     * <li> Phase field at previous time step 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 ) const;
    /**
     * What to do before solving the system. Since the solving is always part of an inner iteration loop, we
     * save the previous iterate so we can check convergence.
     */
    virtual void SolveSystemPreffix();
    /**
     * 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 initial data. In this case $\mathbf{u}_0 = 0$.
     */
    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_velocity; ///< velocity
};

#endif // _NSE_H_
