parsed_convergence_table.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2019 - 2021 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_base_parsed_convergence_table_h
#define dealii_base_parsed_convergence_table_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/convergence_table.h>
#include <deal.II/base/function.h>
#include <deal.II/base/function_parser.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/utilities.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/fe/mapping_q.h>
 
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
 
#include <deal.II/lac/vector.h>
 
#include <deal.II/numerics/vector_tools_integrate_difference.h>
 
DEAL_II_NAMESPACE_OPEN
 
class ParsedConvergenceTable
{
public:
  ParsedConvergenceTable(
    const std::vector<std::string> &                    component_names = {"u"},
    const std::vector<std::set<VectorTools::NormType>> &list_of_error_norms = {
      {VectorTools::H1_norm, VectorTools::L2_norm, VectorTools::Linfty_norm}});
 
  ParsedConvergenceTable(
    const std::vector<std::string> &                    component_names,
    const std::vector<std::set<VectorTools::NormType>> &list_of_error_norms,
    const double                                        exponent,
    const std::set<std::string> &                       extra_columns,
    const std::string &                                 rate_key,
    const std::string &                                 rate_mode,
    const std::string &                                 error_file_name,
    const unsigned int                                  precision,
    const bool                                          compute_error);
 
  void
  add_parameters(ParameterHandler &prm);
 
  template <int dim, int spacedim, typename VectorType>
  void
  error_from_exact(const DoFHandler<dim, spacedim> &vspace,
                   const VectorType &               solution,
                   const Function<spacedim> &       exact,
                   const Function<spacedim> *       weight = nullptr);
 
  template <int dim, int spacedim, typename VectorType>
  void
  error_from_exact(const Mapping<dim, spacedim> &   mapping,
                   const DoFHandler<dim, spacedim> &vspace,
                   const VectorType &               solution,
                   const Function<spacedim> &       exact,
                   const Function<spacedim> *       weight = nullptr);
 
  void
  add_extra_column(const std::string &            column_name,
                   const std::function<double()> &custom_function,
                   const bool                     compute_rate = true);
 
  template <int dim, int spacedim, typename VectorType>
  void
  difference(const DoFHandler<dim, spacedim> &,
             const VectorType &,
             const VectorType &,
             const Function<spacedim> *weight = nullptr);
 
  template <int dim, int spacedim, typename VectorType>
  void
  difference(const Mapping<dim, spacedim> &mapping,
             const DoFHandler<dim, spacedim> &,
             const VectorType &,
             const VectorType &,
             const Function<spacedim> *weight = nullptr);
 
  void
  output_table(std::ostream &out);
 
  void
  output_table();
 
private:
  void
  prepare_table_for_output();
 
  const std::vector<std::string> component_names;
 
  const std::vector<std::string> unique_component_names;
 
  const std::vector<ComponentMask> unique_component_masks;
 
  std::map<std::string, std::pair<std::function<double()>, bool>>
    extra_column_functions;
 
  std::vector<std::set<VectorTools::NormType>> norms_per_unique_component;
 
  double exponent;
 
  ConvergenceTable table;
 
  std::set<std::string> extra_columns;
 
  std::string rate_key;
 
  std::string rate_mode;
 
  unsigned int precision;
 
  std::string error_file_name;
 
  bool compute_error;
};
 
 
 
#ifndef DOXYGEN
// ============================================================
// Template functions
// ============================================================
template <int dim, int spacedim, typename VectorType>
void
ParsedConvergenceTable::difference(const DoFHandler<dim, spacedim> &dh,
                                   const VectorType &               solution1,
                                   const VectorType &               solution2,
                                   const Function<spacedim> *       weight)
{
  AssertThrow(solution1.size() == solution2.size(),
              ExcDimensionMismatch(solution1.size(), solution2.size()));
  VectorType solution(solution1);
  solution -= solution2;
  error_from_exact(
    dh,
    solution,
    Functions::ConstantFunction<spacedim>(0, component_names.size()),
    weight);
}
 
 
 
template <int dim, int spacedim, typename VectorType>
void
ParsedConvergenceTable::difference(const Mapping<dim, spacedim> &   mapping,
                                   const DoFHandler<dim, spacedim> &dh,
                                   const VectorType &               solution1,
                                   const VectorType &               solution2,
                                   const Function<spacedim> *       weight)
{
  AssertThrow(solution1.size() == solution2.size(),
              ExcDimensionMismatch(solution1.size(), solution2.size()));
  VectorType solution(solution1);
  solution -= solution2;
  error_from_exact(
    mapping,
    dh,
    solution,
    Functions::ConstantFunction<spacedim>(0, component_names.size()),
    weight);
}
 
 
 
