error_estimator.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__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/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 typename FunctionMap<spacedim>::type &neumann_bc,
   const InputVector                          &solution,
   Vector<float>                              &error,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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 typename FunctionMap<spacedim>::type &neumann_bc,
   const InputVector                          &solution,
   Vector<float>                              &error,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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 typename FunctionMap<spacedim>::type &neumann_bc,
   const std::vector<const InputVector *>     &solutions,
   std::vector<Vector<float>*>                &errors,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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 typename FunctionMap<spacedim>::type &neumann_bc,
   const std::vector<const InputVector *>     &solutions,
   std::vector<Vector<float>*>                &errors,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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 typename FunctionMap<spacedim>::type &neumann_bc,
   const InputVector                          &solution,
   Vector<float>                              &error,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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 typename FunctionMap<spacedim>::type &neumann_bc,
   const InputVector                          &solution,
   Vector<float>                              &error,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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 typename FunctionMap<spacedim>::type &neumann_bc,
   const std::vector<const InputVector *>     &solutions,
   std::vector<Vector<float>*>                &errors,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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 typename FunctionMap<spacedim>::type &neumann_bc,
   const std::vector<const InputVector *>     &solutions,
   std::vector<Vector<float>*>                &errors,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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:
  template <typename InputVector, typename DoFHandlerType>
  static void estimate
  (const Mapping<1,spacedim>                  &mapping,
   const DoFHandlerType                       &dof,
   const Quadrature<0>                        &quadrature,
   const typename FunctionMap<spacedim>::type &neumann_bc,
   const InputVector                          &solution,
   Vector<float>                              &error,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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);

  template <typename InputVector, typename DoFHandlerType>
  static void estimate
  (const DoFHandlerType                       &dof,
   const Quadrature<0>                        &quadrature,
   const typename FunctionMap<spacedim>::type &neumann_bc,
   const InputVector                          &solution,
   Vector<float>                              &error,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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);

  template <typename InputVector, typename DoFHandlerType>
  static void estimate
  (const Mapping<1,spacedim>                  &mapping,
   const DoFHandlerType                       &dof,
   const Quadrature<0>                        &quadrature,
   const typename FunctionMap<spacedim>::type &neumann_bc,
   const std::vector<const InputVector *>     &solutions,
   std::vector<Vector<float>*>                &errors,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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);

  template <typename InputVector, typename DoFHandlerType>
  static void estimate
  (const DoFHandlerType                       &dof,
   const Quadrature<0>                        &quadrature,
   const typename FunctionMap<spacedim>::type &neumann_bc,
   const std::vector<const InputVector *>     &solutions,
   std::vector<Vector<float>*>                &errors,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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);


  template <typename InputVector, typename DoFHandlerType>
  static void estimate
  (const Mapping<1,spacedim>                  &mapping,
   const DoFHandlerType                       &dof,
   const hp::QCollection<0>                   &quadrature,
   const typename FunctionMap<spacedim>::type &neumann_bc,
   const InputVector                          &solution,
   Vector<float>                              &error,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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);


  template <typename InputVector, typename DoFHandlerType>
  static void estimate
  (const DoFHandlerType                       &dof,
   const hp::QCollection<0>                   &quadrature,
   const typename FunctionMap<spacedim>::type &neumann_bc,
   const InputVector                          &solution,
   Vector<float>                              &error,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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);


  template <typename InputVector, typename DoFHandlerType>
  static void estimate
  (const Mapping<1,spacedim>                  &mapping,
   const DoFHandlerType                       &dof,
   const hp::QCollection<0>                   &quadrature,
   const typename FunctionMap<spacedim>::type &neumann_bc,
   const std::vector<const InputVector *>     &solutions,
   std::vector<Vector<float>*>                &errors,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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);


  template <typename InputVector, typename DoFHandlerType>
  static void estimate
  (const DoFHandlerType                       &dof,
   const hp::QCollection<0>                   &quadrature,
   const typename FunctionMap<spacedim>::type &neumann_bc,
   const std::vector<const InputVector *>     &solutions,
   std::vector<Vector<float>*>                &errors,
   const ComponentMask                        &component_mask = ComponentMask(),
   const Function<spacedim>                   *coefficients   = 0,
   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);

  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