From: David Wells <drwells@email.unc.edu>
Date: Mon, 25 Jul 2022 15:28:18 +0000 (-0400)
Subject: Run clang-format on everything
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F116%2Fhead;p=code-gallery.git

Run clang-format on everything
---

diff --git a/cdr/common/include/deal.II-cdr/parameters.h b/cdr/common/include/deal.II-cdr/parameters.h
index 83e21bf..4a88d3b 100644
--- a/cdr/common/include/deal.II-cdr/parameters.h
+++ b/cdr/common/include/deal.II-cdr/parameters.h
@@ -23,22 +23,25 @@ namespace CDR
 
     double diffusion_coefficient;
     double reaction_coefficient;
-    bool time_dependent_forcing;
+    bool   time_dependent_forcing;
 
     unsigned int initial_refinement_level;
     unsigned int max_refinement_level;
     unsigned int fe_order;
 
-    double start_time;
-    double stop_time;
+    double       start_time;
+    double       stop_time;
     unsigned int n_time_steps;
 
     unsigned int save_interval;
     unsigned int patch_level;
 
-    void read_parameter_file(const std::string &file_name);
+    void
+    read_parameter_file(const std::string &file_name);
+
   private:
-    void configure_parameter_handler(ParameterHandler &parameter_handler);
+    void
+    configure_parameter_handler(ParameterHandler &parameter_handler);
   };
-}
+} // namespace CDR
 #endif
diff --git a/cdr/common/include/deal.II-cdr/system_matrix.h b/cdr/common/include/deal.II-cdr/system_matrix.h
index 41fd50b..489497f 100644
--- a/cdr/common/include/deal.II-cdr/system_matrix.h
+++ b/cdr/common/include/deal.II-cdr/system_matrix.h
@@ -18,23 +18,25 @@ namespace CDR
 {
   using namespace dealii;
 
-  template<int dim, typename MatrixType>
-  void create_system_matrix
-  (const DoFHandler<dim>                                 &dof_handler,
-   const QGauss<dim>                                     &quad,
-   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
-   const CDR::Parameters                                 &parameters,
-   const double                                           time_step,
-   MatrixType                                            &system_matrix);
+  template <int dim, typename MatrixType>
+  void
+  create_system_matrix(
+    const DoFHandler<dim> &                                dof_handler,
+    const QGauss<dim> &                                    quad,
+    const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+    const CDR::Parameters &                                parameters,
+    const double                                           time_step,
+    MatrixType &                                           system_matrix);
 
-  template<int dim, typename MatrixType>
-  void create_system_matrix
-  (const DoFHandler<dim>                                 &dof_handler,
-   const QGauss<dim>                                     &quad,
-   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
-   const CDR::Parameters                                 &parameters,
-   const double                                           time_step,
-   const AffineConstraints<double>                       &constraints,
-   MatrixType                                            &system_matrix);
-}
+  template <int dim, typename MatrixType>
+  void
+  create_system_matrix(
+    const DoFHandler<dim> &                                dof_handler,
+    const QGauss<dim> &                                    quad,
+    const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+    const CDR::Parameters &                                parameters,
+    const double                                           time_step,
+    const AffineConstraints<double> &                      constraints,
+    MatrixType &                                           system_matrix);
+} // namespace CDR
 #endif
diff --git a/cdr/common/include/deal.II-cdr/system_matrix.templates.h b/cdr/common/include/deal.II-cdr/system_matrix.templates.h
index 12340b5..8e4cabe 100644
--- a/cdr/common/include/deal.II-cdr/system_matrix.templates.h
+++ b/cdr/common/include/deal.II-cdr/system_matrix.templates.h
@@ -4,11 +4,11 @@
 
 #include <deal.II/dofs/dof_handler.h>
 
-#include <deal.II/lac/affine_constraints.h>
-
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_values.h>
 
+#include <deal.II/lac/affine_constraints.h>
+
 #include <deal.II-cdr/parameters.h>
 #include <deal.II-cdr/system_matrix.h>
 
@@ -22,20 +22,23 @@ namespace CDR
   // This is the actual implementation of the <code>create_system_matrix</code>
   // function described in the header file. It is similar to the system matrix
   // assembly routine in step-40.
-  template<int dim, typename UpdateFunction>
-  void internal_create_system_matrix
-  (const DoFHandler<dim>                                 &dof_handler,
-   const QGauss<dim>                                     &quad,
-   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
-   const CDR::Parameters                                 &parameters,
-   const double                                           time_step,
-   UpdateFunction                                         update_system_matrix)
+  template <int dim, typename UpdateFunction>
+  void
+  internal_create_system_matrix(
+    const DoFHandler<dim> &                                dof_handler,
+    const QGauss<dim> &                                    quad,
+    const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+    const CDR::Parameters &                                parameters,
+    const double                                           time_step,
+    UpdateFunction                                         update_system_matrix)
   {
-    auto &fe = dof_handler.get_fe();
-    const auto dofs_per_cell = fe.dofs_per_cell;
+    auto &             fe            = dof_handler.get_fe();
+    const auto         dofs_per_cell = fe.dofs_per_cell;
     FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
-    FEValues<dim> fe_values(fe, quad, update_values | update_gradients |
-                            update_quadrature_points | update_JxW_values);
+    FEValues<dim>      fe_values(fe,
+                            quad,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
 
     std::vector<types::global_dof_index> local_indices(dofs_per_cell);
 
@@ -55,19 +58,23 @@ namespace CDR
                   {
                     for (unsigned int j = 0; j < dofs_per_cell; ++j)
                       {
-                        const auto convection_contribution = current_convection
-                                                             *fe_values.shape_grad(j, q);
-                        cell_matrix(i, j) += fe_values.JxW(q)*
-                                             // Here are the time step, mass, and reaction parts:
-                                             ((1.0 + time_step/2.0*parameters.reaction_coefficient)
-                                              *fe_values.shape_value(i, q)*fe_values.shape_value(j, q)
-                                              + time_step/2.0*
-                                              // and the convection part:
-                                              (fe_values.shape_value(i, q)*convection_contribution
-                                               // and, finally, the diffusion part:
-                                               + parameters.diffusion_coefficient
-                                               *(fe_values.shape_grad(i, q)*fe_values.shape_grad(j, q)))
-                                             );
+                        const auto convection_contribution =
+                          current_convection * fe_values.shape_grad(j, q);
+                        cell_matrix(i, j) +=
+                          fe_values.JxW(q) *
+                          // Here are the time step, mass, and reaction parts:
+                          ((1.0 +
+                            time_step / 2.0 * parameters.reaction_coefficient) *
+                             fe_values.shape_value(i, q) *
+                             fe_values.shape_value(j, q) +
+                           time_step / 2.0 *
+                             // and the convection part:
+                             (fe_values.shape_value(i, q) *
+                                convection_contribution
+                              // and, finally, the diffusion part:
+                              + parameters.diffusion_coefficient *
+                                  (fe_values.shape_grad(i, q) *
+                                   fe_values.shape_grad(j, q))));
                       }
                   }
               }
@@ -76,42 +83,53 @@ namespace CDR
       }
   }
 
