Reference documentation for deal.II version GIT 5743fb25e5 2023-09-22 04:10: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.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_index_set_h
17 #define dealii_index_set_h
18 
19 #include <deal.II/base/config.h>
20 
22 #include <deal.II/base/mpi_stub.h>
23 #include <deal.II/base/mutex.h>
24 
25 #include <boost/container/small_vector.hpp>
26 
27 #include <algorithm>
28 #include <vector>
29 
30 
31 #ifdef DEAL_II_WITH_TRILINOS
32 # include <Epetra_Map.h>
33 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
34 # include <Tpetra_Map.hpp>
35 # endif
36 #endif
37 
38 #ifdef DEAL_II_WITH_PETSC
39 # include <petscis.h>
40 #endif
41 
43 
66 class IndexSet
67 {
68 public:
69  // forward declarations:
70  class ElementIterator;
71  class IntervalIterator;
72 
78 
96  using value_type = signed int;
97 
98 
102  IndexSet();
103 
107  explicit IndexSet(const size_type size);
108 
112  IndexSet(const IndexSet &) = default;
113 
117  IndexSet &
118  operator=(const IndexSet &) = default;
119 
124  IndexSet(IndexSet &&is) noexcept;
125 
130  IndexSet &
131  operator=(IndexSet &&is) noexcept;
132 
133 #ifdef DEAL_II_WITH_TRILINOS
137  explicit IndexSet(const Epetra_BlockMap &map);
138 #endif
139 
144  void
145  clear();
146 
153  void
154  set_size(const size_type size);
155 
163  size_type
164  size() const;
165 
172  void
173  add_range(const size_type begin, const size_type end);
174 
183  void
184  add_index(const size_type index);
185 
210  template <typename ForwardIterator>
211  void
212  add_indices(const ForwardIterator &begin, const ForwardIterator &end);
213 
229  void
230  add_indices(const IndexSet &other, const size_type offset = 0);
231 
235  bool
236  is_element(const size_type index) const;
237 
242  bool
243  is_contiguous() const;
244 
249  bool
250  is_empty() const;
251 
261  bool
262  is_ascending_and_one_to_one(const MPI_Comm communicator) const;
263 
267  size_type
268  n_elements() const;
269 
275  size_type
276  nth_index_in_set(const size_type local_index) const;
277 
284  size_type
286 
295  unsigned int
296  n_intervals() const;
297 
309  size_type
311 
316  void
317  compress() const;
318 
333  bool
334  operator==(const IndexSet &is) const;
335 
350  bool
351  operator!=(const IndexSet &is) const;
352 
359  IndexSet
360  operator&(const IndexSet &is) const;
361 
375  IndexSet
376  get_view(const size_type begin, const size_type end) const;
377 
383  std::vector<IndexSet>
385  const std::vector<types::global_dof_index> &n_indices_per_block) const;
386 
392  void
393  subtract_set(const IndexSet &other);
394 
411  IndexSet
412  tensor_product(const IndexSet &other) const;
413 
418  size_type
419  pop_back();
420 
425  size_type
426  pop_front();
427 
436  std::vector<size_type>
437  get_index_vector() const;
438 
446  void
447  fill_index_vector(std::vector<size_type> &indices) const;
448 
461  template <typename VectorType>
462  void
463  fill_binary_vector(VectorType &vector) const;
464 
472  bool
473  is_subset_of(const IndexSet &other) const;
474 
479  template <typename StreamType>
480  void
481  print(StreamType &out) const;
482 
487  void
488  write(std::ostream &out) const;
489 
494  void
495  read(std::istream &in);
496 
501  void
502  block_write(std::ostream &out) const;
503 
508  void
509  block_read(std::istream &in);
510 
511 #ifdef DEAL_II_WITH_TRILINOS
540  Epetra_Map
541  make_trilinos_map(const MPI_Comm communicator = MPI_COMM_WORLD,
542  const bool overlapping = false) const;
543 
544 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
545  Tpetra::Map<int, types::signed_global_dof_index>
546  make_tpetra_map(const MPI_Comm communicator = MPI_COMM_WORLD,
547  const bool overlapping = false) const;
548 # endif
549 #endif
550 
551 #ifdef DEAL_II_WITH_PETSC
552  IS
553  make_petsc_is(const MPI_Comm communicator = MPI_COMM_WORLD) const;
554 #endif
555 
556 
561  std::size_t
562  memory_consumption() const;
563 
565  size_type,
566  << "The global index " << arg1
567  << " is not an element of this set.");
568 
574  template <class Archive>
575  void
576  serialize(Archive &ar, const unsigned int version);
577 
578 
590  {
591  public:
596  IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
597 
601  explicit IntervalAccessor(const IndexSet *idxset);
602 
606  size_type
607  n_elements() const;
608 
612  bool
613  is_valid() const;
614 
619  begin() const;
620 
626  end() const;
627 
631  size_type
632  last() const;
633 
634  private:
638  IntervalAccessor(const IntervalAccessor &other);
643  operator=(const IntervalAccessor &other);
644 
648  bool
649  operator==(const IntervalAccessor &other) const;
653  bool
654  operator<(const IntervalAccessor &other) const;
659  void
660  advance();
665 
671 
672  friend class IntervalIterator;
673  };
674 
680  {
681  public:
686  IntervalIterator(const IndexSet *idxset, const size_type range_idx);
687 
691  explicit IntervalIterator(const IndexSet *idxset);
692 
697 
701  IntervalIterator(const IntervalIterator &other) = default;
702 
707  operator=(const IntervalIterator &other) = default;
708 
713  operator++();
714 
719  operator++(int);
720 
724  const IntervalAccessor &
725  operator*() const;
726 
730  const IntervalAccessor *
731  operator->() const;
732 
736  bool
737  operator==(const IntervalIterator &) const;
738 
742  bool
743  operator!=(const IntervalIterator &) const;
744 
748  bool
749  operator<(const IntervalIterator &) const;
750 
757  int
758  operator-(const IntervalIterator &p) const;
759 
765  using iterator_category = std::forward_iterator_tag;
767  using difference_type = std::ptrdiff_t;
770 
771  private:
776  };
777 
783  {
784  public:
789  ElementIterator(const IndexSet *idxset,
790  const size_type range_idx,
791  const size_type index);
792 
796  explicit ElementIterator(const IndexSet *idxset);
797 
802  size_type
803  operator*() const;
804 
808  bool
809  is_valid() const;
810 
815  operator++();
816 
821  operator++(int);
822 
826  bool
827  operator==(const ElementIterator &) const;
828 
832  bool
833  operator!=(const ElementIterator &) const;
834 
838  bool
839  operator<(const ElementIterator &) const;
840 
849  std::ptrdiff_t
850  operator-(const ElementIterator &p) const;
851 
857  using iterator_category = std::forward_iterator_tag;
859  using difference_type = std::ptrdiff_t;
860  using pointer = size_type *;
861  using reference = size_type &;
862 
863  private:
867  void
868  advance();
869 
882  };
883 
889  begin() const;
890 
906  at(const size_type global_index) const;
907 
913  end() const;
914 
919  begin_intervals() const;
920 
926  end_intervals() const;
927 
932 private:
941  struct Range
942  {
945 
947 
956  Range();
957 
965  Range(const size_type i1, const size_type i2);
966 
967  friend inline bool
968  operator<(const Range &range_1, const Range &range_2)
969  {
970  return (
971  (range_1.begin < range_2.begin) ||
972  ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
973  }
974 
975  static bool
977  {
978  return x.end < y.end;
979  }
980 
981  static bool
983  {
984  return (x.nth_index_in_set + (x.end - x.begin) <
985  y.nth_index_in_set + (y.end - y.begin));
986  }
987 
988  friend inline bool
989  operator==(const Range &range_1, const Range &range_2)
990  {
991  return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
992  }
993 
994  static std::size_t
996  {
997  return sizeof(Range);
998  }
999 
1005  template <class Archive>
1006  void
1007  serialize(Archive &ar, const unsigned int version);
1008  };
1009 
1018  mutable std::vector<Range> ranges;
1019 
1028  mutable bool is_compressed;
1029 
1035 
1046 
1052 
1056  void
1057  do_compress() const;
1058 
1065  bool
1066  is_element_binary_search(const size_type local_index) const;
1067 
1074  size_type
1075  nth_index_in_set_binary_search(const size_type local_index) const;
1076 
1083  size_type
1085 
1091  void
1092  add_range_lower_bound(const Range &range);
1093 
1097  void
1099  boost::container::small_vector<std::pair<size_type, size_type>, 200>
1100  &tmp_ranges,
1101  const bool ranges_are_sorted);
1102 };
1103 
1104 
1105 
1123 inline IndexSet
1125 {
1126  IndexSet is(N);
1127  is.add_range(0, N);
1128  is.compress();
1129  return is;
1130 }
1131 
1132 /* ------------------ inline functions ------------------ */
1133 
1134 
1135 /* IntervalAccessor */
1136 
1138  const IndexSet *idxset,
1139  const IndexSet::size_type range_idx)
1140  : index_set(idxset)
1141  , range_idx(range_idx)
1142 {
1143  Assert(range_idx < idxset->n_intervals(),
1144  ExcInternalError("Invalid range index"));
1145 }
1146 
1147 
1148 
1150  : index_set(idxset)
1151  , range_idx(numbers::invalid_dof_index)
1152 {}
1153 
1154 
1155 
1157  const IndexSet::IntervalAccessor &other)
1158  : index_set(other.index_set)
1159  , range_idx(other.range_idx)
1160 {
1162  ExcMessage("invalid iterator"));
1163 }
1164 
1165 
1166 
1167 inline IndexSet::size_type
1169 {
1170  Assert(is_valid(), ExcMessage("invalid iterator"));
1171  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1172 }
1173 
1174 
1175 
1176 inline bool
1178 {
1179  return index_set != nullptr && range_idx < index_set->n_intervals();
1180 }
1181 
1182 
1183 
1186 {
1187  Assert(is_valid(), ExcMessage("invalid iterator"));
1188  return {index_set, range_idx, index_set->ranges[range_idx].begin};
1189 }
1190 
1191 
1192 
1195 {
1196  Assert(is_valid(), ExcMessage("invalid iterator"));
1197 
1198  // point to first index in next interval unless we are the last interval.
1199  if (range_idx < index_set->ranges.size() - 1)
1200  return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1201  else
1202  return index_set->end();
1203 }
1204 
1205 
1206 
1207 inline IndexSet::size_type
1209 {
1210  Assert(is_valid(), ExcMessage("invalid iterator"));
1211 
1212  return index_set->ranges[range_idx].end - 1;
1213 }
1214 
1215 
1216 
1219 {
1220  index_set = other.index_set;
1221  range_idx = other.range_idx;
1222  Assert(range_idx == numbers::invalid_dof_index || is_valid(),
1223  ExcMessage("invalid iterator"));
1224  return *this;
1225 }
1226 
1227 
1228 
1229 inline bool
1231  const IndexSet::IntervalAccessor &other) const
1232 {
1233  Assert(index_set == other.index_set,
1234  ExcMessage(
1235  "Can not compare accessors pointing to different IndexSets"));
1236  return range_idx == other.range_idx;
1237 }
1238 
1239 
1240 
1241 inline bool
1243  const IndexSet::IntervalAccessor &other) const
1244 {
1245  Assert(index_set == other.index_set,
1246  ExcMessage(
1247  "Can not compare accessors pointing to different IndexSets"));
1248  return range_idx < other.range_idx;
1249 }
1250 
1251 
1252 
1253 inline void
1255 {
1256  Assert(
1257  is_valid(),
1258  ExcMessage(
1259  "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1260  ++range_idx;
1261 
1262  // set ourselves to invalid if we walk off the end
1263  if (range_idx >= index_set->ranges.size())
1264  range_idx = numbers::invalid_dof_index;
1265 }
1266 
1267 
1268 /* IntervalIterator */
1269 
1271  const IndexSet *idxset,
1272  const IndexSet::size_type range_idx)
1273  : accessor(idxset, range_idx)
1274 {}
1275 
1276 
1277 
1279  : accessor(nullptr)
1280 {}
1281 
1282 
1283 
1285  : accessor(idxset)
1286 {}
1287 
1288 
1289 
1292 {
1293  accessor.advance();
1294  return *this;
1295 }
1296 
1297 
1298 
1301 {
1302  const IndexSet::IntervalIterator iter = *this;
1303  accessor.advance();
1304  return iter;
1305 }
1306 
1307 
1308 
1309 inline const IndexSet::IntervalAccessor &
1311 {
1312  return accessor;
1313 }
1314 
1315 
1316 
1317 inline const IndexSet::IntervalAccessor *
1319 {
1320  return &accessor;
1321 }
1322 
1323 
1324 
1325 inline bool
1327  const IndexSet::IntervalIterator &other) const
1328 {
1329  return accessor == other.accessor;
1330 }
1331 
1332 
1333 
1334 inline bool
1336  const IndexSet::IntervalIterator &other) const
1337 {
1338  return !(*this == other);
1339 }
1340 
1341 
1342 
1343 inline bool
1345  const IndexSet::IntervalIterator &other) const
1346 {
1347  return accessor < other.accessor;
1348 }
1349 
1350 
1351 
1352 inline int
1354  const IndexSet::IntervalIterator &other) const
1355 {
1356  Assert(accessor.index_set == other.accessor.index_set,
1357  ExcMessage(
1358  "Can not compare iterators belonging to different IndexSets"));
1359 
1360  const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index) ?
1361  accessor.index_set->ranges.size() :
1362  accessor.range_idx;
1363  const size_type rhs =
1365  accessor.index_set->ranges.size() :
1366  other.accessor.range_idx;
1367 
1368  if (lhs > rhs)
1369  return static_cast<int>(lhs - rhs);
1370  else
1371  return -static_cast<int>(rhs - lhs);
1372 }
1373 
1374 
1375 
1376 /* ElementIterator */
1377 
1379  const IndexSet *idxset,
1380  const IndexSet::size_type range_idx,
1381  const IndexSet::size_type index)
1382  : index_set(idxset)
1383  , range_idx(range_idx)
1384  , idx(index)
1385 {
1386  Assert(range_idx < index_set->ranges.size(),
1387  ExcMessage(
1388  "Invalid range index for IndexSet::ElementIterator constructor."));
1389  Assert(
1390  idx >= index_set->ranges[range_idx].begin &&
1391  idx < index_set->ranges[range_idx].end,
1393  "Invalid index argument for IndexSet::ElementIterator constructor."));
1394 }
1395 
1396 
1397 
1399  : index_set(idxset)
1400  , range_idx(numbers::invalid_dof_index)
1401  , idx(numbers::invalid_dof_index)
1402 {}
1403 
1404 
1405 
1406 inline bool
1408 {
1409  Assert((range_idx == numbers::invalid_dof_index &&
1410  idx == numbers::invalid_dof_index) ||
1411  (range_idx < index_set->ranges.size() &&
1412  idx < index_set->ranges[range_idx].end),
1413  ExcInternalError("Invalid ElementIterator state."));
1414 
1415  return (range_idx < index_set->ranges.size() &&
1416  idx < index_set->ranges[range_idx].end);
1417 }
1418 
1419 
1420 
1421 inline IndexSet::size_type
1423 {
1424  Assert(
1425  is_valid(),
1426  ExcMessage(
1427  "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1428  return idx;
1429 }
1430 
1431 
1432 
1433 inline bool
1435  const IndexSet::ElementIterator &other) const
1436 {
1437  Assert(index_set == other.index_set,
1438  ExcMessage(
1439  "Can not compare iterators belonging to different IndexSets"));
1440  return range_idx == other.range_idx && idx == other.idx;
1441 }
1442 
1443 
1444 
1445 inline void
1447 {
1448  Assert(
1449  is_valid(),
1450  ExcMessage(
1451  "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1452  if (idx < index_set->ranges[range_idx].end)
1453  ++idx;
1454  // end of this range?
1455  if (idx == index_set->ranges[range_idx].end)
1456  {
1457  // point to first element in next interval if possible
1458  if (range_idx < index_set->ranges.size() - 1)
1459  {
1460  ++range_idx;
1461  idx = index_set->ranges[range_idx].begin;
1462  }
1463  else
1464  {
1465  // we just fell off the end, set to invalid:
1466  range_idx = numbers::invalid_dof_index;
1468  }
1469  }
1470 }
1471 
1472 
1473 
1476 {
1477  advance();
1478  return *this;
1479 }
1480 
1481 
1482 
1485 {
1486  const IndexSet::ElementIterator it = *this;
1487  advance();
1488  return it;
1489 }
1490 
1491 
1492 
1493 inline bool
1495  const IndexSet::ElementIterator &other) const
1496 {
1497  return !(*this == other);
1498 }
1499 
1500 
1501 
1502 inline bool
1504  const IndexSet::ElementIterator &other) const
1505 {
1506  Assert(index_set == other.index_set,
1507  ExcMessage(
1508  "Can not compare iterators belonging to different IndexSets"));
1509  return range_idx < other.range_idx ||
1510  (range_idx == other.range_idx && idx < other.idx);
1511 }
1512 
1513 
1514 
1515 inline std::ptrdiff_t
1517  const IndexSet::ElementIterator &other) const
1518 {
1519  Assert(index_set == other.index_set,
1520  ExcMessage(
1521  "Can not compare iterators belonging to different IndexSets"));
1522  if (*this == other)
1523  return 0;
1524  if (!(*this < other))
1525  return -(other - *this);
1526 
1527  // only other can be equal to end() because of the checks above.
1528  Assert(is_valid(), ExcInternalError());
1529 
1530  // Note: we now compute how far advance *this in "*this < other" to get other,
1531  // so we need to return -c at the end.
1532 
1533  // first finish the current range:
1534  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1535 
1536  // now walk in steps of ranges (need to start one behind our current one):
1537  for (size_type range = range_idx + 1;
1538  range < index_set->ranges.size() && range <= other.range_idx;
1539  ++range)
1540  c += index_set->ranges[range].end - index_set->ranges[range].begin;
1541 
1542  Assert(
1543  other.range_idx < index_set->ranges.size() ||
1545  ExcMessage(
1546  "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1547 
1548  // We might have walked too far because we went until the end of
1549  // other.range_idx, so walk backwards to other.idx:
1551  c -= index_set->ranges[other.range_idx].end - other.idx;
1552 
1553  return -c;
1554 }
1555 
1556 
1557 /* Range */
1558 
1563 {}
1564 
1565 
1566 
1567 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1568  : begin(i1)
1569  , end(i2)
1571 {}
1572 
1573 
1574 
1575 /* IndexSet itself */
1576 
1578  : is_compressed(true)
1579  , index_space_size(0)
1581 {}
1582 
1583 
1584 
1585 inline IndexSet::IndexSet(const size_type size)
1586  : is_compressed(true)
1587  , index_space_size(size)
1588  , largest_range(numbers::invalid_unsigned_int)
1589 {}
1590 
1591 
1592 
1593 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1594  : ranges(std::move(is.ranges))
1595  , is_compressed(is.is_compressed)
1596  , index_space_size(is.index_space_size)
1597  , largest_range(is.largest_range)
1598 {
1599  is.ranges.clear();
1600  is.is_compressed = true;
1601  is.index_space_size = 0;
1602  is.largest_range = numbers::invalid_unsigned_int;
1603 
1604  compress();
1605 }
1606 
1607 
1608 
1609 inline IndexSet &
1611 {
1612  ranges = std::move(is.ranges);
1613  is_compressed = is.is_compressed;
1614  index_space_size = is.index_space_size;
1615  largest_range = is.largest_range;
1616 
1617  is.ranges.clear();
1618  is.is_compressed = true;
1619  is.index_space_size = 0;
1620  is.largest_range = numbers::invalid_unsigned_int;
1621 
1622  compress();
1623 
1624  return *this;
1625 }
1626 
1627 
1628 
1631 {
1632  compress();
1633  if (ranges.size() > 0)
1634  return {this, 0, ranges[0].begin};
1635  else
1636  return end();
1637 }
1638 
1639 
1640 
1643 {
1644  compress();
1645  return IndexSet::ElementIterator(this);
1646 }
1647 
1648 
1649 
1652 {
1653  compress();
1654  if (ranges.size() > 0)
1655  return IndexSet::IntervalIterator(this, 0);
1656  else
1657  return end_intervals();
1658 }
1659 
1660 
1661 
1664 {
1665  compress();
1666  return IndexSet::IntervalIterator(this);
1667 }
1668 
1669 
1670 
1671 inline void
1673 {
1674  // reset so that there are no indices in the set any more; however,
1675  // as documented, the index set retains its size
1676  ranges.clear();
1677  is_compressed = true;
1679 }
1680 
1681 
1682 
1683 inline void
1685 {
1686  Assert(ranges.empty(),
1687  ExcMessage("This function can only be called if the current "
1688  "object does not yet contain any elements."));
1689  index_space_size = sz;
1690  is_compressed = true;
1691 }
1692 
1693 
1694 
1695 inline IndexSet::size_type
1697 {
1698  return index_space_size;
1699 }
1700 
1701 
1702 
1703 inline void
1705 {
1706  if (is_compressed == true)
1707  return;
1708 
1709  do_compress();
1710 }
1711 
1712 
1713 
1714 inline void
1716 {
1717  add_range(index, index + 1);
1718 }
1719 
1720 
1721 
1722 inline void
1724 {
1726  ((begin == index_space_size) && (end == index_space_size)),
1727  ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1729  ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1730  AssertIndexRange(begin, end + 1);
1731 
1732  if (begin != end)
1733  {
1734  // the new index might be larger than the last index present in the
1735  // ranges. Then we can skip the binary search
1736  if (ranges.empty() || begin > ranges.back().end)
1737  ranges.emplace_back(begin, end);
1738  else if (begin == ranges.back().end)
1739  ranges.back().end = end;
1740  else
1742 
1743  is_compressed = false;
1744  }
1745 }
1746 
1747 
1748 
1749 template <typename ForwardIterator>
1750 inline void
1751 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1752 {
1753  if (begin == end)
1754  return;
1755 
1756  // identify ranges in the given iterator range by checking whether some
1757  // indices happen to be consecutive. to avoid quadratic complexity when
1758  // calling add_range many times (as add_range() going into the middle of an
1759  // already existing range must shift entries around), we first collect a
1760  // vector of ranges.
1761  boost::container::small_vector<std::pair<size_type, size_type>, 200>
1762  tmp_ranges;
1763  bool ranges_are_sorted = true;
1764  for (ForwardIterator p = begin; p != end;)
1765  {
1766  // Starting with the current iterator 'p', find an iterator
1767  // 'q' so that the indices pointed to by the iterators in
1768  // the range [p,q) are consecutive. These indices then form
1769  // a range that is contiguous, and that can be added all
1770  // at once.
1771  const size_type begin_index = *p;
1772  size_type end_index = begin_index + 1;
1773 
1774  // Start looking at the position after 'p', and keep iterating while
1775  // 'q' points to a duplicate of 'p':
1776  ForwardIterator q = p;
1777  ++q;
1778  while ((q != end) && (*q == *p))
1779  ++q;
1780 
1781  // Now we know that 'q' is either past the end, or points to a value
1782  // other than 'p'. If it points to 'end_index', we are still good with
1783  // a contiguous range; then increment the end index of that range, and
1784  // move to the next iterator that is not a duplicate of what
1785  // we were just looking at:
1786  while ((q != end) && (static_cast<size_type>(*q) == end_index))
1787  {
1788  ++q;
1789  while ((q != end) && (static_cast<size_type>(*q) == end_index))
1790  ++q;
1791 
1792  ++end_index;
1793  }
1794 
1795  // Add this range:
1796  tmp_ranges.emplace_back(begin_index, end_index);
1797 
1798  // Then move on to the next element in the input range.
1799  // If the starting index of the next go-around of the for loop is less
1800  // than the end index of the one just identified, then we will have at
1801  // least one pair of ranges that are not sorted, and consequently the
1802  // whole collection of ranges is not sorted.
1803  p = q;
1804  if ((p != end) && (static_cast<size_type>(*p) < end_index))
1805  ranges_are_sorted = false;
1806  }
1807 
1808  add_ranges_internal(tmp_ranges, ranges_are_sorted);
1809 }
1810 
1811 
1812 
1813 inline bool
1815 {
1816  if (ranges.empty() == false)
1817  {
1818  compress();
1819 
1820  // fast check whether the index is in the largest range
1822  if (index >= ranges[largest_range].begin &&
1823  index < ranges[largest_range].end)
1824  return true;
1825  else if (ranges.size() > 1)
1826  return is_element_binary_search(index);
1827  else
1828  return false;
1829  }
1830  else
1831  return false;
1832 }
1833 
1834 
1835 
1836 inline bool
1838 {
1839  compress();
1840  return (ranges.size() <= 1);
1841 }
1842 
1843 
1844 
1845 inline bool
1847 {
1848  return ranges.empty();
1849 }
1850 
1851 
1852 
1853 inline IndexSet::size_type
1855 {
1856  // make sure we have non-overlapping ranges
1857  compress();
1858 
1859  size_type v = 0;
1860  if (!ranges.empty())
1861  {
1862  const Range &r = ranges.back();
1863  v = r.nth_index_in_set + r.end - r.begin;
1864  }
1865 
1866 #ifdef DEBUG
1867  size_type s = 0;
1868  for (const auto &range : ranges)
1869  s += (range.end - range.begin);
1870  Assert(s == v, ExcInternalError());
1871 #endif
1872 
1873  return v;
1874 }
1875 
1876 
1877 
1878 inline unsigned int
1880 {
1881  compress();
1882  return ranges.size();
1883 }
1884 
1885 
1886 
1887 inline IndexSet::size_type
1889 {
1890  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1891 
1892  compress();
1893  const std::vector<Range>::const_iterator main_range =
1894  ranges.begin() + largest_range;
1895 
1896  return main_range->nth_index_in_set;
1897 }
1898 
1899 
1900 
1901 inline IndexSet::size_type
1903 {
1905 
1906  compress();
1907 
1908  // first check whether the index is in the largest range
1910  const auto main_range = ranges.begin() + largest_range;
1911  if (n >= main_range->nth_index_in_set &&
1912  n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1913  return main_range->begin + (n - main_range->nth_index_in_set);
1914  else
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  if (n >= ranges[largest_range].begin && n < ranges[largest_range].end)
1936  return (n - ranges[largest_range].begin) +
1937  ranges[largest_range].nth_index_in_set;
1938  else if (ranges.size() > 1)
1940  else
1942 }
1943 
1944 
1945 
1946 inline bool
1948 {
1949  // If one of the two index sets has size zero, the other one has to
1950  // have size zero as well:
1951  if (size() == 0)
1952  return (is.size() == 0);
1953  if (is.size() == 0)
1954  return (size() == 0);
1955 
1956  // Otherwise, they must have the same size (see the documentation):
1957  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1958 
1959  compress();
1960  is.compress();
1961 
1962  return (ranges == is.ranges);
1963 }
1964 
1965 
1966 
1967 inline bool
1969 {
1970  // If one of the two index sets has size zero, the other one has to
1971  // have a non-zero size for inequality:
1972  if (size() == 0)
1973  return (is.size() != 0);
1974  if (is.size() == 0)
1975  return (size() != 0);
1976 
1977  // Otherwise, they must have the same size (see the documentation):
1978  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1979 
1980  compress();
1981  is.compress();
1982 
1983  return (ranges != is.ranges);
1984 }
1985 
1986 
1987 
1988 template <typename Vector>
1989 void
1991 {
1992  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
1993 
1994  compress();
1995  // first fill all elements of the vector with zeroes.
1996  std::fill(vector.begin(), vector.end(), 0);
1997 
1998  // then write ones into the elements whose indices are contained in the
1999  // index set
2000  for (const auto &range : ranges)
2001  for (size_type i = range.begin; i < range.end; ++i)
2002  vector[i] = 1;
2003 }
2004 
2005 
2006 
2007 template <typename StreamType>
2008 inline void
2009 IndexSet::print(StreamType &out) const
2010 {
2011  compress();
2012  out << '{';
2013  std::vector<Range>::const_iterator p;
2014  for (p = ranges.begin(); p != ranges.end(); ++p)
2015  {
2016  if (p->end - p->begin == 1)
2017  out << p->begin;
2018  else
2019  out << '[' << p->begin << ',' << p->end - 1 << ']';
2020 
2021  if (p != --ranges.end())
2022  out << ", ";
2023  }
2024  out << '}' << std::endl;
2025 }
2026 
2027 
2028 
2029 template <class Archive>
2030 inline void
2031 IndexSet::Range::serialize(Archive &ar, const unsigned int)
2032 {
2033  ar &begin &end &nth_index_in_set;
2034 }
2035 
2036 
2037 
2038 template <class Archive>
2039 inline void
2040 IndexSet::serialize(Archive &ar, const unsigned int)
2041 {
2043 }
2044 
2046 
2047 #endif
size_type operator*() const
Definition: index_set.h:1422
bool operator==(const ElementIterator &) const
Definition: index_set.h:1434
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition: index_set.h:1378
std::ptrdiff_t difference_type
Definition: index_set.h:859
ElementIterator & operator++()
Definition: index_set.h:1475
bool operator<(const ElementIterator &) const
Definition: index_set.h:1503
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1494
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1516
std::forward_iterator_tag iterator_category
Definition: index_set.h:857
const IndexSet * index_set
Definition: index_set.h:873
size_type last() const
Definition: index_set.h:1208
ElementIterator end() const
Definition: index_set.h:1194
size_type n_elements() const
Definition: index_set.h:1168
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1230
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1218
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1242
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:1137
const IndexSet * index_set
Definition: index_set.h:664
ElementIterator begin() const
Definition: index_set.h:1185
std::ptrdiff_t difference_type
Definition: index_set.h:767
const IntervalAccessor & operator*() const
Definition: index_set.h:1310
IntervalIterator(const IntervalIterator &other)=default
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1326
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1353
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1344
std::forward_iterator_tag iterator_category
Definition: index_set.h:765
IntervalIterator & operator++()
Definition: index_set.h:1291
IntervalIterator & operator=(const IntervalIterator &other)=default
const IntervalAccessor * operator->() const
Definition: index_set.h:1318
IntervalAccessor accessor
Definition: index_set.h:775
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1335
bool is_subset_of(const IndexSet &other) const
Definition: index_set.cc:521
bool is_element_binary_search(const size_type local_index) const
Definition: index_set.cc:609
size_type index_within_set_binary_search(const size_type global_index) const
Definition: index_set.cc:662
size_type largest_range
Definition: index_set.h:1045
IS make_petsc_is(const MPI_Comm communicator=MPI_COMM_WORLD) const
Definition: index_set.cc:902
bool is_ascending_and_one_to_one(const MPI_Comm communicator) const
Definition: index_set.cc:922
bool is_contiguous() const
Definition: index_set.h:1837
unsigned int n_intervals() const
Definition: index_set.h:1879
IndexSet(const IndexSet &)=default
IndexSet & operator=(const IndexSet &)=default
IntervalIterator end_intervals() const
Definition: index_set.h:1663
void do_compress() const
Definition: index_set.cc:100
ElementIterator at(const size_type global_index) const
Definition: index_set.cc:684
size_type size() const
Definition: index_set.h:1696
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:753
bool is_empty() const
Definition: index_set.h:1846
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1921
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
Definition: index_set.cc:261
size_type n_elements() const
Definition: index_set.h:1854
bool operator==(const IndexSet &is) const
Definition: index_set.h:1947
void add_range_lower_bound(const Range &range)
Definition: index_set.cc:409
bool is_element(const size_type index) const
Definition: index_set.h:1814
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2040
ElementIterator begin() const
Definition: index_set.h:1630
size_type pop_front()
Definition: index_set.cc:387
signed int value_type
Definition: index_set.h:96
void set_size(const size_type size)
Definition: index_set.h:1684
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1968
void read(std::istream &in)
Definition: index_set.cc:553
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:833
void clear()
Definition: index_set.h:1672
size_type largest_range_starting_index() const
Definition: index_set.h:1888
IndexSet tensor_product(const IndexSet &other) const
Definition: index_set.cc:357
void add_index(const size_type index)
Definition: index_set.h:1715
void write(std::ostream &out) const
Definition: index_set.cc:538
void block_read(std::istream &in)
Definition: index_set.cc:589
bool is_compressed
Definition: index_set.h:1028
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:425
IntervalIterator begin_intervals() const
Definition: index_set.h:1651
Tpetra::Map< int, types::signed_global_dof_index > make_tpetra_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:764
void fill_binary_vector(VectorType &vector) const
std::vector< Range > ranges
Definition: index_set.h:1018
void subtract_set(const IndexSet &other)
Definition: index_set.cc:288
ElementIterator end() const
Definition: index_set.h:1642
Threads::Mutex compress_mutex
Definition: index_set.h:1051
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1124
size_type index_space_size
Definition: index_set.h:1034
void block_write(std::ostream &out) const
Definition: index_set.cc:575
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:234
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1723
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1902
std::size_t memory_consumption() const
Definition: index_set.cc:988
void print(StreamType &out) const
Definition: index_set.h:2009
size_type nth_index_in_set_binary_search(const size_type local_index) const
Definition: index_set.cc:645
void compress() const
Definition: index_set.h:1704
size_type pop_back()
Definition: index_set.cc:369
std::vector< size_type > get_index_vector() const
Definition: index_set.cc:734
types::global_dof_index size_type
Definition: index_set.h:77
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1751
IndexSet operator&(const IndexSet &is) const
Definition: vector.h:110
virtual size_type size() const override
iterator end()
iterator begin()
#define DEAL_II_DEPRECATED
Definition: config.h:178
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcIndexNotPresent(size_type arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:512
static ::ExceptionBase & ExcMessage(std::string arg1)
static const char N
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)
std::string compress(const std::string &input)
Definition: utilities.cc:390
static const unsigned int invalid_unsigned_int
Definition: types.h:213
const types::global_dof_index invalid_dof_index
Definition: types.h:233
unsigned int global_dof_index
Definition: types.h:82
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:976
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:982
friend bool operator==(const Range &range_1, const Range &range_2)
Definition: index_set.h:989
size_type end
Definition: index_set.h:944
size_type nth_index_in_set
Definition: index_set.h:946
static std::size_t memory_consumption()
Definition: index_set.h:995
size_type begin
Definition: index_set.h:943
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2031
friend bool operator<(const Range &range_1, const Range &range_2)
Definition: index_set.h:968
void advance(std::tuple< I1, I2 > &t, const unsigned int n)