vector_tools.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__vector_tools_h
#define dealii__vector_tools_h


#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/point.h>
#include <deal.II/base/std_cxx11/function.h>
#include <deal.II/dofs/function_map.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/mapping_collection.h>

#include <map>
#include <vector>
#include <set>

DEAL_II_NAMESPACE_OPEN

template <int dim, typename Number> class Function;
template <int dim> class Quadrature;
template <int dim> class QGauss;
template <int dim, typename number> class MatrixFree;

template <typename number> class Vector;
template <typename number> class FullMatrix;
template <int dim, int spacedim> class Mapping;
template <typename gridtype> class InterGridMap;
namespace hp
{
  template <int dim> class QCollection;
}
class ConstraintMatrix;


//TODO: Move documentation of functions to the functions!

namespace VectorTools
{
  enum NormType
  {
    mean,

    L1_norm,

    L2_norm,

    Lp_norm,

    Linfty_norm,

    H1_seminorm,

    Hdiv_seminorm,

    H1_norm,

    W1p_seminorm,

    W1p_norm,

    W1infty_seminorm,

    W1infty_norm

  };
  template <int dim, int spacedim, typename VectorType,
            template <int, int> class DoFHandlerType>
  void interpolate (const Mapping<dim,spacedim>        &mapping,
                    const DoFHandlerType<dim,spacedim> &dof,
                    const Function<spacedim,typename VectorType::value_type> &function,
                    VectorType                         &vec,
                    const ComponentMask                &component_mask = ComponentMask());

  template <int dim, int spacedim, typename VectorType,
            template <int, int> class DoFHandlerType>
  void interpolate
  (const DoFHandlerType<dim,spacedim>                       &dof,
   const Function<spacedim,typename VectorType::value_type> &function,
   VectorType                                               &vec,
   const ComponentMask                                      &component_mask = ComponentMask());

  template <int dim, class InVector, class OutVector, int spacedim>
  void interpolate (const DoFHandler<dim,spacedim> &dof_1,
                    const DoFHandler<dim,spacedim> &dof_2,
                    const FullMatrix<double>       &transfer,
                    const InVector                 &data_1,
                    OutVector                      &data_2);

  template <int dim, int spacedim, typename VectorType,
            template <int,int> class DoFHandlerType>
  void
  interpolate_based_on_material_id
  (const Mapping<dim,spacedim>        &mapping,
   const DoFHandlerType<dim,spacedim> &dof_handler,
   const std::map<types::material_id,
   const Function<spacedim, typename VectorType::value_type> *> &function_map,
   VectorType                         &dst,
   const ComponentMask                &component_mask = ComponentMask());

  template <int dim, int spacedim, typename VectorType,
            template <int, int> class DoFHandlerType>
  void
  interpolate_to_different_mesh (const DoFHandlerType<dim, spacedim> &dof1,
                                 const VectorType                    &u1,
                                 const DoFHandlerType<dim, spacedim> &dof2,
                                 VectorType                          &u2);

  template <int dim, int spacedim, typename VectorType,
            template <int, int> class DoFHandlerType>
  void
  interpolate_to_different_mesh (const DoFHandlerType<dim, spacedim> &dof1,
                                 const VectorType                    &u1,
                                 const DoFHandlerType<dim, spacedim> &dof2,
                                 const ConstraintMatrix              &constraints,
                                 VectorType                          &u2);

  template <int dim, int spacedim, typename VectorType,
            template <int, int> class DoFHandlerType>
  void
  interpolate_to_different_mesh
  (const InterGridMap<DoFHandlerType<dim, spacedim> > &intergridmap,
   const VectorType                                   &u1,
   const ConstraintMatrix                             &constraints,
   VectorType                                         &u2);

  template <int dim, typename VectorType, int spacedim>
  void project (const Mapping<dim, spacedim>    &mapping,
                const DoFHandler<dim,spacedim>  &dof,
                const ConstraintMatrix          &constraints,
                const Quadrature<dim>           &quadrature,
                const Function<spacedim,typename VectorType::value_type> &function,
                VectorType                      &vec,
                const bool                      enforce_zero_boundary = false,
                const Quadrature<dim-1>         &q_boundary = (dim > 1 ?
                                                              QGauss<dim-1>(2) :
                                                              Quadrature<dim-1>(0)),
                const bool                      project_to_boundary_first = false);

