#ifndef _PROBLEM_H_
#define _PROBLEM_H_

#include <deal.II/base/point.h>
#include <deal.II/base/function.h>
#include <deal.II/base/multithread_info.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/work_stream.h>
#include <deal.II/base/parallel.h>

#include <deal.II/grid/tria.h>

#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/fe_collection.h>

#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/function_map.h>

#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_tools.h>

#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/compressed_sparsity_pattern.h>
#include <deal.II/lac/constraint_matrix.h>

#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/solution_transfer.h>

#include "Material.h"
#include "AsFunction.h"
#include "GroupedIterators.h"

using namespace dealii;

/**
 * Once a problem is discretized in time and space, its solution basically amounts to
 * defining a finite element space, assembling a matrix and solving a linear system.
 * This class is an attempt at maximizing the reuse of the common parts of this process.
 *
 * The way we try to maximize this is by wrapping around all the common features of this
 * procedures in this class. There is, however, some caveats. It is not possible -- or at least
 * we do not know how --  to fully abstract some parts of this process. For instace, different
 * solvers might be needed depending on whether the problem at hand is symmetric or not. For
 * this reason, we chose to set these as pure virtual functions and the child classes should
 * implement them.
 */
template<int dim> class Problem{
  public:
    /**
     * Copy the data to members and/or initialize them
     * @param tria : A reference to the triangulation where this problem is defined
     * @param _n_transfers : We want to be able to use mesh adaptation. Since most likely
     * our problem is time dependent we want to be able to keep the solution when doing this.
     * In addition, there might be some auxiliary quantities that we want to keep as well.
     * This will indicate how many of these extra quantities we have.
     * @param _eps : The tolerance for the iterative solver.
     * @param _prec_needs_reinit : A flag that indicates if the matrix of this problem 
     * will change during the existence of the object. If it does then we will need to 
     * reinitialize the preconditioner every once in a while.
     * @param _update_prec : If the preconditioner needs reinitialization, then this indicates
     * what exactly we mean by "once in a while".
     */
    Problem( const Triangulation<dim> &tria, const unsigned _n_transfers, const double _eps,
             bool _prec_needs_reinit = false, const unsigned _update_prec = 10 );
    /**
     * Clear data
     */
    virtual ~Problem();
    /**
     * Set the time step
     */
    void set_dt( const double _dt );
    /**
     * Initialize the problem. This requires:
     * <ol>
     *  <li> Set up the Degrees of Freedom and constraints </li>
     *  <li> Initialize the linear algebra data </li>
     *  <li> Set the initial data </li>
     *  <li> Connect PreRefinement() and PostRefinement() to the signals emmited 
     *    by Triangulation before and after refinement, respectively.
     * </ol>
     */
    void init();
    /**
     * We assemble and solve the system.
     * @param data : The views from other problems this one uses as data or coefficients.
     * @param time : The current time
     */
    virtual void solve( std::vector< AsFunction<dim> *> &data, const double time );
    /**
     * Returns a constant reference to the used FiniteElement object
     */
    const hp::FECollection<dim> &get_fe() const;
    /**
     * Return the number of degrees of freedom
     */
    unsigned size() const;
  protected:
    const Triangulation<dim> &tri; ///< Reference to the triangulation
    bool mesh_has_changed; ///< Flag to indicate that the mesh changed so that we can reinit the prec
    hp::DoFHandler<dim> dh; ///< The DoF handler object
    ConstraintMatrix constraints; ///< The hanging node (and other) constraints
    hp::FECollection<dim> fe; ///< the finite element
    unsigned n_dofs; ///< number of degrees of freedom
    double dt; ///< the time-step
    unsigned steps; ///< the number of steps
    TrilinosWrappers::Vector sol, ///< The solution to the problem
                             rhs; ///< The RHS
    TrilinosWrappers::SparseMatrix K; ///< The system matrix
    const double eps; ///< Tolerance for the iterative solver
    const unsigned update_prec; ///< How often to update the preconditoner
    const bool prec_needs_reinit; ///< Flag to indicate if the preconditioner needs reinitialization
    /**
     * We distribute the degrees of freedom, set the number of degrees of freedom
     * and make the hanging node constraints.
     *
     */
    void SetupDoFs();
    /**
     * It is possible that some constraints must be added to the problem at hand. For instance,
     * if the problem has Dirichlet or Neumann boundary conditions. It is also possible that
     * something else needs to be done. We will call this method just from SetupDoFs() just
     * before closing the ConstraintMatrix so that the child classes can do what they
     * want there.
     */
    virtual void SetupDoFsSuffix();
    /**
     * We initialize the members that are related to the linear algebra part of the solution
     * process. That is, we create the sparsity pattern, which in turn is needed to
     * initialize the system matrix. Then we set the solution and right hand side vectors
     * to the correct size.
     */
    void InitLAData();
    /**
     * In a particular problem there might be more linear algebra related objects that 
     * need initialization. For instace if we are keeping extrapolations of the solution
     * or the solution at previous time-steps. This cannot be known in advance and so,
     * we add a suffix that will be called at the end InitLAData() and that the 
     * child classes will overload.
     */
    virtual void InitLADataSuffix();
    /**
     * The data for assembling the problem. This is the PerTaskData in the sense of the
     * parallel module in deal.II. In other words, this only needs to bundle the number
     * of degrees of freedom per cell, the local matrix and vector.
     */
    struct PerTaskData{
      unsigned dpc; /// The number of degrees of freedom per cell
      std::vector<unsigned> ldi; ///< The local to global map for the DoFs
      FullMatrix<double> loc_m; ///< The local matrix
      Vector<double> loc_rhs; ///< The local vector
      PerTaskData( const unsigned _dpc ) : dpc( _dpc ), ldi( dpc ), loc_m( dpc, dpc ),
                loc_rhs( dpc ) {}
      PerTaskData( const PerTaskData &data ) : dpc( data.dpc ), ldi( dpc ),
                loc_m( dpc, dpc ), loc_rhs( dpc ) {}
    };
    /**
     * Assemble the system. There is a lot of common structure in the assembly of the system,
     * since the only thing that is really problem-specific is the expression to assemble
     * the local matrices and vectors. However, due to the structure used by the 
     * deal.II parallel routines, i.e. division in Scratch and PerTask data, this method
     * needs to be made pure virtual, because there is now way of knowing beforehand what kind
     * of scratch data might be needed for the assembly of the local problems.
     * @param data : The AsFunction objects that the problem uses as coefficients and or
     *      right hand side.
     * @param time : The current time
     * @param m_threaded : Indicates if the assembly routine must be performed in parallel
     */
    virtual void AssembleSystem( std::vector< AsFunction<dim> *> &data, const double time,
                                  const bool m_threaded = true ) = 0;
    /**
     * This <b>IS</b> the assembly procedure. That is, it adds the local matrices and right hand
     * sides to the global ones, taking into account the constraints in the system.
     * @param data : This bundles the local matrix, vector and the local to global map.
     */
    void CopyToGlob( const PerTaskData &data );
    /**
     * Do all the bookkeeping needed before solving the system, actually solve it and do
     * all the bookkeeping that is needed after the system is solved. In other words, this
     * method does the following:
     * <ol>
     *    <li> Call SolveSystemPreffix()</li>
     *    <li> Check if the preconditioner needs to be updated and if it does, reinitialize it </li>
     *    <li> Create the SolverControl object that will drive the solution process </li>
     *    <li> Do the actual solution </li>
     *    <li> Call SolveSystemSuffix() </li>
     * </ol>
     */
    void SolveSystem();
    /**
     * What to do before solving the system. For instance we might need to save the
     * old solution.
     */
    virtual void SolveSystemPreffix();
    /**
     * Reinit the preconditioner. The type of preconditioner to be used
     * <b>must be </b> problem specific, and so it should be its initialization.
     * This justifies this function being pure virtual.
     */
    virtual void ReinitPrec() = 0;
    /**
     * Call the actual solver. The solution procedure <b>must be</b> problem specific,
     * so this is pure virtual
     */
    virtual void DoSolve( SolverControl &control ) = 0;
    /**
     * What to do after the system is solved. For instance, we might need to 
     * update extrapolations.
     */
    virtual void SolveSystemSuffix();
    /**
     * Set the initial data
     */
    virtual void SetInitialData() = 0;
    const unsigned n_of_transfers; ///< Indicates how many vectors we want to transfer on each mesh
    std::vector<TrilinosWrappers::Vector> x_sol; ///< The vectors that are going to be transferred
    SolutionTransfer<dim, TrilinosWrappers::Vector, hp::DoFHandler<dim> > transfer; ///< The transferring mechanism
    /**
     * Executed before the refinement
     */
    void PreRefinement();
    /**
     * Since we do not know how to fill x_sol, we deferr it to the child classes
     */
    virtual void PreRefinementPreffix() = 0;
    /**
     * Executed after the refinement. We call SetupDofs() and InitLAData(). Then we check if
     * we are at the beginning of time. If we are we call SetInitialData(), otherwise 
     * we transfer the vectors to the new mesh.
     */
    void PostRefinement();
    /**
     * We do not know where x_sol came from, so its transfer is deferred to the child classes
     */
    virtual void PostRefinementSuffix( const std::vector<TrilinosWrappers::Vector> &sol_tmp ) = 0;
    /**
     * In our application, special care must be taken at the interface between the fluid
     * and the plates. This method determines wether a face is there.
     * @param c : The cell
     * @param face : The face
     */
    bool isNeeded( const typename hp::DoFHandler<dim>::active_cell_iterator &c,
                   const unsigned face ) const;
};

#endif
