vector_tools.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2016 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/dofs/function_map.h>
#include <deal.II/fe/mapping_q.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, typename Number> struct FunctionMap;
template <int dim> class Quadrature;
template <int dim> class QGauss;

template <typename number> class Vector;
template <typename number> class FullMatrix;
template <int dim, int spacedim> class Mapping;
template <int dim, int spacedim> class DoFHandler;
template <typename gridtype> class InterGridMap;
namespace hp
{
  template <int dim, int spacedim> class DoFHandler;
  template <int dim, int spacedim> class MappingCollection;
  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 <typename VectorType, int dim, int spacedim, template <int, int> class DoFHandlerType>
  void interpolate (const Mapping<dim,spacedim>        &mapping,
                    const DoFHandlerType<dim,spacedim> &dof,
                    const Function<spacedim,double>    &function,
                    VectorType                         &vec);

  template <typename VectorType, typename DoFHandlerType>
  void interpolate (const DoFHandlerType                                   &dof,
                    const Function<DoFHandlerType::space_dimension,double> &function,
                    VectorType                                             &vec);

  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 <typename VectorType, typename DoFHandlerType>
  void
  interpolate_based_on_material_id
  (const Mapping<DoFHandlerType::dimension, DoFHandlerType::space_dimension> &mapping,
   const DoFHandlerType                                                  &dof_handler,
   const std::map<types::material_id, const Function<DoFHandlerType::space_dimension, double> *> &function_map,
   VectorType                                                            &dst,
   const ComponentMask                                                   &component_mask = ComponentMask());

  template <int dim, int spacedim,
            template <int, int> class DoFHandlerType,
            typename VectorType>
  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,
            template <int, int> class DoFHandlerType,
            typename VectorType>
  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,
            template <int, int> class DoFHandlerType,
            typename VectorType>
  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,double> &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,double> &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,double>            &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,double>    &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 <typename DoFHandlerType>
  void
  interpolate_boundary_values
  (const Mapping<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &mapping,
   const DoFHandlerType                                                     &dof,
   const typename FunctionMap<DoFHandlerType::space_dimension>::type        &function_map,
   std::map<types::global_dof_index,double>                                 &boundary_values,
   const ComponentMask                                                      &component_mask = ComponentMask());

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

  template <typename DoFHandlerType>
  void
  interpolate_boundary_values
  (const Mapping<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &mapping,
   const DoFHandlerType                                                     &dof,
   const types::boundary_id                                                  boundary_component,
   const Function<DoFHandlerType::space_dimension,double>                   &boundary_function,
   std::map<types::global_dof_index,double>                                 &boundary_values,
   const ComponentMask                                                      &component_mask = ComponentMask());

  template <typename DoFHandlerType>
  void
  interpolate_boundary_values
  (const DoFHandlerType                                   &dof,
   const types::boundary_id                                boundary_component,
   const Function<DoFHandlerType::space_dimension,double> &boundary_function,
   std::map<types::global_dof_index,double>               &boundary_values,
   const ComponentMask                                    &component_mask = ComponentMask());


  template <typename DoFHandlerType>
  void
  interpolate_boundary_values
  (const DoFHandlerType                                              &dof,
   const typename FunctionMap<DoFHandlerType::space_dimension>::type &function_map,
   std::map<types::global_dof_index,double>                          &boundary_values,
   const ComponentMask                                               &component_mask = ComponentMask());


  template <typename DoFHandlerType>
  void
  interpolate_boundary_values
  (const Mapping<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &mapping,
   const DoFHandlerType                                                     &dof,
   const typename FunctionMap<DoFHandlerType::space_dimension>::type        &function_map,
   ConstraintMatrix                                                         &constraints,
   const ComponentMask                                                      &component_mask = ComponentMask());

  template <typename DoFHandlerType>
  void
  interpolate_boundary_values
  (const Mapping<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &mapping,
   const DoFHandlerType                                                     &dof,
   const types::boundary_id                                                  boundary_component,
   const Function<DoFHandlerType::space_dimension,double>                   &boundary_function,
   ConstraintMatrix                                                         &constraints,
   const ComponentMask                                                      &component_mask = ComponentMask());

