error_estimator.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2019 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_error_estimator_h
#define dealii_error_estimator_h


#include <deal.II/base/config.h>

#include <deal.II/base/exceptions.h>
#include <deal.II/base/function.h>

#include <deal.II/dofs/deprecated_function_map.h>

#include <deal.II/fe/component_mask.h>

#include <map>

DEAL_II_NAMESPACE_OPEN


template <int, int>
class Mapping;
template <int>
class Quadrature;

namespace hp
{
  template <int>
  class QCollection;
}



template <int dim, int spacedim = dim>
class KellyErrorEstimator
{
public:
  enum Strategy
  {
    cell_diameter_over_24 = 0,
    face_diameter_over_twice_max_degree,
    cell_diameter
  };

  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const Mapping<dim, spacedim> &mapping,
    const DoFHandlerType &        dof,
    const Quadrature<dim - 1> &   quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                       neumann_bc,
    const InputVector &       solution,
    Vector<float> &           error,
    const ComponentMask &     component_mask = ComponentMask(),
    const Function<spacedim> *coefficients   = nullptr,
    const unsigned int        n_threads      = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id   = numbers::invalid_subdomain_id,
    const types::material_id  material_id    = numbers::invalid_material_id,
    const Strategy            strategy       = cell_diameter_over_24);

  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const DoFHandlerType &     dof,
    const Quadrature<dim - 1> &quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                       neumann_bc,
    const InputVector &       solution,
    Vector<float> &           error,
    const ComponentMask &     component_mask = ComponentMask(),
    const Function<spacedim> *coefficients   = nullptr,
    const unsigned int        n_threads      = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id   = numbers::invalid_subdomain_id,
    const types::material_id  material_id    = numbers::invalid_material_id,
    const Strategy            strategy       = cell_diameter_over_24);

  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const Mapping<dim, spacedim> &mapping,
    const DoFHandlerType &        dof,
    const Quadrature<dim - 1> &   quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                                     neumann_bc,
    const std::vector<const InputVector *> &solutions,
    std::vector<Vector<float> *> &          errors,
    const ComponentMask &                   component_mask = ComponentMask(),
    const Function<spacedim> *              coefficients   = nullptr,
    const unsigned int        n_threads    = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
    const types::material_id  material_id  = numbers::invalid_material_id,
    const Strategy            strategy     = cell_diameter_over_24);

  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const DoFHandlerType &     dof,
    const Quadrature<dim - 1> &quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                                     neumann_bc,
    const std::vector<const InputVector *> &solutions,
    std::vector<Vector<float> *> &          errors,
    const ComponentMask &                   component_mask = ComponentMask(),
    const Function<spacedim> *              coefficients   = nullptr,
    const unsigned int        n_threads    = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
    const types::material_id  material_id  = numbers::invalid_material_id,
    const Strategy            strategy     = cell_diameter_over_24);


  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const Mapping<dim, spacedim> &  mapping,
    const DoFHandlerType &          dof,
    const hp::QCollection<dim - 1> &quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                       neumann_bc,
    const InputVector &       solution,
    Vector<float> &           error,
    const ComponentMask &     component_mask = ComponentMask(),
    const Function<spacedim> *coefficients   = nullptr,
    const unsigned int        n_threads      = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id   = numbers::invalid_subdomain_id,
    const types::material_id  material_id    = numbers::invalid_material_id,
    const Strategy            strategy       = cell_diameter_over_24);


  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const DoFHandlerType &          dof,
    const hp::QCollection<dim - 1> &quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                       neumann_bc,
    const InputVector &       solution,
    Vector<float> &           error,
    const ComponentMask &     component_mask = ComponentMask(),
    const Function<spacedim> *coefficients   = nullptr,
    const unsigned int        n_threads      = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id   = numbers::invalid_subdomain_id,
    const types::material_id  material_id    = numbers::invalid_material_id,
    const Strategy            strategy       = cell_diameter_over_24);


  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const Mapping<dim, spacedim> &  mapping,
    const DoFHandlerType &          dof,
    const hp::QCollection<dim - 1> &quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                                     neumann_bc,
    const std::vector<const InputVector *> &solutions,
    std::vector<Vector<float> *> &          errors,
    const ComponentMask &                   component_mask = ComponentMask(),
    const Function<spacedim> *              coefficients   = nullptr,
    const unsigned int        n_threads    = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
    const types::material_id  material_id  = numbers::invalid_material_id,
    const Strategy            strategy     = cell_diameter_over_24);


  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const DoFHandlerType &          dof,
    const hp::QCollection<dim - 1> &quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                                     neumann_bc,
    const std::vector<const InputVector *> &solutions,
    std::vector<Vector<float> *> &          errors,
    const ComponentMask &                   component_mask = ComponentMask(),
    const Function<spacedim> *              coefficients   = nullptr,
    const unsigned int        n_threads    = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
    const types::material_id  material_id  = numbers::invalid_material_id,
    const Strategy            strategy     = cell_diameter_over_24);

  DeclExceptionMsg(ExcInvalidComponentMask,
                   "You provided a ComponentMask argument that is invalid. "
                   "Component masks need to be either default constructed "
                   "(in which case they indicate that every component is "
                   "selected) or need to have a length equal to the number "
                   "of vector components of the finite element in use "
                   "by the DoFHandler object. In the latter case, at "
                   "least one component needs to be selected.");
  DeclExceptionMsg(
    ExcInvalidCoefficient,
    "If you do specify the argument for a (possibly "
    "spatially variable) coefficient function for this function, "
    "then it needs to refer to a coefficient that is either "
    "scalar (has one vector component) or has as many vector "
    "components as there are in the finite element used by "
    "the DoFHandler argument.");
  DeclException3(ExcInvalidBoundaryFunction,
                 types::boundary_id,
                 int,
                 int,
                 << "You provided a function map that for boundary indicator "
                 << arg1 << " specifies a function with " << arg2
                 << " vector components. However, the finite "
                    "element in use has "
                 << arg2
                 << " components, and these two numbers need to match.");
  DeclException2(ExcIncompatibleNumberOfElements,
                 int,
                 int,
                 << "The number of input vectors, " << arg1
                 << " needs to be equal to the number of output vectors, "
                 << arg2
                 << ". This is not the case in your call of this function.");
  DeclExceptionMsg(ExcNoSolutions,
                   "You need to specify at least one solution vector as "
                   "input.");
};



