#ifndef _AS_FUNCTION_H_
#define _AS_FUNCTION_H_

#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/numerics/data_out.h>

using namespace dealii;

/**
 * When we have more than one coupled PDEs, we usually need to use the computed solution of a subproblem
 * as right hand side or coefficient of the other subproblems.
 * This class implements this technology, by storing a <code>FEValues</code> object,
 * being able to loop over all cells and, more importantly, call the <code>get_function_*</code> methods
 * with a fixed global vector which, usually, is the solution of a given problem or an easily derived
 * quantity.
 */
template<int dim> class AsFunction{
  public:
    /**
     * In deal.II the <code>FEValuesExtractors</code> do not have a common type
     * therefore we will need to keep track of the type of function (scalar/vector valued)
     * and treat the <code>get_function_*</code> methods accordingly
     */ 
    enum FunctionType{
      isScalarValued, ///< The function is scalar valued
      isVectorValued ///< The function is vector valued
    };
    /**
     * @param _dh : The <code>DoFHandler</code> object that is needed to interpret vectors as functions.
     * @param _sol : The coefficient vector that represent the discrete function.
     * @param _name : The name of this function. Used for plotting.
     * @param _type : Indicates if the function must be interpreted as scalar or vector valued
     * @param fc : First component. This is used if, for instance, the finite element space is vector valued
     * and we want only one of the components of the space.
     */
    AsFunction( const hp::DoFHandler<dim> &_dh, const TrilinosWrappers::Vector &_sol, const char* _name,
                 const FunctionType _type = isScalarValued, const unsigned fc = 0 );
    /**
     * Copy constructor
     */
    AsFunction( const AsFunction<dim> &_asfunction );
    /**
     * Clear allocated memory
     */
    ~AsFunction();
    /**
     * Reset the object,
     * @param q : is the new quadrature
     * @param f : the update flags
     */
    void reset( const hp::QCollection<dim> &q, const UpdateFlags f );
    /**
     * Clear the object
     */
    void reset();
    /**
     * Reinit to a new cell
     */
    void reinit( const typename hp::DoFHandler<dim>::active_cell_iterator &cell );
    /**
     * Reset the face object
     * @param q : is the new quadrature
     * @param f : the update flags
     */
    void reset_face( const hp::QCollection<dim-1> &q, const UpdateFlags f );
    /**
     * Clears the face object
     */
    void reset_face();
    /**
     * Reinit to a new face.
     * @param cell : the cell
     * @param face : the face
     */
    void reinit( const typename hp::DoFHandler<dim>::active_cell_iterator &cell, const unsigned face );
    /**
     * Maps to <code>FEValues.get_function_values</code> for scalars
     */
    void get_function_values( std::vector<double> &values ) const;
    /**
     * Maps to <code>FEValues.get_function_gradients</code> for scalars
     */
    void get_function_gradients( std::vector< Tensor<1, dim> > &values ) const;
    /**
     * Maps to <code>FEValues.get_function_values</code> for vectors
     */
    void get_function_values( std::vector< Tensor<1, dim> > &values ) const;
    /**
     * Maps to <code>FEValues.get_function_gradients</code> for vectors
     */
    void get_function_gradients( std::vector< Tensor<2, dim> > &values ) const;
    /**
     * Maps to <code>FEValues.get_function_divergences</code> (valid only for vectors)
     */
    void get_function_divergences( std::vector<double> &values ) const;
    /**
     * Maps to <code>FEFaceValues.get_function_values</code> for scalars
     */
    void get_function_face_values( std::vector<double> &values ) const;
    /**
     * Maps to <code>FEFaceValues.get_function_gradients</code> for scalars
     */
    void get_function_face_gradients( std::vector< Tensor<1, dim> > &values ) const;
    /**
     * Maps to <code>FEFaceValues.get_function_values</code> for vectors
     */
    void get_function_face_values( std::vector< Tensor<1, dim> > &values ) const;
    /**
     * Maps to <code>FEFaceValues.get_function_gradients</code> for vectors
     */
    void get_function_face_gradients( std::vector< Tensor<2, dim> > &values ) const;
    /**
     * Raw access to the data vector
     */
    const TrilinosWrappers::Vector & get_vec() const;
    /**
     * Raw access to an entry of the data vector
     */
    const double operator()( unsigned i ) const;
    /**
     * Get the associated <code>DoFHandler</code> object
     */
    const hp::DoFHandler<dim> & get_dh() const;
    /**
     * Get the beginning of the cell iterator associated with this <code>DoFHandler</code>
     */
    typename hp::DoFHandler<dim>::active_cell_iterator get_begin() const;
    /**
     * Get the end of the cell iterator associated with this <code>DoFHandler</code>
     */
    typename hp::DoFHandler<dim>::active_cell_iterator get_end() const;
  protected:
    /**
     * The internal cell quadrature rule.
     */
    mutable hp::QCollection<dim> *quad;
    /**
     * The internal cell update flags
     */
    mutable UpdateFlags flags;
    /**
     * The internal <code>FEValues</code> object
     */
    mutable hp::FEValues<dim> *fe_val;
    /**
     * The internal face quadrature rule
     */
    mutable hp::QCollection<dim-1> *f_quad;
    /**
     * The internal face update flags
     */
    mutable UpdateFlags f_flags;
    /**
     * The internal <code>FEFaceValues</code> object
     */
    mutable hp::FEFaceValues<dim> *fe_face_val;
    /**
     * Points to <code>DoFHandler</code>
     */
    SmartPointer< const hp::DoFHandler<dim> > dh;
    /**
     * Points to the coefficient vector
     */
    SmartPointer< const TrilinosWrappers::Vector > sol;
    /**
     * The type of function that we have (scalar/vector)
     */
    const FunctionType type;
    /**
     * The first component of the function (meaningful if the <code>DoFHandler</code> is vector valued)
     */
    const unsigned first_component;
  public:
    /**
     * The function name
     */
    const std::string name;
    /**
     * The interpetation of the function (scalar/vector)
     */
    const DataComponentInterpretation::DataComponentInterpretation interpretation;
  private:
    /**
     * Disable default constructor
     */
    AsFunction();
};

#endif