  template <int dim, typename VectorType, int spacedim>
  void project (const DoFHandler<dim,spacedim>  &dof,
                const ConstraintMatrix          &constraints,
                const Quadrature<dim>           &quadrature,
                const Function<spacedim,typename VectorType::value_type> &function,
                VectorType                      &vec,
                const bool                       enforce_zero_boundary = false,
                const Quadrature<dim-1>         &q_boundary = (dim > 1 ?
                                                              QGauss<dim-1>(2) :
                                                              Quadrature<dim-1>(0)),
                const bool                       project_to_boundary_first = false);

  template <int dim, typename VectorType, int spacedim>
  void project (const hp::MappingCollection<dim, spacedim> &mapping,
                const hp::DoFHandler<dim,spacedim>         &dof,
                const ConstraintMatrix                     &constraints,
                const hp::QCollection<dim>                 &quadrature,
                const Function<spacedim,typename VectorType::value_type>            &function,
                VectorType                                 &vec,
                const bool                                  enforce_zero_boundary = false,
                const hp::QCollection<dim-1> &q_boundary = hp::QCollection<dim-1>(dim > 1 ?
                                                           QGauss<dim-1>(2) :
                                                           Quadrature<dim-1>(0)),
                const bool                                  project_to_boundary_first = false);

  template <int dim, typename VectorType, int spacedim>
  void project (const hp::DoFHandler<dim,spacedim> &dof,
                const ConstraintMatrix             &constraints,
                const hp::QCollection<dim>         &quadrature,
                const Function<spacedim,typename VectorType::value_type> &function,
                VectorType                         &vec,
                const bool                          enforce_zero_boundary = false,
                const hp::QCollection<dim-1>       &q_boundary = hp::QCollection<dim-1>(dim > 1 ?
                    QGauss<dim-1>(2) :
                    Quadrature<dim-1>(0)),
                const bool                          project_to_boundary_first = false);

  template <int dim, typename VectorType, int spacedim>
  void project (const Mapping<dim, spacedim>   &mapping,
                const DoFHandler<dim,spacedim> &dof,
                const ConstraintMatrix         &constraints,
                const Quadrature<dim>          &quadrature,
                const std_cxx11::function< typename VectorType::value_type (const typename DoFHandler<dim, spacedim>::active_cell_iterator &, const unsigned int)> func,
                VectorType                     &vec_result);

  template <int dim, typename VectorType>
  void project (std_cxx11::shared_ptr<const MatrixFree<dim,typename VectorType::value_type> > data,
                const ConstraintMatrix &constraints,
                const unsigned int n_q_points_1d,
                const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> func,
                VectorType &vec_result);

  template <int dim, typename VectorType>
  void project (std_cxx11::shared_ptr<const MatrixFree<dim,typename VectorType::value_type> > data,
                const ConstraintMatrix &constraints,
                const std_cxx11::function< VectorizedArray<typename VectorType::value_type> (const unsigned int, const unsigned int)> func,
                VectorType &vec_result);

  template <int dim, int spacedim, template <int,int> class DoFHandlerType,
            typename number>
  void
  interpolate_boundary_values
  (const Mapping<dim,spacedim>                                          &mapping,
   const DoFHandlerType<dim,spacedim>                                   &dof,
   const std::map<types::boundary_id, const Function<spacedim,number>*> &function_map,
   std::map<types::global_dof_index,number>                             &boundary_values,
   const ComponentMask &component_mask = ComponentMask());

  template <int dim, int spacedim, typename number>
  void
  interpolate_boundary_values
  (const hp::MappingCollection<dim,spacedim>                            &mapping,
   const hp::DoFHandler<dim,spacedim>                                   &dof,
   const std::map<types::boundary_id, const Function<spacedim,number>*> &function_map,
   std::map<types::global_dof_index,number>                             &boundary_values,
   const ComponentMask &component_mask = ComponentMask());

