]> https://gitweb.dealii.org/ - dealii.git/commitdiff
AffineConstraints::close(): Simplify combination of constraints 16136/head
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Fri, 13 Oct 2023 18:18:56 +0000 (20:18 +0200)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Sun, 15 Oct 2023 12:38:05 +0000 (14:38 +0200)
include/deal.II/lac/affine_constraints.templates.h

index 75bf4a95064184a2595d5fdfe28634df41a1a857..3cc940ff661611cf0e911b47c944d3370e106589 100644 (file)
@@ -782,10 +782,8 @@ AffineConstraints<number>::close()
 #ifdef DEBUG
   size_type iteration = 0;
 #endif
-  bool                                chained_constraint_replaced = false;
-  std::vector<bool>                   line_finalized(lines.size(), false);
-  std::vector<const ConstraintLine *> sub_constraints;
-  std::vector<unsigned int>           entries_to_delete;
+  bool              chained_constraint_replaced = false;
+  std::vector<bool> line_finalized(lines.size(), false);
   do
     {
       chained_constraint_replaced = false;
@@ -805,128 +803,112 @@ AffineConstraints<number>::close()
             // 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.
-            sub_constraints.resize(line.entries.size());
-            bool has_sub_constraints = false;
-            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)]];
-                  has_sub_constraints = true;
-                }
-              else
-                sub_constraints[entry] = nullptr;
-
-            // If none of the entries in the current line refer to DoFs that are
-            // themselves constrained, then we can move on:
-            if (has_sub_constraints == false)
-              {
-                line_finalized[line_index] = true;
-                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();
-            entries_to_delete.clear();
+            const unsigned int n_original_entries  = line.entries.size();
+            bool               has_sub_constraints = false;
+            const size_type    lines_cache_size    = lines_cache.size();
             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
-                  [[maybe_unused]] 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 size_type dof_index =
+                  calculate_line_index(line.entries[entry].first);
 
-                  const ConstraintLine &constrained_line =
-                    *sub_constraints[entry];
-                  Assert(constrained_line.index == dof_index,
-                         ExcInternalError());
+                if (dof_index < lines_cache_size &&
+                    lines_cache[dof_index] != numbers::invalid_size_type)
+                  {
+                    // Now we found an index that is further constrained, so
+                    // we need to resolve the 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. There is no harm in
+                    // doing so: we will simply treat those in the next round
+                    // around of the outer iteration.
+                    has_sub_constraints = true;
+
+                    const number weight = line.entries[entry].second;
+                    Assert(line.entries[entry].first != line.index,
+                           ExcMessage("Cycle in constraints detected!"));
+
+                    const ConstraintLine &constrained_line =
+                      lines[lines_cache[dof_index]];
+                    Assert(constrained_line.index == line.entries[entry].first,
+                           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(line.entries[entry].first !=
+                                   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<size_type, number>(
+                          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);
 
-                  // 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,
+#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!"));
+#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. we do not want to change the loop length above we
+                      // do so by setting the 'first' entry to
+                      // invalid_size_type here to finally remove entries in a
+                      // second loop
+                      {
+                        line.entries[entry].first = numbers::invalid_size_type;
+                      }
 
-                      // replace first entry, then tack the rest to the end
-                      // of the list
-                      line.entries[entry] = std::pair<size_type, number>(
-                        constrained_line.entries[0].first,
-                        constrained_line.entries[0].second * weight);
+                    line.inhomogeneity +=
+                      constrained_line.inhomogeneity * 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);
+            // If none of the entries in the current line refer to DoFs that are
+            // themselves constrained, then we can move on:
+            if (has_sub_constraints == false)
+              line_finalized[line_index] = true;
+            else
+              {
+                chained_constraint_replaced = true;
 
-#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!"));
-#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
+                // Now delete the elements we have marked for deletion.
+                auto remaining_entries = line.entries.begin();
+                for (const auto &entry : line.entries)
+                  if (entry.first != numbers::invalid_size_type)
                     {
-                      entries_to_delete.emplace_back(entry);
+                      *remaining_entries = entry;
+                      ++remaining_entries;
                     }
-
-                  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).
-            std::sort(entries_to_delete.begin(), entries_to_delete.end());
-            for (auto it = entries_to_delete.rbegin();
-                 it != entries_to_delete.rend();
-                 ++it)
-              line.entries.erase(line.entries.begin() + *it);
+                line.entries.erase(remaining_entries, line.entries.end());
+              }
           }
 
 #ifdef DEBUG

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.