]> https://gitweb.dealii.org/ - code-gallery.git/commitdiff
Clean a little and fix indentation 114/head
authorMarco Feder <marco.feder@sissa.it>
Wed, 20 Jul 2022 17:38:15 +0000 (19:38 +0200)
committerMarco Feder <marco.feder@sissa.it>
Wed, 20 Jul 2022 21:57:30 +0000 (23:57 +0200)
advection_reaction_estimator/README.md
advection_reaction_estimator/include/DG_advection_reaction.h
advection_reaction_estimator/main.cc
advection_reaction_estimator/source/DG_advection_reaction.cc

index b33d7940355741324c08f6e446b703e6ef42c580..f243a2ab5e8ffaed9db629a4cca58cd23c4ea80f 100644 (file)
@@ -133,15 +133,15 @@ This solution has an internal layer along the line $y=\frac{1}{2} -x$, hence we
 
 The next image is the 3D view of the numerical solution:
 
-![Screenshot](doc/images/warp_by_scalar_solution_layer.png)
+![Screenshot](./doc/images/warp_by_scalar_solution_layer.png)
 
 More interestingly, we see that the estimator has been able to capture the layer. Here a bulk-chasing criterion is used, with bottom fraction ´0.5´ and no coarsening. This mesh is obtained after 12 refinement cycles.
-![Screenshot](doc/images/refined_mesh_internal_layer.png)
+![Screenshot](./doc/images/refined_mesh_internal_layer.png)
 
 
 If we look at the decrease of the energy norm of the error in the globally refined case and in the adaptively case, with respect to the DoFs, we obtain:
 
