]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Avoid the use of MeshWorker::LocalIntegrator and MeshWorker::integration_loop(). 17065/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Tue, 21 May 2024 02:41:42 +0000 (20:41 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Thu, 23 May 2024 03:18:39 +0000 (21:18 -0600)
examples/step-39/step-39.cc

index 105f692e555a532be6f4e33f883b8267eff43956..b423668346796c88adaac3f5988c78603e918fc5 100644 (file)
@@ -76,35 +76,26 @@ namespace Step39
 
   // @sect3{The local integrators}
 
-  // MeshWorker separates local integration from the loops over cells and
-  // faces. Thus, we have to write local integration classes for generating
-  // matrices, the right hand side and the error estimator.
-
-  // All these classes have the same three functions for integrating over
-  // cells, boundary faces and interior faces, respectively. All the
-  // information needed for the local integration is provided by
-  // MeshWorker::IntegrationInfo<dim>. Note that the signature of the
-  // functions cannot be changed, because it is expected by
-  // MeshWorker::integration_loop().
-
-  // The first class defining local integrators is responsible for computing
-  // cell and face matrices. It is used to assemble the global matrix as well
-  // as the level matrices.
-  template <int dim>
-  class MatrixIntegrator : public MeshWorker::LocalIntegrator<dim>
-  {
-  public:
-    void cell(MeshWorker::DoFInfo<dim>         &dinfo,
-              MeshWorker::IntegrationInfo<dim> &info) const override;
-    void boundary(MeshWorker::DoFInfo<dim>         &dinfo,
-                  MeshWorker::IntegrationInfo<dim> &info) const override;
-    void face(MeshWorker::DoFInfo<dim>         &dinfo1,
-              MeshWorker::DoFInfo<dim>         &dinfo2,
-              MeshWorker::IntegrationInfo<dim> &info1,
-              MeshWorker::IntegrationInfo<dim> &info2) const override;
-  };
-
-
+  // The MeshWorker::loop() function separates what needs to be done for
+  // local integration, from the loops over cells and
+  // faces. It does this by calling functions that integrate over a cell,
+  // a boundary face, or an interior face, and letting them create the
+  // local contributions and then in a separate step calling a function
+  // that moves these local contributions into the global objects.
+  // We will use this approach for computing the
+  // matrices, the right hand side, the error estimator, and the actual
+  // error computation in the functions below. For each of these operations,
+  // we provide a namespace that contains a set of functions for cell, boundary,
+  // and interior face contributions.
+  //
+  // All the information needed for these local integration is provided by
+  // MeshWorker::DoFInfo<dim> and MeshWorker::IntegrationInfo<dim>. In each
+  // case, the functions' signatures is fixed: MeshWorker::loop() wants to call
+  // functions with a specific set of arguments, so the signature of the
+  // functions cannot be changed.
+
+  // The first namespace defining local integrators is responsible for
+  // assembling the global matrix as well as the level matrices.
   // On each cell, we integrate the Dirichlet form as well as the
   // Nitsche boundary conditions and the interior penalty fluxes between
   // cells.
@@ -121,317 +112,292 @@ namespace Step39
   // cells are refined, in which case we are integrating over a non-active
   // face and no adjustment is necessary. Finally, we return the average
   // of the two penalty values.
-  template <int dim>
-  double ip_penalty_factor(const MeshWorker::DoFInfo<dim> &dinfo1,
-                           const MeshWorker::DoFInfo<dim> &dinfo2,
-                           unsigned int                    deg1,
-                           unsigned int                    deg2)
-  {
-    const unsigned int normal1 =
-      GeometryInfo<dim>::unit_normal_direction[dinfo1.face_number];
-    const unsigned int normal2 =
-      GeometryInfo<dim>::unit_normal_direction[dinfo2.face_number];
-    const unsigned int deg1sq = (deg1 == 0) ? 1 : deg1 * (deg1 + 1);
-    const unsigned int deg2sq = (deg2 == 0) ? 1 : deg2 * (deg2 + 1);
-
-    double penalty1 = deg1sq / dinfo1.cell->extent_in_direction(normal1);
-    double penalty2 = deg2sq / dinfo2.cell->extent_in_direction(normal2);
-    if (dinfo1.cell->has_children() && !dinfo2.cell->has_children())
-      penalty1 *= 2;
-    else if (!dinfo1.cell->has_children() && dinfo2.cell->has_children())
-      penalty2 *= 2;
-
-    const double penalty = 0.5 * (penalty1 + penalty2);
-    return penalty;
-  }
-
-
-  template <int dim>
-  void MatrixIntegrator<dim>::cell(MeshWorker::DoFInfo<dim>         &dinfo,
-                                   MeshWorker::IntegrationInfo<dim> &info) const
-  {
-    FullMatrix<double> &M = dinfo.matrix(0, false).matrix;
-
-    for (unsigned int k = 0; k < info.fe_values().n_quadrature_points; ++k)
-      {
-        const double dx = info.fe_values().JxW(k);
-
-        for (unsigned int i = 0; i < info.fe_values().dofs_per_cell; ++i)
-          {
-            const double Mii = (info.fe_values().shape_grad(i, k) *
-                                info.fe_values().shape_grad(i, k) * dx);
-
-            M(i, i) += Mii;
-
-            for (unsigned int j = i + 1; j < info.fe_values().dofs_per_cell;
-                 ++j)
-              {
-                const double Mij = info.fe_values().shape_grad(j, k) *
-                                   info.fe_values().shape_grad(i, k) * dx;
-
-                M(i, j) += Mij;
-                M(j, i) += Mij;
-              }
-          }
-      }
-  }
-
-
-  // Boundary faces use the Nitsche method to impose boundary values:
-  template <int dim>
-  void
-  MatrixIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim>         &dinfo,
-                                  MeshWorker::IntegrationInfo<dim> &info) const
+  namespace MatrixIntegrator
   {
-    const FEValuesBase<dim> &fe_face_values = info.fe_values(0);
-
-    FullMatrix<double> &M = dinfo.matrix(0, false).matrix;
-    AssertDimension(M.n(), fe_face_values.dofs_per_cell);
-    AssertDimension(M.m(), fe_face_values.dofs_per_cell);
-
-    const unsigned int polynomial_degree =
-      info.fe_values(0).get_fe().tensor_degree();
+    template <int dim>
+    double ip_penalty_factor(const MeshWorker::DoFInfo<dim> &dinfo1,
+                             const MeshWorker::DoFInfo<dim> &dinfo2,
+                             unsigned int                    deg1,
+                             unsigned int                    deg2)
+    {
+      const unsigned int normal1 =
+        GeometryInfo<dim>::unit_normal_direction[dinfo1.face_number];
+      const unsigned int normal2 =
+        GeometryInfo<dim>::unit_normal_direction[dinfo2.face_number];
+      const unsigned int deg1sq = (deg1 == 0) ? 1 : deg1 * (deg1 + 1);
+      const unsigned int deg2sq = (deg2 == 0) ? 1 : deg2 * (deg2 + 1);
+
+      double penalty1 = deg1sq / dinfo1.cell->extent_in_direction(normal1);
+      double penalty2 = deg2sq / dinfo2.cell->extent_in_direction(normal2);
+      if (dinfo1.cell->has_children() && !dinfo2.cell->has_children())
+        penalty1 *= 2;
+      else if (!dinfo1.cell->has_children() && dinfo2.cell->has_children())
+        penalty2 *= 2;
+
+      const double penalty = 0.5 * (penalty1 + penalty2);
+      return penalty;
+    }
 
-    const double ip_penalty =
-      ip_penalty_factor(dinfo, dinfo, polynomial_degree, polynomial_degree);
 
-    for (unsigned int k = 0; k < fe_face_values.n_quadrature_points; ++k)
-      {
-        const double          dx = fe_face_values.JxW(k);
-        const Tensor<1, dim> &n  = fe_face_values.normal_vector(k);
-
-        for (unsigned int i = 0; i < fe_face_values.dofs_per_cell; ++i)
-          for (unsigned int j = 0; j < fe_face_values.dofs_per_cell; ++j)
-            M(i, j) += (2. * fe_face_values.shape_value(i, k) * ip_penalty *
-                          fe_face_values.shape_value(j, k) -
-                        (n * fe_face_values.shape_grad(i, k)) *
-                          fe_face_values.shape_value(j, k) -
-                        (n * fe_face_values.shape_grad(j, k)) *
-                          fe_face_values.shape_value(i, k)) *
-                       dx;
-      }
-  }
-
-  // Interior faces use the interior penalty method:
-  template <int dim>
-  void
-  MatrixIntegrator<dim>::face(MeshWorker::DoFInfo<dim>         &dinfo1,
-                              MeshWorker::DoFInfo<dim>         &dinfo2,
-                              MeshWorker::IntegrationInfo<dim> &info1,
-                              MeshWorker::IntegrationInfo<dim> &info2) const
-  {
-    const FEValuesBase<dim> &fe_face_values_1 = info1.fe_values(0);
-    const FEValuesBase<dim> &fe_face_values_2 = info2.fe_values(0);
-
-    FullMatrix<double> &M11 = dinfo1.matrix(0, false).matrix;
-    FullMatrix<double> &M12 = dinfo1.matrix(0, true).matrix;
-    FullMatrix<double> &M21 = dinfo2.matrix(0, true).matrix;
-    FullMatrix<double> &M22 = dinfo2.matrix(0, false).matrix;
-
-    AssertDimension(M11.n(), fe_face_values_1.dofs_per_cell);
-    AssertDimension(M11.m(), fe_face_values_1.dofs_per_cell);
-    AssertDimension(M12.n(), fe_face_values_1.dofs_per_cell);
-    AssertDimension(M12.m(), fe_face_values_1.dofs_per_cell);
-    AssertDimension(M21.n(), fe_face_values_1.dofs_per_cell);
-    AssertDimension(M21.m(), fe_face_values_1.dofs_per_cell);
-    AssertDimension(M22.n(), fe_face_values_1.dofs_per_cell);
-    AssertDimension(M22.m(), fe_face_values_1.dofs_per_cell);
-
-    const unsigned int polynomial_degree =
-      info1.fe_values(0).get_fe().tensor_degree();
-    const double ip_penalty =
-      ip_penalty_factor(dinfo1, dinfo2, polynomial_degree, polynomial_degree);
-
-    const double nui = 1.;
-    const double nue = 1.;
-    const double nu  = .5 * (nui + nue);
-
-    for (unsigned int k = 0; k < fe_face_values_1.n_quadrature_points; ++k)
-      {
-        const double          dx = fe_face_values_1.JxW(k);
-        const Tensor<1, dim> &n  = fe_face_values_1.normal_vector(k);
+    template <int dim>
+    void cell(MeshWorker::DoFInfo<dim>         &dinfo,
+              MeshWorker::IntegrationInfo<dim> &info)
+    {
+      FullMatrix<double> &M = dinfo.matrix(0, false).matrix;
+
+      for (unsigned int k = 0; k < info.fe_values().n_quadrature_points; ++k)
+        {
+          const double dx = info.fe_values().JxW(k);
+
+          for (unsigned int i = 0; i < info.fe_values().dofs_per_cell; ++i)
+            {
+              const double Mii = (info.fe_values().shape_grad(i, k) *
+                                  info.fe_values().shape_grad(i, k) * dx);
+
+              M(i, i) += Mii;
+
+              for (unsigned int j = i + 1; j < info.fe_values().dofs_per_cell;
+                   ++j)
+                {
+                  const double Mij = info.fe_values().shape_grad(j, k) *
+                                     info.fe_values().shape_grad(i, k) * dx;
+
+                  M(i, j) += Mij;
+                  M(j, i) += Mij;
+                }
+            }
+        }
+    }
 
-        for (unsigned int i = 0; i < fe_face_values_1.dofs_per_cell; ++i)
-          {
-            for (unsigned int j = 0; j < fe_face_values_1.dofs_per_cell; ++j)
-              {
-                const double vi   = fe_face_values_1.shape_value(i, k);
-                const double dnvi = n * fe_face_values_1.shape_grad(i, k);
-                const double ve   = fe_face_values_2.shape_value(i, k);
-                const double dnve = n * fe_face_values_2.shape_grad(i, k);
-                const double ui   = fe_face_values_1.shape_value(j, k);
-                const double dnui = n * fe_face_values_1.shape_grad(j, k);
-                const double ue   = fe_face_values_2.shape_value(j, k);
-                const double dnue = n * fe_face_values_2.shape_grad(j, k);
-
-                M11(i, j) += (-.5 * nui * dnvi * ui - .5 * nui * dnui * vi +
-                              nu * ip_penalty * ui * vi) *
-                             dx;
-                M12(i, j) += (.5 * nui * dnvi * ue - .5 * nue * dnue * vi -
-                              nu * ip_penalty * vi * ue) *
-                             dx;
-                M21(i, j) += (-.5 * nue * dnve * ui + .5 * nui * dnui * ve -
-                              nu * ip_penalty * ui * ve) *
-                             dx;
-                M22(i, j) += (.5 * nue * dnve * ue + .5 * nue * dnue * ve +
-                              nu * ip_penalty * ue * ve) *
-                             dx;
-              }
-          }
-      }
-  }
 
