Reference documentation for deal.II version Git 6d836ad036 2021-06-16 20:57:45 +0200
\(\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 - 2021 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 
511  template <class Archive>
512  void
513  serialize(Archive &ar, const unsigned int version);
514 
515 
527  {
528  public:
533  IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
534 
538  explicit IntervalAccessor(const IndexSet *idxset);
539 
543  size_type
544  n_elements() const;
545 
549  bool
550  is_valid() const;
551 
556  begin() const;
557 
563  end() const;
564 
568  size_type
569  last() const;
570 
571  private:
575  IntervalAccessor(const IntervalAccessor &other);
580  operator=(const IntervalAccessor &other);
581 
585  bool
586  operator==(const IntervalAccessor &other) const;
590  bool
591  operator<(const IntervalAccessor &other) const;
596  void
597  advance();
602 
608 
609  friend class IntervalIterator;
610  };
611 
617  {
618  public:
623  IntervalIterator(const IndexSet *idxset, const size_type range_idx);
624 
628  explicit IntervalIterator(const IndexSet *idxset);
629 
634 
638  IntervalIterator(const IntervalIterator &other) = default;
639 
644  operator=(const IntervalIterator &other) = default;
645 
650  operator++();
651 
656  operator++(int);
657 
661  const IntervalAccessor &operator*() const;
662 
666  const IntervalAccessor *operator->() const;
667 
671  bool
672  operator==(const IntervalIterator &) const;
673 
677  bool
678  operator!=(const IntervalIterator &) const;
679 
683  bool
684  operator<(const IntervalIterator &) const;
685 
692  int
693  operator-(const IntervalIterator &p) const;
694 
700  using iterator_category = std::forward_iterator_tag;
702  using difference_type = std::ptrdiff_t;
705 
706  private:
711  };
712 
718  {
719  public:
724  ElementIterator(const IndexSet *idxset,
725  const size_type range_idx,
726  const size_type index);
727 
731  explicit ElementIterator(const IndexSet *idxset);
732 
737  size_type operator*() const;
738 
742  bool
743  is_valid() const;
744 
749  operator++();
750 
755  operator++(int);
756 
760  bool
761  operator==(const ElementIterator &) const;
762 
766  bool
767  operator!=(const ElementIterator &) const;
768 
772  bool
773  operator<(const ElementIterator &) const;
774 
783  std::ptrdiff_t
784  operator-(const ElementIterator &p) const;
785 
791  using iterator_category = std::forward_iterator_tag;
793  using difference_type = std::ptrdiff_t;
794  using pointer = size_type *;
795  using reference = size_type &;
796 
797  private:
801  void
802  advance();
803 
816  };
817 
823  begin() const;
824 
840  at(const size_type global_index) const;
841 
847  end() const;
848 
853  begin_intervals() const;
854 
860  end_intervals() const;
861 
866 private:
875  struct Range
876  {
879 
881 
890  Range();
891 
899  Range(const size_type i1, const size_type i2);
900 
901  friend inline bool
902  operator<(const Range &range_1, const Range &range_2)
903  {
904  return (
905  (range_1.begin < range_2.begin) ||
906  ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
907  }
908 
909  static bool
911  {
912  return x.end < y.end;
913  }
914 
915  static bool
917  {
918  return (x.nth_index_in_set + (x.end - x.begin) <
919  y.nth_index_in_set + (y.end - y.begin));
920  }
921 
922  friend inline bool
923  operator==(const Range &range_1, const Range &range_2)
924  {
925  return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
926  }
927 
928  static std::size_t
930  {
931  return sizeof(Range);
932  }
933 
939  template <class Archive>
940  void
941  serialize(Archive &ar, const unsigned int version);
942  };
943 
952  mutable std::vector<Range> ranges;
953 
962  mutable bool is_compressed;
963 
969 
980 
986 
990  void
991  do_compress() const;
992 };
993 
994 
1012 inline IndexSet
1014 {
1015  IndexSet is(N);
1016  is.add_range(0, N);
1017  is.compress();
1018  return is;
1019 }
1020 
1021 /* ------------------ inline functions ------------------ */
1022 
1023 
1024 /* IntervalAccessor */
1025 
1027  const IndexSet * idxset,
1029  : index_set(idxset)
1030  , range_idx(range_idx)
1031 {
1032  Assert(range_idx < idxset->n_intervals(),
1033  ExcInternalError("Invalid range index"));
1034 }
1035 
1036 
1037 
1039  : index_set(idxset)
1041 {}
1042 
1043 
1044 
1046  const IndexSet::IntervalAccessor &other)
1047  : index_set(other.index_set)
1048  , range_idx(other.range_idx)
1049 {
1051  ExcMessage("invalid iterator"));
1052 }
1053 
1054 
1055 
1056 inline IndexSet::size_type
1058 {
1059  Assert(is_valid(), ExcMessage("invalid iterator"));
1060  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1061 }
1062 
1063 
1064 
1065 inline bool
1067 {
1068  return index_set != nullptr && range_idx < index_set->n_intervals();
1069 }
1070 
1071 
1072 
1075 {
1076  Assert(is_valid(), ExcMessage("invalid iterator"));
1077  return {index_set, range_idx, index_set->ranges[range_idx].begin};
1078 }
1079 
1080 
1081 
1084 {
1085  Assert(is_valid(), ExcMessage("invalid iterator"));
1086 
1087  // point to first index in next interval unless we are the last interval.
1088  if (range_idx < index_set->ranges.size() - 1)
1089  return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1090  else
1091  return index_set->end();
1092 }
1093 
1094 
1095 
1096 inline IndexSet::size_type
1098 {
1099  Assert(is_valid(), ExcMessage("invalid iterator"));
1100 
1101  return index_set->ranges[range_idx].end - 1;
1102 }
1103 
1104 
1105 
1108 {
1109  index_set = other.index_set;
1110  range_idx = other.range_idx;
1112  ExcMessage("invalid iterator"));
1113  return *this;
1114 }
1115 
1116 
1117 
1118 inline bool
1121 {
1122  Assert(index_set == other.index_set,
1123  ExcMessage(
1124  "Can not compare accessors pointing to different IndexSets"));
1125  return range_idx == other.range_idx;
1126 }
1127 
1128 
1129 
1130 inline bool
1133 {
1134  Assert(index_set == other.index_set,
1135  ExcMessage(
1136  "Can not compare accessors pointing to different IndexSets"));
1137  return range_idx < other.range_idx;
1138 }
1139 
1140 
1141 
1142 inline void
1144 {
1145  Assert(
1146  is_valid(),
1147  ExcMessage(
1148  "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1149  ++range_idx;
1150 
1151  // set ourselves to invalid if we walk off the end
1152  if (range_idx >= index_set->ranges.size())
1154 }
1155 
1156 
1157 /* IntervalIterator */
1158 
1160  const IndexSet * idxset,
1162  : accessor(idxset, range_idx)
1163 {}
1164 
1165 
1166 
1168  : accessor(nullptr)
1169 {}
1170 
1171 
1172 
1174  : accessor(idxset)
1175 {}
1176 
1177 
1178 
1181 {
1182  accessor.advance();
1183  return *this;
1184 }
1185 
1186 
1187 
1190 {
1191  const IndexSet::IntervalIterator iter = *this;
1192  accessor.advance();
1193  return iter;
1194 }
1195 
1196 
1197 
1199  operator*() const
1200 {
1201  return accessor;
1202 }
1203 
1204 
1205 
1208 {
1209  return &accessor;
1210 }
1211 
1212 
1213 
1214 inline bool
1217 {
1218  return accessor == other.accessor;
1219 }
1220 
1221 
1222 
1223 inline bool
1226 {
1227  return !(*this == other);
1228 }
1229 
1230 
1231 
1232 inline bool
1235 {
1236  return accessor < other.accessor;
1237 }
1238 
1239 
1240 
1241 inline int
1244 {
1246  ExcMessage(
1247  "Can not compare iterators belonging to different IndexSets"));
1248 
1250  accessor.index_set->ranges.size() :
1252  const size_type rhs =
1254  accessor.index_set->ranges.size() :
1255  other.accessor.range_idx;
1256 
1257  if (lhs > rhs)
1258  return static_cast<int>(lhs - rhs);
1259  else
1260  return -static_cast<int>(rhs - lhs);
1261 }
1262 
1263 
1264 
1265 /* ElementIterator */
1266 
1268  const IndexSet * idxset,
1269  const IndexSet::size_type range_idx,
1270  const IndexSet::size_type index)
1271  : index_set(idxset)
1272  , range_idx(range_idx)
1273  , idx(index)
1274 {
1275  Assert(range_idx < index_set->ranges.size(),
1276  ExcMessage(
1277  "Invalid range index for IndexSet::ElementIterator constructor."));
1278  Assert(
1279  idx >= index_set->ranges[range_idx].begin &&
1280  idx < index_set->ranges[range_idx].end,
1282  "Invalid index argument for IndexSet::ElementIterator constructor."));
1283 }
1284 
1285 
1286 
1288  : index_set(idxset)
1291 {}
1292 
1293 
1294 
1295 inline bool
1297 {
1300  (range_idx < index_set->ranges.size() &&
1301  idx < index_set->ranges[range_idx].end),
1302  ExcInternalError("Invalid ElementIterator state."));
1303 
1304  return (range_idx < index_set->ranges.size() &&
1305  idx < index_set->ranges[range_idx].end);
1306 }
1307 
1308 
1309 
1311 {
1312  Assert(
1313  is_valid(),
1314  ExcMessage(
1315  "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1316  return idx;
1317 }
1318 
1319 
1320 
1321 inline bool
1324 {
1325  Assert(index_set == other.index_set,
1326  ExcMessage(
1327  "Can not compare iterators belonging to different IndexSets"));
1328  return range_idx == other.range_idx && idx == other.idx;
1329 }
1330 
1331 
1332 
1333 inline void
1335 {
1336  Assert(
1337  is_valid(),
1338  ExcMessage(
1339  "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1340  if (idx < index_set->ranges[range_idx].end)
1341  ++idx;
1342  // end of this range?
1343  if (idx == index_set->ranges[range_idx].end)
1344  {
1345  // point to first element in next interval if possible
1346  if (range_idx < index_set->ranges.size() - 1)
1347  {
1348  ++range_idx;
1349  idx = index_set->ranges[range_idx].begin;
1350  }
1351  else
1352  {
1353  // we just fell off the end, set to invalid:
1356  }
1357  }
1358 }
1359 
1360 
1361 
1364 {
1365  advance();
1366  return *this;
1367 }
1368 
1369 
1370 
1373 {
1374  const IndexSet::ElementIterator it = *this;
1375  advance();
1376  return it;
1377 }
1378 
1379 
1380 
1381 inline bool
1384 {
1385  return !(*this == other);
1386 }
1387 
1388 
1389 
1390 inline bool
1393 {
1394  Assert(index_set == other.index_set,
1395  ExcMessage(
1396  "Can not compare iterators belonging to different IndexSets"));
1397  return range_idx < other.range_idx ||
1398  (range_idx == other.range_idx && idx < other.idx);
1399 }
1400 
1401 
1402 
1403 inline std::ptrdiff_t
1406 {
1407  Assert(index_set == other.index_set,
1408  ExcMessage(
1409  "Can not compare iterators belonging to different IndexSets"));
1410  if (*this == other)
1411  return 0;
1412  if (!(*this < other))
1413  return -(other - *this);
1414 
1415  // only other can be equal to end() because of the checks above.
1417 
1418  // Note: we now compute how far advance *this in "*this < other" to get other,
1419  // so we need to return -c at the end.
1420 
1421  // first finish the current range:
1422  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1423 
1424  // now walk in steps of ranges (need to start one behind our current one):
1425  for (size_type range = range_idx + 1;
1426  range < index_set->ranges.size() && range <= other.range_idx;
1427  ++range)
1428  c += index_set->ranges[range].end - index_set->ranges[range].begin;
1429 
1430  Assert(
1431  other.range_idx < index_set->ranges.size() ||
1433  ExcMessage(
1434  "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1435 
1436  // We might have walked too far because we went until the end of
1437  // other.range_idx, so walk backwards to other.idx:
1439  c -= index_set->ranges[other.range_idx].end - other.idx;
1440 
1441  return -c;
1442 }
1443 
1444 
1445 /* Range */
1446 
1451 {}
1452 
1453 
1454 
1455 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1456  : begin(i1)
1457  , end(i2)
1459 {}
1460 
1461 
1462 
1463 /* IndexSet itself */
1464 
1466  : is_compressed(true)
1467  , index_space_size(0)
1469 {}
1470 
1471 
1472 
1474  : is_compressed(true)
1475  , index_space_size(size)
1477 {}
1478 
1479 
1480 
1481 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1482  : ranges(std::move(is.ranges))
1483  , is_compressed(is.is_compressed)
1484  , index_space_size(is.index_space_size)
1485  , largest_range(is.largest_range)
1486 {
1487  is.ranges.clear();
1488  is.is_compressed = true;
1489  is.index_space_size = 0;
1490  is.largest_range = numbers::invalid_unsigned_int;
1491 
1492  compress();
1493 }
1494 
1495 
1496 
1497 inline IndexSet &
1499 {
1500  ranges = std::move(is.ranges);
1501  is_compressed = is.is_compressed;
1502  index_space_size = is.index_space_size;
1503  largest_range = is.largest_range;
1504 
1505  is.ranges.clear();
1506  is.is_compressed = true;
1507  is.index_space_size = 0;
1508  is.largest_range = numbers::invalid_unsigned_int;
1509 
1510  compress();
1511 
1512  return *this;
1513 }
1514 
1515 
1516 
1519 {
1520  compress();
1521  if (ranges.size() > 0)
1522  return {this, 0, ranges[0].begin};
1523  else
1524  return end();
1525 }
1526 
1527 
1528 
1531 {
1532  compress();
1533  AssertIndexRange(global_index, size());
1534 
1535  if (ranges.empty())
1536  return end();
1537 
1538  std::vector<Range>::const_iterator main_range =
1539  ranges.begin() + largest_range;
1540 
1541  Range r(global_index, global_index + 1);
1542  // This optimization makes the bounds for lower_bound smaller by checking
1543  // the largest range first.
1544  std::vector<Range>::const_iterator range_begin, range_end;
1545  if (global_index < main_range->begin)
1546  {
1547  range_begin = ranges.begin();
1548  range_end = main_range;
1549  }
1550  else
1551  {
1552  range_begin = main_range;
1553  range_end = ranges.end();
1554  }
1555 
1556  // This will give us the first range p=[a,b[ with b>=global_index using
1557  // a binary search
1558  const std::vector<Range>::const_iterator p =
1559  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1560 
1561  // We couldn't find a range, which means we have no range that contains
1562  // global_index and also no range behind it, meaning we need to return end().
1563  if (p == ranges.end())
1564  return end();
1565 
1566  // Finally, we can have two cases: Either global_index is not in [a,b[,
1567  // which means we need to return an iterator to a because global_index, ...,
1568  // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
1569  // [a,b[ and we will return an iterator pointing directly at global_index
1570  // (else branch).
1571  if (global_index < p->begin)
1572  return {this, static_cast<size_type>(p - ranges.begin()), p->begin};
1573  else
1574  return {this, static_cast<size_type>(p - ranges.begin()), global_index};
1575 }
1576 
1577 
1578 
1581 {
1582  compress();
1583  return IndexSet::ElementIterator(this);
1584 }
1585 
1586 
1587 
1590 {
1591  compress();
1592  if (ranges.size() > 0)
1593  return IndexSet::IntervalIterator(this, 0);
1594  else
1595  return end_intervals();
1596 }
1597 
1598 
1599 
1602 {
1603  compress();
1604  return IndexSet::IntervalIterator(this);
1605 }
1606 
1607 
1608 
1609 inline void
1611 {
1612  // reset so that there are no indices in the set any more; however,
1613  // as documented, the index set retains its size
1614  ranges.clear();
1615  is_compressed = true;
1617 }
1618 
1619 
1620 
1621 inline void
1623 {
1624  Assert(ranges.empty(),
1625  ExcMessage("This function can only be called if the current "
1626  "object does not yet contain any elements."));
1627  index_space_size = sz;
1628  is_compressed = true;
1629 }
1630 
1631 
1632 
1633 inline IndexSet::size_type
1635 {
1636  return index_space_size;
1637 }
1638 
1639 
1640 
1641 inline void
1643 {
1644  if (is_compressed == true)
1645  return;
1646 
1647  do_compress();
1648 }
1649 
1650 
1651 
1652 inline void
1654 {
1656 
1657  const Range new_range(index, index + 1);
1658  if (ranges.size() == 0 || index > ranges.back().end)
1659  ranges.push_back(new_range);
1660  else if (index == ranges.back().end)
1661  ranges.back().end++;
1662  else
1663  ranges.insert(Utilities::lower_bound(ranges.begin(),
1664  ranges.end(),
1665  new_range),
1666  new_range);
1667  is_compressed = false;
1668 }
1669 
1670 
1671 
1672 inline void
1674 {
1675  Assert((begin < index_space_size) ||
1676  ((begin == index_space_size) && (end == index_space_size)),
1677  ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1678  Assert(end <= index_space_size,
1679  ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1680  AssertIndexRange(begin, end + 1);
1681 
1682  if (begin != end)
1683  {
1684  const Range new_range(begin, end);
1685 
1686  // the new index might be larger than the last index present in the
1687  // ranges. Then we can skip the binary search
1688  if (ranges.size() == 0 || begin > ranges.back().end)
1689  ranges.push_back(new_range);
1690  else
1691  ranges.insert(Utilities::lower_bound(ranges.begin(),
1692  ranges.end(),
1693  new_range),
1694  new_range);
1695  is_compressed = false;
1696  }
1697 }
1698 
1699 
1700 
1701 template <typename ForwardIterator>
1702 inline void
1703 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1704 {
1705  if (begin == end)
1706  return;
1707 
1708  // identify ranges in the given iterator range by checking whether some
1709  // indices happen to be consecutive. to avoid quadratic complexity when
1710  // calling add_range many times (as add_range() going into the middle of an
1711  // already existing range must shift entries around), we first collect a
1712  // vector of ranges.
1713  std::vector<std::pair<size_type, size_type>> tmp_ranges;
1714  bool ranges_are_sorted = true;
1715  for (ForwardIterator p = begin; p != end;)
1716  {
1717  const size_type begin_index = *p;
1718  size_type end_index = begin_index + 1;
1719  ForwardIterator q = p;
1720  ++q;
1721  while ((q != end) && (*q == end_index))
1722  {
1723  ++end_index;
1724  ++q;
1725  }
1726 
1727  tmp_ranges.emplace_back(begin_index, end_index);
1728  p = q;
1729 
1730  // if the starting index of the next go-around of the for loop is less
1731  // than the end index of the one just identified, then we will have at
1732  // least one pair of ranges that are not sorted, and consequently the
1733  // whole collection of ranges is not sorted.
1734  if (p != end && *p < end_index)
1735  ranges_are_sorted = false;
1736  }
1737 
1738  if (!ranges_are_sorted)
1739  std::sort(tmp_ranges.begin(), tmp_ranges.end());
1740 
1741  // if we have many ranges, we first construct a temporary index set (where
1742  // we add ranges in a consecutive way, so fast), otherwise, we work with
1743  // add_range(). the number 9 is chosen heuristically given the fact that
1744  // there are typically up to 8 independent ranges when adding the degrees of
1745  // freedom on a 3D cell or 9 when adding degrees of freedom of faces. if
1746  // doing cell-by-cell additions, we want to avoid repeated calls to
1747  // IndexSet::compress() which gets called upon merging two index sets, so we
1748  // want to be in the other branch then.
1749  if (tmp_ranges.size() > 9)
1750  {
1751  IndexSet tmp_set(size());
1752  tmp_set.ranges.reserve(tmp_ranges.size());
1753  for (const auto &i : tmp_ranges)
1754  tmp_set.add_range(i.first, i.second);
1755  this->add_indices(tmp_set);
1756  }
1757  else
1758  for (const auto &i : tmp_ranges)
1759  add_range(i.first, i.second);
1760 }
1761 
1762 
1763 
1764 inline bool
1766 {
1767  if (ranges.empty() == false)
1768  {
1769  compress();
1770 
1771  // fast check whether the index is in the largest range
1773  if (index >= ranges[largest_range].begin &&
1774  index < ranges[largest_range].end)
1775  return true;
1776 
1777  // get the element after which we would have to insert a range that
1778  // consists of all elements from this element to the end of the index
1779  // range plus one. after this call we know that if p!=end() then
1780  // p->begin<=index unless there is no such range at all
1781  //
1782  // if the searched for element is an element of this range, then we're
1783  // done. otherwise, the element can't be in one of the following ranges
1784  // because otherwise p would be a different iterator
1785  //
1786  // since we already know the position relative to the largest range (we
1787  // called compress!), we can perform the binary search on ranges with
1788  // lower/higher number compared to the largest range
1789  std::vector<Range>::const_iterator p = std::upper_bound(
1790  ranges.begin() +
1791  (index < ranges[largest_range].begin ? 0 : largest_range + 1),
1792  index < ranges[largest_range].begin ? ranges.begin() + largest_range :
1793  ranges.end(),
1794  Range(index, size() + 1));
1795 
1796  if (p == ranges.begin())
1797  return ((index >= p->begin) && (index < p->end));
1798 
1799  Assert((p == ranges.end()) || (p->begin > index), ExcInternalError());
1800 
1801  // now move to that previous range
1802  --p;
1803  Assert(p->begin <= index, ExcInternalError());
1804 
1805  return (p->end > index);
1806  }
1807 
1808  // didn't find this index, so it's not in the set
1809  return false;
1810 }
1811 
1812 
1813 
1814 inline bool
1816 {
1817  compress();
1818  return (ranges.size() <= 1);
1819 }
1820 
1821 
1822 
1823 inline bool
1825 {
1826  return ranges.empty();
1827 }
1828 
1829 
1830 
1831 inline IndexSet::size_type
1833 {
1834  // make sure we have non-overlapping ranges
1835  compress();
1836 
1837  size_type v = 0;
1838  if (!ranges.empty())
1839  {
1840  Range &r = ranges.back();
1841  v = r.nth_index_in_set + r.end - r.begin;
1842  }
1843 
1844 #ifdef DEBUG
1845  size_type s = 0;
1846  for (const auto &range : ranges)
1847  s += (range.end - range.begin);
1848  Assert(s == v, ExcInternalError());
1849 #endif
1850 
1851  return v;
1852 }
1853 
1854 
1855 
1856 inline unsigned int
1858 {
1859  compress();
1860  return ranges.size();
1861 }
1862 
1863 
1864 
1865 inline IndexSet::size_type
1867 {
1868  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1869 
1870  compress();
1871  const std::vector<Range>::const_iterator main_range =
1872  ranges.begin() + largest_range;
1873 
1874  return main_range->nth_index_in_set;
1875 }
1876 
1877 
1878 
1879 inline IndexSet::size_type
1881 {
1883 
1884  compress();
1885 
1886  // first check whether the index is in the largest range
1888  std::vector<Range>::const_iterator main_range =
1889  ranges.begin() + largest_range;
1890  if (n >= main_range->nth_index_in_set &&
1891  n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1892  return main_range->begin + (n - main_range->nth_index_in_set);
1893 
1894  // find out which chunk the local index n belongs to by using a binary
1895  // search. the comparator is based on the end of the ranges. Use the
1896  // position relative to main_range to subdivide the ranges
1897  Range r(n, n + 1);
1898  r.nth_index_in_set = n;
1899  std::vector<Range>::const_iterator range_begin, range_end;
1900  if (n < main_range->nth_index_in_set)
1901  {
1902  range_begin = ranges.begin();
1903  range_end = main_range;
1904  }
1905  else
1906  {
1907  range_begin = main_range + 1;
1908  range_end = ranges.end();
1909  }
1910 
1911  const std::vector<Range>::const_iterator p =
1912  Utilities::lower_bound(range_begin, range_end, r, Range::nth_index_compare);
1913 
1914  Assert(p != ranges.end(), ExcInternalError());
1915  return p->begin + (n - p->nth_index_in_set);
1916 }
1917 
1918 
1919 
1920 inline IndexSet::size_type
1922 {
1923  // to make this call thread-safe, compress() must not be called through this
1924  // function
1925  Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1926  AssertIndexRange(n, size());
1927 
1928  // return immediately if the index set is empty
1929  if (is_empty())
1931 
1932  // check whether the index is in the largest range. use the result to
1933  // perform a one-sided binary search afterward
1935  std::vector<Range>::const_iterator main_range =
1936  ranges.begin() + largest_range;
1937  if (n >= main_range->begin && n < main_range->end)
1938  return (n - main_range->begin) + main_range->nth_index_in_set;
1939 
1940  Range r(n, n);
1941  std::vector<Range>::const_iterator range_begin, range_end;
1942  if (n < main_range->begin)
1943  {
1944  range_begin = ranges.begin();
1945  range_end = main_range;
1946  }
1947  else
1948  {
1949  range_begin = main_range + 1;
1950  range_end = ranges.end();
1951  }
1952 
1953  std::vector<Range>::const_iterator p =
1954  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1955 
1956  // if n is not in this set
1957  if (p == range_end || p->end == n || p->begin > n)
1959 
1960  Assert(p != ranges.end(), ExcInternalError());
1961  Assert(p->begin <= n, ExcInternalError());
1962  Assert(n < p->end, ExcInternalError());
1963  return (n - p->begin) + p->nth_index_in_set;
1964 }
1965 
1966 
1967 
1968 inline bool
1970 {
1971  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1972 
1973  compress();
1974  is.compress();
1975 
1976  return ranges == is.ranges;
1977 }
1978 
1979 
1980 
1981 inline bool
1983 {
1984  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1985 
1986  compress();
1987  is.compress();
1988 
1989  return ranges != is.ranges;
1990 }
1991 
1992 
1993 
1994 template <typename Vector>
1995 void
1997 {
1998  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
1999 
2000  compress();
2001  // first fill all elements of the vector with zeroes.
2002  std::fill(vector.begin(), vector.end(), 0);
2003 
2004  // then write ones into the elements whose indices are contained in the
2005  // index set
2006  for (const auto &range : ranges)
2007  for (size_type i = range.begin; i < range.end; ++i)
2008  vector[i] = 1;
2009 }
2010 
2011 
2012 
2013 template <class StreamType>
2014 inline void
2015 IndexSet::print(StreamType &out) const
2016 {
2017  compress();
2018  out << "{";
2019  std::vector<Range>::const_iterator p;
2020  for (p = ranges.begin(); p != ranges.end(); ++p)
2021  {
2022  if (p->end - p->begin == 1)
2023  out << p->begin;
2024  else
2025  out << "[" << p->begin << "," << p->end - 1 << "]";
2026 
2027  if (p != --ranges.end())
2028  out << ", ";
2029  }
2030  out << "}" << std::endl;
2031 }
2032 
2033 
2034 
2035 template <class Archive>
2036 inline void
2037 IndexSet::Range::serialize(Archive &ar, const unsigned int)
2038 {
2039  ar &begin &end &nth_index_in_set;
2040 }
2041 
2042 
2043 
2044 template <class Archive>
2045 inline void
2046 IndexSet::serialize(Archive &ar, const unsigned int)
2047 {
2049 }
2050 
2052 
2053 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1111
friend bool operator==(const Range &range_1, const Range &range_2)
Definition: index_set.h:923
void fill_binary_vector(VectorType &vector) const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
ElementIterator end() const
Definition: index_set.h:1580
void add_index(const size_type index)
Definition: index_set.h:1653
const IntervalAccessor * operator->() const
Definition: index_set.h:1207
const IntervalAccessor & operator*() const
Definition: index_set.h:1199
const IndexSet * index_set
Definition: index_set.h:807
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:1866
#define AssertIndexRange(index, range)
Definition: exceptions.h:1721
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1405
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1703
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:1267
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1234
std::forward_iterator_tag iterator_category
Definition: index_set.h:700
std::size_t memory_consumption() const
Definition: index_set.cc:732
size_type end
Definition: index_set.h:878
size_type n_elements() const
Definition: index_set.h:1057
size_type index_space_size
Definition: index_set.h:968
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:1634
void print(StreamType &out) const
Definition: index_set.h:2015
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2037
void clear()
Definition: index_set.h:1610
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2046
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:952
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:1478
friend class IntervalIterator
Definition: index_set.h:609
std::forward_iterator_tag iterator_category
Definition: index_set.h:791
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
bool operator<(const ElementIterator &) const
Definition: index_set.h:1392
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1921
bool operator==(const IndexSet &is) const
Definition: index_set.h:1969
signed int value_type
Definition: index_set.h:102
IntervalAccessor accessor
Definition: index_set.h:710
void do_compress() const
Definition: index_set.cc:98
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:399
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1107
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:910
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1216
size_type last() const
Definition: index_set.h:1097
const IndexSet * index_set
Definition: index_set.h:601
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1120
unsigned int n_intervals() const
Definition: index_set.h:1857
IndexSet tensor_product(const IndexSet &other) const
Definition: index_set.cc:330
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1243
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1982
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:1673
IndexSet operator &(const IndexSet &is) const
void set_size(const size_type size)
Definition: index_set.h:1622
ElementIterator & operator++()
Definition: index_set.h:1363
unsigned int global_dof_index
Definition: types.h:76
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:1026
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:1642
ElementIterator end() const
Definition: index_set.h:1083
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:916
static std::size_t memory_consumption()
Definition: index_set.h:929
IntervalIterator begin_intervals() const
Definition: index_set.h:1589
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1013
size_type begin
Definition: index_set.h:877
bool is_contiguous() const
Definition: index_set.h:1815
Threads::Mutex compress_mutex
Definition: index_set.h:985
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:398
size_type pop_front()
Definition: index_set.cc:360
static const char N
IntervalIterator end_intervals() const
Definition: index_set.h:1601
bool is_compressed
Definition: index_set.h:962
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1132
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1383
ElementIterator begin() const
Definition: index_set.h:1074
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:1765
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:1323
ElementIterator at(const size_type global_index) const
Definition: index_set.h:1530
const types::global_dof_index invalid_dof_index
Definition: types.h:211
std::ptrdiff_t difference_type
Definition: index_set.h:793
std::ptrdiff_t difference_type
Definition: index_set.h:702
ElementIterator begin() const
Definition: index_set.h:1518
void read(std::istream &in)
Definition: index_set.cc:451
size_type n_elements() const
Definition: index_set.h:1832
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1880
bool is_empty() const
Definition: index_set.h:1824
size_type largest_range
Definition: index_set.h:979
IntervalIterator & operator++()
Definition: index_set.h:1180
size_type nth_index_in_set
Definition: index_set.h:880
static ::ExceptionBase & ExcInternalError()
size_type operator*() const
Definition: index_set.h:1310
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1225