-  template<int dim, typename MatrixType>
-  void create_system_matrix
-  (const DoFHandler<dim>                                 &dof_handler,
-   const QGauss<dim>                                     &quad,
-   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
-   const CDR::Parameters                                 &parameters,
-   const double                                           time_step,
-   const AffineConstraints<double>                       &constraints,
-   MatrixType                                            &system_matrix)
+  template <int dim, typename MatrixType>
+  void
+  create_system_matrix(
+    const DoFHandler<dim> &                                dof_handler,
+    const QGauss<dim> &                                    quad,
+    const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+    const CDR::Parameters &                                parameters,
+    const double                                           time_step,
+    const AffineConstraints<double> &                      constraints,
+    MatrixType &                                           system_matrix)
   {
-    internal_create_system_matrix<dim>
-    (dof_handler, quad, convection_function, parameters, time_step,
-     [&constraints, &system_matrix](const std::vector<types::global_dof_index> &local_indices,
-                                    const FullMatrix<double> &cell_matrix)
-    {
-      constraints.distribute_local_to_global
-      (cell_matrix, local_indices, system_matrix);
-    });
+    internal_create_system_matrix<dim>(
+      dof_handler,
+      quad,
+      convection_function,
+      parameters,
+      time_step,
+      [&constraints, &system_matrix](
+        const std::vector<types::global_dof_index> &local_indices,
+        const FullMatrix<double> &                  cell_matrix) {
+        constraints.distribute_local_to_global(cell_matrix,
+                                               local_indices,
+                                               system_matrix);
+      });
   }
 
-  template<int dim, typename MatrixType>
-  void create_system_matrix
-  (const DoFHandler<dim>                                 &dof_handler,
-   const QGauss<dim>                                     &quad,
-   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
-   const CDR::Parameters                                 &parameters,
-   const double                                           time_step,
-   MatrixType                                            &system_matrix)
+  template <int dim, typename MatrixType>
+  void
+  create_system_matrix(
+    const DoFHandler<dim> &                                dof_handler,
+    const QGauss<dim> &                                    quad,
+    const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+    const CDR::Parameters &                                parameters,
+    const double                                           time_step,
+    MatrixType &                                           system_matrix)
   {
-    internal_create_system_matrix<dim>
-    (dof_handler, quad, convection_function, parameters, time_step,
-     [&system_matrix](const std::vector<types::global_dof_index> &local_indices,
-                      const FullMatrix<double> &cell_matrix)
-    {
-      system_matrix.add(local_indices, cell_matrix);
-    });
+    internal_create_system_matrix<dim>(
+      dof_handler,
+      quad,
+      convection_function,
+      parameters,
+      time_step,
+      [&system_matrix](
+        const std::vector<types::global_dof_index> &local_indices,
+        const FullMatrix<double> &                  cell_matrix) {
+        system_matrix.add(local_indices, cell_matrix);
+      });
   }
-}
+} // namespace CDR
 #endif
diff --git a/cdr/common/include/deal.II-cdr/system_rhs.h b/cdr/common/include/deal.II-cdr/system_rhs.h
index 986b5a1..a21d043 100644
--- a/cdr/common/include/deal.II-cdr/system_rhs.h
+++ b/cdr/common/include/deal.II-cdr/system_rhs.h
@@ -18,16 +18,17 @@ namespace CDR
 {
   using namespace dealii;
 
-  template<int dim, typename VectorType>
-  void create_system_rhs
-  (const DoFHandler<dim>                                 &dof_handler,
-   const QGauss<dim>                                     &quad,
-   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
-   const std::function<double(double, const Point<dim>)> &forcing_function,
-   const CDR::Parameters                                 &parameters,
-   const VectorType                                      &previous_solution,
-   const AffineConstraints<double>                       &constraints,
-   const double                                           current_time,
-   VectorType                                            &system_rhs);
-}
+  template <int dim, typename VectorType>
+  void
+  create_system_rhs(
+    const DoFHandler<dim> &                                dof_handler,
+    const QGauss<dim> &                                    quad,
+    const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+    const std::function<double(double, const Point<dim>)> &forcing_function,
+    const CDR::Parameters &                                parameters,
+    const VectorType &                                     previous_solution,
+    const AffineConstraints<double> &                      constraints,
+    const double                                           current_time,
+    VectorType &                                           system_rhs);
+} // namespace CDR
 #endif
diff --git a/cdr/common/include/deal.II-cdr/system_rhs.templates.h b/cdr/common/include/deal.II-cdr/system_rhs.templates.h
index e7e2ae6..34f1378 100644
--- a/cdr/common/include/deal.II-cdr/system_rhs.templates.h
+++ b/cdr/common/include/deal.II-cdr/system_rhs.templates.h
@@ -22,32 +22,35 @@ namespace CDR
 {
   using namespace dealii;
 
-  template<int dim, typename VectorType>
-  void create_system_rhs
-  (const DoFHandler<dim>                                 &dof_handler,
-   const QGauss<dim>                                     &quad,
-   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
-   const std::function<double(double, const Point<dim>)> &forcing_function,
-   const CDR::Parameters                                 &parameters,
-   const VectorType                                      &previous_solution,
-   const AffineConstraints<double>                       &constraints,
-   const double                                           current_time,
-   VectorType                                            &system_rhs)
+  template <int dim, typename VectorType>
+  void
+  create_system_rhs(
+    const DoFHandler<dim> &                                dof_handler,
+    const QGauss<dim> &                                    quad,
+    const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
+    const std::function<double(double, const Point<dim>)> &forcing_function,
+    const CDR::Parameters &                                parameters,
+    const VectorType &                                     previous_solution,
+    const AffineConstraints<double> &                      constraints,
+    const double                                           current_time,
+    VectorType &                                           system_rhs)
   {
-    auto &fe = dof_handler.get_fe();
-    const auto dofs_per_cell = fe.dofs_per_cell;
-    const double time_step = (parameters.stop_time - parameters.start_time)
-                             /parameters.n_time_steps;
-    FEValues<dim> fe_values(fe, quad, update_values | update_gradients |
-                            update_quadrature_points | update_JxW_values);
+    auto &       fe            = dof_handler.get_fe();
+    const auto   dofs_per_cell = fe.dofs_per_cell;
+    const double time_step =
+      (parameters.stop_time - parameters.start_time) / parameters.n_time_steps;
+    FEValues<dim> fe_values(fe,
+                            quad,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
 
-    Vector<double> cell_rhs(dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
     FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 
-    Vector<double> current_fe_coefficients(dofs_per_cell);
+    Vector<double>                       current_fe_coefficients(dofs_per_cell);
     std::vector<types::global_dof_index> local_indices(dofs_per_cell);
 
-    const double previous_time {current_time - time_step};
+    const double previous_time{current_time - time_step};
 
     for (const auto &cell : dof_handler.active_cell_iterators())
       {
@@ -58,7 +61,8 @@ namespace CDR
             cell->get_dof_indices(local_indices);
             for (unsigned int i = 0; i < dofs_per_cell; ++i)
               {
-                current_fe_coefficients[i] = previous_solution[local_indices[i]];
+                current_fe_coefficients[i] =
+                  previous_solution[local_indices[i]];
               }
 
             for (unsigned int q = 0; q < quad.size(); ++q)
@@ -66,38 +70,46 @@ namespace CDR
                 const auto current_convection =
                   convection_function(fe_values.quadrature_point(q));
 
-                const double current_forcing = forcing_function
-                                               (current_time, fe_values.quadrature_point(q));
-                const double previous_forcing = forcing_function
-                                                (previous_time, fe_values.quadrature_point(q));
+                const double current_forcing =
+                  forcing_function(current_time, fe_values.quadrature_point(q));
+                const double previous_forcing =
+                  forcing_function(previous_time,
+                                   fe_values.quadrature_point(q));
                 for (unsigned int i = 0; i < dofs_per_cell; ++i)
                   {
                     for (unsigned int j = 0; j < dofs_per_cell; ++j)
                       {
-                        const auto convection_contribution = current_convection
-                                                             *fe_values.shape_grad(j, q);
+                        const auto convection_contribution =
+                          current_convection * fe_values.shape_grad(j, q);
 
-                        cell_rhs(i) += fe_values.JxW(q)*
-                                       // Here are the mass and reaction part:
-                                       (((1.0 - time_step/2.0*parameters.reaction_coefficient)
-                                         *fe_values.shape_value(i, q)*fe_values.shape_value(j, q)
-                                         - time_step/2.0*
-                                         // the convection part:
-                                         (fe_values.shape_value(i, q)*convection_contribution
-                                          // the diffusion part:
-                                          + parameters.diffusion_coefficient
-                                          *(fe_values.shape_grad(i, q)*fe_values.shape_grad(j, q))))
-                                        *current_fe_coefficients[j]
-                                        // and, finally, the forcing function part:
-                                        + time_step/2.0*
-                                        (current_forcing + previous_forcing)
-                                        *fe_values.shape_value(i, q));
+                        cell_rhs(i) +=
+                          fe_values.JxW(q) *
+                          // Here are the mass and reaction part:
+                          (((1.0 - time_step / 2.0 *
+                                     parameters.reaction_coefficient) *
+                              fe_values.shape_value(i, q) *
+                              fe_values.shape_value(j, q) -
+                            time_step / 2.0 *
+                              // the convection part:
+                              (fe_values.shape_value(i, q) *
+                                 convection_contribution
+                               // the diffusion part:
+                               + parameters.diffusion_coefficient *
+                                   (fe_values.shape_grad(i, q) *
+                                    fe_values.shape_grad(j, q)))) *
+                             current_fe_coefficients[j]
+                           // and, finally, the forcing function part:
+                           + time_step / 2.0 *
+                               (current_forcing + previous_forcing) *
+                               fe_values.shape_value(i, q));
                       }
                   }
               }
-            constraints.distribute_local_to_global(cell_rhs, local_indices, system_rhs);
+            constraints.distribute_local_to_global(cell_rhs,
+                                                   local_indices,
+                                                   system_rhs);
           }
       }
   }
-}
+} // namespace CDR
 #endif
