Reference documentation for deal.II version GIT relicensing-1291-gf3a46e514e 2024-07-24 11:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
index_set.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2009 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
17#include <deal.II/base/mpi.h>
18
20
21#include <vector>
22
23#ifdef DEAL_II_WITH_TRILINOS
24# ifdef DEAL_II_WITH_MPI
25# include <Epetra_MpiComm.h>
26# endif
27# include <Epetra_Map.h>
28# include <Epetra_SerialComm.h>
29# ifdef DEAL_II_TRILINOS_WITH_TPETRA
30# include <Tpetra_Map.hpp>
31# endif
32#endif
33
35
36
37
38#ifdef DEAL_II_WITH_TRILINOS
39
40# ifdef DEAL_II_TRILINOS_WITH_TPETRA
41
43 const Teuchos::RCP<const Tpetra::Map<int, types::signed_global_dof_index>>
44 &map)
45 : is_compressed(true)
46 , index_space_size(1 + map->getMaxAllGlobalIndex())
47 , largest_range(numbers::invalid_unsigned_int)
48{
49 Assert(map->getMinAllGlobalIndex() == 0,
51 "The Tpetra::Map does not contain the global index 0, "
52 "which means some entries are not present on any processor."));
53
54 // For a contiguous map, we do not need to go through the whole data...
55 if (map->isContiguous())
56 add_range(size_type(map->getMinGlobalIndex()),
57 size_type(map->getMaxGlobalIndex() + 1));
58 else
59 {
60# if DEAL_II_TRILINOS_VERSION_GTE(13, 4, 0)
61 const size_type n_indices = map->getLocalNumElements();
62# else
63 const size_type n_indices = map->getNodeNumElements();
64# endif
65 const types::signed_global_dof_index *indices =
66 map->getMyGlobalIndices().data();
67 add_indices(indices, indices + n_indices);
68 }
69 compress();
70}
71
72# endif // DEAL_II_TRILINOS_WITH_TPETRA
73
74
75// the 64-bit path uses a few different names, so put that into a separate
76// implementation
77
78# ifdef DEAL_II_WITH_64BIT_INDICES
79
80IndexSet::IndexSet(const Epetra_BlockMap &map)
81 : is_compressed(true)
82 , index_space_size(1 + map.MaxAllGID64())
83 , largest_range(numbers::invalid_unsigned_int)
84{
85 Assert(map.MinAllGID64() == 0,
87 "The Epetra_BlockMap does not contain the global index 0, "
88 "which means some entries are not present on any processor."));
89
90 // For a contiguous map, we do not need to go through the whole data...
91 if (map.LinearMap())
92 add_range(size_type(map.MinMyGID64()), size_type(map.MaxMyGID64() + 1));
93 else
94 {
95 const size_type n_indices = map.NumMyElements();
96 size_type *indices =
97 reinterpret_cast<size_type *>(map.MyGlobalElements64());
98 add_indices(indices, indices + n_indices);
99 }
100 compress();
101}
102
103# else
104
105// this is the standard 32-bit implementation
106
107IndexSet::IndexSet(const Epetra_BlockMap &map)
108 : is_compressed(true)
109 , index_space_size(1 + map.MaxAllGID())
110 , largest_range(numbers::invalid_unsigned_int)
111{
112 Assert(map.MinAllGID() == 0,
114 "The Epetra_BlockMap does not contain the global index 0, "
115 "which means some entries are not present on any processor."));
116
117 // For a contiguous map, we do not need to go through the whole data...
118 if (map.LinearMap())
119 add_range(size_type(map.MinMyGID()), size_type(map.MaxMyGID() + 1));
120 else
121 {
122 const size_type n_indices = map.NumMyElements();
123 unsigned int *indices =
124 reinterpret_cast<unsigned int *>(map.MyGlobalElements());
125 add_indices(indices, indices + n_indices);
126 }
127 compress();
128}
129
130# endif
131
132#endif // ifdef DEAL_II_WITH_TRILINOS
133
134
135
136void
138{
139 {
140 // we will, in the following, modify mutable variables. this can only
141 // work in multithreaded applications if we lock the data structures
142 // via a mutex, so that users can call 'const' functions from threads
143 // in parallel (and these 'const' functions can then call compress()
144 // which itself calls the current function)
145 std::lock_guard<std::mutex> lock(compress_mutex);
146
147 // see if any of the contiguous ranges can be merged. do not use
148 // std::vector::erase in-place as it is quadratic in the number of
149 // ranges. since the ranges are sorted by their first index, determining
150 // overlap isn't all that hard
151 std::vector<Range>::iterator store = ranges.begin();
152 for (std::vector<Range>::iterator i = ranges.begin(); i != ranges.end();)
153 {
154 std::vector<Range>::iterator next = i;
155 ++next;
156
157 size_type first_index = i->begin;
158 size_type last_index = i->end;
159
160 // see if we can merge any of the following ranges
161 while (next != ranges.end() && (next->begin <= last_index))
162 {
163 last_index = std::max(last_index, next->end);
164 ++next;
165 }
166 i = next;
167
168 // store the new range in the slot we last occupied
169 *store = Range(first_index, last_index);
170 ++store;
171 }
172 // use a compact array with exactly the right amount of storage
173 if (store != ranges.end())
174 {
175 std::vector<Range> new_ranges(ranges.begin(), store);
176 ranges.swap(new_ranges);
177 }
178
179 // now compute indices within set and the range with most elements
180 size_type next_index = 0, largest_range_size = 0;
181 for (std::vector<Range>::iterator i = ranges.begin(); i != ranges.end();
182 ++i)
183 {
184 Assert(i->begin < i->end, ExcInternalError());
185
186 i->nth_index_in_set = next_index;
187 next_index += (i->end - i->begin);
188 if (i->end - i->begin > largest_range_size)
189 {
190 largest_range_size = i->end - i->begin;
191 largest_range = i - ranges.begin();
192 }
193 }
194 is_compressed = true;
195
196 // check that next_index is correct. needs to be after the previous
197 // statement because we otherwise will get into an endless loop
198 Assert(next_index == n_elements(), ExcInternalError());
199 }
200
201#ifdef DEBUG
202 // A consistency check: We should only ever have added indices
203 // that are within the range of the index set. Instead of doing
204 // this in every one of the many functions that add indices,
205 // do this in the current, central location
206 for (const auto &range : ranges)
207 Assert((range.begin < index_space_size) && (range.end <= index_space_size),
208 ExcMessage("In the process of creating the current IndexSet "
209 "object, you added indices beyond the size of the index "
210 "space. Specifically, you added elements that form the "
211 "range [" +
212 std::to_string(range.begin) + "," +
213 std::to_string(range.end) +
214 "), but the size of the index space is only " +
215 std::to_string(index_space_size) + "."));
216#endif
217}
218
219
220
221#ifndef DOXYGEN
223IndexSet::operator&(const IndexSet &is) const
224{
225 Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
226
227 compress();
228 is.compress();
229
230 std::vector<Range>::const_iterator r1 = ranges.begin(),
231 r2 = is.ranges.begin();
232 IndexSet result(size());
233
234 while ((r1 != ranges.end()) && (r2 != is.ranges.end()))
235 {
236 // if r1 and r2 do not overlap at all, then move the pointer that sits
237 // to the left of the other up by one
238 if (r1->end <= r2->begin)
239 ++r1;
240 else if (r2->end <= r1->begin)
241 ++r2;
242 else
243 {
244 // the ranges must overlap somehow
245 Assert(((r1->begin <= r2->begin) && (r1->end > r2->begin)) ||
246 ((r2->begin <= r1->begin) && (r2->end > r1->begin)),
248
249 // add the overlapping range to the result
250 result.add_range(std::max(r1->begin, r2->begin),
251 std::min(r1->end, r2->end));
252
253 // now move that iterator that ends earlier one up. note that it has
254 // to be this one because a subsequent range may still have a chance
255 // of overlapping with the range that ends later
256 if (r1->end <= r2->end)
257 ++r1;
258 else
259 ++r2;
260 }
261 }
262
263 result.compress();
264 return result;
265}
266#endif
267
268
269
271IndexSet::get_view(const size_type begin, const size_type end) const
272{
273 Assert(begin <= end,
274 ExcMessage("End index needs to be larger or equal to begin index!"));
275 Assert(end <= size(),
276 ExcMessage("You are asking for a view into an IndexSet object "
277 "that would cover the sub-range [" +
278 std::to_string(begin) + ',' + std::to_string(end) +
279 "). But this is not a subset of the range "
280 "of the current object, which is [0," +
281 std::to_string(size()) + ")."));
282
283 IndexSet result(end - begin);
284 std::vector<Range>::const_iterator r1 = ranges.begin();
285
286 while (r1 != ranges.end())
287 {
288 if ((r1->end > begin) && (r1->begin < end))
289 {
290 result.add_range(std::max(r1->begin, begin) - begin,
291 std::min(r1->end, end) - begin);
292 }
293 else if (r1->begin >= end)
294 break;
295
296 ++r1;
297 }
298
299 result.compress();
300 return result;
301}
302
303
304
306IndexSet::get_view(const IndexSet &mask) const
307{
308 Assert(size() == mask.size(),
309 ExcMessage("The mask must have the same size index space "
310 "as the index set it is applied to."));
311
312 // If 'other' is an empty set, then the view is also empty:
313 if (mask == IndexSet())
314 return {};
315
316 // For everything, it is more efficient to work on compressed sets:
317 compress();
318 mask.compress();
319
320 // If 'other' has a single range, then we can just defer to the
321 // previous function
322 if (mask.ranges.size() == 1)
323 return get_view(mask.ranges[0].begin, mask.ranges[0].end);
324
325 // For the general case where the mask is an arbitrary set,
326 // the situation is slightly more complicated. We need to walk
327 // the ranges of the two index sets in parallel and search for
328 // overlaps, and then appropriately shift
329
330 // we save all new ranges to our IndexSet in an temporary vector and
331 // add all of them in one go at the end.
332 std::vector<Range> new_ranges;
333
334 std::vector<Range>::iterator own_it = ranges.begin();
335 std::vector<Range>::iterator mask_it = mask.ranges.begin();
336
337 while ((own_it != ranges.end()) && (mask_it != mask.ranges.end()))
338 {
339 // If our own range lies completely ahead of the current
340 // range in the mask, move forward and start the loop body
341 // anew. If this was the last range, the 'while' loop above
342 // will terminate, so we don't have to check for end iterators
343 if (own_it->end <= mask_it->begin)
344 {
345 ++own_it;
346 continue;
347 }
348
349 // Do the same if the current mask range lies completely ahead of
350 // the current range of the this object:
351 if (mask_it->end <= own_it->begin)
352 {
353 ++mask_it;
354 continue;
355 }
356
357 // Now own_it and other_it overlap. Check that that is true by
358 // enumerating the cases that can happen. This is
359 // surprisingly tricky because the two intervals can intersect in
360 // a number of different ways, but there really are only the four
361 // following possibilities:
362
363 // Case 1: our interval overlaps the left end of the other interval
364 //
365 // So we need to add the elements from the first element of the mask's
366 // interval to the end of our own interval. But we need to shift the
367 // indices so that they correspond to the how many'th element within the
368 // mask this is; fortunately (because we compressed the mask), this
369 // is recorded in the mask's ranges.
370 if ((own_it->begin <= mask_it->begin) && (own_it->end <= mask_it->end))
371 {
372 new_ranges.emplace_back(mask_it->begin - mask_it->nth_index_in_set,
373 own_it->end - mask_it->nth_index_in_set);
374 }
375 else
376 // Case 2:our interval overlaps the tail end of the other interval
377 if ((mask_it->begin <= own_it->begin) && (mask_it->end <= own_it->end))
378 {
379 const size_type offset_within_mask_interval =
380 own_it->begin - mask_it->begin;
381 new_ranges.emplace_back(mask_it->nth_index_in_set +
382 offset_within_mask_interval,
383 mask_it->nth_index_in_set +
384 (mask_it->end - mask_it->begin));
385 }
386 else
387 // Case 3: Our own interval completely encloses the other interval
388 if ((own_it->begin <= mask_it->begin) &&
389 (own_it->end >= mask_it->end))
390 {
391 new_ranges.emplace_back(mask_it->begin -
392 mask_it->nth_index_in_set,
393 mask_it->end - mask_it->nth_index_in_set);
394 }
395 else
396 // Case 3: The other interval completely encloses our own interval
397 if ((mask_it->begin <= own_it->begin) &&
398 (mask_it->end >= own_it->end))
399 {
400 const size_type offset_within_mask_interval =
401 own_it->begin - mask_it->begin;
402 new_ranges.emplace_back(mask_it->nth_index_in_set +
403 offset_within_mask_interval,
404 mask_it->nth_index_in_set +
405 offset_within_mask_interval +
406 (own_it->end - own_it->begin));
407 }
408 else
410
411 // We considered the overlap of these two intervals. It may of course
412 // be that one of them overlaps with another one, but that can only
413 // be the case for the interval that extends further to the right. So
414 // we can safely move on from the interval that terminates earlier:
415 if (own_it->end < mask_it->end)
416 ++own_it;
417 else if (mask_it->end < own_it->end)
418 ++mask_it;
419 else
420 {
421 // The intervals ended at the same point. We can move on from both.
422 // (The algorithm would also work if we only moved on from one,
423 // but we can micro-optimize here without too much effort.)
424 ++own_it;
425 ++mask_it;
426 }
427 }
428
429 // Now turn the ranges of overlap we have accumulated into an IndexSet in
430 // its own right:
431 IndexSet result(mask.n_elements());
432 for (const auto &range : new_ranges)
433 result.add_range(range.begin, range.end);
434 result.compress();
435
436 return result;
437}
438
439
440
441std::vector<IndexSet>
443 const std::vector<types::global_dof_index> &n_indices_per_block) const
444{
445 std::vector<IndexSet> partitioned;
446 const unsigned int n_blocks = n_indices_per_block.size();
447
448 partitioned.reserve(n_blocks);
449 types::global_dof_index start = 0;
450 for (const auto n_block_indices : n_indices_per_block)
451 {
452 partitioned.push_back(this->get_view(start, start + n_block_indices));
453 start += n_block_indices;
454 }
455
456#ifdef DEBUG
458 for (const auto &partition : partitioned)
459 {
460 sum += partition.size();
461 }
462 AssertDimension(sum, this->size());
463#endif
464
465 return partitioned;
466}
467
468
469
470void
472{
473 compress();
474 other.compress();
475 is_compressed = false;
476
477
478 // we save all new ranges to our IndexSet in an temporary vector and
479 // add all of them in one go at the end.
480 std::vector<Range> new_ranges;
481
482 std::vector<Range>::iterator own_it = ranges.begin();
483 std::vector<Range>::iterator other_it = other.ranges.begin();
484
485 while (own_it != ranges.end() && other_it != other.ranges.end())
486 {
487 // advance own iterator until we get an overlap
488 if (own_it->end <= other_it->begin)
489 {
490 new_ranges.push_back(*own_it);
491 ++own_it;
492 continue;
493 }
494 // we are done with other_it, so advance
495 if (own_it->begin >= other_it->end)
496 {
497 ++other_it;
498 continue;
499 }
500
501 // Now own_it and other_it overlap. First save the part of own_it that
502 // is before other_it (if not empty).
503 if (own_it->begin < other_it->begin)
504 {
505 Range r(own_it->begin, other_it->begin);
506 r.nth_index_in_set = 0; // fix warning of unused variable
507 new_ranges.push_back(r);
508 }
509 // change own_it to the sub range behind other_it. Do not delete own_it
510 // in any case. As removal would invalidate iterators, we just shrink
511 // the range to an empty one.
512 own_it->begin = other_it->end;
513 if (own_it->begin > own_it->end)
514 {
515 own_it->begin = own_it->end;
516 ++own_it;
517 }
518
519 // continue without advancing iterators, the right one will be advanced
520 // next.
521 }
522
523 // make sure to take over the remaining ranges
524 for (; own_it != ranges.end(); ++own_it)
525 new_ranges.push_back(*own_it);
526
527 ranges.clear();
528
529 // done, now add the temporary ranges
530 const std::vector<Range>::iterator end = new_ranges.end();
531 for (std::vector<Range>::iterator it = new_ranges.begin(); it != end; ++it)
532 add_range(it->begin, it->end);
533
534 compress();
535}
536
537
538
541{
542 IndexSet set(this->size() * other.size());
543 for (const auto el : *this)
544 set.add_indices(other, el * other.size());
545 set.compress();
546 return set;
547}
548
549
550
553{
554 Assert(is_empty() == false,
556 "pop_back() failed, because this IndexSet contains no entries."));
557
558 const size_type index = ranges.back().end - 1;
559 --ranges.back().end;
560
561 if (ranges.back().begin == ranges.back().end)
562 ranges.pop_back();
563
564 return index;
565}
566
567
568
571{
572 Assert(is_empty() == false,
574 "pop_front() failed, because this IndexSet contains no entries."));
575
576 const size_type index = ranges.front().begin;
577 ++ranges.front().begin;
578
579 if (ranges.front().begin == ranges.front().end)
580 ranges.erase(ranges.begin());
581
582 // We have to set this in any case, because nth_index_in_set is no longer
583 // up to date for all but the first range
584 is_compressed = false;
585
586 return index;
587}
588
589
590
591void
593{
594 // if the inserted range is already within the range we find by lower_bound,
595 // there is no need to do anything; we do not try to be clever here and
596 // leave all other work to compress().
597 const auto insert_position =
598 Utilities::lower_bound(ranges.begin(), ranges.end(), new_range);
599 if (insert_position == ranges.end() ||
600 insert_position->begin > new_range.begin ||
601 insert_position->end < new_range.end)
602 ranges.insert(insert_position, new_range);
603}
604
605
606
607void
609 boost::container::small_vector<std::pair<size_type, size_type>, 200>
610 &tmp_ranges,
611 const bool ranges_are_sorted)
612{
613 if (!ranges_are_sorted)
614 std::sort(tmp_ranges.begin(), tmp_ranges.end());
615
616 // if we have many ranges, we first construct a temporary index set (where
617 // we add ranges in a consecutive way, so fast), otherwise, we work with
618 // add_range(). the number 9 is chosen heuristically given the fact that
619 // there are typically up to 8 independent ranges when adding the degrees of
620 // freedom on a 3d cell or 9 when adding degrees of freedom of faces. if
621 // doing cell-by-cell additions, we want to avoid repeated calls to
622 // IndexSet::compress() which gets called upon merging two index sets, so we
623 // want to be in the other branch then.
624 if (tmp_ranges.size() > 9)
625 {
626 IndexSet tmp_set(size());
627 tmp_set.ranges.reserve(tmp_ranges.size());
628 for (const auto &i : tmp_ranges)
629 tmp_set.add_range(i.first, i.second);
630
631 // Case if we have zero or just one range: Add into the other set with
632 // its indices, as that is cheaper
633 if (this->ranges.size() <= 1)
634 {
635 if (this->ranges.size() == 1)
636 tmp_set.add_range(ranges[0].begin, ranges[0].end);
637 std::swap(*this, tmp_set);
638 }
639 else
640 this->add_indices(tmp_set);
641 }
642 else
643 for (const auto &i : tmp_ranges)
644 add_range(i.first, i.second);
645}
646
647
648
649void
650IndexSet::add_indices(const IndexSet &other, const size_type offset)
651{
652 if ((this == &other) && (offset == 0))
653 return;
654
655 if (other.ranges.size() != 0)
656 {
657 AssertIndexRange(other.ranges.back().end - 1, index_space_size);
658 }
659
660 compress();
661 other.compress();
662
663 std::vector<Range>::const_iterator r1 = ranges.begin(),
664 r2 = other.ranges.begin();
665
666 std::vector<Range> new_ranges;
667 // just get the start and end of the ranges right in this method, everything
668 // else will be done in compress()
669 while (r1 != ranges.end() || r2 != other.ranges.end())
670 {
671 // the two ranges do not overlap or we are at the end of one of the
672 // ranges
673 if (r2 == other.ranges.end() ||
674 (r1 != ranges.end() && r1->end < (r2->begin + offset)))
675 {
676 new_ranges.push_back(*r1);
677 ++r1;
678 }
679 else if (r1 == ranges.end() || (r2->end + offset) < r1->begin)
680 {
681 new_ranges.emplace_back(r2->begin + offset, r2->end + offset);
682 ++r2;
683 }
684 else
685 {
686 // ok, we do overlap, so just take the combination of the current
687 // range (do not bother to merge with subsequent ranges)
688 Range next(std::min(r1->begin, r2->begin + offset),
689 std::max(r1->end, r2->end + offset));
690 new_ranges.push_back(next);
691 ++r1;
692 ++r2;
693 }
694 }
695 ranges.swap(new_ranges);
696
697 is_compressed = false;
698 compress();
699}
700
701
702
703bool
705{
706 Assert(size() == other.size(),
707 ExcMessage("One index set can only be a subset of another if they "
708 "describe index spaces of the same size. The ones in "
709 "question here have sizes " +
710 std::to_string(size()) + " and " +
711 std::to_string(other.size()) + "."));
712
713 // See whether there are indices in the current set that are not in 'other'.
714 // If so, then this is clearly not a subset of 'other'.
715 IndexSet A_minus_B = *this;
716 A_minus_B.subtract_set(other);
717 if (A_minus_B.n_elements() > 0)
718 return false;
719 else
720 // Else, every index in 'this' is also in 'other', since we ended up
721 // with an empty set upon subtraction. This means that we have a subset:
722 return true;
723}
724
725
726
727void
728IndexSet::write(std::ostream &out) const
729{
730 compress();
731 out << size() << " ";
732 out << ranges.size() << std::endl;
733 std::vector<Range>::const_iterator r = ranges.begin();
734 for (; r != ranges.end(); ++r)
735 {
736 out << r->begin << " " << r->end << std::endl;
737 }
738}
739
740
741
742void
743IndexSet::read(std::istream &in)
744{
745 AssertThrow(in.fail() == false, ExcIO());
746
747 size_type s;
748 unsigned int n_ranges;
749
750 in >> s >> n_ranges;
751 ranges.clear();
752 set_size(s);
753 for (unsigned int i = 0; i < n_ranges; ++i)
754 {
755 AssertThrow(in.fail() == false, ExcIO());
756
757 size_type b, e;
758 in >> b >> e;
759 add_range(b, e);
760 }
761}
762
763
764void
765IndexSet::block_write(std::ostream &out) const
766{
767 AssertThrow(out.fail() == false, ExcIO());
768 out.write(reinterpret_cast<const char *>(&index_space_size),
769 sizeof(index_space_size));
770 std::size_t n_ranges = ranges.size();
771 out.write(reinterpret_cast<const char *>(&n_ranges), sizeof(n_ranges));
772 if (ranges.empty() == false)
773 out.write(reinterpret_cast<const char *>(&*ranges.begin()),
774 ranges.size() * sizeof(Range));
775 AssertThrow(out.fail() == false, ExcIO());
776}
777
778void
779IndexSet::block_read(std::istream &in)
780{
782 std::size_t n_ranges;
783 in.read(reinterpret_cast<char *>(&size), sizeof(size));
784 in.read(reinterpret_cast<char *>(&n_ranges), sizeof(n_ranges));
785 // we have to clear ranges first
786 ranges.clear();
787 set_size(size);
788 ranges.resize(n_ranges, Range(0, 0));
789 if (n_ranges != 0u)
790 in.read(reinterpret_cast<char *>(&*ranges.begin()),
791 ranges.size() * sizeof(Range));
792
793 do_compress(); // needed so that largest_range can be recomputed
794}
795
796
797
798bool
800{
801 // get the element after which we would have to insert a range that
802 // consists of all elements from this element to the end of the index
803 // range plus one. after this call we know that if p!=end() then
804 // p->begin<=index unless there is no such range at all
805 //
806 // if the searched for element is an element of this range, then we're
807 // done. otherwise, the element can't be in one of the following ranges
808 // because otherwise p would be a different iterator
809 //
810 // since we already know the position relative to the largest range (we
811 // called compress!), we can perform the binary search on ranges with
812 // lower/higher number compared to the largest range
813 std::vector<Range>::const_iterator p = std::upper_bound(
814 ranges.begin() +
815 (index < ranges[largest_range].begin ? 0 : largest_range + 1),
816 index < ranges[largest_range].begin ? ranges.begin() + largest_range :
817 ranges.end(),
818 Range(index, size() + 1));
819
820 if (p == ranges.begin())
821 return ((index >= p->begin) && (index < p->end));
822
823 Assert((p == ranges.end()) || (p->begin > index), ExcInternalError());
824
825 // now move to that previous range
826 --p;
827 Assert(p->begin <= index, ExcInternalError());
828
829 return (p->end > index);
830}
831
832
833
836{
837 // find out which chunk the local index n belongs to by using a binary
838 // search. the comparator is based on the end of the ranges.
839 Range r(n, n + 1);
840 r.nth_index_in_set = n;
841
842 const std::vector<Range>::const_iterator p = Utilities::lower_bound(
843 ranges.begin(), ranges.end(), r, Range::nth_index_compare);
844
845 Assert(p != ranges.end(), ExcInternalError());
846 return p->begin + (n - p->nth_index_in_set);
847}
848
849
850
853{
854 // we could try to use the main range for splitting up the search range, but
855 // since we only come here when the largest range did not contain the index,
856 // there is little gain from doing a first step manually.
857 Range r(n, n);
858 std::vector<Range>::const_iterator p =
860
861 // if n is not in this set
862 if (p == ranges.end() || p->end == n || p->begin > n)
864
865 Assert(p != ranges.end(), ExcInternalError());
866 Assert(p->begin <= n, ExcInternalError());
867 Assert(n < p->end, ExcInternalError());
868 return (n - p->begin) + p->nth_index_in_set;
869}
870
871
872
874IndexSet::at(const size_type global_index) const
875{
876 compress();
877 AssertIndexRange(global_index, size());
878
879 if (ranges.empty())
880 return end();
881
882 std::vector<Range>::const_iterator main_range =
883 ranges.begin() + largest_range;
884
885 Range r(global_index, global_index + 1);
886 // This optimization makes the bounds for lower_bound smaller by checking
887 // the largest range first.
888 std::vector<Range>::const_iterator range_begin, range_end;
889 if (global_index < main_range->begin)
890 {
891 range_begin = ranges.begin();
892 range_end = main_range;
893 }
894 else
895 {
896 range_begin = main_range;
897 range_end = ranges.end();
898 }
899
900 // This will give us the first range p=[a,b[ with b>=global_index using
901 // a binary search
902 const std::vector<Range>::const_iterator p =
903 Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
904
905 // We couldn't find a range, which means we have no range that contains
906 // global_index and also no range behind it, meaning we need to return end().
907 if (p == ranges.end())
908 return end();
909
910 // Finally, we can have two cases: Either global_index is not in [a,b[,
911 // which means we need to return an iterator to a because global_index, ...,
912 // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
913 // [a,b[ and we will return an iterator pointing directly at global_index
914 // (else branch).
915 if (global_index < p->begin)
916 return {this, static_cast<size_type>(p - ranges.begin()), p->begin};
917 else
918 return {this, static_cast<size_type>(p - ranges.begin()), global_index};
919}
920
921
922
923std::vector<IndexSet::size_type>
925{
926 compress();
927
928 std::vector<size_type> indices;
929 indices.reserve(n_elements());
930
931 for (const auto &range : ranges)
932 for (size_type entry = range.begin; entry < range.end; ++entry)
933 indices.push_back(entry);
934
935 Assert(indices.size() == n_elements(), ExcInternalError());
936
937 return indices;
938}
939
940
941
942void
943IndexSet::fill_index_vector(std::vector<size_type> &indices) const
944{
945 indices = get_index_vector();
946}
947
948
949
950#ifdef DEAL_II_WITH_TRILINOS
951# ifdef DEAL_II_TRILINOS_WITH_TPETRA
952
953Tpetra::Map<int, types::signed_global_dof_index>
955 const bool overlapping) const
956{
957 return *make_tpetra_map_rcp(communicator, overlapping);
958}
959
960
961
962Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index>>
964 const bool overlapping) const
965{
966 compress();
967 (void)communicator;
968
969# ifdef DEBUG
970 if (!overlapping)
971 {
972 const size_type n_global_elements =
973 Utilities::MPI::sum(n_elements(), communicator);
974 Assert(n_global_elements == size(),
975 ExcMessage("You are trying to create an Tpetra::Map object "
976 "that partitions elements of an index set "
977 "between processors. However, the union of the "
978 "index sets on different processors does not "
979 "contain all indices exactly once: the sum of "
980 "the number of entries the various processors "
981 "want to store locally is " +
982 std::to_string(n_global_elements) +
983 " whereas the total size of the object to be "
984 "allocated is " +
985 std::to_string(size()) +
986 ". In other words, there are "
987 "either indices that are not spoken for "
988 "by any processor, or there are indices that are "
989 "claimed by multiple processors."));
990 }
991# endif
992
993 // Find out if the IndexSet is ascending and 1:1. This corresponds to a
994 // linear Tpetra::Map. Overlapping IndexSets are never 1:1.
995 const bool linear =
996 overlapping ? false : is_ascending_and_one_to_one(communicator);
997 if (linear)
999 Tpetra::Map<int, types::signed_global_dof_index>>(
1000 size(),
1001 n_elements(),
1002 0,
1003# ifdef DEAL_II_WITH_MPI
1004 Utilities::Trilinos::internal::make_rcp<Teuchos::MpiComm<int>>(
1005 communicator)
1006# else
1007 Utilities::Trilinos::internal::make_rcp<Teuchos::Comm<int>>()
1008# endif // DEAL_II_WITH_MPI
1009 );
1010 else
1011 {
1012 const std::vector<size_type> indices = get_index_vector();
1013 std::vector<types::signed_global_dof_index> int_indices(indices.size());
1014 std::copy(indices.begin(), indices.end(), int_indices.begin());
1015 const Teuchos::ArrayView<types::signed_global_dof_index> arr_view(
1016 int_indices);
1017
1019 Tpetra::Map<int, types::signed_global_dof_index>>(
1020 size(),
1021 arr_view,
1022 0,
1023# ifdef DEAL_II_WITH_MPI
1024 Utilities::Trilinos::internal::make_rcp<Teuchos::MpiComm<int>>(
1025 communicator)
1026# else
1027 Utilities::Trilinos::internal::make_rcp<Teuchos::Comm<int>>()
1028# endif // DEAL_II_WITH_MPI
1029 );
1030 }
1031}
1032# endif
1033
1034
1035
1036Epetra_Map
1038 const bool overlapping) const
1039{
1040 compress();
1041 (void)communicator;
1042
1043# ifdef DEBUG
1044 if (!overlapping)
1045 {
1046 const size_type n_global_elements =
1047 Utilities::MPI::sum(n_elements(), communicator);
1048 Assert(n_global_elements == size(),
1049 ExcMessage("You are trying to create an Epetra_Map object "
1050 "that partitions elements of an index set "
1051 "between processors. However, the union of the "
1052 "index sets on different processors does not "
1053 "contain all indices exactly once: the sum of "
1054 "the number of entries the various processors "
1055 "want to store locally is " +
1056 std::to_string(n_global_elements) +
1057 " whereas the total size of the object to be "
1058 "allocated is " +
1059 std::to_string(size()) +
1060 ". In other words, there are "
1061 "either indices that are not spoken for "
1062 "by any processor, or there are indices that are "
1063 "claimed by multiple processors."));
1064 }
1065# endif
1066
1067 // Find out if the IndexSet is ascending and 1:1. This corresponds to a
1068 // linear EpetraMap. Overlapping IndexSets are never 1:1.
1069 const bool linear =
1070 overlapping ? false : is_ascending_and_one_to_one(communicator);
1071
1072 if (linear)
1073 return Epetra_Map(TrilinosWrappers::types::int_type(size()),
1075 0,
1076# ifdef DEAL_II_WITH_MPI
1077 Epetra_MpiComm(communicator)
1078# else
1079 Epetra_SerialComm()
1080# endif
1081 );
1082 else
1083 {
1084 const std::vector<size_type> indices = get_index_vector();
1085 return Epetra_Map(
1088 (n_elements() > 0 ?
1089 reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1090 indices.data()) :
1091 nullptr),
1092 0,
1093# ifdef DEAL_II_WITH_MPI
1094 Epetra_MpiComm(communicator)
1095# else
1096 Epetra_SerialComm()
1097# endif
1098 );
1099 }
1100}
1101#endif
1102
1103
1104#ifdef DEAL_II_WITH_PETSC
1105IS
1106IndexSet::make_petsc_is(const MPI_Comm communicator) const
1107{
1108 std::vector<size_type> indices;
1109 fill_index_vector(indices);
1110
1111 PetscInt n = indices.size();
1112 std::vector<PetscInt> pindices(indices.begin(), indices.end());
1113
1114 IS is;
1115 PetscErrorCode ierr =
1116 ISCreateGeneral(communicator, n, pindices.data(), PETSC_COPY_VALUES, &is);
1117 AssertThrow(ierr == 0, ExcPETScError(ierr));
1118
1119 return is;
1120}
1121#endif
1122
1123
1124
1125bool
1127{
1128 // If the sum of local elements does not add up to the total size,
1129 // the IndexSet can't be complete.
1130 const size_type n_global_elements =
1131 Utilities::MPI::sum(n_elements(), communicator);
1132 if (n_global_elements != size())
1133 return false;
1134
1135 if (n_global_elements == 0)
1136 return true;
1137
1138#ifdef DEAL_II_WITH_MPI
1139 // Non-contiguous IndexSets can't be linear.
1140 const bool all_contiguous =
1141 (Utilities::MPI::min(is_contiguous() ? 1 : 0, communicator) == 1);
1142 if (!all_contiguous)
1143 return false;
1144
1145 bool is_globally_ascending = true;
1146 // we know that there is only one interval
1147 types::global_dof_index first_local_dof = (n_elements() > 0) ?
1148 *(begin_intervals()->begin()) :
1150
1151 const unsigned int my_rank = Utilities::MPI::this_mpi_process(communicator);
1152 const std::vector<types::global_dof_index> global_dofs =
1153 Utilities::MPI::gather(communicator, first_local_dof, 0);
1154
1155 if (my_rank == 0)
1156 {
1157 // find out if the received std::vector is ascending
1158 types::global_dof_index index = 0;
1159 while (global_dofs[index] == numbers::invalid_dof_index)
1160 ++index;
1161 types::global_dof_index old_dof = global_dofs[index++];
1162 for (; index < global_dofs.size(); ++index)
1163 {
1164 const types::global_dof_index new_dof = global_dofs[index];
1165 if (new_dof != numbers::invalid_dof_index)
1166 {
1167 if (new_dof <= old_dof)
1168 {
1169 is_globally_ascending = false;
1170 break;
1171 }
1172 else
1173 old_dof = new_dof;
1174 }
1175 }
1176 }
1177
1178 // now broadcast the result
1179 int is_ascending = is_globally_ascending ? 1 : 0;
1180 int ierr = MPI_Bcast(&is_ascending, 1, MPI_INT, 0, communicator);
1181 AssertThrowMPI(ierr);
1182
1183 return (is_ascending == 1);
1184#else
1185 return true;
1186#endif // DEAL_II_WITH_MPI
1187}
1188
1189
1190
1191std::size_t
1199
1200
1201
ElementIterator begin() const
Definition index_set.h:1255
bool is_subset_of(const IndexSet &other) const
Definition index_set.cc:704
bool is_element_binary_search(const size_type local_index) const
Definition index_set.cc:799
size_type index_within_set_binary_search(const size_type global_index) const
Definition index_set.cc:852
size_type largest_range
Definition index_set.h:1115
IS make_petsc_is(const MPI_Comm communicator=MPI_COMM_WORLD) const
bool is_ascending_and_one_to_one(const MPI_Comm communicator) const
bool is_contiguous() const
Definition index_set.h:1907
void do_compress() const
Definition index_set.cc:137
ElementIterator at(const size_type global_index) const
Definition index_set.cc:874
size_type size() const
Definition index_set.h:1766
void fill_index_vector(std::vector< size_type > &indices) const
Definition index_set.cc:943
bool is_empty() const
Definition index_set.h:1916
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
Definition index_set.cc:442
size_type n_elements() const
Definition index_set.h:1924
void add_range_lower_bound(const Range &range)
Definition index_set.cc:592
ElementIterator begin() const
Definition index_set.h:1700
size_type pop_front()
Definition index_set.cc:570
void set_size(const size_type size)
Definition index_set.h:1754
void read(std::istream &in)
Definition index_set.cc:743
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
IndexSet tensor_product(const IndexSet &other) const
Definition index_set.cc:540
void write(std::ostream &out) const
Definition index_set.cc:728
void block_read(std::istream &in)
Definition index_set.cc:779
bool is_compressed
Definition index_set.h:1098
void add_ranges_internal(boost::container::small_vector< std::pair< size_type, size_type >, 200 > &tmp_ranges, const bool ranges_are_sorted)
Definition index_set.cc:608
IntervalIterator begin_intervals() const
Definition index_set.h:1721
Tpetra::Map< int, types::signed_global_dof_index > make_tpetra_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition index_set.cc:954
std::vector< Range > ranges
Definition index_set.h:1088
void subtract_set(const IndexSet &other)
Definition index_set.cc:471
ElementIterator end() const
Definition index_set.h:1712
Threads::Mutex compress_mutex
Definition index_set.h:1121
size_type index_space_size
Definition index_set.h:1104
void block_write(std::ostream &out) const
Definition index_set.cc:765
IndexSet get_view(const size_type begin, const size_type end) const
Definition index_set.cc:271
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1793
std::size_t memory_consumption() const
size_type nth_index_in_set_binary_search(const size_type local_index) const
Definition index_set.cc:835
void compress() const
Definition index_set.h:1774
Teuchos::RCP< Tpetra::Map< int, types::signed_global_dof_index > > make_tpetra_map_rcp(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition index_set.cc:963
size_type pop_back()
Definition index_set.cc:552
std::vector< size_type > get_index_vector() const
Definition index_set.cc:924
types::global_dof_index size_type
Definition index_set.h:77
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1821
IndexSet operator&(const IndexSet &is) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
static ::ExceptionBase & ExcIO()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_ASSERT_UNREACHABLE()
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
T sum(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:107
std::vector< T > gather(const MPI_Comm comm, const T &object_to_send, const unsigned int root_process=0)
Teuchos::RCP< T > make_rcp(Args &&...args)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition utilities.h:1032
const types::global_dof_index invalid_dof_index
Definition types.h:252
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition index_set.h:1046
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition index_set.h:1052
size_type end
Definition index_set.h:1014
size_type nth_index_in_set
Definition index_set.h:1016
size_type begin
Definition index_set.h:1013