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