-  // The second local integrator builds the right hand side. In our example,
-  // the right hand side function is zero, such that only the boundary
-  // condition is set here in weak form.
-  template <int dim>
-  class RHSIntegrator : public MeshWorker::LocalIntegrator<dim>
-  {
-  public:
-    void cell(MeshWorker::DoFInfo<dim>         &dinfo,
-              MeshWorker::IntegrationInfo<dim> &info) const override;
+    // Boundary faces use the Nitsche method to impose boundary values:
+    template <int dim>
     void boundary(MeshWorker::DoFInfo<dim>         &dinfo,
-                  MeshWorker::IntegrationInfo<dim> &info) const override;
+                  MeshWorker::IntegrationInfo<dim> &info)
+    {
+      const FEValuesBase<dim> &fe_face_values = info.fe_values(0);
+
+      FullMatrix<double> &M = dinfo.matrix(0, false).matrix;
+      AssertDimension(M.n(), fe_face_values.dofs_per_cell);
+      AssertDimension(M.m(), fe_face_values.dofs_per_cell);
+
+      const unsigned int polynomial_degree =
+        info.fe_values(0).get_fe().tensor_degree();
+
+      const double ip_penalty =
+        ip_penalty_factor(dinfo, dinfo, polynomial_degree, polynomial_degree);
+
+      for (unsigned int k = 0; k < fe_face_values.n_quadrature_points; ++k)
+        {
+          const double          dx = fe_face_values.JxW(k);
+          const Tensor<1, dim> &n  = fe_face_values.normal_vector(k);
+
+          for (unsigned int i = 0; i < fe_face_values.dofs_per_cell; ++i)
+            for (unsigned int j = 0; j < fe_face_values.dofs_per_cell; ++j)
+              M(i, j) += (2. * fe_face_values.shape_value(i, k) * ip_penalty *
+                            fe_face_values.shape_value(j, k) -
+                          (n * fe_face_values.shape_grad(i, k)) *
+                            fe_face_values.shape_value(j, k) -
+                          (n * fe_face_values.shape_grad(j, k)) *
+                            fe_face_values.shape_value(i, k)) *
+                         dx;
+        }
+    }
+
+    // Interior faces use the interior penalty method:
+    template <int dim>
     void face(MeshWorker::DoFInfo<dim>         &dinfo1,
               MeshWorker::DoFInfo<dim>         &dinfo2,
               MeshWorker::IntegrationInfo<dim> &info1,
-              MeshWorker::IntegrationInfo<dim> &info2) const override;
-  };
-
-
-  template <int dim>
-  void RHSIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> &,
-                                MeshWorker::IntegrationInfo<dim> &) const
-  {}
-
+              MeshWorker::IntegrationInfo<dim> &info2)
+    {
+      const FEValuesBase<dim> &fe_face_values_1 = info1.fe_values(0);
+      const FEValuesBase<dim> &fe_face_values_2 = info2.fe_values(0);
+
+      FullMatrix<double> &M11 = dinfo1.matrix(0, false).matrix;
+      FullMatrix<double> &M12 = dinfo1.matrix(0, true).matrix;
+      FullMatrix<double> &M21 = dinfo2.matrix(0, true).matrix;
+      FullMatrix<double> &M22 = dinfo2.matrix(0, false).matrix;
+
+      AssertDimension(M11.n(), fe_face_values_1.dofs_per_cell);
+      AssertDimension(M11.m(), fe_face_values_1.dofs_per_cell);
+      AssertDimension(M12.n(), fe_face_values_1.dofs_per_cell);
+      AssertDimension(M12.m(), fe_face_values_1.dofs_per_cell);
+      AssertDimension(M21.n(), fe_face_values_1.dofs_per_cell);
+      AssertDimension(M21.m(), fe_face_values_1.dofs_per_cell);
+      AssertDimension(M22.n(), fe_face_values_1.dofs_per_cell);
+      AssertDimension(M22.m(), fe_face_values_1.dofs_per_cell);
+
+      const unsigned int polynomial_degree =
+        info1.fe_values(0).get_fe().tensor_degree();
+      const double ip_penalty =
+        ip_penalty_factor(dinfo1, dinfo2, polynomial_degree, polynomial_degree);
+
+      const double nui = 1.;
+      const double nue = 1.;
+      const double nu  = .5 * (nui + nue);
+
+      for (unsigned int k = 0; k < fe_face_values_1.n_quadrature_points; ++k)
+        {
+          const double          dx = fe_face_values_1.JxW(k);
+          const Tensor<1, dim> &n  = fe_face_values_1.normal_vector(k);
+
+          for (unsigned int i = 0; i < fe_face_values_1.dofs_per_cell; ++i)
+            {
+              for (unsigned int j = 0; j < fe_face_values_1.dofs_per_cell; ++j)
+                {
+                  const double vi   = fe_face_values_1.shape_value(i, k);
+                  const double dnvi = n * fe_face_values_1.shape_grad(i, k);
+                  const double ve   = fe_face_values_2.shape_value(i, k);
+                  const double dnve = n * fe_face_values_2.shape_grad(i, k);
+                  const double ui   = fe_face_values_1.shape_value(j, k);
+                  const double dnui = n * fe_face_values_1.shape_grad(j, k);
+                  const double ue   = fe_face_values_2.shape_value(j, k);
+                  const double dnue = n * fe_face_values_2.shape_grad(j, k);
+
+                  M11(i, j) += (-.5 * nui * dnvi * ui - .5 * nui * dnui * vi +
+                                nu * ip_penalty * ui * vi) *
+                               dx;
+                  M12(i, j) += (.5 * nui * dnvi * ue - .5 * nue * dnue * vi -
+                                nu * ip_penalty * vi * ue) *
+                               dx;
+                  M21(i, j) += (-.5 * nue * dnve * ui + .5 * nui * dnui * ve -
+                                nu * ip_penalty * ui * ve) *
+                               dx;
+                  M22(i, j) += (.5 * nue * dnve * ue + .5 * nue * dnue * ve +
+                                nu * ip_penalty * ue * ve) *
+                               dx;
+                }
+            }
+        }
+    }
+  } // namespace MatrixIntegrator
 
