From 8071f1141a05b42ca95c8fc924283e0e5e3d1d5c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 20 May 2024 20:41:42 -0600 Subject: [PATCH] Avoid the use of MeshWorker::LocalIntegrator and MeshWorker::integration_loop(). --- examples/step-39/step-39.cc | 838 +++++++++++++++++------------------- 1 file changed, 398 insertions(+), 440 deletions(-) diff --git a/examples/step-39/step-39.cc b/examples/step-39/step-39.cc index 105f692e55..b423668346 100644 --- a/examples/step-39/step-39.cc +++ b/examples/step-39/step-39.cc @@ -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. 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 - class MatrixIntegrator : public MeshWorker::LocalIntegrator - { - public: - void cell(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const override; - void boundary(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const override; - void face(MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - MeshWorker::IntegrationInfo &info1, - MeshWorker::IntegrationInfo &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 and MeshWorker::IntegrationInfo. 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 - double ip_penalty_factor(const MeshWorker::DoFInfo &dinfo1, - const MeshWorker::DoFInfo &dinfo2, - unsigned int deg1, - unsigned int deg2) - { - const unsigned int normal1 = - GeometryInfo::unit_normal_direction[dinfo1.face_number]; - const unsigned int normal2 = - GeometryInfo::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 - void MatrixIntegrator::cell(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const - { - FullMatrix &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 - void - MatrixIntegrator::boundary(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const + namespace MatrixIntegrator { - const FEValuesBase &fe_face_values = info.fe_values(0); - - FullMatrix &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 + double ip_penalty_factor(const MeshWorker::DoFInfo &dinfo1, + const MeshWorker::DoFInfo &dinfo2, + unsigned int deg1, + unsigned int deg2) + { + const unsigned int normal1 = + GeometryInfo::unit_normal_direction[dinfo1.face_number]; + const unsigned int normal2 = + GeometryInfo::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 - void - MatrixIntegrator::face(MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - MeshWorker::IntegrationInfo &info1, - MeshWorker::IntegrationInfo &info2) const - { - const FEValuesBase &fe_face_values_1 = info1.fe_values(0); - const FEValuesBase &fe_face_values_2 = info2.fe_values(0); - - FullMatrix &M11 = dinfo1.matrix(0, false).matrix; - FullMatrix &M12 = dinfo1.matrix(0, true).matrix; - FullMatrix &M21 = dinfo2.matrix(0, true).matrix; - FullMatrix &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 + void cell(MeshWorker::DoFInfo &dinfo, + MeshWorker::IntegrationInfo &info) + { + FullMatrix &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 - class RHSIntegrator : public MeshWorker::LocalIntegrator - { - public: - void cell(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const override; + // Boundary faces use the Nitsche method to impose boundary values: + template void boundary(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const override; + MeshWorker::IntegrationInfo &info) + { + const FEValuesBase &fe_face_values = info.fe_values(0); + + FullMatrix &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 void face(MeshWorker::DoFInfo &dinfo1, MeshWorker::DoFInfo &dinfo2, MeshWorker::IntegrationInfo &info1, - MeshWorker::IntegrationInfo &info2) const override; - }; - - - template - void RHSIntegrator::cell(MeshWorker::DoFInfo &, - MeshWorker::IntegrationInfo &) const - {} - + MeshWorker::IntegrationInfo &info2) + { + const FEValuesBase &fe_face_values_1 = info1.fe_values(0); + const FEValuesBase &fe_face_values_2 = info2.fe_values(0); + + FullMatrix &M11 = dinfo1.matrix(0, false).matrix; + FullMatrix &M12 = dinfo1.matrix(0, true).matrix; + FullMatrix &M21 = dinfo2.matrix(0, true).matrix; + FullMatrix &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 - void - RHSIntegrator::boundary(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &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 &fe = info.fe_values(); - Vector &local_vector = dinfo.vector(0).block(0); - - std::vector 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 + void cell(MeshWorker::DoFInfo &, MeshWorker::IntegrationInfo &) + {} - template - void RHSIntegrator::face(MeshWorker::DoFInfo &, - MeshWorker::DoFInfo &, - MeshWorker::IntegrationInfo &, - MeshWorker::IntegrationInfo &) const - {} + template + void boundary(MeshWorker::DoFInfo &dinfo, + MeshWorker::IntegrationInfo &info) + { + const FEValuesBase &fe = info.fe_values(); + Vector &local_vector = dinfo.vector(0).block(0); + + std::vector 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 + void face(MeshWorker::DoFInfo &, + MeshWorker::DoFInfo &, + MeshWorker::IntegrationInfo &, + MeshWorker::IntegrationInfo &) + {} + } // 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 - class Estimator : public MeshWorker::LocalIntegrator - { - public: - void cell(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const override; - void boundary(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const override; - void face(MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - MeshWorker::IntegrationInfo &info1, - MeshWorker::IntegrationInfo &info2) const override; - }; - - // The cell contribution is the Laplacian of the discrete solution, since // the right hand side is zero. - template - void Estimator::cell(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const + namespace Estimator { - const FEValuesBase &fe = info.fe_values(); - - const std::vector> &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 + void cell(MeshWorker::DoFInfo &dinfo, + MeshWorker::IntegrationInfo &info) + { + const FEValuesBase &fe = info.fe_values(); + + const std::vector> &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 - void Estimator::boundary(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const - { - const FEValuesBase &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 + void boundary(MeshWorker::DoFInfo &dinfo, + MeshWorker::IntegrationInfo &info) + { + const FEValuesBase &fe = info.fe_values(); - std::vector boundary_values(fe.n_quadrature_points); - exact_solution.value_list(fe.get_quadrature_points(), boundary_values); + std::vector boundary_values(fe.n_quadrature_points); + exact_solution.value_list(fe.get_quadrature_points(), boundary_values); - const std::vector &uh = info.values[0][0]; + const std::vector &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 - void Estimator::face(MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - MeshWorker::IntegrationInfo &info1, - MeshWorker::IntegrationInfo &info2) const - { - const FEValuesBase &fe = info1.fe_values(); - const std::vector &uh1 = info1.values[0][0]; - const std::vector &uh2 = info2.values[0][0]; - const std::vector> &Duh1 = info1.gradients[0][0]; - const std::vector> &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 + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + MeshWorker::IntegrationInfo &info1, + MeshWorker::IntegrationInfo &info2) + { + const FEValuesBase &fe = info1.fe_values(); + const std::vector &uh1 = info1.values[0][0]; + const std::vector &uh2 = info2.values[0][0]; + const std::vector> &Duh1 = info1.gradients[0][0]; + const std::vector> &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 - class ErrorIntegrator : public MeshWorker::LocalIntegrator - { - public: - void cell(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const override; - void boundary(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const override; - void face(MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - MeshWorker::IntegrationInfo &info1, - MeshWorker::IntegrationInfo &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 L2-error in the same loop. Obviously, // this one does not have any jump terms and only appears in the integration // on cells. - template - void ErrorIntegrator::cell(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const - { - const FEValuesBase &fe = info.fe_values(); - std::vector> exact_gradients(fe.n_quadrature_points); - std::vector 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> &Duh = info.gradients[0][0]; - const std::vector &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 - void - ErrorIntegrator::boundary(MeshWorker::DoFInfo &dinfo, - MeshWorker::IntegrationInfo &info) const + namespace ErrorIntegrator { - const FEValuesBase &fe = info.fe_values(); + template + void cell(MeshWorker::DoFInfo &dinfo, + MeshWorker::IntegrationInfo &info) + { + const FEValuesBase &fe = info.fe_values(); + std::vector> exact_gradients(fe.n_quadrature_points); + std::vector 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> &Duh = info.gradients[0][0]; + const std::vector &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 exact_values(fe.n_quadrature_points); - exact_solution.value_list(fe.get_quadrature_points(), exact_values); - const std::vector &uh = info.values[0][0]; + template + void boundary(MeshWorker::DoFInfo &dinfo, + MeshWorker::IntegrationInfo &info) + { + const FEValuesBase &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 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 &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 - void ErrorIntegrator::face(MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - MeshWorker::IntegrationInfo &info1, - MeshWorker::IntegrationInfo &info2) const - { - const FEValuesBase &fe = info1.fe_values(); - const std::vector &uh1 = info1.values[0][0]; - const std::vector &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 + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + MeshWorker::IntegrationInfo &info1, + MeshWorker::IntegrationInfo &info2) + { + const FEValuesBase &fe = info1.fe_values(); + const std::vector &uh1 = info1.values[0][0]; + const std::vector &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 dof_info(dof_handler); @@ -735,20 +690,19 @@ namespace Step39 MeshWorker::Assembler::MatrixSimple> assembler; assembler.initialize(matrix); - // Now comes the part we coded ourselves, the local - // integrator. This is the only part which is problem dependent. - MatrixIntegrator integrator; - // Now, we throw everything into a MeshWorker::loop(), which here + // Now, we throw everything into a MeshWorker::loop(), 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 // dof_handler here in order to use the global numbering of // degrees of freedom. - MeshWorker::integration_loop(dof_handler.begin_active(), - dof_handler.end(), - dof_info, - info_box, - integrator, - assembler); + MeshWorker::loop(dof_handler.begin_active(), + dof_handler.end(), + dof_info, + info_box, + &MatrixIntegrator::cell, + &MatrixIntegrator::boundary, + &MatrixIntegrator::face, + assembler); } @@ -771,17 +725,18 @@ namespace Step39 assembler.initialize(mg_matrix); assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down); - MatrixIntegrator 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 _mg since we need the degrees of freedom // on each level, not the global numbering. - MeshWorker::integration_loop(dof_handler.begin_mg(), - dof_handler.end_mg(), - dof_info, - info_box, - integrator, - assembler); + MeshWorker::loop(dof_handler.begin_mg(), + dof_handler.end_mg(), + dof_info, + info_box, + &MatrixIntegrator::cell, + &MatrixIntegrator::boundary, + &MatrixIntegrator::face, + assembler); } @@ -808,13 +763,14 @@ namespace Step39 data.add *>(&right_hand_side, "RHS"); assembler.initialize(data); - RHSIntegrator integrator; - MeshWorker::integration_loop(dof_handler.begin_active(), - dof_handler.end(), - dof_info, - info_box, - integrator, - assembler); + MeshWorker::loop(dof_handler.begin_active(), + dof_handler.end(), + dof_info, + info_box, + &RHSIntegrator::cell, + &RHSIntegrator::boundary, + &RHSIntegrator::face, + assembler); right_hand_side *= -1.; } @@ -946,13 +902,14 @@ namespace Step39 out_data.add *>(&estimates, "cells"); assembler.initialize(out_data, false); - Estimator integrator; - MeshWorker::integration_loop(dof_handler.begin_active(), - dof_handler.end(), - dof_info, - info_box, - integrator, - assembler); + MeshWorker::loop(dof_handler.begin_active(), + dof_handler.end(), + dof_info, + info_box, + &Estimator::cell, + &Estimator::boundary, + &Estimator::face, + 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 *>(&errors, "cells"); assembler.initialize(out_data, false); - ErrorIntegrator integrator; - MeshWorker::integration_loop(dof_handler.begin_active(), - dof_handler.end(), - dof_info, - info_box, - integrator, - assembler); + MeshWorker::loop(dof_handler.begin_active(), + dof_handler.end(), + dof_info, + info_box, + &ErrorIntegrator::cell, + &ErrorIntegrator::boundary, + &ErrorIntegrator::face, + assembler); triangulation.load_user_indices(old_user_indices); deallog << "energy-error: " << errors.block(0).l2_norm() << std::endl; -- 2.39.5