template <int spacedim>
class KellyErrorEstimator<1, spacedim>
{
public:
  enum Strategy
  {
    cell_diameter_over_24 = 0,
    face_diameter_over_twice_max_degree,
    cell_diameter
  };

  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const Mapping<1, spacedim> &mapping,
    const DoFHandlerType &      dof,
    const Quadrature<0> &       quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                       neumann_bc,
    const InputVector &       solution,
    Vector<float> &           error,
    const ComponentMask &     component_mask = ComponentMask(),
    const Function<spacedim> *coefficient    = nullptr,
    const unsigned int        n_threads      = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id   = numbers::invalid_subdomain_id,
    const types::material_id  material_id    = numbers::invalid_material_id,
    const Strategy            strategy       = cell_diameter_over_24);

  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const DoFHandlerType &dof,
    const Quadrature<0> & quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                       neumann_bc,
    const InputVector &       solution,
    Vector<float> &           error,
    const ComponentMask &     component_mask = ComponentMask(),
    const Function<spacedim> *coefficients   = nullptr,
    const unsigned int        n_threads      = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id   = numbers::invalid_subdomain_id,
    const types::material_id  material_id    = numbers::invalid_material_id,
    const Strategy            strategy       = cell_diameter_over_24);

  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const Mapping<1, spacedim> &mapping,
    const DoFHandlerType &      dof,
    const Quadrature<0> &       quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                                     neumann_bc,
    const std::vector<const InputVector *> &solutions,
    std::vector<Vector<float> *> &          errors,
    const ComponentMask &                   component_mask = ComponentMask(),
    const Function<spacedim> *              coefficients   = nullptr,
    const unsigned int        n_threads    = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
    const types::material_id  material_id  = numbers::invalid_material_id,
    const Strategy            strategy     = cell_diameter_over_24);

  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const DoFHandlerType &dof,
    const Quadrature<0> & quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                                     neumann_bc,
    const std::vector<const InputVector *> &solutions,
    std::vector<Vector<float> *> &          errors,
    const ComponentMask &                   component_mask = ComponentMask(),
    const Function<spacedim> *              coefficients   = nullptr,
    const unsigned int        n_threads    = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
    const types::material_id  material_id  = numbers::invalid_material_id,
    const Strategy            strategy     = cell_diameter_over_24);


  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const Mapping<1, spacedim> &mapping,
    const DoFHandlerType &      dof,
    const hp::QCollection<0> &  quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                       neumann_bc,
    const InputVector &       solution,
    Vector<float> &           error,
    const ComponentMask &     component_mask = ComponentMask(),
    const Function<spacedim> *coefficients   = nullptr,
    const unsigned int        n_threads      = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id   = numbers::invalid_subdomain_id,
    const types::material_id  material_id    = numbers::invalid_material_id,
    const Strategy            strategy       = cell_diameter_over_24);


  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const DoFHandlerType &    dof,
    const hp::QCollection<0> &quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                       neumann_bc,
    const InputVector &       solution,
    Vector<float> &           error,
    const ComponentMask &     component_mask = ComponentMask(),
    const Function<spacedim> *coefficients   = nullptr,
    const unsigned int        n_threads      = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id   = numbers::invalid_subdomain_id,
    const types::material_id  material_id    = numbers::invalid_material_id,
    const Strategy            strategy       = cell_diameter_over_24);


  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const Mapping<1, spacedim> &mapping,
    const DoFHandlerType &      dof,
    const hp::QCollection<0> &  quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                                     neumann_bc,
    const std::vector<const InputVector *> &solutions,
    std::vector<Vector<float> *> &          errors,
    const ComponentMask &                   component_mask = ComponentMask(),
    const Function<spacedim> *              coefficients   = nullptr,
    const unsigned int        n_threads    = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
    const types::material_id  material_id  = numbers::invalid_material_id,
    const Strategy            strategy     = cell_diameter_over_24);


  template <typename InputVector, typename DoFHandlerType>
  static void
  estimate(
    const DoFHandlerType &    dof,
    const hp::QCollection<0> &quadrature,
    const std::map<types::boundary_id,
                   const Function<spacedim, typename InputVector::value_type> *>
      &                                     neumann_bc,
    const std::vector<const InputVector *> &solutions,
    std::vector<Vector<float> *> &          errors,
    const ComponentMask &                   component_mask = ComponentMask(),
    const Function<spacedim> *              coefficients   = nullptr,
    const unsigned int        n_threads    = numbers::invalid_unsigned_int,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id,
    const types::material_id  material_id  = numbers::invalid_material_id,
    const Strategy            strategy     = cell_diameter_over_24);

  DeclExceptionMsg(ExcInvalidComponentMask,
                   "You provided a ComponentMask argument that is invalid. "
                   "Component masks need to be either default constructed "
                   "(in which case they indicate that every component is "
                   "selected) or need to have a length equal to the number "
                   "of vector components of the finite element in use "
                   "by the DoFHandler object. In the latter case, at "
                   "least one component needs to be selected.");
  DeclExceptionMsg(
    ExcInvalidCoefficient,
    "If you do specify the argument for a (possibly "
    "spatially variable) coefficient function for this function, "
    "then it needs to refer to a coefficient that is either "
    "scalar (has one vector component) or has as many vector "
    "components as there are in the finite element used by "
    "the DoFHandler argument.");
  DeclException3(ExcInvalidBoundaryFunction,
                 types::boundary_id,
                 int,
                 int,
                 << "You provided a function map that for boundary indicator "
                 << arg1 << " specifies a function with " << arg2
                 << " vector components. However, the finite "
                    "element in use has "
                 << arg3
                 << " components, and these two numbers need to match.");
  DeclException2(ExcIncompatibleNumberOfElements,
                 int,
                 int,
                 << "The number of input vectors, " << arg1
                 << " needs to be equal to the number of output vectors, "
                 << arg2
                 << ". This is not the case in your call of this function.");
  DeclExceptionMsg(ExcNoSolutions,
                   "You need to specify at least one solution vector as "
                   "input.");
};



DEAL_II_NAMESPACE_CLOSE

#endif