1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15
16 #ifndef dealii_index_set_h
17 #define dealii_index_set_h
18
19 #include <deal.II/base/config.h>
20
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/thread_management.h>
23 #include <deal.II/base/utilities.h>
24
25 #include <boost/serialization/vector.hpp>
26
27 #include <algorithm>
28 #include <iterator>
29 #include <vector>
30
31
32 #ifdef DEAL_II_WITH_TRILINOS
33 # include <Epetra_Map.h>
34 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
35 # include <Tpetra_Map.hpp>
36 # endif
37 #endif
38
39 #if defined(DEAL_II_WITH_MPI) || defined(DEAL_II_WITH_PETSC)
40 # include <mpi.h>
41 #else
42 using MPI_Comm = int;
43 # ifndef MPI_COMM_WORLD
44 # define MPI_COMM_WORLD 0
45 # endif
46 #endif
47
48 DEAL_II_NAMESPACE_OPEN
49
50 /**
51 * A class that represents a subset of indices among a larger set. For
52 * example, it can be used to denote the set of degrees of freedom within the
53 * range $[0,\mathrm{dof\_handler.n\_dofs()})$ that belongs to a particular
54 * subdomain, or those among all degrees of freedom that are stored on a
55 * particular processor in a distributed parallel computation.
56 *
57 * This class can represent a collection of half-open ranges of indices as
58 * well as individual elements. For practical purposes it also stores the
59 * overall range these indices can assume. In other words, you need to specify
60 * the size of the index space $[0,\text{size})$ of which objects of this
61 * class are a subset.
62 *
63 * There are two ways to iterate over the IndexSets: First, begin() and end()
64 * allow iteration over individual indices in the set. Second,
65 * begin_interval() and end_interval() allow iteration over the half-open
66 * ranges as described above.
67 *
68 * The data structures used in this class along with a rationale can be found
69 * in the
70 * @ref distributed_paper "Distributed Computing paper".
71 */
72 class IndexSet
73 {
74 public:
75 // forward declarations:
76 class ElementIterator;
77 class IntervalIterator;
78
79 /**
80 * @p size_type is the type used for storing the size and the individual
81 * entries in the IndexSet.
82 */
83 using size_type = types::global_dof_index;
84
85 /**
86 * One can see an IndexSet as a container of size size(), where the elements
87 * of the containers are bool values that are either false or true,
88 * depending on whether a particular index is an element of the IndexSet or
89 * not. In other words, an IndexSet is a bit like a vector in which the
90 * elements we store are booleans. In this view, the correct local alias
91 * indicating the type of the elements of the vector would then be @p bool.
92 *
93 * On the other hand, @p bool has the disadvantage that it is not a
94 * numerical type that, for example, allows multiplication with a @p double.
95 * In other words, one can not easily use a vector of booleans in a place
96 * where other vectors are allowed. Consequently, we declare the type of the
97 * elements of such a vector as a signed integer. This uses the fact that in
98 * the C++ language, booleans are implicitly convertible to integers. In
99 * other words, declaring the type of the elements of the vector as a signed
100 * integer is only a small lie, but it is a useful one.
101 */
102 using value_type = signed int;
103
104
105 /**
106 * Default constructor.
107 */
108 IndexSet();
109
110 /**
111 * Constructor that also sets the overall size of the index range.
112 */
113 explicit IndexSet(const size_type size);
114
115 /**
116 * Copy constructor.
117 */
118 IndexSet(const IndexSet &) = default;
119
120 /**
121 * Copy assignment operator.
122 */
123 IndexSet &
124 operator=(const IndexSet &) = default;
125
126 /**
127 * Move constructor. Create a new IndexSet by transferring the internal data
128 * of the input set.
129 */
130 IndexSet(IndexSet &&is) noexcept;
131
132 /**
133 * Move assignment operator. Transfer the internal data of the input set into
134 * the current one.
135 */
136 IndexSet &
137 operator=(IndexSet &&is) noexcept;
138
139 #ifdef DEAL_II_WITH_TRILINOS
140 /**
141 * Constructor from a Trilinos Epetra_BlockMap.
142 */
143 explicit IndexSet(const Epetra_BlockMap &map);
144 #endif
145
146 /**
147 * Remove all indices from this index set. The index set retains its size,
148 * however.
149 */
150 void
151 clear();
152
153 /**
154 * Set the maximal size of the indices upon which this object operates.
155 *
156 * This function can only be called if the index set does not yet contain
157 * any elements. This can be achieved by calling clear(), for example.
158 */
159 void
160 set_size(const size_type size);
161
162 /**
163 * Return the size of the index space of which this index set is a subset
164 * of.
165 *
166 * Note that the result is not equal to the number of indices within this
167 * set. The latter information is returned by n_elements().
168 */
169 size_type
170 size() const;
171
172 /**
173 * Add the half-open range $[\text{begin},\text{end})$ to the set of indices
174 * represented by this class.
175 * @param[in] begin The first element of the range to be added.
176 * @param[in] end The past-the-end element of the range to be added.
177 */
178 void
179 add_range(const size_type begin, const size_type end);
180
181 /**
182 * Add an individual index to the set of indices.
183 */
184 void
185 add_index(const size_type index);
186
187 /**
188 * Add a whole set of indices described by dereferencing every element of
189 * the iterator range <code>[begin,end)</code>.
190 *
191 * @param[in] begin Iterator to the first element of range of indices to be
192 * added
193 * @param[in] end The past-the-end iterator for the range of elements to be
194 * added. @pre The condition <code>begin@<=end</code> needs to be satisfied.
195 */
196 template <typename ForwardIterator>
197 void
198 add_indices(const ForwardIterator &begin, const ForwardIterator &end);
199
200 /**
201 * Add the given IndexSet @p other to the current one, constructing the
202 * union of *this and @p other.
203 *
204 * If the @p offset argument is nonzero, then every index in @p other is
205 * shifted by @p offset before being added to the current index set. This
206 * allows to construct, for example, one index set from several others that
207 * are supposed to represent index sets corresponding to different ranges
208 * (e.g., when constructing the set of nonzero entries of a block vector
209 * from the sets of nonzero elements of the individual blocks of a vector).
210 *
211 * This function will generate an exception if any of the (possibly shifted)
212 * indices of the @p other index set lie outside the range
213 * <code>[0,size())</code> represented by the current object.
214 */
215 void
216 add_indices(const IndexSet &other, const size_type offset = 0);
217
218 /**
219 * Return whether the specified index is an element of the index set.
220 */
221 bool
222 is_element(const size_type index) const;
223
224 /**
225 * Return whether the index set stored by this object defines a contiguous
226 * range. This is true also if no indices are stored at all.
227 */
228 bool
229 is_contiguous() const;
230
231 /**
232 * Return whether the index set stored by this object contains no elements.
233 * This is similar, but faster than checking <code>n_elements() == 0</code>.
234 */
235 bool
236 is_empty() const;
237
238 /**
239 * Return whether the IndexSets are ascending with respect to MPI process
240 * number and 1:1, i.e., each index is contained in exactly one IndexSet
241 * (among those stored on the different processes), each process stores
242 * contiguous subset of indices, and the index set on process $p+1$ starts
243 * at the index one larger than the last one stored on process $p$.
244 * In case there is only one MPI process, this just means that the IndexSet
245 * is complete.
246 */
247 bool
248 is_ascending_and_one_to_one(const MPI_Comm &communicator) const;
249
250 /**
251 * Return the number of elements stored in this index set.
252 */
253 size_type
254 n_elements() const;
255
256 /**
257 * Return the global index of the local index with number @p local_index
258 * stored in this index set. @p local_index obviously needs to be less than
259 * n_elements().
260 */
261 size_type
262 nth_index_in_set(const size_type local_index) const;
263
264 /**
265 * Return the how-manyth element of this set (counted in ascending order) @p
266 * global_index is. @p global_index needs to be less than the size(). This
267 * function returns numbers::invalid_dof_index if the index @p global_index is not actually
268 * a member of this index set, i.e. if is_element(global_index) is false.
269 */
270 size_type
271 index_within_set(const size_type global_index) const;
272
273 /**
274 * Each index set can be represented as the union of a number of contiguous
275 * intervals of indices, where if necessary intervals may only consist of
276 * individual elements to represent isolated members of the index set.
277 *
278 * This function returns the minimal number of such intervals that are
279 * needed to represent the index set under consideration.
280 */
281 unsigned int
282 n_intervals() const;
283
284 /**
285 * This function returns the local index of the beginning of the largest
286 * range.
287 *
288 * In other words, the return value is nth_index_in_set(x), where x is the
289 * first index of the largest contiguous range of indices in the
290 * IndexSet. The return value is therefore equal to the number of elements
291 * in the set that come before the largest range.
292 *
293 * This call assumes that the IndexSet is nonempty.
294 */
295 size_type
296 largest_range_starting_index() const;
297
298 /**
299 * Compress the internal representation by merging individual elements with
300 * contiguous ranges, etc. This function does not have any external effect.
301 */
302 void
303 compress() const;
304
305 /**
306 * Comparison for equality of index sets. This operation is only allowed if
307 * the size of the two sets is the same (though of course they do not have
308 * to have the same number of indices).
309 */
310 bool
311 operator==(const IndexSet &is) const;
312
313 /**
314 * Comparison for inequality of index sets. This operation is only allowed
315 * if the size of the two sets is the same (though of course they do not
316 * have to have the same number of indices).
317 */
318 bool
319 operator!=(const IndexSet &is) const;
320
321 /**
322 * Return the intersection of the current index set and the argument given,
323 * i.e. a set of indices that are elements of both index sets. The two index
324 * sets must have the same size (though of course they do not have to have
325 * the same number of indices).
326 */
327 IndexSet operator&(const IndexSet &is) const;
328
329 /**
330 * This command takes an interval <tt>[begin, end)</tt> and returns the
331 * intersection of the current index set with the interval, shifted to the
332 * range <tt>[0, end-begin)</tt>.
333 *
334 * In other words, the result of this operation is the intersection of the
335 * set represented by the current object and the interval <tt>[begin,
336 * end)</tt>, as seen <i>within the interval <tt>[begin, end)</tt></i> by
337 * shifting the result of the intersection operation to the left by
338 * <tt>begin</tt>. This corresponds to the notion of a <i>view</i>: The
339 * interval <tt>[begin, end)</tt> is a <i>window</i> through which we see
340 * the set represented by the current object.
341 */
342 IndexSet
343 get_view(const size_type begin, const size_type end) const;
344
345 /**
346 * Split the set indices represented by this object into blocks given by the
347 * @p n_indices_per_block structure. The sum of its entries must match the
348 * global size of the current object.
349 */
350 std::vector<IndexSet>
351 split_by_block(
352 const std::vector<types::global_dof_index> &n_indices_per_block) const;
353
354 /**
355 * Remove all elements contained in @p other from this set. In other words,
356 * if $x$ is the current object and $o$ the argument, then we compute $x
357 * \leftarrow x \backslash o$.
358 */
359 void
360 subtract_set(const IndexSet &other);
361
362 /**
363 * Return a new IndexSet, with global size equal to
364 * `this->size()*other.size()`, containing for every element `n` of this
365 * IndexSet, the entries in the half open range `[n*other.size(),
366 * (n+1)*other.size())` of the @p other IndexSet.
367 *
368 * The name results from the perspective that one starts with an IndexSet and
369 * takes the tensor product with another IndexSet with `other.size()`
370 * elements; this results in a matrix of size `this->size()` times
371 * `other.size()` that has ones in exactly the rows for which this IndexSet
372 * contained an index and in the columns for which the @p other IndexSet
373 * contained an index. This matrix is then "unrolled" again by going through
374 * each row one by one and reindexing the entries of the matrix in consecutive
375 * order. A one in the matrix then corresponds to an entry in the reindexed
376 * IndexSet that is returned by this function.
377 */
378 IndexSet
379 tensor_product(const IndexSet &other) const;
380
381 /**
382 * Remove and return the last element of the last range.
383 * This function throws an exception if the IndexSet is empty.
384 */
385 size_type
386 pop_back();
387
388 /**
389 * Remove and return the first element of the first range.
390 * This function throws an exception if the IndexSet is empty.
391 */
392 size_type
393 pop_front();
394
395 /**
396 * Fill the given vector with all indices contained in this IndexSet.
397 */
398 void
399 fill_index_vector(std::vector<size_type> &indices) const;
400
401 /**
402 * Fill the given vector with either zero or one elements, providing a
403 * binary representation of this index set. The given vector is assumed to
404 * already have the correct size.
405 *
406 * The given argument is filled with integer values zero and one, using
407 * <code>vector.operator[]</code>. Thus, any object that has such an
408 * operator can be used as long as it allows conversion of integers zero and
409 * one to elements of the vector. Specifically, this is the case for classes
410 * Vector, BlockVector, but also std::vector@<bool@>, std::vector@<int@>,
411 * and std::vector@<double@>.
412 */
413 template <typename VectorType>
414 void
415 fill_binary_vector(VectorType &vector) const;
416
417 /**
418 * Output a text representation of this IndexSet to the given stream. Used
419 * for testing.
420 */
421 template <class StreamType>
422 void
423 print(StreamType &out) const;
424
425 /**
426 * Write the IndexSet into a text based file format, that can be read in
427 * again using the read() function.
428 */
429 void
430 write(std::ostream &out) const;
431
432 /**
433 * Construct the IndexSet from a text based representation given by the
434 * stream @p in written by the write() function.
435 */
436 void
437 read(std::istream &in);
438
439 /**
440 * Write the IndexSet into a binary, compact representation, that can be
441 * read in again using the block_read() function.
442 */
443 void
444 block_write(std::ostream &out) const;
445
446 /**
447 * Construct the IndexSet from a binary representation given by the stream
448 * @p in written by the write_block() function.
449 */
450 void
451 block_read(std::istream &in);
452
453 #ifdef DEAL_II_WITH_TRILINOS
454 /**
455 * Given an MPI communicator, create a Trilinos map object that represents a
456 * distribution of vector elements or matrix rows in which we will locally
457 * store those elements or rows for which we store the index in the current
458 * index set, and all the other elements/rows elsewhere on one of the other
459 * MPI processes.
460 *
461 * The last argument only plays a role if the communicator is a parallel
462 * one, distributing computations across multiple processors. In that case,
463 * if the last argument is false, then it is assumed that the index sets
464 * this function is called with on all processors are mutually exclusive but
465 * together enumerate each index exactly once. In other words, if you call
466 * this function on two processors, then the index sets this function is
467 * called with must together have all possible indices from zero to
468 * size()-1, and no index must appear in both index sets. This corresponds,
469 * for example, to the case where we want to split the elements of vectors
470 * into unique subsets to be stored on different processors -- no element
471 * should be owned by more than one processor, but each element must be
472 * owned by one.
473 *
474 * On the other hand, if the second argument is true, then the index sets
475 * can be overlapping, and they also do not need to span the whole index
476 * set. This is a useful operation if we want to create vectors that not
477 * only contain the locally owned indices, but for example also the elements
478 * that correspond to degrees of freedom located on ghost cells. Another
479 * application of this method is to select a subset of the elements of a
480 * vector, e.g. for extracting only certain solution components.
481 */
482 Epetra_Map
483 make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
484 const bool overlapping = false) const;
485
486 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
487 Tpetra::Map<int, types::global_dof_index>
488 make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
489 const bool overlapping = false) const;
490 # endif
491 #endif
492
493
494 /**
495 * Determine an estimate for the memory consumption (in bytes) of this
496 * object.
497 */
498 std::size_t
499 memory_consumption() const;
500
501 DeclException1(ExcIndexNotPresent,
502 size_type,
503 << "The global index " << arg1
504 << " is not an element of this set.");
505
506 /**
507 * Write or read the data of this object to or from a stream for the purpose
508 * of serialization
509 */
510 template <class Archive>
511 void
512 serialize(Archive &ar, const unsigned int version);
513
514
515 /**
516 * @name Iterators
517 * @{
518 */
519
520 /**
521 * Dereferencing an IntervalIterator will return a reference to an object of
522 * this type. It allows access to a contiguous interval $[a,b[$ (also called
523 * a range) of the IndexSet being iterated over.
524 */
525 class IntervalAccessor
526 {
527 public:
528 /**
529 * Construct a valid accessor given an IndexSet and the index @p range_idx
530 * of the range to point to.
531 */
532 IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
533
534 /**
535 * Construct an invalid accessor for the IndexSet.
536 */
537 explicit IntervalAccessor(const IndexSet *idxset);
538
539 /**
540 * Number of elements in this interval.
541 */
542 size_type
543 n_elements() const;
544
545 /**
546 * If true, we are pointing at a valid interval in the IndexSet.
547 */
548 bool
549 is_valid() const;
550
551 /**
552 * Return an iterator pointing at the first index in this interval.
553 */
554 ElementIterator
555 begin() const;
556
557 /**
558 * Return an iterator pointing directly after the last index in this
559 * interval.
560 */
561 ElementIterator
562 end() const;
563
564 /**
565 * Return the index of the last index in this interval.
566 */
567 size_type
568 last() const;
569
570 private:
571 /**
572 * Private copy constructor.
573 */
574 IntervalAccessor(const IntervalAccessor &other);
575 /**
576 * Private copy operator.
577 */
578 IntervalAccessor &
579 operator=(const IntervalAccessor &other);
580
581 /**
582 * Test for equality, used by IntervalIterator.
583 */
584 bool
585 operator==(const IntervalAccessor &other) const;
586 /**
587 * Smaller-than operator, used by IntervalIterator.
588 */
589 bool
590 operator<(const IntervalAccessor &other) const;
591 /**
592 * Advance this accessor to point to the next interval in the @p
593 * index_set.
594 */
595 void
596 advance();
597 /**
598 * Reference to the IndexSet.
599 */
600 const IndexSet *index_set;
601
602 /**
603 * Index into index_set.ranges[]. Set to numbers::invalid_dof_index if
604 * invalid or the end iterator.
605 */
606 size_type range_idx;
607
608 friend class IntervalIterator;
609 };
610
611 /**
612 * Class that represents an iterator pointing to a contiguous interval
613 * $[a,b[$ as returned by IndexSet::begin_interval().
614 */
615 class IntervalIterator
616 {
617 public:
618 /**
619 * Construct a valid iterator pointing to the interval with index @p
620 * range_idx.
621 */
622 IntervalIterator(const IndexSet *idxset, const size_type range_idx);
623
624 /**
625 * Construct an invalid iterator (used as end()).
626 */
627 explicit IntervalIterator(const IndexSet *idxset);
628
629 /**
630 * Construct an empty iterator.
631 */
632 IntervalIterator();
633
634 /**
635 * Copy constructor from @p other iterator.
636 */
637 IntervalIterator(const IntervalIterator &other) = default;
638
639 /**
640 * Assignment of another iterator.
641 */
642 IntervalIterator &
643 operator=(const IntervalIterator &other) = default;
644
645 /**
646 * Prefix increment.
647 */
648 IntervalIterator &
649 operator++();
650
651 /**
652 * Postfix increment.
653 */
654 IntervalIterator
655 operator++(int);
656
657 /**
658 * Dereferencing operator, returns an IntervalAccessor.
659 */
660 const IntervalAccessor &operator*() const;
661
662 /**
663 * Dereferencing operator, returns a pointer to an IntervalAccessor.
664 */
665 const IntervalAccessor *operator->() const;
666
667 /**
668 * Comparison.
669 */
670 bool
671 operator==(const IntervalIterator &) const;
672
673 /**
674 * Inverse of <tt>==</tt>.
675 */
676 bool
677 operator!=(const IntervalIterator &) const;
678
679 /**
680 * Comparison operator.
681 */
682 bool
683 operator<(const IntervalIterator &) const;
684
685 /**
686 * Return the distance between the current iterator and the argument. The
687 * distance is given by how many times one has to apply operator++ to the
688 * current iterator to get the argument (for a positive return value), or
689 * operator-- (for a negative return value).
690 */
691 int
692 operator-(const IntervalIterator &p) const;
693
694 /**
695 * Mark the class as forward iterator and declare some alias which are
696 * standard for iterators and are used by algorithms to enquire about the
697 * specifics of the iterators they work on.
698 */
699 using iterator_category = std::forward_iterator_tag;
700 using value_type = IntervalAccessor;
701 using difference_type = std::ptrdiff_t;
702 using pointer = IntervalAccessor *;
703 using reference = IntervalAccessor &;
704
705 private:
706 /**
707 * Accessor that contains what IndexSet and interval we are pointing at.
708 */
709 IntervalAccessor accessor;
710 };
711
712 /**
713 * Class that represents an iterator pointing to a single element in the
714 * IndexSet as returned by IndexSet::begin().
715 */
716 class ElementIterator
717 {
718 public:
719 /**
720 * Construct an iterator pointing to the global index @p index in the
721 * interval @p range_idx
722 */
723 ElementIterator(const IndexSet *idxset,
724 const size_type range_idx,
725 const size_type index);
726
727 /**
728 * Construct an iterator pointing to the end of the IndexSet.
729 */
730 explicit ElementIterator(const IndexSet *idxset);
731
732 /**
733 * Dereferencing operator. The returned value is the index of the element
734 * inside the IndexSet.
735 */
736 size_type operator*() const;
737
738 /**
739 * Does this iterator point to an existing element?
740 */
741 bool
742 is_valid() const;
743
744 /**
745 * Prefix increment.
746 */
747 ElementIterator &
748 operator++();
749
750 /**
751 * Postfix increment.
752 */
753 ElementIterator
754 operator++(int);
755
756 /**
757 * Comparison.
758 */
759 bool
760 operator==(const ElementIterator &) const;
761
762 /**
763 * Inverse of <tt>==</tt>.
764 */
765 bool
766 operator!=(const ElementIterator &) const;
767
768 /**
769 * Comparison operator.
770 */
771 bool
772 operator<(const ElementIterator &) const;
773
774 /**
775 * Return the distance between the current iterator and the argument. In
776 * the expression <code>it_left-it_right</code> the distance is given by
777 * how many times one has to apply operator++ to the right operand @p
778 * it_right to get the left operand @p it_left (for a positive return
779 * value), or to @p it_left to get the @p it_right (for a negative return
780 * value).
781 */
782 std::ptrdiff_t
783 operator-(const ElementIterator &p) const;
784
785 /**
786 * Mark the class as forward iterator and declare some alias which are
787 * standard for iterators and are used by algorithms to enquire about the
788 * specifics of the iterators they work on.
789 */
790 using iterator_category = std::forward_iterator_tag;
791 using value_type = size_type;
792 using difference_type = std::ptrdiff_t;
793 using pointer = size_type *;
794 using reference = size_type &;
795
796 private:
797 /**
798 * Advance iterator by one.
799 */
800 void
801 advance();
802
803 /**
804 * The parent IndexSet.
805 */
806 const IndexSet *index_set;
807 /**
808 * Index into index_set.ranges.
809 */
810 size_type range_idx;
811 /**
812 * The global index this iterator is pointing at.
813 */
814 size_type idx;
815 };
816
817 /**
818 * Return an iterator that points at the first index that is contained in
819 * this IndexSet.
820 */
821 ElementIterator
822 begin() const;
823
824 /**
825 * Return an element iterator pointing to the element with global index
826 * @p global_index or the next larger element if the index is not in the
827 * set. This is equivalent to
828 * @code
829 * auto p = begin();
830 * while (*p<global_index)
831 * ++p;
832 * return p;
833 * @endcode
834 *
835 * If there is no element in this IndexSet at or behind @p global_index,
836 * this method will return end().
837 */
838 ElementIterator
839 at(const size_type global_index) const;
840
841 /**
842 * Return an iterator that points one after the last index that is contained
843 * in this IndexSet.
844 */
845 ElementIterator
846 end() const;
847
848 /**
849 * Return an Iterator that points at the first interval of this IndexSet.
850 */
851 IntervalIterator
852 begin_intervals() const;
853
854 /**
855 * Return an Iterator that points one after the last interval of this
856 * IndexSet.
857 */
858 IntervalIterator
859 end_intervals() const;
860
861 /**
862 * @}
863 */
864
865 private:
866 /**
867 * A type that denotes the half open index range <code>[begin,end)</code>.
868 *
869 * The nth_index_in_set denotes the how many-th index within this IndexSet
870 * the first element of the current range is. This information is only
871 * accurate if IndexSet::compress() has been called after the last
872 * insertion.
873 */
874 struct Range
875 {
876 size_type begin;
877 size_type end;
878
879 size_type nth_index_in_set;
880
881 /**
882 * Default constructor. Since there is no useful choice for a default
883 * constructed interval, this constructor simply creates something that
884 * resembles an invalid range. We need this constructor for serialization
885 * purposes, but the invalid range should be filled with something read
886 * from the archive before it is used, so we should hopefully never get to
887 * see an invalid range in the wild.
888 */
889 Range();
890
891 /**
892 * Constructor. Create a half-open interval with the given indices.
893 *
894 * @param i1 Left end point of the interval.
895 * @param i2 First index greater than the last index of the indicated
896 * range.
897 */
898 Range(const size_type i1, const size_type i2);
899
900 friend inline bool
901 operator<(const Range &range_1, const Range &range_2)
902 {
903 return (
904 (range_1.begin < range_2.begin) ||
905 ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
906 }
907
908 static bool
end_compareRange909 end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
910 {
911 return x.end < y.end;
912 }
913
914 static bool
nth_index_compareRange915 nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
916 {
917 return (x.nth_index_in_set + (x.end - x.begin) <
918 y.nth_index_in_set + (y.end - y.begin));
919 }
920
921 friend inline bool
922 operator==(const Range &range_1, const Range &range_2)
923 {
924 return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
925 }
926
927 static std::size_t
memory_consumptionRange928 memory_consumption()
929 {
930 return sizeof(Range);
931 }
932
933 /**
934 * Write or read the data of this object to or from a stream for the
935 * purpose of serialization
936 */
937 template <class Archive>
938 void
939 serialize(Archive &ar, const unsigned int version);
940 };
941
942 /**
943 * A set of contiguous ranges of indices that make up (part of) this index
944 * set. This variable is always kept sorted.
945 *
946 * The variable is marked "mutable" so that it can be changed by compress(),
947 * though this of course doesn't change anything about the external
948 * representation of this index set.
949 */
950 mutable std::vector<Range> ranges;
951
952 /**
953 * True if compress() has been called after the last change in the set of
954 * indices.
955 *
956 * The variable is marked "mutable" so that it can be changed by compress(),
957 * though this of course doesn't change anything about the external
958 * representation of this index set.
959 */
960 mutable bool is_compressed;
961
962 /**
963 * The overall size of the index range. Elements of this index set have to
964 * have a smaller number than this value.
965 */
966 size_type index_space_size;
967
968 /**
969 * This integer caches the index of the largest range in @p ranges. This
970 * gives <tt>O(1)</tt> access to the range with most elements, while general
971 * access costs <tt>O(log(n_ranges))</tt>. The largest range is needed for
972 * the methods @p is_element(), @p index_within_set(), @p nth_index_in_set.
973 * In many applications, the largest range contains most elements (the
974 * locally owned range), whereas there are only a few other elements
975 * (ghosts).
976 */
977 mutable size_type largest_range;
978
979 /**
980 * A mutex that is used to synchronize operations of the do_compress()
981 * function that is called from many 'const' functions via compress().
982 */
983 mutable Threads::Mutex compress_mutex;
984
985 /**
986 * Actually perform the compress() operation.
987 */
988 void
989 do_compress() const;
990 };
991
992
993 /**
994 * Create and return an index set of size $N$ that contains every single index
995 * within this range. In essence, this function returns an index set created
996 * by
997 * @code
998 * IndexSet is (N);
999 * is.add_range(0, N);
1000 * @endcode
1001 * This function exists so that one can create and initialize index sets that
1002 * are complete in one step, or so one can write code like
1003 * @code
1004 * if (my_index_set == complete_index_set(my_index_set.size())
1005 * ...
1006 * @endcode
1007 *
1008 * @relatesalso IndexSet
1009 */
1010 inline IndexSet
complete_index_set(const IndexSet::size_type N)1011 complete_index_set(const IndexSet::size_type N)
1012 {
1013 IndexSet is(N);
1014 is.add_range(0, N);
1015 is.compress();
1016 return is;
1017 }
1018
1019 /* ------------------ inline functions ------------------ */
1020
1021
1022 /* IntervalAccessor */
1023
IntervalAccessor(const IndexSet * idxset,const IndexSet::size_type range_idx)1024 inline IndexSet::IntervalAccessor::IntervalAccessor(
1025 const IndexSet * idxset,
1026 const IndexSet::size_type range_idx)
1027 : index_set(idxset)
1028 , range_idx(range_idx)
1029 {
1030 Assert(range_idx < idxset->n_intervals(),
1031 ExcInternalError("Invalid range index"));
1032 }
1033
1034
1035
IntervalAccessor(const IndexSet * idxset)1036 inline IndexSet::IntervalAccessor::IntervalAccessor(const IndexSet *idxset)
1037 : index_set(idxset)
1038 , range_idx(numbers::invalid_dof_index)
1039 {}
1040
1041
1042
IntervalAccessor(const IndexSet::IntervalAccessor & other)1043 inline IndexSet::IntervalAccessor::IntervalAccessor(
1044 const IndexSet::IntervalAccessor &other)
1045 : index_set(other.index_set)
1046 , range_idx(other.range_idx)
1047 {
1048 Assert(range_idx == numbers::invalid_dof_index || is_valid(),
1049 ExcMessage("invalid iterator"));
1050 }
1051
1052
1053
1054 inline IndexSet::size_type
n_elements()1055 IndexSet::IntervalAccessor::n_elements() const
1056 {
1057 Assert(is_valid(), ExcMessage("invalid iterator"));
1058 return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1059 }
1060
1061
1062
1063 inline bool
is_valid()1064 IndexSet::IntervalAccessor::is_valid() const
1065 {
1066 return index_set != nullptr && range_idx < index_set->n_intervals();
1067 }
1068
1069
1070
1071 inline IndexSet::ElementIterator
begin()1072 IndexSet::IntervalAccessor::begin() const
1073 {
1074 Assert(is_valid(), ExcMessage("invalid iterator"));
1075 return {index_set, range_idx, index_set->ranges[range_idx].begin};
1076 }
1077
1078
1079
1080 inline IndexSet::ElementIterator
end()1081 IndexSet::IntervalAccessor::end() const
1082 {
1083 Assert(is_valid(), ExcMessage("invalid iterator"));
1084
1085 // point to first index in next interval unless we are the last interval.
1086 if (range_idx < index_set->ranges.size() - 1)
1087 return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1088 else
1089 return index_set->end();
1090 }
1091
1092
1093
1094 inline IndexSet::size_type
last()1095 IndexSet::IntervalAccessor::last() const
1096 {
1097 Assert(is_valid(), ExcMessage("invalid iterator"));
1098
1099 return index_set->ranges[range_idx].end - 1;
1100 }
1101
1102
1103
1104 inline IndexSet::IntervalAccessor &
1105 IndexSet::IntervalAccessor::operator=(const IndexSet::IntervalAccessor &other)
1106 {
1107 index_set = other.index_set;
1108 range_idx = other.range_idx;
1109 Assert(range_idx == numbers::invalid_dof_index || is_valid(),
1110 ExcMessage("invalid iterator"));
1111 return *this;
1112 }
1113
1114
1115
1116 inline bool
1117 IndexSet::IntervalAccessor::
1118 operator==(const IndexSet::IntervalAccessor &other) const
1119 {
1120 Assert(index_set == other.index_set,
1121 ExcMessage(
1122 "Can not compare accessors pointing to different IndexSets"));
1123 return range_idx == other.range_idx;
1124 }
1125
1126
1127
1128 inline bool
1129 IndexSet::IntervalAccessor::
1130 operator<(const IndexSet::IntervalAccessor &other) const
1131 {
1132 Assert(index_set == other.index_set,
1133 ExcMessage(
1134 "Can not compare accessors pointing to different IndexSets"));
1135 return range_idx < other.range_idx;
1136 }
1137
1138
1139
1140 inline void
advance()1141 IndexSet::IntervalAccessor::advance()
1142 {
1143 Assert(
1144 is_valid(),
1145 ExcMessage(
1146 "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1147 ++range_idx;
1148
1149 // set ourselves to invalid if we walk off the end
1150 if (range_idx >= index_set->ranges.size())
1151 range_idx = numbers::invalid_dof_index;
1152 }
1153
1154
1155 /* IntervalIterator */
1156
IntervalIterator(const IndexSet * idxset,const IndexSet::size_type range_idx)1157 inline IndexSet::IntervalIterator::IntervalIterator(
1158 const IndexSet * idxset,
1159 const IndexSet::size_type range_idx)
1160 : accessor(idxset, range_idx)
1161 {}
1162
1163
1164
IntervalIterator()1165 inline IndexSet::IntervalIterator::IntervalIterator()
1166 : accessor(nullptr)
1167 {}
1168
1169
1170
IntervalIterator(const IndexSet * idxset)1171 inline IndexSet::IntervalIterator::IntervalIterator(const IndexSet *idxset)
1172 : accessor(idxset)
1173 {}
1174
1175
1176
1177 inline IndexSet::IntervalIterator &
1178 IndexSet::IntervalIterator::operator++()
1179 {
1180 accessor.advance();
1181 return *this;
1182 }
1183
1184
1185
1186 inline IndexSet::IntervalIterator
1187 IndexSet::IntervalIterator::operator++(int)
1188 {
1189 const IndexSet::IntervalIterator iter = *this;
1190 accessor.advance();
1191 return iter;
1192 }
1193
1194
1195
1196 inline const IndexSet::IntervalAccessor &IndexSet::IntervalIterator::
1197 operator*() const
1198 {
1199 return accessor;
1200 }
1201
1202
1203
1204 inline const IndexSet::IntervalAccessor *IndexSet::IntervalIterator::
1205 operator->() const
1206 {
1207 return &accessor;
1208 }
1209
1210
1211
1212 inline bool
1213 IndexSet::IntervalIterator::
1214 operator==(const IndexSet::IntervalIterator &other) const
1215 {
1216 return accessor == other.accessor;
1217 }
1218
1219
1220
1221 inline bool
1222 IndexSet::IntervalIterator::
1223 operator!=(const IndexSet::IntervalIterator &other) const
1224 {
1225 return !(*this == other);
1226 }
1227
1228
1229
1230 inline bool
1231 IndexSet::IntervalIterator::
1232 operator<(const IndexSet::IntervalIterator &other) const
1233 {
1234 return accessor < other.accessor;
1235 }
1236
1237
1238
1239 inline int
1240 IndexSet::IntervalIterator::
1241 operator-(const IndexSet::IntervalIterator &other) const
1242 {
1243 Assert(accessor.index_set == other.accessor.index_set,
1244 ExcMessage(
1245 "Can not compare iterators belonging to different IndexSets"));
1246
1247 const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index) ?
1248 accessor.index_set->ranges.size() :
1249 accessor.range_idx;
1250 const size_type rhs =
1251 (other.accessor.range_idx == numbers::invalid_dof_index) ?
1252 accessor.index_set->ranges.size() :
1253 other.accessor.range_idx;
1254
1255 if (lhs > rhs)
1256 return static_cast<int>(lhs - rhs);
1257 else
1258 return -static_cast<int>(rhs - lhs);
1259 }
1260
1261
1262
1263 /* ElementIterator */
1264
ElementIterator(const IndexSet * idxset,const IndexSet::size_type range_idx,const IndexSet::size_type index)1265 inline IndexSet::ElementIterator::ElementIterator(
1266 const IndexSet * idxset,
1267 const IndexSet::size_type range_idx,
1268 const IndexSet::size_type index)
1269 : index_set(idxset)
1270 , range_idx(range_idx)
1271 , idx(index)
1272 {
1273 Assert(range_idx < index_set->ranges.size(),
1274 ExcMessage(
1275 "Invalid range index for IndexSet::ElementIterator constructor."));
1276 Assert(
1277 idx >= index_set->ranges[range_idx].begin &&
1278 idx < index_set->ranges[range_idx].end,
1279 ExcInternalError(
1280 "Invalid index argument for IndexSet::ElementIterator constructor."));
1281 }
1282
1283
1284
ElementIterator(const IndexSet * idxset)1285 inline IndexSet::ElementIterator::ElementIterator(const IndexSet *idxset)
1286 : index_set(idxset)
1287 , range_idx(numbers::invalid_dof_index)
1288 , idx(numbers::invalid_dof_index)
1289 {}
1290
1291
1292
1293 inline bool
is_valid()1294 IndexSet::ElementIterator::is_valid() const
1295 {
1296 Assert((range_idx == numbers::invalid_dof_index &&
1297 idx == numbers::invalid_dof_index) ||
1298 (range_idx < index_set->ranges.size() &&
1299 idx < index_set->ranges[range_idx].end),
1300 ExcInternalError("Invalid ElementIterator state."));
1301
1302 return (range_idx < index_set->ranges.size() &&
1303 idx < index_set->ranges[range_idx].end);
1304 }
1305
1306
1307
1308 inline IndexSet::size_type IndexSet::ElementIterator::operator*() const
1309 {
1310 Assert(
1311 is_valid(),
1312 ExcMessage(
1313 "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1314 return idx;
1315 }
1316
1317
1318
1319 inline bool
1320 IndexSet::ElementIterator::
1321 operator==(const IndexSet::ElementIterator &other) const
1322 {
1323 Assert(index_set == other.index_set,
1324 ExcMessage(
1325 "Can not compare iterators belonging to different IndexSets"));
1326 return range_idx == other.range_idx && idx == other.idx;
1327 }
1328
1329
1330
1331 inline void
advance()1332 IndexSet::ElementIterator::advance()
1333 {
1334 Assert(
1335 is_valid(),
1336 ExcMessage(
1337 "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1338 if (idx < index_set->ranges[range_idx].end)
1339 ++idx;
1340 // end of this range?
1341 if (idx == index_set->ranges[range_idx].end)
1342 {
1343 // point to first element in next interval if possible
1344 if (range_idx < index_set->ranges.size() - 1)
1345 {
1346 ++range_idx;
1347 idx = index_set->ranges[range_idx].begin;
1348 }
1349 else
1350 {
1351 // we just fell off the end, set to invalid:
1352 range_idx = numbers::invalid_dof_index;
1353 idx = numbers::invalid_dof_index;
1354 }
1355 }
1356 }
1357
1358
1359
1360 inline IndexSet::ElementIterator &
1361 IndexSet::ElementIterator::operator++()
1362 {
1363 advance();
1364 return *this;
1365 }
1366
1367
1368
1369 inline IndexSet::ElementIterator
1370 IndexSet::ElementIterator::operator++(int)
1371 {
1372 const IndexSet::ElementIterator it = *this;
1373 advance();
1374 return it;
1375 }
1376
1377
1378
1379 inline bool
1380 IndexSet::ElementIterator::
1381 operator!=(const IndexSet::ElementIterator &other) const
1382 {
1383 return !(*this == other);
1384 }
1385
1386
1387
1388 inline bool
1389 IndexSet::ElementIterator::
1390 operator<(const IndexSet::ElementIterator &other) const
1391 {
1392 Assert(index_set == other.index_set,
1393 ExcMessage(
1394 "Can not compare iterators belonging to different IndexSets"));
1395 return range_idx < other.range_idx ||
1396 (range_idx == other.range_idx && idx < other.idx);
1397 }
1398
1399
1400
1401 inline std::ptrdiff_t
1402 IndexSet::ElementIterator::
1403 operator-(const IndexSet::ElementIterator &other) const
1404 {
1405 Assert(index_set == other.index_set,
1406 ExcMessage(
1407 "Can not compare iterators belonging to different IndexSets"));
1408 if (*this == other)
1409 return 0;
1410 if (!(*this < other))
1411 return -(other - *this);
1412
1413 // only other can be equal to end() because of the checks above.
1414 Assert(is_valid(), ExcInternalError());
1415
1416 // Note: we now compute how far advance *this in "*this < other" to get other,
1417 // so we need to return -c at the end.
1418
1419 // first finish the current range:
1420 std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1421
1422 // now walk in steps of ranges (need to start one behind our current one):
1423 for (size_type range = range_idx + 1;
1424 range < index_set->ranges.size() && range <= other.range_idx;
1425 ++range)
1426 c += index_set->ranges[range].end - index_set->ranges[range].begin;
1427
1428 Assert(
1429 other.range_idx < index_set->ranges.size() ||
1430 other.range_idx == numbers::invalid_dof_index,
1431 ExcMessage(
1432 "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1433
1434 // We might have walked too far because we went until the end of
1435 // other.range_idx, so walk backwards to other.idx:
1436 if (other.range_idx != numbers::invalid_dof_index)
1437 c -= index_set->ranges[other.range_idx].end - other.idx;
1438
1439 return -c;
1440 }
1441
1442
1443 /* Range */
1444
Range()1445 inline IndexSet::Range::Range()
1446 : begin(numbers::invalid_dof_index)
1447 , end(numbers::invalid_dof_index)
1448 , nth_index_in_set(numbers::invalid_dof_index)
1449 {}
1450
1451
1452
Range(const size_type i1,const size_type i2)1453 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1454 : begin(i1)
1455 , end(i2)
1456 , nth_index_in_set(numbers::invalid_dof_index)
1457 {}
1458
1459
1460
1461 /* IndexSet itself */
1462
IndexSet()1463 inline IndexSet::IndexSet()
1464 : is_compressed(true)
1465 , index_space_size(0)
1466 , largest_range(numbers::invalid_unsigned_int)
1467 {}
1468
1469
1470
IndexSet(const size_type size)1471 inline IndexSet::IndexSet(const size_type size)
1472 : is_compressed(true)
1473 , index_space_size(size)
1474 , largest_range(numbers::invalid_unsigned_int)
1475 {}
1476
1477
1478
IndexSet(IndexSet && is)1479 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1480 : ranges(std::move(is.ranges))
1481 , is_compressed(is.is_compressed)
1482 , index_space_size(is.index_space_size)
1483 , largest_range(is.largest_range)
1484 {
1485 is.ranges.clear();
1486 is.is_compressed = true;
1487 is.index_space_size = 0;
1488 is.largest_range = numbers::invalid_unsigned_int;
1489
1490 compress();
1491 }
1492
1493
1494
1495 inline IndexSet &
1496 IndexSet::operator=(IndexSet &&is) noexcept
1497 {
1498 ranges = std::move(is.ranges);
1499 is_compressed = is.is_compressed;
1500 index_space_size = is.index_space_size;
1501 largest_range = is.largest_range;
1502
1503 is.ranges.clear();
1504 is.is_compressed = true;
1505 is.index_space_size = 0;
1506 is.largest_range = numbers::invalid_unsigned_int;
1507
1508 compress();
1509
1510 return *this;
1511 }
1512
1513
1514
1515 inline IndexSet::ElementIterator
begin()1516 IndexSet::begin() const
1517 {
1518 compress();
1519 if (ranges.size() > 0)
1520 return {this, 0, ranges[0].begin};
1521 else
1522 return end();
1523 }
1524
1525
1526
1527 inline IndexSet::ElementIterator
at(const size_type global_index)1528 IndexSet::at(const size_type global_index) const
1529 {
1530 compress();
1531 AssertIndexRange(global_index, size());
1532
1533 if (ranges.empty())
1534 return end();
1535
1536 std::vector<Range>::const_iterator main_range =
1537 ranges.begin() + largest_range;
1538
1539 Range r(global_index, global_index + 1);
1540 // This optimization makes the bounds for lower_bound smaller by checking
1541 // the largest range first.
1542 std::vector<Range>::const_iterator range_begin, range_end;
1543 if (global_index < main_range->begin)
1544 {
1545 range_begin = ranges.begin();
1546 range_end = main_range;
1547 }
1548 else
1549 {
1550 range_begin = main_range;
1551 range_end = ranges.end();
1552 }
1553
1554 // This will give us the first range p=[a,b[ with b>=global_index using
1555 // a binary search
1556 const std::vector<Range>::const_iterator p =
1557 Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1558
1559 // We couldn't find a range, which means we have no range that contains
1560 // global_index and also no range behind it, meaning we need to return end().
1561 if (p == ranges.end())
1562 return end();
1563
1564 // Finally, we can have two cases: Either global_index is not in [a,b[,
1565 // which means we need to return an iterator to a because global_index, ...,
1566 // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
1567 // [a,b[ and we will return an iterator pointing directly at global_index
1568 // (else branch).
1569 if (global_index < p->begin)
1570 return {this, static_cast<size_type>(p - ranges.begin()), p->begin};
1571 else
1572 return {this, static_cast<size_type>(p - ranges.begin()), global_index};
1573 }
1574
1575
1576
1577 inline IndexSet::ElementIterator
end()1578 IndexSet::end() const
1579 {
1580 compress();
1581 return IndexSet::ElementIterator(this);
1582 }
1583
1584
1585
1586 inline IndexSet::IntervalIterator
begin_intervals()1587 IndexSet::begin_intervals() const
1588 {
1589 compress();
1590 if (ranges.size() > 0)
1591 return IndexSet::IntervalIterator(this, 0);
1592 else
1593 return end_intervals();
1594 }
1595
1596
1597
1598 inline IndexSet::IntervalIterator
end_intervals()1599 IndexSet::end_intervals() const
1600 {
1601 compress();
1602 return IndexSet::IntervalIterator(this);
1603 }
1604
1605
1606
1607 inline void
clear()1608 IndexSet::clear()
1609 {
1610 // reset so that there are no indices in the set any more; however,
1611 // as documented, the index set retains its size
1612 ranges.clear();
1613 is_compressed = true;
1614 largest_range = numbers::invalid_unsigned_int;
1615 }
1616
1617
1618
1619 inline void
set_size(const size_type sz)1620 IndexSet::set_size(const size_type sz)
1621 {
1622 Assert(ranges.empty(),
1623 ExcMessage("This function can only be called if the current "
1624 "object does not yet contain any elements."));
1625 index_space_size = sz;
1626 is_compressed = true;
1627 }
1628
1629
1630
1631 inline IndexSet::size_type
size()1632 IndexSet::size() const
1633 {
1634 return index_space_size;
1635 }
1636
1637
1638
1639 inline void
compress()1640 IndexSet::compress() const
1641 {
1642 if (is_compressed == true)
1643 return;
1644
1645 do_compress();
1646 }
1647
1648
1649
1650 inline void
add_index(const size_type index)1651 IndexSet::add_index(const size_type index)
1652 {
1653 AssertIndexRange(index, index_space_size);
1654
1655 const Range new_range(index, index + 1);
1656 if (ranges.size() == 0 || index > ranges.back().end)
1657 ranges.push_back(new_range);
1658 else if (index == ranges.back().end)
1659 ranges.back().end++;
1660 else
1661 ranges.insert(Utilities::lower_bound(ranges.begin(),
1662 ranges.end(),
1663 new_range),
1664 new_range);
1665 is_compressed = false;
1666 }
1667
1668
1669
1670 inline void
add_range(const size_type begin,const size_type end)1671 IndexSet::add_range(const size_type begin, const size_type end)
1672 {
1673 Assert((begin < index_space_size) ||
1674 ((begin == index_space_size) && (end == index_space_size)),
1675 ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1676 Assert(end <= index_space_size,
1677 ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1678 AssertIndexRange(begin, end + 1);
1679
1680 if (begin != end)
1681 {
1682 const Range new_range(begin, end);
1683
1684 // the new index might be larger than the last index present in the
1685 // ranges. Then we can skip the binary search
1686 if (ranges.size() == 0 || begin > ranges.back().end)
1687 ranges.push_back(new_range);
1688 else
1689 ranges.insert(Utilities::lower_bound(ranges.begin(),
1690 ranges.end(),
1691 new_range),
1692 new_range);
1693 is_compressed = false;
1694 }
1695 }
1696
1697
1698
1699 template <typename ForwardIterator>
1700 inline void
add_indices(const ForwardIterator & begin,const ForwardIterator & end)1701 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1702 {
1703 if (begin == end)
1704 return;
1705
1706 // identify ranges in the given iterator range by checking whether some
1707 // indices happen to be consecutive. to avoid quadratic complexity when
1708 // calling add_range many times (as add_range() going into the middle of an
1709 // already existing range must shift entries around), we first collect a
1710 // vector of ranges.
1711 std::vector<std::pair<size_type, size_type>> tmp_ranges;
1712 bool ranges_are_sorted = true;
1713 for (ForwardIterator p = begin; p != end;)
1714 {
1715 const size_type begin_index = *p;
1716 size_type end_index = begin_index + 1;
1717 ForwardIterator q = p;
1718 ++q;
1719 while ((q != end) && (*q == end_index))
1720 {
1721 ++end_index;
1722 ++q;
1723 }
1724
1725 tmp_ranges.emplace_back(begin_index, end_index);
1726 p = q;
1727
1728 // if the starting index of the next go-around of the for loop is less
1729 // than the end index of the one just identified, then we will have at
1730 // least one pair of ranges that are not sorted, and consequently the
1731 // whole collection of ranges is not sorted.
1732 if (p != end && *p < end_index)
1733 ranges_are_sorted = false;
1734 }
1735
1736 if (!ranges_are_sorted)
1737 std::sort(tmp_ranges.begin(), tmp_ranges.end());
1738
1739 // if we have many ranges, we first construct a temporary index set (where
1740 // we add ranges in a consecutive way, so fast), otherwise, we work with
1741 // add_range(). the number 9 is chosen heuristically given the fact that
1742 // there are typically up to 8 independent ranges when adding the degrees of
1743 // freedom on a 3D cell or 9 when adding degrees of freedom of faces. if
1744 // doing cell-by-cell additions, we want to avoid repeated calls to
1745 // IndexSet::compress() which gets called upon merging two index sets, so we
1746 // want to be in the other branch then.
1747 if (tmp_ranges.size() > 9)
1748 {
1749 IndexSet tmp_set(size());
1750 tmp_set.ranges.reserve(tmp_ranges.size());
1751 for (const auto &i : tmp_ranges)
1752 tmp_set.add_range(i.first, i.second);
1753 this->add_indices(tmp_set);
1754 }
1755 else
1756 for (const auto &i : tmp_ranges)
1757 add_range(i.first, i.second);
1758 }
1759
1760
1761
1762 inline bool
is_element(const size_type index)1763 IndexSet::is_element(const size_type index) const
1764 {
1765 if (ranges.empty() == false)
1766 {
1767 compress();
1768
1769 // fast check whether the index is in the largest range
1770 Assert(largest_range < ranges.size(), ExcInternalError());
1771 if (index >= ranges[largest_range].begin &&
1772 index < ranges[largest_range].end)
1773 return true;
1774
1775 // get the element after which we would have to insert a range that
1776 // consists of all elements from this element to the end of the index
1777 // range plus one. after this call we know that if p!=end() then
1778 // p->begin<=index unless there is no such range at all
1779 //
1780 // if the searched for element is an element of this range, then we're
1781 // done. otherwise, the element can't be in one of the following ranges
1782 // because otherwise p would be a different iterator
1783 //
1784 // since we already know the position relative to the largest range (we
1785 // called compress!), we can perform the binary search on ranges with
1786 // lower/higher number compared to the largest range
1787 std::vector<Range>::const_iterator p = std::upper_bound(
1788 ranges.begin() +
1789 (index < ranges[largest_range].begin ? 0 : largest_range + 1),
1790 index < ranges[largest_range].begin ? ranges.begin() + largest_range :
1791 ranges.end(),
1792 Range(index, size() + 1));
1793
1794 if (p == ranges.begin())
1795 return ((index >= p->begin) && (index < p->end));
1796
1797 Assert((p == ranges.end()) || (p->begin > index), ExcInternalError());
1798
1799 // now move to that previous range
1800 --p;
1801 Assert(p->begin <= index, ExcInternalError());
1802
1803 return (p->end > index);
1804 }
1805
1806 // didn't find this index, so it's not in the set
1807 return false;
1808 }
1809
1810
1811
1812 inline bool
is_contiguous()1813 IndexSet::is_contiguous() const
1814 {
1815 compress();
1816 return (ranges.size() <= 1);
1817 }
1818
1819
1820
1821 inline bool
is_empty()1822 IndexSet::is_empty() const
1823 {
1824 return ranges.empty();
1825 }
1826
1827
1828
1829 inline IndexSet::size_type
n_elements()1830 IndexSet::n_elements() const
1831 {
1832 // make sure we have non-overlapping ranges
1833 compress();
1834
1835 size_type v = 0;
1836 if (!ranges.empty())
1837 {
1838 Range &r = ranges.back();
1839 v = r.nth_index_in_set + r.end - r.begin;
1840 }
1841
1842 #ifdef DEBUG
1843 size_type s = 0;
1844 for (const auto &range : ranges)
1845 s += (range.end - range.begin);
1846 Assert(s == v, ExcInternalError());
1847 #endif
1848
1849 return v;
1850 }
1851
1852
1853
1854 inline unsigned int
n_intervals()1855 IndexSet::n_intervals() const
1856 {
1857 compress();
1858 return ranges.size();
1859 }
1860
1861
1862
1863 inline IndexSet::size_type
largest_range_starting_index()1864 IndexSet::largest_range_starting_index() const
1865 {
1866 Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1867
1868 compress();
1869 const std::vector<Range>::const_iterator main_range =
1870 ranges.begin() + largest_range;
1871
1872 return main_range->nth_index_in_set;
1873 }
1874
1875
1876
1877 inline IndexSet::size_type
nth_index_in_set(const size_type n)1878 IndexSet::nth_index_in_set(const size_type n) const
1879 {
1880 AssertIndexRange(n, n_elements());
1881
1882 compress();
1883
1884 // first check whether the index is in the largest range
1885 Assert(largest_range < ranges.size(), ExcInternalError());
1886 std::vector<Range>::const_iterator main_range =
1887 ranges.begin() + largest_range;
1888 if (n >= main_range->nth_index_in_set &&
1889 n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1890 return main_range->begin + (n - main_range->nth_index_in_set);
1891
1892 // find out which chunk the local index n belongs to by using a binary
1893 // search. the comparator is based on the end of the ranges. Use the
1894 // position relative to main_range to subdivide the ranges
1895 Range r(n, n + 1);
1896 r.nth_index_in_set = n;
1897 std::vector<Range>::const_iterator range_begin, range_end;
1898 if (n < main_range->nth_index_in_set)
1899 {
1900 range_begin = ranges.begin();
1901 range_end = main_range;
1902 }
1903 else
1904 {
1905 range_begin = main_range + 1;
1906 range_end = ranges.end();
1907 }
1908
1909 const std::vector<Range>::const_iterator p =
1910 Utilities::lower_bound(range_begin, range_end, r, Range::nth_index_compare);
1911
1912 Assert(p != ranges.end(), ExcInternalError());
1913 return p->begin + (n - p->nth_index_in_set);
1914 }
1915
1916
1917
1918 inline IndexSet::size_type
index_within_set(const size_type n)1919 IndexSet::index_within_set(const size_type n) const
1920 {
1921 // to make this call thread-safe, compress() must not be called through this
1922 // function
1923 Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1924 AssertIndexRange(n, size());
1925
1926 // return immediately if the index set is empty
1927 if (is_empty())
1928 return numbers::invalid_dof_index;
1929
1930 // check whether the index is in the largest range. use the result to
1931 // perform a one-sided binary search afterward
1932 Assert(largest_range < ranges.size(), ExcInternalError());
1933 std::vector<Range>::const_iterator main_range =
1934 ranges.begin() + largest_range;
1935 if (n >= main_range->begin && n < main_range->end)
1936 return (n - main_range->begin) + main_range->nth_index_in_set;
1937
1938 Range r(n, n);
1939 std::vector<Range>::const_iterator range_begin, range_end;
1940 if (n < main_range->begin)
1941 {
1942 range_begin = ranges.begin();
1943 range_end = main_range;
1944 }
1945 else
1946 {
1947 range_begin = main_range + 1;
1948 range_end = ranges.end();
1949 }
1950
1951 std::vector<Range>::const_iterator p =
1952 Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1953
1954 // if n is not in this set
1955 if (p == range_end || p->end == n || p->begin > n)
1956 return numbers::invalid_dof_index;
1957
1958 Assert(p != ranges.end(), ExcInternalError());
1959 Assert(p->begin <= n, ExcInternalError());
1960 Assert(n < p->end, ExcInternalError());
1961 return (n - p->begin) + p->nth_index_in_set;
1962 }
1963
1964
1965
1966 inline bool
1967 IndexSet::operator==(const IndexSet &is) const
1968 {
1969 Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1970
1971 compress();
1972 is.compress();
1973
1974 return ranges == is.ranges;
1975 }
1976
1977
1978
1979 inline bool
1980 IndexSet::operator!=(const IndexSet &is) const
1981 {
1982 Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1983
1984 compress();
1985 is.compress();
1986
1987 return ranges != is.ranges;
1988 }
1989
1990
1991
1992 template <typename Vector>
1993 void
fill_binary_vector(Vector & vector)1994 IndexSet::fill_binary_vector(Vector &vector) const
1995 {
1996 Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
1997
1998 compress();
1999 // first fill all elements of the vector with zeroes.
2000 std::fill(vector.begin(), vector.end(), 0);
2001
2002 // then write ones into the elements whose indices are contained in the
2003 // index set
2004 for (const auto &range : ranges)
2005 for (size_type i = range.begin; i < range.end; ++i)
2006 vector[i] = 1;
2007 }
2008
2009
2010
2011 template <class StreamType>
2012 inline void
print(StreamType & out)2013 IndexSet::print(StreamType &out) const
2014 {
2015 compress();
2016 out << "{";
2017 std::vector<Range>::const_iterator p;
2018 for (p = ranges.begin(); p != ranges.end(); ++p)
2019 {
2020 if (p->end - p->begin == 1)
2021 out << p->begin;
2022 else
2023 out << "[" << p->begin << "," << p->end - 1 << "]";
2024
2025 if (p != --ranges.end())
2026 out << ", ";
2027 }
2028 out << "}" << std::endl;
2029 }
2030
2031
2032
2033 template <class Archive>
2034 inline void
serialize(Archive & ar,const unsigned int)2035 IndexSet::Range::serialize(Archive &ar, const unsigned int)
2036 {
2037 ar &begin &end &nth_index_in_set;
2038 }
2039
2040
2041
2042 template <class Archive>
2043 inline void
serialize(Archive & ar,const unsigned int)2044 IndexSet::serialize(Archive &ar, const unsigned int)
2045 {
2046 ar &ranges &is_compressed &index_space_size &largest_range;
2047 }
2048
2049 DEAL_II_NAMESPACE_CLOSE
2050
2051 #endif
2052