#include "../include/AsFunction.h"


// @sect2{The AsFunction class template}
// In the constructor we just store a reference to the <code>DoFHandler</code> object we want
// to use to interpret the <code>Vector</code> as a function

/**
 * Constructor.
 * 
 */
template<int dim> AsFunction<dim>::AsFunction( const hp::DoFHandler<dim> &_dh, const TrilinosWrappers::Vector &_sol,
                                               const char* _name, const FunctionType _type, const unsigned fc ):
    quad( NULL ),fe_val( NULL ),
    f_quad( NULL ), fe_face_val( NULL ),
    dh( &_dh ), sol( &_sol ),
    type( _type ), first_component( fc ), name( _name ),
    interpretation( type == isScalarValued ? DataComponentInterpretation::component_is_scalar
                     : DataComponentInterpretation::component_is_part_of_vector ) {}

template<int dim> AsFunction<dim>::AsFunction( const AsFunction<dim> &_asfunction ):
    quad( NULL ), flags( _asfunction.flags ),
    fe_val( NULL ),
    f_quad( NULL ), f_flags( _asfunction.f_flags ),
    fe_face_val( NULL ),
    dh( _asfunction.dh ), sol( _asfunction.sol ),
    type( _asfunction.type ), first_component( _asfunction.first_component ), name( _asfunction.name ),
    interpretation( _asfunction.interpretation ) {
  if( _asfunction.quad )
    quad = new hp::QCollection<dim>( *_asfunction.quad );
  if( quad )
    fe_val = new hp::FEValues<dim>( dh->get_fe(), *quad, flags );
  if( _asfunction.f_quad )
    f_quad = new hp::QCollection<dim-1>( *_asfunction.f_quad );
  if( f_quad )
    fe_face_val = new hp::FEFaceValues<dim>( dh->get_fe(), *f_quad, f_flags );
}

template<int dim> AsFunction<dim>::~AsFunction(){
  dh = NULL;
  sol = NULL;
  reset();
  reset_face();
}

template<int dim> void AsFunction<dim>::reset(){
  if( fe_val ){
    delete fe_val;
    fe_val = NULL;
  }
  if( quad ){
    delete quad;
    quad = NULL;
  }
}

template<int dim> void AsFunction<dim>::reset( const hp::QCollection<dim> &q, const UpdateFlags f ){
  reset();
  quad = new hp::QCollection<dim>( q );
  flags = f;
  fe_val = new hp::FEValues<dim>( dh->get_fe(), *quad, flags );
}

template<int dim> void AsFunction<dim>::reinit( const typename hp::DoFHandler<dim>::active_cell_iterator &cell ){
  Assert( fe_val, ExcInvalidState() );
  fe_val->reinit( cell );
}

template<int dim> void AsFunction<dim>::reset_face(){
  if( fe_face_val ){
    delete fe_face_val;
    fe_face_val = NULL;
  }
  if( f_quad ){
    delete f_quad;
    f_quad = NULL;
  }
}

template<int dim> void AsFunction<dim>::reset_face( const hp::QCollection <dim-1> &q, const UpdateFlags f ){
  reset_face();
  f_quad = new hp::QCollection<dim-1>( q );
  f_flags = f;
  fe_face_val = new hp::FEFaceValues<dim>( dh->get_fe(), *f_quad, f_flags );
}

template<int dim> void AsFunction<dim>::reinit( const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
                                                const unsigned face ){
  Assert( fe_face_val, ExcInvalidState() );
  fe_face_val->reinit( cell, face );
}

template<int dim> void AsFunction<dim>::get_function_values( std::vector<double> &values ) const{
  Assert( type == isScalarValued, ExcInvalidState() );
  FEValuesExtractors::Scalar aa( first_component );
  fe_val->get_present_fe_values()[aa].get_function_values( *sol, values );
}

template<int dim> void AsFunction<dim>::get_function_values( std::vector< Tensor<1, dim> > &values ) const{
  Assert( type == isVectorValued, ExcInvalidState() );
  FEValuesExtractors::Vector aa( first_component );
  fe_val->get_present_fe_values()[aa].get_function_values( *sol, values );
}

template<int dim> void AsFunction<dim>::get_function_gradients( std::vector< Tensor<1, dim> > &values ) const{
  Assert( type == isScalarValued, ExcInvalidState() );
  FEValuesExtractors::Scalar aa( first_component );
  fe_val->get_present_fe_values()[aa].get_function_gradients( *sol, values );
}

template<int dim> void AsFunction<dim>::get_function_gradients( std::vector< Tensor<2, dim> > &values ) const{
  Assert( type == isVectorValued, ExcInvalidState() );
  FEValuesExtractors::Vector aa( first_component );
  fe_val->get_present_fe_values()[aa].get_function_gradients( *sol, values );
}

template<int dim> void AsFunction<dim>::get_function_divergences( std::vector<double> &values ) const{
  Assert( type == isVectorValued, ExcInvalidState() );
  FEValuesExtractors::Vector aa( first_component );
  fe_val->get_present_fe_values()[aa].get_function_divergences( *sol, values );
}


template<int dim> void AsFunction<dim>::get_function_face_values( std::vector<double> &values ) const{
  Assert( type == isScalarValued, ExcInvalidState() );
  FEValuesExtractors::Scalar aa( first_component );
  fe_face_val->get_present_fe_values()[aa].get_function_values( *sol, values );
}

template<int dim> void AsFunction<dim>::get_function_face_values( std::vector< Tensor<1, dim> > &values ) const{
  Assert( type == isVectorValued, ExcInvalidState() );
  FEValuesExtractors::Vector aa( first_component );
  fe_face_val->get_present_fe_values()[aa].get_function_values( *sol, values );
}

template<int dim> void AsFunction<dim>::get_function_face_gradients( std::vector< Tensor<1, dim> > &values ) const{
  Assert( type == isScalarValued, ExcInvalidState() );
  FEValuesExtractors::Scalar aa( first_component );
  fe_face_val->get_present_fe_values()[aa].get_function_gradients( *sol, values );
}

template<int dim> void AsFunction<dim>::get_function_face_gradients( std::vector< Tensor<2, dim> > &values ) const{
  Assert( type == isVectorValued, ExcInvalidState() );
  FEValuesExtractors::Vector aa( first_component );
  fe_face_val->get_present_fe_values()[aa].get_function_gradients( *sol, values );
}

template<int dim> const TrilinosWrappers::Vector &AsFunction<dim>::get_vec() const{
  return *sol;
}

template<int dim> const double AsFunction<dim>::operator()( unsigned i ) const{
  return ( *sol )( i );
}

template<int dim> const hp::DoFHandler<dim> &AsFunction<dim>::get_dh() const{
  return *dh;
}

template<int dim> typename hp::DoFHandler<dim>::active_cell_iterator AsFunction<dim>::get_begin() const{
  return dh->begin_active();
}

template<int dim> typename hp::DoFHandler<dim>::active_cell_iterator AsFunction<dim>::get_end() const{
  return dh->end();
}

template class AsFunction<DIM>;
