Reference documentation for deal.II version GIT 5dd0dfed0b 2022-12-08 21:45: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\}}\)
index_set.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2022 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 <vector>
21 
22 #ifdef DEAL_II_WITH_TRILINOS
23 # ifdef DEAL_II_WITH_MPI
24 # include <Epetra_MpiComm.h>
25 # endif
26 # include <Epetra_Map.h>
27 # include <Epetra_SerialComm.h>
28 #endif
29 
31 
32 
33 
34 #ifdef DEAL_II_WITH_TRILINOS
35 
36 // the 64-bit path uses a few different names, so put that into a separate
37 // implementation
38 
39 # ifdef DEAL_II_WITH_64BIT_INDICES
40 
41 IndexSet::IndexSet(const Epetra_BlockMap &map)
42  : is_compressed(true)
43  , index_space_size(1 + map.MaxAllGID64())
44  , largest_range(numbers::invalid_unsigned_int)
45 {
46  Assert(map.MinAllGID64() == 0,
47  ExcMessage(
48  "The Epetra_BlockMap does not contain the global index 0, "
49  "which means some entries are not present on any processor."));
50 
51  // For a contiguous map, we do not need to go through the whole data...
52  if (map.LinearMap())
53  add_range(size_type(map.MinMyGID64()), size_type(map.MaxMyGID64() + 1));
54  else
55  {
56  const size_type n_indices = map.NumMyElements();
57  size_type * indices =
58  reinterpret_cast<size_type *>(map.MyGlobalElements64());
59  add_indices(indices, indices + n_indices);
60  }
61  compress();
62 }
63 
64 # else
65 
66 // this is the standard 32-bit implementation
67 
68 IndexSet::IndexSet(const Epetra_BlockMap &map)
69  : is_compressed(true)
70  , index_space_size(1 + map.MaxAllGID())
71  , largest_range(numbers::invalid_unsigned_int)
72 {
73  Assert(map.MinAllGID() == 0,
74  ExcMessage(
75  "The Epetra_BlockMap does not contain the global index 0, "
76  "which means some entries are not present on any processor."));
77 
78  // For a contiguous map, we do not need to go through the whole data...
79  if (map.LinearMap())
80  add_range(size_type(map.MinMyGID()), size_type(map.MaxMyGID() + 1));
81  else
82  {
83  const size_type n_indices = map.NumMyElements();
84  unsigned int * indices =
85  reinterpret_cast<unsigned int *>(map.MyGlobalElements());
86  add_indices(indices, indices + n_indices);
87  }
88  compress();
89 }
90 
91 # endif
92 
93 #endif // ifdef DEAL_II_WITH_TRILINOS
94 
95 
96 
97 void
99 {
100  // we will, in the following, modify mutable variables. this can only
101  // work in multithreaded applications if we lock the data structures
102  // via a mutex, so that users can call 'const' functions from threads
103  // in parallel (and these 'const' functions can then call compress()
104  // which itself calls the current function)
105  std::lock_guard<std::mutex> lock(compress_mutex);
106 
107  // see if any of the contiguous ranges can be merged. do not use
108  // std::vector::erase in-place as it is quadratic in the number of
109  // ranges. since the ranges are sorted by their first index, determining
110  // overlap isn't all that hard
111  std::vector<Range>::iterator store = ranges.begin();
112  for (std::vector<Range>::iterator i = ranges.begin(); i != ranges.end();)
113  {
114  std::vector<Range>::iterator next = i;
115  ++next;
116 
117  size_type first_index = i->begin;
118  size_type last_index = i->end;
119 
120  // see if we can merge any of the following ranges
121  while (next != ranges.end() && (next->begin <= last_index))
122  {
123  last_index = std::max(last_index, next->end);
124  ++next;
125  }
126  i = next;
127 
128  // store the new range in the slot we last occupied
129  *store = Range(first_index, last_index);
130  ++store;
131  }
132  // use a compact array with exactly the right amount of storage
133  if (store != ranges.end())
134  {
135  std::vector<Range> new_ranges(ranges.begin(), store);
136  ranges.swap(new_ranges);
137  }
138 
139  // now compute indices within set and the range with most elements
140  size_type next_index = 0, largest_range_size = 0;
141  for (std::vector<Range>::iterator i = ranges.begin(); i != ranges.end(); ++i)
142  {
143  Assert(i->begin < i->end, ExcInternalError());
144 
145  i->nth_index_in_set = next_index;
146  next_index += (i->end - i->begin);
147  if (i->end - i->begin > largest_range_size)
148  {
149  largest_range_size = i->end - i->begin;
150  largest_range = i - ranges.begin();
151  }
152  }
153  is_compressed = true;
154 
155  // check that next_index is correct. needs to be after the previous
156  // statement because we otherwise will get into an endless loop
157  Assert(next_index == n_elements(), ExcInternalError());
158 }
159 
160 
161 
162 #ifndef DOXYGEN
163 IndexSet
164 IndexSet::operator&(const IndexSet &is) const
165 {
166  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
167 
168  compress();
169  is.compress();
170 
171  std::vector<Range>::const_iterator r1 = ranges.begin(),
172  r2 = is.ranges.begin();
173  IndexSet result(size());
174 
175  while ((r1 != ranges.end()) && (r2 != is.ranges.end()))
176  {
177  // if r1 and r2 do not overlap at all, then move the pointer that sits
178  // to the left of the other up by one
179  if (r1->end <= r2->begin)
180  ++r1;
181  else if (r2->end <= r1->begin)
182  ++r2;
183  else
184  {
185  // the ranges must overlap somehow
186  Assert(((r1->begin <= r2->begin) && (r1->end > r2->begin)) ||
187  ((r2->begin <= r1->begin) && (r2->end > r1->begin)),
188  ExcInternalError());
189 
190  // add the overlapping range to the result
191  result.add_range(std::max(r1->begin, r2->begin),
192  std::min(r1->end, r2->end));
193 
194  // now move that iterator that ends earlier one up. note that it has
195  // to be this one because a subsequent range may still have a chance
196  // of overlapping with the range that ends later
197  if (r1->end <= r2->end)
198  ++r1;
199  else
200  ++r2;
201  }
202  }
203 
204  result.compress();
205  return result;
206 }
207 #endif
208 
209 
210 
211 IndexSet
213 {
214  Assert(begin <= end,
215  ExcMessage("End index needs to be larger or equal to begin index!"));
216  Assert(end <= size(), ExcMessage("Given range exceeds index set dimension"));
217 
218  IndexSet result(end - begin);
219  std::vector<Range>::const_iterator r1 = ranges.begin();
220 
221  while (r1 != ranges.end())
222  {
223  if ((r1->end > begin) && (r1->begin < end))
224  {
225  result.add_range(std::max(r1->begin, begin) - begin,
226  std::min(r1->end, end) - begin);
227  }
228  else if (r1->begin >= end)
229  break;
230 
231  ++r1;
232  }
233 
234  result.compress();
235  return result;
236 }
237 
238 std::vector<IndexSet>
240  const std::vector<types::global_dof_index> &n_indices_per_block) const
241 {
242  std::vector<IndexSet> partitioned;
243  const unsigned int n_blocks = n_indices_per_block.size();
244 
245  partitioned.reserve(n_blocks);
246  types::global_dof_index start = 0;
247  for (const auto n_block_indices : n_indices_per_block)
248  {
249  partitioned.push_back(this->get_view(start, start + n_block_indices));
250  start += n_block_indices;
251  }
252 
253 #ifdef DEBUG
255  for (const auto &partition : partitioned)
256  {
257  sum += partition.size();
258  }
259  AssertDimension(sum, this->size());
260 #endif
261 
262  return partitioned;
263 }
264 
265 void
267 {
268  compress();
269  other.compress();
270  is_compressed = false;
271 
272 
273  // we save all new ranges to our IndexSet in an temporary vector and
274  // add all of them in one go at the end.
275  std::vector<Range> new_ranges;
276 
277  std::vector<Range>::iterator own_it = ranges.begin();
278  std::vector<Range>::iterator other_it = other.ranges.begin();
279 
280  while (own_it != ranges.end() && other_it != other.ranges.end())
281  {
282  // advance own iterator until we get an overlap
283  if (own_it->end <= other_it->begin)
284  {
285  new_ranges.push_back(*own_it);
286  ++own_it;
287  continue;
288  }
289  // we are done with other_it, so advance
290  if (own_it->begin >= other_it->end)
291  {
292  ++other_it;
293  continue;
294  }
295 
296  // Now own_it and other_it overlap. First save the part of own_it that
297  // is before other_it (if not empty).
298  if (own_it->begin < other_it->begin)
299  {
300  Range r(own_it->begin, other_it->begin);
301  r.nth_index_in_set = 0; // fix warning of unused variable
302  new_ranges.push_back(r);
303  }
304  // change own_it to the sub range behind other_it. Do not delete own_it
305  // in any case. As removal would invalidate iterators, we just shrink
306  // the range to an empty one.
307  own_it->begin = other_it->end;
308  if (own_it->begin > own_it->end)
309  {
310  own_it->begin = own_it->end;
311  ++own_it;
312  }
313 
314  // continue without advancing iterators, the right one will be advanced
315  // next.
316  }
317 
318  // make sure to take over the remaining ranges
319  for (; own_it != ranges.end(); ++own_it)
320  new_ranges.push_back(*own_it);
321 
322  ranges.clear();
323 
324  // done, now add the temporary ranges
325  const std::vector<Range>::iterator end = new_ranges.end();
326  for (std::vector<Range>::iterator it = new_ranges.begin(); it != end; ++it)
327  add_range(it->begin, it->end);
328 
329  compress();
330 }
331 
332 
333 
334 IndexSet
336 {
337  IndexSet set(this->size() * other.size());
338  for (const auto el : *this)
339  set.add_indices(other, el * other.size());
340  set.compress();
341  return set;
342 }
343 
344 
345 
348 {
349  Assert(is_empty() == false,
350  ExcMessage(
351  "pop_back() failed, because this IndexSet contains no entries."));
352 
353  const size_type index = ranges.back().end - 1;
354  --ranges.back().end;
355 
356  if (ranges.back().begin == ranges.back().end)
357  ranges.pop_back();
358 
359  return index;
360 }
361 
362 
363 
366 {
367  Assert(is_empty() == false,
368  ExcMessage(
369  "pop_front() failed, because this IndexSet contains no entries."));
370 
371  const size_type index = ranges.front().begin;
372  ++ranges.front().begin;
373 
374  if (ranges.front().begin == ranges.front().end)
375  ranges.erase(ranges.begin());
376 
377  // We have to set this in any case, because nth_index_in_set is no longer
378  // up to date for all but the first range
379  is_compressed = false;
380 
381  return index;
382 }
383 
384 
385 
386 void
388 {
389  ranges.insert(Utilities::lower_bound(ranges.begin(), ranges.end(), new_range),
390  new_range);
391 }
392 
393 
394 
395 void
397  boost::container::small_vector<std::pair<size_type, size_type>, 200>
398  & tmp_ranges,
399  const bool ranges_are_sorted)
400 {
401  if (!ranges_are_sorted)
402  std::sort(tmp_ranges.begin(), tmp_ranges.end());
403 
404  // if we have many ranges, we first construct a temporary index set (where
405  // we add ranges in a consecutive way, so fast), otherwise, we work with
406  // add_range(). the number 9 is chosen heuristically given the fact that
407  // there are typically up to 8 independent ranges when adding the degrees of
408  // freedom on a 3D cell or 9 when adding degrees of freedom of faces. if
409  // doing cell-by-cell additions, we want to avoid repeated calls to
410  // IndexSet::compress() which gets called upon merging two index sets, so we
411  // want to be in the other branch then.
412  if (tmp_ranges.size() > 9)
413  {
414  IndexSet tmp_set(size());
415  tmp_set.ranges.reserve(tmp_ranges.size());
416  for (const auto &i : tmp_ranges)
417  tmp_set.add_range(i.first, i.second);
418 
419  // Case if we have zero or just one range: Add into the other set with
420  // its indices, as that is cheaper
421  if (this->ranges.size() <= 1)
422  {
423  if (this->ranges.size() == 1)
424  tmp_set.add_range(ranges[0].begin, ranges[0].end);
425  std::swap(*this, tmp_set);
426  }
427  else
428  this->add_indices(tmp_set);
429  }
430  else
431  for (const auto &i : tmp_ranges)
432  add_range(i.first, i.second);
433 }
434 
435 
436 
437 void
438 IndexSet::add_indices(const IndexSet &other, const size_type offset)
439 {
440  if ((this == &other) && (offset == 0))
441  return;
442 
443  if (other.ranges.size() != 0)
444  {
445  AssertIndexRange(other.ranges.back().end - 1, index_space_size);
446  }
447 
448  compress();
449  other.compress();
450 
451  std::vector<Range>::const_iterator r1 = ranges.begin(),
452  r2 = other.ranges.begin();
453 
454  std::vector<Range> new_ranges;
455  // just get the start and end of the ranges right in this method, everything
456  // else will be done in compress()
457  while (r1 != ranges.end() || r2 != other.ranges.end())
458  {
459  // the two ranges do not overlap or we are at the end of one of the
460  // ranges
461  if (r2 == other.ranges.end() ||
462  (r1 != ranges.end() && r1->end < (r2->begin + offset)))
463  {
464  new_ranges.push_back(*r1);
465  ++r1;
466  }
467  else if (r1 == ranges.end() || (r2->end + offset) < r1->begin)
468  {
469  new_ranges.emplace_back(r2->begin + offset, r2->end + offset);
470  ++r2;
471  }
472  else
473  {
474  // ok, we do overlap, so just take the combination of the current
475  // range (do not bother to merge with subsequent ranges)
476  Range next(std::min(r1->begin, r2->begin + offset),
477  std::max(r1->end, r2->end + offset));
478  new_ranges.push_back(next);
479  ++r1;
480  ++r2;
481  }
482  }
483  ranges.swap(new_ranges);
484 
485  is_compressed = false;
486  compress();
487 }
488 
489 
490 
491 void
492 IndexSet::write(std::ostream &out) const
493 {
494  compress();
495  out << size() << " ";
496  out << ranges.size() << std::endl;
497  std::vector<Range>::const_iterator r = ranges.begin();
498  for (; r != ranges.end(); ++r)
499  {
500  out << r->begin << " " << r->end << std::endl;
501  }
502 }
503 
504 
505 
506 void
507 IndexSet::read(std::istream &in)
508 {
509  AssertThrow(in.fail() == false, ExcIO());
510 
511  size_type s;
512  unsigned int n_ranges;
513 
514  in >> s >> n_ranges;
515  ranges.clear();
516  set_size(s);
517  for (unsigned int i = 0; i < n_ranges; ++i)
518  {
519  AssertThrow(in.fail() == false, ExcIO());
520 
521  size_type b, e;
522  in >> b >> e;
523  add_range(b, e);
524  }
525 }
526 
527 
528 void
529 IndexSet::block_write(std::ostream &out) const
530 {
531  AssertThrow(out.fail() == false, ExcIO());
532  out.write(reinterpret_cast<const char *>(&index_space_size),
533  sizeof(index_space_size));
534  std::size_t n_ranges = ranges.size();
535  out.write(reinterpret_cast<const char *>(&n_ranges), sizeof(n_ranges));
536  if (ranges.empty() == false)
537  out.write(reinterpret_cast<const char *>(&*ranges.begin()),
538  ranges.size() * sizeof(Range));
539  AssertThrow(out.fail() == false, ExcIO());
540 }
541 
542 void
543 IndexSet::block_read(std::istream &in)
544 {
545  size_type size;
546  std::size_t n_ranges;
547  in.read(reinterpret_cast<char *>(&size), sizeof(size));
548  in.read(reinterpret_cast<char *>(&n_ranges), sizeof(n_ranges));
549  // we have to clear ranges first
550  ranges.clear();
551  set_size(size);
552  ranges.resize(n_ranges, Range(0, 0));
553  if (n_ranges != 0u)
554  in.read(reinterpret_cast<char *>(&*ranges.begin()),
555  ranges.size() * sizeof(Range));
556 
557  do_compress(); // needed so that largest_range can be recomputed
558 }
559 
560 
561 
562 bool
564 {
565  // get the element after which we would have to insert a range that
566  // consists of all elements from this element to the end of the index
567  // range plus one. after this call we know that if p!=end() then
568  // p->begin<=index unless there is no such range at all
569  //
570  // if the searched for element is an element of this range, then we're
571  // done. otherwise, the element can't be in one of the following ranges
572  // because otherwise p would be a different iterator
573  //
574  // since we already know the position relative to the largest range (we
575  // called compress!), we can perform the binary search on ranges with
576  // lower/higher number compared to the largest range
577  std::vector<Range>::const_iterator p = std::upper_bound(
578  ranges.begin() +
579  (index < ranges[largest_range].begin ? 0 : largest_range + 1),
580  index < ranges[largest_range].begin ? ranges.begin() + largest_range :
581  ranges.end(),
582  Range(index, size() + 1));
583 
584  if (p == ranges.begin())
585  return ((index >= p->begin) && (index < p->end));
586 
587  Assert((p == ranges.end()) || (p->begin > index), ExcInternalError());
588 
589  // now move to that previous range
590  --p;
591  Assert(p->begin <= index, ExcInternalError());
592 
593  return (p->end > index);
594 }
595 
596 
597 
600 {
601  // find out which chunk the local index n belongs to by using a binary
602  // search. the comparator is based on the end of the ranges.
603  Range r(n, n + 1);
604  r.nth_index_in_set = n;
605 
606  const std::vector<Range>::const_iterator p = Utilities::lower_bound(
607  ranges.begin(), ranges.end(), r, Range::nth_index_compare);
608 
609  Assert(p != ranges.end(), ExcInternalError());
610  return p->begin + (n - p->nth_index_in_set);
611 }
612 
613 
614 
617 {
618  // we could try to use the main range for splitting up the search range, but
619  // since we only come here when the largest range did not contain the index,
620  // there is little gain from doing a first step manually.
621  Range r(n, n);
622  std::vector<Range>::const_iterator p =
624 
625  // if n is not in this set
626  if (p == ranges.end() || p->end == n || p->begin > n)
628 
629  Assert(p != ranges.end(), ExcInternalError());
630  Assert(p->begin <= n, ExcInternalError());
631  Assert(n < p->end, ExcInternalError());
632  return (n - p->begin) + p->nth_index_in_set;
633 }
634 
635 
636 
639 {
640  compress();
642 
643  if (ranges.empty())
644  return end();
645 
646  std::vector<Range>::const_iterator main_range =
647  ranges.begin() + largest_range;
648 
650  // This optimization makes the bounds for lower_bound smaller by checking
651  // the largest range first.
652  std::vector<Range>::const_iterator range_begin, range_end;
653  if (global_index < main_range->begin)
654  {
655  range_begin = ranges.begin();
656  range_end = main_range;
657  }
658  else
659  {
660  range_begin = main_range;
661  range_end = ranges.end();
662  }
663 
664  // This will give us the first range p=[a,b[ with b>=global_index using
665  // a binary search
666  const std::vector<Range>::const_iterator p =
667  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
668 
669  // We couldn't find a range, which means we have no range that contains
670  // global_index and also no range behind it, meaning we need to return end().
671  if (p == ranges.end())
672  return end();
673 
674  // Finally, we can have two cases: Either global_index is not in [a,b[,
675  // which means we need to return an iterator to a because global_index, ...,
676  // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
677  // [a,b[ and we will return an iterator pointing directly at global_index
678  // (else branch).
679  if (global_index < p->begin)
680  return {this, static_cast<size_type>(p - ranges.begin()), p->begin};
681  else
682  return {this, static_cast<size_type>(p - ranges.begin()), global_index};
683 }
684 
685 
686 
687 void
688 IndexSet::fill_index_vector(std::vector<size_type> &indices) const
689 {
690  compress();
691 
692  indices.clear();
693  indices.reserve(n_elements());
694 
695  for (const auto &range : ranges)
696  for (size_type entry = range.begin; entry < range.end; ++entry)
697  indices.push_back(entry);
698 
699  Assert(indices.size() == n_elements(), ExcInternalError());
700 }
701 
702 
703 
704 #ifdef DEAL_II_WITH_TRILINOS
705 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
706 
707 Tpetra::Map<int, types::global_dof_index>
708 IndexSet::make_tpetra_map(const MPI_Comm &communicator,
709  const bool overlapping) const
710 {
711  compress();
712  (void)communicator;
713 
714 # ifdef DEBUG
715  if (!overlapping)
716  {
718  Utilities::MPI::sum(n_elements(), communicator);
720  ExcMessage("You are trying to create an Tpetra::Map object "
721  "that partitions elements of an index set "
722  "between processors. However, the union of the "
723  "index sets on different processors does not "
724  "contain all indices exactly once: the sum of "
725  "the number of entries the various processors "
726  "want to store locally is " +
728  " whereas the total size of the object to be "
729  "allocated is " +
730  std::to_string(size()) +
731  ". In other words, there are "
732  "either indices that are not spoken for "
733  "by any processor, or there are indices that are "
734  "claimed by multiple processors."));
735  }
736 # endif
737 
738  // Find out if the IndexSet is ascending and 1:1. This corresponds to a
739  // linear Tpetra::Map. Overlapping IndexSets are never 1:1.
740  const bool linear =
741  overlapping ? false : is_ascending_and_one_to_one(communicator);
742  if (linear)
743  return Tpetra::Map<int, types::global_dof_index>(
744  size(),
745  n_elements(),
746  0,
747 # ifdef DEAL_II_WITH_MPI
748  Teuchos::rcp(new Teuchos::MpiComm<int>(communicator))
749 # else
750  Teuchos::rcp(new Teuchos::Comm<int>())
751 # endif
752  );
753  else
754  {
755  std::vector<size_type> indices;
756  fill_index_vector(indices);
757  std::vector<types::global_dof_index> int_indices(indices.size());
758  std::copy(indices.begin(), indices.end(), int_indices.begin());
759  const Teuchos::ArrayView<types::global_dof_index> arr_view(int_indices);
760  return Tpetra::Map<int, types::global_dof_index>(
761  size(),
762  arr_view,
763  0,
764 # ifdef DEAL_II_WITH_MPI
765  Teuchos::rcp(new Teuchos::MpiComm<int>(communicator))
766 # else
767  Teuchos::rcp(new Teuchos::Comm<int>())
768 # endif
769  );
770  }
771 }
772 # endif
773 
774 
775 
776 Epetra_Map
777 IndexSet::make_trilinos_map(const MPI_Comm &communicator,
778  const bool overlapping) const
779 {
780  compress();
781  (void)communicator;
782 
783 # ifdef DEBUG
784  if (!overlapping)
785  {
787  Utilities::MPI::sum(n_elements(), communicator);
789  ExcMessage("You are trying to create an Epetra_Map object "
790  "that partitions elements of an index set "
791  "between processors. However, the union of the "
792  "index sets on different processors does not "
793  "contain all indices exactly once: the sum of "
794  "the number of entries the various processors "
795  "want to store locally is " +
797  " whereas the total size of the object to be "
798  "allocated is " +
799  std::to_string(size()) +
800  ". In other words, there are "
801  "either indices that are not spoken for "
802  "by any processor, or there are indices that are "
803  "claimed by multiple processors."));
804  }
805 # endif
806 
807  // Find out if the IndexSet is ascending and 1:1. This corresponds to a
808  // linear EpetraMap. Overlapping IndexSets are never 1:1.
809  const bool linear =
810  overlapping ? false : is_ascending_and_one_to_one(communicator);
811 
812  if (linear)
813  return Epetra_Map(TrilinosWrappers::types::int_type(size()),
815  0,
816 # ifdef DEAL_II_WITH_MPI
817  Epetra_MpiComm(communicator)
818 # else
819  Epetra_SerialComm()
820 # endif
821  );
822  else
823  {
824  std::vector<size_type> indices;
825  fill_index_vector(indices);
826  return Epetra_Map(
829  (n_elements() > 0 ?
830  reinterpret_cast<TrilinosWrappers::types::int_type *>(
831  indices.data()) :
832  nullptr),
833  0,
834 # ifdef DEAL_II_WITH_MPI
835  Epetra_MpiComm(communicator)
836 # else
837  Epetra_SerialComm()
838 # endif
839  );
840  }
841 }
842 #endif
843 
844 
845 
846 bool
847 IndexSet::is_ascending_and_one_to_one(const MPI_Comm &communicator) const
848 {
849  // If the sum of local elements does not add up to the total size,
850  // the IndexSet can't be complete.
852  Utilities::MPI::sum(n_elements(), communicator);
853  if (n_global_elements != size())
854  return false;
855 
856  if (n_global_elements == 0)
857  return true;
858 
859 #ifdef DEAL_II_WITH_MPI
860  // Non-contiguous IndexSets can't be linear.
861  const bool all_contiguous =
862  (Utilities::MPI::min(is_contiguous() ? 1 : 0, communicator) == 1);
863  if (!all_contiguous)
864  return false;
865 
866  bool is_globally_ascending = true;
867  // we know that there is only one interval
868  types::global_dof_index first_local_dof = (n_elements() > 0) ?
869  *(begin_intervals()->begin()) :
871 
872  const unsigned int my_rank = Utilities::MPI::this_mpi_process(communicator);
873  const std::vector<types::global_dof_index> global_dofs =
874  Utilities::MPI::gather(communicator, first_local_dof, 0);
875 
876  if (my_rank == 0)
877  {
878  // find out if the received std::vector is ascending
879  types::global_dof_index index = 0;
880  while (global_dofs[index] == numbers::invalid_dof_index)
881  ++index;
882  types::global_dof_index old_dof = global_dofs[index++];
883  for (; index < global_dofs.size(); ++index)
884  {
885  const types::global_dof_index new_dof = global_dofs[index];
886  if (new_dof != numbers::invalid_dof_index)
887  {
888  if (new_dof <= old_dof)
889  {
890  is_globally_ascending = false;
891  break;
892  }
893  else
894  old_dof = new_dof;
895  }
896  }
897  }
898 
899  // now broadcast the result
900  int is_ascending = is_globally_ascending ? 1 : 0;
901  int ierr = MPI_Bcast(&is_ascending, 1, MPI_INT, 0, communicator);
902  AssertThrowMPI(ierr);
903 
904  return (is_ascending == 1);
905 #else
906  return true;
907 #endif // DEAL_II_WITH_MPI
908 }
909 
910 
911 
912 std::size_t
914 {
918  sizeof(compress_mutex));
919 }
920 
921 
922 
ElementIterator begin() const
Definition: index_set.h:1115
bool is_element_binary_search(const size_type local_index) const
Definition: index_set.cc:563
size_type index_within_set_binary_search(const size_type global_index) const
Definition: index_set.cc:616
size_type largest_range
Definition: index_set.h:975
bool is_contiguous() const
Definition: index_set.h:1757
void do_compress() const
Definition: index_set.cc:98
ElementIterator at(const size_type global_index) const
Definition: index_set.cc:638
size_type size() const
Definition: index_set.h:1626
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:688
bool is_empty() const
Definition: index_set.h:1766
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
Definition: index_set.cc:239
size_type n_elements() const
Definition: index_set.h:1774
void add_range_lower_bound(const Range &range)
Definition: index_set.cc:387
ElementIterator begin() const
Definition: index_set.h:1560
size_type pop_front()
Definition: index_set.cc:365
void set_size(const size_type size)
Definition: index_set.h:1614
void read(std::istream &in)
Definition: index_set.cc:507
Tpetra::Map< int, types::global_dof_index > make_tpetra_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:708
IndexSet tensor_product(const IndexSet &other) const
Definition: index_set.cc:335
void write(std::ostream &out) const
Definition: index_set.cc:492
void block_read(std::istream &in)
Definition: index_set.cc:543
bool is_compressed
Definition: index_set.h:958
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:396
IntervalIterator begin_intervals() const
Definition: index_set.h:1581
std::vector< Range > ranges
Definition: index_set.h:948
void subtract_set(const IndexSet &other)
Definition: index_set.cc:266
ElementIterator end() const
Definition: index_set.h:1572
Threads::Mutex compress_mutex
Definition: index_set.h:981
size_type index_space_size
Definition: index_set.h:964
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:777
void block_write(std::ostream &out) const
Definition: index_set.cc:529
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:212
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1662
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:847
std::size_t memory_consumption() const
Definition: index_set.cc:913
size_type nth_index_in_set_binary_search(const size_type local_index) const
Definition: index_set.cc:599
void compress() const
Definition: index_set.h:1634
size_type pop_back()
Definition: index_set.cc:347
types::global_dof_index size_type
Definition: index_set.h:75
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1692
IndexSet operator&(const IndexSet &is) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
__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:1501
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1818
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
types::global_dof_index size_type
Definition: cuda_kernels.h:45
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< T >::value, std::size_t > memory_consumption(const T &t)
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
std::string to_string(const T &t)
Definition: patterns.h:2394
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)
std::vector< T > gather(const MPI_Comm &comm, const T &object_to_send, const unsigned int root_process=0)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:157
T min(const T &t, const MPI_Comm &mpi_communicator)
T sum(const T &t, const MPI_Comm &mpi_communicator)
std::string compress(const std::string &input)
Definition: utilities.cc:392
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:999
void copy(const T *begin, const T *end, U *dest)
static const unsigned int invalid_unsigned_int
Definition: types.h:206
const types::global_dof_index invalid_dof_index
Definition: types.h:226
unsigned int global_dof_index
Definition: types.h:81
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:906
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:912
size_type nth_index_in_set
Definition: index_set.h:876