-  template <int dim>
-  void
-  RHSIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim>         &dinfo,
-                               MeshWorker::IntegrationInfo<dim> &info) const
+  // The second set of local integrators builds the right hand side. In our
+  // example, the right hand side function is zero, such that only the boundary
+  // condition is set here in weak form.
+  namespace RHSIntegrator
   {
-    const FEValuesBase<dim> &fe           = info.fe_values();
-    Vector<double>          &local_vector = dinfo.vector(0).block(0);
-
-    std::vector<double> boundary_values(fe.n_quadrature_points);
-    exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
-
-    const unsigned int degree = fe.get_fe().tensor_degree();
-    const double penalty = 2. * degree * (degree + 1) * dinfo.face->measure() /
-                           dinfo.cell->measure();
-
-    for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
-      for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
-        local_vector(i) +=
-          (-penalty * fe.shape_value(i, k)              // (-sigma * v_i(x_k)
-           + fe.normal_vector(k) * fe.shape_grad(i, k)) // + n * grad v_i(x_k))
-          * boundary_values[k] * fe.JxW(k);             // u^D(x_k) * dx
-  }
+    template <int dim>
+    void cell(MeshWorker::DoFInfo<dim> &, MeshWorker::IntegrationInfo<dim> &)
+    {}
 
 
-  template <int dim>
-  void RHSIntegrator<dim>::face(MeshWorker::DoFInfo<dim> &,
-                                MeshWorker::DoFInfo<dim> &,
-                                MeshWorker::IntegrationInfo<dim> &,
-                                MeshWorker::IntegrationInfo<dim> &) const
-  {}
+    template <int dim>
+    void boundary(MeshWorker::DoFInfo<dim>         &dinfo,
+                  MeshWorker::IntegrationInfo<dim> &info)
+    {
+      const FEValuesBase<dim> &fe           = info.fe_values();
+      Vector<double>          &local_vector = dinfo.vector(0).block(0);
+
+      std::vector<double> boundary_values(fe.n_quadrature_points);
+      exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
+
+      const unsigned int degree  = fe.get_fe().tensor_degree();
+      const double       penalty = 2. * degree * (degree + 1) *
+                             dinfo.face->measure() / dinfo.cell->measure();
+
+      for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
+        for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+          local_vector(i) +=
+            (-penalty * fe.shape_value(i, k) // (-sigma * v_i(x_k)
+             +
+             fe.normal_vector(k) * fe.shape_grad(i, k)) // + n * grad v_i(x_k))
+            * boundary_values[k] * fe.JxW(k);           // u^D(x_k) * dx
+    }
 
 
+    template <int dim>
+    void face(MeshWorker::DoFInfo<dim> &,
+              MeshWorker::DoFInfo<dim> &,
+              MeshWorker::IntegrationInfo<dim> &,
+              MeshWorker::IntegrationInfo<dim> &)
+    {}
+  } // namespace RHSIntegrator
+
   // The third local integrator is responsible for the contributions to the
   // error estimate. This is the standard energy estimator due to Karakashian
   // and Pascal (2003).
