Reference documentation for deal.II version GIT 196b8d1939 2022-07-03 13:35:01+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 
23 #include <deal.II/base/utilities.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 
40 #if defined(DEAL_II_WITH_MPI) || defined(DEAL_II_WITH_PETSC)
41 # include <mpi.h>
42 #else
43 using MPI_Comm = int;
44 # ifndef MPI_COMM_WORLD
45 # define MPI_COMM_WORLD 0
46 # endif
47 #endif
48 
50 
73 class IndexSet
74 {
75 public:
76  // forward declarations:
77  class ElementIterator;
78  class IntervalIterator;
79 
85 
103  using value_type = signed int;
104 
105 
109  IndexSet();
110 
114  explicit IndexSet(const size_type size);
115 
119  IndexSet(const IndexSet &) = default;
120 
124  IndexSet &
125  operator=(const IndexSet &) = default;
126 
131  IndexSet(IndexSet &&is) noexcept;
132 
137  IndexSet &
138  operator=(IndexSet &&is) noexcept;
139 
140 #ifdef DEAL_II_WITH_TRILINOS
144  explicit IndexSet(const Epetra_BlockMap &map);
145 #endif
146 
151  void
152  clear();
153 
160  void
161  set_size(const size_type size);
162 
170  size_type
171  size() const;
172 
179  void
180  add_range(const size_type begin, const size_type end);
181 
185  void
186  add_index(const size_type index);
187 
197  template <typename ForwardIterator>
198  void
199  add_indices(const ForwardIterator &begin, const ForwardIterator &end);
200 
216  void
217  add_indices(const IndexSet &other, const size_type offset = 0);
218 
222  bool
223  is_element(const size_type index) const;
224 
229  bool
230  is_contiguous() const;
231 
236  bool
237  is_empty() const;
238 
248  bool
249  is_ascending_and_one_to_one(const MPI_Comm &communicator) const;
250 
254  size_type
255  n_elements() const;
256 
262  size_type
263  nth_index_in_set(const size_type local_index) const;
264 
271  size_type
273 
282  unsigned int
283  n_intervals() const;
284 
296  size_type
298 
303  void
304  compress() const;
305 
311  bool
312  operator==(const IndexSet &is) const;
313 
319  bool
320  operator!=(const IndexSet &is) const;
321 
328  IndexSet
329  operator&(const IndexSet &is) const;
330 
344  IndexSet
345  get_view(const size_type begin, const size_type end) const;
346 
352  std::vector<IndexSet>
354  const std::vector<types::global_dof_index> &n_indices_per_block) const;
355 
361  void
362  subtract_set(const IndexSet &other);
363 
380  IndexSet
381  tensor_product(const IndexSet &other) const;
382 
387  size_type
388  pop_back();
389 
394  size_type
395  pop_front();
396 
400  void
401  fill_index_vector(std::vector<size_type> &indices) const;
402 
415  template <typename VectorType>
416  void
417  fill_binary_vector(VectorType &vector) const;
418 
423  template <class StreamType>
424  void
425  print(StreamType &out) const;
426 
431  void
432  write(std::ostream &out) const;
433 
438  void
439  read(std::istream &in);
440 
445  void
446  block_write(std::ostream &out) const;
447 
452  void
453  block_read(std::istream &in);
454 
455 #ifdef DEAL_II_WITH_TRILINOS
484  Epetra_Map
485  make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
486  const bool overlapping = false) const;
487 
488 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
489  Tpetra::Map<int, types::global_dof_index>
490  make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
491  const bool overlapping = false) const;
492 # endif
493 #endif
494 
495 
500  std::size_t
501  memory_consumption() const;
502 
504  size_type,
505  << "The global index " << arg1
506  << " is not an element of this set.");
507 
513  template <class Archive>
514  void
515  serialize(Archive &ar, const unsigned int version);
516 
517 
529  {
530  public:
535  IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
536 
540  explicit IntervalAccessor(const IndexSet *idxset);
541 
545  size_type
546  n_elements() const;
547 
551  bool
552  is_valid() const;
553 
558  begin() const;
559 
565  end() const;
566 
570  size_type
571  last() const;
572 
573  private:
577  IntervalAccessor(const IntervalAccessor &other);
582  operator=(const IntervalAccessor &other);
583 
587  bool
588  operator==(const IntervalAccessor &other) const;
592  bool
593  operator<(const IntervalAccessor &other) const;
598  void
599  advance();
604 
610 
611  friend class IntervalIterator;
612  };
613 
619  {
620  public:
625  IntervalIterator(const IndexSet *idxset, const size_type range_idx);
626 
630  explicit IntervalIterator(const IndexSet *idxset);
631 
636 
640  IntervalIterator(const IntervalIterator &other) = default;
641 
646  operator=(const IntervalIterator &other) = default;
647 
652  operator++();
653 
658  operator++(int);
659 
663  const IntervalAccessor &
664  operator*() const;
665 
669  const IntervalAccessor *
670  operator->() const;
671 
675  bool
676  operator==(const IntervalIterator &) const;
677 
681  bool
682  operator!=(const IntervalIterator &) const;
683 
687  bool
688  operator<(const IntervalIterator &) const;
689 
696  int
697  operator-(const IntervalIterator &p) const;
698 
704  using iterator_category = std::forward_iterator_tag;
706  using difference_type = std::ptrdiff_t;
709 
710  private:
715  };
716 
722  {
723  public:
728  ElementIterator(const IndexSet *idxset,
729  const size_type range_idx,
730  const size_type index);
731 
735  explicit ElementIterator(const IndexSet *idxset);
736 
741  size_type
742  operator*() const;
743 
747  bool
748  is_valid() const;
749 
754  operator++();
755 
760  operator++(int);
761 
765  bool
766  operator==(const ElementIterator &) const;
767 
771  bool
772  operator!=(const ElementIterator &) const;
773 
777  bool
778  operator<(const ElementIterator &) const;
779 
788  std::ptrdiff_t
789  operator-(const ElementIterator &p) const;
790 
796  using iterator_category = std::forward_iterator_tag;
798  using difference_type = std::ptrdiff_t;
799  using pointer = size_type *;
800  using reference = size_type &;
801 
802  private:
806  void
807  advance();
808 
821  };
822 
828  begin() const;
829 
845  at(const size_type global_index) const;
846 
852  end() const;
853 
858  begin_intervals() const;
859 
865  end_intervals() const;
866 
871 private:
880  struct Range
881  {
884 
886 
895  Range();
896 
904  Range(const size_type i1, const size_type i2);
905 
906  friend inline bool
907  operator<(const Range &range_1, const Range &range_2)
908  {
909  return (
910  (range_1.begin < range_2.begin) ||
911  ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
912  }
913 
914  static bool
916  {
917  return x.end < y.end;
918  }
919 
920  static bool
922  {
923  return (x.nth_index_in_set + (x.end - x.begin) <
924  y.nth_index_in_set + (y.end - y.begin));
925  }
926 
927  friend inline bool
928  operator==(const Range &range_1, const Range &range_2)
929  {
930  return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
931  }
932 
933  static std::size_t
935  {
936  return sizeof(Range);
937  }
938 
944  template <class Archive>
945  void
946  serialize(Archive &ar, const unsigned int version);
947  };
948 
957  mutable std::vector<Range> ranges;
958 
967  mutable bool is_compressed;
968 
974 
985 
991 
995  void
996  do_compress() const;
997 
1004  bool
1005  is_element_binary_search(const size_type local_index) const;
1006 
1013  size_type
1014  nth_index_in_set_binary_search(const size_type local_index) const;
1015 
1022  size_type
1024 
1030  void
1031  add_range_lower_bound(const Range &range);
1032 
1036  void
1038  boost::container::small_vector<std::pair<size_type, size_type>, 200>
1039  & tmp_ranges,
1040  const bool ranges_are_sorted);
1041 };
1042 
1043 
1044 
1062 inline IndexSet
1064 {
1065  IndexSet is(N);
1066  is.add_range(0, N);
1067  is.compress();
1068  return is;
1069 }
1070 
1071 /* ------------------ inline functions ------------------ */
1072 
1073 
1074 /* IntervalAccessor */
1075 
1077  const IndexSet * idxset,
1078  const IndexSet::size_type range_idx)
1079  : index_set(idxset)
1080  , range_idx(range_idx)
1081 {
1082  Assert(range_idx < idxset->n_intervals(),
1083  ExcInternalError("Invalid range index"));
1084 }
1085 
1086 
1087 
1089  : index_set(idxset)
1090  , range_idx(numbers::invalid_dof_index)
1091 {}
1092 
1093 
1094 
1096  const IndexSet::IntervalAccessor &other)
1097  : index_set(other.index_set)
1098  , range_idx(other.range_idx)
1099 {
1101  ExcMessage("invalid iterator"));
1102 }
1103 
1104 
1105 
1106 inline IndexSet::size_type
1108 {
1109  Assert(is_valid(), ExcMessage("invalid iterator"));
1110  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1111 }
1112 
1113 
1114 
1115 inline bool
1117 {
1118  return index_set != nullptr && range_idx < index_set->n_intervals();
1119 }
1120 
1121 
1122 
1125 {
1126  Assert(is_valid(), ExcMessage("invalid iterator"));
1127  return {index_set, range_idx, index_set->ranges[range_idx].begin};
1128 }
1129 
1130 
1131 
1134 {
1135  Assert(is_valid(), ExcMessage("invalid iterator"));
1136 
1137  // point to first index in next interval unless we are the last interval.
1138  if (range_idx < index_set->ranges.size() - 1)
1139  return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1140  else
1141  return index_set->end();
1142 }
1143 
1144 
1145 
1146 inline IndexSet::size_type
1148 {
1149  Assert(is_valid(), ExcMessage("invalid iterator"));
1150 
1151  return index_set->ranges[range_idx].end - 1;
1152 }
1153 
1154 
1155 
1158 {
1159  index_set = other.index_set;
1160  range_idx = other.range_idx;
1161  Assert(range_idx == numbers::invalid_dof_index || is_valid(),
1162  ExcMessage("invalid iterator"));
1163  return *this;
1164 }
1165 
1166 
1167 
1168 inline bool
1170  const IndexSet::IntervalAccessor &other) const
1171 {
1172  Assert(index_set == other.index_set,
1173  ExcMessage(
1174  "Can not compare accessors pointing to different IndexSets"));
1175  return range_idx == other.range_idx;
1176 }
1177 
1178 
1179 
1180 inline bool
1182  const IndexSet::IntervalAccessor &other) const
1183 {
1184  Assert(index_set == other.index_set,
1185  ExcMessage(
1186  "Can not compare accessors pointing to different IndexSets"));
1187  return range_idx < other.range_idx;
1188 }
1189 
1190 
1191 
1192 inline void
1194 {
1195  Assert(
1196  is_valid(),
1197  ExcMessage(
1198  "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1199  ++range_idx;
1200 
1201  // set ourselves to invalid if we walk off the end
1202  if (range_idx >= index_set->ranges.size())
1203  range_idx = numbers::invalid_dof_index;
1204 }
1205 
1206 
1207 /* IntervalIterator */
1208 
1210  const IndexSet * idxset,
1211  const IndexSet::size_type range_idx)
1212  : accessor(idxset, range_idx)
1213 {}
1214 
1215 
1216 
1218  : accessor(nullptr)
1219 {}
1220 
1221 
1222 
1224  : accessor(idxset)
1225 {}
1226 
1227 
1228 
1231 {
1232  accessor.advance();
1233  return *this;
1234 }
1235 
1236 
1237 
1240 {
1241  const IndexSet::IntervalIterator iter = *this;
1242  accessor.advance();
1243  return iter;
1244 }
1245 
1246 
1247 
1248 inline const IndexSet::IntervalAccessor &
1250 {
1251  return accessor;
1252 }
1253 
1254 
1255 
1256 inline const IndexSet::IntervalAccessor *
1258 {
1259  return &accessor;
1260 }
1261 
1262 
1263 
1264 inline bool
1266  const IndexSet::IntervalIterator &other) const
1267 {
1268  return accessor == other.accessor;
1269 }
1270 
1271 
1272 
1273 inline bool
1275  const IndexSet::IntervalIterator &other) const
1276 {
1277  return !(*this == other);
1278 }
1279 
1280 
1281 
1282 inline bool
1284  const IndexSet::IntervalIterator &other) const
1285 {
1286  return accessor < other.accessor;
1287 }
1288 
1289 
1290 
1291 inline int
1293  const IndexSet::IntervalIterator &other) const
1294 {
1295  Assert(accessor.index_set == other.accessor.index_set,
1296  ExcMessage(
1297  "Can not compare iterators belonging to different IndexSets"));
1298 
1299  const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index) ?
1300  accessor.index_set->ranges.size() :
1301  accessor.range_idx;
1302  const size_type rhs =
1304  accessor.index_set->ranges.size() :
1305  other.accessor.range_idx;
1306 
1307  if (lhs > rhs)
1308  return static_cast<int>(lhs - rhs);
1309  else
1310  return -static_cast<int>(rhs - lhs);
1311 }
1312 
1313 
1314 
1315 /* ElementIterator */
1316 
1318  const IndexSet * idxset,
1319  const IndexSet::size_type range_idx,
1320  const IndexSet::size_type index)
1321  : index_set(idxset)
1322  , range_idx(range_idx)
1323  , idx(index)
1324 {
1325  Assert(range_idx < index_set->ranges.size(),
1326  ExcMessage(
1327  "Invalid range index for IndexSet::ElementIterator constructor."));
1328  Assert(
1329  idx >= index_set->ranges[range_idx].begin &&
1330  idx < index_set->ranges[range_idx].end,
1332  "Invalid index argument for IndexSet::ElementIterator constructor."));
1333 }
1334 
1335 
1336 
1338  : index_set(idxset)
1339  , range_idx(numbers::invalid_dof_index)
1340  , idx(numbers::invalid_dof_index)
1341 {}
1342 
1343 
1344 
1345 inline bool
1347 {
1348  Assert((range_idx == numbers::invalid_dof_index &&
1349  idx == numbers::invalid_dof_index) ||
1350  (range_idx < index_set->ranges.size() &&
1351  idx < index_set->ranges[range_idx].end),
1352  ExcInternalError("Invalid ElementIterator state."));
1353 
1354  return (range_idx < index_set->ranges.size() &&
1355  idx < index_set->ranges[range_idx].end);
1356 }
1357 
1358 
1359 
1360 inline IndexSet::size_type
1362 {
1363  Assert(
1364  is_valid(),
1365  ExcMessage(
1366  "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1367  return idx;
1368 }
1369 
1370 
1371 
1372 inline bool
1374  const IndexSet::ElementIterator &other) const
1375 {
1376  Assert(index_set == other.index_set,
1377  ExcMessage(
1378  "Can not compare iterators belonging to different IndexSets"));
1379  return range_idx == other.range_idx && idx == other.idx;
1380 }
1381 
1382 
1383 
1384 inline void
1386 {
1387  Assert(
1388  is_valid(),
1389  ExcMessage(
1390  "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1391  if (idx < index_set->ranges[range_idx].end)
1392  ++idx;
1393  // end of this range?
1394  if (idx == index_set->ranges[range_idx].end)
1395  {
1396  // point to first element in next interval if possible
1397  if (range_idx < index_set->ranges.size() - 1)
1398  {
1399  ++range_idx;
1400  idx = index_set->ranges[range_idx].begin;
1401  }
1402  else
1403  {
1404  // we just fell off the end, set to invalid:
1405  range_idx = numbers::invalid_dof_index;
1407  }
1408  }
1409 }
1410 
1411 
1412 
1415 {
1416  advance();
1417  return *this;
1418 }
1419 
1420 
1421 
1424 {
1425  const IndexSet::ElementIterator it = *this;
1426  advance();
1427  return it;
1428 }
1429 
1430 
1431 
1432 inline bool
1434  const IndexSet::ElementIterator &other) const
1435 {
1436  return !(*this == other);
1437 }
1438 
1439 
1440 
1441 inline bool
1443  const IndexSet::ElementIterator &other) const
1444 {
1445  Assert(index_set == other.index_set,
1446  ExcMessage(
1447  "Can not compare iterators belonging to different IndexSets"));
1448  return range_idx < other.range_idx ||
1449  (range_idx == other.range_idx && idx < other.idx);
1450 }
1451 
1452 
1453 
1454 inline std::ptrdiff_t
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  if (*this == other)
1462  return 0;
1463  if (!(*this < other))
1464  return -(other - *this);
1465 
1466  // only other can be equal to end() because of the checks above.
1467  Assert(is_valid(), ExcInternalError());
1468 
1469  // Note: we now compute how far advance *this in "*this < other" to get other,
1470  // so we need to return -c at the end.
1471 
1472  // first finish the current range:
1473  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1474 
1475  // now walk in steps of ranges (need to start one behind our current one):
1476  for (size_type range = range_idx + 1;
1477  range < index_set->ranges.size() && range <= other.range_idx;
1478  ++range)
1479  c += index_set->ranges[range].end - index_set->ranges[range].begin;
1480 
1481  Assert(
1482  other.range_idx < index_set->ranges.size() ||
1484  ExcMessage(
1485  "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1486 
1487  // We might have walked too far because we went until the end of
1488  // other.range_idx, so walk backwards to other.idx:
1490  c -= index_set->ranges[other.range_idx].end - other.idx;
1491 
1492  return -c;
1493 }
1494 
1495 
1496 /* Range */
1497 
1502 {}
1503 
1504 
1505 
1506 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1507  : begin(i1)
1508  , end(i2)
1510 {}
1511 
1512 
1513 
1514 /* IndexSet itself */
1515 
1517  : is_compressed(true)
1518  , index_space_size(0)
1520 {}
1521 
1522 
1523 
1524 inline IndexSet::IndexSet(const size_type size)
1525  : is_compressed(true)
1526  , index_space_size(size)
1527  , largest_range(numbers::invalid_unsigned_int)
1528 {}
1529 
1530 
1531 
1532 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1533  : ranges(std::move(is.ranges))
1534  , is_compressed(is.is_compressed)
1535  , index_space_size(is.index_space_size)
1536  , largest_range(is.largest_range)
1537 {
1538  is.ranges.clear();
1539  is.is_compressed = true;
1540  is.index_space_size = 0;
1541  is.largest_range = numbers::invalid_unsigned_int;
1542 
1543  compress();
1544 }
1545 
1546 
1547 
1548 inline IndexSet &
1550 {
1551  ranges = std::move(is.ranges);
1552  is_compressed = is.is_compressed;
1553  index_space_size = is.index_space_size;
1554  largest_range = is.largest_range;
1555 
1556  is.ranges.clear();
1557  is.is_compressed = true;
1558  is.index_space_size = 0;
1559  is.largest_range = numbers::invalid_unsigned_int;
1560 
1561  compress();
1562 
1563  return *this;
1564 }
1565 
1566 
1567 
1570 {
1571  compress();
1572  if (ranges.size() > 0)
1573  return {this, 0, ranges[0].begin};
1574  else
1575  return end();
1576 }
1577 
1578 
1579 
1582 {
1583  compress();
1584  return IndexSet::ElementIterator(this);
1585 }
1586 
1587 
1588 
1591 {
1592  compress();
1593  if (ranges.size() > 0)
1594  return IndexSet::IntervalIterator(this, 0);
1595  else
1596  return end_intervals();
1597 }
1598 
1599 
1600 
1603 {
1604  compress();
1605  return IndexSet::IntervalIterator(this);
1606 }
1607 
1608 
1609 
1610 inline void
1612 {
1613  // reset so that there are no indices in the set any more; however,
1614  // as documented, the index set retains its size
1615  ranges.clear();
1616  is_compressed = true;
1618 }
1619 
1620 
1621 
1622 inline void
1624 {
1625  Assert(ranges.empty(),
1626  ExcMessage("This function can only be called if the current "
1627  "object does not yet contain any elements."));
1628  index_space_size = sz;
1629  is_compressed = true;
1630 }
1631 
1632 
1633 
1634 inline IndexSet::size_type
1636 {
1637  return index_space_size;
1638 }
1639 
1640 
1641 
1642 inline void
1644 {
1645  if (is_compressed == true)
1646  return;
1647 
1648  do_compress();
1649 }
1650 
1651 
1652 
1653 inline void
1655 {
1657 
1658  const Range new_range(index, index + 1);
1659  if (ranges.size() == 0 || index > ranges.back().end)
1660  ranges.push_back(new_range);
1661  else if (index == ranges.back().end)
1662  ranges.back().end++;
1663  else
1664  add_range_lower_bound(new_range);
1665  is_compressed = false;
1666 }
1667 
1668 
1669 
1670 inline void
1672 {
1674  ((begin == index_space_size) && (end == index_space_size)),
1675  ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1677  ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1678  AssertIndexRange(begin, end + 1);
1679 
1680  if (begin != end)
1681  {
1682  const Range new_range(begin, end);
1683 
1684  // the new index might be larger than the last index present in the
1685  // ranges. Then we can skip the binary search
1686  if (ranges.size() == 0 || begin > ranges.back().end)
1687  ranges.push_back(new_range);
1688  else if (begin == ranges.back().end)
1689  ranges.back().end = end;
1690  else
1691  add_range_lower_bound(new_range);
1692 
1693  is_compressed = false;
1694  }
1695 }
1696 
1697 
1698 
1699 template <typename ForwardIterator>
1700 inline void
1701 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1702 {
1703  if (begin == end)
1704  return;
1705 
1706  // identify ranges in the given iterator range by checking whether some
1707  // indices happen to be consecutive. to avoid quadratic complexity when
1708  // calling add_range many times (as add_range() going into the middle of an
1709  // already existing range must shift entries around), we first collect a
1710  // vector of ranges.
1711  boost::container::small_vector<std::pair<size_type, size_type>, 200>
1712  tmp_ranges;
1713  bool ranges_are_sorted = true;
1714  for (ForwardIterator p = begin; p != end;)
1715  {
1716  const size_type begin_index = *p;
1717  size_type end_index = begin_index + 1;
1718  ForwardIterator q = p;
1719  ++q;
1720  while ((q != end) && (static_cast<size_type>(*q) == end_index))
1721  {
1722  ++end_index;
1723  ++q;
1724  }
1725 
1726  tmp_ranges.emplace_back(begin_index, end_index);
1727  p = q;
1728 
1729  // if the starting index of the next go-around of the for loop is less
1730  // than the end index of the one just identified, then we will have at
1731  // least one pair of ranges that are not sorted, and consequently the
1732  // whole collection of ranges is not sorted.
1733  if (p != end && static_cast<size_type>(*p) < end_index)
1734  ranges_are_sorted = false;
1735  }
1736 
1737  add_ranges_internal(tmp_ranges, ranges_are_sorted);
1738 }
1739 
1740 
1741 
1742 inline bool
1744 {
1745  if (ranges.empty() == false)
1746  {
1747  compress();
1748 
1749  // fast check whether the index is in the largest range
1751  if (index >= ranges[largest_range].begin &&
1752  index < ranges[largest_range].end)
1753  return true;
1754  else if (ranges.size() > 1)
1755  return is_element_binary_search(index);
1756  else
1757  return false;
1758  }
1759  else
1760  return false;
1761 }
1762 
1763 
1764 
1765 inline bool
1767 {
1768  compress();
1769  return (ranges.size() <= 1);
1770 }
1771 
1772 
1773 
1774 inline bool
1776 {
1777  return ranges.empty();
1778 }
1779 
1780 
1781 
1782 inline IndexSet::size_type
1784 {
1785  // make sure we have non-overlapping ranges
1786  compress();
1787 
1788  size_type v = 0;
1789  if (!ranges.empty())
1790  {
1791  Range &r = ranges.back();
1792  v = r.nth_index_in_set + r.end - r.begin;
1793  }
1794 
1795 #ifdef DEBUG
1796  size_type s = 0;
1797  for (const auto &range : ranges)
1798  s += (range.end - range.begin);
1799  Assert(s == v, ExcInternalError());
1800 #endif
1801 
1802  return v;
1803 }
1804 
1805 
1806 
1807 inline unsigned int
1809 {
1810  compress();
1811  return ranges.size();
1812 }
1813 
1814 
1815 
1816 inline IndexSet::size_type
1818 {
1819  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1820 
1821  compress();
1822  const std::vector<Range>::const_iterator main_range =
1823  ranges.begin() + largest_range;
1824 
1825  return main_range->nth_index_in_set;
1826 }
1827 
1828 
1829 
1830 inline IndexSet::size_type
1832 {
1834 
1835  compress();
1836 
1837  // first check whether the index is in the largest range
1839  const auto main_range = ranges.begin() + largest_range;
1840  if (n >= main_range->nth_index_in_set &&
1841  n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1842  return main_range->begin + (n - main_range->nth_index_in_set);
1843  else
1845 }
1846 
1847 
1848 
1849 inline IndexSet::size_type
1851 {
1852  // to make this call thread-safe, compress() must not be called through this
1853  // function
1854  Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1855  AssertIndexRange(n, size());
1856 
1857  // return immediately if the index set is empty
1858  if (is_empty())
1860 
1861  // check whether the index is in the largest range. use the result to
1862  // perform a one-sided binary search afterward
1864  if (n >= ranges[largest_range].begin && n < ranges[largest_range].end)
1865  return (n - ranges[largest_range].begin) +
1866  ranges[largest_range].nth_index_in_set;
1867  else if (ranges.size() > 1)
1869  else
1871 }
1872 
1873 
1874 
1875 inline bool
1877 {
1878  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1879 
1880  compress();
1881  is.compress();
1882 
1883  return ranges == is.ranges;
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 template <typename Vector>
1902 void
1904 {
1905  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
1906 
1907  compress();
1908  // first fill all elements of the vector with zeroes.
1909  std::fill(vector.begin(), vector.end(), 0);
1910 
1911  // then write ones into the elements whose indices are contained in the
1912  // index set
1913  for (const auto &range : ranges)
1914  for (size_type i = range.begin; i < range.end; ++i)
1915  vector[i] = 1;
1916 }
1917 
1918 
1919 
1920 template <class StreamType>
1921 inline void
1922 IndexSet::print(StreamType &out) const
1923 {
1924  compress();
1925  out << '{';
1926  std::vector<Range>::const_iterator p;
1927  for (p = ranges.begin(); p != ranges.end(); ++p)
1928  {
1929  if (p->end - p->begin == 1)
1930  out << p->begin;
1931  else
1932  out << '[' << p->begin << ',' << p->end - 1 << ']';
1933 
1934  if (p != --ranges.end())
1935  out << ", ";
1936  }
1937  out << '}' << std::endl;
1938 }
1939 
1940 
1941 
1942 template <class Archive>
1943 inline void
1944 IndexSet::Range::serialize(Archive &ar, const unsigned int)
1945 {
1946  ar &begin &end &nth_index_in_set;
1947 }
1948 
1949 
1950 
1951 template <class Archive>
1952 inline void
1953 IndexSet::serialize(Archive &ar, const unsigned int)
1954 {
1956 }
1957 
1959 
1960 #endif
size_type operator*() const
Definition: index_set.h:1361
bool operator==(const ElementIterator &) const
Definition: index_set.h:1373
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition: index_set.h:1317
std::ptrdiff_t difference_type
Definition: index_set.h:798
ElementIterator & operator++()
Definition: index_set.h:1414
bool operator<(const ElementIterator &) const
Definition: index_set.h:1442
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1433
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1455
std::forward_iterator_tag iterator_category
Definition: index_set.h:796
const IndexSet * index_set
Definition: index_set.h:812
size_type last() const
Definition: index_set.h:1147
ElementIterator end() const
Definition: index_set.h:1133
size_type n_elements() const
Definition: index_set.h:1107
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1169
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1157
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1181
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:1076
const IndexSet * index_set
Definition: index_set.h:603
ElementIterator begin() const
Definition: index_set.h:1124
std::ptrdiff_t difference_type
Definition: index_set.h:706
const IntervalAccessor & operator*() const
Definition: index_set.h:1249
IntervalIterator(const IntervalIterator &other)=default
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1265
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1292
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1283
std::forward_iterator_tag iterator_category
Definition: index_set.h:704
IntervalIterator & operator++()
Definition: index_set.h:1230
IntervalIterator & operator=(const IntervalIterator &other)=default
const IntervalAccessor * operator->() const
Definition: index_set.h:1257
IntervalAccessor accessor
Definition: index_set.h:714
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1274
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:984
bool is_contiguous() const
Definition: index_set.h:1766
unsigned int n_intervals() const
Definition: index_set.h:1808
IndexSet(const IndexSet &)=default
IndexSet & operator=(const IndexSet &)=default
IntervalIterator end_intervals() const
Definition: index_set.h:1602
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:1635
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:688
bool is_empty() const
Definition: index_set.h:1775
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1850
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:1783
bool operator==(const IndexSet &is) const
Definition: index_set.h:1876
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:1743
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1953
ElementIterator begin() const
Definition: index_set.h:1569
size_type pop_front()
Definition: index_set.cc:365
signed int value_type
Definition: index_set.h:103
void set_size(const size_type size)
Definition: index_set.h:1623
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1889
void read(std::istream &in)
Definition: index_set.cc:507
void clear()
Definition: index_set.h:1611
size_type largest_range_starting_index() const
Definition: index_set.h:1817
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:1654
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:967
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:1590
void fill_binary_vector(VectorType &vector) const
std::vector< Range > ranges
Definition: index_set.h:957
void subtract_set(const IndexSet &other)
Definition: index_set.cc:266
ElementIterator end() const
Definition: index_set.h:1581
Threads::Mutex compress_mutex
Definition: index_set.h:990
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1063
size_type index_space_size
Definition: index_set.h:973
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:1671
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1831
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:1922
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:1643
size_type pop_back()
Definition: index_set.cc:347
types::global_dof_index size_type
Definition: index_set.h:84
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1701
IndexSet operator&(const IndexSet &is) const
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:456
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:495
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:1473
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcMessage(std::string arg1)
iterator end()
size_type size() const
iterator begin()
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:201
const types::global_dof_index invalid_dof_index
Definition: types.h:216
unsigned int global_dof_index
Definition: types.h:76
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:915
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:921
friend bool operator==(const Range &range_1, const Range &range_2)
Definition: index_set.h:928
size_type end
Definition: index_set.h:883
size_type nth_index_in_set
Definition: index_set.h:885
static std::size_t memory_consumption()
Definition: index_set.h:934
size_type begin
Definition: index_set.h:882
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1944
friend bool operator<(const Range &range_1, const Range &range_2)
Definition: index_set.h:907
void advance(std::tuple< I1, I2 > &t, const unsigned int n)