Reference documentation for deal.II version GIT 29f9da0a34 2023-12-07 10:00:01+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\}}\)
index_set.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2023 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 #include <deal.II/base/index_set.h>
18 #include <deal.II/base/mpi.h>
19 
20 #include <deal.II/lac/exceptions.h>
21 
22 #include <vector>
23 
24 #ifdef DEAL_II_WITH_TRILINOS
25 # ifdef DEAL_II_WITH_MPI
26 # include <Epetra_MpiComm.h>
27 # endif
28 # include <Epetra_Map.h>
29 # include <Epetra_SerialComm.h>
30 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
31 # include <Tpetra_Map.hpp>
32 # endif
33 #endif
34 
36 
37 
38 
39 #ifdef DEAL_II_WITH_TRILINOS
40 
41 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
42 
44  Teuchos::RCP<const Tpetra::Map<int, types::signed_global_dof_index>> 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,
50  ExcMessage(
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 
80 IndexSet::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,
86  ExcMessage(
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 
107 IndexSet::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,
113  ExcMessage(
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 
136 void
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 " +
216 #endif
217 }
218 
219 
220 
221 #ifndef DOXYGEN
222 IndexSet
223 IndexSet::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)),
247  ExcInternalError());
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 
270 IndexSet
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 [" +
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 
305 IndexSet
306 IndexSet::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
409  Assert(false, ExcInternalError());
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 
441 std::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 
470 void
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 
539 IndexSet
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,
555  ExcMessage(
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,
573  ExcMessage(
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 
591 void
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 
607 void
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 
649 void
650 IndexSet::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 
703 bool
704 IndexSet::is_subset_of(const IndexSet &other) const
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 
727 void
728 IndexSet::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 
742 void
743 IndexSet::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 
764 void
765 IndexSet::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 
778 void
779 IndexSet::block_read(std::istream &in)
780 {
781  size_type size;
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 
798 bool
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 
875 {
876  compress();
878 
879  if (ranges.empty())
880  return end();
881 
882  std::vector<Range>::const_iterator main_range =
883  ranges.begin() + largest_range;
884 
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 
923 std::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 
942 void
943 IndexSet::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 
953 Tpetra::Map<int, types::signed_global_dof_index>
954 IndexSet::make_tpetra_map(const MPI_Comm communicator,
955  const bool overlapping) const
956 {
957  return *make_tpetra_map_rcp(communicator, overlapping);
958 }
959 
960 
961 
962 Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index>>
963 IndexSet::make_tpetra_map_rcp(const MPI_Comm communicator,
964  const bool overlapping) const
965 {
966  compress();
967  (void)communicator;
968 
969 # ifdef DEBUG
970  if (!overlapping)
971  {
973  Utilities::MPI::sum(n_elements(), communicator);
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 " +
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_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 
1036 Epetra_Map
1037 IndexSet::make_trilinos_map(const MPI_Comm communicator,
1038  const bool overlapping) const
1039 {
1040  compress();
1041  (void)communicator;
1042 
1043 # ifdef DEBUG
1044  if (!overlapping)
1045  {
1047  Utilities::MPI::sum(n_elements(), communicator);
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 " +
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
1105 IS
1106 IndexSet::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 
1125 bool
1126 IndexSet::is_ascending_and_one_to_one(const MPI_Comm communicator) const
1127 {
1128  // If the sum of local elements does not add up to the total size,
1129  // the IndexSet can't be complete.
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 
1191 std::size_t
1193 {
1197  sizeof(compress_mutex));
1198 }
1199 
1200 
1201 
ElementIterator begin() const
Definition: index_set.h:1250
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:1110
IS make_petsc_is(const MPI_Comm communicator=MPI_COMM_WORLD) const
Definition: index_set.cc:1106
bool is_ascending_and_one_to_one(const MPI_Comm communicator) const
Definition: index_set.cc:1126
bool is_contiguous() const
Definition: index_set.h:1902
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:1761
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:943
bool is_empty() const
Definition: index_set.h:1911
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:1919
void add_range_lower_bound(const Range &range)
Definition: index_set.cc:592
ElementIterator begin() const
Definition: index_set.h:1695
size_type pop_front()
Definition: index_set.cc:570
void set_size(const size_type size)
Definition: index_set.h:1749
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
Definition: index_set.cc:1037
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:1093
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:1716
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:1083
void subtract_set(const IndexSet &other)
Definition: index_set.cc:471
ElementIterator end() const
Definition: index_set.h:1707
Threads::Mutex compress_mutex
Definition: index_set.h:1116
size_type index_space_size
Definition: index_set.h:1099
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:1788
std::size_t memory_consumption() const
Definition: index_set.cc:1192
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:1769
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:78
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1816
IndexSet operator&(const IndexSet &is) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1947
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
void swap(MemorySpaceData< T, MemorySpace > &u, MemorySpaceData< T, MemorySpace > &v)
std::string to_string(const T &t)
Definition: patterns.h:2391
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
TrilinosWrappers::types::int64_type n_global_elements(const Epetra_BlockMap &map)
T sum(const T &t, const MPI_Comm mpi_communicator)
std::vector< T > gather(const MPI_Comm comm, const T &object_to_send, const unsigned int root_process=0)
T min(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition: mpi.cc:164
Teuchos::RCP< T > make_rcp(Args &&...args)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1008
void copy(const T *begin, const T *end, U *dest)
static const unsigned int invalid_unsigned_int
Definition: types.h:221
const types::global_dof_index invalid_dof_index
Definition: types.h:253
int signed_global_dof_index
Definition: types.h:93
unsigned int global_dof_index
Definition: types.h:82
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:1041
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:1047
size_type end
Definition: index_set.h:1009
size_type nth_index_in_set
Definition: index_set.h:1011
size_type begin
Definition: index_set.h:1008