  template <int dim, int spacedim, template <int,int> class DoFHandlerType,
            typename number>
  void
  interpolate_boundary_values
  (const Mapping<dim,spacedim>              &mapping,
   const DoFHandlerType<dim,spacedim>       &dof,
   const types::boundary_id                  boundary_component,
   const Function<spacedim,number>          &boundary_function,
   std::map<types::global_dof_index,number> &boundary_values,
   const ComponentMask                      &component_mask = ComponentMask());

  template <int dim, int spacedim, template <int,int> class DoFHandlerType,
            typename number>
  void
  interpolate_boundary_values
  (const DoFHandlerType<dim,spacedim>       &dof,
   const types::boundary_id                  boundary_component,
   const Function<spacedim,number>          &boundary_function,
   std::map<types::global_dof_index,number> &boundary_values,
   const ComponentMask                      &component_mask = ComponentMask());


  template <int dim, int spacedim, template <int,int> class DoFHandlerType,
            typename number>
  void
  interpolate_boundary_values
  (const DoFHandlerType<dim,spacedim>                                   &dof,
   const std::map<types::boundary_id, const Function<spacedim,number>*> &function_map,
   std::map<types::global_dof_index,number>                             &boundary_values,
   const ComponentMask &component_mask = ComponentMask());


  template <int dim, int spacedim, template <int,int> class DoFHandlerType,
            typename number>
  void
  interpolate_boundary_values
  (const Mapping<dim,spacedim>                                          &mapping,
   const DoFHandlerType<dim,spacedim>                                   &dof,
   const std::map<types::boundary_id, const Function<spacedim,number>*> &function_map,
   ConstraintMatrix                                                     &constraints,
   const ComponentMask &component_mask = ComponentMask());

  template <int dim, int spacedim, template <int,int> class DoFHandlerType,
            typename number>
  void
  interpolate_boundary_values
  (const Mapping<dim,spacedim>        &mapping,
   const DoFHandlerType<dim,spacedim> &dof,
   const types::boundary_id            boundary_component,
   const Function<spacedim,number>    &boundary_function,
   ConstraintMatrix                   &constraints,
   const ComponentMask                &component_mask = ComponentMask());

  template <int dim, int spacedim, template <int,int> class DoFHandlerType,
            typename number>
  void
  interpolate_boundary_values
  (const DoFHandlerType<dim,spacedim> &dof,
   const types::boundary_id             boundary_component,
   const Function<spacedim,number>     &boundary_function,
   ConstraintMatrix                    &constraints,
   const ComponentMask                 &component_mask = ComponentMask());


  template <int dim, int spacedim, template <int,int> class DoFHandlerType,
            typename number>
  void
  interpolate_boundary_values
  (const DoFHandlerType<dim,spacedim>                                   &dof,
   const std::map<types::boundary_id, const Function<spacedim,number>*> &function_map,
   ConstraintMatrix                                                     &constraints,
   const ComponentMask &component_mask = ComponentMask());


  template <int dim, int spacedim, typename number>
  void project_boundary_values
  (const Mapping<dim, spacedim>                                         &mapping,
   const DoFHandler<dim,spacedim>                                       &dof,
   const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_functions,
   const Quadrature<dim-1>                                              &q,
   std::map<types::global_dof_index,number>                             &boundary_values,
   std::vector<unsigned int> component_mapping = std::vector<unsigned int>());

  template <int dim, int spacedim, typename number>
  void project_boundary_values
  (const DoFHandler<dim,spacedim>                                       &dof,
   const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_function,
   const Quadrature<dim-1>                                              &q,
   std::map<types::global_dof_index,number>                             &boundary_values,
   std::vector<unsigned int> component_mapping = std::vector<unsigned int>());

  template <int dim, int spacedim, typename number>
  void project_boundary_values
  (const hp::MappingCollection<dim, spacedim>                           &mapping,
   const hp::DoFHandler<dim,spacedim>                                   &dof,
   const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_functions,
   const hp::QCollection<dim-1>                                         &q,
   std::map<types::global_dof_index,number>                             &boundary_values,
   std::vector<unsigned int> component_mapping = std::vector<unsigned int>());