diff --git a/cdr/common/include/deal.II-cdr/write_pvtu_output.h b/cdr/common/include/deal.II-cdr/write_pvtu_output.h
index 5b0b1c2..9530d88 100644
--- a/cdr/common/include/deal.II-cdr/write_pvtu_output.h
+++ b/cdr/common/include/deal.II-cdr/write_pvtu_output.h
@@ -12,14 +12,16 @@ namespace CDR
   public:
     WritePVTUOutput(const unsigned int patch_level);
 
-    template<int dim, typename VectorType>
-    void write_output(const DoFHandler<dim> &dof_handler,
-                      const VectorType      &solution,
-                      const unsigned int     time_step_n,
-                      const double           current_time);
+    template <int dim, typename VectorType>
+    void
+    write_output(const DoFHandler<dim> &dof_handler,
+                 const VectorType &     solution,
+                 const unsigned int     time_step_n,
+                 const double           current_time);
+
   private:
     const unsigned int patch_level;
     const unsigned int this_mpi_process;
   };
-}
+} // namespace CDR
 #endif
diff --git a/cdr/common/include/deal.II-cdr/write_pvtu_output.templates.h b/cdr/common/include/deal.II-cdr/write_pvtu_output.templates.h
index 15e5154..a0af9ba 100644
--- a/cdr/common/include/deal.II-cdr/write_pvtu_output.templates.h
+++ b/cdr/common/include/deal.II-cdr/write_pvtu_output.templates.h
@@ -11,8 +11,8 @@
 
 #include <deal.II-cdr/write_pvtu_output.h>
 
-#include <string>
 #include <fstream>
+#include <string>
 #include <vector>
 
 // Here is the implementation of the important function. This is similar to
@@ -21,18 +21,19 @@ namespace CDR
 {
   using namespace dealii;
 
-  template<int dim, typename VectorType>
-  void WritePVTUOutput::write_output(const DoFHandler<dim> &dof_handler,
-                                     const VectorType      &solution,
-                                     const unsigned int     time_step_n,
-                                     const double           current_time)
+  template <int dim, typename VectorType>
+  void
+  WritePVTUOutput::write_output(const DoFHandler<dim> &dof_handler,
+                                const VectorType &     solution,
+                                const unsigned int     time_step_n,
+                                const double           current_time)
   {
     DataOut<dim> data_out;
     data_out.attach_dof_handler(dof_handler);
     data_out.add_data_vector(solution, "u");
 
-    const auto &triangulation = dof_handler.get_triangulation();
-    Vector<float> subdomain (triangulation.n_active_cells());
+    const auto &  triangulation = dof_handler.get_triangulation();
+    Vector<float> subdomain(triangulation.n_active_cells());
     for (auto &domain : subdomain)
       {
         domain = triangulation.locally_owned_subdomain();
@@ -44,7 +45,8 @@ namespace CDR
     flags.time = current_time;
     // While the default flag is for the best compression level, using
     // <code>best_speed</code> makes this function much faster.
-    flags.compression_level = DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed;
+    flags.compression_level =
+      DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed;
     data_out.set_flags(flags);
 
     unsigned int subdomain_n;
@@ -57,25 +59,25 @@ namespace CDR
         subdomain_n = triangulation.locally_owned_subdomain();
       }
 
-    std::ofstream output
-    ("solution-" + Utilities::int_to_string(time_step_n) + "."
-     + Utilities::int_to_string(subdomain_n, 4)
-     + ".vtu");
+    std::ofstream output("solution-" + Utilities::int_to_string(time_step_n) +
+                         "." + Utilities::int_to_string(subdomain_n, 4) +
+                         ".vtu");
 
     data_out.write_vtu(output);
 
     if (this_mpi_process == 0)
       {
         std::vector<std::string> filenames;
-        for (unsigned int i = 0; i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
              ++i)
-          filenames.push_back
-          ("solution-" + Utilities::int_to_string (time_step_n) + "."
-           + Utilities::int_to_string (i, 4) + ".vtu");
-        std::ofstream master_output
-        ("solution-" + Utilities::int_to_string(time_step_n) + ".pvtu");
+          filenames.push_back("solution-" +
+                              Utilities::int_to_string(time_step_n) + "." +
+                              Utilities::int_to_string(i, 4) + ".vtu");
+        std::ofstream master_output(
+          "solution-" + Utilities::int_to_string(time_step_n) + ".pvtu");
         data_out.write_pvtu_record(master_output, filenames);
       }
   }
-}
+} // namespace CDR
 #endif