-![Screenshot](doc/images/adaptive_vs_global_refinement.png)
+![Screenshot](./doc/images/adaptive_vs_global_refinement.png)
 
 ## References 
 * [1] Emmanuil H. Georgoulis, Edward Hall and Charalambos Makridakis (2013), Error Control for Discontinuous Galerkin Methods for First Order Hyperbolic Problems. DOI: [10.1007/978-3-319-01818-8_8
index 0146eab2f5766ef4a9440e7c5c3e15fde57495e2..24c5235de7b0afd2be3364d2bf1a9ca30f265f06 100644 (file)
 #ifndef INCLUDE_DG_UPWIND_H_
 #define INCLUDE_DG_UPWIND_H_
 
-// The first few files have already been covered in  tutorials and will
-// thus not be further commented on:
-#include <deal.II/base/quadrature_lib.h>
+
 #include <deal.II/base/function.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_out.h>
-#include <deal.II/grid/grid_refinement.h>
-#include <deal.II/fe/fe_values.h>
+#include <deal.II/base/quadrature_lib.h>
+
 #include <deal.II/dofs/dof_handler.h>
-#include <deal.II/numerics/vector_tools.h>
 #include <deal.II/dofs/dof_tools.h>
-#include <deal.II/numerics/data_out.h>
-#include <deal.II/fe/mapping_q1.h>
 
-#include <deal.II/lac/full_matrix.h>
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
 // This header is needed for FEInterfaceValues to compute integrals on
 // interfaces:
 #include <deal.II/fe/fe_interface_values.h>
-//Solver
-#include <deal.II/lac/solver_richardson.h>
+// Solver
 #include <deal.II/lac/precondition_block.h>
+#include <deal.II/lac/solver_richardson.h>
 #include <deal.II/lac/sparse_direct.h>
 // We are going to use gradients as refinement indicator.
 #include <deal.II/numerics/derivative_approximation.h>
 // Using using the mesh_loop from the MeshWorker framework
-#include <deal.II/meshworker/mesh_loop.h>
-
 #include <deal.II/base/convergence_table.h>
 
-//To enable parameter handling
+#include <deal.II/meshworker/mesh_loop.h>
+
+// To enable parameter handling
 #include <deal.II/base/function_parser.h>
 #include <deal.II/base/parameter_acceptor.h>
-#include <deal.II/base/parsed_convergence_table.h>
 #include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/parsed_convergence_table.h>
 #include <deal.II/base/symbolic_function.h>
 
 #include <deal.II/meshworker/copy_data.h>
 #include <deal.II/meshworker/mesh_loop.h>
 #include <deal.II/meshworker/scratch_data.h>
 
-#include <iostream>
 #include <fstream>
+#include <iostream>
+using namespace dealii;
 
-//This is a struct used only for throwing an exception when theta parameter is not okay.
+// This is a struct used only for throwing an exception when theta parameter is
+// not okay.
 struct theta_exc
 {
-       std::string message;
-       theta_exc(std::string &&s) : message{std::move(s)} {};
-       const char *what() const { return message.c_str(); }
+  std::string message;
+  theta_exc(std::string &&s)
+    : message{std::move(s)} {};
+  const char *
+  what() const
+  {
+    return message.c_str();
+  }
 };
 
-using namespace dealii;
-// @sect3{Class declaration}
-// In the following we have the declaration of the functions used in the program. As we want to use
-// parameter files, we need to derive our class from `ParameterAcceptor`.
+
 template <int dim>
 class AdvectionReaction : ParameterAcceptor
 {
 public:
-       AdvectionReaction();
-       void initialize_params(const std::string &filename);
-       void run();
+  AdvectionReaction();
+  void
+  initialize_params(const std::string &filename);
+  void
+  run();
 
 private:
-       using Iterator = typename DoFHandler<dim>::active_cell_iterator;
-       void parse_string(const std::string &parameters);
-       void setup_system();
-       void assemble_system();
-       void solve();
-       void refine_grid();
-       void output_results(const unsigned int cycle) const;
-       void compute_error();
-       double compute_energy_norm();
-       void compute_local_projection_and_estimate();
-
-       Triangulation<dim> triangulation;
-       const MappingQ1<dim> mapping;
-
-       // Furthermore we want to use DG elements.
-       std::unique_ptr<FE_DGQ<dim>> fe;
-       DoFHandler<dim> dof_handler;
-
-       SparsityPattern sparsity_pattern;
-       SparseMatrix<double> system_matrix;
-
-       Vector<double> solution;
-       Vector<double> right_hand_side;
-       Vector<double> energy_norm_square_per_cell;
-       Vector<double> error_indicator_per_cell;
-
-       // So far we declared the usual objects. Hereafter we declare `FunctionParser<dim>` objects
-       FunctionParser<dim> exact_solution;
-       FunctionParser<dim> boundary_conditions;
-       FunctionParser<dim> rhs;
-       FunctionParser<dim> advection_coeff;
-
-       unsigned int fe_degree = 1;
-
-       // and then we define default values that will be parsed from the following strings
-       std::string exact_solution_expression = "tanh(100*(x+y-0.5))"; //internal layer solution
-       std::string rhs_expression = "-200*tanh(100*x + 100*y - 50.0)^2 + tanh(100*x + 100*y - 50.0) + 200";
-       std::string advection_coefficient_expression = "1.0";
-       std::string boundary_conditions_expression = "tanh(100*x + 100*y - 50.0)";
-       std::string refinement = "residual";
-       std::string output_filename = "DG_estimator";
-       std::map<std::string, double> constants;
-       ParsedConvergenceTable error_table;
-
-       bool use_direct_solver = true;
-       unsigned int n_refinement_cycles = 14;
-       unsigned int n_global_refinements = 3;
-       double theta = 0.5; //default is 0.5 so that I have classical upwind flux
+  using Iterator = typename DoFHandler<dim>::active_cell_iterator;
+  void
+  parse_string(const std::string &parameters);
+  void
+  setup_system();
+  void
+  assemble_system();
+  void
+  solve();
+  void
+  refine_grid();
+  void
+  output_results(const unsigned int cycle) const;
+  void
+  compute_error();
+  double
+  compute_energy_norm();
+  void
+  compute_local_projection_and_estimate();
+
+  Triangulation<dim>   triangulation;
+  const MappingQ1<dim> mapping;
+
+  // Furthermore we want to use DG elements.
+  std::unique_ptr<FE_DGQ<dim>> fe;
+  DoFHandler<dim>              dof_handler;
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+
+  Vector<double> solution;
+  Vector<double> right_hand_side;
+  Vector<double> energy_norm_square_per_cell;
+  Vector<double> error_indicator_per_cell;
+
+  // So far we declared the usual objects. Hereafter we declare
+  // `FunctionParser<dim>` objects
+  FunctionParser<dim> exact_solution;
+  FunctionParser<dim> boundary_conditions;
+  FunctionParser<dim> rhs;
+  FunctionParser<dim> advection_coeff;
+
+  unsigned int fe_degree = 1;
+
+  // and then we define default values that will be parsed from the following
+  // strings
+  std::string exact_solution_expression =
+    "tanh(100*(x+y-0.5))"; // internal layer solution
+  std::string rhs_expression =
+    "-200*tanh(100*x + 100*y - 50.0)^2 + tanh(100*x + 100*y - 50.0) + 200";
+  std::string advection_coefficient_expression = "1.0";
+  std::string boundary_conditions_expression   = "tanh(100*x + 100*y - 50.0)";
+  std::string refinement                       = "residual";
+  std::string output_filename = "DG_advection_reaction_estimator";
+  std::map<std::string, double> constants;
+  ParsedConvergenceTable        error_table;
+
+  bool         use_direct_solver    = true;
+  unsigned int n_refinement_cycles  = 8;
+  unsigned int n_global_refinements = 3;
+  double theta = 0.5; // default is 0.5 so that I have classical upwind flux
 };
 
 #endif /* INCLUDE_DG_UPWIND_H_ */
index 91eac652ecc0f4bc6c2f8bdd793a6e36e24e4c11..97407d32649d0212913c1b2f79f82d8e33d72864 100644 (file)
@@ -1,62 +1,62 @@
 #include "include/DG_advection_reaction.h"
 
-int main(int argc, char **argv)
+int
+main(int argc, char **argv)
 {
-
   try
-  {
-    std::string par_name = "";
-    if (argc > 1)
     {
-      par_name = argv[1];
+      std::string par_name = "";
+      if (argc > 1)
+        {
+          par_name = argv[1];
+        }
+      else
+        {
+          par_name = "parameters.prm";
+        }
+      deallog.depth_console(2);
+      AdvectionReaction<2> problem;
+      problem.initialize_params(par_name);
+      problem.run();
     }
-    deallog.depth_console(2); //solver infos
-    AdvectionReaction<2> dgmethod;
-    if (par_name != "")
+  catch (std::exception &exc)
     {
-      dgmethod.initialize_params(par_name);
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
     }
-
-    dgmethod.run();
-  }
-  catch (std::exception &exc)
-  {
-    std::cerr << std::endl
-              << std::endl
-              << "----------------------------------------------------"
-              << std::endl;
-    std::cerr << "Exception on processing: " << std::endl
-              << exc.what() << std::endl
-              << "Aborting!" << std::endl
-              << "----------------------------------------------------"
-              << std::endl;
-    return 1;
-  }
   catch (const theta_exc &theta_range)
-  {
-    std::cerr << std::endl
-              << std::endl
-              << "----------------------------------------------------"
-              << std::endl;
-    std::cerr << "Exception on processing: " << std::endl
-              << theta_range.what() << std::endl
-              << "Aborting!" << std::endl
-              << "----------------------------------------------------"
-              << std::endl;
-    return 1;
-  }
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << theta_range.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
   catch (...)
-  {
-    std::cerr << std::endl
-              << std::endl
-              << "----------------------------------------------------"
-              << std::endl;
-    std::cerr << "Unknown exception!" << std::endl
-              << "Aborting!" << std::endl
-              << "----------------------------------------------------"
-              << std::endl;
-    return 1;
-  }
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
 
   return 0;
 }
index 7a69bd73a6ec585ef41e376729db6fa81c3765a9..33cc9390576def74b00f6f6060b01cea24278d2e 100644 (file)
  * the top level directory of deal.II.
  *
  * ---------------------------------------------------------------------
-
-
-
  *
  * Author: Marco Feder, SISSA, 2021
- *         
+ *
  */
 
 #include "../include/DG_advection_reaction.h"
 
-//Compute and returns the wind field b
+// Compute and returns the wind field b
 template <int dim>
-Tensor<1, dim> beta(const Point<dim> &p)
+Tensor<1, dim>
+beta(const Point<dim> &p)
 {
-       Assert(dim >= 2, ExcNotImplemented());
-       (void)p; //suppress warnings from p
-       Tensor<1, dim> wind_field;
-       wind_field[0] = 1.0;
-       wind_field[1] = 1.0;
+  Assert(dim > 1, ExcNotImplemented());
+  (void)p; // suppress warnings from p
+  Tensor<1, dim> wind_field;
+  wind_field[0] = 1.0;
+  wind_field[1] = 1.0;
 
-       return wind_field;
+  return wind_field;
 }
 
 // @sect3{The ScratchData and CopyData classes}
@@ -45,781 +43,882 @@ Tensor<1, dim> beta(const Point<dim> &p)
 template <int dim>
 struct ScratchData
 {
-       ScratchData(const Mapping<dim> &mapping, const FiniteElement<dim> &fe,
-                               const Quadrature<dim> &quadrature,
-                               const Quadrature<dim - 1> &quadrature_face,
-                               const UpdateFlags update_flags = update_values | update_gradients | update_quadrature_points | update_JxW_values,
-                               const UpdateFlags interface_update_flags = update_values | update_gradients | update_quadrature_points | update_JxW_values | update_normal_vectors) : fe_values(mapping, fe, quadrature, update_flags), fe_interface_values(mapping, fe, quadrature_face, interface_update_flags)
-       {
-       }
-
-       ScratchData(const ScratchData<dim> &scratch_data) : fe_values(scratch_data.fe_values.get_mapping(),
-                                                                                                                                 scratch_data.fe_values.get_fe(),
-                                                                                                                                 scratch_data.fe_values.get_quadrature(),
-                                                                                                                                 scratch_data.fe_values.get_update_flags()),
-                                                                                                               fe_interface_values(
-                                                                                                                       scratch_data.fe_interface_values.get_mapping(),
-                                                                                                                       scratch_data.fe_interface_values.get_fe(),
-                                                                                                                       scratch_data.fe_interface_values.get_quadrature(),
-                                                                                                                       scratch_data.fe_interface_values.get_update_flags())
-       {
-       }
-
-       FEValues<dim> fe_values;
-       FEInterfaceValues<dim> fe_interface_values;
+  ScratchData(const Mapping<dim>        &mapping,
+              const FiniteElement<dim>  &fe,
+              const Quadrature<dim>     &quadrature,
+              const Quadrature<dim - 1> &quadrature_face,
+              const UpdateFlags          update_flags = update_values |
+                                               update_gradients |
+                                               update_quadrature_points |
+                                               update_JxW_values,
+              const UpdateFlags interface_update_flags =
+                update_values | update_gradients | update_quadrature_points |
+                update_JxW_values | update_normal_vectors)
+    : fe_values(mapping, fe, quadrature, update_flags)
+    , fe_interface_values(mapping, fe, quadrature_face, interface_update_flags)
+  {}
+
+  ScratchData(const ScratchData<dim> &scratch_data)
+    : fe_values(scratch_data.fe_values.get_mapping(),
+                scratch_data.fe_values.get_fe(),
+                scratch_data.fe_values.get_quadrature(),
+                scratch_data.fe_values.get_update_flags())
+    , fe_interface_values(scratch_data.fe_interface_values.get_mapping(),
+                          scratch_data.fe_interface_values.get_fe(),
+                          scratch_data.fe_interface_values.get_quadrature(),
+                          scratch_data.fe_interface_values.get_update_flags())
+  {}
+
+  FEValues<dim>          fe_values;
+  FEInterfaceValues<dim> fe_interface_values;
 };
 
+
+
 struct CopyDataFace
 {
-       FullMatrix<double> cell_matrix;
-       std::vector<types::global_dof_index> joint_dof_indices;
-       std::array<double, 2> values;
-       std::array<unsigned int, 2> cell_indices;
+  FullMatrix<double>                   cell_matrix;
+  std::vector<types::global_dof_index> joint_dof_indices;
+  std::array<double, 2>                values;
+  std::array<unsigned int, 2>          cell_indices;
 };
 
+
+
 struct CopyData
 {
-       FullMatrix<double> cell_matrix;
-       Vector<double> cell_rhs;
-       std::vector<types::global_dof_index> local_dof_indices;
-       std::vector<CopyDataFace> face_data;
-
-       double value;
-       double value_estimator;
-       unsigned int cell_index;
-
-       FullMatrix<double> cell_mass_matrix;
-       Vector<double> cell_mass_rhs;
-
-       template <class Iterator>
-       void reinit(const Iterator &cell, unsigned int dofs_per_cell)
-       {
-               cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
-               cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
-
-               cell_rhs.reinit(dofs_per_cell);
-               cell_mass_rhs.reinit(dofs_per_cell);
-
-               local_dof_indices.resize(dofs_per_cell);
-               cell->get_dof_indices(local_dof_indices);
-       }
+  FullMatrix<double>                   cell_matrix;
+  Vector<double>                       cell_rhs;
+  std::vector<types::global_dof_index> local_dof_indices;
+  std::vector<CopyDataFace>            face_data;
+
+  double       value;
+  double       value_estimator;
+  unsigned int cell_index;
+
+  FullMatrix<double> cell_mass_matrix;
+  Vector<double>     cell_mass_rhs;
+
+  template <class Iterator>
+  void
+  reinit(const Iterator &cell, unsigned int dofs_per_cell)
+  {
+    cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+    cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
+
+    cell_rhs.reinit(dofs_per_cell);
+    cell_mass_rhs.reinit(dofs_per_cell);
+
+    local_dof_indices.resize(dofs_per_cell);
+    cell->get_dof_indices(local_dof_indices);
+  }
 };
 
+
+
 // @sect3{Auxiliary function}
 // This auxiliary function is taken from step-74 and it's used to
 // compute the jump of the finite element function $u_h$ on a face.
 template <int dim>
-void get_function_jump(const FEInterfaceValues<dim> &fe_iv,
-                                          const Vector<double> &solution,
-                                          std::vector<double> &jump)
+void
+get_function_jump(const FEInterfaceValues<dim> &fe_iv,
+                  const Vector<double>         &solution,
+                  std::vector<double>          &jump)
 {
-       const unsigned int n_q = fe_iv.n_quadrature_points;
-       std::array<std::vector<double>, 2> face_values;
-       jump.resize(n_q);
-       for (unsigned int i = 0; i < 2; ++i)
-       {
-               face_values[i].resize(n_q);
-               fe_iv.get_fe_face_values(i).get_function_values(solution,
-                                                                                                               face_values[i]);
-       }
-       for (unsigned int q = 0; q < n_q; ++q)
-               jump[q] = face_values[0][q] - face_values[1][q];
+  const unsigned int                 n_q = fe_iv.n_quadrature_points;
+  std::array<std::vector<double>, 2> face_values;
+  jump.resize(n_q);
+  for (unsigned int i = 0; i < 2; ++i)
+    {
+      face_values[i].resize(n_q);
+      fe_iv.get_fe_face_values(i).get_function_values(solution, face_values[i]);
+    }
+  for (unsigned int q = 0; q < n_q; ++q)
+    jump[q] = face_values[0][q] - face_values[1][q];
 }
 
-// We start with the constructor. The 1 in the constructor call of
-// <code>fe</code> is the polynomial degree.
+
+
 template <int dim>
-AdvectionReaction<dim>::AdvectionReaction() : mapping(),
-                                                                                         dof_handler(triangulation)
+AdvectionReaction<dim>::AdvectionReaction()
+  : mapping()
+  , dof_handler(triangulation)
 {
-
-       add_parameter("Finite element degree", fe_degree);
-       add_parameter("Problem constants", constants);
-       add_parameter("Output filename", output_filename);
-       add_parameter("Use direct solver", use_direct_solver);
-       add_parameter("Number of refinement cycles", n_refinement_cycles);
-       add_parameter("Number of global refinement", n_global_refinements);
-       add_parameter("Refinement", refinement);
-       add_parameter("Exact solution expression", exact_solution_expression);
-       add_parameter("Boundary conditions expression", boundary_conditions_expression);
-       add_parameter("Theta", theta);
-       add_parameter("Advection coefficient expression", advection_coefficient_expression);
-       add_parameter("Right hand side expression", rhs_expression);
-
-       //
-       this->prm.enter_subsection("Error table");
-       error_table.add_parameters(this->prm);
-       this->prm.leave_subsection();
+  Assert(dim > 1, ExcMessage("Not implemented in 1D."));
+  add_parameter("Finite element degree", fe_degree);
+  add_parameter("Problem constants", constants);
+  add_parameter("Output filename", output_filename);
+  add_parameter("Use direct solver", use_direct_solver);
+  add_parameter("Number of refinement cycles", n_refinement_cycles);
+  add_parameter("Number of global refinement", n_global_refinements);
+  add_parameter("Refinement", refinement);
+  add_parameter("Exact solution expression", exact_solution_expression);
+  add_parameter("Boundary conditions expression",
+                boundary_conditions_expression);
+  add_parameter("Theta", theta);
+  add_parameter("Advection coefficient expression",
+                advection_coefficient_expression);
+  add_parameter("Right hand side expression", rhs_expression);
+
+  this->prm.enter_subsection("Error table");
+  error_table.add_parameters(this->prm);
+  this->prm.leave_subsection();
 }
 
+
+
 template <int dim>
-void AdvectionReaction<dim>::initialize_params(const std::string &filename)
+void
+AdvectionReaction<dim>::initialize_params(const std::string &filename)
 {
-
-       ParameterAcceptor::initialize(filename, "last_used_parameters.prm", ParameterHandler::Short);
-       if (theta < 0.0 || theta > 10.0 || std::abs(theta) < 1e-12)
-       {
-               throw(theta_exc("Theta parameter is not in a suitable range: see paper by Brezzi, Marini, Suli for an extended discussion"));
-       }
+  ParameterAcceptor::initialize(filename,
+                                "last_used_parameters.prm",
+                                ParameterHandler::Short);
+  if (theta < 0.0 || theta > 10.0 || std::abs(theta) < 1e-12)
+    {
+      throw(
+        theta_exc("Theta parameter is not in a suitable range: see paper by "
+                  "Brezzi, Marini, Suli for an extended discussion"));
+    }
 }
 
+
+
 template <int dim>
-void AdvectionReaction<dim>::parse_string(const std::string &parameters)
+void
+AdvectionReaction<dim>::parse_string(const std::string &parameters)
 {
-       ParameterAcceptor::prm.parse_input_from_string(parameters);
-       ParameterAcceptor::parse_all_parameters();
+  ParameterAcceptor::prm.parse_input_from_string(parameters);
+  ParameterAcceptor::parse_all_parameters();
 }
 
+
+
 template <int dim>
-void AdvectionReaction<dim>::setup_system()
+void
+AdvectionReaction<dim>::setup_system()
 {
-
-       // first need to distribute the DoFs.
-       if (!fe)
-       {
-               fe = std::make_unique<FE_DGQ<dim>>(fe_degree);
-               const auto vars = dim == 2 ? "x,y" : "x,y,z";
-               exact_solution.initialize(vars, exact_solution_expression, constants);
-               rhs.initialize(vars, rhs_expression, constants);
-               advection_coeff.initialize(vars, advection_coefficient_expression, constants);
-               boundary_conditions.initialize(vars, boundary_conditions_expression, constants);
-       }
-       dof_handler.distribute_dofs(*fe);
-
-       // To build the sparsity pattern for DG discretizations, we can call the
-       // function analogue to DoFTools::make_sparsity_pattern, which is called
-       // DoFTools::make_flux_sparsity_pattern:
-       DynamicSparsityPattern dsp(dof_handler.n_dofs());
-       DoFTools::make_flux_sparsity_pattern(dof_handler, dsp); //DG sparsity pattern generator
-       sparsity_pattern.copy_from(dsp);
-
-       // Finally, we set up the structure of all components of the linear system.
-       system_matrix.reinit(sparsity_pattern);
-       solution.reinit(dof_handler.n_dofs());
-       right_hand_side.reinit(dof_handler.n_dofs());
+  // first need to distribute the DoFs.
+  if (!fe)
+    {
+      fe              = std::make_unique<FE_DGQ<dim>>(fe_degree);
+      const auto vars = dim == 2 ? "x,y" : "x,y,z";
+      exact_solution.initialize(vars, exact_solution_expression, constants);
+      rhs.initialize(vars, rhs_expression, constants);
+      advection_coeff.initialize(vars,
+                                 advection_coefficient_expression,
+                                 constants);
+      boundary_conditions.initialize(vars,
+                                     boundary_conditions_expression,
+                                     constants);
+    }
+  dof_handler.distribute_dofs(*fe);
+
+  // To build the sparsity pattern for DG discretizations, we can call the
+  // function analogue to DoFTools::make_sparsity_pattern, which is called
+  // DoFTools::make_flux_sparsity_pattern:
+  DynamicSparsityPattern dsp(dof_handler.n_dofs());
+  DoFTools::make_flux_sparsity_pattern(dof_handler,
+                                       dsp); // DG sparsity pattern generator
+  sparsity_pattern.copy_from(dsp);
+
+  // Finally, we set up the structure of all components of the linear system.
+  system_matrix.reinit(sparsity_pattern);
+  solution.reinit(dof_handler.n_dofs());
+  right_hand_side.reinit(dof_handler.n_dofs());
 }
 
-//in the call to  MeshWorker::mesh_loop() we only need to specify what should happen on
-// each cell, each boundary face, and each interior face. These three tasks
-// are handled by the lambda functions inside the function below.
+
+
+// in the call to  MeshWorker::mesh_loop() we only need to specify what should
+// happen on
+//  each cell, each boundary face, and each interior face. These three tasks
+//  are handled by the lambda functions inside the function below.
 
 template <int dim>
-void AdvectionReaction<dim>::assemble_system()
+void
+AdvectionReaction<dim>::assemble_system()
 {
-
-       using Iterator = typename DoFHandler<dim>::active_cell_iterator;
-
-       const QGauss<dim> quadrature = fe->tensor_degree() + 1;
-       const QGauss<dim - 1> quadrature_face = fe->tensor_degree() + 1;
-
-       // This is the function that will be executed for each cell.
-       const auto cell_worker = [&](const Iterator &cell,
-                                                                ScratchData<dim> &scratch_data, CopyData &copy_data)
-       {
-               FEValues<dim> fe_values_continuous(*fe,
-                                                                                  quadrature,
-                                                                                  update_values | update_gradients |
-                                                                                          update_quadrature_points | update_JxW_values);
-
-               const unsigned int n_dofs = scratch_data.fe_values.get_fe().n_dofs_per_cell();
-               copy_data.reinit(cell, n_dofs);
-               scratch_data.fe_values.reinit(cell);
-
-               const auto &q_points = scratch_data.fe_values.get_quadrature_points();
-
-               const FEValues<dim> &fe_v = scratch_data.fe_values;
-               const std::vector<double> &JxW = fe_v.get_JxW_values();
-
-               for (unsigned int point = 0; point < fe_v.n_quadrature_points;
-                        ++point)
-               {
-                       auto beta_q = beta(q_points[point]);
-                       for (unsigned int i = 0; i < n_dofs; ++i)
-                       {
-                               for (unsigned int j = 0; j < n_dofs; ++j)
-                               {
-                                       copy_data.cell_matrix(i, j) += (-beta_q                                                  // -\beta
-                                                                                                               * fe_v.shape_grad(i, point)      // \nabla \phi_i
-                                                                                                               * fe_v.shape_value(j, point) // \phi_j
-                                                                                                       +
-                                                                                                       advection_coeff.value(q_points[point]) * //gamma
-                                                                                                               fe_v.shape_value(i, point)                       //phi_i
-                                                                                                               * fe_v.shape_value(j, point)             //phi_j
-                                                                                                       ) *
-                                                                                                  JxW[point]; // dx
-                               }
-                               copy_data.cell_rhs(i) +=
-                                       rhs.value(q_points[point])       // f(x_q)
-                                       * fe_v.shape_value(i, point) //phi_i(x_q)
-                                       * JxW[point];                            //dx
-                       }
-               }
-       };
-
-       // This is the function called for boundary faces and consists of a normal
-       // integration using FEFaceValues. New is the logic to decide if the term
-       // goes into the system matrix (outflow) or the right-hand side (inflow).
-       const auto boundary_worker = [&](const Iterator &cell,
-                                                                        const unsigned int &face_no, ScratchData<dim> &scratch_data,
-                                                                        CopyData &copy_data)
-       {
-               scratch_data.fe_interface_values.reinit(cell, face_no);
-               const FEFaceValuesBase<dim> &fe_face =
-                       scratch_data.fe_interface_values.get_fe_face_values(0);
-
-               const auto &q_points = fe_face.get_quadrature_points();
-
-               const unsigned int n_facet_dofs = fe_face.get_fe().n_dofs_per_cell();
-               const std::vector<double> &JxW = fe_face.get_JxW_values();
-               const std::vector<Tensor<1, dim>> &normals =
-                       fe_face.get_normal_vectors();
-
-               std::vector<double> g(q_points.size());
-               exact_solution.value_list(q_points, g);
-
-               for (unsigned int point = 0; point < q_points.size(); ++point)
-               {
-                       const double beta_dot_n = beta(q_points[point]) * normals[point];
-
-                       if (beta_dot_n > 0)
-                       {
-                               for (unsigned int i = 0; i < n_facet_dofs; ++i)
-                                       for (unsigned int j = 0; j < n_facet_dofs; ++j)
-                                               copy_data.cell_matrix(i, j) += fe_face.shape_value(i,
-                                                                                                                                                  point)          // \phi_i
-                                                                                                          * fe_face.shape_value(j, point) // \phi_j
-                                                                                                          * beta_dot_n                                    // \beta . n
-                                                                                                          * JxW[point];                                   // dx
-                       }
-                       else
-                               for (unsigned int i = 0; i < n_facet_dofs; ++i)
-                                       copy_data.cell_rhs(i) += -fe_face.shape_value(i, point) // \phi_i
-                                                                                        * g[point]                                             // g*/
-                                                                                        * beta_dot_n                                   // \beta . n
-                                                                                        * JxW[point];                                  // dx
-               }
-       };
-
-       // This is the function called on interior faces. The arguments specify
-       // cells, face and subface indices (for adaptive refinement). We just pass
-       // them along to the reinit() function of FEInterfaceValues.
-       const auto face_worker = [&](const Iterator &cell, const unsigned int &f,
-                                                                const unsigned int &sf, const Iterator &ncell,
-                                                                const unsigned int &nf, const unsigned int &nsf,
-                                                                ScratchData<dim> &scratch_data, CopyData &copy_data)
-       {
-               FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
-               fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
-               const auto &q_points = fe_iv.get_quadrature_points();
-
-               copy_data.face_data.emplace_back();
-               CopyDataFace &copy_data_face = copy_data.face_data.back();
-
-               const unsigned int n_dofs = fe_iv.n_current_interface_dofs();
-               copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
-
-               copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
-
-               const std::vector<double> &JxW = fe_iv.get_JxW_values();
-               const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
-
-               for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
-               {
-                       const double beta_dot_n = beta(q_points[qpoint]) * normals[qpoint];
-                       for (unsigned int i = 0; i < n_dofs; ++i)
-                       {
-                               for (unsigned int j = 0; j < n_dofs; ++j)
-                               {
-                                       copy_data_face.cell_matrix(i, j) += (beta(q_points[qpoint]) * normals[qpoint] * fe_iv.average_of_shape_values(j, qpoint) * fe_iv.jump_in_shape_values(i, qpoint) +
-                                                                                                                theta * std::abs(beta_dot_n) * fe_iv.jump_in_shape_values(j, qpoint) * fe_iv.jump_in_shape_values(i, qpoint)) *
-                                                                                                               JxW[qpoint];
-                               }
-                       }
-               }
-       };
-
-       // The following lambda function will handle copying the data from the
-       // cell and face assembly into the global matrix and right-hand side.
-       //
-       // While we would not need an AffineConstraints object, because there are
-       // no hanging node constraints in DG discretizations, we use an empty
-       // object here as this allows us to use its `copy_local_to_global`
-       // functionality.
-       const AffineConstraints<double> constraints;
-
-       const auto copier = [&](const CopyData &c)
-       {
-               constraints.distribute_local_to_global(c.cell_matrix, c.cell_rhs,
-                                                                                          c.local_dof_indices, system_matrix, right_hand_side);
-
-               for (auto &cdf : c.face_data)
-               {
-                       constraints.distribute_local_to_global(cdf.cell_matrix,
-                                                                                                  cdf.joint_dof_indices, system_matrix);
-               }
-       };
-
-       ScratchData<dim> scratch_data(mapping, *fe, quadrature, quadrature_face);
-       CopyData copy_data;
-
-       // Here, we finally handle the assembly. We pass in ScratchData and
-       // CopyData objects, the lambda functions from above, an specify that we
-       // want to assemble interior faces once.
-       MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(),
-                                                 cell_worker, copier, scratch_data, copy_data,
-                                                 MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces | MeshWorker::assemble_own_interior_faces_once,
-                                                 boundary_worker, face_worker);
+  using Iterator = typename DoFHandler<dim>::active_cell_iterator;
+
+  const QGauss<dim>     quadrature      = fe->tensor_degree() + 1;
+  const QGauss<dim - 1> quadrature_face = fe->tensor_degree() + 1;
+
+  // This is the function that will be executed for each cell.
+  const auto cell_worker = [&](const Iterator   &cell,
+                               ScratchData<dim> &scratch_data,
+                               CopyData         &copy_data) {
+    FEValues<dim> fe_values_continuous(*fe,
+                                       quadrature,
+                                       update_values | update_gradients |
+                                         update_quadrature_points |
+                                         update_JxW_values);
+
+    const unsigned int n_dofs =
+      scratch_data.fe_values.get_fe().n_dofs_per_cell();
+    copy_data.reinit(cell, n_dofs);
+    scratch_data.fe_values.reinit(cell);
+
+    const auto &q_points = scratch_data.fe_values.get_quadrature_points();
+
+    const FEValues<dim>       &fe_v = scratch_data.fe_values;
+    const std::vector<double> &JxW  = fe_v.get_JxW_values();
+
+    for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
+      {
+        auto beta_q = beta(q_points[point]);
+        for (unsigned int i = 0; i < n_dofs; ++i)
+          {
+            for (unsigned int j = 0; j < n_dofs; ++j)
+              {
+                copy_data.cell_matrix(i, j) +=
+                  (-beta_q                                    // -\beta
+                     * fe_v.shape_grad(i, point)              // \nabla \phi_i
+                     * fe_v.shape_value(j, point)             // \phi_j
+                   + advection_coeff.value(q_points[point]) * // gamma
+                       fe_v.shape_value(i, point)             // phi_i
+                       * fe_v.shape_value(j, point)           // phi_j
+                   ) *
+                  JxW[point]; // dx
+              }
+            copy_data.cell_rhs(i) += rhs.value(q_points[point])   // f(x_q)
+                                     * fe_v.shape_value(i, point) // phi_i(x_q)
+                                     * JxW[point];                // dx
+          }
+      }
+  };
+
+  // This is the function called for boundary faces and consists of a normal
+  // integration using FEFaceValues. New is the logic to decide if the term
+  // goes into the system matrix (outflow) or the right-hand side (inflow).
+  const auto boundary_worker = [&](const Iterator     &cell,
+                                   const unsigned int &face_no,
+                                   ScratchData<dim>   &scratch_data,
+                                   CopyData           &copy_data) {
+    scratch_data.fe_interface_values.reinit(cell, face_no);
+    const FEFaceValuesBase<dim> &fe_face =
+      scratch_data.fe_interface_values.get_fe_face_values(0);
+
+    const auto &q_points = fe_face.get_quadrature_points();
+
+    const unsigned int n_facet_dofs        = fe_face.get_fe().n_dofs_per_cell();
+    const std::vector<double>         &JxW = fe_face.get_JxW_values();
+    const std::vector<Tensor<1, dim>> &normals = fe_face.get_normal_vectors();
+
+    std::vector<double> g(q_points.size());
+    exact_solution.value_list(q_points, g);
+
+    for (unsigned int point = 0; point < q_points.size(); ++point)
+      {
+        const double beta_dot_n = beta(q_points[point]) * normals[point];
+
+        if (beta_dot_n > 0)
+          {
+            for (unsigned int i = 0; i < n_facet_dofs; ++i)
+              for (unsigned int j = 0; j < n_facet_dofs; ++j)
+                copy_data.cell_matrix(i, j) +=
+                  fe_face.shape_value(i,
+                                      point)      // \phi_i
+                  * fe_face.shape_value(j, point) // \phi_j
+                  * beta_dot_n                    // \beta . n
+                  * JxW[point];                   // dx
+          }
+        else
+          for (unsigned int i = 0; i < n_facet_dofs; ++i)
+            copy_data.cell_rhs(i) += -fe_face.shape_value(i, point) // \phi_i
+                                     * g[point]                     // g*/
+                                     * beta_dot_n                   // \beta . n
+                                     * JxW[point];                  // dx
+      }
+  };
+
+  // This is the function called on interior faces. The arguments specify
+  // cells, face and subface indices (for adaptive refinement). We just pass
+  // them along to the reinit() function of FEInterfaceValues.
+  const auto face_worker = [&](const Iterator     &cell,
+                               const unsigned int &f,
+                               const unsigned int &sf,
+                               const Iterator     &ncell,
+                               const unsigned int &nf,
+                               const unsigned int &nsf,
+                               ScratchData<dim>   &scratch_data,
+                               CopyData           &copy_data) {
+    FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
+    fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
+    const auto &q_points = fe_iv.get_quadrature_points();
+
+    copy_data.face_data.emplace_back();
+    CopyDataFace &copy_data_face = copy_data.face_data.back();
+
+    const unsigned int n_dofs        = fe_iv.n_current_interface_dofs();
+    copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
+
+    copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
+
+    const std::vector<double>         &JxW     = fe_iv.get_JxW_values();
+    const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
+
+    for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+      {
+        const double beta_dot_n = beta(q_points[qpoint]) * normals[qpoint];
+        for (unsigned int i = 0; i < n_dofs; ++i)
+          {
+            for (unsigned int j = 0; j < n_dofs; ++j)
+              {
+                copy_data_face.cell_matrix(i, j) +=
+                  (beta(q_points[qpoint]) * normals[qpoint] *
+                     fe_iv.average_of_shape_values(j, qpoint) *
+                     fe_iv.jump_in_shape_values(i, qpoint) +
+                   theta * std::abs(beta_dot_n) *
+                     fe_iv.jump_in_shape_values(j, qpoint) *
+                     fe_iv.jump_in_shape_values(i, qpoint)) *
+                  JxW[qpoint];
+              }
+          }
+      }
+  };
+
+  // The following lambda function will handle copying the data from the
+  // cell and face assembly into the global matrix and right-hand side.
+  //
+  // While we would not need an AffineConstraints object, because there are
+  // no hanging node constraints in DG discretizations, we use an empty
+  // object here as this allows us to use its `copy_local_to_global`
+  // functionality.
+  const AffineConstraints<double> constraints;
+
+  const auto copier = [&](const CopyData &c) {
+    constraints.distribute_local_to_global(c.cell_matrix,
+                                           c.cell_rhs,
+                                           c.local_dof_indices,
+                                           system_matrix,
+                                           right_hand_side);
+
+    for (auto &cdf : c.face_data)
+      {
+        constraints.distribute_local_to_global(cdf.cell_matrix,
+                                               cdf.joint_dof_indices,
+                                               system_matrix);
+      }
+  };
+
+  ScratchData<dim> scratch_data(mapping, *fe, quadrature, quadrature_face);
+  CopyData         copy_data;
+
+  // Here, we finally handle the assembly. We pass in ScratchData and
+  // CopyData objects, the lambda functions from above, an specify that we
+  // want to assemble interior faces once.
+  MeshWorker::mesh_loop(dof_handler.begin_active(),
+                        dof_handler.end(),
+                        cell_worker,
+                        copier,
+                        scratch_data,
+                        copy_data,
+                        MeshWorker::assemble_own_cells |
+                          MeshWorker::assemble_boundary_faces |
+                          MeshWorker::assemble_own_interior_faces_once,
+                        boundary_worker,
+                        face_worker);
 }
 