  template <int dim, int spacedim, typename number>
  void project_boundary_values
  (const hp::DoFHandler<dim,spacedim>                                   &dof,
   const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_function,
   const hp::QCollection<dim-1>                                         &q,
   std::map<types::global_dof_index,number>                             &boundary_values,
   std::vector<unsigned int> component_mapping = std::vector<unsigned int>());

  template <int dim, int spacedim, typename number>
  void project_boundary_values
  (const Mapping<dim, spacedim>                                         &mapping,
   const DoFHandler<dim,spacedim>                                       &dof,
   const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_functions,
   const Quadrature<dim-1>                                              &q,
   ConstraintMatrix                                                     &constraints,
   std::vector<unsigned int> component_mapping = std::vector<unsigned int>());

  template <int dim, int spacedim, typename number>
  void project_boundary_values
  (const DoFHandler<dim,spacedim>                                       &dof,
   const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_function,
   const Quadrature<dim-1>                                              &q,
   ConstraintMatrix                                                     &constraints,
   std::vector<unsigned int> component_mapping = std::vector<unsigned int>());


  template <int dim>
  void project_boundary_values_curl_conforming
  (const DoFHandler<dim>      &dof_handler,
   const unsigned int          first_vector_component,
   const Function<dim,double> &boundary_function,
   const types::boundary_id    boundary_component,
   ConstraintMatrix           &constraints,
   const Mapping<dim>         &mapping = StaticMappingQ1<dim>::mapping);

  template <int dim>
  void project_boundary_values_curl_conforming
  (const hp::DoFHandler<dim>             &dof_handler,
   const unsigned int                     first_vector_component,
   const Function<dim,double>            &boundary_function,
   const types::boundary_id               boundary_component,
   ConstraintMatrix                      &constraints,
   const hp::MappingCollection<dim, dim> &mapping_collection = hp::StaticMappingQ1<dim>::mapping_collection);

  template <int dim>
  void project_boundary_values_curl_conforming_l2
  (const DoFHandler<dim>      &dof_handler,
   const unsigned int          first_vector_component,
   const Function<dim,double> &boundary_function,
   const types::boundary_id    boundary_component,
   ConstraintMatrix           &constraints,
   const Mapping<dim>         &mapping = StaticMappingQ1<dim>::mapping);


  template <int dim>
  void project_boundary_values_curl_conforming_l2
  (const hp::DoFHandler<dim>  &dof_handler,
   const unsigned int          first_vector_component,
   const Function<dim,double> &boundary_function,
   const types::boundary_id    boundary_component,
   ConstraintMatrix           &constraints,
   const hp::MappingCollection<dim, dim> &mapping_collection = hp::StaticMappingQ1<dim>::mapping_collection);


  template<int dim>
  void project_boundary_values_div_conforming
  (const DoFHandler<dim>      &dof_handler,
   const unsigned int          first_vector_component,
   const Function<dim,double> &boundary_function,
   const types::boundary_id    boundary_component,
   ConstraintMatrix           &constraints,
   const Mapping<dim>         &mapping = StaticMappingQ1<dim>::mapping);

  template<int dim>
  void project_boundary_values_div_conforming
  (const hp::DoFHandler<dim>             &dof_handler,
   const unsigned int                     first_vector_component,
   const Function<dim,double>            &boundary_function,
   const types::boundary_id               boundary_component,
   ConstraintMatrix                      &constraints,
   const hp::MappingCollection<dim, dim> &mapping_collection
   = hp::StaticMappingQ1<dim>::mapping_collection);


  template <int dim, int spacedim, template <int, int> class DoFHandlerType>
  void
  compute_nonzero_normal_flux_constraints
  (const DoFHandlerType<dim,spacedim>   &dof_handler,
   const unsigned int                    first_vector_component,
   const std::set<types::boundary_id>   &boundary_ids,
   typename FunctionMap<spacedim>::type &function_map,
   ConstraintMatrix                     &constraints,
   const Mapping<dim, spacedim>         &mapping = StaticMappingQ1<dim>::mapping);

