From 7025c0c3ad2fc12c1b795157863bcc0ffbbe0b72 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 12 Oct 2023 16:37:35 -0600 Subject: [PATCH] Use AffineConstraints::add_constraint() throughout. --- examples/step-16/step-16.cc | 11 +- examples/step-37/step-37.cc | 5 +- examples/step-42/step-42.cc | 6 +- examples/step-50/step-50.cc | 15 +- examples/step-56/step-56.cc | 22 ++- examples/step-63/step-63.cc | 11 +- examples/step-66/step-66.cc | 6 +- examples/step-79/step-79.cc | 17 +- .../lac/affine_constraints.templates.h | 32 ++-- .../deal.II/multigrid/mg_constrained_dofs.h | 11 +- .../vector_tools_constraints.templates.h | 180 ++++++++++-------- source/dofs/dof_tools_constraints.cc | 170 +++++++++++------ source/grid/grid_tools.cc | 4 +- 13 files changed, 287 insertions(+), 203 deletions(-) diff --git a/examples/step-16/step-16.cc b/examples/step-16/step-16.cc index 09c4302d7d..a79e2f1fde 100644 --- a/examples/step-16/step-16.cc +++ b/examples/step-16/step-16.cc @@ -436,10 +436,13 @@ namespace Step16 boundary_constraints[level].reinit( dof_handler.locally_owned_mg_dofs(level), DoFTools::extract_locally_relevant_level_dofs(dof_handler, level)); - boundary_constraints[level].add_lines( - mg_constrained_dofs.get_refinement_edge_indices(level)); - boundary_constraints[level].add_lines( - mg_constrained_dofs.get_boundary_indices(level)); + + for (const types::global_dof_index dof_index : + mg_constrained_dofs.get_refinement_edge_indices(level)) + boundary_constraints[level].add_constraint(dof_index, {}, 0.); + for (const types::global_dof_index dof_index : + mg_constrained_dofs.get_boundary_indices(level)) + boundary_constraints[level].add_constraint(dof_index, {}, 0.); boundary_constraints[level].close(); } diff --git a/examples/step-37/step-37.cc b/examples/step-37/step-37.cc index f428beb6c6..a4e98ab431 100644 --- a/examples/step-37/step-37.cc +++ b/examples/step-37/step-37.cc @@ -852,8 +852,9 @@ namespace Step37 AffineConstraints level_constraints( dof_handler.locally_owned_mg_dofs(level), DoFTools::extract_locally_relevant_level_dofs(dof_handler, level)); - level_constraints.add_lines( - mg_constrained_dofs.get_boundary_indices(level)); + for (const types::global_dof_index dof_index : + mg_constrained_dofs.get_boundary_indices(level)) + level_constraints.add_constraint(dof_index, {}, 0.); level_constraints.close(); typename MatrixFree::AdditionalData additional_data; diff --git a/examples/step-42/step-42.cc b/examples/step-42/step-42.cc index 3fcd08064d..9d1374c245 100644 --- a/examples/step-42/step-42.cc +++ b/examples/step-42/step-42.cc @@ -1275,9 +1275,9 @@ namespace Step42 0) && !constraints_hanging_nodes.is_constrained(index_z)) { - all_constraints.add_line(index_z); - all_constraints.set_inhomogeneity(index_z, - undeformed_gap); + all_constraints.add_constraint(index_z, + {}, + undeformed_gap); distributed_solution(index_z) = undeformed_gap; active_set.add_index(index_z); diff --git a/examples/step-50/step-50.cc b/examples/step-50/step-50.cc index 6492d805cf..30d0327651 100644 --- a/examples/step-50/step-50.cc +++ b/examples/step-50/step-50.cc @@ -596,8 +596,9 @@ void LaplaceProblem::setup_multigrid() dof_handler.locally_owned_mg_dofs(level), DoFTools::extract_locally_relevant_level_dofs(dof_handler, level)); - level_constraints.add_lines( - mg_constrained_dofs.get_boundary_indices(level)); + for (const types::global_dof_index dof_index : + mg_constrained_dofs.get_boundary_indices(level)) + level_constraints.add_constraint(dof_index, {}, 0.); level_constraints.close(); typename MatrixFree::AdditionalData additional_data; @@ -826,11 +827,13 @@ void LaplaceProblem::assemble_multigrid() boundary_constraints[level].reinit( dof_handler.locally_owned_mg_dofs(level), DoFTools::extract_locally_relevant_level_dofs(dof_handler, level)); - boundary_constraints[level].add_lines( - mg_constrained_dofs.get_refinement_edge_indices(level)); - boundary_constraints[level].add_lines( - mg_constrained_dofs.get_boundary_indices(level)); + for (const types::global_dof_index dof_index : + mg_constrained_dofs.get_refinement_edge_indices(level)) + boundary_constraints[level].add_constraint(dof_index, {}, 0.); + for (const types::global_dof_index dof_index : + mg_constrained_dofs.get_boundary_indices(level)) + boundary_constraints[level].add_constraint(dof_index, {}, 0.); boundary_constraints[level].close(); } diff --git a/examples/step-56/step-56.cc b/examples/step-56/step-56.cc index 8231689a7b..1dca4bc59e 100644 --- a/examples/step-56/step-56.cc +++ b/examples/step-56/step-56.cc @@ -572,7 +572,7 @@ namespace Step56 // this here by marking the first pressure dof, which has index n_u as a // constrained dof. if (solver_type == SolverType::UMFPACK) - constraints.add_line(n_u); + constraints.add_constraint(n_u, {}, 0.); constraints.close(); } @@ -732,16 +732,22 @@ namespace Step56 triangulation.n_levels()); for (unsigned int level = 0; level < triangulation.n_levels(); ++level) { - boundary_constraints[level].add_lines( - mg_constrained_dofs.get_refinement_edge_indices(level)); - boundary_constraints[level].add_lines( - mg_constrained_dofs.get_boundary_indices(level)); + for (const types::global_dof_index dof_index : + mg_constrained_dofs.get_refinement_edge_indices(level)) + boundary_constraints[level].add_constraint(dof_index, {}, 0.); + for (const types::global_dof_index dof_index : + mg_constrained_dofs.get_boundary_indices(level)) + boundary_constraints[level].add_constraint(dof_index, {}, 0.); boundary_constraints[level].close(); - IndexSet idx = mg_constrained_dofs.get_refinement_edge_indices(level) & - mg_constrained_dofs.get_boundary_indices(level); + const IndexSet idx = + mg_constrained_dofs.get_refinement_edge_indices(level) & + mg_constrained_dofs.get_boundary_indices(level); - boundary_interface_constraints[level].add_lines(idx); + for (const types::global_dof_index dof_index : idx) + boundary_interface_constraints[level].add_constraint(dof_index, + {}, + 0.); boundary_interface_constraints[level].close(); } diff --git a/examples/step-63/step-63.cc b/examples/step-63/step-63.cc index bfead4ce5f..81c59e9f64 100644 --- a/examples/step-63/step-63.cc +++ b/examples/step-63/step-63.cc @@ -817,10 +817,13 @@ namespace Step63 boundary_constraints[level].reinit( dof_handler.locally_owned_mg_dofs(level), DoFTools::extract_locally_relevant_level_dofs(dof_handler, level)); - boundary_constraints[level].add_lines( - mg_constrained_dofs.get_refinement_edge_indices(level)); - boundary_constraints[level].add_lines( - mg_constrained_dofs.get_boundary_indices(level)); + + for (const types::global_dof_index dof_index : + mg_constrained_dofs.get_refinement_edge_indices(level)) + boundary_constraints[level].add_constraint(dof_index, {}, 0.); + for (const types::global_dof_index dof_index : + mg_constrained_dofs.get_boundary_indices(level)) + boundary_constraints[level].add_constraint(dof_index, {}, 0.); boundary_constraints[level].close(); } diff --git a/examples/step-66/step-66.cc b/examples/step-66/step-66.cc index d5c2585063..504042a676 100644 --- a/examples/step-66/step-66.cc +++ b/examples/step-66/step-66.cc @@ -593,8 +593,10 @@ namespace Step66 AffineConstraints level_constraints( dof_handler.locally_owned_mg_dofs(level), DoFTools::extract_locally_relevant_level_dofs(dof_handler, level)); - level_constraints.add_lines( - mg_constrained_dofs.get_boundary_indices(level)); + + for (const types::global_dof_index dof_index : + mg_constrained_dofs.get_boundary_indices(level)) + level_constraints.add_constraint(dof_index, {}, 0.); level_constraints.close(); typename MatrixFree::AdditionalData additional_data; diff --git a/examples/step-79/step-79.cc b/examples/step-79/step-79.cc index faf392acd5..7d4a69c333 100644 --- a/examples/step-79/step-79.cc +++ b/examples/step-79/step-79.cc @@ -536,13 +536,16 @@ namespace SAND types::global_dof_index last_density_dof = density_dofs.nth_index_in_set(density_dofs.n_elements() - 1); constraints.clear(); - constraints.add_line(last_density_dof); - for (unsigned int i = 0; i < density_dofs.n_elements() - 1; ++i) - constraints.add_entry(last_density_dof, - density_dofs.nth_index_in_set(i), - -1); - constraints.set_inhomogeneity(last_density_dof, 0); - + { + std::vector> + constraint_entries; + constraint_entries.reserve(density_dofs.n_elements() - 1); + for (const types::global_dof_index dof_index : density_dofs) + if (dof_index != last_density_dof) + constraint_entries.emplace_back(dof_index, -1.); + + constraints.add_constraint(last_density_dof, constraint_entries, 0.); + } constraints.close(); // We can now finally create the sparsity pattern for the diff --git a/include/deal.II/lac/affine_constraints.templates.h b/include/deal.II/lac/affine_constraints.templates.h index 83d9847583..a6e7abaf1d 100644 --- a/include/deal.II/lac/affine_constraints.templates.h +++ b/include/deal.II/lac/affine_constraints.templates.h @@ -568,18 +568,7 @@ AffineConstraints::make_consistent_in_parallel( // 3) refill this constraint matrix for (const auto &line : temporal_constraint_matrix) - { - // ... line - this->add_line(line.index); - - // ... inhomogeneity - if (line.inhomogeneity != number()) - this->set_inhomogeneity(line.index, line.inhomogeneity); - - // ... entries - if (!line.entries.empty()) - this->add_entries(line.index, line.entries); - } + this->add_constraint(line.index, line.entries, line.inhomogeneity); #ifdef DEBUG Assert(this->is_consistent_in_parallel( @@ -1193,8 +1182,7 @@ AffineConstraints::get_view(const IndexSet &mask) const for (const ConstraintLine &line : lines) if (mask.is_element(line.index)) { - const size_type row = mask.index_within_set(line.index); - view.add_line(row); +#ifdef DEBUG for (const std::pair &entry : line.entries) { Assert( @@ -1203,18 +1191,24 @@ AffineConstraints::get_view(const IndexSet &mask) const "In creating a view of an AffineConstraints " "object, the constraint on degree of freedom " + std::to_string(line.index) + " (which corresponds to the " + - std::to_string(row) + + std::to_string(mask.index_within_set(line.index)) + "th degree of freedom selected in the mask) " "is constrained against degree of freedom " + std::to_string(entry.first) + ", but this degree of freedom is not listed in the mask and " "consequently cannot be transcribed into the index space " "of the output object.")); - view.add_entry(row, - mask.index_within_set(entry.first), - entry.second); } - view.set_inhomogeneity(row, line.inhomogeneity); +#endif + + std::vector> translated_entries = + line.entries; + for (auto &entry : translated_entries) + entry.first = mask.index_within_set(entry.first); + + view.add_constraint(mask.index_within_set(line.index), + translated_entries, + line.inhomogeneity); } view.close(); diff --git a/include/deal.II/multigrid/mg_constrained_dofs.h b/include/deal.II/multigrid/mg_constrained_dofs.h index dd3e9a2b39..e176b208e1 100644 --- a/include/deal.II/multigrid/mg_constrained_dofs.h +++ b/include/deal.II/multigrid/mg_constrained_dofs.h @@ -356,8 +356,9 @@ MGConstrainedDoFs::initialize( !level_constraints[l].is_constrained(dofs_2[i]) && !level_constraints[l].is_constrained(dofs_1[i])) { - level_constraints[l].add_line(dofs_2[i]); - level_constraints[l].add_entry(dofs_2[i], dofs_1[i], 1.); + level_constraints[l].add_constraint(dofs_2[i], + {{dofs_1[i], 1.}}, + 0.); } } level_constraints[l].close(); @@ -624,10 +625,12 @@ MGConstrainedDoFs::merge_constraints(AffineConstraints &constraints, // merge constraints if (add_boundary_indices && this->have_boundary_indices()) - constraints.add_lines(this->get_boundary_indices(level)); + for (const auto i : this->get_boundary_indices(level)) + constraints.add_constraint(i, {}, 0.); if (add_refinement_edge_indices) - constraints.add_lines(this->get_refinement_edge_indices(level)); + for (const auto i : this->get_refinement_edge_indices(level)) + constraints.add_constraint(i, {}, 0.); if (add_level_constraints) constraints.merge(this->get_level_constraints(level), diff --git a/include/deal.II/numerics/vector_tools_constraints.templates.h b/include/deal.II/numerics/vector_tools_constraints.templates.h index 7a28a5efb8..b5a2f68186 100644 --- a/include/deal.II/numerics/vector_tools_constraints.templates.h +++ b/include/deal.II/numerics/vector_tools_constraints.templates.h @@ -170,21 +170,24 @@ namespace VectorTools if (!constraints.is_constrained(dof_indices.dof_indices[0]) && constraints.can_store_line(dof_indices.dof_indices[0])) { - constraints.add_line(dof_indices.dof_indices[0]); + const double normalized_inhomogeneity = + (std::fabs(inhomogeneity / constraining_vector[0]) > + std::numeric_limits::epsilon() ? + inhomogeneity / constraining_vector[0] : + 0); if (std::fabs(constraining_vector[1] / constraining_vector[0]) > std::numeric_limits::epsilon()) - constraints.add_entry(dof_indices.dof_indices[0], - dof_indices.dof_indices[1], - -constraining_vector[1] / - constraining_vector[0]); - - if (std::fabs(inhomogeneity / constraining_vector[0]) > - std::numeric_limits::epsilon()) - constraints.set_inhomogeneity( - dof_indices.dof_indices[0], - inhomogeneity / constraining_vector[0]); + constraints.add_constraint(dof_indices.dof_indices[0], + {{dof_indices.dof_indices[1], + -constraining_vector[1] / + constraining_vector[0]}}, + normalized_inhomogeneity); + else + constraints.add_constraint(dof_indices.dof_indices[0], + {}, + normalized_inhomogeneity); } } else @@ -192,21 +195,24 @@ namespace VectorTools if (!constraints.is_constrained(dof_indices.dof_indices[1]) && constraints.can_store_line(dof_indices.dof_indices[1])) { - constraints.add_line(dof_indices.dof_indices[1]); + const double normalized_inhomogeneity = + (std::fabs(inhomogeneity / constraining_vector[1]) > + std::numeric_limits::epsilon() ? + inhomogeneity / constraining_vector[1] : + 0); if (std::fabs(constraining_vector[0] / constraining_vector[1]) > std::numeric_limits::epsilon()) - constraints.add_entry(dof_indices.dof_indices[1], - dof_indices.dof_indices[0], - -constraining_vector[0] / - constraining_vector[1]); - - if (std::fabs(inhomogeneity / constraining_vector[1]) > - std::numeric_limits::epsilon()) - constraints.set_inhomogeneity( - dof_indices.dof_indices[1], - inhomogeneity / constraining_vector[1]); + constraints.add_constraint(dof_indices.dof_indices[1], + {{dof_indices.dof_indices[0], + -constraining_vector[0] / + constraining_vector[1]}}, + normalized_inhomogeneity); + else + constraints.add_constraint(dof_indices.dof_indices[1], + {}, + normalized_inhomogeneity); } } break; @@ -214,6 +220,16 @@ namespace VectorTools case 3: { + // Store constraints. There are at most 2 entries per + // constraint, so a boost::container::small_vector with + // static size 2 will do just fine. + // + // TODO: This could use std::inplace_vector once available (in + // C++26?) + boost::container:: + small_vector, 2> + constraint_entries; + if ((std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[1]) + 1e-10) && (std::fabs(constraining_vector[0]) >= @@ -222,29 +238,29 @@ namespace VectorTools if (!constraints.is_constrained(dof_indices.dof_indices[0]) && constraints.can_store_line(dof_indices.dof_indices[0])) { - constraints.add_line(dof_indices.dof_indices[0]); - if (std::fabs(constraining_vector[1] / constraining_vector[0]) > std::numeric_limits::epsilon()) - constraints.add_entry(dof_indices.dof_indices[0], - dof_indices.dof_indices[1], - -constraining_vector[1] / - constraining_vector[0]); + constraint_entries.emplace_back( + dof_indices.dof_indices[1], + -constraining_vector[1] / constraining_vector[0]); if (std::fabs(constraining_vector[2] / constraining_vector[0]) > std::numeric_limits::epsilon()) - constraints.add_entry(dof_indices.dof_indices[0], - dof_indices.dof_indices[2], - -constraining_vector[2] / - constraining_vector[0]); + constraint_entries.emplace_back( + dof_indices.dof_indices[2], + -constraining_vector[2] / constraining_vector[0]); - if (std::fabs(inhomogeneity / constraining_vector[0]) > - std::numeric_limits::epsilon()) - constraints.set_inhomogeneity( - dof_indices.dof_indices[0], - inhomogeneity / constraining_vector[0]); + const double normalized_inhomogeneity = + (std::fabs(inhomogeneity / constraining_vector[0]) > + std::numeric_limits::epsilon() ? + inhomogeneity / constraining_vector[0] : + 0); + + constraints.add_constraint(dof_indices.dof_indices[0], + constraint_entries, + normalized_inhomogeneity); } } else if ((std::fabs(constraining_vector[1]) + 1e-10 >= @@ -255,29 +271,29 @@ namespace VectorTools if (!constraints.is_constrained(dof_indices.dof_indices[1]) && constraints.can_store_line(dof_indices.dof_indices[1])) { - constraints.add_line(dof_indices.dof_indices[1]); - if (std::fabs(constraining_vector[0] / constraining_vector[1]) > std::numeric_limits::epsilon()) - constraints.add_entry(dof_indices.dof_indices[1], - dof_indices.dof_indices[0], - -constraining_vector[0] / - constraining_vector[1]); + constraint_entries.emplace_back( + dof_indices.dof_indices[0], + -constraining_vector[0] / constraining_vector[1]); if (std::fabs(constraining_vector[2] / constraining_vector[1]) > std::numeric_limits::epsilon()) - constraints.add_entry(dof_indices.dof_indices[1], - dof_indices.dof_indices[2], - -constraining_vector[2] / - constraining_vector[1]); + constraint_entries.emplace_back( + dof_indices.dof_indices[2], + -constraining_vector[2] / constraining_vector[1]); - if (std::fabs(inhomogeneity / constraining_vector[1]) > - std::numeric_limits::epsilon()) - constraints.set_inhomogeneity( - dof_indices.dof_indices[1], - inhomogeneity / constraining_vector[1]); + const double normalized_inhomogeneity = + (std::fabs(inhomogeneity / constraining_vector[1]) > + std::numeric_limits::epsilon() ? + inhomogeneity / constraining_vector[1] : + 0); + + constraints.add_constraint(dof_indices.dof_indices[1], + constraint_entries, + normalized_inhomogeneity); } } else @@ -285,29 +301,29 @@ namespace VectorTools if (!constraints.is_constrained(dof_indices.dof_indices[2]) && constraints.can_store_line(dof_indices.dof_indices[2])) { - constraints.add_line(dof_indices.dof_indices[2]); - if (std::fabs(constraining_vector[0] / constraining_vector[2]) > std::numeric_limits::epsilon()) - constraints.add_entry(dof_indices.dof_indices[2], - dof_indices.dof_indices[0], - -constraining_vector[0] / - constraining_vector[2]); + constraint_entries.emplace_back( + dof_indices.dof_indices[0], + -constraining_vector[0] / constraining_vector[2]); if (std::fabs(constraining_vector[1] / constraining_vector[2]) > std::numeric_limits::epsilon()) - constraints.add_entry(dof_indices.dof_indices[2], - dof_indices.dof_indices[1], - -constraining_vector[1] / - constraining_vector[2]); + constraint_entries.emplace_back( + dof_indices.dof_indices[1], + -constraining_vector[1] / constraining_vector[2]); - if (std::fabs(inhomogeneity / constraining_vector[2]) > - std::numeric_limits::epsilon()) - constraints.set_inhomogeneity( - dof_indices.dof_indices[2], - inhomogeneity / constraining_vector[2]); + const double normalized_inhomogeneity = + (std::fabs(inhomogeneity / constraining_vector[2]) > + std::numeric_limits::epsilon() ? + inhomogeneity / constraining_vector[2] : + 0); + + constraints.add_constraint(dof_indices.dof_indices[2], + constraint_entries, + normalized_inhomogeneity); } } @@ -363,25 +379,29 @@ namespace VectorTools if (!constraints.is_constrained(dof_indices.dof_indices[d]) && constraints.can_store_line(dof_indices.dof_indices[d])) { - constraints.add_line(dof_indices.dof_indices[d]); - - if (std::fabs(tangent_vector[d] / - tangent_vector[largest_component]) > - std::numeric_limits::epsilon()) - constraints.add_entry( - dof_indices.dof_indices[d], - dof_indices.dof_indices[largest_component], - tangent_vector[d] / tangent_vector[largest_component]); - const double inhomogeneity = (b_values(d) * tangent_vector[largest_component] - b_values(largest_component) * tangent_vector[d]) / tangent_vector[largest_component]; - if (std::fabs(inhomogeneity) > + const double normalized_inhomogeneity = + (std::fabs(inhomogeneity) > + std::numeric_limits::epsilon() ? + inhomogeneity : + 0); + + if (std::fabs(tangent_vector[d] / + tangent_vector[largest_component]) > std::numeric_limits::epsilon()) - constraints.set_inhomogeneity(dof_indices.dof_indices[d], - inhomogeneity); + constraints.add_constraint( + dof_indices.dof_indices[d], + {{dof_indices.dof_indices[largest_component], + tangent_vector[d] / tangent_vector[largest_component]}}, + normalized_inhomogeneity); + else + constraints.add_constraint(dof_indices.dof_indices[d], + {}, + normalized_inhomogeneity); } } diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index 716de8eca8..6ec936f7b2 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -509,8 +509,14 @@ namespace DoFTools Assert(primary_dofs[col] != numbers::invalid_dof_index, ExcInternalError()); - std::vector< - std::pair::size_type, number2>> + // Build constraints in a vector of pairs that can be + // arbitrarily large, but that holds up to 25 elements without + // external memory allocation. This is good enough for hanging + // node constraints of Q4 elements in 3d, so covers most + // common cases. + boost::container::small_vector< + std::pair::size_type, number2>, + 25> entries; entries.reserve(n_primary_dofs); for (unsigned int row = 0; row != n_dependent_dofs; ++row) @@ -650,6 +656,16 @@ namespace DoFTools std::vector dofs_on_mother; std::vector dofs_on_children; + // Build constraints in a vector of pairs that can be + // arbitrarily large, but that holds up to 25 elements without + // external memory allocation. This is good enough for hanging + // node constraints of Q4 elements in 3d, so covers most + // common cases. + boost::container::small_vector< + std::pair::size_type, number>, + 25> + constraint_entries; + // loop over all lines; only on lines there can be constraints. We do so // by looping over all active cells and checking whether any of the faces // are refined which can only be from the neighboring cell because this @@ -658,10 +674,7 @@ namespace DoFTools // note that even though we may visit a face twice if the neighboring // cells are equally refined, we can only visit each face with hanging // nodes once - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) + for (const auto &cell : dof_handler.active_cell_iterators()) { // artificial cells can at best neighbor ghost cells, but we're not // interested in these interfaces @@ -760,13 +773,15 @@ namespace DoFTools for (unsigned int row = 0; row != dofs_on_children.size(); ++row) { - constraints.add_line(dofs_on_children[row]); + constraint_entries.clear(); + constraint_entries.reserve(dofs_on_mother.size()); for (unsigned int i = 0; i != dofs_on_mother.size(); ++i) - constraints.add_entry(dofs_on_children[row], - dofs_on_mother[i], - fe.constraints()(row, i)); + constraint_entries.emplace_back(dofs_on_mother[i], + fe.constraints()(row, i)); - constraints.set_inhomogeneity(dofs_on_children[row], 0.); + constraints.add_constraint(dofs_on_children[row], + constraint_entries, + 0.); } } else @@ -802,6 +817,16 @@ namespace DoFTools std::vector dofs_on_mother; std::vector dofs_on_children; + // Build constraints in a vector of pairs that can be + // arbitrarily large, but that holds up to 25 elements without + // external memory allocation. This is good enough for hanging + // node constraints of Q4 elements in 3d, so covers most + // common cases. + boost::container::small_vector< + std::pair::size_type, number>, + 25> + constraint_entries; + // loop over all quads; only on quads there can be constraints. We do so // by looping over all active cells and checking whether any of the faces // are refined which can only be from the neighboring cell because this @@ -810,10 +835,7 @@ namespace DoFTools // note that even though we may visit a face twice if the neighboring // cells are equally refined, we can only visit each face with hanging // nodes once - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) + for (const auto &cell : dof_handler.active_cell_iterators()) { // artificial cells can at best neighbor ghost cells, but we're not // interested in these interfaces @@ -1004,13 +1026,15 @@ namespace DoFTools for (unsigned int row = 0; row != dofs_on_children.size(); ++row) { - constraints.add_line(dofs_on_children[row]); + constraint_entries.clear(); + constraint_entries.reserve(dofs_on_mother.size()); for (unsigned int i = 0; i != dofs_on_mother.size(); ++i) - constraints.add_entry(dofs_on_children[row], - dofs_on_mother[i], - fe.constraints()(row, i)); + constraint_entries.emplace_back(dofs_on_mother[i], + fe.constraints()(row, i)); - constraints.set_inhomogeneity(dofs_on_children[row], 0.); + constraints.add_constraint(dofs_on_children[row], + constraint_entries, + 0.); } } else @@ -1993,11 +2017,20 @@ namespace DoFTools cell_to_rotated_face_index[cell_index] = i; } + // Build constraints in a vector of pairs that can be + // arbitrarily large, but that holds up to 25 elements without + // external memory allocation. This is good enough for hanging + // node constraints of Q4 elements in 3d, so covers most + // common cases. + boost::container::small_vector< + std::pair::size_type, number>, + 25> + constraint_entries; + // // Loop over all dofs on face 2 and constrain them against all // matching dofs on face 1: // - for (unsigned int i = 0; i < dofs_per_face; ++i) { // Obey the component mask @@ -2055,8 +2088,8 @@ namespace DoFTools if (affine_constraints.is_constrained(dofs_2[i])) continue; - // Enter the constraint piece by piece: - affine_constraints.add_line(dofs_2[i]); + constraint_entries.clear(); + constraint_entries.reserve(dofs_per_face); for (unsigned int jj = 0; jj < dofs_per_face; ++jj) { @@ -2067,11 +2100,16 @@ namespace DoFTools jj, 0, face_orientation, face_flip, face_rotation)]; if (std::abs(transformation(i, jj)) > eps) - affine_constraints.add_entry(dofs_2[i], - dofs_1[j], - transformation(i, jj)); + constraint_entries.emplace_back(dofs_1[j], + transformation(i, jj)); } + // Enter the constraint:: + affine_constraints.add_constraint(dofs_2[i], + constraint_entries, + 0.); + + // Continue with next dof. continue; } @@ -2149,14 +2187,12 @@ namespace DoFTools if (constraints_are_cyclic) { if (std::abs(cycle_constraint_factor - number(1.)) > eps) - affine_constraints.add_line(dof_left); + affine_constraints.add_constraint(dof_left, {}, 0.); } else { - affine_constraints.add_line(dof_left); - affine_constraints.add_entry(dof_left, - dof_right, - constraint_factor); + affine_constraints.add_constraint( + dof_left, {{dof_right, constraint_factor}}, 0.); // The number 1e10 in the assert below is arbitrary. If the // absolute value of constraint_factor is too large, then probably // the absolute value of periodicity_factor is too large or too @@ -3339,8 +3375,6 @@ namespace DoFTools // otherwise enter all constraints - constraints.add_line(global_dof); - constraint_line.clear(); for (types::global_dof_index row = first_used_row; row < n_coarse_dofs; @@ -3352,7 +3386,7 @@ namespace DoFTools constraint_line.emplace_back(representants[row], j->second); } - constraints.add_entries(global_dof, constraint_line); + constraints.add_constraint(global_dof, constraint_line, 0.); } } @@ -3483,10 +3517,14 @@ namespace DoFTools std::vector cell_dofs; cell_dofs.reserve(dof.get_fe_collection().max_dofs_per_cell()); - typename DoFHandler::active_cell_iterator - cell = dof.begin_active(), - endc = dof.end(); - for (; cell != endc; ++cell) + // In looping over faces, we will encounter some DoFs multiple + // times (namely, the ones on vertices and (in 3d) edges shared + // between multiple boundary faces. Keep track of which DoFs we + // have already encountered, so that we do not have to consider + // them a second time. + std::set dofs_already_treated; + + for (const auto &cell : dof.active_cell_iterators()) if (!cell->is_artificial() && cell->at_boundary()) { const FiniteElement &fe = cell->get_fe(); @@ -3513,30 +3551,38 @@ namespace DoFTools // enter those dofs into the list that match the component // signature. for (const types::global_dof_index face_dof : face_dofs) - { - // Find out if a dof has a contribution in this component, - // and if so, add it to the list - const std::vector::iterator - it_index_on_cell = std::find(cell_dofs.begin(), - cell_dofs.end(), - face_dof); - Assert(it_index_on_cell != cell_dofs.end(), - ExcInvalidIterator()); - const unsigned int index_on_cell = - std::distance(cell_dofs.begin(), it_index_on_cell); - const ComponentMask &nonzero_component_array = - cell->get_fe().get_nonzero_components(index_on_cell); - bool nonzero = false; - for (unsigned int c = 0; c < n_components; ++c) - if (nonzero_component_array[c] && component_mask[c]) - { - nonzero = true; - break; - } - - if (nonzero) - zero_boundary_constraints.add_line(face_dof); - } + if (dofs_already_treated.find(face_dof) == + dofs_already_treated.end()) + { + // Find out if a dof has a contribution in this + // component, and if so, add it to the list + const std::vector::iterator + it_index_on_cell = std::find(cell_dofs.begin(), + cell_dofs.end(), + face_dof); + Assert(it_index_on_cell != cell_dofs.end(), + ExcInvalidIterator()); + const unsigned int index_on_cell = + std::distance(cell_dofs.begin(), it_index_on_cell); + const ComponentMask &nonzero_component_array = + cell->get_fe().get_nonzero_components(index_on_cell); + + bool nonzero = false; + for (unsigned int c = 0; c < n_components; ++c) + if (nonzero_component_array[c] && component_mask[c]) + { + nonzero = true; + break; + } + + if (nonzero) + zero_boundary_constraints.add_constraint(face_dof, + {}, + 0.); + // We already dealt with this DoF. Make sure we + // don't touch it again. + dofs_already_treated.insert(face_dof); + } } } } diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 6abda211d9..b62a661b90 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2297,9 +2297,9 @@ namespace GridTools if (map_iter != map_end) for (unsigned int i = 0; i < dim; ++i) { - constraints[i].add_line(cell->vertex_dof_index(vertex_no, 0)); - constraints[i].set_inhomogeneity( + constraints[i].add_constraint( cell->vertex_dof_index(vertex_no, 0), + {}, (solve_for_absolute_positions ? map_iter->second(i) : map_iter->second(i) - vertex_point[i])); -- 2.39.5