-  template <int dim>
-  class Estimator : public MeshWorker::LocalIntegrator<dim>
-  {
-  public:
-    void cell(MeshWorker::DoFInfo<dim>         &dinfo,
-              MeshWorker::IntegrationInfo<dim> &info) const override;
-    void boundary(MeshWorker::DoFInfo<dim>         &dinfo,
-                  MeshWorker::IntegrationInfo<dim> &info) const override;
-    void face(MeshWorker::DoFInfo<dim>         &dinfo1,
-              MeshWorker::DoFInfo<dim>         &dinfo2,
-              MeshWorker::IntegrationInfo<dim> &info1,
-              MeshWorker::IntegrationInfo<dim> &info2) const override;
-  };
-
-
   // The cell contribution is the Laplacian of the discrete solution, since
   // the right hand side is zero.
-  template <int dim>
-  void Estimator<dim>::cell(MeshWorker::DoFInfo<dim>         &dinfo,
-                            MeshWorker::IntegrationInfo<dim> &info) const
+  namespace Estimator
   {
-    const FEValuesBase<dim> &fe = info.fe_values();
-
-    const std::vector<Tensor<2, dim>> &DDuh = info.hessians[0][0];
-    for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
-      {
-        const double t = dinfo.cell->diameter() * trace(DDuh[k]);
-        dinfo.value(0) += t * t * fe.JxW(k);
-      }
-    dinfo.value(0) = std::sqrt(dinfo.value(0));
-  }
+    template <int dim>
+    void cell(MeshWorker::DoFInfo<dim>         &dinfo,
+              MeshWorker::IntegrationInfo<dim> &info)
+    {
+      const FEValuesBase<dim> &fe = info.fe_values();
+
+      const std::vector<Tensor<2, dim>> &DDuh = info.hessians[0][0];
+      for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
+        {
+          const double t = dinfo.cell->diameter() * trace(DDuh[k]);
+          dinfo.value(0) += t * t * fe.JxW(k);
+        }
+      dinfo.value(0) = std::sqrt(dinfo.value(0));
+    }
 