  template <int dim, int spacedim, template <int, int> class DoFHandlerType>
  void
  compute_no_normal_flux_constraints
  (const DoFHandlerType<dim,spacedim> &dof_handler,
   const unsigned int                  first_vector_component,
   const std::set<types::boundary_id> &boundary_ids,
   ConstraintMatrix                   &constraints,
   const Mapping<dim, spacedim>       &mapping = StaticMappingQ1<dim>::mapping);

  template <int dim, int spacedim, template <int, int> class DoFHandlerType>
  void
  compute_nonzero_tangential_flux_constraints
  (const DoFHandlerType<dim,spacedim>   &dof_handler,
   const unsigned int                    first_vector_component,
   const std::set<types::boundary_id>   &boundary_ids,
   typename FunctionMap<spacedim>::type &function_map,
   ConstraintMatrix                     &constraints,
   const Mapping<dim, spacedim>         &mapping = StaticMappingQ1<dim>::mapping);

  template <int dim, int spacedim, template <int, int> class DoFHandlerType>
  void
  compute_normal_flux_constraints
  (const DoFHandlerType<dim,spacedim> &dof_handler,
   const unsigned int                  first_vector_component,
   const std::set<types::boundary_id> &boundary_ids,
   ConstraintMatrix                   &constraints,
   const Mapping<dim, spacedim>       &mapping = StaticMappingQ1<dim>::mapping);




  template <int dim, int spacedim, typename VectorType>
  void create_right_hand_side (const Mapping<dim, spacedim>    &mapping,
                               const DoFHandler<dim,spacedim>  &dof,
                               const Quadrature<dim>           &q,
                               const Function<spacedim,typename VectorType::value_type> &rhs,
                               VectorType                      &rhs_vector,
                               const ConstraintMatrix          &constraints=ConstraintMatrix());

  template <int dim, int spacedim, typename VectorType>
  void create_right_hand_side (const DoFHandler<dim,spacedim>  &dof,
                               const Quadrature<dim>           &q,
                               const Function<spacedim,typename VectorType::value_type> &rhs,
                               VectorType                      &rhs_vector,
                               const ConstraintMatrix          &constraints=ConstraintMatrix());

  template <int dim, int spacedim, typename VectorType>
  void create_right_hand_side (const hp::MappingCollection<dim,spacedim> &mapping,
                               const hp::DoFHandler<dim,spacedim>        &dof,
                               const hp::QCollection<dim>                &q,
                               const Function<spacedim,typename VectorType::value_type> &rhs,
                               VectorType                                &rhs_vector,
                               const ConstraintMatrix                    &constraints=ConstraintMatrix());

  template <int dim, int spacedim, typename VectorType>
  void create_right_hand_side (const hp::DoFHandler<dim,spacedim> &dof,
                               const hp::QCollection<dim>         &q,
                               const Function<spacedim,typename VectorType::value_type> &rhs,
                               VectorType                         &rhs_vector,
                               const ConstraintMatrix             &constraints=ConstraintMatrix());

  template <int dim, int spacedim>
  void create_point_source_vector(const Mapping<dim,spacedim>    &mapping,
                                  const DoFHandler<dim,spacedim> &dof,
                                  const Point<spacedim>          &p,
                                  Vector<double>                 &rhs_vector);

  template <int dim, int spacedim>
  void create_point_source_vector(const DoFHandler<dim,spacedim> &dof,
                                  const Point<spacedim>          &p,
                                  Vector<double>                 &rhs_vector);

  template <int dim, int spacedim>
  void create_point_source_vector(const hp::MappingCollection<dim,spacedim> &mapping,
                                  const hp::DoFHandler<dim,spacedim>        &dof,
                                  const Point<spacedim>                     &p,
                                  Vector<double>                            &rhs_vector);

  template <int dim, int spacedim>
  void create_point_source_vector(const hp::DoFHandler<dim,spacedim> &dof,
                                  const Point<spacedim>              &p,
                                  Vector<double>                     &rhs_vector);

