vector_tools.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2018 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#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/function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>

#include <deal.II/dofs/deprecated_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 <functional>
#include <map>
#include <set>
#include <vector>

DEAL_II_NAMESPACE_OPEN

template <int dim, typename RangeNumberType>
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;
}
template <typename number>
class AffineConstraints;


// 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 AffineConstraints<typename VectorType::value_type> &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 AffineConstraints<typename VectorType::value_type> &constraints,
    VectorType &                                              u2);

  template <int dim, typename VectorType, int spacedim>
  void
  project(const Mapping<dim, spacedim> &                            mapping,
          const DoFHandler<dim, spacedim> &                         dof,
          const AffineConstraints<typename VectorType::value_type> &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 AffineConstraints<typename VectorType::value_type> &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 AffineConstraints<typename VectorType::value_type> &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 AffineConstraints<typename VectorType::value_type> &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 AffineConstraints<typename VectorType::value_type> &constraints,
          const Quadrature<dim> &                                   quadrature,
          const std::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::shared_ptr<
            const MatrixFree<dim, typename VectorType::value_type>> data,
          const AffineConstraints<typename VectorType::value_type> &constraints,
          const unsigned int      n_q_points_1d,
          const std::function<VectorizedArray<typename VectorType::value_type>(
            const unsigned int,
            const unsigned int)> &func,
          VectorType &            vec_result,
          const unsigned int      fe_component = 0);

  template <int dim, typename VectorType>
  void
  project(std::shared_ptr<
            const MatrixFree<dim, typename VectorType::value_type>> data,
          const AffineConstraints<typename VectorType::value_type> &constraints,
          const std::function<VectorizedArray<typename VectorType::value_type>(
            const unsigned int,
            const unsigned int)> &                                  func,
          VectorType &                                              vec_result,
          const unsigned int fe_component = 0);

  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,
    AffineConstraints<number> &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,
    AffineConstraints<number> &          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,
    AffineConstraints<number> &          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,
    AffineConstraints<number> &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,
    AffineConstraints<number> &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,
    AffineConstraints<number> &constraints,
    std::vector<unsigned int>  component_mapping = std::vector<unsigned int>());


  template <int dim>
  DEAL_II_DEPRECATED 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,
    AffineConstraints<double> &  constraints,
    const Mapping<dim> &         mapping = StaticMappingQ1<dim>::mapping);

  template <int dim>
  DEAL_II_DEPRECATED 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,
    AffineConstraints<double> &            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,
    AffineConstraints<double> &  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,
    AffineConstraints<double> &            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,
    AffineConstraints<double> &  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,
    AffineConstraints<double> &            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,
    const std::map<types::boundary_id, const Function<spacedim> *>
      &                           function_map,
    AffineConstraints<double> &   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,
    AffineConstraints<double> &          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,
    const std::map<types::boundary_id, const Function<spacedim> *>
      &                           function_map,
    AffineConstraints<double> &   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,
    AffineConstraints<double> &          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 AffineConstraints<typename VectorType::value_type> & constraints =
      AffineConstraints<typename VectorType::value_type>());

  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 AffineConstraints<typename VectorType::value_type> & constraints =
      AffineConstraints<typename VectorType::value_type>());

  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 AffineConstraints<typename VectorType::value_type> & constraints =
      AffineConstraints<typename VectorType::value_type>());

  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 AffineConstraints<typename VectorType::value_type> & constraints =
      AffineConstraints<typename VectorType::value_type>());

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

  template <int dim, int spacedim>
  void
  create_point_source_vector(const DoFHandler<dim, spacedim> &dof_handler,
                             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_handler,
    const Point<spacedim> &                     p,
    Vector<double> &                            rhs_vector);

  template <int dim, int spacedim>
  void
  create_point_source_vector(const hp::DoFHandler<dim, spacedim> &dof_handler,
                             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_handler,
                             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_handler,
                             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_handler,
    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_handler,
                             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   = nullptr,
                       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   = nullptr,
                       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   = nullptr,
                       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   = nullptr,
                       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(ExcPointNotAvailableHere,
                   "The given point is inside a cell of a "
                   "parallel::distributed::Triangulation that is not "
                   "locally owned by this processor.");
} // namespace VectorTools


DEAL_II_NAMESPACE_CLOSE

#endif