Reference documentation for deal.II version GIT 194dd8bb02 2022-12-03 08:20: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 - 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 
26 #include <boost/container/small_vector.hpp>
28 
29 #include <algorithm>
30 #include <vector>
31 
32 
33 #ifdef DEAL_II_WITH_TRILINOS
34 # include <Epetra_Map.h>
35 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
36 # include <Tpetra_Map.hpp>
37 # endif
38 #endif
39 
41 
64 class IndexSet
65 {
66 public:
67  // forward declarations:
68  class ElementIterator;
69  class IntervalIterator;
70 
76 
94  using value_type = signed int;
95 
96 
100  IndexSet();
101 
105  explicit IndexSet(const size_type size);
106 
110  IndexSet(const IndexSet &) = default;
111 
115  IndexSet &
116  operator=(const IndexSet &) = default;
117 
122  IndexSet(IndexSet &&is) noexcept;
123 
128  IndexSet &
129  operator=(IndexSet &&is) noexcept;
130 
131 #ifdef DEAL_II_WITH_TRILINOS
135  explicit IndexSet(const Epetra_BlockMap &map);
136 #endif
137 
142  void
143  clear();
144 
151  void
152  set_size(const size_type size);
153 
161  size_type
162  size() const;
163 
170  void
171  add_range(const size_type begin, const size_type end);
172 
176  void
177  add_index(const size_type index);
178 
188  template <typename ForwardIterator>
189  void
190  add_indices(const ForwardIterator &begin, const ForwardIterator &end);
191 
207  void
208  add_indices(const IndexSet &other, const size_type offset = 0);
209 
213  bool
214  is_element(const size_type index) const;
215 
220  bool
221  is_contiguous() const;
222 
227  bool
228  is_empty() const;
229 
239  bool
240  is_ascending_and_one_to_one(const MPI_Comm &communicator) const;
241 
245  size_type
246  n_elements() const;
247 
253  size_type
254  nth_index_in_set(const size_type local_index) const;
255 
262  size_type
264 
273  unsigned int
274  n_intervals() const;
275 
287  size_type
289 
294  void
295  compress() const;
296 
302  bool
303  operator==(const IndexSet &is) const;
304 
310  bool
311  operator!=(const IndexSet &is) const;
312 
319  IndexSet
320  operator&(const IndexSet &is) const;
321 
335  IndexSet
336  get_view(const size_type begin, const size_type end) const;
337 
343  std::vector<IndexSet>
345  const std::vector<types::global_dof_index> &n_indices_per_block) const;
346 
352  void
353  subtract_set(const IndexSet &other);
354 
371  IndexSet
372  tensor_product(const IndexSet &other) const;
373 
378  size_type
379  pop_back();
380 
385  size_type
386  pop_front();
387 
391  void
392  fill_index_vector(std::vector<size_type> &indices) const;
393 
406  template <typename VectorType>
407  void
408  fill_binary_vector(VectorType &vector) const;
409 
414  template <class StreamType>
415  void
416  print(StreamType &out) const;
417 
422  void
423  write(std::ostream &out) const;
424 
429  void
430  read(std::istream &in);
431 
436  void
437  block_write(std::ostream &out) const;
438 
443  void
444  block_read(std::istream &in);
445 
446 #ifdef DEAL_II_WITH_TRILINOS
475  Epetra_Map
476  make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
477  const bool overlapping = false) const;
478 
479 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
480  Tpetra::Map<int, types::global_dof_index>
481  make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
482  const bool overlapping = false) const;
483 # endif
484 #endif
485 
486 
491  std::size_t
492  memory_consumption() const;
493 
495  size_type,
496  << "The global index " << arg1
497  << " is not an element of this set.");
498 
504  template <class Archive>
505  void
506  serialize(Archive &ar, const unsigned int version);
507 
508 
520  {
521  public:
526  IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
527 
531  explicit IntervalAccessor(const IndexSet *idxset);
532 
536  size_type
537  n_elements() const;
538 
542  bool
543  is_valid() const;
544 
549  begin() const;
550 
556  end() const;
557 
561  size_type
562  last() const;
563 
564  private:
568  IntervalAccessor(const IntervalAccessor &other);
573  operator=(const IntervalAccessor &other);
574 
578  bool
579  operator==(const IntervalAccessor &other) const;
583  bool
584  operator<(const IntervalAccessor &other) const;
589  void
590  advance();
595 
601 
602  friend class IntervalIterator;
603  };
604 
610  {
611  public:
616  IntervalIterator(const IndexSet *idxset, const size_type range_idx);
617 
621  explicit IntervalIterator(const IndexSet *idxset);
622 
627 
631  IntervalIterator(const IntervalIterator &other) = default;
632 
637  operator=(const IntervalIterator &other) = default;
638 
643  operator++();
644 
649  operator++(int);
650 
654  const IntervalAccessor &
655  operator*() const;
656 
660  const IntervalAccessor *
661  operator->() const;
662 
666  bool
667  operator==(const IntervalIterator &) const;
668 
672  bool
673  operator!=(const IntervalIterator &) const;
674 
678  bool
679  operator<(const IntervalIterator &) const;
680 
687  int
688  operator-(const IntervalIterator &p) const;
689 
695  using iterator_category = std::forward_iterator_tag;
697  using difference_type = std::ptrdiff_t;
700 
701  private:
706  };
707 
713  {
714  public:
719  ElementIterator(const IndexSet *idxset,
720  const size_type range_idx,
721  const size_type index);
722 
726  explicit ElementIterator(const IndexSet *idxset);
727 
732  size_type
733  operator*() const;
734 
738  bool
739  is_valid() const;
740 
745  operator++();
746 
751  operator++(int);
752 
756  bool
757  operator==(const ElementIterator &) const;
758 
762  bool
763  operator!=(const ElementIterator &) const;
764 
768  bool
769  operator<(const ElementIterator &) const;
770 
779  std::ptrdiff_t
780  operator-(const ElementIterator &p) const;
781 
787  using iterator_category = std::forward_iterator_tag;
789  using difference_type = std::ptrdiff_t;
790  using pointer = size_type *;
791  using reference = size_type &;
792 
793  private:
797  void
798  advance();
799 
812  };
813 
819  begin() const;
820 
836  at(const size_type global_index) const;
837 
843  end() const;
844 
849  begin_intervals() const;
850 
856  end_intervals() const;
857 
862 private:
871  struct Range
872  {
875 
877 
886  Range();
887 
895  Range(const size_type i1, const size_type i2);
896 
897  friend inline bool
898  operator<(const Range &range_1, const Range &range_2)
899  {
900  return (
901  (range_1.begin < range_2.begin) ||
902  ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
903  }
904 
905  static bool
907  {
908  return x.end < y.end;
909  }
910 
911  static bool
913  {
914  return (x.nth_index_in_set + (x.end - x.begin) <
915  y.nth_index_in_set + (y.end - y.begin));
916  }
917 
918  friend inline bool
919  operator==(const Range &range_1, const Range &range_2)
920  {
921  return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
922  }
923 
924  static std::size_t
926  {
927  return sizeof(Range);
928  }
929 
935  template <class Archive>
936  void
937  serialize(Archive &ar, const unsigned int version);
938  };
939 
948  mutable std::vector<Range> ranges;
949 
958  mutable bool is_compressed;
959 
965 
976 
982 
986  void
987  do_compress() const;
988 
995  bool
996  is_element_binary_search(const size_type local_index) const;
997 
1004  size_type
1005  nth_index_in_set_binary_search(const size_type local_index) const;
1006 
1013  size_type
1015 
1021  void
1022  add_range_lower_bound(const Range &range);
1023 
1027  void
1029  boost::container::small_vector<std::pair<size_type, size_type>, 200>
1030  & tmp_ranges,
1031  const bool ranges_are_sorted);
1032 };
1033 
1034 
1035 
1053 inline IndexSet
1055 {
1056  IndexSet is(N);
1057  is.add_range(0, N);
1058  is.compress();
1059  return is;
1060 }
1061 
1062 /* ------------------ inline functions ------------------ */
1063 
1064 
1065 /* IntervalAccessor */
1066 
1068  const IndexSet * idxset,
1069  const IndexSet::size_type range_idx)
1070  : index_set(idxset)
1071  , range_idx(range_idx)
1072 {
1073  Assert(range_idx < idxset->n_intervals(),
1074  ExcInternalError("Invalid range index"));
1075 }
1076 
1077 
1078 
1080  : index_set(idxset)
1081  , range_idx(numbers::invalid_dof_index)
1082 {}
1083 
1084 
1085 
1087  const IndexSet::IntervalAccessor &other)
1088  : index_set(other.index_set)
1089  , range_idx(other.range_idx)
1090 {
1092  ExcMessage("invalid iterator"));
1093 }
1094 
1095 
1096 
1097 inline IndexSet::size_type
1099 {
1100  Assert(is_valid(), ExcMessage("invalid iterator"));
1101  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1102 }
1103 
1104 
1105 
1106 inline bool
1108 {
1109  return index_set != nullptr && range_idx < index_set->n_intervals();
1110 }
1111 
1112 
1113 
1116 {
1117  Assert(is_valid(), ExcMessage("invalid iterator"));
1118  return {index_set, range_idx, index_set->ranges[range_idx].begin};
1119 }
1120 
1121 
1122 
1125 {
1126  Assert(is_valid(), ExcMessage("invalid iterator"));
1127 
1128  // point to first index in next interval unless we are the last interval.
1129  if (range_idx < index_set->ranges.size() - 1)
1130  return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1131  else
1132  return index_set->end();
1133 }
1134 
1135 
1136 
1137 inline IndexSet::size_type
1139 {
1140  Assert(is_valid(), ExcMessage("invalid iterator"));
1141 
1142  return index_set->ranges[range_idx].end - 1;
1143 }
1144 
1145 
1146 
1149 {
1150  index_set = other.index_set;
1151  range_idx = other.range_idx;
1152  Assert(range_idx == numbers::invalid_dof_index || is_valid(),
1153  ExcMessage("invalid iterator"));
1154  return *this;
1155 }
1156 
1157 
1158 
1159 inline bool
1161  const IndexSet::IntervalAccessor &other) const
1162 {
1163  Assert(index_set == other.index_set,
1164  ExcMessage(
1165  "Can not compare accessors pointing to different IndexSets"));
1166  return range_idx == other.range_idx;
1167 }
1168 
1169 
1170 
1171 inline bool
1173  const IndexSet::IntervalAccessor &other) const
1174 {
1175  Assert(index_set == other.index_set,
1176  ExcMessage(
1177  "Can not compare accessors pointing to different IndexSets"));
1178  return range_idx < other.range_idx;
1179 }
1180 
1181 
1182 
1183 inline void
1185 {
1186  Assert(
1187  is_valid(),
1188  ExcMessage(
1189  "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1190  ++range_idx;
1191 
1192  // set ourselves to invalid if we walk off the end
1193  if (range_idx >= index_set->ranges.size())
1194  range_idx = numbers::invalid_dof_index;
1195 }
1196 
1197 
1198 /* IntervalIterator */
1199 
1201  const IndexSet * idxset,
1202  const IndexSet::size_type range_idx)
1203  : accessor(idxset, range_idx)
1204 {}
1205 
1206 
1207 
1209  : accessor(nullptr)
1210 {}
1211 
1212 
1213 
1215  : accessor(idxset)
1216 {}
1217 
1218 
1219 
1222 {
1223  accessor.advance();
1224  return *this;
1225 }
1226 
1227 
1228 
1231 {
1232  const IndexSet::IntervalIterator iter = *this;
1233  accessor.advance();
1234  return iter;
1235 }
1236 
1237 
1238 
1239 inline const IndexSet::IntervalAccessor &
1241 {
1242  return accessor;
1243 }
1244 
1245 
1246 
1247 inline const IndexSet::IntervalAccessor *
1249 {
1250  return &accessor;
1251 }
1252 
1253 
1254 
1255 inline bool
1257  const IndexSet::IntervalIterator &other) const
1258 {
1259  return accessor == other.accessor;
1260 }
1261 
1262 
1263 
1264 inline bool
1266  const IndexSet::IntervalIterator &other) const
1267 {
1268  return !(*this == other);
1269 }
1270 
1271 
1272 
1273 inline bool
1275  const IndexSet::IntervalIterator &other) const
1276 {
1277  return accessor < other.accessor;
1278 }
1279 
1280 
1281 
1282 inline int
1284  const IndexSet::IntervalIterator &other) const
1285 {
1286  Assert(accessor.index_set == other.accessor.index_set,
1287  ExcMessage(
1288  "Can not compare iterators belonging to different IndexSets"));
1289 
1290  const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index) ?
1291  accessor.index_set->ranges.size() :
1292  accessor.range_idx;
1293  const size_type rhs =
1295  accessor.index_set->ranges.size() :
1296  other.accessor.range_idx;
1297 
1298  if (lhs > rhs)
1299  return static_cast<int>(lhs - rhs);
1300  else
1301  return -static_cast<int>(rhs - lhs);
1302 }
1303 
1304 
1305 
1306 /* ElementIterator */
1307 
1309  const IndexSet * idxset,
1310  const IndexSet::size_type range_idx,
1311  const IndexSet::size_type index)
1312  : index_set(idxset)
1313  , range_idx(range_idx)
1314  , idx(index)
1315 {
1316  Assert(range_idx < index_set->ranges.size(),
1317  ExcMessage(
1318  "Invalid range index for IndexSet::ElementIterator constructor."));
1319  Assert(
1320  idx >= index_set->ranges[range_idx].begin &&
1321  idx < index_set->ranges[range_idx].end,
1323  "Invalid index argument for IndexSet::ElementIterator constructor."));
1324 }
1325 
1326 
1327 
1329  : index_set(idxset)
1330  , range_idx(numbers::invalid_dof_index)
1331  , idx(numbers::invalid_dof_index)
1332 {}
1333 
1334 
1335 
1336 inline bool
1338 {
1339  Assert((range_idx == numbers::invalid_dof_index &&
1340  idx == numbers::invalid_dof_index) ||
1341  (range_idx < index_set->ranges.size() &&
1342  idx < index_set->ranges[range_idx].end),
1343  ExcInternalError("Invalid ElementIterator state."));
1344 
1345  return (range_idx < index_set->ranges.size() &&
1346  idx < index_set->ranges[range_idx].end);
1347 }
1348 
1349 
1350 
1351 inline IndexSet::size_type
1353 {
1354  Assert(
1355  is_valid(),
1356  ExcMessage(
1357  "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1358  return idx;
1359 }
1360 
1361 
1362 
1363 inline bool
1365  const IndexSet::ElementIterator &other) const
1366 {
1367  Assert(index_set == other.index_set,
1368  ExcMessage(
1369  "Can not compare iterators belonging to different IndexSets"));
1370  return range_idx == other.range_idx && idx == other.idx;
1371 }
1372 
1373 
1374 
1375 inline void
1377 {
1378  Assert(
1379  is_valid(),
1380  ExcMessage(
1381  "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1382  if (idx < index_set->ranges[range_idx].end)
1383  ++idx;
1384  // end of this range?
1385  if (idx == index_set->ranges[range_idx].end)
1386  {
1387  // point to first element in next interval if possible
1388  if (range_idx < index_set->ranges.size() - 1)
1389  {
1390  ++range_idx;
1391  idx = index_set->ranges[range_idx].begin;
1392  }
1393  else
1394  {
1395  // we just fell off the end, set to invalid:
1396  range_idx = numbers::invalid_dof_index;
1398  }
1399  }
1400 }
1401 
1402 
1403 
1406 {
1407  advance();
1408  return *this;
1409 }
1410 
1411 
1412 
1415 {
1416  const IndexSet::ElementIterator it = *this;
1417  advance();
1418  return it;
1419 }
1420 
1421 
1422 
1423 inline bool
1425  const IndexSet::ElementIterator &other) const
1426 {
1427  return !(*this == other);
1428 }
1429 
1430 
1431 
1432 inline bool
1434  const IndexSet::ElementIterator &other) const
1435 {
1436  Assert(index_set == other.index_set,
1437  ExcMessage(
1438  "Can not compare iterators belonging to different IndexSets"));
1439  return range_idx < other.range_idx ||
1440  (range_idx == other.range_idx && idx < other.idx);
1441 }
1442 
1443 
1444 
1445 inline std::ptrdiff_t
1447  const IndexSet::ElementIterator &other) const
1448 {
1449  Assert(index_set == other.index_set,
1450  ExcMessage(
1451  "Can not compare iterators belonging to different IndexSets"));
1452  if (*this == other)
1453  return 0;
1454  if (!(*this < other))
1455  return -(other - *this);
1456 
1457  // only other can be equal to end() because of the checks above.
1458  Assert(is_valid(), ExcInternalError());
1459 
1460  // Note: we now compute how far advance *this in "*this < other" to get other,
1461  // so we need to return -c at the end.
1462 
1463  // first finish the current range:
1464  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1465 
1466  // now walk in steps of ranges (need to start one behind our current one):
1467  for (size_type range = range_idx + 1;
1468  range < index_set->ranges.size() && range <= other.range_idx;
1469  ++range)
1470  c += index_set->ranges[range].end - index_set->ranges[range].begin;
1471 
1472  Assert(
1473  other.range_idx < index_set->ranges.size() ||
1475  ExcMessage(
1476  "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1477 
1478  // We might have walked too far because we went until the end of
1479  // other.range_idx, so walk backwards to other.idx:
1481  c -= index_set->ranges[other.range_idx].end - other.idx;
1482 
1483  return -c;
1484 }
1485 
1486 
1487 /* Range */
1488 
1493 {}
1494 
1495 
1496 
1497 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1498  : begin(i1)
1499  , end(i2)
1501 {}
1502 
1503 
1504 
1505 /* IndexSet itself */
1506 
1508  : is_compressed(true)
1509  , index_space_size(0)
1511 {}
1512 
1513 
1514 
1515 inline IndexSet::IndexSet(const size_type size)
1516  : is_compressed(true)
1517  , index_space_size(size)
1518  , largest_range(numbers::invalid_unsigned_int)
1519 {}
1520 
1521 
1522 
1523 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1524  : ranges(std::move(is.ranges))
1525  , is_compressed(is.is_compressed)
1526  , index_space_size(is.index_space_size)
1527  , largest_range(is.largest_range)
1528 {
1529  is.ranges.clear();
1530  is.is_compressed = true;
1531  is.index_space_size = 0;
1532  is.largest_range = numbers::invalid_unsigned_int;
1533 
1534  compress();
1535 }
1536 
1537 
1538 
1539 inline IndexSet &
1541 {
1542  ranges = std::move(is.ranges);
1543  is_compressed = is.is_compressed;
1544  index_space_size = is.index_space_size;
1545  largest_range = is.largest_range;
1546 
1547  is.ranges.clear();
1548  is.is_compressed = true;
1549  is.index_space_size = 0;
1550  is.largest_range = numbers::invalid_unsigned_int;
1551 
1552  compress();
1553 
1554  return *this;
1555 }
1556 
1557 
1558 
1561 {
1562  compress();
1563  if (ranges.size() > 0)
1564  return {this, 0, ranges[0].begin};
1565  else
1566  return end();
1567 }
1568 
1569 
1570 
1573 {
1574  compress();
1575  return IndexSet::ElementIterator(this);
1576 }
1577 
1578 
1579 
1582 {
1583  compress();
1584  if (ranges.size() > 0)
1585  return IndexSet::IntervalIterator(this, 0);
1586  else
1587  return end_intervals();
1588 }
1589 
1590 
1591 
1594 {
1595  compress();
1596  return IndexSet::IntervalIterator(this);
1597 }
1598 
1599 
1600 
1601 inline void
1603 {
1604  // reset so that there are no indices in the set any more; however,
1605  // as documented, the index set retains its size
1606  ranges.clear();
1607  is_compressed = true;
1609 }
1610 
1611 
1612 
1613 inline void
1615 {
1616  Assert(ranges.empty(),
1617  ExcMessage("This function can only be called if the current "
1618  "object does not yet contain any elements."));
1619  index_space_size = sz;
1620  is_compressed = true;
1621 }
1622 
1623 
1624 
1625 inline IndexSet::size_type
1627 {
1628  return index_space_size;
1629 }
1630 
1631 
1632 
1633 inline void
1635 {
1636  if (is_compressed == true)
1637  return;
1638 
1639  do_compress();
1640 }
1641 
1642 
1643 
1644 inline void
1646 {
1648 
1649  const Range new_range(index, index + 1);
1650  if (ranges.size() == 0 || index > ranges.back().end)
1651  ranges.push_back(new_range);
1652  else if (index == ranges.back().end)
1653  ranges.back().end++;
1654  else
1655  add_range_lower_bound(new_range);
1656  is_compressed = false;
1657 }
1658 
1659 
1660 
1661 inline void
1663 {
1665  ((begin == index_space_size) && (end == index_space_size)),
1666  ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1668  ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1669  AssertIndexRange(begin, end + 1);
1670 
1671  if (begin != end)
1672  {
1673  const Range new_range(begin, end);
1674 
1675  // the new index might be larger than the last index present in the
1676  // ranges. Then we can skip the binary search
1677  if (ranges.size() == 0 || begin > ranges.back().end)
1678  ranges.push_back(new_range);
1679  else if (begin == ranges.back().end)
1680  ranges.back().end = end;
1681  else
1682  add_range_lower_bound(new_range);
1683 
1684  is_compressed = false;
1685  }
1686 }
1687 
1688 
1689 
1690 template <typename ForwardIterator>
1691 inline void
1692 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1693 {
1694  if (begin == end)
1695  return;
1696 
1697  // identify ranges in the given iterator range by checking whether some
1698  // indices happen to be consecutive. to avoid quadratic complexity when
1699  // calling add_range many times (as add_range() going into the middle of an
1700  // already existing range must shift entries around), we first collect a
1701  // vector of ranges.
1702  boost::container::small_vector<std::pair<size_type, size_type>, 200>
1703  tmp_ranges;
1704  bool ranges_are_sorted = true;
1705  for (ForwardIterator p = begin; p != end;)
1706  {
1707  const size_type begin_index = *p;
1708  size_type end_index = begin_index + 1;
1709  ForwardIterator q = p;
1710  ++q;
1711  while ((q != end) && (static_cast<size_type>(*q) == end_index))
1712  {
1713  ++end_index;
1714  ++q;
1715  }
1716 
1717  tmp_ranges.emplace_back(begin_index, end_index);
1718  p = q;
1719 
1720  // if the starting index of the next go-around of the for loop is less
1721  // than the end index of the one just identified, then we will have at
1722  // least one pair of ranges that are not sorted, and consequently the
1723  // whole collection of ranges is not sorted.
1724  if (p != end && static_cast<size_type>(*p) < end_index)
1725  ranges_are_sorted = false;
1726  }
1727 
1728  add_ranges_internal(tmp_ranges, ranges_are_sorted);
1729 }
1730 
1731 
1732 
1733 inline bool
1735 {
1736  if (ranges.empty() == false)
1737  {
1738  compress();
1739 
1740  // fast check whether the index is in the largest range
1742  if (index >= ranges[largest_range].begin &&
1743  index < ranges[largest_range].end)
1744  return true;
1745  else if (ranges.size() > 1)
1746  return is_element_binary_search(index);
1747  else
1748  return false;
1749  }
1750  else
1751  return false;
1752 }
1753 
1754 
1755 
1756 inline bool
1758 {
1759  compress();
1760  return (ranges.size() <= 1);
1761 }
1762 
1763 
1764 
1765 inline bool
1767 {
1768  return ranges.empty();
1769 }
1770 
1771 
1772 
1773 inline IndexSet::size_type
1775 {
1776  // make sure we have non-overlapping ranges
1777  compress();
1778 
1779  size_type v = 0;
1780  if (!ranges.empty())
1781  {
1782  Range &r = ranges.back();
1783  v = r.nth_index_in_set + r.end - r.begin;
1784  }
1785 
1786 #ifdef DEBUG
1787  size_type s = 0;
1788  for (const auto &range : ranges)
1789  s += (range.end - range.begin);
1790  Assert(s == v, ExcInternalError());
1791 #endif
1792 
1793  return v;
1794 }
1795 
1796 
1797 
1798 inline unsigned int
1800 {
1801  compress();
1802  return ranges.size();
1803 }
1804 
1805 
1806 
1807 inline IndexSet::size_type
1809 {
1810  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1811 
1812  compress();
1813  const std::vector<Range>::const_iterator main_range =
1814  ranges.begin() + largest_range;
1815 
1816  return main_range->nth_index_in_set;
1817 }
1818 
1819 
1820 
1821 inline IndexSet::size_type
1823 {
1825 
1826  compress();
1827 
1828  // first check whether the index is in the largest range
1830  const auto main_range = ranges.begin() + largest_range;
1831  if (n >= main_range->nth_index_in_set &&
1832  n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1833  return main_range->begin + (n - main_range->nth_index_in_set);
1834  else
1836 }
1837 
1838 
1839 
1840 inline IndexSet::size_type
1842 {
1843  // to make this call thread-safe, compress() must not be called through this
1844  // function
1845  Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1846  AssertIndexRange(n, size());
1847 
1848  // return immediately if the index set is empty
1849  if (is_empty())
1851 
1852  // check whether the index is in the largest range. use the result to
1853  // perform a one-sided binary search afterward
1855  if (n >= ranges[largest_range].begin && n < ranges[largest_range].end)
1856  return (n - ranges[largest_range].begin) +
1857  ranges[largest_range].nth_index_in_set;
1858  else if (ranges.size() > 1)
1860  else
1862 }
1863 
1864 
1865 
1866 inline bool
1868 {
1869  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1870 
1871  compress();
1872  is.compress();
1873 
1874  return ranges == is.ranges;
1875 }
1876 
1877 
1878 
1879 inline bool
1881 {
1882  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1883 
1884  compress();
1885  is.compress();
1886 
1887  return ranges != is.ranges;
1888 }
1889 
1890 
1891 
1892 template <typename Vector>
1893 void
1895 {
1896  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
1897 
1898  compress();
1899  // first fill all elements of the vector with zeroes.
1900  std::fill(vector.begin(), vector.end(), 0);
1901 
1902  // then write ones into the elements whose indices are contained in the
1903  // index set
1904  for (const auto &range : ranges)
1905  for (size_type i = range.begin; i < range.end; ++i)
1906  vector[i] = 1;
1907 }
1908 
1909 
1910 
1911 template <class StreamType>
1912 inline void
1913 IndexSet::print(StreamType &out) const
1914 {
1915  compress();
1916  out << '{';
1917  std::vector<Range>::const_iterator p;
1918  for (p = ranges.begin(); p != ranges.end(); ++p)
1919  {
1920  if (p->end - p->begin == 1)
1921  out << p->begin;
1922  else
1923  out << '[' << p->begin << ',' << p->end - 1 << ']';
1924 
1925  if (p != --ranges.end())
1926  out << ", ";
1927  }
1928  out << '}' << std::endl;
1929 }
1930 
1931 
1932 
1933 template <class Archive>
1934 inline void
1935 IndexSet::Range::serialize(Archive &ar, const unsigned int)
1936 {
1937  ar &begin &end &nth_index_in_set;
1938 }
1939 
1940 
1941 
1942 template <class Archive>
1943 inline void
1944 IndexSet::serialize(Archive &ar, const unsigned int)
1945 {
1947 }
1948 
1950 
1951 #endif
size_type operator*() const
Definition: index_set.h:1352
bool operator==(const ElementIterator &) const
Definition: index_set.h:1364
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition: index_set.h:1308
std::ptrdiff_t difference_type
Definition: index_set.h:789
ElementIterator & operator++()
Definition: index_set.h:1405
bool operator<(const ElementIterator &) const
Definition: index_set.h:1433
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1424
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1446
std::forward_iterator_tag iterator_category
Definition: index_set.h:787
const IndexSet * index_set
Definition: index_set.h:803
size_type last() const
Definition: index_set.h:1138
ElementIterator end() const
Definition: index_set.h:1124
size_type n_elements() const
Definition: index_set.h:1098
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1160
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1148
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1172
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:1067
const IndexSet * index_set
Definition: index_set.h:594
ElementIterator begin() const
Definition: index_set.h:1115
std::ptrdiff_t difference_type
Definition: index_set.h:697
const IntervalAccessor & operator*() const
Definition: index_set.h:1240
IntervalIterator(const IntervalIterator &other)=default
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1256
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1283
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1274
std::forward_iterator_tag iterator_category
Definition: index_set.h:695
IntervalIterator & operator++()
Definition: index_set.h:1221
IntervalIterator & operator=(const IntervalIterator &other)=default
const IntervalAccessor * operator->() const
Definition: index_set.h:1248
IntervalAccessor accessor
Definition: index_set.h:705
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1265
bool is_element_binary_search(const size_type local_index) const
Definition: index_set.cc:563
size_type index_within_set_binary_search(const size_type global_index) const
Definition: index_set.cc:616
size_type largest_range
Definition: index_set.h:975
bool is_contiguous() const
Definition: index_set.h:1757
unsigned int n_intervals() const
Definition: index_set.h:1799
IndexSet(const IndexSet &)=default
IndexSet & operator=(const IndexSet &)=default
IntervalIterator end_intervals() const
Definition: index_set.h:1593
void do_compress() const
Definition: index_set.cc:98
ElementIterator at(const size_type global_index) const
Definition: index_set.cc:638
size_type size() const
Definition: index_set.h:1626
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:688
bool is_empty() const
Definition: index_set.h:1766
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1841
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
Definition: index_set.cc:239
size_type n_elements() const
Definition: index_set.h:1774
bool operator==(const IndexSet &is) const
Definition: index_set.h:1867
void add_range_lower_bound(const Range &range)
Definition: index_set.cc:387
bool is_element(const size_type index) const
Definition: index_set.h:1734
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1944
ElementIterator begin() const
Definition: index_set.h:1560
size_type pop_front()
Definition: index_set.cc:365
signed int value_type
Definition: index_set.h:94
void set_size(const size_type size)
Definition: index_set.h:1614
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1880
void read(std::istream &in)
Definition: index_set.cc:507
void clear()
Definition: index_set.h:1602
size_type largest_range_starting_index() const
Definition: index_set.h:1808
Tpetra::Map< int, types::global_dof_index > make_tpetra_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:708
IndexSet tensor_product(const IndexSet &other) const
Definition: index_set.cc:335
void add_index(const size_type index)
Definition: index_set.h:1645
void write(std::ostream &out) const
Definition: index_set.cc:492
void block_read(std::istream &in)
Definition: index_set.cc:543
bool is_compressed
Definition: index_set.h:958
void add_ranges_internal(boost::container::small_vector< std::pair< size_type, size_type >, 200 > &tmp_ranges, const bool ranges_are_sorted)
Definition: index_set.cc:396
IntervalIterator begin_intervals() const
Definition: index_set.h:1581
void fill_binary_vector(VectorType &vector) const
std::vector< Range > ranges
Definition: index_set.h:948
void subtract_set(const IndexSet &other)
Definition: index_set.cc:266
ElementIterator end() const
Definition: index_set.h:1572
Threads::Mutex compress_mutex
Definition: index_set.h:981
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1054
size_type index_space_size
Definition: index_set.h:964
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:777
void block_write(std::ostream &out) const
Definition: index_set.cc:529
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:212
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1662
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1822
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:847
std::size_t memory_consumption() const
Definition: index_set.cc:913
void print(StreamType &out) const
Definition: index_set.h:1913
size_type nth_index_in_set_binary_search(const size_type local_index) const
Definition: index_set.cc:599
void compress() const
Definition: index_set.h:1634
size_type pop_back()
Definition: index_set.cc:347
types::global_dof_index size_type
Definition: index_set.h:75
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1692
IndexSet operator&(const IndexSet &is) const
Definition: vector.h:109
iterator end()
size_type size() const
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:512
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:1501
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
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:392
static const unsigned int invalid_unsigned_int
Definition: types.h:206
const types::global_dof_index invalid_dof_index
Definition: types.h:226
unsigned int global_dof_index
Definition: types.h:81
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:906
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:912
friend bool operator==(const Range &range_1, const Range &range_2)
Definition: index_set.h:919
size_type end
Definition: index_set.h:874
size_type nth_index_in_set
Definition: index_set.h:876
static std::size_t memory_consumption()
Definition: index_set.h:925
size_type begin
Definition: index_set.h:873
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1935
friend bool operator<(const Range &range_1, const Range &range_2)
Definition: index_set.h:898
void advance(std::tuple< I1, I2 > &t, const unsigned int n)