diff --git a/cdr/common/source/parameters.cc b/cdr/common/source/parameters.cc
index daefef0..16be678 100644
--- a/cdr/common/source/parameters.cc
+++ b/cdr/common/source/parameters.cc
@@ -5,64 +5,90 @@
 
 namespace CDR
 {
-  void Parameters::configure_parameter_handler(ParameterHandler &parameter_handler)
+  void
+  Parameters::configure_parameter_handler(ParameterHandler &parameter_handler)
   {
     parameter_handler.enter_subsection("Geometry");
     {
-      parameter_handler.declare_entry
-      ("inner_radius", "1.0", Patterns::Double(0.0), "Inner radius.");
-      parameter_handler.declare_entry
-      ("outer_radius", "2.0", Patterns::Double(0.0), "Outer radius.");
+      parameter_handler.declare_entry("inner_radius",
+                                      "1.0",
+                                      Patterns::Double(0.0),
+                                      "Inner radius.");
+      parameter_handler.declare_entry("outer_radius",
+                                      "2.0",
+                                      Patterns::Double(0.0),
+                                      "Outer radius.");
     }
     parameter_handler.leave_subsection();
 
     parameter_handler.enter_subsection("Physical Parameters");
     {
-      parameter_handler.declare_entry
-      ("diffusion_coefficient", "1.0", Patterns::Double(0.0), "Diffusion coefficient.");
-      parameter_handler.declare_entry
-      ("reaction_coefficient", "1.0", Patterns::Double(0.0), "Reaction coefficient.");
-      parameter_handler.declare_entry
-      ("time_dependent_forcing", "true", Patterns::Bool(), "Whether or not "
-       "the forcing function depends on time.");
+      parameter_handler.declare_entry("diffusion_coefficient",
+                                      "1.0",
+                                      Patterns::Double(0.0),
+                                      "Diffusion coefficient.");
+      parameter_handler.declare_entry("reaction_coefficient",
+                                      "1.0",
+                                      Patterns::Double(0.0),
+                                      "Reaction coefficient.");
+      parameter_handler.declare_entry("time_dependent_forcing",
+                                      "true",
+                                      Patterns::Bool(),
+                                      "Whether or not "
+                                      "the forcing function depends on time.");
     }
     parameter_handler.leave_subsection();
 
     parameter_handler.enter_subsection("Finite Element");
     {
-      parameter_handler.declare_entry
-      ("initial_refinement_level", "1", Patterns::Integer(1),
-       "Initial number of levels in the mesh.");
-      parameter_handler.declare_entry
-      ("max_refinement_level", "1", Patterns::Integer(1),
-       "Maximum number of levels in the mesh.");
-      parameter_handler.declare_entry
-      ("fe_order", "1", Patterns::Integer(1), "Finite element order.");
+      parameter_handler.declare_entry("initial_refinement_level",
+                                      "1",
+                                      Patterns::Integer(1),
+                                      "Initial number of levels in the mesh.");
+      parameter_handler.declare_entry("max_refinement_level",
+                                      "1",
+                                      Patterns::Integer(1),
+                                      "Maximum number of levels in the mesh.");
+      parameter_handler.declare_entry("fe_order",
+                                      "1",
+                                      Patterns::Integer(1),
+                                      "Finite element order.");
     }
     parameter_handler.leave_subsection();
 
     parameter_handler.enter_subsection("Time Step");
     {
-      parameter_handler.declare_entry
-      ("start_time", "0.0", Patterns::Double(0.0), "Start time.");
-      parameter_handler.declare_entry
-      ("stop_time", "1.0", Patterns::Double(1.0), "Stop time.");
-      parameter_handler.declare_entry
-      ("n_time_steps", "1", Patterns::Integer(1), "Number of time steps.");
+      parameter_handler.declare_entry("start_time",
+                                      "0.0",
+                                      Patterns::Double(0.0),
+                                      "Start time.");
+      parameter_handler.declare_entry("stop_time",
+                                      "1.0",
+                                      Patterns::Double(1.0),
+                                      "Stop time.");
+      parameter_handler.declare_entry("n_time_steps",
+                                      "1",
+                                      Patterns::Integer(1),
+                                      "Number of time steps.");
     }
     parameter_handler.leave_subsection();
 
     parameter_handler.enter_subsection("Output");
     {
-      parameter_handler.declare_entry
-      ("save_interval", "10", Patterns::Integer(1), "Save interval.");
-      parameter_handler.declare_entry
-      ("patch_level", "2", Patterns::Integer(0), "Patch level.");
+      parameter_handler.declare_entry("save_interval",
+                                      "10",
+                                      Patterns::Integer(1),
+                                      "Save interval.");
+      parameter_handler.declare_entry("patch_level",
+                                      "2",
+                                      Patterns::Integer(0),
+                                      "Patch level.");
     }
     parameter_handler.leave_subsection();
   }
 
-  void Parameters::read_parameter_file(const std::string &file_name)
+  void
+  Parameters::read_parameter_file(const std::string &file_name)
   {
     ParameterHandler parameter_handler;
     {
@@ -80,25 +106,30 @@ namespace CDR
 
     parameter_handler.enter_subsection("Physical Parameters");
     {
-      diffusion_coefficient = parameter_handler.get_double("diffusion_coefficient");
-      reaction_coefficient = parameter_handler.get_double("reaction_coefficient");
-      time_dependent_forcing = parameter_handler.get_bool("time_dependent_forcing");
+      diffusion_coefficient =
+        parameter_handler.get_double("diffusion_coefficient");
+      reaction_coefficient =
+        parameter_handler.get_double("reaction_coefficient");
+      time_dependent_forcing =
+        parameter_handler.get_bool("time_dependent_forcing");
     }
     parameter_handler.leave_subsection();
 
 
     parameter_handler.enter_subsection("Finite Element");
     {
-      initial_refinement_level = parameter_handler.get_integer("initial_refinement_level");
-      max_refinement_level = parameter_handler.get_integer("max_refinement_level");
+      initial_refinement_level =
+        parameter_handler.get_integer("initial_refinement_level");
+      max_refinement_level =
+        parameter_handler.get_integer("max_refinement_level");
       fe_order = parameter_handler.get_integer("fe_order");
     }
     parameter_handler.leave_subsection();
 
     parameter_handler.enter_subsection("Time Step");
     {
-      start_time = parameter_handler.get_double("start_time");
-      stop_time = parameter_handler.get_double("stop_time");
+      start_time   = parameter_handler.get_double("start_time");
+      stop_time    = parameter_handler.get_double("stop_time");
       n_time_steps = parameter_handler.get_integer("n_time_steps");
     }
     parameter_handler.leave_subsection();
@@ -106,8 +137,8 @@ namespace CDR
     parameter_handler.enter_subsection("Output");
     {
       save_interval = parameter_handler.get_integer("save_interval");
-      patch_level = parameter_handler.get_integer("patch_level");
+      patch_level   = parameter_handler.get_integer("patch_level");
     }
     parameter_handler.leave_subsection();
   }
-}
+} // namespace CDR
diff --git a/cdr/common/source/system_matrix.cc b/cdr/common/source/system_matrix.cc
index 4e18902..bc03457 100644
--- a/cdr/common/source/system_matrix.cc
+++ b/cdr/common/source/system_matrix.cc
@@ -13,79 +13,79 @@ namespace CDR
 {
   using namespace dealii;
 
-  template
-  void create_system_matrix<2, SparseMatrix<double>>
-  (const DoFHandler<2>                               &dof_handler,
-   const QGauss<2>                                   &quad,
-   const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
-   const CDR::Parameters                             &parameters,
-   const double                                       time_step,
-   SparseMatrix<double>                              &system_matrix);
+  template void
+  create_system_matrix<2, SparseMatrix<double>>(
+    const DoFHandler<2> &                              dof_handler,
+    const QGauss<2> &                                  quad,
+    const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
+    const CDR::Parameters &                            parameters,
+    const double                                       time_step,
+    SparseMatrix<double> &                             system_matrix);
 
-  template
-  void create_system_matrix<3, SparseMatrix<double>>
-  (const DoFHandler<3>                               &dof_handler,
-   const QGauss<3>                                   &quad,
-   const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
-   const CDR::Parameters                             &parameters,
-   const double                                       time_step,
-   SparseMatrix<double>                              &system_matrix);
+  template void
+  create_system_matrix<3, SparseMatrix<double>>(
+    const DoFHandler<3> &                              dof_handler,
+    const QGauss<3> &                                  quad,
+    const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
+    const CDR::Parameters &                            parameters,
+    const double                                       time_step,
+    SparseMatrix<double> &                             system_matrix);
 
-  template
-  void create_system_matrix<2, SparseMatrix<double>>
-  (const DoFHandler<2>                               &dof_handler,
-   const QGauss<2>                                   &quad,
-   const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
-   const CDR::Parameters                             &parameters,
-   const double                                       time_step,
-   const AffineConstraints<double>                   &constraints,
-   SparseMatrix<double>                              &system_matrix);
+  template void
+  create_system_matrix<2, SparseMatrix<double>>(
+    const DoFHandler<2> &                              dof_handler,
+    const QGauss<2> &                                  quad,
+    const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
+    const CDR::Parameters &                            parameters,
+    const double                                       time_step,
+    const AffineConstraints<double> &                  constraints,
+    SparseMatrix<double> &                             system_matrix);
 
-  template
-  void create_system_matrix<3, SparseMatrix<double>>
-  (const DoFHandler<3>                               &dof_handler,
-   const QGauss<3>                                   &quad,
-   const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
-   const CDR::Parameters                             &parameters,
-   const double                                       time_step,
-   const AffineConstraints<double>                   &constraints,
-   SparseMatrix<double>                              &system_matrix);
+  template void
+  create_system_matrix<3, SparseMatrix<double>>(
+    const DoFHandler<3> &                              dof_handler,
+    const QGauss<3> &                                  quad,
+    const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
+    const CDR::Parameters &                            parameters,
+    const double                                       time_step,
+    const AffineConstraints<double> &                  constraints,
+    SparseMatrix<double> &                             system_matrix);
 
-  template
-  void create_system_matrix<2, TrilinosWrappers::SparseMatrix>
-  (const DoFHandler<2>                               &dof_handler,
-   const QGauss<2>                                   &quad,
-   const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
-   const CDR::Parameters                             &parameters,
-   const double                                       time_step,
-   TrilinosWrappers::SparseMatrix                    &system_matrix);
+  template void
+  create_system_matrix<2, TrilinosWrappers::SparseMatrix>(
+    const DoFHandler<2> &                              dof_handler,
+    const QGauss<2> &                                  quad,
+    const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
+    const CDR::Parameters &                            parameters,
+    const double                                       time_step,
+    TrilinosWrappers::SparseMatrix &                   system_matrix);
 
-  template
-  void create_system_matrix<3, TrilinosWrappers::SparseMatrix>
-  (const DoFHandler<3>                               &dof_handler,
-   const QGauss<3>                                   &quad,
-   const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
-   const CDR::Parameters                             &parameters,
-   const double                                       time_step,
-   TrilinosWrappers::SparseMatrix                    &system_matrix);
+  template void
+  create_system_matrix<3, TrilinosWrappers::SparseMatrix>(
+    const DoFHandler<3> &                              dof_handler,
+    const QGauss<3> &                                  quad,
+    const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
+    const CDR::Parameters &                            parameters,
+    const double                                       time_step,
+    TrilinosWrappers::SparseMatrix &                   system_matrix);
 
-  template
-  void create_system_matrix<2, TrilinosWrappers::SparseMatrix>
-  (const DoFHandler<2>                               &dof_handler,
-   const QGauss<2>                                   &quad,
-   const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
-   const CDR::Parameters                             &parameters,
-   const double                                       time_step,
-   const AffineConstraints<double>                   &constraints,
-   TrilinosWrappers::SparseMatrix                    &system_matrix);
+  template void
+  create_system_matrix<2, TrilinosWrappers::SparseMatrix>(
+    const DoFHandler<2> &                              dof_handler,
+    const QGauss<2> &                                  quad,
+    const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
+    const CDR::Parameters &                            parameters,
+    const double                                       time_step,
+    const AffineConstraints<double> &                  constraints,
+    TrilinosWrappers::SparseMatrix &                   system_matrix);
 
-  template
-  void create_system_matrix<3, TrilinosWrappers::SparseMatrix>
-  (const DoFHandler<3>                               &dof_handler,
-   const QGauss<3>                                   &quad,
-   const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
-   const CDR::Parameters                             &parameters,
-   const double                                       time_step,
-   const AffineConstraints<double>                   &constraints,
-   TrilinosWrappers::SparseMatrix                    &system_matrix);
-}
+  template void
+  create_system_matrix<3, TrilinosWrappers::SparseMatrix>(
+    const DoFHandler<3> &                              dof_handler,
+    const QGauss<3> &                                  quad,
+    const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
+    const CDR::Parameters &                            parameters,
+    const double                                       time_step,
+    const AffineConstraints<double> &                  constraints,
+    TrilinosWrappers::SparseMatrix &                   system_matrix);
+} // namespace CDR
diff --git a/cdr/common/source/system_rhs.cc b/cdr/common/source/system_rhs.cc
index df70485..d70e86e 100644
--- a/cdr/common/source/system_rhs.cc
+++ b/cdr/common/source/system_rhs.cc
@@ -11,51 +11,51 @@ namespace CDR
 {
   using namespace dealii;
 
-  template
-  void create_system_rhs<2, Vector<double>>
-  (const DoFHandler<2>                                 &dof_handler,
-   const QGauss<2>                                     &quad,
-   const std::function<Tensor<1, 2>(const Point<2>)>   &convection_function,
-   const std::function<double(double, const Point<2>)> &forcing_function,
-   const CDR::Parameters                               &parameters,
-   const Vector<double>                                &previous_solution,
-   const AffineConstraints<double>                     &constraints,
-   const double                                         current_time,
-   Vector<double>                                      &system_rhs);
+  template void
+  create_system_rhs<2, Vector<double>>(
+    const DoFHandler<2> &                                dof_handler,
+    const QGauss<2> &                                    quad,
+    const std::function<Tensor<1, 2>(const Point<2>)> &  convection_function,
+    const std::function<double(double, const Point<2>)> &forcing_function,
+    const CDR::Parameters &                              parameters,
+    const Vector<double> &                               previous_solution,
+    const AffineConstraints<double> &                    constraints,
+    const double                                         current_time,
+    Vector<double> &                                     system_rhs);
 
-  template
-  void create_system_rhs<3, Vector<double>>
-  (const DoFHandler<3>                                 &dof_handler,
-   const QGauss<3>                                     &quad,
-   const std::function<Tensor<1, 3>(const Point<3>)>   &convection_function,
-   const std::function<double(double, const Point<3>)> &forcing_function,
-   const CDR::Parameters                               &parameters,
-   const Vector<double>                                &previous_solution,
-   const AffineConstraints<double>                     &constraints,
-   const double                                         current_time,
-   Vector<double>                                      &system_rhs);
+  template void
+  create_system_rhs<3, Vector<double>>(
+    const DoFHandler<3> &                                dof_handler,
+    const QGauss<3> &                                    quad,
+    const std::function<Tensor<1, 3>(const Point<3>)> &  convection_function,
+    const std::function<double(double, const Point<3>)> &forcing_function,
+    const CDR::Parameters &                              parameters,
+    const Vector<double> &                               previous_solution,
+    const AffineConstraints<double> &                    constraints,
+    const double                                         current_time,
+    Vector<double> &                                     system_rhs);
 
-  template
-  void create_system_rhs<2, TrilinosWrappers::MPI::Vector>
-  (const DoFHandler<2>                                 &dof_handler,
-   const QGauss<2>                                     &quad,
-   const std::function<Tensor<1, 2>(const Point<2>)>   &convection_function,
-   const std::function<double(double, const Point<2>)> &forcing_function,
-   const CDR::Parameters                               &parameters,
-   const TrilinosWrappers::MPI::Vector                 &previous_solution,
-   const AffineConstraints<double>                     &constraints,
-   const double                                         current_time,
-   TrilinosWrappers::MPI::Vector                       &system_rhs);
+  template void
+  create_system_rhs<2, TrilinosWrappers::MPI::Vector>(
+    const DoFHandler<2> &                                dof_handler,
+    const QGauss<2> &                                    quad,
+    const std::function<Tensor<1, 2>(const Point<2>)> &  convection_function,
+    const std::function<double(double, const Point<2>)> &forcing_function,
+    const CDR::Parameters &                              parameters,
+    const TrilinosWrappers::MPI::Vector &                previous_solution,
+    const AffineConstraints<double> &                    constraints,
+    const double                                         current_time,
+    TrilinosWrappers::MPI::Vector &                      system_rhs);
 
-  template
-  void create_system_rhs<3, TrilinosWrappers::MPI::Vector>
-  (const DoFHandler<3>                                 &dof_handler,
-   const QGauss<3>                                     &quad,
-   const std::function<Tensor<1, 3>(const Point<3>)>   &convection_function,
-   const std::function<double(double, const Point<3>)> &forcing_function,
-   const CDR::Parameters                               &parameters,
-   const TrilinosWrappers::MPI::Vector                 &previous_solution,
-   const AffineConstraints<double>                     &constraints,
-   const double                                         current_time,
-   TrilinosWrappers::MPI::Vector                       &system_rhs);
-}
+  template void
+  create_system_rhs<3, TrilinosWrappers::MPI::Vector>(
+    const DoFHandler<3> &                                dof_handler,
+    const QGauss<3> &                                    quad,
+    const std::function<Tensor<1, 3>(const Point<3>)> &  convection_function,
+    const std::function<double(double, const Point<3>)> &forcing_function,
+    const CDR::Parameters &                              parameters,
+    const TrilinosWrappers::MPI::Vector &                previous_solution,
+    const AffineConstraints<double> &                    constraints,
+    const double                                         current_time,
+    TrilinosWrappers::MPI::Vector &                      system_rhs);
+} // namespace CDR
diff --git a/cdr/common/source/write_pvtu_output.cc b/cdr/common/source/write_pvtu_output.cc
index f0efe08..b918ebb 100644
--- a/cdr/common/source/write_pvtu_output.cc
+++ b/cdr/common/source/write_pvtu_output.cc
@@ -10,31 +10,31 @@ namespace CDR
   using namespace dealii;
 
   WritePVTUOutput::WritePVTUOutput(const unsigned int patch_level)
-    : patch_level {patch_level},
-  this_mpi_process {Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)}
+    : patch_level{patch_level}
+    , this_mpi_process{Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)}
   {}
 
-  template
-  void WritePVTUOutput::write_output(const DoFHandler<2>  &dof_handler,
-                                     const Vector<double> &solution,
-                                     const unsigned int    time_step_n,
-                                     const double          current_time);
+  template void
+  WritePVTUOutput::write_output(const DoFHandler<2> & dof_handler,
+                                const Vector<double> &solution,
+                                const unsigned int    time_step_n,
+                                const double          current_time);
 
-  template
-  void WritePVTUOutput::write_output(const DoFHandler<3>  &dof_handler,
-                                     const Vector<double> &solution,
-                                     const unsigned int    time_step_n,
-                                     const double          current_time);
+  template void
+  WritePVTUOutput::write_output(const DoFHandler<3> & dof_handler,
+                                const Vector<double> &solution,
+                                const unsigned int    time_step_n,
+                                const double          current_time);
 
-  template
-  void WritePVTUOutput::write_output(const DoFHandler<2>                 &dof_handler,
-                                     const TrilinosWrappers::MPI::Vector &solution,
-                                     const unsigned int                   time_step_n,
-                                     const double                         current_time);
+  template void
+  WritePVTUOutput::write_output(const DoFHandler<2> &dof_handler,
+                                const TrilinosWrappers::MPI::Vector &solution,
+                                const unsigned int time_step_n,
+                                const double       current_time);
 
-  template
-  void WritePVTUOutput::write_output(const DoFHandler<3>                 &dof_handler,
-                                     const TrilinosWrappers::MPI::Vector &solution,
-                                     const unsigned int                   time_step_n,
-                                     const double                         current_time);
-}
+  template void
+  WritePVTUOutput::write_output(const DoFHandler<3> &dof_handler,
+                                const TrilinosWrappers::MPI::Vector &solution,
+                                const unsigned int time_step_n,
+                                const double       current_time);
+} // namespace CDR
diff --git a/cdr/solver/cdr.cc b/cdr/solver/cdr.cc
index f5b35d3..4b62f1e 100644
--- a/cdr/solver/cdr.cc
+++ b/cdr/solver/cdr.cc
@@ -9,59 +9,62 @@
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/manifold_lib.h>
 
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
 #include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
 
 #include <deal.II/numerics/error_estimator.h>
 
 // These headers are for distributed computations:
-#include <deal.II/base/utilities.h>
 #include <deal.II/base/index_set.h>
-#include <deal.II/distributed/tria.h>
+#include <deal.II/base/utilities.h>
+
 #include <deal.II/distributed/grid_refinement.h>
 #include <deal.II/distributed/solution_transfer.h>
-#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/distributed/tria.h>
 
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
 #include <deal.II/lac/trilinos_solver.h>
 #include <deal.II/lac/trilinos_sparse_matrix.h>
-#include <deal.II/lac/trilinos_precondition.h>
 #include <deal.II/lac/trilinos_vector.h>
 
-#include <chrono>
-#include <functional>
-#include <iostream>
-
+#include <deal.II-cdr/parameters.h>
 #include <deal.II-cdr/system_matrix.h>
 #include <deal.II-cdr/system_rhs.h>
-#include <deal.II-cdr/parameters.h>
 #include <deal.II-cdr/write_pvtu_output.h>
 
+#include <chrono>
+#include <functional>
+#include <iostream>
+
 using namespace dealii;
 
-constexpr int manifold_id {0};
+constexpr int manifold_id{0};
 
 // This is the actual solver class which performs time iteration and calls the
 // appropriate library functions to do it.
-template<int dim>
+template <int dim>
 class CDRProblem
 {
 public:
   CDRProblem(const CDR::Parameters &parameters);
-  void run();
+  void
+  run();
+
 private:
   const CDR::Parameters parameters;
-  const double time_step;
-  double current_time;
+  const double          time_step;
+  double                current_time;
 
-  MPI_Comm mpi_communicator;
+  MPI_Comm           mpi_communicator;
   const unsigned int n_mpi_processes;
   const unsigned int this_mpi_process;
 
-  FE_Q<dim> fe;
-  QGauss<dim> quad;
-  const SphericalManifold<dim> boundary_description;
+  FE_Q<dim>                                 fe;
+  QGauss<dim>                               quad;
+  const SphericalManifold<dim>              boundary_description;
   parallel::distributed::Triangulation<dim> triangulation;
-  DoFHandler<dim> dof_handler;
+  DoFHandler<dim>                           dof_handler;
 
   const std::function<Tensor<1, dim>(const Point<dim>)> convection_function;
   const std::function<double(double, const Point<dim>)> forcing_function;
@@ -70,7 +73,7 @@ private:
   IndexSet locally_relevant_dofs;
 
   AffineConstraints<double> constraints;
-  bool first_run;
+  bool                      first_run;
 
   // As is usual in parallel programs, I keep two copies of parts of the
   // complete solution: <code>locally_relevant_solution</code> contains both
@@ -78,64 +81,70 @@ private:
   // boundary (the @ref GlossGhostCell "ghost cells") while
   // <code>completely_distributed_solution</code> only contains the parts of
   // the solution computed on the current @ref GlossMPIProcess "MPI process".
-  TrilinosWrappers::MPI::Vector locally_relevant_solution;
-  TrilinosWrappers::MPI::Vector completely_distributed_solution;
-  TrilinosWrappers::MPI::Vector system_rhs;
-  TrilinosWrappers::SparseMatrix system_matrix;
+  TrilinosWrappers::MPI::Vector     locally_relevant_solution;
+  TrilinosWrappers::MPI::Vector     completely_distributed_solution;
+  TrilinosWrappers::MPI::Vector     system_rhs;
+  TrilinosWrappers::SparseMatrix    system_matrix;
   TrilinosWrappers::PreconditionAMG preconditioner;
 
   ConditionalOStream pcout;
 
-  void setup_geometry();
-  void setup_system();
-  void setup_dofs();
-  void refine_mesh();
-  void time_iterate();
+  void
+  setup_geometry();
+  void
+  setup_system();
+  void
+  setup_dofs();
+  void
+  refine_mesh();
+  void
+  time_iterate();
 };
 
 
-template<int dim>
-CDRProblem<dim>::CDRProblem(const CDR::Parameters &parameters) :
-  parameters(parameters),
-  time_step {(parameters.stop_time - parameters.start_time)
-  /parameters.n_time_steps
-},
-current_time {parameters.start_time},
-mpi_communicator (MPI_COMM_WORLD),
-n_mpi_processes {Utilities::MPI::n_mpi_processes(mpi_communicator)},
-this_mpi_process {Utilities::MPI::this_mpi_process(mpi_communicator)},
-fe(parameters.fe_order),
-quad(parameters.fe_order + 2),
-boundary_description(Point<dim>()),
-triangulation(mpi_communicator, typename Triangulation<dim>::MeshSmoothing
-              (Triangulation<dim>::smoothing_on_refinement |
-               Triangulation<dim>::smoothing_on_coarsening)),
-dof_handler(triangulation),
-convection_function
-{
-  [](const Point<dim> p) -> Tensor<1, dim>
-  {Tensor<1, dim> v; v[0] = -p[1]; v[1] = p[0]; return v;}
-},
-forcing_function
-{
-  [](double t, const Point<dim> p) -> double
-  {
-    return std::exp(-8*t)*std::exp(-40*Utilities::fixed_power<6>(p[0] - 1.5))
-    *std::exp(-40*Utilities::fixed_power<6>(p[1]));
-  }
-},
-first_run {true},
-pcout (std::cout, this_mpi_process == 0)
+template <int dim>
+CDRProblem<dim>::CDRProblem(const CDR::Parameters &parameters)
+  : parameters(parameters)
+  , time_step{(parameters.stop_time - parameters.start_time) /
+              parameters.n_time_steps}
+  , current_time{parameters.start_time}
+  , mpi_communicator(MPI_COMM_WORLD)
+  , n_mpi_processes{Utilities::MPI::n_mpi_processes(mpi_communicator)}
+  , this_mpi_process{Utilities::MPI::this_mpi_process(mpi_communicator)}
+  , fe(parameters.fe_order)
+  , quad(parameters.fe_order + 2)
+  , boundary_description(Point<dim>())
+  , triangulation(mpi_communicator,
+                  typename Triangulation<dim>::MeshSmoothing(
+                    Triangulation<dim>::smoothing_on_refinement |
+                    Triangulation<dim>::smoothing_on_coarsening))
+  , dof_handler(triangulation)
+  , convection_function{[](const Point<dim> p) -> Tensor<1, dim> {
+    Tensor<1, dim> v;
+    v[0] = -p[1];
+    v[1] = p[0];
+    return v;
+  }}
+  , forcing_function{[](double t, const Point<dim> p) -> double {
+    return std::exp(-8 * t) *
+           std::exp(-40 * Utilities::fixed_power<6>(p[0] - 1.5)) *
+           std::exp(-40 * Utilities::fixed_power<6>(p[1]));
+  }}
+  , first_run{true}
+  , pcout(std::cout, this_mpi_process == 0)
 {
   Assert(dim == 2, ExcNotImplemented());
 }
 
 