+
+
 template <int dim>
-void AdvectionReaction<dim>::solve()
+void
+AdvectionReaction<dim>::solve()
 {
-
-       if (use_direct_solver)
-       {
-
-               SparseDirectUMFPACK system_matrix_inverse;
-               system_matrix_inverse.initialize(system_matrix);
-               system_matrix_inverse.vmult(solution, right_hand_side);
-       }
-       else
-       {
-               // Here we have a classic iterative solver, as done in many tutorials:
-               SolverControl solver_control(1000, 1e-15);
-               SolverRichardson<Vector<double>> solver(solver_control);
-               PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
-               preconditioner.initialize(system_matrix, fe->n_dofs_per_cell());
-               solver.solve(system_matrix, solution, right_hand_side, preconditioner);
-               std::cout << "  Solver converged in " << solver_control.last_step()
-                                 << " iterations." << std::endl;
-       }
+  if (use_direct_solver)
+    {
+      SparseDirectUMFPACK system_matrix_inverse;
+      system_matrix_inverse.initialize(system_matrix);
+      system_matrix_inverse.vmult(solution, right_hand_side);
+    }
+  else
+    {
+      // Here we have a classic iterative solver, as done in many tutorials:
+      SolverControl                               solver_control(1000, 1e-15);
+      SolverRichardson<Vector<double>>            solver(solver_control);
+      PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
+      preconditioner.initialize(system_matrix, fe->n_dofs_per_cell());
+      solver.solve(system_matrix, solution, right_hand_side, preconditioner);
+      std::cout << "  Solver converged in " << solver_control.last_step()
+                << " iterations." << std::endl;
+    }
 }
 
+
+
 // @sect3{Mesh refinement}
-// We refine the grid according the proposed estimator or with an approximation to the gradient of the solution.
-// The first option is the default one (you can see it in the header file)
+// We refine the grid according the proposed estimator or with an approximation
+// to the gradient of the solution. The first option is the default one (you can
+// see it in the header file)
 template <int dim>
-void AdvectionReaction<dim>::refine_grid()
+void
+AdvectionReaction<dim>::refine_grid()
 {
-
-       if (refinement == "residual")
-       {
-
-               //If the `refinement` string is `"residual"`, then we first compute the local projection
-               compute_local_projection_and_estimate();
-               //We then set the refinement fraction and as usual we execute the refinement.
-               const double refinement_fraction = 0.6;
-               GridRefinement::refine_and_coarsen_fixed_fraction(triangulation, error_indicator_per_cell, refinement_fraction, 0.0);
-               triangulation.execute_coarsening_and_refinement();
-       }
-       else if (refinement == "gradient")
-       {
-
-               Vector<float> gradient_indicator(triangulation.n_active_cells());
-
-               // Now the approximate gradients are computed
-               DerivativeApproximation::approximate_gradient(mapping, dof_handler,
-                                                                                                         solution, gradient_indicator);
-
-               // and they are cell-wise scaled by the factor $h^{1+d/2}$
-               unsigned int cell_no = 0;
-               for (const auto &cell : dof_handler.active_cell_iterators())
-                       gradient_indicator(cell_no++) *= std::pow(cell->diameter(),
-                                                                                                         1 + 1.0 * dim / 2);
-
-               // Finally they serve as refinement indicator.
-               GridRefinement::refine_and_coarsen_fixed_fraction(triangulation,
-                                                                                                                 gradient_indicator, 0.25, 0.0);
-
-               triangulation.execute_coarsening_and_refinement();
-               std::cout << gradient_indicator.l2_norm() << '\n';
-       }
-       else if (refinement == "global")
-       {
-               triangulation.refine_global(1); //just for testing on uniformly refined meshes
-       }
-       else
-       {
-               Assert(false, ExcInternalError());
-       }
+  if (refinement == "residual")
+    {
+      // If the `refinement` string is `"residual"`, then we first compute the
+      // local projection
+      compute_local_projection_and_estimate();
+      // We then set the refinement fraction and as usual we execute the
+      // refinement.
+      const double refinement_fraction = 0.6;
+      GridRefinement::refine_and_coarsen_fixed_fraction(
+        triangulation, error_indicator_per_cell, refinement_fraction, 0.0);
+      triangulation.execute_coarsening_and_refinement();
+    }
+  else if (refinement == "gradient")
+    {
+      Vector<float> gradient_indicator(triangulation.n_active_cells());
+
+      // Now the approximate gradients are computed
+      DerivativeApproximation::approximate_gradient(mapping,
+                                                    dof_handler,
+                                                    solution,
+                                                    gradient_indicator);
+
+      // and they are cell-wise scaled by the factor $h^{1+d/2}$
+      unsigned int cell_no = 0;
+      for (const auto &cell : dof_handler.active_cell_iterators())
+        gradient_indicator(cell_no++) *=
+          std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
+
+      // Finally they serve as refinement indicator.
+      GridRefinement::refine_and_coarsen_fixed_fraction(triangulation,
+                                                        gradient_indicator,
+                                                        0.25,
+                                                        0.0);
+
+      triangulation.execute_coarsening_and_refinement();
+      std::cout << gradient_indicator.l2_norm() << '\n';
+    }
+  else if (refinement == "global")
+    {
+      triangulation.refine_global(
+        1); // just for testing on uniformly refined meshes
+    }
+  else
+    {
+      Assert(false, ExcInternalError());
+    }
 }
+
+
+
 // The output of this program consists of a vtk file of the adaptively
 // refined grids and the numerical solutions.
 template <int dim>
-void AdvectionReaction<dim>::output_results(const unsigned int cycle) const
+void
+AdvectionReaction<dim>::output_results(const unsigned int cycle) const
 {
-       const std::string filename = "solution-" + std::to_string(cycle) + ".vtk";
-       std::cout << "  Writing solution to <" << filename << ">" << std::endl;
-       std::ofstream output(filename);
-
-       DataOut<dim> data_out;
-       data_out.attach_dof_handler(dof_handler);
-       data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
-
-       data_out.build_patches(mapping);
-
-       data_out.write_vtk(output);
+  const std::string filename = "solution-" + std::to_string(cycle) + ".vtk";
+  std::cout << "  Writing solution to <" << filename << ">" << std::endl;
+  std::ofstream output(filename);
+
+  DataOut<dim> data_out;
+  data_out.attach_dof_handler(dof_handler);
+  data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
+  data_out.build_patches(mapping);
+  data_out.write_vtk(output);
 }
 
 template <int dim>
-void AdvectionReaction<dim>::compute_error()
+void
+AdvectionReaction<dim>::compute_error()
 {
-       error_table.error_from_exact(mapping, dof_handler, solution, exact_solution); //be careful: a FD approximation of the gradient is used to compute the H^1 norm if you're not relying on SymbolicFunction class
-                                                                                                                                                                 //    error_table.error_from_exact(mapping, dof_handler, solution, Solution<dim>()); //provided that Solution<dim> implements the Gradient function
+  error_table.error_from_exact(
+    mapping,
+    dof_handler,
+    solution,
+    exact_solution); // be careful: a FD approximation of the gradient is used
+                     // to compute the H^1 norm if Solution<dim> doesn't
+                     // implements the Gradient function
 }
 
+
+
 // @sect3{Compute the energy norm}
-// The energy norm is defined as $ |||\cdot ||| = \Bigl(||\cdot||_{0,\Omega}^2 + \sum_{F \in \mathbb{F}}||c_F^{\frac{1}{2}}[\cdot] ||_{0,F}^2 \Bigr)^{\frac{1}{2}}$
-// Notice that in the current case we have $c_f = \frac{|b \cdot n|}{2}$
-// Like in the assembly, all the contributions are handled separately by using ScratchData and CopyData objects.
+// The energy norm is defined as $ |||\cdot ||| = \Bigl(||\cdot||_{0,\Omega}^2 +
+// \sum_{F \in \mathbb{F}}||c_F^{\frac{1}{2}}[\cdot] ||_{0,F}^2
+// \Bigr)^{\frac{1}{2}}$ Notice that in the current case we have $c_f = \frac{|b
+// \cdot n|}{2}$ Like in the assembly, all the contributions are handled
+// separately by using ScratchData and CopyData objects.
 template <int dim>
-double AdvectionReaction<dim>::compute_energy_norm()
+double
+AdvectionReaction<dim>::compute_energy_norm()
 {
-
-       energy_norm_square_per_cell.reinit(triangulation.n_active_cells());
-
-       using Iterator = typename DoFHandler<dim>::active_cell_iterator;
-
-       // We start off by adding cell contributions
-       const auto cell_worker = [&](const Iterator &cell,
-                                                                ScratchData<dim> &scratch_data, CopyData &copy_data)
-       {
-               const unsigned int n_dofs =
-                       scratch_data.fe_values.get_fe().n_dofs_per_cell();
-               copy_data.reinit(cell, n_dofs);
-               scratch_data.fe_values.reinit(cell);
-
-               copy_data.cell_index = cell->active_cell_index();
-
-               const auto &q_points = scratch_data.fe_values.get_quadrature_points();
-               const FEValues<dim> &fe_v = scratch_data.fe_values;
-               const std::vector<double> &JxW = fe_v.get_JxW_values();
-
-               double error_square_norm{0.0};
-               std::vector<double> sol_u(fe_v.n_quadrature_points);
-               fe_v.get_function_values(solution, sol_u);
-
-               for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
-               {
-                       const double diff = (sol_u[point] - exact_solution.value(q_points[point]));
-                       error_square_norm += diff * diff * JxW[point];
-               }
-               copy_data.value = error_square_norm;
-       };
-
-       // Here we add contributions coming from the internal faces
-       const auto face_worker = [&](const Iterator &cell,
-                                                                const unsigned int &f,
-                                                                const unsigned int &sf,
-                                                                const Iterator &ncell,
-                                                                const unsigned int &nf,
-                                                                const unsigned int &nsf,
-                                                                ScratchData<dim> &scratch_data,
-                                                                CopyData &copy_data)
-       {
-               FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
-               fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
-
-               copy_data.face_data.emplace_back();
-               CopyDataFace &copy_data_face = copy_data.face_data.back();
-               copy_data_face.cell_indices[0] = cell->active_cell_index();
-               copy_data_face.cell_indices[1] = ncell->active_cell_index();
-
-               const auto &q_points = fe_iv.get_quadrature_points();
-               const unsigned n_q_points = q_points.size();
-               const std::vector<double> &JxW = fe_iv.get_JxW_values();
-               std::vector<double> g(n_q_points);
-
-               std::vector<double> jump(n_q_points);
-               get_function_jump(fe_iv, solution, jump);
-
-               const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
-
-               double error_jump_square{0.0};
-               for (unsigned int point = 0; point < n_q_points; ++point)
-               {
-                       const double beta_dot_n = theta * std::abs(beta(q_points[point]) * normals[point]);
-                       error_jump_square += beta_dot_n * jump[point] * jump[point] * JxW[point];
-               }
-
-               copy_data.value = error_jump_square;
-       };
-
-       // Finally, we add the boundary contributions
-       const auto boundary_worker = [&](const Iterator &cell,
-                                                                        const unsigned int &face_no,
-                                                                        ScratchData<dim> &scratch_data,
-                                                                        CopyData &copy_data)
-       {
-               scratch_data.fe_interface_values.reinit(cell, face_no);
-               const FEFaceValuesBase<dim> &fe_fv = scratch_data.fe_interface_values.get_fe_face_values(0);
-               const auto &q_points = fe_fv.get_quadrature_points();
-               const unsigned n_q_points = q_points.size();
-               const std::vector<double> &JxW = fe_fv.get_JxW_values();
-
-               std::vector<double> g(n_q_points);
-
-               std::vector<double> sol_u(n_q_points);
-               fe_fv.get_function_values(solution, sol_u);
-
-               const std::vector<Tensor<1, dim>> &normals = fe_fv.get_normal_vectors();
-
-               double difference_norm_square = 0.;
-               for (unsigned int point = 0; point < q_points.size(); ++point)
-               {
-                       const double beta_dot_n = theta * std::abs(beta(q_points[point]) * normals[point]);
-                       const double diff = (boundary_conditions.value(q_points[point]) - sol_u[point]);
-                       difference_norm_square += beta_dot_n * diff * diff * JxW[point];
-               }
-               copy_data.value = difference_norm_square;
-       };
-
-       const auto copier = [&](const auto &copy_data)
-       {
-               if (copy_data.cell_index != numbers::invalid_unsigned_int)
-               {
-                       energy_norm_square_per_cell[copy_data.cell_index] += copy_data.value;
-               }
-               for (auto &cdf : copy_data.face_data)
-                       for (unsigned int j = 0; j < 2; ++j)
-                               energy_norm_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
-       };
-
-       ScratchData<dim> scratch_data(mapping, *fe, QGauss<dim>{fe->tensor_degree() + 1},
-                                                                 QGauss<dim - 1>{fe->tensor_degree() + 1});
-
-       CopyData copy_data;
-
-       MeshWorker::mesh_loop(dof_handler.begin_active(),
-                                                 dof_handler.end(),
-                                                 cell_worker,
-                                                 copier,
-                                                 scratch_data,
-                                                 copy_data,
-                                                 MeshWorker::assemble_own_cells |
-                                                         MeshWorker::assemble_own_interior_faces_once |
-                                                         MeshWorker::assemble_boundary_faces,
-                                                 boundary_worker,
-                                                 face_worker);
-
-       const double energy_error = std::sqrt(energy_norm_square_per_cell.l1_norm());
-       return energy_error;
+  energy_norm_square_per_cell.reinit(triangulation.n_active_cells());
+
+  using Iterator = typename DoFHandler<dim>::active_cell_iterator;
+
+  // We start off by adding cell contributions
+  const auto cell_worker = [&](const Iterator   &cell,
+                               ScratchData<dim> &scratch_data,
+                               CopyData         &copy_data) {
+    const unsigned int n_dofs =
+      scratch_data.fe_values.get_fe().n_dofs_per_cell();
+    copy_data.reinit(cell, n_dofs);
+    scratch_data.fe_values.reinit(cell);
+
+    copy_data.cell_index = cell->active_cell_index();
+
+    const auto &q_points = scratch_data.fe_values.get_quadrature_points();
+    const FEValues<dim>       &fe_v = scratch_data.fe_values;
+    const std::vector<double> &JxW  = fe_v.get_JxW_values();
+
+    double              error_square_norm{0.0};
+    std::vector<double> sol_u(fe_v.n_quadrature_points);
+    fe_v.get_function_values(solution, sol_u);
+
+    for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
+      {
+        const double diff =
+          (sol_u[point] - exact_solution.value(q_points[point]));
+        error_square_norm += diff * diff * JxW[point];
+      }
+    copy_data.value = error_square_norm;
+  };
+
+  // Here we add contributions coming from the internal faces
+  const auto face_worker = [&](const Iterator     &cell,
+                               const unsigned int &f,
+                               const unsigned int &sf,
+                               const Iterator     &ncell,
+                               const unsigned int &nf,
+                               const unsigned int &nsf,
+                               ScratchData<dim>   &scratch_data,
+                               CopyData           &copy_data) {
+    FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
+    fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
+
+    copy_data.face_data.emplace_back();
+    CopyDataFace &copy_data_face   = copy_data.face_data.back();
+    copy_data_face.cell_indices[0] = cell->active_cell_index();
+    copy_data_face.cell_indices[1] = ncell->active_cell_index();
+
+    const auto                &q_points   = fe_iv.get_quadrature_points();
+    const unsigned             n_q_points = q_points.size();
+    const std::vector<double> &JxW        = fe_iv.get_JxW_values();
+    std::vector<double>        g(n_q_points);
+
+    std::vector<double> jump(n_q_points);
+    get_function_jump(fe_iv, solution, jump);
+
+    const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
+
+    double error_jump_square{0.0};
+    for (unsigned int point = 0; point < n_q_points; ++point)
+      {
+        const double beta_dot_n =
+          theta * std::abs(beta(q_points[point]) * normals[point]);
+        error_jump_square +=
+          beta_dot_n * jump[point] * jump[point] * JxW[point];
+      }
+
+    copy_data.value = error_jump_square;
+  };
+
+  // Finally, we add the boundary contributions
+  const auto boundary_worker = [&](const Iterator     &cell,
+                                   const unsigned int &face_no,
+                                   ScratchData<dim>   &scratch_data,
+                                   CopyData           &copy_data) {
+    scratch_data.fe_interface_values.reinit(cell, face_no);
+    const FEFaceValuesBase<dim> &fe_fv =
+      scratch_data.fe_interface_values.get_fe_face_values(0);
+    const auto                &q_points   = fe_fv.get_quadrature_points();
+    const unsigned             n_q_points = q_points.size();
+    const std::vector<double> &JxW        = fe_fv.get_JxW_values();
+
+    std::vector<double> g(n_q_points);
+
+    std::vector<double> sol_u(n_q_points);
+    fe_fv.get_function_values(solution, sol_u);
+
+    const std::vector<Tensor<1, dim>> &normals = fe_fv.get_normal_vectors();
+
+    double difference_norm_square = 0.;
+    for (unsigned int point = 0; point < q_points.size(); ++point)
+      {
+        const double beta_dot_n =
+          theta * std::abs(beta(q_points[point]) * normals[point]);
+        const double diff =
+          (boundary_conditions.value(q_points[point]) - sol_u[point]);
+        difference_norm_square += beta_dot_n * diff * diff * JxW[point];
+      }
+    copy_data.value = difference_norm_square;
+  };
+
+  const auto copier = [&](const auto &copy_data) {
+    if (copy_data.cell_index != numbers::invalid_unsigned_int)
+      {
+        energy_norm_square_per_cell[copy_data.cell_index] += copy_data.value;
+      }
+    for (auto &cdf : copy_data.face_data)
+      for (unsigned int j = 0; j < 2; ++j)
+        energy_norm_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
+  };
+
+  ScratchData<dim> scratch_data(mapping,
+                                *fe,
+                                QGauss<dim>{fe->tensor_degree() + 1},
+                                QGauss<dim - 1>{fe->tensor_degree() + 1});
+
+  CopyData copy_data;
+
+  MeshWorker::mesh_loop(dof_handler.begin_active(),
+                        dof_handler.end(),
+                        cell_worker,
+                        copier,
+                        scratch_data,
+                        copy_data,
+                        MeshWorker::assemble_own_cells |
+                          MeshWorker::assemble_own_interior_faces_once |
+                          MeshWorker::assemble_boundary_faces,
+                        boundary_worker,
+                        face_worker);
+
+  const double energy_error = std::sqrt(energy_norm_square_per_cell.l1_norm());
+  return energy_error;
 }
 
+
+
 // @sect3{Computing the estimator}
-// In the estimator, we have to compute the term $||f- c u_h - \Pi(f- c u_h)||_{T}^{2}$ over a generic cell $T$. To achieve this, we first need to
-// compute the projection involving the finite element function $u_h$. Using the definition of orthogonal projection, we're required to solve cellwise
-// $(v_h,f-c u_h)_T = (v_h,\Pi)_T \qquad \forall v_h \in V_h$ for $\Pi$, which means that we have to build a mass-matrix on each cell.
-// Once we have the projection, which is a finite element function, we can add its contribution in the <code>cell_worker</code> lambda.
-// As done in step-74, the square of the error indicator is computed.
+// In the estimator, we have to compute the term $||f- c u_h - \Pi(f- c
+// u_h)||_{T}^{2}$ over a generic cell $T$. To achieve this, we first need to
+// compute the projection involving the finite element function $u_h$. Using the
+// definition of orthogonal projection, we're required to solve cellwise
+// $(v_h,f-c u_h)_T = (v_h,\Pi)_T \qquad \forall v_h \in V_h$ for $\Pi$, which
+// means that we have to build a mass-matrix on each cell. Once we have the
+// projection, which is a finite element function, we can add its contribution
+// in the <code>cell_worker</code> lambda. As done in step-74, the square of the
+// error indicator is computed.
 //
 template <int dim>
-void AdvectionReaction<dim>::compute_local_projection_and_estimate()
+void
+AdvectionReaction<dim>::compute_local_projection_and_estimate()
 {
-
-       // Compute the term $||f-c u_h - \Pi(f- cu_h)||_T^2$
-       using Iterator = typename DoFHandler<dim>::active_cell_iterator;
-       error_indicator_per_cell.reinit(triangulation.n_active_cells());
-
-       const auto cell_worker = [&](const Iterator &cell,
-                                                                ScratchData<dim> &scratch_data, CopyData &copy_data)
-       {
-               const unsigned int n_dofs = scratch_data.fe_values.get_fe().n_dofs_per_cell();
-
-               copy_data.reinit(cell, n_dofs);
-               scratch_data.fe_values.reinit(cell);
-               copy_data.cell_index = cell->active_cell_index();
-
-               const auto &q_points = scratch_data.fe_values.get_quadrature_points();
-               const unsigned n_q_points = q_points.size();
-
-               const FEValues<dim> &fe_v = scratch_data.fe_values;
-               const std::vector<double> &JxW = fe_v.get_JxW_values();
-
-               std::vector<double> sol_u_at_quadrature_points(fe_v.n_quadrature_points);
-               fe_v.get_function_values(solution, sol_u_at_quadrature_points);
-
-               //Compute local L^2 projection of  $f- c u_h$ over the local finite element space
-               for (unsigned int point = 0; point < n_q_points; ++point)
-               {
-                       for (unsigned int i = 0; i < n_dofs; ++i)
-                       {
-                               for (unsigned int j = 0; j < n_dofs; ++j)
-                               {
-
-                                       copy_data.cell_mass_matrix(i, j) += fe_v.shape_value(i, point) * //phi_i(x_q)
-                                                                                                               fe_v.shape_value(j, point) * //phi_j(x_q)
-                                                                                                               JxW[point];                                      // dx(x_q)
-                               }
-                               copy_data.cell_mass_rhs(i) +=
-                                       (rhs.value(q_points[point]) *   // f(x_q)
-                                                fe_v.shape_value(i, point) //phi_i(x_q)
-                                        -
-                                        advection_coeff.value(q_points[point]) *
-                                                fe_v.shape_value(i, point) *             //c*phi_i(x_q)
-                                                sol_u_at_quadrature_points[point]) * //u_h(x_q)
-                                       JxW[point];                                                               //dx
-                       }
-               }
-
-               FullMatrix<double> inverse(fe_v.n_quadrature_points, fe_v.n_quadrature_points);
-               inverse.invert(copy_data.cell_mass_matrix);
-               Vector<double> proj(fe_v.n_quadrature_points); //projection of (f-c*U_h) on the local fe_space
-               inverse.vmult(proj, copy_data.cell_mass_rhs);  //M^{-1}*rhs = proj
-
-               double square_norm_over_cell = 0.0;
-               for (unsigned int point = 0; point < n_q_points; ++point)
-               {
-                       const double diff = rhs.value(q_points[point]) - sol_u_at_quadrature_points[point] - proj[point];
-                       square_norm_over_cell += diff * diff * JxW[point];
-               }
-               copy_data.value_estimator = square_norm_over_cell;
-       };
-
-       // Finally we have the boundary term with $||\beta (g-u_h^+)||^2$
-       const auto boundary_worker = [&](const Iterator &cell,
-                                                                        const unsigned int &face_no,
-                                                                        ScratchData<dim> &scratch_data,
-                                                                        CopyData &copy_data)
-       {
-               scratch_data.fe_interface_values.reinit(cell, face_no);
-               const FEFaceValuesBase<dim> &fe_fv = scratch_data.fe_interface_values.get_fe_face_values(0);
-               const auto &q_points = fe_fv.get_quadrature_points();
-               const unsigned n_q_points = q_points.size();
-               const std::vector<double> &JxW = fe_fv.get_JxW_values();
-
-               std::vector<double> g(n_q_points);
-               exact_solution.value_list(q_points, g);
-
-               std::vector<double> sol_u(n_q_points);
-               fe_fv.get_function_values(solution, sol_u);
-
-               const std::vector<Tensor<1, dim>> &normals = fe_fv.get_normal_vectors();
-
-               double square_norm_over_bdary_face = 0.;
-               for (unsigned int point = 0; point < q_points.size(); ++point)
-               {
-                       const double beta_dot_n = beta(q_points[point]) * normals[point];
-
-                       if (beta_dot_n < 0) //\partial_{-T} \cap \partial_{- \Omega}
-                       {
-                               const double diff = std::abs(beta_dot_n) * (g[point] - sol_u[point]);
-                               square_norm_over_bdary_face += diff * diff * JxW[point];
-                       }
-               }
-               copy_data.value_estimator += square_norm_over_bdary_face;
-       };
-
-       // Then compute the interior face terms with $|| \sqrt{b \cdot n}[u_h]||^2$
-       const auto face_worker = [&](const Iterator &cell,
-                                                                const unsigned int &f,
-                                                                const unsigned int &sf,
-                                                                const Iterator &ncell,
-                                                                const unsigned int &nf,
-                                                                const unsigned int &nsf,
-                                                                ScratchData<dim> &scratch_data,
-                                                                CopyData &copy_data)
-       {
-               FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
-               fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
-
-               copy_data.face_data.emplace_back();
-               CopyDataFace &copy_data_face = copy_data.face_data.back();
-               copy_data_face.cell_indices[0] = cell->active_cell_index();
-               copy_data_face.cell_indices[1] = ncell->active_cell_index();
-
-               const auto &q_points = fe_iv.get_quadrature_points();
-               const unsigned n_q_points = q_points.size();
-
-               const std::vector<double> &JxW = fe_iv.get_JxW_values();
-               std::vector<double> g(n_q_points);
-
-               std::vector<double> jump(n_q_points);
-               get_function_jump(fe_iv, solution, jump);
-
-               const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
-
-               double error_jump_square{0.0};
-               for (unsigned int point = 0; point < n_q_points; ++point)
-               {
-                       const double beta_dot_n = beta(q_points[point]) * normals[point];
-                       if (beta_dot_n < 0)
-                       {
-                               error_jump_square += std::abs(beta_dot_n) * jump[point] * jump[point] * JxW[point];
-                       }
-               }
-
-               copy_data_face.values[0] = error_jump_square;
-               copy_data_face.values[1] = copy_data_face.values[0];
-       };
-
-       ScratchData<dim> scratch_data(mapping, *fe, QGauss<dim>{fe->tensor_degree() + 1},
-                                                                 QGauss<dim - 1>{fe->tensor_degree() + 1});
-
-       const auto copier = [&](const auto &copy_data)
-       {
-               if (copy_data.cell_index != numbers::invalid_unsigned_int)
-               {
-                       error_indicator_per_cell[copy_data.cell_index] += copy_data.value_estimator;
-               }
-               for (auto &cdf : copy_data.face_data)
-               {
-                       for (unsigned int j = 0; j < 2; ++j)
-                       {
-                               error_indicator_per_cell[cdf.cell_indices[j]] += cdf.values[j];
-                       }
-               }
-       };
-
-       // Here, we finally handle the assembly of the Mass matrix (M)_{ij} = (\phi_j, \phi_i)_T. We pass in ScratchData and
-       // CopyData objects
-       CopyData copy_data;
-       MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(),
-                                                 cell_worker, copier, scratch_data, copy_data,
-                                                 MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces | MeshWorker::assemble_own_interior_faces_once,
-                                                 boundary_worker, face_worker);
+  // Compute the term $||f-c u_h - \Pi(f- cu_h)||_T^2$
+  using Iterator = typename DoFHandler<dim>::active_cell_iterator;
+  error_indicator_per_cell.reinit(triangulation.n_active_cells());
+
+  const auto cell_worker = [&](const Iterator   &cell,
+                               ScratchData<dim> &scratch_data,
+                               CopyData         &copy_data) {
+    const unsigned int n_dofs =
+      scratch_data.fe_values.get_fe().n_dofs_per_cell();
+
+    copy_data.reinit(cell, n_dofs);
+    scratch_data.fe_values.reinit(cell);
+    copy_data.cell_index = cell->active_cell_index();
+
+    const auto    &q_points   = scratch_data.fe_values.get_quadrature_points();
+    const unsigned n_q_points = q_points.size();
+
+    const FEValues<dim>       &fe_v = scratch_data.fe_values;
+    const std::vector<double> &JxW  = fe_v.get_JxW_values();
+
+    std::vector<double> sol_u_at_quadrature_points(fe_v.n_quadrature_points);
+    fe_v.get_function_values(solution, sol_u_at_quadrature_points);
+
+    // Compute local L^2 projection of  $f- c u_h$ over the local finite element
+    // space
+    for (unsigned int point = 0; point < n_q_points; ++point)
+      {
+        for (unsigned int i = 0; i < n_dofs; ++i)
+          {
+            for (unsigned int j = 0; j < n_dofs; ++j)
+              {
+                copy_data.cell_mass_matrix(i, j) +=
+                  fe_v.shape_value(i, point) * // phi_i(x_q)
+                  fe_v.shape_value(j, point) * // phi_j(x_q)
+                  JxW[point];                  // dx(x_q)
+              }
+            copy_data.cell_mass_rhs(i) +=
+              (rhs.value(q_points[point]) * // f(x_q)
+                 fe_v.shape_value(i, point) // phi_i(x_q)
+               - advection_coeff.value(q_points[point]) *
+                   fe_v.shape_value(i, point) *         // c*phi_i(x_q)
+                   sol_u_at_quadrature_points[point]) * // u_h(x_q)
+              JxW[point];                               // dx
+          }
+      }
+
+    FullMatrix<double> inverse(fe_v.n_quadrature_points,
+                               fe_v.n_quadrature_points);
+    inverse.invert(copy_data.cell_mass_matrix);
+    Vector<double> proj(fe_v.n_quadrature_points); // projection of (f-c*U_h) on
+                                                   // the local fe_space
+    inverse.vmult(proj, copy_data.cell_mass_rhs);  // M^{-1}*rhs = proj
+
+    double square_norm_over_cell = 0.0;
+    for (unsigned int point = 0; point < n_q_points; ++point)
+      {
+        const double diff = rhs.value(q_points[point]) -
+                            sol_u_at_quadrature_points[point] - proj[point];
+        square_norm_over_cell += diff * diff * JxW[point];
+      }
+    copy_data.value_estimator = square_norm_over_cell;
+  };
+
+  // Finally we have the boundary term with $||\beta (g-u_h^+)||^2$
+  const auto boundary_worker = [&](const Iterator     &cell,
+                                   const unsigned int &face_no,
+                                   ScratchData<dim>   &scratch_data,
+                                   CopyData           &copy_data) {
+    scratch_data.fe_interface_values.reinit(cell, face_no);
+    const FEFaceValuesBase<dim> &fe_fv =
+      scratch_data.fe_interface_values.get_fe_face_values(0);
+    const auto                &q_points   = fe_fv.get_quadrature_points();
+    const unsigned             n_q_points = q_points.size();
+    const std::vector<double> &JxW        = fe_fv.get_JxW_values();
+
+    std::vector<double> g(n_q_points);
+    exact_solution.value_list(q_points, g);
+
+    std::vector<double> sol_u(n_q_points);
+    fe_fv.get_function_values(solution, sol_u);
+
+    const std::vector<Tensor<1, dim>> &normals = fe_fv.get_normal_vectors();
+
+    double square_norm_over_bdary_face = 0.;
+    for (unsigned int point = 0; point < q_points.size(); ++point)
+      {
+        const double beta_dot_n = beta(q_points[point]) * normals[point];
+
+        if (beta_dot_n < 0) //\partial_{-T} \cap \partial_{- \Omega}
+          {
+            const double diff =
+              std::abs(beta_dot_n) * (g[point] - sol_u[point]);
+            square_norm_over_bdary_face += diff * diff * JxW[point];
+          }
+      }
+    copy_data.value_estimator += square_norm_over_bdary_face;
+  };
+
+  // Then compute the interior face terms with $|| \sqrt{b \cdot n}[u_h]||^2$
+  const auto face_worker = [&](const Iterator     &cell,
+                               const unsigned int &f,
+                               const unsigned int &sf,
+                               const Iterator     &ncell,
+                               const unsigned int &nf,
+                               const unsigned int &nsf,
+                               ScratchData<dim>   &scratch_data,
+                               CopyData           &copy_data) {
+    FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
+    fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
+
+    copy_data.face_data.emplace_back();
+    CopyDataFace &copy_data_face   = copy_data.face_data.back();
+    copy_data_face.cell_indices[0] = cell->active_cell_index();
+    copy_data_face.cell_indices[1] = ncell->active_cell_index();
+
+    const auto    &q_points   = fe_iv.get_quadrature_points();
+    const unsigned n_q_points = q_points.size();
+
+    const std::vector<double> &JxW = fe_iv.get_JxW_values();
+    std::vector<double>        g(n_q_points);
+
+    std::vector<double> jump(n_q_points);
+    get_function_jump(fe_iv, solution, jump);
+
+    const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
+
+    double error_jump_square{0.0};
+    for (unsigned int point = 0; point < n_q_points; ++point)
+      {
+        const double beta_dot_n = beta(q_points[point]) * normals[point];
+        if (beta_dot_n < 0)
+          {
+            error_jump_square +=
+              std::abs(beta_dot_n) * jump[point] * jump[point] * JxW[point];
+          }
+      }
+
+    copy_data_face.values[0] = error_jump_square;
+    copy_data_face.values[1] = copy_data_face.values[0];
+  };
+
+  ScratchData<dim> scratch_data(mapping,
+                                *fe,
+                                QGauss<dim>{fe->tensor_degree() + 1},
+                                QGauss<dim - 1>{fe->tensor_degree() + 1});
+
+  const auto copier = [&](const auto &copy_data) {
+    if (copy_data.cell_index != numbers::invalid_unsigned_int)
+      {
+        error_indicator_per_cell[copy_data.cell_index] +=
+          copy_data.value_estimator;
+      }
+    for (auto &cdf : copy_data.face_data)
+      {
+        for (unsigned int j = 0; j < 2; ++j)
+          {
+            error_indicator_per_cell[cdf.cell_indices[j]] += cdf.values[j];
+          }
+      }
+  };
+
+  // Here, we finally handle the assembly of the Mass matrix (M)_{ij} = (\phi_j,
+  // \phi_i)_T. We pass in ScratchData and CopyData objects
+  CopyData copy_data;
+  MeshWorker::mesh_loop(dof_handler.begin_active(),
+                        dof_handler.end(),
+                        cell_worker,
+                        copier,
+                        scratch_data,
+                        copy_data,
+                        MeshWorker::assemble_own_cells |
+                          MeshWorker::assemble_boundary_faces |
+                          MeshWorker::assemble_own_interior_faces_once,
+                        boundary_worker,
+                        face_worker);
 }
 