-  // At the boundary, we use simply a weighted form of the boundary residual,
-  // namely the norm of the difference between the finite element solution and
-  // the correct boundary condition.
-  template <int dim>
-  void Estimator<dim>::boundary(MeshWorker::DoFInfo<dim>         &dinfo,
-                                MeshWorker::IntegrationInfo<dim> &info) const
-  {
-    const FEValuesBase<dim> &fe = info.fe_values();
+    // At the boundary, we use simply a weighted form of the boundary residual,
+    // namely the norm of the difference between the finite element solution and
+    // the correct boundary condition.
+    template <int dim>
+    void boundary(MeshWorker::DoFInfo<dim>         &dinfo,
+                  MeshWorker::IntegrationInfo<dim> &info)
+    {
+      const FEValuesBase<dim> &fe = info.fe_values();
 
-    std::vector<double> boundary_values(fe.n_quadrature_points);
-    exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
+      std::vector<double> boundary_values(fe.n_quadrature_points);
+      exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
 
-    const std::vector<double> &uh = info.values[0][0];
+      const std::vector<double> &uh = info.values[0][0];
 
-    const unsigned int degree = fe.get_fe().tensor_degree();
-    const double penalty = 2. * degree * (degree + 1) * dinfo.face->measure() /
-                           dinfo.cell->measure();
+      const unsigned int degree  = fe.get_fe().tensor_degree();
+      const double       penalty = 2. * degree * (degree + 1) *
+                             dinfo.face->measure() / dinfo.cell->measure();
 
-    for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
-      {
-        const double diff = boundary_values[k] - uh[k];
-        dinfo.value(0) += penalty * diff * diff * fe.JxW(k);
-      }
-    dinfo.value(0) = std::sqrt(dinfo.value(0));
-  }
+      for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
+        {
+          const double diff = boundary_values[k] - uh[k];
+          dinfo.value(0) += penalty * diff * diff * fe.JxW(k);
+        }
+      dinfo.value(0) = std::sqrt(dinfo.value(0));
+    }
 
 
-  // Finally, on interior faces, the estimator consists of the jumps of the
-  // solution and its normal derivative, weighted appropriately.
-  template <int dim>
-  void Estimator<dim>::face(MeshWorker::DoFInfo<dim>         &dinfo1,
-                            MeshWorker::DoFInfo<dim>         &dinfo2,
-                            MeshWorker::IntegrationInfo<dim> &info1,
-                            MeshWorker::IntegrationInfo<dim> &info2) const
-  {
-    const FEValuesBase<dim>           &fe   = info1.fe_values();
-    const std::vector<double>         &uh1  = info1.values[0][0];
-    const std::vector<double>         &uh2  = info2.values[0][0];
-    const std::vector<Tensor<1, dim>> &Duh1 = info1.gradients[0][0];
-    const std::vector<Tensor<1, dim>> &Duh2 = info2.gradients[0][0];
-
-    const unsigned int degree = fe.get_fe().tensor_degree();
-    const double       penalty1 =
-      degree * (degree + 1) * dinfo1.face->measure() / dinfo1.cell->measure();
-    const double penalty2 =
-      degree * (degree + 1) * dinfo2.face->measure() / dinfo2.cell->measure();
-    const double penalty = penalty1 + penalty2;
-    const double h       = dinfo1.face->measure();
-
-    for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
-      {
-        const double diff1 = uh1[k] - uh2[k];
-        const double diff2 =
-          fe.normal_vector(k) * Duh1[k] - fe.normal_vector(k) * Duh2[k];
-        dinfo1.value(0) +=
-          (penalty * diff1 * diff1 + h * diff2 * diff2) * fe.JxW(k);
-      }
-    dinfo1.value(0) = std::sqrt(dinfo1.value(0));
-    dinfo2.value(0) = dinfo1.value(0);
-  }
+    // Finally, on interior faces, the estimator consists of the jumps of the
+    // solution and its normal derivative, weighted appropriately.
+    template <int dim>
+    void face(MeshWorker::DoFInfo<dim>         &dinfo1,
+              MeshWorker::DoFInfo<dim>         &dinfo2,
+              MeshWorker::IntegrationInfo<dim> &info1,
+              MeshWorker::IntegrationInfo<dim> &info2)
+    {
+      const FEValuesBase<dim>           &fe   = info1.fe_values();
+      const std::vector<double>         &uh1  = info1.values[0][0];
+      const std::vector<double>         &uh2  = info2.values[0][0];
+      const std::vector<Tensor<1, dim>> &Duh1 = info1.gradients[0][0];
+      const std::vector<Tensor<1, dim>> &Duh2 = info2.gradients[0][0];
+
+      const unsigned int degree = fe.get_fe().tensor_degree();
+      const double       penalty1 =
+        degree * (degree + 1) * dinfo1.face->measure() / dinfo1.cell->measure();
+      const double penalty2 =
+        degree * (degree + 1) * dinfo2.face->measure() / dinfo2.cell->measure();
+      const double penalty = penalty1 + penalty2;
+      const double h       = dinfo1.face->measure();
+
+      for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
+        {
+          const double diff1 = uh1[k] - uh2[k];
+          const double diff2 =
+            fe.normal_vector(k) * Duh1[k] - fe.normal_vector(k) * Duh2[k];
+          dinfo1.value(0) +=
+            (penalty * diff1 * diff1 + h * diff2 * diff2) * fe.JxW(k);
+        }
+      dinfo1.value(0) = std::sqrt(dinfo1.value(0));
+      dinfo2.value(0) = dinfo1.value(0);
+    }
+  } // namespace Estimator
 
   // Finally we have an integrator for the error. Since the energy norm for
   // discontinuous Galerkin problems not only involves the difference of the