template <int dim, int spacedim, typename VectorType>
void
ParsedConvergenceTable::error_from_exact(const DoFHandler<dim, spacedim> &dh,
                                         const VectorType &        solution,
                                         const Function<spacedim> &exact,
                                         const Function<spacedim> *weight)
{
  error_from_exact(get_default_linear_mapping(dh.get_triangulation()),
                   dh,
                   solution,
                   exact,
                   weight);
}
 
 
 
template <int dim, int spacedim, typename VectorType>
void
ParsedConvergenceTable::error_from_exact(const Mapping<dim, spacedim> &mapping,
                                         const DoFHandler<dim, spacedim> &dh,
                                         const VectorType &        solution,
                                         const Function<spacedim> &exact,
                                         const Function<spacedim> *weight)
{
  const auto n_components = component_names.size();
 
  if (compute_error)
    {
      AssertDimension(exact.n_components, n_components);
      AssertDimension(dh.get_fe().n_components(), n_components);
 
      const types::global_cell_index n_active_cells =
        dh.get_triangulation().n_global_active_cells();
      const unsigned int n_dofs = dh.n_dofs();
 
      for (const auto &col : extra_columns)
        if (col == "cells")
          {
            table.add_value("cells", n_active_cells);
            table.set_tex_caption("cells", "\\# cells");
            table.set_tex_format("cells", "r");
          }
        else if (col == "dofs")
          {
            table.add_value("dofs", n_dofs);
            table.set_tex_caption("dofs", "\\# dofs");
            table.set_tex_format("dofs", "r");
          }
 
      // A vector of zero std::functions with n_components components
      const std::vector<std::function<double(const Point<spacedim> &)>>
        zero_components(n_components,
                        [](const Point<spacedim> &) { return 0.0; });
 
      // The default weight function, with n_components components
      std::vector<std::function<double(const Point<spacedim> &)>>
        weight_components(n_components,
                          [](const Point<spacedim> &) { return 1.0; });
 
      if (weight != nullptr)
        {
          if (weight->n_components == 1)
            {
              for (auto &f : weight_components)
                f = [&](const Point<spacedim> &p) { return weight->value(p); };
            }
          else
            {
              AssertDimension(weight->n_components, n_components);
              for (unsigned int i = 0; i < n_components; ++i)
                weight_components[i] = [&](const Point<spacedim> &p) {
                  return weight->value(p, i);
                };
            }
        }
 
      for (unsigned int i = 0; i < norms_per_unique_component.size(); ++i)
        {
          std::map<VectorTools::NormType, double> errors;
 
          const auto &norms = norms_per_unique_component[i];
          const auto &mask  = unique_component_masks[i];
 
          // Simple case first
          if (norms.empty())
            continue;
 
          auto components_expr = zero_components;
          for (unsigned int i = 0; i < n_components; ++i)
            if (mask[i] == true)
              components_expr[i] = weight_components[i];
 
          FunctionFromFunctionObjects<spacedim> select_component(
            components_expr);
 
          Vector<float> difference_per_cell(
            dh.get_triangulation().n_global_active_cells());
 
          QGauss<dim> q_gauss((dh.get_fe().degree + 1) * 2);
 
          for (const auto &norm : norms)
            {
              difference_per_cell = 0;
              VectorTools::integrate_difference(mapping,
                                                dh,
                                                solution,
                                                exact,
                                                difference_per_cell,
                                                q_gauss,
                                                norm,
                                                &select_component,
                                                exponent);
 
              errors[norm] = VectorTools::compute_global_error(
                dh.get_triangulation(), difference_per_cell, norm, exponent);
 
              std::string name = unique_component_names[i] + "_" +
                                 Patterns::Tools::to_string(norm);
              std::string latex_name = "@f$\\| " + unique_component_names[i] +
                                       " - " + unique_component_names[i] +
                                       "_h \\|_{" +
                                       Patterns::Tools::to_string(norm) + "}@f$";
 
              table.add_value(name, errors[norm]);
              table.set_precision(name, precision);
              table.set_scientific(name, true);
              table.set_tex_caption(name, latex_name);
            }
        }
 
      for (const auto &extra_col : extra_column_functions)
        {
          const double custom_error = extra_col.second.first();
 
          std::string name = extra_col.first;
          table.add_value(name, custom_error);
          table.set_precision(name, precision);
          table.set_scientific(name, true);
        }
    }
}
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif