Reference documentation for deal.II version Git 5a2787e538 2021-09-21 14:55:10 -0600
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
index_set.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_index_set_h
17 #define dealii_index_set_h
18 
19 #include <deal.II/base/config.h>
20 
23 #include <deal.II/base/utilities.h>
24 
25 #include <boost/serialization/vector.hpp>
26 
27 #include <algorithm>
28 #include <iterator>
29 #include <vector>
30 
31 
32 #ifdef DEAL_II_WITH_TRILINOS
33 # include <Epetra_Map.h>
34 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
35 # include <Tpetra_Map.hpp>
36 # endif
37 #endif
38 
39 #if defined(DEAL_II_WITH_MPI) || defined(DEAL_II_WITH_PETSC)
40 # include <mpi.h>
41 #else
42 using MPI_Comm = int;
43 # ifndef MPI_COMM_WORLD
44 # define MPI_COMM_WORLD 0
45 # endif
46 #endif
47 
49 
72 class IndexSet
73 {
74 public:
75  // forward declarations:
76  class ElementIterator;
77  class IntervalIterator;
78 
84 
102  using value_type = signed int;
103 
104 
108  IndexSet();
109 
113  explicit IndexSet(const size_type size);
114 
118  IndexSet(const IndexSet &) = default;
119 
123  IndexSet &
124  operator=(const IndexSet &) = default;
125 
130  IndexSet(IndexSet &&is) noexcept;
131 
136  IndexSet &
137  operator=(IndexSet &&is) noexcept;
138 
139 #ifdef DEAL_II_WITH_TRILINOS
140 
143  explicit IndexSet(const Epetra_BlockMap &map);
144 #endif
145 
150  void
151  clear();
152 
159  void
160  set_size(const size_type size);
161 
169  size_type
170  size() const;
171 
178  void
179  add_range(const size_type begin, const size_type end);
180 
184  void
185  add_index(const size_type index);
186 
196  template <typename ForwardIterator>
197  void
198  add_indices(const ForwardIterator &begin, const ForwardIterator &end);
199 
215  void
216  add_indices(const IndexSet &other, const size_type offset = 0);
217 
221  bool
222  is_element(const size_type index) const;
223 
228  bool
229  is_contiguous() const;
230 
235  bool
236  is_empty() const;
237 
247  bool
248  is_ascending_and_one_to_one(const MPI_Comm &communicator) const;
249 
253  size_type
254  n_elements() const;
255 
261  size_type
262  nth_index_in_set(const size_type local_index) const;
263 
270  size_type
272 
281  unsigned int
282  n_intervals() const;
283 
295  size_type
297 
302  void
303  compress() const;
304 
310  bool
311  operator==(const IndexSet &is) const;
312 
318  bool
319  operator!=(const IndexSet &is) const;
320 
327  IndexSet
328  operator&(const IndexSet &is) const;
329 
343  IndexSet
344  get_view(const size_type begin, const size_type end) const;
345 
351  std::vector<IndexSet>
353  const std::vector<types::global_dof_index> &n_indices_per_block) const;
354 
360  void
361  subtract_set(const IndexSet &other);
362 
379  IndexSet
380  tensor_product(const IndexSet &other) const;
381 
386  size_type
387  pop_back();
388 
393  size_type
394  pop_front();
395 
399  void
400  fill_index_vector(std::vector<size_type> &indices) const;
401 
414  template <typename VectorType>
415  void
416  fill_binary_vector(VectorType &vector) const;
417 
422  template <class StreamType>
423  void
424  print(StreamType &out) const;
425 
430  void
431  write(std::ostream &out) const;
432 
437  void
438  read(std::istream &in);
439 
444  void
445  block_write(std::ostream &out) const;
446 
451  void
452  block_read(std::istream &in);
453 
454 #ifdef DEAL_II_WITH_TRILINOS
455 
483  Epetra_Map
484  make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
485  const bool overlapping = false) const;
486 
487 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
488  Tpetra::Map<int, types::global_dof_index>
489  make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
490  const bool overlapping = false) const;
491 # endif
492 #endif
493 
494 
499  std::size_t
500  memory_consumption() const;
501 
503  size_type,
504  << "The global index " << arg1
505  << " is not an element of this set.");
506 
512  template <class Archive>
513  void
514  serialize(Archive &ar, const unsigned int version);
515 
516 
528  {
529  public:
534  IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
535 
539  explicit IntervalAccessor(const IndexSet *idxset);
540 
544  size_type
545  n_elements() const;
546 
550  bool
551  is_valid() const;
552 
557  begin() const;
558 
564  end() const;
565 
569  size_type
570  last() const;
571 
572  private:
576  IntervalAccessor(const IntervalAccessor &other);
581  operator=(const IntervalAccessor &other);
582 
586  bool
587  operator==(const IntervalAccessor &other) const;
591  bool
592  operator<(const IntervalAccessor &other) const;
597  void
598  advance();
603 
609 
610  friend class IntervalIterator;
611  };
612 
618  {
619  public:
624  IntervalIterator(const IndexSet *idxset, const size_type range_idx);
625 
629  explicit IntervalIterator(const IndexSet *idxset);
630 
635 
639  IntervalIterator(const IntervalIterator &other) = default;
640 
645  operator=(const IntervalIterator &other) = default;
646 
651  operator++();
652 
657  operator++(int);
658 
662  const IntervalAccessor &
663  operator*() const;
664 
668  const IntervalAccessor *
669  operator->() const;
670 
674  bool
675  operator==(const IntervalIterator &) const;
676 
680  bool
681  operator!=(const IntervalIterator &) const;
682 
686  bool
687  operator<(const IntervalIterator &) const;
688 
695  int
696  operator-(const IntervalIterator &p) const;
697 
703  using iterator_category = std::forward_iterator_tag;
705  using difference_type = std::ptrdiff_t;
708 
709  private:
714  };
715 
721  {
722  public:
727  ElementIterator(const IndexSet *idxset,
728  const size_type range_idx,
729  const size_type index);
730 
734  explicit ElementIterator(const IndexSet *idxset);
735 
740  size_type
741  operator*() const;
742 
746  bool
747  is_valid() const;
748 
753  operator++();
754 
759  operator++(int);
760 
764  bool
765  operator==(const ElementIterator &) const;
766 
770  bool
771  operator!=(const ElementIterator &) const;
772 
776  bool
777  operator<(const ElementIterator &) const;
778 
787  std::ptrdiff_t
788  operator-(const ElementIterator &p) const;
789 
795  using iterator_category = std::forward_iterator_tag;
797  using difference_type = std::ptrdiff_t;
798  using pointer = size_type *;
799  using reference = size_type &;
800 
801  private:
805  void
806  advance();
807 
820  };
821 
827  begin() const;
828 
844  at(const size_type global_index) const;
845 
851  end() const;
852 
857  begin_intervals() const;
858 
864  end_intervals() const;
865 
870 private:
879  struct Range
880  {
883 
885 
894  Range();
895 
903  Range(const size_type i1, const size_type i2);
904 
905  friend inline bool
906  operator<(const Range &range_1, const Range &range_2)
907  {
908  return (
909  (range_1.begin < range_2.begin) ||
910  ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
911  }
912 
913  static bool
915  {
916  return x.end < y.end;
917  }
918 
919  static bool
921  {
922  return (x.nth_index_in_set + (x.end - x.begin) <
923  y.nth_index_in_set + (y.end - y.begin));
924  }
925 
926  friend inline bool
927  operator==(const Range &range_1, const Range &range_2)
928  {
929  return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
930  }
931 
932  static std::size_t
934  {
935  return sizeof(Range);
936  }
937 
943  template <class Archive>
944  void
945  serialize(Archive &ar, const unsigned int version);
946  };
947 
956  mutable std::vector<Range> ranges;
957 
966  mutable bool is_compressed;
967 
973 
984 
990 
994  void
995  do_compress() const;
996 };
997 
998 
1016 inline IndexSet
1018 {
1019  IndexSet is(N);
1020  is.add_range(0, N);
1021  is.compress();
1022  return is;
1023 }
1024 
1025 /* ------------------ inline functions ------------------ */
1026 
1027 
1028 /* IntervalAccessor */
1029 
1031  const IndexSet * idxset,
1033  : index_set(idxset)
1034  , range_idx(range_idx)
1035 {
1036  Assert(range_idx < idxset->n_intervals(),
1037  ExcInternalError("Invalid range index"));
1038 }
1039 
1040 
1041 
1043  : index_set(idxset)
1045 {}
1046 
1047 
1048 
1050  const IndexSet::IntervalAccessor &other)
1051  : index_set(other.index_set)
1052  , range_idx(other.range_idx)
1053 {
1055  ExcMessage("invalid iterator"));
1056 }
1057 
1058 
1059 
1060 inline IndexSet::size_type
1062 {
1063  Assert(is_valid(), ExcMessage("invalid iterator"));
1064  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1065 }
1066 
1067 
1068 
1069 inline bool
1071 {
1072  return index_set != nullptr && range_idx < index_set->n_intervals();
1073 }
1074 
1075 
1076 
1079 {
1080  Assert(is_valid(), ExcMessage("invalid iterator"));
1081  return {index_set, range_idx, index_set->ranges[range_idx].begin};
1082 }
1083 
1084 
1085 
1088 {
1089  Assert(is_valid(), ExcMessage("invalid iterator"));
1090 
1091  // point to first index in next interval unless we are the last interval.
1092  if (range_idx < index_set->ranges.size() - 1)
1093  return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1094  else
1095  return index_set->end();
1096 }
1097 
1098 
1099 
1100 inline IndexSet::size_type
1102 {
1103  Assert(is_valid(), ExcMessage("invalid iterator"));
1104 
1105  return index_set->ranges[range_idx].end - 1;
1106 }
1107 
1108 
1109 
1112 {
1113  index_set = other.index_set;
1114  range_idx = other.range_idx;
1116  ExcMessage("invalid iterator"));
1117  return *this;
1118 }
1119 
1120 
1121 
1122 inline bool
1124  const IndexSet::IntervalAccessor &other) const
1125 {
1126  Assert(index_set == other.index_set,
1127  ExcMessage(
1128  "Can not compare accessors pointing to different IndexSets"));
1129  return range_idx == other.range_idx;
1130 }
1131 
1132 
1133 
1134 inline bool
1136  const IndexSet::IntervalAccessor &other) const
1137 {
1138  Assert(index_set == other.index_set,
1139  ExcMessage(
1140  "Can not compare accessors pointing to different IndexSets"));
1141  return range_idx < other.range_idx;
1142 }
1143 
1144 
1145 
1146 inline void
1148 {
1149  Assert(
1150  is_valid(),
1151  ExcMessage(
1152  "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1153  ++range_idx;
1154 
1155  // set ourselves to invalid if we walk off the end
1156  if (range_idx >= index_set->ranges.size())
1158 }
1159 
1160 
1161 /* IntervalIterator */
1162 
1164  const IndexSet * idxset,
1166  : accessor(idxset, range_idx)
1167 {}
1168 
1169 
1170 
1172  : accessor(nullptr)
1173 {}
1174 
1175 
1176 
1178  : accessor(idxset)
1179 {}
1180 
1181 
1182 
1185 {
1186  accessor.advance();
1187  return *this;
1188 }
1189 
1190 
1191 
1194 {
1195  const IndexSet::IntervalIterator iter = *this;
1196  accessor.advance();
1197  return iter;
1198 }
1199 
1200 
1201 
1202 inline const IndexSet::IntervalAccessor &
1204 {
1205  return accessor;
1206 }
1207 
1208 
1209 
1210 inline const IndexSet::IntervalAccessor *
1212 {
1213  return &accessor;
1214 }
1215 
1216 
1217 
1218 inline bool
1220  const IndexSet::IntervalIterator &other) const
1221 {
1222  return accessor == other.accessor;
1223 }
1224 
1225 
1226 
1227 inline bool
1229  const IndexSet::IntervalIterator &other) const
1230 {
1231  return !(*this == other);
1232 }
1233 
1234 
1235 
1236 inline bool
1238  const IndexSet::IntervalIterator &other) const
1239 {
1240  return accessor < other.accessor;
1241 }
1242 
1243 
1244 
1245 inline int
1247  const IndexSet::IntervalIterator &other) const
1248 {
1250  ExcMessage(
1251  "Can not compare iterators belonging to different IndexSets"));
1252 
1254  accessor.index_set->ranges.size() :
1256  const size_type rhs =
1258  accessor.index_set->ranges.size() :
1259  other.accessor.range_idx;
1260 
1261  if (lhs > rhs)
1262  return static_cast<int>(lhs - rhs);
1263  else
1264  return -static_cast<int>(rhs - lhs);
1265 }
1266 
1267 
1268 
1269 /* ElementIterator */
1270 
1272  const IndexSet * idxset,
1273  const IndexSet::size_type range_idx,
1274  const IndexSet::size_type index)
1275  : index_set(idxset)
1276  , range_idx(range_idx)
1277  , idx(index)
1278 {
1279  Assert(range_idx < index_set->ranges.size(),
1280  ExcMessage(
1281  "Invalid range index for IndexSet::ElementIterator constructor."));
1282  Assert(
1283  idx >= index_set->ranges[range_idx].begin &&
1284  idx < index_set->ranges[range_idx].end,
1286  "Invalid index argument for IndexSet::ElementIterator constructor."));
1287 }
1288 
1289 
1290 
1292  : index_set(idxset)
1295 {}
1296 
1297 
1298 
1299 inline bool
1301 {
1304  (range_idx < index_set->ranges.size() &&
1305  idx < index_set->ranges[range_idx].end),
1306  ExcInternalError("Invalid ElementIterator state."));
1307 
1308  return (range_idx < index_set->ranges.size() &&
1309  idx < index_set->ranges[range_idx].end);
1310 }
1311 
1312 
1313 
1314 inline IndexSet::size_type
1316 {
1317  Assert(
1318  is_valid(),
1319  ExcMessage(
1320  "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1321  return idx;
1322 }
1323 
1324 
1325 
1326 inline bool
1328  const IndexSet::ElementIterator &other) const
1329 {
1330  Assert(index_set == other.index_set,
1331  ExcMessage(
1332  "Can not compare iterators belonging to different IndexSets"));
1333  return range_idx == other.range_idx && idx == other.idx;
1334 }
1335 
1336 
1337 
1338 inline void
1340 {
1341  Assert(
1342  is_valid(),
1343  ExcMessage(
1344  "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1345  if (idx < index_set->ranges[range_idx].end)
1346  ++idx;
1347  // end of this range?
1348  if (idx == index_set->ranges[range_idx].end)
1349  {
1350  // point to first element in next interval if possible
1351  if (range_idx < index_set->ranges.size() - 1)
1352  {
1353  ++range_idx;
1354  idx = index_set->ranges[range_idx].begin;
1355  }
1356  else
1357  {
1358  // we just fell off the end, set to invalid:
1361  }
1362  }
1363 }
1364 
1365 
1366 
1369 {
1370  advance();
1371  return *this;
1372 }
1373 
1374 
1375 
1378 {
1379  const IndexSet::ElementIterator it = *this;
1380  advance();
1381  return it;
1382 }
1383 
1384 
1385 
1386 inline bool
1388  const IndexSet::ElementIterator &other) const
1389 {
1390  return !(*this == other);
1391 }
1392 
1393 
1394 
1395 inline bool
1397  const IndexSet::ElementIterator &other) const
1398 {
1399  Assert(index_set == other.index_set,
1400  ExcMessage(
1401  "Can not compare iterators belonging to different IndexSets"));
1402  return range_idx < other.range_idx ||
1403  (range_idx == other.range_idx && idx < other.idx);
1404 }
1405 
1406 
1407 
1408 inline std::ptrdiff_t
1410  const IndexSet::ElementIterator &other) const
1411 {
1412  Assert(index_set == other.index_set,
1413  ExcMessage(
1414  "Can not compare iterators belonging to different IndexSets"));
1415  if (*this == other)
1416  return 0;
1417  if (!(*this < other))
1418  return -(other - *this);
1419 
1420  // only other can be equal to end() because of the checks above.
1422 
1423  // Note: we now compute how far advance *this in "*this < other" to get other,
1424  // so we need to return -c at the end.
1425 
1426  // first finish the current range:
1427  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1428 
1429  // now walk in steps of ranges (need to start one behind our current one):
1430  for (size_type range = range_idx + 1;
1431  range < index_set->ranges.size() && range <= other.range_idx;
1432  ++range)
1433  c += index_set->ranges[range].end - index_set->ranges[range].begin;
1434 
1435  Assert(
1436  other.range_idx < index_set->ranges.size() ||
1438  ExcMessage(
1439  "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1440 
1441  // We might have walked too far because we went until the end of
1442  // other.range_idx, so walk backwards to other.idx:
1444  c -= index_set->ranges[other.range_idx].end - other.idx;
1445 
1446  return -c;
1447 }
1448 
1449 
1450 /* Range */
1451 
1456 {}
1457 
1458 
1459 
1460 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1461  : begin(i1)
1462  , end(i2)
1464 {}
1465 
1466 
1467 
1468 /* IndexSet itself */
1469 
1471  : is_compressed(true)
1472  , index_space_size(0)
1474 {}
1475 
1476 
1477 
1479  : is_compressed(true)
1480  , index_space_size(size)
1482 {}
1483 
1484 
1485 
1486 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1487  : ranges(std::move(is.ranges))
1488  , is_compressed(is.is_compressed)
1489  , index_space_size(is.index_space_size)
1490  , largest_range(is.largest_range)
1491 {
1492  is.ranges.clear();
1493  is.is_compressed = true;
1494  is.index_space_size = 0;
1495  is.largest_range = numbers::invalid_unsigned_int;
1496 
1497  compress();
1498 }
1499 
1500 
1501 
1502 inline IndexSet &
1504 {
1505  ranges = std::move(is.ranges);
1506  is_compressed = is.is_compressed;
1507  index_space_size = is.index_space_size;
1508  largest_range = is.largest_range;
1509 
1510  is.ranges.clear();
1511  is.is_compressed = true;
1512  is.index_space_size = 0;
1513  is.largest_range = numbers::invalid_unsigned_int;
1514 
1515  compress();
1516 
1517  return *this;
1518 }
1519 
1520 
1521 
1524 {
1525  compress();
1526  if (ranges.size() > 0)
1527  return {this, 0, ranges[0].begin};
1528  else
1529  return end();
1530 }
1531 
1532 
1533 
1536 {
1537  compress();
1538  AssertIndexRange(global_index, size());
1539 
1540  if (ranges.empty())
1541  return end();
1542 
1543  std::vector<Range>::const_iterator main_range =
1544  ranges.begin() + largest_range;
1545 
1546  Range r(global_index, global_index + 1);
1547  // This optimization makes the bounds for lower_bound smaller by checking
1548  // the largest range first.
1549  std::vector<Range>::const_iterator range_begin, range_end;
1550  if (global_index < main_range->begin)
1551  {
1552  range_begin = ranges.begin();
1553  range_end = main_range;
1554  }
1555  else
1556  {
1557  range_begin = main_range;
1558  range_end = ranges.end();
1559  }
1560 
1561  // This will give us the first range p=[a,b[ with b>=global_index using
1562  // a binary search
1563  const std::vector<Range>::const_iterator p =
1564  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1565 
1566  // We couldn't find a range, which means we have no range that contains
1567  // global_index and also no range behind it, meaning we need to return end().
1568  if (p == ranges.end())
1569  return end();
1570 
1571  // Finally, we can have two cases: Either global_index is not in [a,b[,
1572  // which means we need to return an iterator to a because global_index, ...,
1573  // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
1574  // [a,b[ and we will return an iterator pointing directly at global_index
1575  // (else branch).
1576  if (global_index < p->begin)
1577  return {this, static_cast<size_type>(p - ranges.begin()), p->begin};
1578  else
1579  return {this, static_cast<size_type>(p - ranges.begin()), global_index};
1580 }
1581 
1582 
1583 
1586 {
1587  compress();
1588  return IndexSet::ElementIterator(this);
1589 }
1590 
1591 
1592 
1595 {
1596  compress();
1597  if (ranges.size() > 0)
1598  return IndexSet::IntervalIterator(this, 0);
1599  else
1600  return end_intervals();
1601 }
1602 
1603 
1604 
1607 {
1608  compress();
1609  return IndexSet::IntervalIterator(this);
1610 }
1611 
1612 
1613 
1614 inline void
1616 {
1617  // reset so that there are no indices in the set any more; however,
1618  // as documented, the index set retains its size
1619  ranges.clear();
1620  is_compressed = true;
1622 }
1623 
1624 
1625 
1626 inline void
1628 {
1629  Assert(ranges.empty(),
1630  ExcMessage("This function can only be called if the current "
1631  "object does not yet contain any elements."));
1632  index_space_size = sz;
1633  is_compressed = true;
1634 }
1635 
1636 
1637 
1638 inline IndexSet::size_type
1640 {
1641  return index_space_size;
1642 }
1643 
1644 
1645 
1646 inline void
1648 {
1649  if (is_compressed == true)
1650  return;
1651 
1652  do_compress();
1653 }
1654 
1655 
1656 
1657 inline void
1659 {
1661 
1662  const Range new_range(index, index + 1);
1663  if (ranges.size() == 0 || index > ranges.back().end)
1664  ranges.push_back(new_range);
1665  else if (index == ranges.back().end)
1666  ranges.back().end++;
1667  else
1668  ranges.insert(Utilities::lower_bound(ranges.begin(),
1669  ranges.end(),
1670  new_range),
1671  new_range);
1672  is_compressed = false;
1673 }
1674 
1675 
1676 
1677 inline void
1679 {
1680  Assert((begin < index_space_size) ||
1681  ((begin == index_space_size) && (end == index_space_size)),
1682  ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1683  Assert(end <= index_space_size,
1684  ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1685  AssertIndexRange(begin, end + 1);
1686 
1687  if (begin != end)
1688  {
1689  const Range new_range(begin, end);
1690 
1691  // the new index might be larger than the last index present in the
1692  // ranges. Then we can skip the binary search
1693  if (ranges.size() == 0 || begin > ranges.back().end)
1694  ranges.push_back(new_range);
1695  else
1696  ranges.insert(Utilities::lower_bound(ranges.begin(),
1697  ranges.end(),
1698  new_range),
1699  new_range);
1700  is_compressed = false;
1701  }
1702 }
1703 
1704 
1705 
1706 template <typename ForwardIterator>
1707 inline void
1708 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1709 {
1710  if (begin == end)
1711  return;
1712 
1713  // identify ranges in the given iterator range by checking whether some
1714  // indices happen to be consecutive. to avoid quadratic complexity when
1715  // calling add_range many times (as add_range() going into the middle of an
1716  // already existing range must shift entries around), we first collect a
1717  // vector of ranges.
1718  std::vector<std::pair<size_type, size_type>> tmp_ranges;
1719  bool ranges_are_sorted = true;
1720  for (ForwardIterator p = begin; p != end;)
1721  {
1722  const size_type begin_index = *p;
1723  size_type end_index = begin_index + 1;
1724  ForwardIterator q = p;
1725  ++q;
1726  while ((q != end) && (*q == end_index))
1727  {
1728  ++end_index;
1729  ++q;
1730  }
1731 
1732  tmp_ranges.emplace_back(begin_index, end_index);
1733  p = q;
1734 
1735  // if the starting index of the next go-around of the for loop is less
1736  // than the end index of the one just identified, then we will have at
1737  // least one pair of ranges that are not sorted, and consequently the
1738  // whole collection of ranges is not sorted.
1739  if (p != end && *p < end_index)
1740  ranges_are_sorted = false;
1741  }
1742 
1743  if (!ranges_are_sorted)
1744  std::sort(tmp_ranges.begin(), tmp_ranges.end());
1745 
1746  // if we have many ranges, we first construct a temporary index set (where
1747  // we add ranges in a consecutive way, so fast), otherwise, we work with
1748  // add_range(). the number 9 is chosen heuristically given the fact that
1749  // there are typically up to 8 independent ranges when adding the degrees of
1750  // freedom on a 3D cell or 9 when adding degrees of freedom of faces. if
1751  // doing cell-by-cell additions, we want to avoid repeated calls to
1752  // IndexSet::compress() which gets called upon merging two index sets, so we
1753  // want to be in the other branch then.
1754  if (tmp_ranges.size() > 9)
1755  {
1756  IndexSet tmp_set(size());
1757  tmp_set.ranges.reserve(tmp_ranges.size());
1758  for (const auto &i : tmp_ranges)
1759  tmp_set.add_range(i.first, i.second);
1760  this->add_indices(tmp_set);
1761  }
1762  else
1763  for (const auto &i : tmp_ranges)
1764  add_range(i.first, i.second);
1765 }
1766 
1767 
1768 
1769 inline bool
1771 {
1772  if (ranges.empty() == false)
1773  {
1774  compress();
1775 
1776  // fast check whether the index is in the largest range
1778  if (index >= ranges[largest_range].begin &&
1779  index < ranges[largest_range].end)
1780  return true;
1781 
1782  // get the element after which we would have to insert a range that
1783  // consists of all elements from this element to the end of the index
1784  // range plus one. after this call we know that if p!=end() then
1785  // p->begin<=index unless there is no such range at all
1786  //
1787  // if the searched for element is an element of this range, then we're
1788  // done. otherwise, the element can't be in one of the following ranges
1789  // because otherwise p would be a different iterator
1790  //
1791  // since we already know the position relative to the largest range (we
1792  // called compress!), we can perform the binary search on ranges with
1793  // lower/higher number compared to the largest range
1794  std::vector<Range>::const_iterator p = std::upper_bound(
1795  ranges.begin() +
1796  (index < ranges[largest_range].begin ? 0 : largest_range + 1),
1797  index < ranges[largest_range].begin ? ranges.begin() + largest_range :
1798  ranges.end(),
1799  Range(index, size() + 1));
1800 
1801  if (p == ranges.begin())
1802  return ((index >= p->begin) && (index < p->end));
1803 
1804  Assert((p == ranges.end()) || (p->begin > index), ExcInternalError());
1805 
1806  // now move to that previous range
1807  --p;
1808  Assert(p->begin <= index, ExcInternalError());
1809 
1810  return (p->end > index);
1811  }
1812 
1813  // didn't find this index, so it's not in the set
1814  return false;
1815 }
1816 
1817 
1818 
1819 inline bool
1821 {
1822  compress();
1823  return (ranges.size() <= 1);
1824 }
1825 
1826 
1827 
1828 inline bool
1830 {
1831  return ranges.empty();
1832 }
1833 
1834 
1835 
1836 inline IndexSet::size_type
1838 {
1839  // make sure we have non-overlapping ranges
1840  compress();
1841 
1842  size_type v = 0;
1843  if (!ranges.empty())
1844  {
1845  Range &r = ranges.back();
1846  v = r.nth_index_in_set + r.end - r.begin;
1847  }
1848 
1849 #ifdef DEBUG
1850  size_type s = 0;
1851  for (const auto &range : ranges)
1852  s += (range.end - range.begin);
1853  Assert(s == v, ExcInternalError());
1854 #endif
1855 
1856  return v;
1857 }
1858 
1859 
1860 
1861 inline unsigned int
1863 {
1864  compress();
1865  return ranges.size();
1866 }
1867 
1868 
1869 
1870 inline IndexSet::size_type
1872 {
1873  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1874 
1875  compress();
1876  const std::vector<Range>::const_iterator main_range =
1877  ranges.begin() + largest_range;
1878 
1879  return main_range->nth_index_in_set;
1880 }
1881 
1882 
1883 
1884 inline IndexSet::size_type
1886 {
1888 
1889  compress();
1890 
1891  // first check whether the index is in the largest range
1893  std::vector<Range>::const_iterator main_range =
1894  ranges.begin() + largest_range;
1895  if (n >= main_range->nth_index_in_set &&
1896  n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1897  return main_range->begin + (n - main_range->nth_index_in_set);
1898 
1899  // find out which chunk the local index n belongs to by using a binary
1900  // search. the comparator is based on the end of the ranges. Use the
1901  // position relative to main_range to subdivide the ranges
1902  Range r(n, n + 1);
1903  r.nth_index_in_set = n;
1904  std::vector<Range>::const_iterator range_begin, range_end;
1905  if (n < main_range->nth_index_in_set)
1906  {
1907  range_begin = ranges.begin();
1908  range_end = main_range;
1909  }
1910  else
1911  {
1912  range_begin = main_range + 1;
1913  range_end = ranges.end();
1914  }
1915 
1916  const std::vector<Range>::const_iterator p =
1917  Utilities::lower_bound(range_begin, range_end, r, Range::nth_index_compare);
1918 
1919  Assert(p != ranges.end(), ExcInternalError());
1920  return p->begin + (n - p->nth_index_in_set);
1921 }
1922 
1923 
1924 
1925 inline IndexSet::size_type
1927 {
1928  // to make this call thread-safe, compress() must not be called through this
1929  // function
1930  Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1931  AssertIndexRange(n, size());
1932 
1933  // return immediately if the index set is empty
1934  if (is_empty())
1936 
1937  // check whether the index is in the largest range. use the result to
1938  // perform a one-sided binary search afterward
1940  std::vector<Range>::const_iterator main_range =
1941  ranges.begin() + largest_range;
1942  if (n >= main_range->begin && n < main_range->end)
1943  return (n - main_range->begin) + main_range->nth_index_in_set;
1944 
1945  Range r(n, n);
1946  std::vector<Range>::const_iterator range_begin, range_end;
1947  if (n < main_range->begin)
1948  {
1949  range_begin = ranges.begin();
1950  range_end = main_range;
1951  }
1952  else
1953  {
1954  range_begin = main_range + 1;
1955  range_end = ranges.end();
1956  }
1957 
1958  std::vector<Range>::const_iterator p =
1959  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1960 
1961  // if n is not in this set
1962  if (p == range_end || p->end == n || p->begin > n)
1964 
1965  Assert(p != ranges.end(), ExcInternalError());
1966  Assert(p->begin <= n, ExcInternalError());
1967  Assert(n < p->end, ExcInternalError());
1968  return (n - p->begin) + p->nth_index_in_set;
1969 }
1970 
1971 
1972 
1973 inline bool
1975 {
1976  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1977 
1978  compress();
1979  is.compress();
1980 
1981  return ranges == is.ranges;
1982 }
1983 
1984 
1985 
1986 inline bool
1988 {
1989  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1990 
1991  compress();
1992  is.compress();
1993 
1994  return ranges != is.ranges;
1995 }
1996 
1997 
1998 
1999 template <typename Vector>
2000 void
2002 {
2003  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
2004 
2005  compress();
2006  // first fill all elements of the vector with zeroes.
2007  std::fill(vector.begin(), vector.end(), 0);
2008 
2009  // then write ones into the elements whose indices are contained in the
2010  // index set
2011  for (const auto &range : ranges)
2012  for (size_type i = range.begin; i < range.end; ++i)
2013  vector[i] = 1;
2014 }
2015 
2016 
2017 
2018 template <class StreamType>
2019 inline void
2020 IndexSet::print(StreamType &out) const
2021 {
2022  compress();
2023  out << "{";
2024  std::vector<Range>::const_iterator p;
2025  for (p = ranges.begin(); p != ranges.end(); ++p)
2026  {
2027  if (p->end - p->begin == 1)
2028  out << p->begin;
2029  else
2030  out << "[" << p->begin << "," << p->end - 1 << "]";
2031 
2032  if (p != --ranges.end())
2033  out << ", ";
2034  }
2035  out << "}" << std::endl;
2036 }
2037 
2038 
2039 
2040 template <class Archive>
2041 inline void
2042 IndexSet::Range::serialize(Archive &ar, const unsigned int)
2043 {
2044  ar &begin &end &nth_index_in_set;
2045 }
2046 
2047 
2048 
2049 template <class Archive>
2050 inline void
2051 IndexSet::serialize(Archive &ar, const unsigned int)
2052 {
2054 }
2055 
2057 
2058 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1111
friend bool operator==(const Range &range_1, const Range &range_2)
Definition: index_set.h:927
void fill_binary_vector(VectorType &vector) const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
ElementIterator end() const
Definition: index_set.h:1585
void add_index(const size_type index)
Definition: index_set.h:1658
const IntervalAccessor * operator->() const
Definition: index_set.h:1211
const IntervalAccessor & operator*() const
Definition: index_set.h:1203
const IndexSet * index_set
Definition: index_set.h:811
IndexSet & operator=(const IndexSet &)=default
void block_read(std::istream &in)
Definition: index_set.cc:495
size_type largest_range_starting_index() const
Definition: index_set.h:1871
#define AssertIndexRange(index, range)
Definition: exceptions.h:1718
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1409
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1708
types::global_dof_index size_type
Definition: index_set.h:83
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition: index_set.h:1271
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1237
std::forward_iterator_tag iterator_category
Definition: index_set.h:703
std::size_t memory_consumption() const
Definition: index_set.cc:740
size_type end
Definition: index_set.h:882
size_type n_elements() const
Definition: index_set.h:1061
size_type index_space_size
Definition: index_set.h:972
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:535
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:674
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:604
size_type size() const
Definition: index_set.h:1639
void print(StreamType &out) const
Definition: index_set.h:2020
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2042
void clear()
Definition: index_set.h:1615
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2051
static ::ExceptionBase & ExcMessage(std::string arg1)
void block_write(std::ostream &out) const
Definition: index_set.cc:481
std::vector< Range > ranges
Definition: index_set.h:956
static ::ExceptionBase & ExcIndexNotPresent(size_type arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
void subtract_set(const IndexSet &other)
Definition: index_set.cc:266
#define Assert(cond, exc)
Definition: exceptions.h:1461
friend class IntervalIterator
Definition: index_set.h:610
std::forward_iterator_tag iterator_category
Definition: index_set.h:795
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
bool operator<(const ElementIterator &) const
Definition: index_set.h:1396
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1926
bool operator==(const IndexSet &is) const
Definition: index_set.h:1974
signed int value_type
Definition: index_set.h:102
IntervalAccessor accessor
Definition: index_set.h:713
void do_compress() const
Definition: index_set.cc:98
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1111
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:914
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1219
size_type last() const
Definition: index_set.h:1101
const IndexSet * index_set
Definition: index_set.h:602
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1123
unsigned int n_intervals() const
Definition: index_set.h:1862
IndexSet tensor_product(const IndexSet &other) const
Definition: index_set.cc:338
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1246
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1987
size_type pop_back()
Definition: index_set.cc:350
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:515
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1678
IndexSet operator &(const IndexSet &is) const
void set_size(const size_type size)
Definition: index_set.h:1627
ElementIterator & operator++()
Definition: index_set.h:1368
unsigned int global_dof_index
Definition: types.h:76
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:1030
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
Definition: index_set.cc:239
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:212
void compress() const
Definition: index_set.h:1647
ElementIterator end() const
Definition: index_set.h:1087
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:920
static std::size_t memory_consumption()
Definition: index_set.h:933
IntervalIterator begin_intervals() const
Definition: index_set.h:1594
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1017
size_type begin
Definition: index_set.h:881
bool is_contiguous() const
Definition: index_set.h:1820
Threads::Mutex compress_mutex
Definition: index_set.h:989
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
size_type pop_front()
Definition: index_set.cc:368
static const char N
IntervalIterator end_intervals() const
Definition: index_set.h:1606
bool is_compressed
Definition: index_set.h:966
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1135
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1387
ElementIterator begin() const
Definition: index_set.h:1078
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
void write(std::ostream &out) const
Definition: index_set.cc:444
bool is_element(const size_type index) const
Definition: index_set.h:1770
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
bool operator==(const ElementIterator &) const
Definition: index_set.h:1327
ElementIterator at(const size_type global_index) const
Definition: index_set.h:1535
const types::global_dof_index invalid_dof_index
Definition: types.h:211
std::ptrdiff_t difference_type
Definition: index_set.h:797
std::ptrdiff_t difference_type
Definition: index_set.h:705
ElementIterator begin() const
Definition: index_set.h:1523
void read(std::istream &in)
Definition: index_set.cc:459
size_type n_elements() const
Definition: index_set.h:1837
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1885
bool is_empty() const
Definition: index_set.h:1829
size_type largest_range
Definition: index_set.h:983
IntervalIterator & operator++()
Definition: index_set.h:1184
size_type nth_index_in_set
Definition: index_set.h:884
static ::ExceptionBase & ExcInternalError()
size_type operator*() const
Definition: index_set.h:1315
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1228