@@ -445,22 +411,9 @@ namespace Step39
   // \sum_{K\in \mathbb T_h} \|\nabla u\|_K^2 + \sum_{F \in F_h^i}
   // 4\sigma_F\|\average{ u \mathbf n}\|^2_F + \sum_{F \in F_h^b}
   // 2\sigma_F\|u\|^2_F @f]
-
-  template <int dim>
-  class ErrorIntegrator : public MeshWorker::LocalIntegrator<dim>
-  {
-  public:
-    void cell(MeshWorker::DoFInfo<dim>         &dinfo,
-              MeshWorker::IntegrationInfo<dim> &info) const override;
-    void boundary(MeshWorker::DoFInfo<dim>         &dinfo,
-                  MeshWorker::IntegrationInfo<dim> &info) const override;
-    void face(MeshWorker::DoFInfo<dim>         &dinfo1,
-              MeshWorker::DoFInfo<dim>         &dinfo2,
-              MeshWorker::IntegrationInfo<dim> &info1,
-              MeshWorker::IntegrationInfo<dim> &info2) const override;
-  };
-
-  // Here we have the integration on cells. There is currently no good
+  //
+  // Below, the first function is, as always, the integration
+  // on cells. There is currently no good
   // interface in MeshWorker that would allow us to access values of regular
   // functions in the quadrature points. Thus, we have to create the vectors
   // for the exact function's values and gradients inside the cell
@@ -472,89 +425,91 @@ namespace Step39
   // and compute the <i>L<sup>2</sup></i>-error in the same loop. Obviously,
   // this one does not have any jump terms and only appears in the integration
   // on cells.
