]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a function to take a view of an index set with regard to an index set.
authorWolfgang Bangerth <bangerth@colostate.edu>
Mon, 2 Oct 2023 17:52:23 +0000 (11:52 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Thu, 5 Oct 2023 02:38:00 +0000 (20:38 -0600)
include/deal.II/base/index_set.h
source/base/index_set.cc

index 5c74bb8e8f57ea4360c6662e4e3cc1687e4fc1b8..be095d38580dc8f2fa6f4b3e9873515a456b114e 100644 (file)
@@ -375,6 +375,30 @@ public:
   IndexSet
   get_view(const size_type begin, const size_type end) const;
 
+  /**
+   * This command takes a "mask", i.e., a second index set of same size as the
+   * current one and returns the intersection of the current index set the mask,
+   * shifted to the index of an entry within the given mask. For example,
+   * if the current object is a an IndexSet object representing an index space
+   * `[0,100)` containing indices `[20,40)`, and if the mask represents
+   * an index space of the same size but containing all 50 *odd* indices in this
+   * range, then the result will be an index set for a space of size 50 that
+   * contains those indices that correspond to the question "the how many'th
+   * entry in the mask are the indices `[20,40)`. This will result in an index
+   * set of size 50 that contains the indices `{11,12,13,14,15,16,17,18,19,20}`
+   * (because, for example, the index 20 in the original set is not in the mask,
+   * but 21 is and corresponds to the 11th entry of the mask -- the mask
+   * contains the elements `{1,3,5,7,9,11,13,15,17,19,21,...}`).
+   *
+   * In other words, the result of this operation is the intersection of the
+   * set represented by the current object and the mask, as seen
+   * <i>within the mask</i>. This corresponds to the notion of a <i>view</i>:
+   * The mask is a <i>window</i> through which we see the set represented by the
+   * current object.
+   */
+  IndexSet
+  get_view(const IndexSet &mask) const;
+
   /**
    * Split the set indices represented by this object into blocks given by the
    * @p n_indices_per_block structure. The sum of its entries must match the
index d9174be6b9008e64fb4a1740eace35da0215917d..383799cd9c334d93eb39f246dd54ef6e76c3ce43 100644 (file)
@@ -257,6 +257,144 @@ IndexSet::get_view(const size_type begin, const size_type end) const
   return result;
 }
 
+
+
+IndexSet
+IndexSet::get_view(const IndexSet &mask) const
+{
+  Assert(size() == mask.size(),
+         ExcMessage("The mask must have the same size index space "
+                    "as the index set it is applied to."));
+
+  // If 'other' is an empty set, then the view is also empty:
+  if (mask == IndexSet())
+    return {};
+
+  // For everything, it is more efficient to work on compressed sets:
+  compress();
+  mask.compress();
+
+  // If 'other' has a single range, then we can just defer to the
+  // previous function
+  if (mask.ranges.size() == 1)
+    return get_view(mask.ranges[0].begin, mask.ranges[0].end);
+
+  // For the general case where the mask is an arbitrary set,
+  // the situation is slightly more complicated. We need to walk
+  // the ranges of the two index sets in parallel and search for
+  // overlaps, and then appropriately shift
+
+  // we save all new ranges to our IndexSet in an temporary vector and
+  // add all of them in one go at the end.
+  std::vector<Range> new_ranges;
+
+  std::vector<Range>::iterator own_it  = ranges.begin();
+  std::vector<Range>::iterator mask_it = mask.ranges.begin();
+
+  while ((own_it != ranges.end()) && (mask_it != mask.ranges.end()))
+    {
+      // If our own range lies completely ahead of the current
+      // range in the mask, move forward and start the loop body
+      // anew. If this was the last range, the 'while' loop above
+      // will terminate, so we don't have to check for end iterators
+      if (own_it->end <= mask_it->begin)
+        {
+          ++own_it;
+          continue;
+        }
+
+      // Do the same if the current mask range lies completely ahead of
+      // the current range of the this object:
+      if (mask_it->end <= own_it->begin)
+        {
+          ++mask_it;
+          continue;
+        }
+
+      // Now own_it and other_it overlap. Check that that is true by
+      // enumerating the cases that can happen. This is
+      // surprisingly tricky because the two intervals can intersect in
+      // a number of different ways, but there really are only the four
+      // following possibilities:
+
+      // Case 1: our interval overlaps the left end of the other interval
+      //
+      // So we need to add the elements from the first element of the mask's
+      // interval to the end of our own interval. But we need to shift the
+      // indices so that they correspond to the how many'th element within the
+      // mask this is; fortunately (because we compressed the mask), this
+      // is recorded in the mask's ranges.
+      if ((own_it->begin <= mask_it->begin) && (own_it->end <= mask_it->end))
+        {
+          new_ranges.emplace_back(mask_it->begin - mask_it->nth_index_in_set,
+                                  own_it->end - mask_it->nth_index_in_set);
+        }
+      else
+        // Case 2:our interval overlaps the tail end of the other interval
+        if ((mask_it->begin <= own_it->begin) && (mask_it->end <= own_it->end))
+          {
+            const size_type offset_within_mask_interval =
+              own_it->begin - mask_it->begin;
+            new_ranges.emplace_back(mask_it->nth_index_in_set +
+                                      offset_within_mask_interval,
+                                    mask_it->nth_index_in_set +
+                                      (mask_it->end - mask_it->begin));
+          }
+        else
+          // Case 3: Our own interval completely encloses the other interval
+          if ((own_it->begin <= mask_it->begin) &&
+              (own_it->end >= mask_it->end))
+            {
+              new_ranges.emplace_back(mask_it->begin -
+                                        mask_it->nth_index_in_set,
+                                      mask_it->end - mask_it->nth_index_in_set);
+            }
+          else
+            // Case 3: The other interval completely encloses our own interval
+            if ((mask_it->begin <= own_it->begin) &&
+                (mask_it->end >= own_it->end))
+              {
+                const size_type offset_within_mask_interval =
+                  own_it->begin - mask_it->begin;
+                new_ranges.emplace_back(mask_it->nth_index_in_set +
+                                          offset_within_mask_interval,
+                                        mask_it->nth_index_in_set +
+                                          offset_within_mask_interval +
+                                          (own_it->end - own_it->begin));
+              }
+            else
+              Assert(false, ExcInternalError());
+
+      // We considered the overlap of these two intervals. It may of course
+      // be that one of them overlaps with another one, but that can only
+      // be the case for the interval that extends further to the right. So
+      // we can safely move on from the interval that terminates earlier:
+      if (own_it->end < mask_it->end)
+        ++own_it;
+      else if (mask_it->end < own_it->end)
+        ++mask_it;
+      else
+        {
+          // The intervals ended at the same point. We can move on from both.
+          // (The algorithm would also work if we only moved on from one,
+          // but we can micro-optimize here without too much effort.)
+          ++own_it;
+          ++mask_it;
+        }
+    }
+
+  // Now turn the ranges of overlap we have accumulated into an IndexSet in
+  // its own right:
+  IndexSet result(mask.n_elements());
+  for (const auto &range : new_ranges)
+    result.add_range(range.begin, range.end);
+  result.compress();
+
+  return result;
+}
+
+
+
 std::vector<IndexSet>
 IndexSet::split_by_block(
   const std::vector<types::global_dof_index> &n_indices_per_block) const
@@ -284,6 +422,8 @@ IndexSet::split_by_block(
   return partitioned;
 }
 
+
+
 void
 IndexSet::subtract_set(const IndexSet &other)
 {

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.