  template <int dim, int spacedim>
  void create_point_source_vector(const Mapping<dim,spacedim>    &mapping,
                                  const DoFHandler<dim,spacedim> &dof,
                                  const Point<spacedim>          &p,
                                  const Point<dim>               &direction,
                                  Vector<double>                 &rhs_vector);

  template <int dim, int spacedim>
  void create_point_source_vector(const DoFHandler<dim,spacedim> &dof,
                                  const Point<spacedim>          &p,
                                  const Point<dim>               &direction,
                                  Vector<double>                 &rhs_vector);

  template <int dim, int spacedim>
  void create_point_source_vector(const hp::MappingCollection<dim,spacedim> &mapping,
                                  const hp::DoFHandler<dim,spacedim>        &dof,
                                  const Point<spacedim>                     &p,
                                  const Point<dim>                          &direction,
                                  Vector<double>                            &rhs_vector);

  template <int dim, int spacedim>
  void create_point_source_vector(const hp::DoFHandler<dim,spacedim> &dof,
                                  const Point<spacedim>              &p,
                                  const Point<dim>                   &direction,
                                  Vector<double>                     &rhs_vector);

  template <int dim, int spacedim, typename VectorType>
  void create_boundary_right_hand_side (const Mapping<dim,spacedim>        &mapping,
                                        const DoFHandler<dim,spacedim>     &dof,
                                        const Quadrature<dim-1>            &q,
                                        const Function<spacedim,typename VectorType::value_type> &rhs,
                                        VectorType                         &rhs_vector,
                                        const std::set<types::boundary_id> &boundary_ids = std::set<types::boundary_id>());

  template <int dim, int spacedim, typename VectorType>
  void create_boundary_right_hand_side
  (const DoFHandler<dim,spacedim>     &dof,
   const Quadrature<dim-1>            &q,
   const Function<spacedim,typename VectorType::value_type> &rhs,
   VectorType                         &rhs_vector,
   const std::set<types::boundary_id> &boundary_ids = std::set<types::boundary_id>());

  template <int dim, int spacedim, typename VectorType>
  void create_boundary_right_hand_side
  (const hp::MappingCollection<dim,spacedim> &mapping,
   const hp::DoFHandler<dim,spacedim>        &dof,
   const hp::QCollection<dim-1>              &q,
   const Function<spacedim,typename VectorType::value_type> &rhs,
   VectorType                                &rhs_vector,
   const std::set<types::boundary_id>        &boundary_ids = std::set<types::boundary_id>());

  template <int dim, int spacedim, typename VectorType>
  void create_boundary_right_hand_side
  (const hp::DoFHandler<dim,spacedim> &dof,
   const hp::QCollection<dim-1>       &q,
   const Function<spacedim,typename VectorType::value_type> &rhs,
   VectorType                         &rhs_vector,
   const std::set<types::boundary_id> &boundary_ids = std::set<types::boundary_id>());



  template <int dim, class InVector, class OutVector, int spacedim>
  void integrate_difference (const Mapping<dim,spacedim>     &mapping,
                             const DoFHandler<dim,spacedim>  &dof,
                             const InVector                  &fe_function,
                             const Function<spacedim,double> &exact_solution,
                             OutVector                       &difference,
                             const Quadrature<dim>           &q,
                             const NormType                  &norm,
                             const Function<spacedim,double> *weight   = 0,
                             const double                     exponent = 2.);

  template <int dim, class InVector, class OutVector, int spacedim>
  void integrate_difference (const DoFHandler<dim,spacedim>  &dof,
                             const InVector                  &fe_function,
                             const Function<spacedim,double> &exact_solution,
                             OutVector                       &difference,
                             const Quadrature<dim>           &q,
                             const NormType                  &norm,
                             const Function<spacedim,double> *weight   = 0,
                             const double                     exponent = 2.);