-  template <int dim>
-  void ErrorIntegrator<dim>::cell(MeshWorker::DoFInfo<dim>         &dinfo,
-                                  MeshWorker::IntegrationInfo<dim> &info) const
-  {
-    const FEValuesBase<dim>    &fe = info.fe_values();
-    std::vector<Tensor<1, dim>> exact_gradients(fe.n_quadrature_points);
-    std::vector<double>         exact_values(fe.n_quadrature_points);
-
-    exact_solution.gradient_list(fe.get_quadrature_points(), exact_gradients);
-    exact_solution.value_list(fe.get_quadrature_points(), exact_values);
-
-    const std::vector<Tensor<1, dim>> &Duh = info.gradients[0][0];
-    const std::vector<double>         &uh  = info.values[0][0];
 
-    for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
-      {
-        double sum = 0;
-        for (unsigned int d = 0; d < dim; ++d)
-          {
-            const double diff = exact_gradients[k][d] - Duh[k][d];
-            sum += diff * diff;
-          }
-        const double diff = exact_values[k] - uh[k];
-        dinfo.value(0) += sum * fe.JxW(k);
-        dinfo.value(1) += diff * diff * fe.JxW(k);
-      }
-    dinfo.value(0) = std::sqrt(dinfo.value(0));
-    dinfo.value(1) = std::sqrt(dinfo.value(1));
-  }
-
-
-  template <int dim>
-  void
-  ErrorIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim>         &dinfo,
-                                 MeshWorker::IntegrationInfo<dim> &info) const
+  namespace ErrorIntegrator
   {
-    const FEValuesBase<dim> &fe = info.fe_values();
+    template <int dim>
+    void cell(MeshWorker::DoFInfo<dim>         &dinfo,
+              MeshWorker::IntegrationInfo<dim> &info)
+    {
+      const FEValuesBase<dim>    &fe = info.fe_values();
+      std::vector<Tensor<1, dim>> exact_gradients(fe.n_quadrature_points);
+      std::vector<double>         exact_values(fe.n_quadrature_points);
+
+      exact_solution.gradient_list(fe.get_quadrature_points(), exact_gradients);
+      exact_solution.value_list(fe.get_quadrature_points(), exact_values);
+
+      const std::vector<Tensor<1, dim>> &Duh = info.gradients[0][0];
+      const std::vector<double>         &uh  = info.values[0][0];
+
+      for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
+        {
+          double sum = 0;
+          for (unsigned int d = 0; d < dim; ++d)
+            {
+              const double diff = exact_gradients[k][d] - Duh[k][d];
+              sum += diff * diff;
+            }
+          const double diff = exact_values[k] - uh[k];
+          dinfo.value(0) += sum * fe.JxW(k);
+          dinfo.value(1) += diff * diff * fe.JxW(k);
+        }
+      dinfo.value(0) = std::sqrt(dinfo.value(0));
+      dinfo.value(1) = std::sqrt(dinfo.value(1));
+    }
 
-    std::vector<double> exact_values(fe.n_quadrature_points);
-    exact_solution.value_list(fe.get_quadrature_points(), exact_values);
 
-    const std::vector<double> &uh = info.values[0][0];
+    template <int dim>
+    void boundary(MeshWorker::DoFInfo<dim>         &dinfo,
+                  MeshWorker::IntegrationInfo<dim> &info)
+    {
+      const FEValuesBase<dim> &fe = info.fe_values();
 
-    const unsigned int degree = fe.get_fe().tensor_degree();
-    const double penalty = 2. * degree * (degree + 1) * dinfo.face->measure() /
-                           dinfo.cell->measure();
+      std::vector<double> exact_values(fe.n_quadrature_points);
+      exact_solution.value_list(fe.get_quadrature_points(), exact_values);
 
-    for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
-      {
-        const double diff = exact_values[k] - uh[k];
-        dinfo.value(0) += penalty * diff * diff * fe.JxW(k);
-      }
-    dinfo.value(0) = std::sqrt(dinfo.value(0));
-  }
+      const std::vector<double> &uh = info.values[0][0];
 
+      const unsigned int degree  = fe.get_fe().tensor_degree();
+      const double       penalty = 2. * degree * (degree + 1) *
+                             dinfo.face->measure() / dinfo.cell->measure();
 
-  template <int dim>
-  void ErrorIntegrator<dim>::face(MeshWorker::DoFInfo<dim>         &dinfo1,
-                                  MeshWorker::DoFInfo<dim>         &dinfo2,
-                                  MeshWorker::IntegrationInfo<dim> &info1,
-                                  MeshWorker::IntegrationInfo<dim> &info2) const
-  {
-    const FEValuesBase<dim>   &fe  = info1.fe_values();
-    const std::vector<double> &uh1 = info1.values[0][0];
-    const std::vector<double> &uh2 = info2.values[0][0];
-
-    const unsigned int degree = fe.get_fe().tensor_degree();
-    const double       penalty1 =
-      degree * (degree + 1) * dinfo1.face->measure() / dinfo1.cell->measure();
-    const double penalty2 =
-      degree * (degree + 1) * dinfo2.face->measure() / dinfo2.cell->measure();
-    const double penalty = penalty1 + penalty2;
-
-    for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
-      {
-        const double diff = uh1[k] - uh2[k];
-        dinfo1.value(0) += (penalty * diff * diff) * fe.JxW(k);
-      }
-    dinfo1.value(0) = std::sqrt(dinfo1.value(0));
-    dinfo2.value(0) = dinfo1.value(0);
-  }
+      for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
+        {
+          const double diff = exact_values[k] - uh[k];
+          dinfo.value(0) += penalty * diff * diff * fe.JxW(k);
+        }
+      dinfo.value(0) = std::sqrt(dinfo.value(0));
+    }
 
 
+    template <int dim>
+    void face(MeshWorker::DoFInfo<dim>         &dinfo1,
+              MeshWorker::DoFInfo<dim>         &dinfo2,
+              MeshWorker::IntegrationInfo<dim> &info1,
+              MeshWorker::IntegrationInfo<dim> &info2)
+    {
+      const FEValuesBase<dim>   &fe  = info1.fe_values();
+      const std::vector<double> &uh1 = info1.values[0][0];
+      const std::vector<double> &uh2 = info2.values[0][0];
+
+      const unsigned int degree = fe.get_fe().tensor_degree();
+      const double       penalty1 =
+        degree * (degree + 1) * dinfo1.face->measure() / dinfo1.cell->measure();
+      const double penalty2 =
+        degree * (degree + 1) * dinfo2.face->measure() / dinfo2.cell->measure();
+      const double penalty = penalty1 + penalty2;
+
+      for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
+        {
+          const double diff = uh1[k] - uh2[k];
+          dinfo1.value(0) += (penalty * diff * diff) * fe.JxW(k);
+        }
+      dinfo1.value(0) = std::sqrt(dinfo1.value(0));
+      dinfo2.value(0) = dinfo1.value(0);
+    }
+  } // namespace ErrorIntegrator
+
 
   // @sect3{The main class}
 
