From: Wolfgang Bangerth Date: Mon, 2 Oct 2023 17:52:23 +0000 (-0600) Subject: Add a function to take a view of an index set with regard to an index set. X-Git-Tag: relicensing~435^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ed7cbfbf02a1b5b07e13aa574b75d9f5370da811;p=dealii.git Add a function to take a view of an index set with regard to an index set. --- diff --git a/include/deal.II/base/index_set.h b/include/deal.II/base/index_set.h index 5c74bb8e8f..be095d3858 100644 --- a/include/deal.II/base/index_set.h +++ b/include/deal.II/base/index_set.h @@ -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 + * within the mask. This corresponds to the notion of a view: + * The mask is a window 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 diff --git a/source/base/index_set.cc b/source/base/index_set.cc index d9174be6b9..383799cd9c 100644 --- a/source/base/index_set.cc +++ b/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 new_ranges; + + std::vector::iterator own_it = ranges.begin(); + std::vector::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::split_by_block( const std::vector &n_indices_per_block) const @@ -284,6 +422,8 @@ IndexSet::split_by_block( return partitioned; } + + void IndexSet::subtract_set(const IndexSet &other) {