Reference documentation for deal.II version Git cd8b20072f 2020-09-24 16:28:23 -0400
\(\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.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_index_set_h
17 #define dealii_index_set_h
18 
19 #include <deal.II/base/config.h>
20 
23 #include <deal.II/base/utilities.h>
24 
25 #include <boost/serialization/vector.hpp>
26 
27 #include <algorithm>
28 #include <iterator>
29 #include <vector>
30 
31 
32 #ifdef DEAL_II_WITH_TRILINOS
33 # include <Epetra_Map.h>
34 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
35 # include <Tpetra_Map.hpp>
36 # endif
37 #endif
38 
39 #if defined(DEAL_II_WITH_MPI) || defined(DEAL_II_WITH_PETSC)
40 # include <mpi.h>
41 #else
42 using MPI_Comm = int;
43 # ifndef MPI_COMM_WORLD
44 # define MPI_COMM_WORLD 0
45 # endif
46 #endif
47 
49 
72 class IndexSet
73 {
74 public:
75  // forward declarations:
76  class ElementIterator;
77  class IntervalIterator;
78 
84 
102  using value_type = signed int;
103 
104 
108  IndexSet();
109 
113  explicit IndexSet(const size_type size);
114 
118  IndexSet(const IndexSet &) = default;
119 
123  IndexSet &
124  operator=(const IndexSet &) = default;
125 
130  IndexSet(IndexSet &&is) noexcept;
131 
136  IndexSet &
137  operator=(IndexSet &&is) noexcept;
138 
139 #ifdef DEAL_II_WITH_TRILINOS
140 
143  explicit IndexSet(const Epetra_BlockMap &map);
144 #endif
145 
150  void
151  clear();
152 
159  void
160  set_size(const size_type size);
161 
169  size_type
170  size() const;
171 
178  void
179  add_range(const size_type begin, const size_type end);
180 
184  void
185  add_index(const size_type index);
186 
196  template <typename ForwardIterator>
197  void
198  add_indices(const ForwardIterator &begin, const ForwardIterator &end);
199 
215  void
216  add_indices(const IndexSet &other, const size_type offset = 0);
217 
221  bool
222  is_element(const size_type index) const;
223 
228  bool
229  is_contiguous() const;
230 
235  bool
236  is_empty() const;
237 
247  bool
248  is_ascending_and_one_to_one(const MPI_Comm &communicator) const;
249 
253  size_type
254  n_elements() const;
255 
261  size_type
262  nth_index_in_set(const size_type local_index) const;
263 
270  size_type
272 
281  unsigned int
282  n_intervals() const;
283 
295  size_type
297 
302  void
303  compress() const;
304 
310  bool
311  operator==(const IndexSet &is) const;
312 
318  bool
319  operator!=(const IndexSet &is) const;
320 
327  IndexSet operator&(const IndexSet &is) const;
328 
342  IndexSet
343  get_view(const size_type begin, const size_type end) const;
344 
350  std::vector<IndexSet>
352  const std::vector<types::global_dof_index> &n_indices_per_block) const;
353 
359  void
360  subtract_set(const IndexSet &other);
361 
378  IndexSet
379  tensor_product(const IndexSet &other) const;
380 
385  size_type
386  pop_back();
387 
392  size_type
393  pop_front();
394 
398  void
399  fill_index_vector(std::vector<size_type> &indices) const;
400 
413  template <typename VectorType>
414  void
415  fill_binary_vector(VectorType &vector) const;
416 
421  template <class StreamType>
422  void
423  print(StreamType &out) const;
424 
429  void
430  write(std::ostream &out) const;
431 
436  void
437  read(std::istream &in);
438 
443  void
444  block_write(std::ostream &out) const;
445 
450  void
451  block_read(std::istream &in);
452 
453 #ifdef DEAL_II_WITH_TRILINOS
454 
482  Epetra_Map
483  make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
484  const bool overlapping = false) const;
485 
486 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
487  Tpetra::Map<int, types::global_dof_index>
488  make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
489  const bool overlapping = false) const;
490 # endif
491 #endif
492 
493 
498  std::size_t
499  memory_consumption() const;
500 
502  size_type,
503  << "The global index " << arg1
504  << " is not an element of this set.");
505 
510  template <class Archive>
511  void
512  serialize(Archive &ar, const unsigned int version);
513 
514 
526  {
527  public:
532  IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
533 
537  explicit IntervalAccessor(const IndexSet *idxset);
538 
542  size_type
543  n_elements() const;
544 
548  bool
549  is_valid() const;
550 
555  begin() const;
556 
562  end() const;
563 
567  size_type
568  last() const;
569 
570  private:
574  IntervalAccessor(const IntervalAccessor &other);
579  operator=(const IntervalAccessor &other);
580 
584  bool
585  operator==(const IntervalAccessor &other) const;
589  bool
590  operator<(const IntervalAccessor &other) const;
595  void
596  advance();
601 
607 
608  friend class IntervalIterator;
609  };
610 
616  {
617  public:
622  IntervalIterator(const IndexSet *idxset, const size_type range_idx);
623 
627  explicit IntervalIterator(const IndexSet *idxset);
628 
633 
637  IntervalIterator(const IntervalIterator &other) = default;
638 
643  operator=(const IntervalIterator &other) = default;
644 
649  operator++();
650 
655  operator++(int);
656 
660  const IntervalAccessor &operator*() const;
661 
665  const IntervalAccessor *operator->() const;
666 
670  bool
671  operator==(const IntervalIterator &) const;
672 
676  bool
677  operator!=(const IntervalIterator &) const;
678 
682  bool
683  operator<(const IntervalIterator &) const;
684 
691  int
692  operator-(const IntervalIterator &p) const;
693 
699  using iterator_category = std::forward_iterator_tag;
701  using difference_type = std::ptrdiff_t;
704 
705  private:
710  };
711 
717  {
718  public:
723  ElementIterator(const IndexSet *idxset,
724  const size_type range_idx,
725  const size_type index);
726 
730  explicit ElementIterator(const IndexSet *idxset);
731 
736  size_type operator*() const;
737 
741  bool
742  is_valid() const;
743 
748  operator++();
749 
754  operator++(int);
755 
759  bool
760  operator==(const ElementIterator &) const;
761 
765  bool
766  operator!=(const ElementIterator &) const;
767 
771  bool
772  operator<(const ElementIterator &) const;
773 
782  std::ptrdiff_t
783  operator-(const ElementIterator &p) const;
784 
790  using iterator_category = std::forward_iterator_tag;
792  using difference_type = std::ptrdiff_t;
793  using pointer = size_type *;
794  using reference = size_type &;
795 
796  private:
800  void
801  advance();
802 
815  };
816 
822  begin() const;
823 
839  at(const size_type global_index) const;
840 
846  end() const;
847 
852  begin_intervals() const;
853 
859  end_intervals() const;
860 
865 private:
874  struct Range
875  {
878 
880 
889  Range();
890 
898  Range(const size_type i1, const size_type i2);
899 
900  friend inline bool
901  operator<(const Range &range_1, const Range &range_2)
902  {
903  return (
904  (range_1.begin < range_2.begin) ||
905  ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
906  }
907 
908  static bool
910  {
911  return x.end < y.end;
912  }
913 
914  static bool
916  {
917  return (x.nth_index_in_set + (x.end - x.begin) <
918  y.nth_index_in_set + (y.end - y.begin));
919  }
920 
921  friend inline bool
922  operator==(const Range &range_1, const Range &range_2)
923  {
924  return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
925  }
926 
927  static std::size_t
929  {
930  return sizeof(Range);
931  }
932 
937  template <class Archive>
938  void
939  serialize(Archive &ar, const unsigned int version);
940  };
941 
950  mutable std::vector<Range> ranges;
951 
960  mutable bool is_compressed;
961 
967 
978 
984 
988  void
989  do_compress() const;
990 };
991 
992 
1010 inline IndexSet
1012 {
1013  IndexSet is(N);
1014  is.add_range(0, N);
1015  is.compress();
1016  return is;
1017 }
1018 
1019 /* ------------------ inline functions ------------------ */
1020 
1021 
1022 /* IntervalAccessor */
1023 
1025  const IndexSet * idxset,
1027  : index_set(idxset)
1028  , range_idx(range_idx)
1029 {
1030  Assert(range_idx < idxset->n_intervals(),
1031  ExcInternalError("Invalid range index"));
1032 }
1033 
1034 
1035 
1037  : index_set(idxset)
1039 {}
1040 
1041 
1042 
1044  const IndexSet::IntervalAccessor &other)
1045  : index_set(other.index_set)
1046  , range_idx(other.range_idx)
1047 {
1049  ExcMessage("invalid iterator"));
1050 }
1051 
1052 
1053 
1054 inline IndexSet::size_type
1056 {
1057  Assert(is_valid(), ExcMessage("invalid iterator"));
1058  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1059 }
1060 
1061 
1062 
1063 inline bool
1065 {
1066  return index_set != nullptr && range_idx < index_set->n_intervals();
1067 }
1068 
1069 
1070 
1073 {
1074  Assert(is_valid(), ExcMessage("invalid iterator"));
1075  return {index_set, range_idx, index_set->ranges[range_idx].begin};
1076 }
1077 
1078 
1079 
1082 {
1083  Assert(is_valid(), ExcMessage("invalid iterator"));
1084 
1085  // point to first index in next interval unless we are the last interval.
1086  if (range_idx < index_set->ranges.size() - 1)
1087  return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1088  else
1089  return index_set->end();
1090 }
1091 
1092 
1093 
1094 inline IndexSet::size_type
1096 {
1097  Assert(is_valid(), ExcMessage("invalid iterator"));
1098 
1099  return index_set->ranges[range_idx].end - 1;
1100 }
1101 
1102 
1103 
1106 {
1107  index_set = other.index_set;
1108  range_idx = other.range_idx;
1110  ExcMessage("invalid iterator"));
1111  return *this;
1112 }
1113 
1114 
1115 
1116 inline bool
1119 {
1120  Assert(index_set == other.index_set,
1121  ExcMessage(
1122  "Can not compare accessors pointing to different IndexSets"));
1123  return range_idx == other.range_idx;
1124 }
1125 
1126 
1127 
1128 inline bool
1131 {
1132  Assert(index_set == other.index_set,
1133  ExcMessage(
1134  "Can not compare accessors pointing to different IndexSets"));
1135  return range_idx < other.range_idx;
1136 }
1137 
1138 
1139 
1140 inline void
1142 {
1143  Assert(
1144  is_valid(),
1145  ExcMessage(
1146  "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1147  ++range_idx;
1148 
1149  // set ourselves to invalid if we walk off the end
1150  if (range_idx >= index_set->ranges.size())
1152 }
1153 
1154 
1155 /* IntervalIterator */
1156 
1158  const IndexSet * idxset,
1160  : accessor(idxset, range_idx)
1161 {}
1162 
1163 
1164 
1166  : accessor(nullptr)
1167 {}
1168 
1169 
1170 
1172  : accessor(idxset)
1173 {}
1174 
1175 
1176 
1179 {
1180  accessor.advance();
1181  return *this;
1182 }
1183 
1184 
1185 
1188 {
1189  const IndexSet::IntervalIterator iter = *this;
1190  accessor.advance();
1191  return iter;
1192 }
1193 
1194 
1195 
1197  operator*() const
1198 {
1199  return accessor;
1200 }
1201 
1202 
1203 
1206 {
1207  return &accessor;
1208 }
1209 
1210 
1211 
1212 inline bool
1215 {
1216  return accessor == other.accessor;
1217 }
1218 
1219 
1220 
1221 inline bool
1224 {
1225  return !(*this == other);
1226 }
1227 
1228 
1229 
1230 inline bool
1233 {
1234  return accessor < other.accessor;
1235 }
1236 
1237 
1238 
1239 inline int
1242 {
1244  ExcMessage(
1245  "Can not compare iterators belonging to different IndexSets"));
1246 
1248  accessor.index_set->ranges.size() :
1250  const size_type rhs =
1252  accessor.index_set->ranges.size() :
1253  other.accessor.range_idx;
1254 
1255  if (lhs > rhs)
1256  return static_cast<int>(lhs - rhs);
1257  else
1258  return -static_cast<int>(rhs - lhs);
1259 }
1260 
1261 
1262 
1263 /* ElementIterator */
1264 
1266  const IndexSet * idxset,
1267  const IndexSet::size_type range_idx,
1268  const IndexSet::size_type index)
1269  : index_set(idxset)
1270  , range_idx(range_idx)
1271  , idx(index)
1272 {
1273  Assert(range_idx < index_set->ranges.size(),
1274  ExcMessage(
1275  "Invalid range index for IndexSet::ElementIterator constructor."));
1276  Assert(
1277  idx >= index_set->ranges[range_idx].begin &&
1278  idx < index_set->ranges[range_idx].end,
1280  "Invalid index argument for IndexSet::ElementIterator constructor."));
1281 }
1282 
1283 
1284 
1286  : index_set(idxset)
1289 {}
1290 
1291 
1292 
1293 inline bool
1295 {
1298  (range_idx < index_set->ranges.size() &&
1299  idx < index_set->ranges[range_idx].end),
1300  ExcInternalError("Invalid ElementIterator state."));
1301 
1302  return (range_idx < index_set->ranges.size() &&
1303  idx < index_set->ranges[range_idx].end);
1304 }
1305 
1306 
1307 
1309 {
1310  Assert(
1311  is_valid(),
1312  ExcMessage(
1313  "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1314  return idx;
1315 }
1316 
1317 
1318 
1319 inline bool
1322 {
1323  Assert(index_set == other.index_set,
1324  ExcMessage(
1325  "Can not compare iterators belonging to different IndexSets"));
1326  return range_idx == other.range_idx && idx == other.idx;
1327 }
1328 
1329 
1330 
1331 inline void
1333 {
1334  Assert(
1335  is_valid(),
1336  ExcMessage(
1337  "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1338  if (idx < index_set->ranges[range_idx].end)
1339  ++idx;
1340  // end of this range?
1341  if (idx == index_set->ranges[range_idx].end)
1342  {
1343  // point to first element in next interval if possible
1344  if (range_idx < index_set->ranges.size() - 1)
1345  {
1346  ++range_idx;
1347  idx = index_set->ranges[range_idx].begin;
1348  }
1349  else
1350  {
1351  // we just fell off the end, set to invalid:
1354  }
1355  }
1356 }
1357 
1358 
1359 
1362 {
1363  advance();
1364  return *this;
1365 }
1366 
1367 
1368 
1371 {
1372  const IndexSet::ElementIterator it = *this;
1373  advance();
1374  return it;
1375 }
1376 
1377 
1378 
1379 inline bool
1382 {
1383  return !(*this == other);
1384 }
1385 
1386 
1387 
1388 inline bool
1391 {
1392  Assert(index_set == other.index_set,
1393  ExcMessage(
1394  "Can not compare iterators belonging to different IndexSets"));
1395  return range_idx < other.range_idx ||
1396  (range_idx == other.range_idx && idx < other.idx);
1397 }
1398 
1399 
1400 
1401 inline std::ptrdiff_t
1404 {
1405  Assert(index_set == other.index_set,
1406  ExcMessage(
1407  "Can not compare iterators belonging to different IndexSets"));
1408  if (*this == other)
1409  return 0;
1410  if (!(*this < other))
1411  return -(other - *this);
1412 
1413  // only other can be equal to end() because of the checks above.
1415 
1416  // Note: we now compute how far advance *this in "*this < other" to get other,
1417  // so we need to return -c at the end.
1418 
1419  // first finish the current range:
1420  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1421 
1422  // now walk in steps of ranges (need to start one behind our current one):
1423  for (size_type range = range_idx + 1;
1424  range < index_set->ranges.size() && range <= other.range_idx;
1425  ++range)
1426  c += index_set->ranges[range].end - index_set->ranges[range].begin;
1427 
1428  Assert(
1429  other.range_idx < index_set->ranges.size() ||
1431  ExcMessage(
1432  "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1433 
1434  // We might have walked too far because we went until the end of
1435  // other.range_idx, so walk backwards to other.idx:
1437  c -= index_set->ranges[other.range_idx].end - other.idx;
1438 
1439  return -c;
1440 }
1441 
1442 
1443 /* Range */
1444 
1449 {}
1450 
1451 
1452 
1453 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1454  : begin(i1)
1455  , end(i2)
1457 {}
1458 
1459 
1460 
1461 /* IndexSet itself */
1462 
1464  : is_compressed(true)
1465  , index_space_size(0)
1467 {}
1468 
1469 
1470 
1472  : is_compressed(true)
1473  , index_space_size(size)
1475 {}
1476 
1477 
1478 
1479 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1480  : ranges(std::move(is.ranges))
1481  , is_compressed(is.is_compressed)
1482  , index_space_size(is.index_space_size)
1483  , largest_range(is.largest_range)
1484 {
1485  is.ranges.clear();
1486  is.is_compressed = true;
1487  is.index_space_size = 0;
1488  is.largest_range = numbers::invalid_unsigned_int;
1489 
1490  compress();
1491 }
1492 
1493 
1494 
1495 inline IndexSet &
1497 {
1498  ranges = std::move(is.ranges);
1499  is_compressed = is.is_compressed;
1500  index_space_size = is.index_space_size;
1501  largest_range = is.largest_range;
1502 
1503  is.ranges.clear();
1504  is.is_compressed = true;
1505  is.index_space_size = 0;
1506  is.largest_range = numbers::invalid_unsigned_int;
1507 
1508  compress();
1509 
1510  return *this;
1511 }
1512 
1513 
1514 
1517 {
1518  compress();
1519  if (ranges.size() > 0)
1520  return {this, 0, ranges[0].begin};
1521  else
1522  return end();
1523 }
1524 
1525 
1526 
1529 {
1530  compress();
1531  AssertIndexRange(global_index, size());
1532 
1533  if (ranges.empty())
1534  return end();
1535 
1536  std::vector<Range>::const_iterator main_range =
1537  ranges.begin() + largest_range;
1538 
1539  Range r(global_index, global_index + 1);
1540  // This optimization makes the bounds for lower_bound smaller by checking
1541  // the largest range first.
1542  std::vector<Range>::const_iterator range_begin, range_end;
1543  if (global_index < main_range->begin)
1544  {
1545  range_begin = ranges.begin();
1546  range_end = main_range;
1547  }
1548  else
1549  {
1550  range_begin = main_range;
1551  range_end = ranges.end();
1552  }
1553 
1554  // This will give us the first range p=[a,b[ with b>=global_index using
1555  // a binary search
1556  const std::vector<Range>::const_iterator p =
1557  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1558 
1559  // We couldn't find a range, which means we have no range that contains
1560  // global_index and also no range behind it, meaning we need to return end().
1561  if (p == ranges.end())
1562  return end();
1563 
1564  // Finally, we can have two cases: Either global_index is not in [a,b[,
1565  // which means we need to return an iterator to a because global_index, ...,
1566  // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
1567  // [a,b[ and we will return an iterator pointing directly at global_index
1568  // (else branch).
1569  if (global_index < p->begin)
1570  return {this, static_cast<size_type>(p - ranges.begin()), p->begin};
1571  else
1572  return {this, static_cast<size_type>(p - ranges.begin()), global_index};
1573 }
1574 
1575 
1576 
1579 {
1580  compress();
1581  return IndexSet::ElementIterator(this);
1582 }
1583 
1584 
1585 
1588 {
1589  compress();
1590  if (ranges.size() > 0)
1591  return IndexSet::IntervalIterator(this, 0);
1592  else
1593  return end_intervals();
1594 }
1595 
1596 
1597 
1600 {
1601  compress();
1602  return IndexSet::IntervalIterator(this);
1603 }
1604 
1605 
1606 
1607 inline void
1609 {
1610  // reset so that there are no indices in the set any more; however,
1611  // as documented, the index set retains its size
1612  ranges.clear();
1613  is_compressed = true;
1615 }
1616 
1617 
1618 
1619 inline void
1621 {
1622  Assert(ranges.empty(),
1623  ExcMessage("This function can only be called if the current "
1624  "object does not yet contain any elements."));
1625  index_space_size = sz;
1626  is_compressed = true;
1627 }
1628 
1629 
1630 
1631 inline IndexSet::size_type
1633 {
1634  return index_space_size;
1635 }
1636 
1637 
1638 
1639 inline void
1641 {
1642  if (is_compressed == true)
1643  return;
1644 
1645  do_compress();
1646 }
1647 
1648 
1649 
1650 inline void
1652 {
1654 
1655  const Range new_range(index, index + 1);
1656  if (ranges.size() == 0 || index > ranges.back().end)
1657  ranges.push_back(new_range);
1658  else if (index == ranges.back().end)
1659  ranges.back().end++;
1660  else
1661  ranges.insert(Utilities::lower_bound(ranges.begin(),
1662  ranges.end(),
1663  new_range),
1664  new_range);
1665  is_compressed = false;
1666 }
1667 
1668 
1669 
1670 inline void
1672 {
1673  Assert((begin < index_space_size) ||
1674  ((begin == index_space_size) && (end == index_space_size)),
1675  ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1676  Assert(end <= index_space_size,
1677  ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1678  AssertIndexRange(begin, end + 1);
1679 
1680  if (begin != end)
1681  {
1682  const Range new_range(begin, end);
1683 
1684  // the new index might be larger than the last index present in the
1685  // ranges. Then we can skip the binary search
1686  if (ranges.size() == 0 || begin > ranges.back().end)
1687  ranges.push_back(new_range);
1688  else
1689  ranges.insert(Utilities::lower_bound(ranges.begin(),
1690  ranges.end(),
1691  new_range),
1692  new_range);
1693  is_compressed = false;
1694  }
1695 }
1696 
1697 
1698 
1699 template <typename ForwardIterator>
1700 inline void
1701 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1702 {
1703  if (begin == end)
1704  return;
1705 
1706  // identify ranges in the given iterator range by checking whether some
1707  // indices happen to be consecutive. to avoid quadratic complexity when
1708  // calling add_range many times (as add_range() going into the middle of an
1709  // already existing range must shift entries around), we first collect a
1710  // vector of ranges.
1711  std::vector<std::pair<size_type, size_type>> tmp_ranges;
1712  bool ranges_are_sorted = true;
1713  for (ForwardIterator p = begin; p != end;)
1714  {
1715  const size_type begin_index = *p;
1716  size_type end_index = begin_index + 1;
1717  ForwardIterator q = p;
1718  ++q;
1719  while ((q != end) && (*q == end_index))
1720  {
1721  ++end_index;
1722  ++q;
1723  }
1724 
1725  tmp_ranges.emplace_back(begin_index, end_index);
1726  p = q;
1727 
1728  // if the starting index of the next go-around of the for loop is less
1729  // than the end index of the one just identified, then we will have at
1730  // least one pair of ranges that are not sorted, and consequently the
1731  // whole collection of ranges is not sorted.
1732  if (p != end && *p < end_index)
1733  ranges_are_sorted = false;
1734  }
1735 
1736  if (!ranges_are_sorted)
1737  std::sort(tmp_ranges.begin(), tmp_ranges.end());
1738 
1739  // if we have many ranges, we first construct a temporary index set (where
1740  // we add ranges in a consecutive way, so fast), otherwise, we work with
1741  // add_range(). the number 9 is chosen heuristically given the fact that
1742  // there are typically up to 8 independent ranges when adding the degrees of
1743  // freedom on a 3D cell or 9 when adding degrees of freedom of faces. if
1744  // doing cell-by-cell additions, we want to avoid repeated calls to
1745  // IndexSet::compress() which gets called upon merging two index sets, so we
1746  // want to be in the other branch then.
1747  if (tmp_ranges.size() > 9)
1748  {
1749  IndexSet tmp_set(size());
1750  tmp_set.ranges.reserve(tmp_ranges.size());
1751  for (const auto &i : tmp_ranges)
1752  tmp_set.add_range(i.first, i.second);
1753  this->add_indices(tmp_set);
1754  }
1755  else
1756  for (const auto &i : tmp_ranges)
1757  add_range(i.first, i.second);
1758 }
1759 
1760 
1761 
1762 inline bool
1764 {
1765  if (ranges.empty() == false)
1766  {
1767  compress();
1768 
1769  // fast check whether the index is in the largest range
1771  if (index >= ranges[largest_range].begin &&
1772  index < ranges[largest_range].end)
1773  return true;
1774 
1775  // get the element after which we would have to insert a range that
1776  // consists of all elements from this element to the end of the index
1777  // range plus one. after this call we know that if p!=end() then
1778  // p->begin<=index unless there is no such range at all
1779  //
1780  // if the searched for element is an element of this range, then we're
1781  // done. otherwise, the element can't be in one of the following ranges
1782  // because otherwise p would be a different iterator
1783  //
1784  // since we already know the position relative to the largest range (we
1785  // called compress!), we can perform the binary search on ranges with
1786  // lower/higher number compared to the largest range
1787  std::vector<Range>::const_iterator p = std::upper_bound(
1788  ranges.begin() +
1789  (index < ranges[largest_range].begin ? 0 : largest_range + 1),
1790  index < ranges[largest_range].begin ? ranges.begin() + largest_range :
1791  ranges.end(),
1792  Range(index, size() + 1));
1793 
1794  if (p == ranges.begin())
1795  return ((index >= p->begin) && (index < p->end));
1796 
1797  Assert((p == ranges.end()) || (p->begin > index), ExcInternalError());
1798 
1799  // now move to that previous range
1800  --p;
1801  Assert(p->begin <= index, ExcInternalError());
1802 
1803  return (p->end > index);
1804  }
1805 
1806  // didn't find this index, so it's not in the set
1807  return false;
1808 }
1809 
1810 
1811 
1812 inline bool
1814 {
1815  compress();
1816  return (ranges.size() <= 1);
1817 }
1818 
1819 
1820 
1821 inline bool
1823 {
1824  return ranges.empty();
1825 }
1826 
1827 
1828 
1829 inline IndexSet::size_type
1831 {
1832  // make sure we have non-overlapping ranges
1833  compress();
1834 
1835  size_type v = 0;
1836  if (!ranges.empty())
1837  {
1838  Range &r = ranges.back();
1839  v = r.nth_index_in_set + r.end - r.begin;
1840  }
1841 
1842 #ifdef DEBUG
1843  size_type s = 0;
1844  for (const auto &range : ranges)
1845  s += (range.end - range.begin);
1846  Assert(s == v, ExcInternalError());
1847 #endif
1848 
1849  return v;
1850 }
1851 
1852 
1853 
1854 inline unsigned int
1856 {
1857  compress();
1858  return ranges.size();
1859 }
1860 
1861 
1862 
1863 inline IndexSet::size_type
1865 {
1866  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1867 
1868  compress();
1869  const std::vector<Range>::const_iterator main_range =
1870  ranges.begin() + largest_range;
1871 
1872  return main_range->nth_index_in_set;
1873 }
1874 
1875 
1876 
1877 inline IndexSet::size_type
1879 {
1881 
1882  compress();
1883 
1884  // first check whether the index is in the largest range
1886  std::vector<Range>::const_iterator main_range =
1887  ranges.begin() + largest_range;
1888  if (n >= main_range->nth_index_in_set &&
1889  n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1890  return main_range->begin + (n - main_range->nth_index_in_set);
1891 
1892  // find out which chunk the local index n belongs to by using a binary
1893  // search. the comparator is based on the end of the ranges. Use the
1894  // position relative to main_range to subdivide the ranges
1895  Range r(n, n + 1);
1896  r.nth_index_in_set = n;
1897  std::vector<Range>::const_iterator range_begin, range_end;
1898  if (n < main_range->nth_index_in_set)
1899  {
1900  range_begin = ranges.begin();
1901  range_end = main_range;
1902  }
1903  else
1904  {
1905  range_begin = main_range + 1;
1906  range_end = ranges.end();
1907  }
1908 
1909  const std::vector<Range>::const_iterator p =
1910  Utilities::lower_bound(range_begin, range_end, r, Range::nth_index_compare);
1911 
1912  Assert(p != ranges.end(), ExcInternalError());
1913  return p->begin + (n - p->nth_index_in_set);
1914 }
1915 
1916 
1917 
1918 inline IndexSet::size_type
1920 {
1921  // to make this call thread-safe, compress() must not be called through this
1922  // function
1923  Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1924  AssertIndexRange(n, size());
1925 
1926  // return immediately if the index set is empty
1927  if (is_empty())
1929 
1930  // check whether the index is in the largest range. use the result to
1931  // perform a one-sided binary search afterward
1933  std::vector<Range>::const_iterator main_range =
1934  ranges.begin() + largest_range;
1935  if (n >= main_range->begin && n < main_range->end)
1936  return (n - main_range->begin) + main_range->nth_index_in_set;
1937 
1938  Range r(n, n);
1939  std::vector<Range>::const_iterator range_begin, range_end;
1940  if (n < main_range->begin)
1941  {
1942  range_begin = ranges.begin();
1943  range_end = main_range;
1944  }
1945  else
1946  {
1947  range_begin = main_range + 1;
1948  range_end = ranges.end();
1949  }
1950 
1951  std::vector<Range>::const_iterator p =
1952  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1953 
1954  // if n is not in this set
1955  if (p == range_end || p->end == n || p->begin > n)
1957 
1958  Assert(p != ranges.end(), ExcInternalError());
1959  Assert(p->begin <= n, ExcInternalError());
1960  Assert(n < p->end, ExcInternalError());
1961  return (n - p->begin) + p->nth_index_in_set;
1962 }
1963 
1964 
1965 
1966 inline bool
1968 {
1969  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1970 
1971  compress();
1972  is.compress();
1973 
1974  return ranges == is.ranges;
1975 }
1976 
1977 
1978 
1979 inline bool
1981 {
1982  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1983 
1984  compress();
1985  is.compress();
1986 
1987  return ranges != is.ranges;
1988 }
1989 
1990 
1991 
1992 template <typename Vector>
1993 void
1995 {
1996  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
1997 
1998  compress();
1999  // first fill all elements of the vector with zeroes.
2000  std::fill(vector.begin(), vector.end(), 0);
2001 
2002  // then write ones into the elements whose indices are contained in the
2003  // index set
2004  for (const auto &range : ranges)
2005  for (size_type i = range.begin; i < range.end; ++i)
2006  vector[i] = 1;
2007 }
2008 
2009 
2010 
2011 template <class StreamType>
2012 inline void
2013 IndexSet::print(StreamType &out) const
2014 {
2015  compress();
2016  out << "{";
2017  std::vector<Range>::const_iterator p;
2018  for (p = ranges.begin(); p != ranges.end(); ++p)
2019  {
2020  if (p->end - p->begin == 1)
2021  out << p->begin;
2022  else
2023  out << "[" << p->begin << "," << p->end - 1 << "]";
2024 
2025  if (p != --ranges.end())
2026  out << ", ";
2027  }
2028  out << "}" << std::endl;
2029 }
2030 
2031 
2032 
2033 template <class Archive>
2034 inline void
2035 IndexSet::Range::serialize(Archive &ar, const unsigned int)
2036 {
2037  ar &begin &end &nth_index_in_set;
2038 }
2039 
2040 
2041 
2042 template <class Archive>
2043 inline void
2044 IndexSet::serialize(Archive &ar, const unsigned int)
2045 {
2047 }
2048 
2050 
2051 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1075
friend bool operator==(const Range &range_1, const Range &range_2)
Definition: index_set.h:922
void fill_binary_vector(VectorType &vector) const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
ElementIterator end() const
Definition: index_set.h:1578
void add_index(const size_type index)
Definition: index_set.h:1651
const IntervalAccessor * operator->() const
Definition: index_set.h:1205
const IntervalAccessor & operator*() const
Definition: index_set.h:1197
const IndexSet * index_set
Definition: index_set.h:806
IndexSet & operator=(const IndexSet &)=default
void block_read(std::istream &in)
Definition: index_set.cc:487
size_type largest_range_starting_index() const
Definition: index_set.h:1864
#define AssertIndexRange(index, range)
Definition: exceptions.h:1636
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1403
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1701
types::global_dof_index size_type
Definition: index_set.h:83
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition: index_set.h:1265
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1232
std::forward_iterator_tag iterator_category
Definition: index_set.h:699
std::size_t memory_consumption() const
Definition: index_set.cc:732
size_type end
Definition: index_set.h:877
size_type n_elements() const
Definition: index_set.h:1055
size_type index_space_size
Definition: index_set.h:966
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:527
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:666
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:596
size_type size() const
Definition: index_set.h:1632
void print(StreamType &out) const
Definition: index_set.h:2013
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2035
void clear()
Definition: index_set.h:1608
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2044
static ::ExceptionBase & ExcMessage(std::string arg1)
void block_write(std::ostream &out) const
Definition: index_set.cc:473
std::vector< Range > ranges
Definition: index_set.h:950
static ::ExceptionBase & ExcIndexNotPresent(size_type arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
void subtract_set(const IndexSet &other)
Definition: index_set.cc:258
#define Assert(cond, exc)
Definition: exceptions.h:1411
friend class IntervalIterator
Definition: index_set.h:608
std::forward_iterator_tag iterator_category
Definition: index_set.h:790
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
bool operator<(const ElementIterator &) const
Definition: index_set.h:1390
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1919
bool operator==(const IndexSet &is) const
Definition: index_set.h:1967
signed int value_type
Definition: index_set.h:102
IntervalAccessor accessor
Definition: index_set.h:709
void do_compress() const
Definition: index_set.cc:98
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1105
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:909
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1214
size_type last() const
Definition: index_set.h:1095
const IndexSet * index_set
Definition: index_set.h:600
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1118
unsigned int n_intervals() const
Definition: index_set.h:1855
IndexSet tensor_product(const IndexSet &other) const
Definition: index_set.cc:330
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1241
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1980
size_type pop_back()
Definition: index_set.cc:342
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:507
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1671
IndexSet operator &(const IndexSet &is) const
void set_size(const size_type size)
Definition: index_set.h:1620
ElementIterator & operator++()
Definition: index_set.h:1361
unsigned int global_dof_index
Definition: types.h:76
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:1024
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
Definition: index_set.cc:238
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:211
void compress() const
Definition: index_set.h:1640
ElementIterator end() const
Definition: index_set.h:1081
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:915
static std::size_t memory_consumption()
Definition: index_set.h:928
IntervalIterator begin_intervals() const
Definition: index_set.h:1587
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1011
size_type begin
Definition: index_set.h:876
bool is_contiguous() const
Definition: index_set.h:1813
Threads::Mutex compress_mutex
Definition: index_set.h:983
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
size_type pop_front()
Definition: index_set.cc:360
static const char N
IntervalIterator end_intervals() const
Definition: index_set.h:1599
bool is_compressed
Definition: index_set.h:960
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1130
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1381
ElementIterator begin() const
Definition: index_set.h:1072
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
void write(std::ostream &out) const
Definition: index_set.cc:436
bool is_element(const size_type index) const
Definition: index_set.h:1763
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
bool operator==(const ElementIterator &) const
Definition: index_set.h:1321
ElementIterator at(const size_type global_index) const
Definition: index_set.h:1528
const types::global_dof_index invalid_dof_index
Definition: types.h:211
std::ptrdiff_t difference_type
Definition: index_set.h:792
std::ptrdiff_t difference_type
Definition: index_set.h:701
ElementIterator begin() const
Definition: index_set.h:1516
void read(std::istream &in)
Definition: index_set.cc:451
size_type n_elements() const
Definition: index_set.h:1830
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1878
bool is_empty() const
Definition: index_set.h:1822
size_type largest_range
Definition: index_set.h:977
IntervalIterator & operator++()
Definition: index_set.h:1178
size_type nth_index_in_set
Definition: index_set.h:879
static ::ExceptionBase & ExcInternalError()
size_type operator*() const
Definition: index_set.h:1308
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1223