-template<int dim>
-void CDRProblem<dim>::setup_geometry()
+template <int dim>
+void
+CDRProblem<dim>::setup_geometry()
 {
   const Point<dim> center;
-  GridGenerator::hyper_shell(triangulation, center, parameters.inner_radius,
+  GridGenerator::hyper_shell(triangulation,
+                             center,
+                             parameters.inner_radius,
                              parameters.outer_radius);
   triangulation.set_manifold(manifold_id, boundary_description);
   for (const auto &cell : triangulation.active_cell_iterators())
@@ -146,12 +155,12 @@ void CDRProblem<dim>::setup_geometry()
 }
 
 
-template<int dim>
-void CDRProblem<dim>::setup_dofs()
+template <int dim>
+void
+CDRProblem<dim>::setup_dofs()
 {
   dof_handler.distribute_dofs(fe);
-  pcout << "Number of degrees of freedom: "
-        << dof_handler.n_dofs()
+  pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
         << std::endl;
   locally_owned_dofs = dof_handler.locally_owned_dofs();
   DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
@@ -159,43 +168,55 @@ void CDRProblem<dim>::setup_dofs()
   constraints.clear();
   constraints.reinit(locally_relevant_dofs);
   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
-  DoFTools::make_zero_boundary_constraints(dof_handler, manifold_id, constraints);
+  DoFTools::make_zero_boundary_constraints(dof_handler,
+                                           manifold_id,
+                                           constraints);
   constraints.close();
 
-  completely_distributed_solution.reinit
-  (locally_owned_dofs, mpi_communicator);
+  completely_distributed_solution.reinit(locally_owned_dofs, mpi_communicator);
 
-  locally_relevant_solution.reinit(locally_owned_dofs, locally_relevant_dofs,
+  locally_relevant_solution.reinit(locally_owned_dofs,
+                                   locally_relevant_dofs,
                                    mpi_communicator);
 }
 
 
-template<int dim>
-void CDRProblem<dim>::setup_system()
+template <int dim>
+void
+CDRProblem<dim>::setup_system()
 {
   DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs());
-  DoFTools::make_sparsity_pattern(dof_handler, dynamic_sparsity_pattern,
-                                  constraints, /*keep_constrained_dofs*/true);
-  SparsityTools::distribute_sparsity_pattern
-  (dynamic_sparsity_pattern, dof_handler.locally_owned_dofs(),
-   mpi_communicator, locally_relevant_dofs);
+  DoFTools::make_sparsity_pattern(dof_handler,
+                                  dynamic_sparsity_pattern,
+                                  constraints,
+                                  /*keep_constrained_dofs*/ true);
+  SparsityTools::distribute_sparsity_pattern(dynamic_sparsity_pattern,
+                                             dof_handler.locally_owned_dofs(),
+                                             mpi_communicator,
+                                             locally_relevant_dofs);
 
   system_rhs.reinit(locally_owned_dofs, mpi_communicator);
-  system_matrix.reinit(locally_owned_dofs, dynamic_sparsity_pattern,
+  system_matrix.reinit(locally_owned_dofs,
+                       dynamic_sparsity_pattern,
                        mpi_communicator);
 
-  CDR::create_system_matrix<dim>
-  (dof_handler, quad, convection_function, parameters, time_step, constraints,
-   system_matrix);
+  CDR::create_system_matrix<dim>(dof_handler,
+                                 quad,
+                                 convection_function,
+                                 parameters,
+                                 time_step,
+                                 constraints,
+                                 system_matrix);
   system_matrix.compress(VectorOperation::add);
   preconditioner.initialize(system_matrix);
 }
 
 
-template<int dim>
-void CDRProblem<dim>::time_iterate()
+template <int dim>
+void
+CDRProblem<dim>::time_iterate()
 {
-  double current_time = parameters.start_time;
+  double               current_time = parameters.start_time;
   CDR::WritePVTUOutput pvtu_output(parameters.patch_level);
   for (unsigned int time_step_n = 0; time_step_n < parameters.n_time_steps;
        ++time_step_n)
@@ -203,25 +224,35 @@ void CDRProblem<dim>::time_iterate()
       current_time += time_step;
 
       system_rhs = 0.0;
-      CDR::create_system_rhs<dim>
-      (dof_handler, quad, convection_function, forcing_function, parameters,
-       locally_relevant_solution, constraints, current_time, system_rhs);
+      CDR::create_system_rhs<dim>(dof_handler,
+                                  quad,
+                                  convection_function,
+                                  forcing_function,
+                                  parameters,
+                                  locally_relevant_solution,
+                                  constraints,
+                                  current_time,
+                                  system_rhs);
       system_rhs.compress(VectorOperation::add);
 
-      SolverControl solver_control(dof_handler.n_dofs(),
-                                   1e-6*system_rhs.l2_norm(),
+      SolverControl                 solver_control(dof_handler.n_dofs(),
+                                   1e-6 * system_rhs.l2_norm(),
                                    /*log_history = */ false,
                                    /*log_result = */ false);
       TrilinosWrappers::SolverGMRES solver(solver_control);
-      solver.solve(system_matrix, completely_distributed_solution, system_rhs,
+      solver.solve(system_matrix,
+                   completely_distributed_solution,
+                   system_rhs,
                    preconditioner);
       constraints.distribute(completely_distributed_solution);
       locally_relevant_solution = completely_distributed_solution;
 
       if (time_step_n % parameters.save_interval == 0)
         {
-          pvtu_output.write_output(dof_handler, locally_relevant_solution,
-                                   time_step_n, current_time);
+          pvtu_output.write_output(dof_handler,
+                                   locally_relevant_solution,
+                                   time_step_n,
+                                   current_time);
         }
 
       refine_mesh();
@@ -229,16 +260,18 @@ void CDRProblem<dim>::time_iterate()
 }
 
 
-template<int dim>
-void CDRProblem<dim>::refine_mesh()
+template <int dim>
+void
+CDRProblem<dim>::refine_mesh()
 {
-  using FunctionMap =
-    std::map<types::boundary_id, const Function<dim> *>;
+  using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
 
   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
-  KellyErrorEstimator<dim>::estimate
-  (dof_handler, QGauss<dim - 1>(fe.degree + 1), FunctionMap(),
-   locally_relevant_solution, estimated_error_per_cell);
+  KellyErrorEstimator<dim>::estimate(dof_handler,
+                                     QGauss<dim - 1>(fe.degree + 1),
+                                     FunctionMap(),
+                                     locally_relevant_solution,
+                                     estimated_error_per_cell);
 
   // This solver uses a crude refinement strategy where cells with relatively
   // high errors are refined and cells with relatively low errors are
@@ -250,7 +283,8 @@ void CDRProblem<dim>::refine_mesh()
         {
           cell->set_refine_flag();
         }
-      else if (std::abs(estimated_error_per_cell[cell->active_cell_index()]) <= 1e-5)
+      else if (std::abs(estimated_error_per_cell[cell->active_cell_index()]) <=
+               1e-5)
         {
           cell->set_coarsen_flag();
         }
@@ -258,8 +292,8 @@ void CDRProblem<dim>::refine_mesh()
 
   if (triangulation.n_levels() > parameters.max_refinement_level)
     {
-      for (const auto &cell :
-           triangulation.cell_iterators_on_level(parameters.max_refinement_level))
+      for (const auto &cell : triangulation.cell_iterators_on_level(
+             parameters.max_refinement_level))
         {
           cell->clear_refine_flag();
         }
@@ -268,11 +302,11 @@ void CDRProblem<dim>::refine_mesh()
   // Transferring the solution between different grids is ultimately just a
   // few function calls but they must be made in exactly the right order.
   parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::Vector>
-  solution_transfer(dof_handler);
+    solution_transfer(dof_handler);
 
   triangulation.prepare_coarsening_and_refinement();
-  solution_transfer.prepare_for_coarsening_and_refinement
-  (locally_relevant_solution);
+  solution_transfer.prepare_for_coarsening_and_refinement(
+    locally_relevant_solution);
   triangulation.execute_coarsening_and_refinement();
 
   setup_dofs();
@@ -281,8 +315,7 @@ void CDRProblem<dim>::refine_mesh()
   // <code>locally_relevant_solution</code>, so when
   // parallel::distributed::SolutionTransfer::interpolate is called it uses
   // those values to populate <code>temporary</code>.
-  TrilinosWrappers::MPI::Vector temporary
-  (locally_owned_dofs, mpi_communicator);
+  TrilinosWrappers::MPI::Vector temporary(locally_owned_dofs, mpi_communicator);
   solution_transfer.interpolate(temporary);
   // After <code>temporary</code> has the correct value, this call correctly
   // populates <code>completely_distributed_solution</code>, which had its
@@ -298,8 +331,9 @@ void CDRProblem<dim>::refine_mesh()
 }
 
 
-template<int dim>
-void CDRProblem<dim>::run()
+template <int dim>
+void
+CDRProblem<dim>::run()
 {
   setup_geometry();
   setup_dofs();
@@ -308,17 +342,18 @@ void CDRProblem<dim>::run()
 }
 
 
-constexpr int dim {2};
+constexpr int dim{2};
 
 
-int main(int argc, char *argv[])
+int
+main(int argc, char *argv[])
 {
   // One of the new features in C++11 is the <code>chrono</code> component of
   // the standard library. This gives us an easy way to time the output.
   auto t0 = std::chrono::high_resolution_clock::now();
 
   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
-  CDR::Parameters parameters;
+  CDR::Parameters                  parameters;
   parameters.read_parameter_file("parameters.prm");
   CDRProblem<dim> cdr_problem(parameters);
   cdr_problem.run();
@@ -327,9 +362,10 @@ int main(int argc, char *argv[])
   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
     {
       std::cout << "time elapsed: "
-                << std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0).count()
-                << " milliseconds."
-                << std::endl;
+                << std::chrono::duration_cast<std::chrono::milliseconds>(t1 -
+                                                                         t0)
+                     .count()
+                << " milliseconds." << std::endl;
     }
 
   return 0;