From a10c3c006719989ef90104445725c2abe073bb07 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 20 Jun 2024 16:46:12 -0600 Subject: [PATCH] Avoid the use of 'auto' in step-74. --- examples/step-74/step-74.cc | 464 +++++++++++++++++++----------------- 1 file changed, 239 insertions(+), 225 deletions(-) diff --git a/examples/step-74/step-74.cc b/examples/step-74/step-74.cc index 98873dcb22..e8584faa31 100644 --- a/examples/step-74/step-74.cc +++ b/examples/step-74/step-74.cc @@ -343,14 +343,17 @@ namespace Step74 void SIPGLaplace::assemble_system() { const auto cell_worker = - [&](const auto &cell, auto &scratch_data, auto ©_data) { + [&](const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch_data, + CopyData ©_data) { const FEValues &fe_v = scratch_data.reinit(cell); const unsigned int dofs_per_cell = fe_v.dofs_per_cell; copy_data.reinit(cell, dofs_per_cell); - const auto &q_points = scratch_data.get_quadrature_points(); - const unsigned int n_q_points = q_points.size(); - const std::vector &JxW = scratch_data.get_JxW_values(); + const std::vector> &q_points = + scratch_data.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + const std::vector &JxW = scratch_data.get_JxW_values(); std::vector rhs(n_q_points); rhs_function->value_list(q_points, rhs); @@ -372,117 +375,120 @@ namespace Step74 }; // Next, we need a function that assembles face integrals on the boundary: - const auto boundary_worker = [&](const auto &cell, - const unsigned int &face_no, - auto &scratch_data, - auto ©_data) { - const FEFaceValuesBase &fe_fv = scratch_data.reinit(cell, face_no); + const auto boundary_worker = + [&](const typename DoFHandler::active_cell_iterator &cell, + const unsigned int &face_no, + ScratchData &scratch_data, + CopyData ©_data) { + const FEFaceValuesBase &fe_fv = scratch_data.reinit(cell, face_no); - const auto &q_points = scratch_data.get_quadrature_points(); - const unsigned int n_q_points = q_points.size(); - const unsigned int dofs_per_cell = fe_fv.dofs_per_cell; + const std::vector> &q_points = + scratch_data.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + const unsigned int dofs_per_cell = fe_fv.dofs_per_cell; - const std::vector &JxW = scratch_data.get_JxW_values(); - const std::vector> &normals = - scratch_data.get_normal_vectors(); + const std::vector &JxW = scratch_data.get_JxW_values(); + const std::vector> &normals = + scratch_data.get_normal_vectors(); - std::vector g(n_q_points); - exact_solution->value_list(q_points, g); + std::vector g(n_q_points); + exact_solution->value_list(q_points, g); - const double extent1 = cell->measure() / cell->face(face_no)->measure(); - const double penalty = get_penalty_factor(degree, extent1, extent1); + const double extent1 = cell->measure() / cell->face(face_no)->measure(); + const double penalty = get_penalty_factor(degree, extent1, extent1); - for (unsigned int point = 0; point < n_q_points; ++point) - { - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int j = 0; j < dofs_per_cell; ++j) - copy_data.cell_matrix(i, j) += - (-diffusion_coefficient * // - nu - fe_fv.shape_value(i, point) * // v_h - (fe_fv.shape_grad(j, point) * // (grad u_h . - normals[point]) // n) + for (unsigned int point = 0; point < n_q_points; ++point) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + copy_data.cell_matrix(i, j) += + (-diffusion_coefficient * // - nu + fe_fv.shape_value(i, point) * // v_h + (fe_fv.shape_grad(j, point) * // (grad u_h . + normals[point]) // n) - - diffusion_coefficient * // - nu - (fe_fv.shape_grad(i, point) * // (grad v_h . - normals[point]) * // n) - fe_fv.shape_value(j, point) // u_h + - diffusion_coefficient * // - nu + (fe_fv.shape_grad(i, point) * // (grad v_h . + normals[point]) * // n) + fe_fv.shape_value(j, point) // u_h - + diffusion_coefficient * penalty * // + nu sigma - fe_fv.shape_value(i, point) * // v_h - fe_fv.shape_value(j, point) // u_h + + diffusion_coefficient * penalty * // + nu sigma + fe_fv.shape_value(i, point) * // v_h + fe_fv.shape_value(j, point) // u_h - ) * - JxW[point]; // dx + ) * + JxW[point]; // dx - for (unsigned int i = 0; i < dofs_per_cell; ++i) - copy_data.cell_rhs(i) += - (-diffusion_coefficient * // - nu - (fe_fv.shape_grad(i, point) * // (grad v_h . - normals[point]) * // n) - g[point] // g + for (unsigned int i = 0; i < dofs_per_cell; ++i) + copy_data.cell_rhs(i) += + (-diffusion_coefficient * // - nu + (fe_fv.shape_grad(i, point) * // (grad v_h . + normals[point]) * // n) + g[point] // g - + diffusion_coefficient * penalty * // + nu sigma - fe_fv.shape_value(i, point) * g[point] // v_h g + + diffusion_coefficient * penalty * // + nu sigma + fe_fv.shape_value(i, point) * g[point] // v_h g - ) * - JxW[point]; // dx - } - }; + ) * + JxW[point]; // dx + } + }; // Finally, a function that assembles face integrals on interior // faces. To reinitialize FEInterfaceValues, we need to pass // cells, face and subface indices (for adaptive refinement) to // the reinit() function of FEInterfaceValues: - const auto face_worker = [&](const auto &cell, - const unsigned int &f, - const unsigned int &sf, - const auto &ncell, - const unsigned int &nf, - const unsigned int &nsf, - auto &scratch_data, - auto ©_data) { - const FEInterfaceValues &fe_iv = - scratch_data.reinit(cell, f, sf, ncell, nf, nsf); - - copy_data.face_data.emplace_back(); - CopyDataFace ©_data_face = copy_data.face_data.back(); - const unsigned int n_dofs_face = fe_iv.n_current_interface_dofs(); - copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices(); - copy_data_face.cell_matrix.reinit(n_dofs_face, n_dofs_face); - - const std::vector &JxW = fe_iv.get_JxW_values(); - const std::vector> &normals = fe_iv.get_normal_vectors(); - - const double extent1 = cell->measure() / cell->face(f)->measure(); - const double extent2 = ncell->measure() / ncell->face(nf)->measure(); - const double penalty = get_penalty_factor(degree, extent1, extent2); - - for (const unsigned int point : fe_iv.quadrature_point_indices()) - { - for (const unsigned int i : fe_iv.dof_indices()) - for (const unsigned int j : fe_iv.dof_indices()) - copy_data_face.cell_matrix(i, j) += - (-diffusion_coefficient * // - nu - fe_iv.jump_in_shape_values(i, point) * // [v_h] - (fe_iv.average_of_shape_gradients(j, // - point) * // ({grad u_h} . - normals[point]) // n) - - - - diffusion_coefficient * // - nu - (fe_iv.average_of_shape_gradients(i, point) * // (grad v_h . - normals[point]) * // n) - fe_iv.jump_in_shape_values(j, point) // [u_h] - - + diffusion_coefficient * penalty * // + nu sigma - fe_iv.jump_in_shape_values(i, point) * // [v_h] - fe_iv.jump_in_shape_values(j, point) // [u_h] - - ) * - JxW[point]; // dx - } - }; + const auto face_worker = + [&](const typename DoFHandler::cell_iterator &cell, + const unsigned int &f, + const unsigned int &sf, + const typename DoFHandler::cell_iterator &ncell, + const unsigned int &nf, + const unsigned int &nsf, + ScratchData &scratch_data, + CopyData ©_data) { + const FEInterfaceValues &fe_iv = + scratch_data.reinit(cell, f, sf, ncell, nf, nsf); + + copy_data.face_data.emplace_back(); + CopyDataFace ©_data_face = copy_data.face_data.back(); + const unsigned int n_dofs_face = fe_iv.n_current_interface_dofs(); + copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices(); + copy_data_face.cell_matrix.reinit(n_dofs_face, n_dofs_face); + + const std::vector &JxW = fe_iv.get_JxW_values(); + const std::vector> &normals = fe_iv.get_normal_vectors(); + + const double extent1 = cell->measure() / cell->face(f)->measure(); + const double extent2 = ncell->measure() / ncell->face(nf)->measure(); + const double penalty = get_penalty_factor(degree, extent1, extent2); + + for (const unsigned int point : fe_iv.quadrature_point_indices()) + { + for (const unsigned int i : fe_iv.dof_indices()) + for (const unsigned int j : fe_iv.dof_indices()) + copy_data_face.cell_matrix(i, j) += + (-diffusion_coefficient * // - nu + fe_iv.jump_in_shape_values(i, point) * // [v_h] + (fe_iv.average_of_shape_gradients(j, // + point) * // ({grad u_h} . + normals[point]) // n) + + - diffusion_coefficient * // - nu + (fe_iv.average_of_shape_gradients(i, + point) * // (grad v_h . + normals[point]) * // n) + fe_iv.jump_in_shape_values(j, point) // [u_h] + + + diffusion_coefficient * penalty * // + nu sigma + fe_iv.jump_in_shape_values(i, point) * // [v_h] + fe_iv.jump_in_shape_values(j, point) // [u_h] + + ) * + JxW[point]; // dx + } + }; // The following lambda function will then copy data into the // global matrix and right-hand side. Though there are no hanging @@ -491,7 +497,7 @@ namespace Step74 // AffineConstraints::distribute_local_to_global() functionality. AffineConstraints constraints; constraints.close(); - const auto copier = [&](const auto &c) { + const auto copier = [&](const CopyData &c) { constraints.distribute_local_to_global(c.cell_matrix, c.cell_rhs, c.local_dof_indices, @@ -499,7 +505,7 @@ namespace Step74 system_rhs); // Copy data from interior face assembly to the global matrix. - for (auto &cdf : c.face_data) + for (const CopyDataFace &cdf : c.face_data) { constraints.distribute_local_to_global(cdf.cell_matrix, cdf.joint_dof_indices, @@ -577,14 +583,16 @@ namespace Step74 void SIPGLaplace::compute_error_estimate() { const auto cell_worker = - [&](const auto &cell, auto &scratch_data, auto ©_data) { + [&](const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch_data, + CopyData ©_data) { const FEValues &fe_v = scratch_data.reinit(cell); copy_data.cell_index = cell->active_cell_index(); - const auto &q_points = fe_v.get_quadrature_points(); - const unsigned int n_q_points = q_points.size(); - const std::vector &JxW = fe_v.get_JxW_values(); + const std::vector> &q_points = fe_v.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + const std::vector &JxW = fe_v.get_JxW_values(); std::vector> hessians(n_q_points); fe_v.get_function_hessians(solution, hessians); @@ -606,95 +614,97 @@ namespace Step74 // Next compute boundary terms $\sum_{f\in \partial K \cap \partial \Omega} // \sigma \left\| [ u_h-g_D ] \right\|_f^2 $: - const auto boundary_worker = [&](const auto &cell, - const unsigned int &face_no, - auto &scratch_data, - auto ©_data) { - const FEFaceValuesBase &fe_fv = scratch_data.reinit(cell, face_no); + const auto boundary_worker = + [&](const typename DoFHandler::active_cell_iterator &cell, + const unsigned int &face_no, + ScratchData &scratch_data, + CopyData ©_data) { + const FEFaceValuesBase &fe_fv = scratch_data.reinit(cell, face_no); - const auto &q_points = fe_fv.get_quadrature_points(); - const unsigned n_q_points = q_points.size(); + const std::vector> &q_points = fe_fv.get_quadrature_points(); + const unsigned n_q_points = q_points.size(); - const std::vector &JxW = fe_fv.get_JxW_values(); + const std::vector &JxW = fe_fv.get_JxW_values(); - std::vector g(n_q_points); - exact_solution->value_list(q_points, g); + std::vector g(n_q_points); + exact_solution->value_list(q_points, g); - std::vector sol_u(n_q_points); - fe_fv.get_function_values(solution, sol_u); + std::vector sol_u(n_q_points); + fe_fv.get_function_values(solution, sol_u); - const double extent1 = cell->measure() / cell->face(face_no)->measure(); - const double penalty = get_penalty_factor(degree, extent1, extent1); + const double extent1 = cell->measure() / cell->face(face_no)->measure(); + const double penalty = get_penalty_factor(degree, extent1, extent1); - double difference_norm_square = 0.; - for (unsigned int point = 0; point < q_points.size(); ++point) - { - const double diff = (g[point] - sol_u[point]); - difference_norm_square += diff * diff * JxW[point]; - } - copy_data.value += penalty * difference_norm_square; - }; + double difference_norm_square = 0.; + for (unsigned int point = 0; point < q_points.size(); ++point) + { + const double diff = (g[point] - sol_u[point]); + difference_norm_square += diff * diff * JxW[point]; + } + copy_data.value += penalty * difference_norm_square; + }; // And finally interior face terms $\sum_{f\in \partial K}\lbrace \sigma // \left\| [u_h] \right\|_f^2 + h_f \left\| [\nu \nabla u_h \cdot // \mathbf n ] \right\|_f^2 \rbrace$: - const auto face_worker = [&](const auto &cell, - const unsigned int &f, - const unsigned int &sf, - const auto &ncell, - const unsigned int &nf, - const unsigned int &nsf, - auto &scratch_data, - auto ©_data) { - const FEInterfaceValues &fe_iv = - scratch_data.reinit(cell, f, sf, ncell, nf, nsf); + const auto face_worker = + [&](const typename DoFHandler::cell_iterator &cell, + const unsigned int &f, + const unsigned int &sf, + const typename DoFHandler::cell_iterator &ncell, + const unsigned int &nf, + const unsigned int &nsf, + ScratchData &scratch_data, + CopyData ©_data) { + const FEInterfaceValues &fe_iv = + scratch_data.reinit(cell, f, sf, ncell, nf, nsf); - copy_data.face_data.emplace_back(); - CopyDataFace ©_data_face = copy_data.face_data.back(); + copy_data.face_data.emplace_back(); + CopyDataFace ©_data_face = copy_data.face_data.back(); - copy_data_face.cell_indices[0] = cell->active_cell_index(); - copy_data_face.cell_indices[1] = ncell->active_cell_index(); + copy_data_face.cell_indices[0] = cell->active_cell_index(); + copy_data_face.cell_indices[1] = ncell->active_cell_index(); - const std::vector &JxW = fe_iv.get_JxW_values(); - const std::vector> &normals = fe_iv.get_normal_vectors(); + const std::vector &JxW = fe_iv.get_JxW_values(); + const std::vector> &normals = fe_iv.get_normal_vectors(); - const auto &q_points = fe_iv.get_quadrature_points(); - const unsigned int n_q_points = q_points.size(); + const std::vector> &q_points = fe_iv.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); - std::vector jump(n_q_points); - fe_iv.get_jump_in_function_values(solution, jump); + std::vector jump(n_q_points); + fe_iv.get_jump_in_function_values(solution, jump); - std::vector> grad_jump(n_q_points); - fe_iv.get_jump_in_function_gradients(solution, grad_jump); + std::vector> grad_jump(n_q_points); + fe_iv.get_jump_in_function_gradients(solution, grad_jump); - const double h = cell->face(f)->diameter(); + const double h = cell->face(f)->diameter(); - const double extent1 = cell->measure() / cell->face(f)->measure(); - const double extent2 = ncell->measure() / ncell->face(nf)->measure(); - const double penalty = get_penalty_factor(degree, extent1, extent2); + const double extent1 = cell->measure() / cell->face(f)->measure(); + const double extent2 = ncell->measure() / ncell->face(nf)->measure(); + const double penalty = get_penalty_factor(degree, extent1, extent2); - double flux_jump_square = 0; - double u_jump_square = 0; - for (unsigned int point = 0; point < n_q_points; ++point) - { - u_jump_square += jump[point] * jump[point] * JxW[point]; - const double flux_jump = grad_jump[point] * normals[point]; - flux_jump_square += - diffusion_coefficient * flux_jump * flux_jump * JxW[point]; - } - copy_data_face.values[0] = - 0.5 * h * (flux_jump_square + penalty * u_jump_square); - copy_data_face.values[1] = copy_data_face.values[0]; - }; + double flux_jump_square = 0; + double u_jump_square = 0; + for (unsigned int point = 0; point < n_q_points; ++point) + { + u_jump_square += jump[point] * jump[point] * JxW[point]; + const double flux_jump = grad_jump[point] * normals[point]; + flux_jump_square += + diffusion_coefficient * flux_jump * flux_jump * JxW[point]; + } + copy_data_face.values[0] = + 0.5 * h * (flux_jump_square + penalty * u_jump_square); + copy_data_face.values[1] = copy_data_face.values[0]; + }; // Having computed local contributions for each cell, we still // need a way to copy these into the global vector that will hold // the error estimators for all cells: - const auto copier = [&](const auto ©_data) { + const auto copier = [&](const CopyData ©_data) { if (copy_data.cell_index != numbers::invalid_unsigned_int) estimated_error_square_per_cell[copy_data.cell_index] += copy_data.value; - for (auto &cdf : copy_data.face_data) + for (const CopyDataFace &cdf : copy_data.face_data) for (unsigned int j = 0; j < 2; ++j) estimated_error_square_per_cell[cdf.cell_indices[j]] += cdf.values[j]; }; @@ -750,14 +760,16 @@ namespace Step74 // Assemble $\sum_{K \in \Gamma_h} \nu\|\nabla (u_h - u) \|_K^2 $. const auto cell_worker = - [&](const auto &cell, auto &scratch_data, auto ©_data) { + [&](const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch_data, + CopyData ©_data) { const FEValues &fe_v = scratch_data.reinit(cell); copy_data.cell_index = cell->active_cell_index(); - const auto &q_points = fe_v.get_quadrature_points(); - const unsigned int n_q_points = q_points.size(); - const std::vector &JxW = fe_v.get_JxW_values(); + const std::vector> &q_points = fe_v.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + const std::vector &JxW = fe_v.get_JxW_values(); std::vector> grad_u(n_q_points); fe_v.get_function_gradients(solution, grad_u); @@ -775,78 +787,80 @@ namespace Step74 }; // Assemble $\sum_{f \in F_b}\sigma \|u_h-g_D\|_f^2$. - const auto boundary_worker = [&](const auto &cell, - const unsigned int &face_no, - auto &scratch_data, - auto ©_data) { - const FEFaceValuesBase &fe_fv = scratch_data.reinit(cell, face_no); + const auto boundary_worker = + [&](const typename DoFHandler::active_cell_iterator &cell, + const unsigned int &face_no, + ScratchData &scratch_data, + CopyData ©_data) { + const FEFaceValuesBase &fe_fv = scratch_data.reinit(cell, face_no); - const auto &q_points = fe_fv.get_quadrature_points(); - const unsigned n_q_points = q_points.size(); + const std::vector> &q_points = fe_fv.get_quadrature_points(); + const unsigned n_q_points = q_points.size(); - const std::vector &JxW = fe_fv.get_JxW_values(); + const std::vector &JxW = fe_fv.get_JxW_values(); - std::vector g(n_q_points); - exact_solution->value_list(q_points, g); + std::vector g(n_q_points); + exact_solution->value_list(q_points, g); - std::vector sol_u(n_q_points); - fe_fv.get_function_values(solution, sol_u); + std::vector sol_u(n_q_points); + fe_fv.get_function_values(solution, sol_u); - const double extent1 = cell->measure() / cell->face(face_no)->measure(); - const double penalty = get_penalty_factor(degree, extent1, extent1); + const double extent1 = cell->measure() / cell->face(face_no)->measure(); + const double penalty = get_penalty_factor(degree, extent1, extent1); - double difference_norm_square = 0.; - for (unsigned int point = 0; point < q_points.size(); ++point) - { - const double diff = (g[point] - sol_u[point]); - difference_norm_square += diff * diff * JxW[point]; - } - copy_data.value += penalty * difference_norm_square; - }; + double difference_norm_square = 0.; + for (unsigned int point = 0; point < q_points.size(); ++point) + { + const double diff = (g[point] - sol_u[point]); + difference_norm_square += diff * diff * JxW[point]; + } + copy_data.value += penalty * difference_norm_square; + }; // Assemble $\sum_{f \in F_i} \sigma \| [ u_h ] \|_f^2$. - const auto face_worker = [&](const auto &cell, - const unsigned int &f, - const unsigned int &sf, - const auto &ncell, - const unsigned int &nf, - const unsigned int &nsf, - auto &scratch_data, - auto ©_data) { - const FEInterfaceValues &fe_iv = - scratch_data.reinit(cell, f, sf, ncell, nf, nsf); + const auto face_worker = + [&](const typename DoFHandler::cell_iterator &cell, + const unsigned int &f, + const unsigned int &sf, + const typename DoFHandler::cell_iterator &ncell, + const unsigned int &nf, + const unsigned int &nsf, + ScratchData &scratch_data, + CopyData ©_data) { + const FEInterfaceValues &fe_iv = + scratch_data.reinit(cell, f, sf, ncell, nf, nsf); - copy_data.face_data.emplace_back(); - CopyDataFace ©_data_face = copy_data.face_data.back(); + copy_data.face_data.emplace_back(); + CopyDataFace ©_data_face = copy_data.face_data.back(); - copy_data_face.cell_indices[0] = cell->active_cell_index(); - copy_data_face.cell_indices[1] = ncell->active_cell_index(); + copy_data_face.cell_indices[0] = cell->active_cell_index(); + copy_data_face.cell_indices[1] = ncell->active_cell_index(); - const std::vector &JxW = fe_iv.get_JxW_values(); + const std::vector &JxW = fe_iv.get_JxW_values(); - const auto &q_points = fe_iv.get_quadrature_points(); - const unsigned int n_q_points = q_points.size(); + const std::vector> &q_points = fe_iv.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); - std::vector jump(n_q_points); - fe_iv.get_jump_in_function_values(solution, jump); + std::vector jump(n_q_points); + fe_iv.get_jump_in_function_values(solution, jump); - const double extent1 = cell->measure() / cell->face(f)->measure(); - const double extent2 = ncell->measure() / ncell->face(nf)->measure(); - const double penalty = get_penalty_factor(degree, extent1, extent2); + const double extent1 = cell->measure() / cell->face(f)->measure(); + const double extent2 = ncell->measure() / ncell->face(nf)->measure(); + const double penalty = get_penalty_factor(degree, extent1, extent2); - double u_jump_square = 0; - for (unsigned int point = 0; point < n_q_points; ++point) - { - u_jump_square += jump[point] * jump[point] * JxW[point]; - } - copy_data_face.values[0] = 0.5 * penalty * u_jump_square; - copy_data_face.values[1] = copy_data_face.values[0]; - }; + double u_jump_square = 0; + for (unsigned int point = 0; point < n_q_points; ++point) + { + u_jump_square += jump[point] * jump[point] * JxW[point]; + } + copy_data_face.values[0] = 0.5 * penalty * u_jump_square; + copy_data_face.values[1] = copy_data_face.values[0]; + }; - const auto copier = [&](const auto ©_data) { + const auto copier = [&](const CopyData ©_data) { if (copy_data.cell_index != numbers::invalid_unsigned_int) energy_norm_square_per_cell[copy_data.cell_index] += copy_data.value; - for (auto &cdf : copy_data.face_data) + for (const CopyDataFace &cdf : copy_data.face_data) for (unsigned int j = 0; j < 2; ++j) energy_norm_square_per_cell[cdf.cell_indices[j]] += cdf.values[j]; }; -- 2.39.5