  template <int dim, class InVector, class OutVector, int spacedim>
  void integrate_difference (const hp::MappingCollection<dim,spacedim> &mapping,
                             const hp::DoFHandler<dim,spacedim>        &dof,
                             const InVector                            &fe_function,
                             const Function<spacedim,double>           &exact_solution,
                             OutVector                                 &difference,
                             const hp::QCollection<dim>                &q,
                             const NormType                            &norm,
                             const Function<spacedim,double>           *weight   = 0,
                             const double                               exponent = 2.);

  template <int dim, class InVector, class OutVector, int spacedim>
  void integrate_difference (const hp::DoFHandler<dim,spacedim> &dof,
                             const InVector                     &fe_function,
                             const Function<spacedim,double>    &exact_solution,
                             OutVector                          &difference,
                             const hp::QCollection<dim>         &q,
                             const NormType                     &norm,
                             const Function<spacedim,double>    *weight   = 0,
                             const double                        exponent = 2.);

  template <int dim, int spacedim, class InVector>
  double compute_global_error(const Triangulation<dim,spacedim> &tria,
                              const InVector                    &cellwise_error,
                              const NormType                    &norm,
                              const double                       exponent = 2.);

  template <int dim, typename VectorType, int spacedim>
  void point_difference
  (const DoFHandler<dim,spacedim>                           &dof,
   const VectorType                                         &fe_function,
   const Function<spacedim,typename VectorType::value_type> &exact_solution,
   Vector<typename VectorType::value_type>                  &difference,
   const Point<spacedim>                                    &point);

  template <int dim, typename VectorType, int spacedim>
  void point_difference
  (const Mapping<dim, spacedim>                             &mapping,
   const DoFHandler<dim,spacedim>                           &dof,
   const VectorType                                         &fe_function,
   const Function<spacedim,typename VectorType::value_type> &exact_solution,
   Vector<typename VectorType::value_type>                  &difference,
   const Point<spacedim>                                    &point);

  template <int dim, typename VectorType, int spacedim>
  void
  point_value (const DoFHandler<dim,spacedim>          &dof,
               const VectorType                        &fe_function,
               const Point<spacedim>                   &point,
               Vector<typename VectorType::value_type> &value);

  template <int dim, typename VectorType, int spacedim>
  void
  point_value (const hp::DoFHandler<dim,spacedim>      &dof,
               const VectorType                        &fe_function,
               const Point<spacedim>                   &point,
               Vector<typename VectorType::value_type> &value);

  template <int dim, typename VectorType, int spacedim>
  typename VectorType::value_type
  point_value (const DoFHandler<dim,spacedim> &dof,
               const VectorType               &fe_function,
               const Point<spacedim>          &point);

  template <int dim, typename VectorType, int spacedim>
  typename VectorType::value_type
  point_value (const hp::DoFHandler<dim,spacedim> &dof,
               const VectorType                   &fe_function,
               const Point<spacedim>              &point);

  template <int dim, typename VectorType, int spacedim>
  void
  point_value (const Mapping<dim, spacedim>            &mapping,
               const DoFHandler<dim,spacedim>          &dof,
               const VectorType                        &fe_function,
               const Point<spacedim>                   &point,
               Vector<typename VectorType::value_type> &value);

  template <int dim, typename VectorType, int spacedim>
  void
  point_value (const hp::MappingCollection<dim, spacedim> &mapping,
               const hp::DoFHandler<dim,spacedim>         &dof,
               const VectorType                           &fe_function,
               const Point<spacedim>                      &point,
               Vector<typename VectorType::value_type>    &value);

  template <int dim, typename VectorType, int spacedim>
  typename VectorType::value_type
  point_value (const Mapping<dim,spacedim>    &mapping,
               const DoFHandler<dim,spacedim> &dof,
               const VectorType               &fe_function,
               const Point<spacedim>          &point);

  template <int dim, typename VectorType, int spacedim>
  typename VectorType::value_type
  point_value (const hp::MappingCollection<dim,spacedim> &mapping,
               const hp::DoFHandler<dim,spacedim>        &dof,
               const VectorType                          &fe_function,
               const Point<spacedim>                     &point);