-//Usual <code>run</code> function, which runs over several refinement cycles
+
+
+// Usual <code>run</code> function, which runs over several refinement cycles
 template <int dim>
-void AdvectionReaction<dim>::run()
+void
+AdvectionReaction<dim>::run()
 {
-       std::vector<double> energy_errors;
-       std::vector<int> dofs_hist;
-       for (unsigned int cycle = 0; cycle < n_refinement_cycles; ++cycle)
-       {
-               std::cout << "Cycle " << cycle << std::endl;
-
-               if (cycle == 0)
-               {
-                       GridGenerator::hyper_cube(triangulation);
-                       triangulation.refine_global(n_global_refinements);
-               }
-               else
-               {
-                       refine_grid();
-               }
-               std::cout << "  Number of active cells:       "
-                                 << triangulation.n_active_cells() << std::endl;
-
-               setup_system();
-
-               std::cout << "  Number of degrees of freedom: " << dof_handler.n_dofs()
-                                 << std::endl;
-
-               assemble_system();
-               solve();
-               compute_error();
-               output_results(cycle);
-               energy_errors.emplace_back(compute_energy_norm());
-               dofs_hist.emplace_back(triangulation.n_active_cells());
-       }
-       error_table.output_table(std::cout);
-
-       for (unsigned int i = 0; i < n_refinement_cycles; ++i)
-               std::cout << "Cycle " << i << "\t" << energy_errors[i] << '\n';
-       {
-       }
+  std::vector<double> energy_errors;
+  std::vector<int>    dofs_hist;
+  for (unsigned int cycle = 0; cycle < n_refinement_cycles; ++cycle)
+    {
+      std::cout << "Cycle " << cycle << std::endl;
+
+      if (cycle == 0)
+        {
+          GridGenerator::hyper_cube(triangulation);
+          triangulation.refine_global(n_global_refinements);
+        }
+      else
+        {
+          refine_grid();
+        }
+      std::cout << "  Number of active cells:       "
+                << triangulation.n_active_cells() << std::endl;
+      std::cout << "  Number of degrees of freedom: " << dof_handler.n_dofs()
+                << std::endl;
+
+      setup_system();
+      assemble_system();
+      solve();
+      compute_error();
+      output_results(cycle);
+
+      energy_errors.emplace_back(compute_energy_norm());
+      dofs_hist.emplace_back(triangulation.n_active_cells());
+    }
+  error_table.output_table(std::cout);
+
+  for (unsigned int i = 0; i < n_refinement_cycles; ++i)
+    std::cout << "Cycle " << i << "\t" << energy_errors[i] << std::endl;
 }
 // Explicit instantiation
+template class AdvectionReaction<1>;
 template class AdvectionReaction<2>;
+template class AdvectionReaction<3>;

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.