]> https://gitweb.dealii.org/ - dealii.git/commitdiff
ParsedConvergenceTable
authorLuca Heltai <luca.heltai@sissa.it>
Tue, 9 Apr 2019 22:58:34 +0000 (00:58 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 30 Apr 2019 15:45:34 +0000 (17:45 +0200)
include/deal.II/base/parsed_convergence_table.h [new file with mode: 0644]
include/deal.II/base/parsed_convergence_table_flags.h [new file with mode: 0644]
source/base/CMakeLists.txt
source/base/parsed_convergence_table.cc [new file with mode: 0644]

diff --git a/include/deal.II/base/parsed_convergence_table.h b/include/deal.II/base/parsed_convergence_table.h
new file mode 100644 (file)
index 0000000..669fd6e
--- /dev/null
@@ -0,0 +1,766 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2019 by the deal.II authors
+//
+//    This file is part of the deal.II library.
+//
+//    The deal2lkit 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 deal2lkit distribution.
+//
+//-----------------------------------------------------------
+
+#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/parsed_convergence_table_flags.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.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * @brief The ParsedConvergenceTable class
+ *
+ * This class simplifies the construction of convergence tables, reading the
+ * options for the generation of the table from a parameter file. It provides a
+ * series of methods that can be used to compute the error given a reference
+ * exact solution, or the difference between two numerical solutions, or any
+ * other custom computation of the error, given via std::function objects.
+ *
+ * An example usage of this class is given by
+ * @code
+ * ParsedConvergenceTable table;
+ *
+ * ParameterHandler prm;
+ * table.add_parameters(prm);
+ *
+ * for(unsigned int i=0; i<n_cycles; ++i) {
+ *   ... // do some computations
+ *   table.error_from_exact(dof_handler, solution, exact_solution);
+ * }
+ * table.output_table(std::cout);
+ * @endcode
+ *
+ * The above code constructs a ParsedConvergenceTable that works for
+ * scalar problems, and will produce an error table with `H1`, `L2`, and
+ * `Linfty` norms of the error.
+ *
+ * Whenever a call to the methods error_from_exact() or difference() is made,
+ * the instance of this class inspects its parameters, and computes all norms
+ * specified by the parameter given at construction time, possibly modified
+ * via a ParameterFile.
+ *
+ * With a small modification, the same code can be used to estimate the errors
+ * of mixed or multi-physics problems, e.g.:
+ * @code
+ * using namespace ParsedConvergenceTableFlags;
+ * ParsedConvergenceTable table({"u,u,p"},{{H1, L2}, {L2}});
+ *
+ * ParameterHandler prm;
+ * table.add_parameters(prm);
+ *
+ * for(unsigned int i=0; i<n_cycles; ++i) {
+ *   ... // do some computations
+ *   table.error_from_exact(dof_handler, solution, exact_solution);
+ * }
+ * table.output_table(std::cout);
+ * @endcode
+ *
+ * The above code assumes that you are solving a Stokes problem with three
+ * components. Two components for the vector velocity field `u`, and one
+ * component for the pressure field `p`, and will produce an error
+ * table with `H1` and `L2` norm of the error in the velocity field (first two,
+ * components) and `L2` error in the pressure field.
+ *
+ * By calling the method add_parameters() passing a ParameterHandler object,
+ * the following options will be defined in the given ParameterHandler object,
+ * and can be modified at run time through a parameter file:
+ * @code
+ * set Enable computation of the errors = true
+ * set Enable output to streams         = true
+ * set Error file name                  =
+ * set Error precision                  = 3
+ * set Extra columns                    = dofs, cells
+ * set List of error norms to compute   = Linfty, L2, H1
+ * set Rate key                         = dofs
+ * @endcode
+ *
+ * When using this class, please cite
+ * @code{.bib}
+ * @article{SartoriGiulianiBardelloni-2018-a,
+ *     Author = {Sartori, Alberto and Giuliani, Nicola and
+ *            Bardelloni, Mauro and Heltai, Luca},
+ *     Journal = {SoftwareX},
+ *     Pages = {318--327},
+ *     Title = {{deal2lkit: A toolkit library for high performance
+ *            programming in deal.II}},
+ *  Doi = {10.1016/j.softx.2018.09.004},
+ *     Volume = {7},
+ *     Year = {2018}}
+ * @endcode
+ *
+ * @author Luca Heltai, 2019.
+ */
+class ParsedConvergenceTable
+{
+public:
+  /**
+   * Minimal constructor for ParsedConvergenceTable objects.
+   *
+   * The number of components must match the number of components of the
+   * finite element space that is used to compute the errors. If
+   * a component name is repeated, than it is interpreted as a vector field,
+   * and the errors of the repeated components are grouped together.
+   *
+   * The size of the vector @p list_of_error_norms must match the number of
+   * unique component names, and may contain zero or more comma separated
+   * strings that identify the norm types to compute for each component.
+   *
+   * For example, the following constructor
+   * @code
+   * using namespace ParsedConvergenceTableFlags;
+   * ParsedConvergenceTable table({"u", "v", "v"}, {{Linfty}, {L2, H1}});
+   * @endcode
+   *
+   * Will produce (if the parameter file is left untouched) the following table:
+   * @code
+   * cells dofs    u_Linfty         v_L2           v_H1
+   * 4     9    1.183e-01 -    5.156e-02 -    2.615e-01 -
+   * 16    25   3.291e-02 2.50 1.333e-02 2.65 1.272e-01 1.41
+   * 64    81   8.449e-03 2.31 3.360e-03 2.34 6.313e-02 1.19
+   * 256   289  2.126e-03 2.17 8.418e-04 2.18 3.150e-02 1.09
+   * 1024  1089 5.325e-04 2.09 2.106e-04 2.09 1.574e-02 1.05
+   * @endcode
+   *
+   * Se the documentation of NormType for a list of
+   * available norms you can compute using this class, and see the other
+   * constructor for a documentation of all the parameters you can change.
+   *
+   * @param component_names Specify the names of the components
+   * @param list_of_error_norms Specify what error norms to compute for each
+   * unique component name
+   */
+  ParsedConvergenceTable(
+    const std::vector<std::string> &component_names = {"u"},
+    const std::vector<std::set<ParsedConvergenceTableFlags::NormType>>
+      &list_of_error_norms = {{ParsedConvergenceTableFlags::H1,
+                               ParsedConvergenceTableFlags::L2,
+                               ParsedConvergenceTableFlags::Linfty}});
+
+  /**
+   * Full constructor for ParsedConvergenceTable.
+   *
+   * @param component_names Names of the components. Repeated consecutive
+   * names are interpreted as components of a vector valued field;
+   * @param list_of_error_norms Specify what error norms to compute for each
+   * unique component name;
+   * @param extra_columns Extra columns to add;
+   * @param rate_key Specify the extra column by which we will compute the
+   * error rates;
+   * @param error_file_name Name of error output file (with extension txt,
+   * gpl, tex, or org). If different from the empty string, than
+   * output_table() also writes in this file in the format deduced from its
+   * extension;
+   * @param precision How many digits to use when writing the error;
+   * @param compute_error If set to `false`, then no error computation is
+   * performed, and no output is produced;
+   * @param output_error If set to `false`, then the call to output_table()
+   * will produce only an output file (if @p error_file_name is not empty),
+   * but no output will be written on the given stream.
+   *
+   * The parameters you specify with this constructor can be written to a
+   * ParameterHandler object by calling the add_parameters() method. Once you
+   * call the add_parameters() method, the following options will be defined in
+   * the given ParameterHandler object, and the parameters of the instance of
+   * this class will follow the modification you make to the ParameterHandler
+   * object at run time:
+   * @code
+   * # Listing of Parameters
+   * # ---------------------
+   * # When set to false, no computations are performed.
+   * set Enable computation of the errors = true
+   *
+   * # When set to false, printing of the convergence table to the stream
+   * # specified as input to output_table() is disabled.
+   * set Enable output to streams         = true
+   *
+   * # Set this to a filename with extension .txt, .gpl, .org, or .tex to enable
+   * # writing the convergence table to a file.
+   * set Error file name                  =
+   *
+   * # Number of digits to use when printing the error.
+   * set Error precision                  = 3
+   *
+   * # Extra columns to add to the table. Availabl options are dofs, cells, and
+   * # dt.
+   * set Extra columns                    = dofs, cells
+   *
+   * # Each component is separated by a semicolon, and each norm by a comma.
+   * # Implemented norms are Linfty, L2, W1infty, H1, and custom. If you want to
+   * # skip a component, use none.
+   * set List of error norms to compute   = Linfty, L2, H1
+   *
+   * # Key to use when computing convergence rates. If this is set to a
+   * # column that is not present, or to none, then no error rates are computed.
+   * set Rate key                         = dofs
+   * @endcode
+   */
+  ParsedConvergenceTable(
+    const std::vector<std::string> &component_names,
+    const std::vector<std::set<ParsedConvergenceTableFlags::NormType>>
+      &list_of_error_norms,
+    const std::set<ParsedConvergenceTableFlags::ExtraColumns> &extra_columns,
+    const ParsedConvergenceTableFlags::ExtraColumns &          rate_key,
+    const std::string &                                        error_file_name,
+    const unsigned int &                                       precision,
+    const bool &                                               compute_error,
+    const bool &                                               output_error);
+
+  /**
+   * Attach all the parameters in this class to entries of the parameter
+   * handler @p prm. Whenever the content of @p prm changes, the parameters
+   * of this class will be updated.
+   */
+  void
+  add_parameters(ParameterHandler &prm);
+
+  /**
+   * Add a row to the error table, containing the error between @p solution and
+   * the @p exact function, in the norm(s) specified in the parameter file.
+   *
+   * If you specify a @p weight expression during this call, then this is used
+   * to compute weighted errors. If the weight is different from "1.0", then
+   * you have to make sure that muparser is installed and enabled, otherwise an
+   * exception will be thrown.
+   */
+  template <typename DoFHandlerType, typename VectorType>
+  void
+  error_from_exact(const DoFHandlerType &                           vspace,
+                   const VectorType &                               solution,
+                   const Function<DoFHandlerType::space_dimension> &exact,
+                   double                                           dt = 0.,
+                   const std::string &weight                           = "1.0");
+
+  /**
+   * Same as above, with a different mapping.
+   */
+  template <typename DoFHandlerType, typename VectorType>
+  void
+  error_from_exact(const Mapping<DoFHandlerType::dimension,
+                                 DoFHandlerType::space_dimension> & mapping,
+                   const DoFHandlerType &                           vspace,
+                   const VectorType &                               solution,
+                   const Function<DoFHandlerType::space_dimension> &exact,
+                   double                                           dt = 0.,
+                   const std::string &weight                           = "1.0");
+
+
+  /**
+   * Call the given custom function to compute a custom error for any component
+   * for which you specify the ConvergenceTableFlags::custom norm. This function
+   * may be called several times to add different columns with different names,
+   * provided you use an appropriate function every time, and
+   * that you update @p column_name everytime.
+   *
+   * Should you wish to do so, then make sure the parameter
+   * @p add_table_extras is set to false for each call, except one, so
+   * that you only add extra columns informations once.
+   *
+   * An example usage of this class with a custom error function is the
+   * following:
+   * @code
+   * using namespace ParsedConvergenceTableFlags;
+   * ParsedConvergenceTable table({"u"}, {{"custom"}});
+   *
+   * for(unsigned int i=0; i<n_cycles; ++i) {
+   *   // ... refine and distribute dofs
+   *   auto cycle_function = [&](const unsigned int component) {
+   *      return i;
+   *   };
+   *   table.custom_error(cycle_function, dof_handler, "cycle", true);
+   * }
+   * table.output_table(std::cout);
+   * @endcode
+   *
+   * Will produce (if the parameter file is left untouched) the following table:
+   * @code
+   * cells dofs  cycle
+   * 4     9     0
+   * 16    25    1
+   * 64    81    2
+   * 256   289   3
+   * 1024  1089  4
+   * @endcode
+   *
+   * Should you wish to add another column, then the above code should be
+   * replaced by:
+   * @code
+   * using namespace ParsedConvergenceTableFlags;
+   * ParsedConvergenceTable table({"u"}, {"custom"});
+   *
+   * for(unsigned int i=0; i<n_cycles; ++i) {
+   *   // ... refine and distribute dofs
+   *   auto cycle_function = [&](const unsigned int) {
+   *      return i;
+   *   };
+   *   auto cycle_square_function = [&](const unsigned int) {
+   *      return i*i;
+   *   };
+   *   table.custom_error(cycle_function, dof_handler, "cycle", false);
+   *   table.custom_error(cycle_square_function, dof_handler,
+   *                      "cycle_square", true);
+   * }
+   * table.output_table(std::cout);
+   * @endcode
+   *
+   * This will produce the table:
+   * @code
+   * cells dofs  cycle cycle_square
+   * 4     9     0     0
+   * 16    25    1     1
+   * 64    81    2     4
+   * 256   289   3     9
+   * 1024  1089  4     16
+   * @endcode
+   */
+  template <typename DoFHandlerType>
+  void
+  custom_error(const std::function<double(const unsigned int component)>
+                 &                   custom_error_function,
+               const DoFHandlerType &dh,
+               const std::string &   column_name      = "custom",
+               const bool            add_table_extras = false,
+               const double          dt               = 0.);
+
+  /**
+   * Difference between two solutions in the same vector space.
+   */
+  template <typename DoFHandlerType, typename VectorType>
+  void
+  difference(const DoFHandlerType &,
+             const VectorType &,
+             const VectorType &,
+             double dt = 0.);
+
+  /**
+   * Same as above, with a non default mapping.
+   */
+  template <typename DoFHandlerType, typename VectorType>
+  void
+  difference(const Mapping<DoFHandlerType::dimension,
+                           DoFHandlerType::space_dimension> &mapping,
+             const DoFHandlerType &,
+             const VectorType &,
+             const VectorType &,
+             double dt = 0.);
+
+  /**
+   * Write the error table to the @p out stream, and (possibly) to the file
+   * stream specified in the parameters.
+   */
+  void
+  output_table(std::ostream &out = std::cout);
+
+private:
+  /**
+   * Names of the solution components.
+   */
+  const std::vector<std::string> component_names;
+
+  /**
+   * Same as above, but containing repeated component names only once.
+   */
+  const std::vector<std::string> unique_component_names;
+
+  /**
+   * Masks for each unique component name.
+   */
+  const std::vector<ComponentMask> unique_component_masks;
+
+  /**
+   * Type of error to compute per components.
+   */
+  std::vector<std::set<ParsedConvergenceTableFlags::NormType>>
+    norms_per_unique_component;
+
+  /**
+   * The actual table
+   */
+  ConvergenceTable table;
+  /**
+   * Extra columns to add to the table.
+   */
+  std::set<ParsedConvergenceTableFlags::ExtraColumns> extra_columns;
+
+  /**
+   * Wether or not to calculate the rates according to the given keys.
+   */
+  ParsedConvergenceTableFlags::ExtraColumns rate_key;
+
+  /**
+   * The precision used to output the table.
+   */
+  unsigned int precision;
+
+  /**
+   * Filename to use when writing to file.
+   */
+  std::string error_file_name;
+
+  /**
+   * Compute the error. If this is false, all functions regarding
+   * errors are disabled and don't do anything.
+   */
+  bool compute_error;
+
+  /**
+   * Output the error file also on screen.
+   */
+  bool output_error;
+};
+
+
+
+#ifndef DOXYGEN
+// ============================================================
+// Template instantiations
+// ============================================================
+template <typename DoFHandlerType, typename VectorType>
+void
+ParsedConvergenceTable::difference(const DoFHandlerType &dh,
+                                   const VectorType &    solution1,
+                                   const VectorType &    solution2,
+                                   double                dt)
+{
+  AssertThrow(solution1.size() == solution2.size(),
+              ExcDimensionMismatch(solution1.size(), solution2.size()));
+  VectorType solution(solution1);
+  solution -= solution2;
+  error_from_exact(dh,
+                   solution,
+                   ConstantFunction<DoFHandlerType::space_dimension>(
+                     0, component_names.size()),
+                   dt);
+}
+
+
+
+template <typename DoFHandlerType, typename VectorType>
+void
+ParsedConvergenceTable::difference(
+  const Mapping<DoFHandlerType::dimension, DoFHandlerType::space_dimension>
+    &                   mapping,
+  const DoFHandlerType &dh,
+  const VectorType &    solution1,
+  const VectorType &    solution2,
+  double                dt)
+{
+  AssertThrow(solution1.size() == solution2.size(),
+              ExcDimensionMismatch(solution1.size(), solution2.size()));
+  VectorType solution(solution1);
+  solution -= solution2;
+  error_from_exact(mapping,
+                   dh,
+                   solution,
+                   ConstantFunction<DoFHandlerType::space_dimension>(
+                     0, component_names.size()),
+                   dt);
+}
+
+
+
+template <typename DoFHandlerType, typename VectorType>
+void
+ParsedConvergenceTable::error_from_exact(
+  const DoFHandlerType &                           dh,
+  const VectorType &                               solution,
+  const Function<DoFHandlerType::space_dimension> &exact,
+  double                                           dt,
+  const std::string &                              weight)
+{
+  error_from_exact(StaticMappingQ1<DoFHandlerType::dimension,
+                                   DoFHandlerType::space_dimension>::mapping,
+                   dh,
+                   solution,
+                   exact,
+                   dt,
+                   weight);
+}
+
+
+
+template <typename DoFHandlerType, typename VectorType>
+void
+ParsedConvergenceTable::error_from_exact(
+  const Mapping<DoFHandlerType::dimension, DoFHandlerType::space_dimension>
+    &                                              mapping,
+  const DoFHandlerType &                           dh,
+  const VectorType &                               solution,
+  const Function<DoFHandlerType::space_dimension> &exact,
+  double                                           dt,
+  const std::string &                              weight)
+{
+  const int dim      = DoFHandlerType::dimension;
+  const int spacedim = DoFHandlerType::space_dimension;
+  if (compute_error)
+    {
+      AssertDimension(exact.n_components, component_names.size());
+
+      const unsigned int 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 == ParsedConvergenceTableFlags::cells)
+          {
+            table.add_value("cells", n_active_cells);
+            table.set_tex_caption("cells", "\\# cells");
+            table.set_tex_format("cells", "r");
+          }
+        else if (col == ParsedConvergenceTableFlags::dofs)
+          {
+            table.add_value("dofs", n_dofs);
+            table.set_tex_caption("dofs", "\\# dofs");
+            table.set_tex_format("dofs", "r");
+          }
+        else if (col == ParsedConvergenceTableFlags::dt)
+          {
+            table.add_value("dt", dt);
+            table.set_tex_caption("dt", "\\Delta t");
+            table.set_tex_format("dt", "r");
+          }
+
+
+      std::map<ParsedConvergenceTableFlags::NormType, unsigned int>
+        norm_to_index;
+      norm_to_index[ParsedConvergenceTableFlags::L2]      = 0;
+      norm_to_index[ParsedConvergenceTableFlags::H1]      = 1;
+      norm_to_index[ParsedConvergenceTableFlags::Linfty]  = 2;
+      norm_to_index[ParsedConvergenceTableFlags::W1infty] = 3;
+      norm_to_index[ParsedConvergenceTableFlags::custom]  = 4;
+
+      const auto n_components = component_names.size();
+      for (unsigned int i = 0; i < unique_component_names.size(); ++i)
+        {
+          std::vector<double> error(5, 0.0);
+          const auto &        norms = norms_per_unique_component[i];
+          const auto &        mask  = unique_component_masks[i];
+
+#  ifndef DEAL_II_WITH_MUPARSER
+          Assert(weight == "1.0",
+                 ExcMessage("If you want to use a weight which is "
+                            "different from 1.0, you have to "
+                            "install deal.II with muparser enabled."));
+          ComponentSelectFunction<spacedim> select_component(
+            {mask.first_selected_component(),
+             mask.first_selected_component() + mask.n_selected_components()},
+            n_components);
+#  else
+          std::string comp_sel;
+          {
+            // Select all required components
+            std::vector<std::string> expr(n_components, "0");
+
+            for (unsigned int i = 0; i < n_components; ++i)
+              if (mask[i] == true)
+                expr[i] = weight;
+
+            comp_sel = expr[0];
+            for (unsigned int i = 1; i < expr.size(); ++i)
+              comp_sel += "; " + expr[i];
+          }
+          FunctionParser<spacedim> select_component(comp_sel);
+#  endif
+
+          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)
+            {
+              if (norm == ParsedConvergenceTableFlags::L2)
+                {
+                  VectorTools::integrate_difference(mapping,
+                                                    dh,
+                                                    solution,
+                                                    exact,
+                                                    difference_per_cell,
+                                                    q_gauss,
+                                                    VectorTools::L2_norm,
+                                                    &select_component);
+                }
+
+              error[norm_to_index[ParsedConvergenceTableFlags::L2]] =
+                VectorTools::compute_global_error(dh.get_triangulation(),
+                                                  difference_per_cell,
+                                                  VectorTools::L2_norm);
+              difference_per_cell = 0;
+
+              if (norm == ParsedConvergenceTableFlags::H1)
+                {
+                  VectorTools::integrate_difference(mapping,
+                                                    dh,
+                                                    solution,
+                                                    exact,
+                                                    difference_per_cell,
+                                                    q_gauss,
+                                                    VectorTools::H1_norm,
+                                                    &select_component);
+                }
+              error[norm_to_index[ParsedConvergenceTableFlags::H1]] =
+                VectorTools::compute_global_error(dh.get_triangulation(),
+                                                  difference_per_cell,
+                                                  VectorTools::H1_norm);
+              difference_per_cell = 0;
+
+              if (norm == ParsedConvergenceTableFlags::W1infty)
+                {
+                  VectorTools::integrate_difference(mapping,
+                                                    dh,
+                                                    solution,
+                                                    exact,
+                                                    difference_per_cell,
+                                                    q_gauss,
+                                                    VectorTools::W1infty_norm,
+                                                    &select_component);
+                }
+
+              error[norm_to_index[ParsedConvergenceTableFlags::W1infty]] =
+                VectorTools::compute_global_error(
+                  dh.get_triangulation(),
+                  difference_per_cell,
+                  VectorTools::W1infty_seminorm);
+              difference_per_cell = 0;
+
+              if (norm == ParsedConvergenceTableFlags::Linfty)
+                {
+                  VectorTools::integrate_difference(mapping,
+                                                    dh,
+                                                    solution,
+                                                    exact,
+                                                    difference_per_cell,
+                                                    q_gauss,
+                                                    VectorTools::Linfty_norm,
+                                                    &select_component);
+                }
+
+              error[norm_to_index[ParsedConvergenceTableFlags::Linfty]] =
+                VectorTools::compute_global_error(dh.get_triangulation(),
+                                                  difference_per_cell,
+                                                  VectorTools::Linfty_norm);
+              difference_per_cell = 0;
+
+              // Now add what we just computed.
+              if (norm != ParsedConvergenceTableFlags::none)
+                {
+                  std::string name = unique_component_names[i] + "_" +
+                                     Patterns::Tools::to_string(norm);
+                  std::string latex_name =
+                    "$\\| " + unique_component_names[i] + " - " +
+                    unique_component_names[i] + "_h \\|_{" +
+                    Patterns::Tools::to_string(norm) + "}$";
+
+                  table.add_value(name, error[norm_to_index[norm]]);
+                  table.set_precision(name, precision);
+                  table.set_scientific(name, true);
+                  table.set_tex_caption(name, latex_name);
+                }
+            }
+        }
+    }
+}
+
+
+template <typename DoFHandlerType>
+void
+ParsedConvergenceTable::custom_error(
+  const std::function<double(const unsigned int component)>
+    &                   custom_error_function,
+  const DoFHandlerType &dh,
+  const std::string &   column_name,
+  const bool            add_table_extras,
+  const double          dt)
+{
+  if (compute_error)
+    {
+      const unsigned int  n_components = norms_per_unique_component.size();
+      std::vector<double> c_error(norms_per_unique_component.size());
+      const unsigned int  n_active_cells =
+        dh.get_triangulation().n_global_active_cells();
+      const unsigned int n_dofs = dh.n_dofs();
+      if (add_table_extras)
+        {
+          for (const auto &col : extra_columns)
+            {
+              if (col == ParsedConvergenceTableFlags::cells)
+                {
+                  table.add_value("cells", n_active_cells);
+                  table.set_tex_caption("cells", "\\# cells");
+                  table.set_tex_format("cells", "r");
+                }
+              if (col == ParsedConvergenceTableFlags::dofs)
+                {
+                  table.add_value("dofs", n_dofs);
+                  table.set_tex_caption("dofs", "\\# dofs");
+                  table.set_tex_format("dofs", "r");
+                }
+              if (col == ParsedConvergenceTableFlags::dt)
+                {
+                  table.add_value("dt", dt);
+                  table.set_tex_caption("dt", "\\Delta t");
+                  table.set_tex_format("dt", "r");
+                }
+            }
+        }
+
+      for (unsigned int j = 0; j < n_components; ++j)
+        {
+          const auto &norms = norms_per_unique_component[j];
+
+          for (const auto &norm : norms)
+            if (norm == ParsedConvergenceTableFlags::custom)
+              {
+                const double custom_error = custom_error_function(j);
+
+                std::string name =
+                  unique_component_names[j] + "_" + column_name;
+                std::string latex_name = "$\\| " + unique_component_names[j] +
+                                         " - " + unique_component_names[j] +
+                                         "_h \\|_{\\text{" + column_name +
+                                         "} $";
+
+                table.add_value(name, custom_error);
+                table.set_precision(name, precision);
+                table.set_scientific(name, true);
+                table.set_tex_caption(name, latex_name);
+              }
+        }
+    }
+}
+#endif
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/include/deal.II/base/parsed_convergence_table_flags.h b/include/deal.II/base/parsed_convergence_table_flags.h
new file mode 100644 (file)
index 0000000..7ff300d
--- /dev/null
@@ -0,0 +1,232 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2019 by the deal.II authors
+//
+//    This file is part of the deal.II library.
+//
+//    The deal2lkit 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 deal2lkit distribution.
+//
+//-----------------------------------------------------------
+
+#ifndef dealii_base_parsed_convergence_table_flags_h
+#define dealii_base_parsed_convergence_table_flags_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/patterns.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace ParsedConvergenceTableFlags
+{
+  /**
+   * Implemented Error norms.
+   */
+  enum NormType
+  {
+    none    = 0x00,
+    Linfty  = 0x01,
+    L2      = 0x02,
+    W1infty = 0x04,
+    H1      = 0x08,
+    custom  = 0x10,
+  };
+
+  /**
+   * The ExtraColumns enum.
+   */
+  enum ExtraColumns
+  {
+    no_columns = 0x0,
+    dofs       = 0x1,
+    cells      = 0x2,
+    dt         = 0x4
+  };
+} // namespace ParsedConvergenceTableFlags
+
+namespace Patterns
+{
+  namespace Tools
+  {
+    template <>
+    struct Convert<ParsedConvergenceTableFlags::NormType>
+    {
+      /**
+       * Return the Correct pattern for NormType.
+       */
+      static std::unique_ptr<Patterns::PatternBase>
+      to_pattern()
+      {
+        return std::make_unique<Patterns::Selection>(
+          "none|Linfty|L2|W1infty|H1|custom");
+      }
+
+
+
+      /**
+       * Convert a NormType to a string.
+       */
+      static std::string
+      to_string(const ParsedConvergenceTableFlags::NormType & s,
+                const std::unique_ptr<Patterns::PatternBase> &p =
+                  Convert<ParsedConvergenceTableFlags::NormType>::to_pattern())
+      {
+        (void)p;
+        if (s == ParsedConvergenceTableFlags::Linfty)
+          return "Linfty";
+        else if (s == ParsedConvergenceTableFlags::L2)
+          return "L2";
+        else if (s == ParsedConvergenceTableFlags::W1infty)
+          return "W1infty";
+        else if (s == ParsedConvergenceTableFlags::H1)
+          return "H1";
+        else if (s == ParsedConvergenceTableFlags::custom)
+          return "custom";
+        else if (s == ParsedConvergenceTableFlags::none)
+          return "none";
+        else
+          {
+            AssertThrow(false, ExcMessage("Didn't recognize a norm type."));
+            return "";
+          }
+      }
+
+
+
+      /**
+       * Convert a string to a NormType.
+       */
+      static ParsedConvergenceTableFlags::NormType
+      to_value(const std::string &                           norm_string,
+               const std::unique_ptr<Patterns::PatternBase> &p =
+                 Convert<ParsedConvergenceTableFlags::NormType>::to_pattern())
+      {
+        (void)p;
+        ParsedConvergenceTableFlags::NormType norm =
+          ParsedConvergenceTableFlags::none;
+
+        if (norm_string == "Linfty")
+          {
+            norm = ParsedConvergenceTableFlags::Linfty;
+          }
+        else if (norm_string == "L2")
+          {
+            norm = ParsedConvergenceTableFlags::L2;
+          }
+        else if (norm_string == "W1infty")
+          {
+            norm = ParsedConvergenceTableFlags::W1infty;
+          }
+        else if (norm_string == "H1")
+          {
+            norm = ParsedConvergenceTableFlags::H1;
+          }
+        else if (norm_string == "custom")
+          {
+            norm = ParsedConvergenceTableFlags::custom;
+          }
+        else if (norm_string == "none")
+          {
+            norm = ParsedConvergenceTableFlags::none;
+          }
+        else
+          {
+            AssertThrow(false, ExcMessage("Didn't recognize a norm type."));
+          }
+
+        return norm;
+      }
+    };
+
+
+
+    template <>
+    struct Convert<ParsedConvergenceTableFlags::ExtraColumns>
+    {
+      /**
+       * Return the Correct pattern for ExtraColumns.
+       */
+      static std::unique_ptr<Patterns::PatternBase>
+      to_pattern()
+      {
+        return std::make_unique<Patterns::Selection>("none|cells|dofs|dt");
+      }
+
+
+
+      /**
+       * Convert an ExtraColumn object to a string.
+       */
+      static std::string
+      to_string(
+        const ParsedConvergenceTableFlags::ExtraColumns &s,
+        const std::unique_ptr<Patterns::PatternBase> &   p =
+          Convert<ParsedConvergenceTableFlags::ExtraColumns>::to_pattern())
+      {
+        (void)p;
+        if (s == ParsedConvergenceTableFlags::cells)
+          return "cells";
+        else if (s == ParsedConvergenceTableFlags::dofs)
+          return "dofs";
+        else if (s == ParsedConvergenceTableFlags::dt)
+          return "dt";
+        else if (s == ParsedConvergenceTableFlags::no_columns)
+          return "none";
+        else
+          {
+            AssertThrow(false, ExcMessage("Didn't recognize a column type."));
+            return "";
+          }
+      }
+
+
+
+      /**
+       * Convert a string to an ExtraColumn object.
+       */
+      static ParsedConvergenceTableFlags::ExtraColumns
+      to_value(
+        const std::string &                           column_string,
+        const std::unique_ptr<Patterns::PatternBase> &p =
+          Convert<ParsedConvergenceTableFlags::ExtraColumns>::to_pattern())
+      {
+        (void)p;
+        ParsedConvergenceTableFlags::ExtraColumns column =
+          ParsedConvergenceTableFlags::no_columns;
+
+        if (column_string == "cells")
+          {
+            column = ParsedConvergenceTableFlags::cells;
+          }
+        else if (column_string == "dofs")
+          {
+            column = ParsedConvergenceTableFlags::dofs;
+          }
+        else if (column_string == "dt")
+          {
+            column = ParsedConvergenceTableFlags::dt;
+          }
+        else if (column_string == "none")
+          {
+            column = ParsedConvergenceTableFlags::no_columns;
+          }
+        else
+          {
+            AssertThrow(false,
+                        ExcMessage("Didn't recognize an extra column type."));
+          }
+
+        return column;
+      }
+    };
+  } // namespace Tools
+} // namespace Patterns
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index eacda918d507ac6e557ca12698f76839c53fa77a..e3c3f8ee0afcb872671703dc24c65c601c4d16eb 100644 (file)
@@ -51,6 +51,7 @@ SET(_unity_include_src
   parallel.cc
   parameter_handler.cc
   parameter_acceptor.cc
+  parsed_convergence_table.cc
   parsed_function.cc
   partitioner.cc
   patterns.cc
diff --git a/source/base/parsed_convergence_table.cc b/source/base/parsed_convergence_table.cc
new file mode 100644 (file)
index 0000000..0139d97
--- /dev/null
@@ -0,0 +1,201 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/parsed_convergence_table.h>
+#include <deal.II/base/patterns.h>
+#include <deal.II/base/utilities.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace
+{
+  std::vector<std::string>
+  get_unique_component_names(const std::vector<std::string> &component_names)
+  {
+    auto elements = component_names;
+    auto last     = std::unique(elements.begin(), elements.end());
+    elements.resize(last - elements.begin());
+    return elements;
+  }
+
+
+
+  std::vector<ComponentMask>
+  get_unique_component_masks(const std::vector<std::string> &component_names)
+  {
+    const auto unique_component_names =
+      get_unique_component_names(component_names);
+
+    std::vector<ComponentMask> masks;
+
+    std::vector<std::vector<bool>> bools(
+      unique_component_names.size(),
+      std::vector<bool>(component_names.size(), false));
+
+    unsigned int j = 0;
+    for (unsigned int i = 0; i < component_names.size(); ++i)
+      {
+        if (unique_component_names[j] != component_names[i])
+          masks.push_back(ComponentMask(bools[j++]));
+        bools[j][i] = true;
+      }
+    masks.push_back(ComponentMask(bools[j++]));
+    AssertDimension(j, unique_component_names.size());
+    return masks;
+  }
+} // namespace
+
+ParsedConvergenceTable::ParsedConvergenceTable(
+  const std::vector<std::string> &solution_names,
+  const std::vector<std::set<ParsedConvergenceTableFlags::NormType>>
+    &list_of_error_norms)
+  : ParsedConvergenceTable(solution_names,
+                           list_of_error_norms,
+                           {ParsedConvergenceTableFlags::cells,
+                            ParsedConvergenceTableFlags::dofs},
+                           ParsedConvergenceTableFlags::dofs,
+                           "",
+                           3,
+                           true,
+                           true)
+{}
+
+
+
+ParsedConvergenceTable::ParsedConvergenceTable(
+  const std::vector<std::string> &component_names,
+  const std::vector<std::set<ParsedConvergenceTableFlags::NormType>>
+    &list_of_error_norms,
+  const std::set<ParsedConvergenceTableFlags::ExtraColumns> &extra_columns,
+  const ParsedConvergenceTableFlags::ExtraColumns &          rate_key,
+  const std::string &                                        error_file_name,
+  const unsigned int &                                       precision,
+  const bool &                                               compute_error,
+  const bool &                                               output_error)
+  : component_names(component_names)
+  , unique_component_names(get_unique_component_names(component_names))
+  , unique_component_masks(get_unique_component_masks(component_names))
+  , norms_per_unique_component(list_of_error_norms)
+  , extra_columns(extra_columns)
+  , rate_key(rate_key)
+  , precision(precision)
+  , error_file_name(error_file_name)
+  , compute_error(compute_error)
+  , output_error(output_error)
+{
+  AssertDimension(unique_component_names.size(), list_of_error_norms.size());
+}
+
+
+
+void
+ParsedConvergenceTable::add_parameters(ParameterHandler &prm)
+{
+  prm.add_parameter("Enable output to streams",
+                    output_error,
+                    "When set to false, printing of the convergence table "
+                    "to the stream specified as input to output_table() "
+                    "is disabled.");
+
+  prm.add_parameter("Enable computation of the errors",
+                    compute_error,
+                    "When set to false, no computations are performed.");
+  prm.add_parameter("Error precision",
+                    precision,
+                    "Number of digits to use when printing the error.",
+                    Patterns::Integer(0));
+
+  prm.add_parameter("Error file name",
+                    error_file_name,
+                    "Set this to a filename with extension .txt, .gpl, .org, "
+                    "or .tex to enable writing the convergence table to a "
+                    "file.");
+
+  prm.add_parameter(
+    "List of error norms to compute",
+    norms_per_unique_component,
+    "Each component is separated by a semicolon, "
+    "and each norm by a comma. Implemented norms are Linfty, L2, "
+    "W1infty, H1, and custom. If you want to skip a component, use none.");
+
+  prm.add_parameter("Extra columns",
+                    extra_columns,
+                    "Extra columns to add to the table. Availabl options "
+                    "are dofs, cells, and dt.");
+
+  prm.add_parameter("Rate key",
+                    rate_key,
+                    "Key to use when computing convergence rates. If "
+                    "this is set to a column that is not present, or to none, "
+                    "then no error rates are computed.");
+}
+
+
+
+void
+ParsedConvergenceTable::output_table(std::ostream &out)
+{
+  if (compute_error)
+    {
+      // Add convergence rates if the rate_key is not empty
+      if (rate_key != ParsedConvergenceTableFlags::no_columns)
+        {
+          bool has_key = false;
+          for (const auto &col : extra_columns)
+            {
+              if (rate_key == col)
+                has_key = true;
+
+              if (col != ParsedConvergenceTableFlags::no_columns)
+                table.omit_column_from_convergence_rate_evaluation(
+                  Patterns::Tools::to_string(col));
+            }
+          if (has_key)
+            table.evaluate_all_convergence_rates(
+              Patterns::Tools::to_string(rate_key),
+              ConvergenceTable::reduction_rate_log2);
+        }
+
+      if (output_error)
+        table.write_text(out);
+
+      if (error_file_name != "")
+        {
+          const std::string error_file_format =
+            error_file_name.substr(error_file_name.find_last_of(".") + 1);
+
+          std::ofstream table_file(error_file_name);
+
+          if (error_file_format == "tex")
+            table.write_tex(table_file);
+          else if (error_file_format == "txt")
+            table.write_text(table_file);
+          else if (error_file_format == "gpl")
+            table.write_text(
+              table_file, TableHandler::table_with_separate_column_description);
+          else if (error_file_format == "org")
+            table.write_text(table_file, TableHandler::org_mode_table);
+          else
+            {
+              AssertThrow(false,
+                          ExcInternalError(
+                            std::string("Unrecognized file format: ") +
+                            error_file_format));
+            }
+          table_file.close();
+        }
+    }
+}
+DEAL_II_NAMESPACE_CLOSE

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.