@@ -724,7 +679,7 @@ namespace Step39
     info_box.initialize(fe, mapping);
 
     // This is the object into which we integrate local data. It is filled by
-    // the local integration routines in MatrixIntegrator and then used by the
+    // the local integration routines in `MatrixIntegrator` and then used by the
     // assembler to distribute the information into the global matrix.
     MeshWorker::DoFInfo<dim> dof_info(dof_handler);
 
@@ -735,20 +690,19 @@ namespace Step39
     MeshWorker::Assembler::MatrixSimple<SparseMatrix<double>> assembler;
     assembler.initialize(matrix);
 
-    // Now comes the part we coded ourselves, the local
-    // integrator. This is the only part which is problem dependent.
-    MatrixIntegrator<dim> integrator;
-    // Now, we throw everything into a MeshWorker::loop(), which here
+    // Now, we throw everything into a MeshWorker::loop<dim, dim>(), which here
     // traverses all active cells of the mesh, computes cell and face matrices
     // and assembles them into the global matrix. We use the variable
     // <tt>dof_handler</tt> here in order to use the global numbering of
     // degrees of freedom.
-    MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
-                                           dof_handler.end(),
-                                           dof_info,
-                                           info_box,
-                                           integrator,
-                                           assembler);
+    MeshWorker::loop<dim, dim>(dof_handler.begin_active(),
+                               dof_handler.end(),
+                               dof_info,
+                               info_box,
+                               &MatrixIntegrator::cell<dim>,
+                               &MatrixIntegrator::boundary<dim>,
+                               &MatrixIntegrator::face<dim>,
+                               assembler);
   }
 
 
@@ -771,17 +725,18 @@ namespace Step39
     assembler.initialize(mg_matrix);
     assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down);
 
-    MatrixIntegrator<dim> integrator;
     // Here is the other difference to the previous function: we run
     // over all cells, not only the active ones. And we use functions
     // ending on <code>_mg</code> since we need the degrees of freedom
     // on each level, not the global numbering.
-    MeshWorker::integration_loop<dim, dim>(dof_handler.begin_mg(),
-                                           dof_handler.end_mg(),
-                                           dof_info,
-                                           info_box,
-                                           integrator,
-                                           assembler);
+    MeshWorker::loop<dim, dim>(dof_handler.begin_mg(),
+                               dof_handler.end_mg(),
+                               dof_info,
+                               info_box,
+                               &MatrixIntegrator::cell<dim>,
+                               &MatrixIntegrator::boundary<dim>,
+                               &MatrixIntegrator::face<dim>,
+                               assembler);
   }
 
 
@@ -808,13 +763,14 @@ namespace Step39
     data.add<Vector<double> *>(&right_hand_side, "RHS");
     assembler.initialize(data);
 
-    RHSIntegrator<dim> integrator;
-    MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
-                                           dof_handler.end(),
-                                           dof_info,
-                                           info_box,
-                                           integrator,
-                                           assembler);
+    MeshWorker::loop<dim, dim>(dof_handler.begin_active(),
+                               dof_handler.end(),
+                               dof_info,
+                               info_box,
+                               &RHSIntegrator::cell<dim>,
+                               &RHSIntegrator::boundary<dim>,
+                               &RHSIntegrator::face<dim>,
+                               assembler);
 
     right_hand_side *= -1.;
   }
@@ -946,13 +902,14 @@ namespace Step39
     out_data.add<BlockVector<double> *>(&estimates, "cells");
     assembler.initialize(out_data, false);
 
-    Estimator<dim> integrator;
-    MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
-                                           dof_handler.end(),
-                                           dof_info,
-                                           info_box,
-                                           integrator,
-                                           assembler);
+    MeshWorker::loop<dim, dim>(dof_handler.begin_active(),
+                               dof_handler.end(),
+                               dof_info,
+                               info_box,
+                               &Estimator::cell<dim>,
+                               &Estimator::boundary<dim>,
+                               &Estimator::face<dim>,
+                               assembler);
 
     // Right before we return the result of the error estimate, we restore the
     // old user indices.
@@ -1006,13 +963,14 @@ namespace Step39
     out_data.add<BlockVector<double> *>(&errors, "cells");
     assembler.initialize(out_data, false);
 
-    ErrorIntegrator<dim> integrator;
-    MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
-                                           dof_handler.end(),
-                                           dof_info,
-                                           info_box,
-                                           integrator,
-                                           assembler);
+    MeshWorker::loop<dim, dim>(dof_handler.begin_active(),
+                               dof_handler.end(),
+                               dof_info,
+                               info_box,
+                               &ErrorIntegrator::cell<dim>,
+                               &ErrorIntegrator::boundary<dim>,
+                               &ErrorIntegrator::face<dim>,
+                               assembler);
     triangulation.load_user_indices(old_user_indices);
 
     deallog << "energy-error: " << errors.block(0).l2_norm() << std::endl;

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.