From bfce59424b679187a1bb5e2161df1d8ad33ac0a3 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 10 Jul 2023 23:08:16 -0600 Subject: [PATCH] Keep track of constraints we have already determined are finalized. --- .../lac/affine_constraints.templates.h | 240 +++++++++--------- 1 file changed, 124 insertions(+), 116 deletions(-) diff --git a/include/deal.II/lac/affine_constraints.templates.h b/include/deal.II/lac/affine_constraints.templates.h index 834384fa66..f9a6e01eb6 100644 --- a/include/deal.II/lac/affine_constraints.templates.h +++ b/include/deal.II/lac/affine_constraints.templates.h @@ -712,139 +712,147 @@ AffineConstraints::close() #ifdef DEBUG size_type iteration = 0; #endif - bool chained_constraint_replaced = false; + bool chained_constraint_replaced = false; + std::vector line_finalized(lines.size(), false); do { chained_constraint_replaced = false; - for (ConstraintLine &line : lines) - { + for (unsigned int line_index = 0; line_index < lines.size(); ++line_index) + if (line_finalized[line_index] == false) + { + ConstraintLine &line = lines[line_index]; + #ifdef DEBUG - // we need to keep track of how many replacements we do in this line, - // because we can end up in a cycle A->B->C->A without the number of - // entries growing. - size_type n_replacements = 0; + // we need to keep track of how many replacements we do in this + // line, because we can end up in a cycle A->B->C->A without the + // number of entries growing. + size_type n_replacements = 0; #endif - // loop over all entries of this line (including ones that we - // have appended in this go around) and see whether they are - // further constrained. ignore elements that we don't store on - // the current processor. - // - // If we find that one of the entries in the current line is indeed - // constrained, then store the address of the line it corresponds to - // for further processing; this avoids having to look it up again, - // which requires an IndexSet lookup. If not, store a nullptr. - std::vector sub_constraints( - line.entries.size(), nullptr); - for (unsigned int entry = 0; entry < line.entries.size(); ++entry) - if (((local_lines.size() == 0) || - (local_lines.is_element(line.entries[entry].first))) && - is_constrained(line.entries[entry].first)) + // loop over all entries of this line (including ones that we + // have appended in this go around) and see whether they are + // further constrained. ignore elements that we don't store on + // the current processor. + // + // If we find that one of the entries in the current line is indeed + // constrained, then store the address of the line it corresponds to + // for further processing; this avoids having to look it up again, + // which requires an IndexSet lookup. If not, store a nullptr. + std::vector sub_constraints( + line.entries.size(), nullptr); + for (unsigned int entry = 0; entry < line.entries.size(); ++entry) + if (((local_lines.size() == 0) || + (local_lines.is_element(line.entries[entry].first))) && + is_constrained(line.entries[entry].first)) + { + const size_type dof_index = line.entries[entry].first; + sub_constraints[entry] = + &lines[lines_cache[calculate_line_index(dof_index)]]; + } + + // If none of the entries in the current line refer to DoFs that are + // themselves constrained, then we can move on: + if (std::none_of(sub_constraints.begin(), + sub_constraints.end(), + [](const auto p) { return (p != nullptr); })) { - const size_type dof_index = line.entries[entry].first; - sub_constraints[entry] = - &lines[lines_cache[calculate_line_index(dof_index)]]; + line_finalized[line_index] = true; + continue; } - // If none of the entries in the current line refer to DoFs that are - // themselves constrained, then we can move on: - if (std::none_of(sub_constraints.begin(), - sub_constraints.end(), - [](const auto p) { return (p != nullptr); })) - continue; + // Now walk through the original entries. We may replace some of + // them, and we may add some, but we only walk through the original + // entries (and also don't touch the replacements) even though in + // principle we could also do the ones we have added or replaced. We + // don't because we don't have information for those in the + // sub_constraints array above. But there is no harm in doing so: + // we will simply treat those in the next round around of the outer + // iteration. + // + // Since we want to delete some entries (see below) but don't want + // to disturb the order of the correspondence between lines.entries + // and sub_constraints, we store up which entries we will later + // have to delete. + const unsigned int n_original_entries = line.entries.size(); + std::set entries_to_delete; + for (unsigned int entry = 0; entry < n_original_entries; ++entry) + if (sub_constraints[entry] != nullptr) + { + // ok, this entry is further constrained: + chained_constraint_replaced = true; - // Now walk through the original entries. We may replace some of them, - // and we may add some, but we only walk through the original entries - // (and also don't touch the replacements) even though in principle - // we could also do the ones we have added or replaced. - // We don't because we don't have information for those in the - // sub_constraints array above. But there is no harm in doing so: - // we will simply treat those in the next round around of the outer - // iteration. - // - // Since we want to delete some entries (see below) but don't want - // to disturb the order of the correspondence between lines.entries - // and sub_constraints, we store up which entries we will later - // have to delete. - const unsigned int n_original_entries = line.entries.size(); - std::set entries_to_delete; - for (unsigned int entry = 0; entry < n_original_entries; ++entry) - if (sub_constraints[entry] != nullptr) - { - // ok, this entry is further constrained: - chained_constraint_replaced = true; - - // look up the chain of constraints for this entry - const size_type dof_index = line.entries[entry].first; - const number weight = line.entries[entry].second; - - Assert(dof_index != line.index, - ExcMessage("Cycle in constraints detected!")); - - const ConstraintLine &constrained_line = - *sub_constraints[entry]; - Assert(constrained_line.index == dof_index, ExcInternalError()); - - // now we have to replace an entry by its expansion. we do - // that by overwriting the entry by the first entry of the - // expansion and adding the remaining ones to the end, - // where we will later process them once more - // - // we can of course only do that if the DoF that we are - // currently handling is constrained by a linear combination - // of other dofs: - if (constrained_line.entries.size() > 0) - { - for (size_type i = 0; i < constrained_line.entries.size(); - ++i) - Assert(dof_index != constrained_line.entries[i].first, - ExcMessage("Cycle in constraints detected!")); + // look up the chain of constraints for this entry + const size_type dof_index = line.entries[entry].first; + const number weight = line.entries[entry].second; - // replace first entry, then tack the rest to the end - // of the list - line.entries[entry] = std::pair( - constrained_line.entries[0].first, - constrained_line.entries[0].second * weight); + Assert(dof_index != line.index, + ExcMessage("Cycle in constraints detected!")); - for (size_type i = 1; i < constrained_line.entries.size(); - ++i) - line.entries.emplace_back( - constrained_line.entries[i].first, - constrained_line.entries[i].second * weight); + const ConstraintLine &constrained_line = + *sub_constraints[entry]; + Assert(constrained_line.index == dof_index, + ExcInternalError()); + + // now we have to replace an entry by its expansion. we do + // that by overwriting the entry by the first entry of the + // expansion and adding the remaining ones to the end, + // where we will later process them once more + // + // we can of course only do that if the DoF that we are + // currently handling is constrained by a linear combination + // of other dofs: + if (constrained_line.entries.size() > 0) + { + for (size_type i = 0; i < constrained_line.entries.size(); + ++i) + Assert(dof_index != constrained_line.entries[i].first, + ExcMessage("Cycle in constraints detected!")); + + // replace first entry, then tack the rest to the end + // of the list + line.entries[entry] = std::pair( + constrained_line.entries[0].first, + constrained_line.entries[0].second * weight); + + for (size_type i = 1; i < constrained_line.entries.size(); + ++i) + line.entries.emplace_back( + constrained_line.entries[i].first, + constrained_line.entries[i].second * weight); #ifdef DEBUG - // keep track of how many entries we replace in this - // line. If we do more than there are constraints or - // dofs in our system, we must have a cycle. - ++n_replacements; - Assert(n_replacements / 2 < largest_idx, - ExcMessage("Cycle in constraints detected!")); + // keep track of how many entries we replace in this + // line. If we do more than there are constraints or + // dofs in our system, we must have a cycle. + ++n_replacements; + Assert(n_replacements / 2 < largest_idx, + ExcMessage("Cycle in constraints detected!")); #endif - } - else - // the DoF that we encountered is not constrained by a - // linear combination of other dofs but is equal to just - // the inhomogeneity (i.e. its chain of entries is - // empty). in that case, we can't just overwrite the - // current entry, but we have to actually eliminate it - { - entries_to_delete.insert(entry); - } + } + else + // the DoF that we encountered is not constrained by a + // linear combination of other dofs but is equal to just + // the inhomogeneity (i.e. its chain of entries is + // empty). in that case, we can't just overwrite the + // current entry, but we have to actually eliminate it + { + entries_to_delete.insert(entry); + } - line.inhomogeneity += constrained_line.inhomogeneity * weight; - } + line.inhomogeneity += constrained_line.inhomogeneity * weight; + } - // Now delete the elements we have marked for deletion. Do - // so in reverse order so that we compress the array walking - // backward from the end without having to keep track that - // we have already erased earlier elements (as we would have to - // do if we walked the list of entries to delete in forward - // order). - for (auto it = entries_to_delete.rbegin(); - it != entries_to_delete.rend(); - ++it) - line.entries.erase(line.entries.begin() + *it); - } + // Now delete the elements we have marked for deletion. Do + // so in reverse order so that we compress the array walking + // backward from the end without having to keep track that + // we have already erased earlier elements (as we would have to + // do if we walked the list of entries to delete in forward + // order). + for (auto it = entries_to_delete.rbegin(); + it != entries_to_delete.rend(); + ++it) + line.entries.erase(line.entries.begin() + *it); + } #ifdef DEBUG // increase iteration count. note that we should not iterate more -- 2.39.5