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
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
return partitioned;
}
+
+
void
IndexSet::subtract_set(const IndexSet &other)
{
projects / dealii.git / commitdiff