Reference documentation for deal.II version GIT 0980a66d4b 2023-03-23 20:20:03+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 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #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 
178  void
179  add_index(const size_type index);
180 
190  template <typename ForwardIterator>
191  void
192  add_indices(const ForwardIterator &begin, const ForwardIterator &end);
193 
209  void
210  add_indices(const IndexSet &other, const size_type offset = 0);
211 
215  bool
216  is_element(const size_type index) const;
217 
222  bool
223  is_contiguous() const;
224 
229  bool
230  is_empty() const;
231 
241  bool
242  is_ascending_and_one_to_one(const MPI_Comm &communicator) const;
243 
247  size_type
248  n_elements() const;
249 
255  size_type
256  nth_index_in_set(const size_type local_index) const;
257 
264  size_type
266 
275  unsigned int
276  n_intervals() const;
277 
289  size_type
291 
296  void
297  compress() const;
298 
304  bool
305  operator==(const IndexSet &is) const;
306 
312  bool
313  operator!=(const IndexSet &is) const;
314 
321  IndexSet
322  operator&(const IndexSet &is) const;
323 
337  IndexSet
338  get_view(const size_type begin, const size_type end) const;
339 
345  std::vector<IndexSet>
347  const std::vector<types::global_dof_index> &n_indices_per_block) const;
348 
354  void
355  subtract_set(const IndexSet &other);
356 
373  IndexSet
374  tensor_product(const IndexSet &other) const;
375 
380  size_type
381  pop_back();
382 
387  size_type
388  pop_front();
389 
398  std::vector<size_type>
399  get_index_vector() const;
400 
407  DEAL_II_DEPRECATED_EARLY
408  void
409  fill_index_vector(std::vector<size_type> &indices) const;
410 
423  template <typename VectorType>
424  void
425  fill_binary_vector(VectorType &vector) const;
426 
431  template <class StreamType>
432  void
433  print(StreamType &out) const;
434 
439  void
440  write(std::ostream &out) const;
441 
446  void
447  read(std::istream &in);
448 
453  void
454  block_write(std::ostream &out) const;
455 
460  void
461  block_read(std::istream &in);
462 
463 #ifdef DEAL_II_WITH_TRILINOS
492  Epetra_Map
493  make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
494  const bool overlapping = false) const;
495 
496 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
497  Tpetra::Map<int, types::signed_global_dof_index>
498  make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
499  const bool overlapping = false) const;
500 # endif
501 #endif
502 
503 #ifdef DEAL_II_WITH_PETSC
504  IS
505  make_petsc_is(const MPI_Comm &communicator = MPI_COMM_WORLD) const;
506 #endif
507 
508 
513  std::size_t
514  memory_consumption() const;
515 
517  size_type,
518  << "The global index " << arg1
519  << " is not an element of this set.");
520 
526  template <class Archive>
527  void
528  serialize(Archive &ar, const unsigned int version);
529 
530 
542  {
543  public:
548  IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
549 
553  explicit IntervalAccessor(const IndexSet *idxset);
554 
558  size_type
559  n_elements() const;
560 
564  bool
565  is_valid() const;
566 
571  begin() const;
572 
578  end() const;
579 
583  size_type
584  last() const;
585 
586  private:
590  IntervalAccessor(const IntervalAccessor &other);
595  operator=(const IntervalAccessor &other);
596 
600  bool
601  operator==(const IntervalAccessor &other) const;
605  bool
606  operator<(const IntervalAccessor &other) const;
611  void
612  advance();
617 
623 
624  friend class IntervalIterator;
625  };
626 
632  {
633  public:
638  IntervalIterator(const IndexSet *idxset, const size_type range_idx);
639 
643  explicit IntervalIterator(const IndexSet *idxset);
644 
649 
653  IntervalIterator(const IntervalIterator &other) = default;
654 
659  operator=(const IntervalIterator &other) = default;
660 
665  operator++();
666 
671  operator++(int);
672 
676  const IntervalAccessor &
677  operator*() const;
678 
682  const IntervalAccessor *
683  operator->() const;
684 
688  bool
689  operator==(const IntervalIterator &) const;
690 
694  bool
695  operator!=(const IntervalIterator &) const;
696 
700  bool
701  operator<(const IntervalIterator &) const;
702 
709  int
710  operator-(const IntervalIterator &p) const;
711 
717  using iterator_category = std::forward_iterator_tag;
719  using difference_type = std::ptrdiff_t;
722 
723  private:
728  };
729 
735  {
736  public:
741  ElementIterator(const IndexSet *idxset,
742  const size_type range_idx,
743  const size_type index);
744 
748  explicit ElementIterator(const IndexSet *idxset);
749 
754  size_type
755  operator*() const;
756 
760  bool
761  is_valid() const;
762 
767  operator++();
768 
773  operator++(int);
774 
778  bool
779  operator==(const ElementIterator &) const;
780 
784  bool
785  operator!=(const ElementIterator &) const;
786 
790  bool
791  operator<(const ElementIterator &) const;
792 
801  std::ptrdiff_t
802  operator-(const ElementIterator &p) const;
803 
809  using iterator_category = std::forward_iterator_tag;
811  using difference_type = std::ptrdiff_t;
812  using pointer = size_type *;
813  using reference = size_type &;
814 
815  private:
819  void
820  advance();
821 
834  };
835 
841  begin() const;
842 
858  at(const size_type global_index) const;
859 
865  end() const;
866 
871  begin_intervals() const;
872 
878  end_intervals() const;
879 
884 private:
893  struct Range
894  {
897 
899 
908  Range();
909 
917  Range(const size_type i1, const size_type i2);
918 
919  friend inline bool
920  operator<(const Range &range_1, const Range &range_2)
921  {
922  return (
923  (range_1.begin < range_2.begin) ||
924  ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
925  }
926 
927  static bool
929  {
930  return x.end < y.end;
931  }
932 
933  static bool
935  {
936  return (x.nth_index_in_set + (x.end - x.begin) <
937  y.nth_index_in_set + (y.end - y.begin));
938  }
939 
940  friend inline bool
941  operator==(const Range &range_1, const Range &range_2)
942  {
943  return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
944  }
945 
946  static std::size_t
948  {
949  return sizeof(Range);
950  }
951 
957  template <class Archive>
958  void
959  serialize(Archive &ar, const unsigned int version);
960  };
961 
970  mutable std::vector<Range> ranges;
971 
980  mutable bool is_compressed;
981 
987 
998 
1004 
1008  void
1009  do_compress() const;
1010 
1017  bool
1018  is_element_binary_search(const size_type local_index) const;
1019 
1026  size_type
1027  nth_index_in_set_binary_search(const size_type local_index) const;
1028 
1035  size_type
1037 
1043  void
1044  add_range_lower_bound(const Range &range);
1045 
1049  void
1051  boost::container::small_vector<std::pair<size_type, size_type>, 200>
1052  & tmp_ranges,
1053  const bool ranges_are_sorted);
1054 };
1055 
1056 
1057 
1075 inline IndexSet
1077 {
1078  IndexSet is(N);
1079  is.add_range(0, N);
1080  is.compress();
1081  return is;
1082 }
1083 
1084 /* ------------------ inline functions ------------------ */
1085 
1086 
1087 /* IntervalAccessor */
1088 
1090  const IndexSet * idxset,
1091  const IndexSet::size_type range_idx)
1092  : index_set(idxset)
1093  , range_idx(range_idx)
1094 {
1095  Assert(range_idx < idxset->n_intervals(),
1096  ExcInternalError("Invalid range index"));
1097 }
1098 
1099 
1100 
1102  : index_set(idxset)
1103  , range_idx(numbers::invalid_dof_index)
1104 {}
1105 
1106 
1107 
1109  const IndexSet::IntervalAccessor &other)
1110  : index_set(other.index_set)
1111  , range_idx(other.range_idx)
1112 {
1114  ExcMessage("invalid iterator"));
1115 }
1116 
1117 
1118 
1119 inline IndexSet::size_type
1121 {
1122  Assert(is_valid(), ExcMessage("invalid iterator"));
1123  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1124 }
1125 
1126 
1127 
1128 inline bool
1130 {
1131  return index_set != nullptr && range_idx < index_set->n_intervals();
1132 }
1133 
1134 
1135 
1138 {
1139  Assert(is_valid(), ExcMessage("invalid iterator"));
1140  return {index_set, range_idx, index_set->ranges[range_idx].begin};
1141 }
1142 
1143 
1144 
1147 {
1148  Assert(is_valid(), ExcMessage("invalid iterator"));
1149 
1150  // point to first index in next interval unless we are the last interval.
1151  if (range_idx < index_set->ranges.size() - 1)
1152  return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1153  else
1154  return index_set->end();
1155 }
1156 
1157 
1158 
1159 inline IndexSet::size_type
1161 {
1162  Assert(is_valid(), ExcMessage("invalid iterator"));
1163 
1164  return index_set->ranges[range_idx].end - 1;
1165 }
1166 
1167 
1168 
1171 {
1172  index_set = other.index_set;
1173  range_idx = other.range_idx;
1174  Assert(range_idx == numbers::invalid_dof_index || is_valid(),
1175  ExcMessage("invalid iterator"));
1176  return *this;
1177 }
1178 
1179 
1180 
1181 inline bool
1183  const IndexSet::IntervalAccessor &other) const
1184 {
1185  Assert(index_set == other.index_set,
1186  ExcMessage(
1187  "Can not compare accessors pointing to different IndexSets"));
1188  return range_idx == other.range_idx;
1189 }
1190 
1191 
1192 
1193 inline bool
1195  const IndexSet::IntervalAccessor &other) const
1196 {
1197  Assert(index_set == other.index_set,
1198  ExcMessage(
1199  "Can not compare accessors pointing to different IndexSets"));
1200  return range_idx < other.range_idx;
1201 }
1202 
1203 
1204 
1205 inline void
1207 {
1208  Assert(
1209  is_valid(),
1210  ExcMessage(
1211  "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1212  ++range_idx;
1213 
1214  // set ourselves to invalid if we walk off the end
1215  if (range_idx >= index_set->ranges.size())
1216  range_idx = numbers::invalid_dof_index;
1217 }
1218 
1219 
1220 /* IntervalIterator */
1221 
1223  const IndexSet * idxset,
1224  const IndexSet::size_type range_idx)
1225  : accessor(idxset, range_idx)
1226 {}
1227 
1228 
1229 
1231  : accessor(nullptr)
1232 {}
1233 
1234 
1235 
1237  : accessor(idxset)
1238 {}
1239 
1240 
1241 
1244 {
1245  accessor.advance();
1246  return *this;
1247 }
1248 
1249 
1250 
1253 {
1254  const IndexSet::IntervalIterator iter = *this;
1255  accessor.advance();
1256  return iter;
1257 }
1258 
1259 
1260 
1261 inline const IndexSet::IntervalAccessor &
1263 {
1264  return accessor;
1265 }
1266 
1267 
1268 
1269 inline const IndexSet::IntervalAccessor *
1271 {
1272  return &accessor;
1273 }
1274 
1275 
1276 
1277 inline bool
1279  const IndexSet::IntervalIterator &other) const
1280 {
1281  return accessor == other.accessor;
1282 }
1283 
1284 
1285 
1286 inline bool
1288  const IndexSet::IntervalIterator &other) const
1289 {
1290  return !(*this == other);
1291 }
1292 
1293 
1294 
1295 inline bool
1297  const IndexSet::IntervalIterator &other) const
1298 {
1299  return accessor < other.accessor;
1300 }
1301 
1302 
1303 
1304 inline int
1306  const IndexSet::IntervalIterator &other) const
1307 {
1308  Assert(accessor.index_set == other.accessor.index_set,
1309  ExcMessage(
1310  "Can not compare iterators belonging to different IndexSets"));
1311 
1312  const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index) ?
1313  accessor.index_set->ranges.size() :
1314  accessor.range_idx;
1315  const size_type rhs =
1317  accessor.index_set->ranges.size() :
1318  other.accessor.range_idx;
1319 
1320  if (lhs > rhs)
1321  return static_cast<int>(lhs - rhs);
1322  else
1323  return -static_cast<int>(rhs - lhs);
1324 }
1325 
1326 
1327 
1328 /* ElementIterator */
1329 
1331  const IndexSet * idxset,
1332  const IndexSet::size_type range_idx,
1333  const IndexSet::size_type index)
1334  : index_set(idxset)
1335  , range_idx(range_idx)
1336  , idx(index)
1337 {
1338  Assert(range_idx < index_set->ranges.size(),
1339  ExcMessage(
1340  "Invalid range index for IndexSet::ElementIterator constructor."));
1341  Assert(
1342  idx >= index_set->ranges[range_idx].begin &&
1343  idx < index_set->ranges[range_idx].end,
1345  "Invalid index argument for IndexSet::ElementIterator constructor."));
1346 }
1347 
1348 
1349 
1351  : index_set(idxset)
1352  , range_idx(numbers::invalid_dof_index)
1353  , idx(numbers::invalid_dof_index)
1354 {}
1355 
1356 
1357 
1358 inline bool
1360 {
1361  Assert((range_idx == numbers::invalid_dof_index &&
1362  idx == numbers::invalid_dof_index) ||
1363  (range_idx < index_set->ranges.size() &&
1364  idx < index_set->ranges[range_idx].end),
1365  ExcInternalError("Invalid ElementIterator state."));
1366 
1367  return (range_idx < index_set->ranges.size() &&
1368  idx < index_set->ranges[range_idx].end);
1369 }
1370 
1371 
1372 
1373 inline IndexSet::size_type
1375 {
1376  Assert(
1377  is_valid(),
1378  ExcMessage(
1379  "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1380  return idx;
1381 }
1382 
1383 
1384 
1385 inline bool
1387  const IndexSet::ElementIterator &other) const
1388 {
1389  Assert(index_set == other.index_set,
1390  ExcMessage(
1391  "Can not compare iterators belonging to different IndexSets"));
1392  return range_idx == other.range_idx && idx == other.idx;
1393 }
1394 
1395 
1396 
1397 inline void
1399 {
1400  Assert(
1401  is_valid(),
1402  ExcMessage(
1403  "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1404  if (idx < index_set->ranges[range_idx].end)
1405  ++idx;
1406  // end of this range?
1407  if (idx == index_set->ranges[range_idx].end)
1408  {
1409  // point to first element in next interval if possible
1410  if (range_idx < index_set->ranges.size() - 1)
1411  {
1412  ++range_idx;
1413  idx = index_set->ranges[range_idx].begin;
1414  }
1415  else
1416  {
1417  // we just fell off the end, set to invalid:
1418  range_idx = numbers::invalid_dof_index;
1420  }
1421  }
1422 }
1423 
1424 
1425 
1428 {
1429  advance();
1430  return *this;
1431 }
1432 
1433 
1434 
1437 {
1438  const IndexSet::ElementIterator it = *this;
1439  advance();
1440  return it;
1441 }
1442 
1443 
1444 
1445 inline bool
1447  const IndexSet::ElementIterator &other) const
1448 {
1449  return !(*this == other);
1450 }
1451 
1452 
1453 
1454 inline bool
1456  const IndexSet::ElementIterator &other) const
1457 {
1458  Assert(index_set == other.index_set,
1459  ExcMessage(
1460  "Can not compare iterators belonging to different IndexSets"));
1461  return range_idx < other.range_idx ||
1462  (range_idx == other.range_idx && idx < other.idx);
1463 }
1464 
1465 
1466 
1467 inline std::ptrdiff_t
1469  const IndexSet::ElementIterator &other) const
1470 {
1471  Assert(index_set == other.index_set,
1472  ExcMessage(
1473  "Can not compare iterators belonging to different IndexSets"));
1474  if (*this == other)
1475  return 0;
1476  if (!(*this < other))
1477  return -(other - *this);
1478 
1479  // only other can be equal to end() because of the checks above.
1480  Assert(is_valid(), ExcInternalError());
1481 
1482  // Note: we now compute how far advance *this in "*this < other" to get other,
1483  // so we need to return -c at the end.
1484 
1485  // first finish the current range:
1486  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1487 
1488  // now walk in steps of ranges (need to start one behind our current one):
1489  for (size_type range = range_idx + 1;
1490  range < index_set->ranges.size() && range <= other.range_idx;
1491  ++range)
1492  c += index_set->ranges[range].end - index_set->ranges[range].begin;
1493 
1494  Assert(
1495  other.range_idx < index_set->ranges.size() ||
1497  ExcMessage(
1498  "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1499 
1500  // We might have walked too far because we went until the end of
1501  // other.range_idx, so walk backwards to other.idx:
1503  c -= index_set->ranges[other.range_idx].end - other.idx;
1504 
1505  return -c;
1506 }
1507 
1508 
1509 /* Range */
1510 
1515 {}
1516 
1517 
1518 
1519 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1520  : begin(i1)
1521  , end(i2)
1523 {}
1524 
1525 
1526 
1527 /* IndexSet itself */
1528 
1530  : is_compressed(true)
1531  , index_space_size(0)
1533 {}
1534 
1535 
1536 
1537 inline IndexSet::IndexSet(const size_type size)
1538  : is_compressed(true)
1539  , index_space_size(size)
1540  , largest_range(numbers::invalid_unsigned_int)
1541 {}
1542 
1543 
1544 
1545 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1546  : ranges(std::move(is.ranges))
1547  , is_compressed(is.is_compressed)
1548  , index_space_size(is.index_space_size)
1549  , largest_range(is.largest_range)
1550 {
1551  is.ranges.clear();
1552  is.is_compressed = true;
1553  is.index_space_size = 0;
1554  is.largest_range = numbers::invalid_unsigned_int;
1555 
1556  compress();
1557 }
1558 
1559 
1560 
1561 inline IndexSet &
1563 {
1564  ranges = std::move(is.ranges);
1565  is_compressed = is.is_compressed;
1566  index_space_size = is.index_space_size;
1567  largest_range = is.largest_range;
1568 
1569  is.ranges.clear();
1570  is.is_compressed = true;
1571  is.index_space_size = 0;
1572  is.largest_range = numbers::invalid_unsigned_int;
1573 
1574  compress();
1575 
1576  return *this;
1577 }
1578 
1579 
1580 
1583 {
1584  compress();
1585  if (ranges.size() > 0)
1586  return {this, 0, ranges[0].begin};
1587  else
1588  return end();
1589 }
1590 
1591 
1592 
1595 {
1596  compress();
1597  return IndexSet::ElementIterator(this);
1598 }
1599 
1600 
1601 
1604 {
1605  compress();
1606  if (ranges.size() > 0)
1607  return IndexSet::IntervalIterator(this, 0);
1608  else
1609  return end_intervals();
1610 }
1611 
1612 
1613 
1616 {
1617  compress();
1618  return IndexSet::IntervalIterator(this);
1619 }
1620 
1621 
1622 
1623 inline void
1625 {
1626  // reset so that there are no indices in the set any more; however,
1627  // as documented, the index set retains its size
1628  ranges.clear();
1629  is_compressed = true;
1631 }
1632 
1633 
1634 
1635 inline void
1637 {
1638  Assert(ranges.empty(),
1639  ExcMessage("This function can only be called if the current "
1640  "object does not yet contain any elements."));
1641  index_space_size = sz;
1642  is_compressed = true;
1643 }
1644 
1645 
1646 
1647 inline IndexSet::size_type
1649 {
1650  return index_space_size;
1651 }
1652 
1653 
1654 
1655 inline void
1657 {
1658  if (is_compressed == true)
1659  return;
1660 
1661  do_compress();
1662 }
1663 
1664 
1665 
1666 inline void
1668 {
1670 
1671  const Range new_range(index, index + 1);
1672  if (ranges.size() == 0 || index > ranges.back().end)
1673  ranges.push_back(new_range);
1674  else if (index == ranges.back().end)
1675  ranges.back().end++;
1676  else
1677  add_range_lower_bound(new_range);
1678  is_compressed = false;
1679 }
1680 
1681 
1682 
1683 inline void
1685 {
1687  ((begin == index_space_size) && (end == index_space_size)),
1688  ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1690  ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1691  AssertIndexRange(begin, end + 1);
1692 
1693  if (begin != end)
1694  {
1695  const Range new_range(begin, end);
1696 
1697  // the new index might be larger than the last index present in the
1698  // ranges. Then we can skip the binary search
1699  if (ranges.size() == 0 || begin > ranges.back().end)
1700  ranges.push_back(new_range);
1701  else if (begin == ranges.back().end)
1702  ranges.back().end = end;
1703  else
1704  add_range_lower_bound(new_range);
1705 
1706  is_compressed = false;
1707  }
1708 }
1709 
1710 
1711 
1712 template <typename ForwardIterator>
1713 inline void
1714 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1715 {
1716  if (begin == end)
1717  return;
1718 
1719  // identify ranges in the given iterator range by checking whether some
1720  // indices happen to be consecutive. to avoid quadratic complexity when
1721  // calling add_range many times (as add_range() going into the middle of an
1722  // already existing range must shift entries around), we first collect a
1723  // vector of ranges.
1724  boost::container::small_vector<std::pair<size_type, size_type>, 200>
1725  tmp_ranges;
1726  bool ranges_are_sorted = true;
1727  for (ForwardIterator p = begin; p != end;)
1728  {
1729  const size_type begin_index = *p;
1730  size_type end_index = begin_index + 1;
1731  ForwardIterator q = p;
1732  ++q;
1733  while ((q != end) && (static_cast<size_type>(*q) == end_index))
1734  {
1735  ++end_index;
1736  ++q;
1737  }
1738 
1739  tmp_ranges.emplace_back(begin_index, end_index);
1740  p = q;
1741 
1742  // if the starting index of the next go-around of the for loop is less
1743  // than the end index of the one just identified, then we will have at
1744  // least one pair of ranges that are not sorted, and consequently the
1745  // whole collection of ranges is not sorted.
1746  if (p != end && static_cast<size_type>(*p) < end_index)
1747  ranges_are_sorted = false;
1748  }
1749 
1750  add_ranges_internal(tmp_ranges, ranges_are_sorted);
1751 }
1752 
1753 
1754 
1755 inline bool
1757 {
1758  if (ranges.empty() == false)
1759  {
1760  compress();
1761 
1762  // fast check whether the index is in the largest range
1764  if (index >= ranges[largest_range].begin &&
1765  index < ranges[largest_range].end)
1766  return true;
1767  else if (ranges.size() > 1)
1768  return is_element_binary_search(index);
1769  else
1770  return false;
1771  }
1772  else
1773  return false;
1774 }
1775 
1776 
1777 
1778 inline bool
1780 {
1781  compress();
1782  return (ranges.size() <= 1);
1783 }
1784 
1785 
1786 
1787 inline bool
1789 {
1790  return ranges.empty();
1791 }
1792 
1793 
1794 
1795 inline IndexSet::size_type
1797 {
1798  // make sure we have non-overlapping ranges
1799  compress();
1800 
1801  size_type v = 0;
1802  if (!ranges.empty())
1803  {
1804  Range &r = ranges.back();
1805  v = r.nth_index_in_set + r.end - r.begin;
1806  }
1807 
1808 #ifdef DEBUG
1809  size_type s = 0;
1810  for (const auto &range : ranges)
1811  s += (range.end - range.begin);
1812  Assert(s == v, ExcInternalError());
1813 #endif
1814 
1815  return v;
1816 }
1817 
1818 
1819 
1820 inline unsigned int
1822 {
1823  compress();
1824  return ranges.size();
1825 }
1826 
1827 
1828 
1829 inline IndexSet::size_type
1831 {
1832  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1833 
1834  compress();
1835  const std::vector<Range>::const_iterator main_range =
1836  ranges.begin() + largest_range;
1837 
1838  return main_range->nth_index_in_set;
1839 }
1840 
1841 
1842 
1843 inline IndexSet::size_type
1845 {
1847 
1848  compress();
1849 
1850  // first check whether the index is in the largest range
1852  const auto main_range = ranges.begin() + largest_range;
1853  if (n >= main_range->nth_index_in_set &&
1854  n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1855  return main_range->begin + (n - main_range->nth_index_in_set);
1856  else
1858 }
1859 
1860 
1861 
1862 inline IndexSet::size_type
1864 {
1865  // to make this call thread-safe, compress() must not be called through this
1866  // function
1867  Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1868  AssertIndexRange(n, size());
1869 
1870  // return immediately if the index set is empty
1871  if (is_empty())
1873 
1874  // check whether the index is in the largest range. use the result to
1875  // perform a one-sided binary search afterward
1877  if (n >= ranges[largest_range].begin && n < ranges[largest_range].end)
1878  return (n - ranges[largest_range].begin) +
1879  ranges[largest_range].nth_index_in_set;
1880  else if (ranges.size() > 1)
1882  else
1884 }
1885 
1886 
1887 
1888 inline bool
1890 {
1891  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1892 
1893  compress();
1894  is.compress();
1895 
1896  return ranges == is.ranges;
1897 }
1898 
1899 
1900 
1901 inline bool
1903 {
1904  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1905 
1906  compress();
1907  is.compress();
1908 
1909  return ranges != is.ranges;
1910 }
1911 
1912 
1913 
1914 template <typename Vector>
1915 void
1917 {
1918  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
1919 
1920  compress();
1921  // first fill all elements of the vector with zeroes.
1922  std::fill(vector.begin(), vector.end(), 0);
1923 
1924  // then write ones into the elements whose indices are contained in the
1925  // index set
1926  for (const auto &range : ranges)
1927  for (size_type i = range.begin; i < range.end; ++i)
1928  vector[i] = 1;
1929 }
1930 
1931 
1932 
1933 template <class StreamType>
1934 inline void
1935 IndexSet::print(StreamType &out) const
1936 {
1937  compress();
1938  out << '{';
1939  std::vector<Range>::const_iterator p;
1940  for (p = ranges.begin(); p != ranges.end(); ++p)
1941  {
1942  if (p->end - p->begin == 1)
1943  out << p->begin;
1944  else
1945  out << '[' << p->begin << ',' << p->end - 1 << ']';
1946 
1947  if (p != --ranges.end())
1948  out << ", ";
1949  }
1950  out << '}' << std::endl;
1951 }
1952 
1953 
1954 
1955 template <class Archive>
1956 inline void
1957 IndexSet::Range::serialize(Archive &ar, const unsigned int)
1958 {
1959  ar &begin &end &nth_index_in_set;
1960 }
1961 
1962 
1963 
1964 template <class Archive>
1965 inline void
1966 IndexSet::serialize(Archive &ar, const unsigned int)
1967 {
1969 }
1970 
1972 
1973 #endif
size_type operator*() const
Definition: index_set.h:1374
bool operator==(const ElementIterator &) const
Definition: index_set.h:1386
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition: index_set.h:1330
std::ptrdiff_t difference_type
Definition: index_set.h:811
ElementIterator & operator++()
Definition: index_set.h:1427
bool operator<(const ElementIterator &) const
Definition: index_set.h:1455
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1446
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1468
std::forward_iterator_tag iterator_category
Definition: index_set.h:809
const IndexSet * index_set
Definition: index_set.h:825
size_type last() const
Definition: index_set.h:1160
ElementIterator end() const
Definition: index_set.h:1146
size_type n_elements() const
Definition: index_set.h:1120
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1182
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1170
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1194
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:1089
const IndexSet * index_set
Definition: index_set.h:616
ElementIterator begin() const
Definition: index_set.h:1137
std::ptrdiff_t difference_type
Definition: index_set.h:719
const IntervalAccessor & operator*() const
Definition: index_set.h:1262
IntervalIterator(const IntervalIterator &other)=default
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1278
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1305
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1296
std::forward_iterator_tag iterator_category
Definition: index_set.h:717
IntervalIterator & operator++()
Definition: index_set.h:1243
IntervalIterator & operator=(const IntervalIterator &other)=default
const IntervalAccessor * operator->() const
Definition: index_set.h:1270
IntervalAccessor accessor
Definition: index_set.h:727
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1287
bool is_element_binary_search(const size_type local_index) const
Definition: index_set.cc:565
size_type index_within_set_binary_search(const size_type global_index) const
Definition: index_set.cc:618
size_type largest_range
Definition: index_set.h:997
bool is_contiguous() const
Definition: index_set.h:1779
unsigned int n_intervals() const
Definition: index_set.h:1821
IndexSet(const IndexSet &)=default
IndexSet & operator=(const IndexSet &)=default
IntervalIterator end_intervals() const
Definition: index_set.h:1615
void do_compress() const
Definition: index_set.cc:100
ElementIterator at(const size_type global_index) const
Definition: index_set.cc:640
size_type size() const
Definition: index_set.h:1648
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:709
bool is_empty() const
Definition: index_set.h:1788
IS make_petsc_is(const MPI_Comm &communicator=MPI_COMM_WORLD) const
Definition: index_set.cc:858
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1863
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
Definition: index_set.cc:241
size_type n_elements() const
Definition: index_set.h:1796
bool operator==(const IndexSet &is) const
Definition: index_set.h:1889
void add_range_lower_bound(const Range &range)
Definition: index_set.cc:389
bool is_element(const size_type index) const
Definition: index_set.h:1756
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1966
ElementIterator begin() const
Definition: index_set.h:1582
size_type pop_front()
Definition: index_set.cc:367
signed int value_type
Definition: index_set.h:96
void set_size(const size_type size)
Definition: index_set.h:1636
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1902
void read(std::istream &in)
Definition: index_set.cc:509
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:720
void clear()
Definition: index_set.h:1624
size_type largest_range_starting_index() const
Definition: index_set.h:1830
IndexSet tensor_product(const IndexSet &other) const
Definition: index_set.cc:337
void add_index(const size_type index)
Definition: index_set.h:1667
void write(std::ostream &out) const
Definition: index_set.cc:494
void block_read(std::istream &in)
Definition: index_set.cc:545
bool is_compressed
Definition: index_set.h:980
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:398
IntervalIterator begin_intervals() const
Definition: index_set.h:1603
void fill_binary_vector(VectorType &vector) const
std::vector< Range > ranges
Definition: index_set.h:970
void subtract_set(const IndexSet &other)
Definition: index_set.cc:268
ElementIterator end() const
Definition: index_set.h:1594
Threads::Mutex compress_mutex
Definition: index_set.h:1003
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1076
size_type index_space_size
Definition: index_set.h:986
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:789
void block_write(std::ostream &out) const
Definition: index_set.cc:531
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:214
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1684
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1844
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:878
std::size_t memory_consumption() const
Definition: index_set.cc:944
void print(StreamType &out) const
Definition: index_set.h:1935
size_type nth_index_in_set_binary_search(const size_type local_index) const
Definition: index_set.cc:601
void compress() const
Definition: index_set.h:1656
size_type pop_back()
Definition: index_set.cc:349
std::vector< size_type > get_index_vector() const
Definition: index_set.cc:690
types::global_dof_index size_type
Definition: index_set.h:77
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1714
IndexSet operator&(const IndexSet &is) const
Definition: vector.h:109
iterator end()
size_type size() const
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
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:1586
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:510
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:928
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:934
friend bool operator==(const Range &range_1, const Range &range_2)
Definition: index_set.h:941
size_type end
Definition: index_set.h:896
size_type nth_index_in_set
Definition: index_set.h:898
static std::size_t memory_consumption()
Definition: index_set.h:947
size_type begin
Definition: index_set.h:895
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1957
friend bool operator<(const Range &range_1, const Range &range_2)
Definition: index_set.h:920
void advance(std::tuple< I1, I2 > &t, const unsigned int n)