Merge pull request #16079 from bangerth/view

Add a function to take a view of an index set with regard to an index set mask.

doc/news/changes/minor/20231002Bangerth

@@ -0,0 +1,4 @@
+New: The function IndexSet::get_view() now has an overload that
+computes the view with regard to an arbitrary mask.
+<br>
+(Wolfgang Bangerth, 2023/10/02)

include/deal.II/base/index_set.h

@@ -371,10 +371,61 @@ class IndexSet
    * <tt>begin</tt>. This corresponds to the notion of a <i>view</i>: The
    * interval <tt>[begin, end)</tt> is a <i>window</i> through which we see
    * the set represented by the current object.
+   *
+   * A more general function of the same name, taking a mask instead
+   * of just an interval to define the view, is below.
    */
   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.
+   *
+   * A typical case where this function is useful is as follows. Say,
+   * you have a block linear system in which you have blocks
+   * corresponding to variables $(u,p,T,c)$ (which you can think of as
+   * velocity, pressure, temperature, and chemical composition -- or
+   * whatever other kind of problem you are currently considering in
+   * your own work). We solve this in parallel, so every MPI process
+   * has its own `locally_owned_dofs` index set that describes which
+   * among all $N_\text{dofs}$ degrees of freedom this process
+   * owns. Let's assume we have developed a linear solver or
+   * preconditioner that first solves the coupled $u$-$T$ system, and
+   * once that is done, solves the $p$-$c$ system. In this case, it is
+   * often useful to set up block vectors with only two components
+   * corresponding to the $u$ and $T$ components, and later for only
+   * the $p$-$c$ components of the solution. The question is which of
+   * the components of these 2-block vectors are locally owned? The
+   * answer is that we need to get a view of the `locally_owned_dofs`
+   * index set in which we apply a mask that corresponds to the
+   * variables we're currently interested in. For the $u$-$T$ system,
+   * we need a mask (corresponding to an index set of size
+   * $N_\text{dofs}$ that contains all indices of $u$ degrees of
+   * freedom as well as all indices of $T$ degrees of freedom. The
+   * resulting view is an index set of size $N_u+N_T$ that contains
+   * the indices of the locally owned $u$ and $T$ degrees of freedom.
+   */
+  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

source/base/index_set.cc

@@ -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)
 {

tests/base/index_set_36.cc

@@ -0,0 +1,63 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2010 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// test IndexSet::get_view(IndexSet) with a mask that is identical to
+// the one used in the _18 test just that it is given as an IndexSet
+// rather than a single range.
+
+#include <deal.II/base/index_set.h>
+
+#include <stdlib.h>
+
+#include "../tests.h"
+
+
+void
+test()
+{
+  IndexSet is1(100);
+
+  // randomly add 90 elements to each
+  // set, some of which may be
+  // repetitions of previous ones
+  for (unsigned int i = 0; i < 9 * is1.size() / 10; ++i)
+    is1.add_index(Testing::rand() % is1.size());
+
+  is1.compress();
+
+  IndexSet mask(100);
+  mask.add_range(20, 50);
+
+  const IndexSet is2 = is1.get_view(mask);
+  Assert(is2.size() == mask.n_elements(), ExcInternalError());
+
+  deallog << "Original index set: " << std::endl;
+  is1.print(deallog);
+  deallog << "View of index set between 20 and 50: " << std::endl;
+  is2.print(deallog);
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  test();
+}

tests/base/index_set_36.output

@@ -0,0 +1,6 @@
+
+DEAL::Original index set: 
+DEAL::{[2,3], 5, 8, 11, [13,15], 19, [21,27], [29,30], [34,37], 40, [42,43], [45,46], [49,51], 54, [56,59], [62,64], [67,70], [72,73], [76,78], [80,84], [86,88], [90,93], [95,96], [98,99]}
+DEAL::View of index set between 20 and 50: 
+DEAL::{[1,7], [9,10], [14,17], 20, [22,23], [25,26], 29}
+DEAL::OK

tests/base/index_set_37.cc

@@ -0,0 +1,64 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2010 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// test IndexSet::get_view(IndexSet) with a mask that is not just a
+// single range.
+
+#include <deal.II/base/index_set.h>
+
+#include <stdlib.h>
+
+#include "../tests.h"
+
+
+void
+test()
+{
+  IndexSet is1(100);
+
+  // randomly add 90 elements to each
+  // set, some of which may be
+  // repetitions of previous ones
+  for (unsigned int i = 0; i < 9 * is1.size() / 10; ++i)
+    is1.add_index(Testing::rand() % is1.size());
+
+  IndexSet mask(100);
+  mask.add_range(20, 40);
+  mask.add_range(60, 80);
+  mask.compress();
+
+  const IndexSet is2 = is1.get_view(mask);
+  Assert(is2.size() == mask.n_elements(), ExcInternalError());
+
+  deallog << "Original index set: " << std::endl;
+  is1.print(deallog);
+  deallog << "Mask: " << std::endl;
+  mask.print(deallog);
+  deallog << "View of index set between 20-40 and 60-80: " << std::endl;
+  is2.print(deallog);
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  test();
+}

tests/base/index_set_37.output

@@ -0,0 +1,8 @@
+
+DEAL::Original index set: 
+DEAL::{[2,3], 5, 8, 11, [13,15], 19, [21,27], [29,30], [34,37], 40, [42,43], [45,46], [49,51], 54, [56,59], [62,64], [67,70], [72,73], [76,78], [80,84], [86,88], [90,93], [95,96], [98,99]}
+DEAL::Mask: 
+DEAL::{[20,39], [60,79]}
+DEAL::View of index set between 20-40 and 60-80: 
+DEAL::{[1,7], [9,10], [14,17], [22,24], [27,30], [32,33], [36,38]}
+DEAL::OK