  template <typename DoFHandlerType>
  void
  interpolate_boundary_values
  (const DoFHandlerType                                   &dof,
   const types::boundary_id                                boundary_component,
   const Function<DoFHandlerType::space_dimension,double> &boundary_function,
   ConstraintMatrix                                       &constraints,
   const ComponentMask                                    &component_mask = ComponentMask());


  template <typename DoFHandlerType>
  void
  interpolate_boundary_values
  (const DoFHandlerType                                              &dof,
   const typename FunctionMap<DoFHandlerType::space_dimension>::type &function_map,
   ConstraintMatrix                                                  &constraints,
   const ComponentMask                                               &component_mask = ComponentMask());


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

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

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

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

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

  template <int dim, int spacedim>
  void project_boundary_values (const DoFHandler<dim,spacedim> &dof,
                                const typename FunctionMap<spacedim>::type &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, template <int, int> class DoFHandlerType, int spacedim>
  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, template <int, int> class DoFHandlerType, int spacedim>
  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, template <int, int> class DoFHandlerType, int spacedim>
  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, template <int, int> class DoFHandlerType, int spacedim>
  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>
  void create_right_hand_side (const Mapping<dim, spacedim>    &mapping,
                               const DoFHandler<dim,spacedim> &dof,
                               const Quadrature<dim> &q,
                               const Function<spacedim,double>   &rhs,
                               Vector<double>        &rhs_vector);

  template <int dim, int spacedim>
  void create_right_hand_side (const DoFHandler<dim,spacedim> &dof,
                               const Quadrature<dim> &q,
                               const Function<spacedim,double>   &rhs,
                               Vector<double>        &rhs_vector);

  template <int dim, int spacedim>
  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,double>   &rhs,
                               Vector<double>        &rhs_vector);

  template <int dim, int spacedim>
  void create_right_hand_side (const hp::DoFHandler<dim,spacedim> &dof,
                               const hp::QCollection<dim> &q,
                               const Function<spacedim,double>   &rhs,
                               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,
                                  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>
  void create_boundary_right_hand_side (const Mapping<dim,spacedim>      &mapping,
                                        const DoFHandler<dim,spacedim>   &dof,
                                        const Quadrature<dim-1> &q,
                                        const Function<spacedim,double>     &rhs,
                                        Vector<double>          &rhs_vector,
                                        const std::set<types::boundary_id> &boundary_ids = std::set<types::boundary_id>());

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

  template <int dim, int spacedim>
  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,double>     &rhs,
                                        Vector<double>          &rhs_vector,
                                        const std::set<types::boundary_id> &boundary_ids = std::set<types::boundary_id>());

  template <int dim, int spacedim>
  void create_boundary_right_hand_side (const hp::DoFHandler<dim,spacedim>   &dof,
                                        const hp::QCollection<dim-1> &q,
                                        const Function<spacedim,double>     &rhs,
                                        Vector<double>          &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, typename VectorType, int spacedim>
  void point_difference (const DoFHandler<dim,spacedim>  &dof,
                         const VectorType                &fe_function,
                         const Function<spacedim,double> &exact_solution,
                         Vector<double>                  &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,double> &exact_solution,
                         Vector<double>                  &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<double>                 &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<double>                     &value);

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

  template <int dim, typename VectorType, int spacedim>
  double
  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<double>                 &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<double>                             &value);

  template <int dim, typename VectorType, int spacedim>
  double
  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>
  double
  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>
  double 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>
  double compute_mean_value (const DoFHandler<dim,spacedim> &dof,
                             const Quadrature<dim>          &quadrature,
                             const VectorType               &v,
                             const unsigned int             component);

  template<typename DoFHandlerType, typename VectorType>
  void get_position_vector(const DoFHandlerType &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.");

  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