  template <int dim, typename VectorType, int spacedim>
  void
  point_gradient (const DoFHandler<dim,spacedim> &dof,
                  const VectorType               &fe_function,
                  const Point<spacedim>          &point,
                  std::vector<Tensor<1, spacedim, typename VectorType::value_type> > &value);

  template <int dim, typename VectorType, int spacedim>
  void
  point_gradient (const hp::DoFHandler<dim,spacedim> &dof,
                  const VectorType                   &fe_function,
                  const Point<spacedim>              &point,
                  std::vector<Tensor<1, spacedim, typename VectorType::value_type> > &value);

  template <int dim, typename VectorType, int spacedim>
  Tensor<1, spacedim, typename VectorType::value_type>
  point_gradient (const DoFHandler<dim,spacedim> &dof,
                  const VectorType               &fe_function,
                  const Point<spacedim>          &point);

  template <int dim, typename VectorType, int spacedim>
  Tensor<1, spacedim, typename VectorType::value_type>
  point_gradient (const hp::DoFHandler<dim,spacedim> &dof,
                  const VectorType                   &fe_function,
                  const Point<spacedim>              &point);

  template <int dim, typename VectorType, int spacedim>
  void
  point_gradient (const Mapping<dim, spacedim>   &mapping,
                  const DoFHandler<dim,spacedim> &dof,
                  const VectorType               &fe_function,
                  const Point<spacedim>          &point,
                  std::vector<Tensor<1, spacedim, typename VectorType::value_type> > &value);

  template <int dim, typename VectorType, int spacedim>
  void
  point_gradient (const hp::MappingCollection<dim, spacedim> &mapping,
                  const hp::DoFHandler<dim,spacedim>         &dof,
                  const VectorType                           &fe_function,
                  const Point<spacedim>                      &point,
                  std::vector<Tensor<1, spacedim, typename VectorType::value_type> > &value);

  template <int dim, typename VectorType, int spacedim>
  Tensor<1, spacedim, typename VectorType::value_type>
  point_gradient (const Mapping<dim,spacedim>    &mapping,
                  const DoFHandler<dim,spacedim> &dof,
                  const VectorType               &fe_function,
                  const Point<spacedim>          &point);

  template <int dim, typename VectorType, int spacedim>
  Tensor<1, spacedim, typename VectorType::value_type>
  point_gradient (const hp::MappingCollection<dim,spacedim> &mapping,
                  const hp::DoFHandler<dim,spacedim>        &dof,
                  const VectorType                          &fe_function,
                  const Point<spacedim>                     &point);



  template <typename VectorType>
  void subtract_mean_value(VectorType              &v,
                           const std::vector<bool> &p_select = std::vector<bool>());


  template <int dim, typename VectorType, int spacedim>
  typename VectorType::value_type
  compute_mean_value (const Mapping<dim, spacedim>   &mapping,
                      const DoFHandler<dim,spacedim> &dof,
                      const Quadrature<dim>          &quadrature,
                      const VectorType               &v,
                      const unsigned int             component);

  template <int dim, typename VectorType, int spacedim>
  typename VectorType::value_type
  compute_mean_value (const DoFHandler<dim,spacedim> &dof,
                      const Quadrature<dim>          &quadrature,
                      const VectorType               &v,
                      const unsigned int             component);

  template <int dim, int spacedim, template <int,int> class DoFHandlerType,
            typename VectorType>
  void get_position_vector(const DoFHandlerType<dim,spacedim> &dh,
                           VectorType           &vector,
                           const ComponentMask  &mask = ComponentMask());


  DeclExceptionMsg (ExcNonInterpolatingFE,
                    "You are attempting an operation that requires the "
                    "finite element involved to be 'interpolating', i.e., "
                    "it needs to have support points. The finite element "
                    "you are using here does not appear to have those for "
                    "the required components.");

  DeclExceptionMsg (ExcPointNotAvailableHere,
                    "The given point is inside a cell of a "
                    "parallel::distributed::Triangulation that is not "
                    "locally owned by this processor.");
}


DEAL_II_NAMESPACE_CLOSE

#endif