Reference documentation for deal.II version Git 423cd11810 2020-05-28 01:07:25 +0200
\(\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 - 2020 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 
74 class IndexSet
75 {
76 public:
77  // forward declarations:
78  class ElementIterator;
79  class IntervalIterator;
80 
86 
104  using value_type = signed int;
105 
106 
110  IndexSet();
111 
115  explicit IndexSet(const size_type size);
116 
120  IndexSet(const IndexSet &) = default;
121 
125  IndexSet &
126  operator=(const IndexSet &) = default;
127 
132  IndexSet(IndexSet &&is) noexcept;
133 
138  IndexSet &
139  operator=(IndexSet &&is) noexcept;
140 
141 #ifdef DEAL_II_WITH_TRILINOS
142 
145  explicit IndexSet(const Epetra_BlockMap &map);
146 #endif
147 
152  void
153  clear();
154 
161  void
162  set_size(const size_type size);
163 
171  size_type
172  size() const;
173 
180  void
181  add_range(const size_type begin, const size_type end);
182 
186  void
187  add_index(const size_type index);
188 
198  template <typename ForwardIterator>
199  void
200  add_indices(const ForwardIterator &begin, const ForwardIterator &end);
201 
217  void
218  add_indices(const IndexSet &other, const size_type offset = 0);
219 
223  bool
224  is_element(const size_type index) const;
225 
230  bool
231  is_contiguous() const;
232 
237  bool
238  is_empty() const;
239 
249  bool
250  is_ascending_and_one_to_one(const MPI_Comm &communicator) const;
251 
255  size_type
256  n_elements() const;
257 
263  size_type
264  nth_index_in_set(const size_type local_index) const;
265 
272  size_type
274 
283  unsigned int
284  n_intervals() const;
285 
297  size_type
299 
304  void
305  compress() const;
306 
312  bool
313  operator==(const IndexSet &is) const;
314 
320  bool
321  operator!=(const IndexSet &is) const;
322 
329  IndexSet operator&(const IndexSet &is) const;
330 
344  IndexSet
345  get_view(const size_type begin, const size_type end) const;
346 
353  std::vector<IndexSet>
355  const std::vector<types::global_dof_index> &n_indices_per_block) const;
356 
362  void
363  subtract_set(const IndexSet &other);
364 
381  IndexSet
382  tensor_product(const IndexSet &other) const;
383 
388  size_type
389  pop_back();
390 
395  size_type
396  pop_front();
397 
401  void
402  fill_index_vector(std::vector<size_type> &indices) const;
403 
416  template <typename VectorType>
417  void
418  fill_binary_vector(VectorType &vector) const;
419 
424  template <class StreamType>
425  void
426  print(StreamType &out) const;
427 
432  void
433  write(std::ostream &out) const;
434 
439  void
440  read(std::istream &in);
441 
446  void
447  block_write(std::ostream &out) const;
448 
453  void
454  block_read(std::istream &in);
455 
456 #ifdef DEAL_II_WITH_TRILINOS
457 
485  Epetra_Map
486  make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
487  const bool overlapping = false) const;
488 
489 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
490  Tpetra::Map<int, types::global_dof_index>
491  make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
492  const bool overlapping = false) const;
493 # endif
494 #endif
495 
496 
501  std::size_t
502  memory_consumption() const;
503 
505  size_type,
506  << "The global index " << arg1
507  << " is not an element of this set.");
508 
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 &operator*() const;
664 
668  const IntervalAccessor *operator->() const;
669 
673  bool
674  operator==(const IntervalIterator &) const;
675 
679  bool
680  operator!=(const IntervalIterator &) const;
681 
685  bool
686  operator<(const IntervalIterator &) const;
687 
694  int
695  operator-(const IntervalIterator &p) const;
696 
702  using iterator_category = std::forward_iterator_tag;
704  using difference_type = std::ptrdiff_t;
707 
708  private:
713  };
714 
720  {
721  public:
726  ElementIterator(const IndexSet *idxset,
727  const size_type range_idx,
728  const size_type index);
729 
733  explicit ElementIterator(const IndexSet *idxset);
734 
739  size_type operator*() const;
740 
744  bool
745  is_valid() const;
746 
751  operator++();
752 
757  operator++(int);
758 
762  bool
763  operator==(const ElementIterator &) const;
764 
768  bool
769  operator!=(const ElementIterator &) const;
770 
774  bool
775  operator<(const ElementIterator &) const;
776 
785  std::ptrdiff_t
786  operator-(const ElementIterator &p) const;
787 
793  using iterator_category = std::forward_iterator_tag;
795  using difference_type = std::ptrdiff_t;
796  using pointer = size_type *;
797  using reference = size_type &;
798 
799  private:
803  void
804  advance();
805 
818  };
819 
825  begin() const;
826 
842  at(const size_type global_index) const;
843 
849  end() const;
850 
855  begin_intervals() const;
856 
862  end_intervals() const;
863 
868 private:
877  struct Range
878  {
881 
883 
892  Range();
893 
901  Range(const size_type i1, const size_type i2);
902 
903  friend inline bool
904  operator<(const Range &range_1, const Range &range_2)
905  {
906  return (
907  (range_1.begin < range_2.begin) ||
908  ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
909  }
910 
911  static bool
913  {
914  return x.end < y.end;
915  }
916 
917  static bool
919  {
920  return (x.nth_index_in_set + (x.end - x.begin) <
921  y.nth_index_in_set + (y.end - y.begin));
922  }
923 
924  friend inline bool
925  operator==(const Range &range_1, const Range &range_2)
926  {
927  return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
928  }
929 
930  static std::size_t
932  {
933  return sizeof(Range);
934  }
935 
940  template <class Archive>
941  void
942  serialize(Archive &ar, const unsigned int version);
943  };
944 
953  mutable std::vector<Range> ranges;
954 
963  mutable bool is_compressed;
964 
970 
981 
987 
991  void
992  do_compress() const;
993 };
994 
995 
1013 inline IndexSet
1015 {
1016  IndexSet is(N);
1017  is.add_range(0, N);
1018  is.compress();
1019  return is;
1020 }
1021 
1022 /* ------------------ inline functions ------------------ */
1023 
1024 
1025 /* IntervalAccessor */
1026 
1028  const IndexSet * idxset,
1030  : index_set(idxset)
1031  , range_idx(range_idx)
1032 {
1033  Assert(range_idx < idxset->n_intervals(),
1034  ExcInternalError("Invalid range index"));
1035 }
1036 
1037 
1038 
1040  : index_set(idxset)
1042 {}
1043 
1044 
1045 
1047  const IndexSet::IntervalAccessor &other)
1048  : index_set(other.index_set)
1049  , range_idx(other.range_idx)
1050 {
1052  ExcMessage("invalid iterator"));
1053 }
1054 
1055 
1056 
1057 inline IndexSet::size_type
1059 {
1060  Assert(is_valid(), ExcMessage("invalid iterator"));
1061  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1062 }
1063 
1064 
1065 
1066 inline bool
1068 {
1069  return index_set != nullptr && range_idx < index_set->n_intervals();
1070 }
1071 
1072 
1073 
1076 {
1077  Assert(is_valid(), ExcMessage("invalid iterator"));
1078  return {index_set, range_idx, index_set->ranges[range_idx].begin};
1079 }
1080 
1081 
1082 
1085 {
1086  Assert(is_valid(), ExcMessage("invalid iterator"));
1087 
1088  // point to first index in next interval unless we are the last interval.
1089  if (range_idx < index_set->ranges.size() - 1)
1090  return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1091  else
1092  return index_set->end();
1093 }
1094 
1095 
1096 
1097 inline IndexSet::size_type
1099 {
1100  Assert(is_valid(), ExcMessage("invalid iterator"));
1101 
1102  return index_set->ranges[range_idx].end - 1;
1103 }
1104 
1105 
1106 
1109 {
1110  index_set = other.index_set;
1111  range_idx = other.range_idx;
1113  ExcMessage("invalid iterator"));
1114  return *this;
1115 }
1116 
1117 
1118 
1119 inline bool
1122 {
1123  Assert(index_set == other.index_set,
1124  ExcMessage(
1125  "Can not compare accessors pointing to different IndexSets"));
1126  return range_idx == other.range_idx;
1127 }
1128 
1129 
1130 
1131 inline bool
1134 {
1135  Assert(index_set == other.index_set,
1136  ExcMessage(
1137  "Can not compare accessors pointing to different IndexSets"));
1138  return range_idx < other.range_idx;
1139 }
1140 
1141 
1142 
1143 inline void
1145 {
1146  Assert(
1147  is_valid(),
1148  ExcMessage(
1149  "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1150  ++range_idx;
1151 
1152  // set ourselves to invalid if we walk off the end
1153  if (range_idx >= index_set->ranges.size())
1155 }
1156 
1157 
1158 /* IntervalIterator */
1159 
1161  const IndexSet * idxset,
1163  : accessor(idxset, range_idx)
1164 {}
1165 
1166 
1167 
1169  : accessor(nullptr)
1170 {}
1171 
1172 
1173 
1175  : accessor(idxset)
1176 {}
1177 
1178 
1179 
1182 {
1183  accessor.advance();
1184  return *this;
1185 }
1186 
1187 
1188 
1191 {
1192  const IndexSet::IntervalIterator iter = *this;
1193  accessor.advance();
1194  return iter;
1195 }
1196 
1197 
1198 
1200  operator*() const
1201 {
1202  return accessor;
1203 }
1204 
1205 
1206 
1209 {
1210  return &accessor;
1211 }
1212 
1213 
1214 
1215 inline bool
1218 {
1219  return accessor == other.accessor;
1220 }
1221 
1222 
1223 
1224 inline bool
1227 {
1228  return !(*this == other);
1229 }
1230 
1231 
1232 
1233 inline bool
1236 {
1237  return accessor < other.accessor;
1238 }
1239 
1240 
1241 
1242 inline int
1245 {
1247  ExcMessage(
1248  "Can not compare iterators belonging to different IndexSets"));
1249 
1251  accessor.index_set->ranges.size() :
1253  const size_type rhs =
1255  accessor.index_set->ranges.size() :
1256  other.accessor.range_idx;
1257 
1258  if (lhs > rhs)
1259  return static_cast<int>(lhs - rhs);
1260  else
1261  return -static_cast<int>(rhs - lhs);
1262 }
1263 
1264 
1265 
1266 /* ElementIterator */
1267 
1269  const IndexSet * idxset,
1270  const IndexSet::size_type range_idx,
1271  const IndexSet::size_type index)
1272  : index_set(idxset)
1273  , range_idx(range_idx)
1274  , idx(index)
1275 {
1276  Assert(range_idx < index_set->ranges.size(),
1277  ExcMessage(
1278  "Invalid range index for IndexSet::ElementIterator constructor."));
1279  Assert(
1280  idx >= index_set->ranges[range_idx].begin &&
1281  idx < index_set->ranges[range_idx].end,
1283  "Invalid index argument for IndexSet::ElementIterator constructor."));
1284 }
1285 
1286 
1287 
1289  : index_set(idxset)
1292 {}
1293 
1294 
1295 
1296 inline bool
1298 {
1301  (range_idx < index_set->ranges.size() &&
1302  idx < index_set->ranges[range_idx].end),
1303  ExcInternalError("Invalid ElementIterator state."));
1304 
1305  return (range_idx < index_set->ranges.size() &&
1306  idx < index_set->ranges[range_idx].end);
1307 }
1308 
1309 
1310 
1312 {
1313  Assert(
1314  is_valid(),
1315  ExcMessage(
1316  "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1317  return idx;
1318 }
1319 
1320 
1321 
1322 inline bool
1325 {
1326  Assert(index_set == other.index_set,
1327  ExcMessage(
1328  "Can not compare iterators belonging to different IndexSets"));
1329  return range_idx == other.range_idx && idx == other.idx;
1330 }
1331 
1332 
1333 
1334 inline void
1336 {
1337  Assert(
1338  is_valid(),
1339  ExcMessage(
1340  "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1341  if (idx < index_set->ranges[range_idx].end)
1342  ++idx;
1343  // end of this range?
1344  if (idx == index_set->ranges[range_idx].end)
1345  {
1346  // point to first element in next interval if possible
1347  if (range_idx < index_set->ranges.size() - 1)
1348  {
1349  ++range_idx;
1350  idx = index_set->ranges[range_idx].begin;
1351  }
1352  else
1353  {
1354  // we just fell off the end, set to invalid:
1357  }
1358  }
1359 }
1360 
1361 
1362 
1365 {
1366  advance();
1367  return *this;
1368 }
1369 
1370 
1371 
1374 {
1375  const IndexSet::ElementIterator it = *this;
1376  advance();
1377  return it;
1378 }
1379 
1380 
1381 
1382 inline bool
1385 {
1386  return !(*this == other);
1387 }
1388 
1389 
1390 
1391 inline bool
1394 {
1395  Assert(index_set == other.index_set,
1396  ExcMessage(
1397  "Can not compare iterators belonging to different IndexSets"));
1398  return range_idx < other.range_idx ||
1399  (range_idx == other.range_idx && idx < other.idx);
1400 }
1401 
1402 
1403 
1404 inline std::ptrdiff_t
1407 {
1408  Assert(index_set == other.index_set,
1409  ExcMessage(
1410  "Can not compare iterators belonging to different IndexSets"));
1411  if (*this == other)
1412  return 0;
1413  if (!(*this < other))
1414  return -(other - *this);
1415 
1416  // only other can be equal to end() because of the checks above.
1418 
1419  // Note: we now compute how far advance *this in "*this < other" to get other,
1420  // so we need to return -c at the end.
1421 
1422  // first finish the current range:
1423  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1424 
1425  // now walk in steps of ranges (need to start one behind our current one):
1426  for (size_type range = range_idx + 1;
1427  range < index_set->ranges.size() && range <= other.range_idx;
1428  ++range)
1429  c += index_set->ranges[range].end - index_set->ranges[range].begin;
1430 
1431  Assert(
1432  other.range_idx < index_set->ranges.size() ||
1434  ExcMessage(
1435  "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1436 
1437  // We might have walked too far because we went until the end of
1438  // other.range_idx, so walk backwards to other.idx:
1440  c -= index_set->ranges[other.range_idx].end - other.idx;
1441 
1442  return -c;
1443 }
1444 
1445 
1446 /* Range */
1447 
1452 {}
1453 
1454 
1455 
1456 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1457  : begin(i1)
1458  , end(i2)
1460 {}
1461 
1462 
1463 
1464 /* IndexSet itself */
1465 
1467  : is_compressed(true)
1468  , index_space_size(0)
1470 {}
1471 
1472 
1473 
1475  : is_compressed(true)
1476  , index_space_size(size)
1478 {}
1479 
1480 
1481 
1482 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1483  : ranges(std::move(is.ranges))
1484  , is_compressed(is.is_compressed)
1485  , index_space_size(is.index_space_size)
1486  , largest_range(is.largest_range)
1487 {
1488  is.ranges.clear();
1489  is.is_compressed = true;
1490  is.index_space_size = 0;
1491  is.largest_range = numbers::invalid_unsigned_int;
1492 
1493  compress();
1494 }
1495 
1496 
1497 
1498 inline IndexSet &
1500 {
1501  ranges = std::move(is.ranges);
1502  is_compressed = is.is_compressed;
1503  index_space_size = is.index_space_size;
1504  largest_range = is.largest_range;
1505 
1506  is.ranges.clear();
1507  is.is_compressed = true;
1508  is.index_space_size = 0;
1509  is.largest_range = numbers::invalid_unsigned_int;
1510 
1511  compress();
1512 
1513  return *this;
1514 }
1515 
1516 
1517 
1520 {
1521  compress();
1522  if (ranges.size() > 0)
1523  return {this, 0, ranges[0].begin};
1524  else
1525  return end();
1526 }
1527 
1528 
1529 
1532 {
1533  compress();
1534  AssertIndexRange(global_index, size());
1535 
1536  if (ranges.empty())
1537  return end();
1538 
1539  std::vector<Range>::const_iterator main_range =
1540  ranges.begin() + largest_range;
1541 
1542  Range r(global_index, global_index + 1);
1543  // This optimization makes the bounds for lower_bound smaller by checking
1544  // the largest range first.
1545  std::vector<Range>::const_iterator range_begin, range_end;
1546  if (global_index < main_range->begin)
1547  {
1548  range_begin = ranges.begin();
1549  range_end = main_range;
1550  }
1551  else
1552  {
1553  range_begin = main_range;
1554  range_end = ranges.end();
1555  }
1556 
1557  // This will give us the first range p=[a,b[ with b>=global_index using
1558  // a binary search
1559  const std::vector<Range>::const_iterator p =
1560  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1561 
1562  // We couldn't find a range, which means we have no range that contains
1563  // global_index and also no range behind it, meaning we need to return end().
1564  if (p == ranges.end())
1565  return end();
1566 
1567  // Finally, we can have two cases: Either global_index is not in [a,b[,
1568  // which means we need to return an iterator to a because global_index, ...,
1569  // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
1570  // [a,b[ and we will return an iterator pointing directly at global_index
1571  // (else branch).
1572  if (global_index < p->begin)
1573  return {this, static_cast<size_type>(p - ranges.begin()), p->begin};
1574  else
1575  return {this, static_cast<size_type>(p - ranges.begin()), global_index};
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  ranges.insert(Utilities::lower_bound(ranges.begin(),
1665  ranges.end(),
1666  new_range),
1667  new_range);
1668  is_compressed = false;
1669 }
1670 
1671 
1672 
1673 inline void
1675 {
1676  Assert((begin < index_space_size) ||
1677  ((begin == index_space_size) && (end == index_space_size)),
1678  ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1679  Assert(end <= index_space_size,
1680  ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1681  AssertIndexRange(begin, end + 1);
1682 
1683  if (begin != end)
1684  {
1685  const Range new_range(begin, end);
1686 
1687  // the new index might be larger than the last index present in the
1688  // ranges. Then we can skip the binary search
1689  if (ranges.size() == 0 || begin > ranges.back().end)
1690  ranges.push_back(new_range);
1691  else
1692  ranges.insert(Utilities::lower_bound(ranges.begin(),
1693  ranges.end(),
1694  new_range),
1695  new_range);
1696  is_compressed = false;
1697  }
1698 }
1699 
1700 
1701 
1702 template <typename ForwardIterator>
1703 inline void
1704 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1705 {
1706  if (begin == end)
1707  return;
1708 
1709  // identify ranges in the given iterator range by checking whether some
1710  // indices happen to be consecutive. to avoid quadratic complexity when
1711  // calling add_range many times (as add_range() going into the middle of an
1712  // already existing range must shift entries around), we first collect a
1713  // vector of ranges.
1714  std::vector<std::pair<size_type, size_type>> tmp_ranges;
1715  bool ranges_are_sorted = true;
1716  for (ForwardIterator p = begin; p != end;)
1717  {
1718  const size_type begin_index = *p;
1719  size_type end_index = begin_index + 1;
1720  ForwardIterator q = p;
1721  ++q;
1722  while ((q != end) && (*q == end_index))
1723  {
1724  ++end_index;
1725  ++q;
1726  }
1727 
1728  tmp_ranges.emplace_back(begin_index, end_index);
1729  p = q;
1730 
1731  // if the starting index of the next go-around of the for loop is less
1732  // than the end index of the one just identified, then we will have at
1733  // least one pair of ranges that are not sorted, and consequently the
1734  // whole collection of ranges is not sorted.
1735  if (p != end && *p < end_index)
1736  ranges_are_sorted = false;
1737  }
1738 
1739  if (!ranges_are_sorted)
1740  std::sort(tmp_ranges.begin(), tmp_ranges.end());
1741 
1742  // if we have many ranges, we first construct a temporary index set (where
1743  // we add ranges in a consecutive way, so fast), otherwise, we work with
1744  // add_range(). the number 9 is chosen heuristically given the fact that
1745  // there are typically up to 8 independent ranges when adding the degrees of
1746  // freedom on a 3D cell or 9 when adding degrees of freedom of faces. if
1747  // doing cell-by-cell additions, we want to avoid repeated calls to
1748  // IndexSet::compress() which gets called upon merging two index sets, so we
1749  // want to be in the other branch then.
1750  if (tmp_ranges.size() > 9)
1751  {
1752  IndexSet tmp_set(size());
1753  tmp_set.ranges.reserve(tmp_ranges.size());
1754  for (const auto &i : tmp_ranges)
1755  tmp_set.add_range(i.first, i.second);
1756  this->add_indices(tmp_set);
1757  }
1758  else
1759  for (const auto &i : tmp_ranges)
1760  add_range(i.first, i.second);
1761 }
1762 
1763 
1764 
1765 inline bool
1767 {
1768  if (ranges.empty() == false)
1769  {
1770  compress();
1771 
1772  // fast check whether the index is in the largest range
1774  if (index >= ranges[largest_range].begin &&
1775  index < ranges[largest_range].end)
1776  return true;
1777 
1778  // get the element after which we would have to insert a range that
1779  // consists of all elements from this element to the end of the index
1780  // range plus one. after this call we know that if p!=end() then
1781  // p->begin<=index unless there is no such range at all
1782  //
1783  // if the searched for element is an element of this range, then we're
1784  // done. otherwise, the element can't be in one of the following ranges
1785  // because otherwise p would be a different iterator
1786  //
1787  // since we already know the position relative to the largest range (we
1788  // called compress!), we can perform the binary search on ranges with
1789  // lower/higher number compared to the largest range
1790  std::vector<Range>::const_iterator p = std::upper_bound(
1791  ranges.begin() +
1792  (index < ranges[largest_range].begin ? 0 : largest_range + 1),
1793  index < ranges[largest_range].begin ? ranges.begin() + largest_range :
1794  ranges.end(),
1795  Range(index, size() + 1));
1796 
1797  if (p == ranges.begin())
1798  return ((index >= p->begin) && (index < p->end));
1799 
1800  Assert((p == ranges.end()) || (p->begin > index), ExcInternalError());
1801 
1802  // now move to that previous range
1803  --p;
1804  Assert(p->begin <= index, ExcInternalError());
1805 
1806  return (p->end > index);
1807  }
1808 
1809  // didn't find this index, so it's not in the set
1810  return false;
1811 }
1812 
1813 
1814 
1815 inline bool
1817 {
1818  compress();
1819  return (ranges.size() <= 1);
1820 }
1821 
1822 
1823 
1824 inline bool
1826 {
1827  return ranges.empty();
1828 }
1829 
1830 
1831 
1832 inline IndexSet::size_type
1834 {
1835  // make sure we have non-overlapping ranges
1836  compress();
1837 
1838  size_type v = 0;
1839  if (!ranges.empty())
1840  {
1841  Range &r = ranges.back();
1842  v = r.nth_index_in_set + r.end - r.begin;
1843  }
1844 
1845 #ifdef DEBUG
1846  size_type s = 0;
1847  for (const auto &range : ranges)
1848  s += (range.end - range.begin);
1849  Assert(s == v, ExcInternalError());
1850 #endif
1851 
1852  return v;
1853 }
1854 
1855 
1856 
1857 inline unsigned int
1859 {
1860  compress();
1861  return ranges.size();
1862 }
1863 
1864 
1865 
1866 inline IndexSet::size_type
1868 {
1869  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1870 
1871  compress();
1872  const std::vector<Range>::const_iterator main_range =
1873  ranges.begin() + largest_range;
1874 
1875  return main_range->nth_index_in_set;
1876 }
1877 
1878 
1879 
1880 inline IndexSet::size_type
1882 {
1884 
1885  compress();
1886 
1887  // first check whether the index is in the largest range
1889  std::vector<Range>::const_iterator main_range =
1890  ranges.begin() + largest_range;
1891  if (n >= main_range->nth_index_in_set &&
1892  n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1893  return main_range->begin + (n - main_range->nth_index_in_set);
1894 
1895  // find out which chunk the local index n belongs to by using a binary
1896  // search. the comparator is based on the end of the ranges. Use the
1897  // position relative to main_range to subdivide the ranges
1898  Range r(n, n + 1);
1899  r.nth_index_in_set = n;
1900  std::vector<Range>::const_iterator range_begin, range_end;
1901  if (n < main_range->nth_index_in_set)
1902  {
1903  range_begin = ranges.begin();
1904  range_end = main_range;
1905  }
1906  else
1907  {
1908  range_begin = main_range + 1;
1909  range_end = ranges.end();
1910  }
1911 
1912  const std::vector<Range>::const_iterator p =
1913  Utilities::lower_bound(range_begin, range_end, r, Range::nth_index_compare);
1914 
1915  Assert(p != ranges.end(), ExcInternalError());
1916  return p->begin + (n - p->nth_index_in_set);
1917 }
1918 
1919 
1920 
1921 inline IndexSet::size_type
1923 {
1924  // to make this call thread-safe, compress() must not be called through this
1925  // function
1926  Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1927  AssertIndexRange(n, size());
1928 
1929  // return immediately if the index set is empty
1930  if (is_empty())
1932 
1933  // check whether the index is in the largest range. use the result to
1934  // perform a one-sided binary search afterward
1936  std::vector<Range>::const_iterator main_range =
1937  ranges.begin() + largest_range;
1938  if (n >= main_range->begin && n < main_range->end)
1939  return (n - main_range->begin) + main_range->nth_index_in_set;
1940 
1941  Range r(n, n);
1942  std::vector<Range>::const_iterator range_begin, range_end;
1943  if (n < main_range->begin)
1944  {
1945  range_begin = ranges.begin();
1946  range_end = main_range;
1947  }
1948  else
1949  {
1950  range_begin = main_range + 1;
1951  range_end = ranges.end();
1952  }
1953 
1954  std::vector<Range>::const_iterator p =
1955  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1956 
1957  // if n is not in this set
1958  if (p == range_end || p->end == n || p->begin > n)
1960 
1961  Assert(p != ranges.end(), ExcInternalError());
1962  Assert(p->begin <= n, ExcInternalError());
1963  Assert(n < p->end, ExcInternalError());
1964  return (n - p->begin) + p->nth_index_in_set;
1965 }
1966 
1967 
1968 
1969 inline bool
1971 {
1972  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1973 
1974  compress();
1975  is.compress();
1976 
1977  return ranges == is.ranges;
1978 }
1979 
1980 
1981 
1982 inline bool
1984 {
1985  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1986 
1987  compress();
1988  is.compress();
1989 
1990  return ranges != is.ranges;
1991 }
1992 
1993 
1994 
1995 template <typename Vector>
1996 void
1998 {
1999  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
2000 
2001  compress();
2002  // first fill all elements of the vector with zeroes.
2003  std::fill(vector.begin(), vector.end(), 0);
2004 
2005  // then write ones into the elements whose indices are contained in the
2006  // index set
2007  for (const auto &range : ranges)
2008  for (size_type i = range.begin; i < range.end; ++i)
2009  vector[i] = 1;
2010 }
2011 
2012 
2013 
2014 template <class StreamType>
2015 inline void
2016 IndexSet::print(StreamType &out) const
2017 {
2018  compress();
2019  out << "{";
2020  std::vector<Range>::const_iterator p;
2021  for (p = ranges.begin(); p != ranges.end(); ++p)
2022  {
2023  if (p->end - p->begin == 1)
2024  out << p->begin;
2025  else
2026  out << "[" << p->begin << "," << p->end - 1 << "]";
2027 
2028  if (p != --ranges.end())
2029  out << ", ";
2030  }
2031  out << "}" << std::endl;
2032 }
2033 
2034 
2035 
2036 template <class Archive>
2037 inline void
2038 IndexSet::Range::serialize(Archive &ar, const unsigned int)
2039 {
2040  ar &begin &end &nth_index_in_set;
2041 }
2042 
2043 
2044 
2045 template <class Archive>
2046 inline void
2047 IndexSet::serialize(Archive &ar, const unsigned int)
2048 {
2050 }
2051 
2053 
2054 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1102
friend bool operator==(const Range &range_1, const Range &range_2)
Definition: index_set.h:925
void fill_binary_vector(VectorType &vector) const
static const unsigned int invalid_unsigned_int
Definition: types.h:191
ElementIterator end() const
Definition: index_set.h:1581
void add_index(const size_type index)
Definition: index_set.h:1654
const IntervalAccessor * operator->() const
Definition: index_set.h:1208
const IntervalAccessor & operator*() const
Definition: index_set.h:1200
const IndexSet * index_set
Definition: index_set.h:809
IndexSet & operator=(const IndexSet &)=default
void block_read(std::istream &in)
Definition: index_set.cc:487
size_type largest_range_starting_index() const
Definition: index_set.h:1867
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1406
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1704
types::global_dof_index size_type
Definition: index_set.h:85
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition: index_set.h:1268
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1235
std::forward_iterator_tag iterator_category
Definition: index_set.h:702
std::size_t memory_consumption() const
Definition: index_set.cc:729
size_type end
Definition: index_set.h:880
size_type n_elements() const
Definition: index_set.h:1058
size_type index_space_size
Definition: index_set.h:969
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:527
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:666
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:596
size_type size() const
Definition: index_set.h:1635
void print(StreamType &out) const
Definition: index_set.h:2016
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2038
void clear()
Definition: index_set.h:1611
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2047
static ::ExceptionBase & ExcMessage(std::string arg1)
void block_write(std::ostream &out) const
Definition: index_set.cc:473
std::vector< Range > ranges
Definition: index_set.h:953
static ::ExceptionBase & ExcIndexNotPresent(size_type arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
void subtract_set(const IndexSet &other)
Definition: index_set.cc:258
#define Assert(cond, exc)
Definition: exceptions.h:1419
friend class IntervalIterator
Definition: index_set.h:611
std::forward_iterator_tag iterator_category
Definition: index_set.h:793
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
bool operator<(const ElementIterator &) const
Definition: index_set.h:1393
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1922
bool operator==(const IndexSet &is) const
Definition: index_set.h:1970
signed int value_type
Definition: index_set.h:104
IntervalAccessor accessor
Definition: index_set.h:712
void do_compress() const
Definition: index_set.cc:98
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:362
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1108
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:912
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1217
size_type last() const
Definition: index_set.h:1098
const IndexSet * index_set
Definition: index_set.h:603
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1121
unsigned int n_intervals() const
Definition: index_set.h:1858
IndexSet tensor_product(const IndexSet &other) const
Definition: index_set.cc:330
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1244
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1983
size_type pop_back()
Definition: index_set.cc:342
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:507
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1674
IndexSet operator &(const IndexSet &is) const
void set_size(const size_type size)
Definition: index_set.h:1623
ElementIterator & operator++()
Definition: index_set.h:1364
unsigned int global_dof_index
Definition: types.h:76
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:1027
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
Definition: index_set.cc:238
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:211
void compress() const
Definition: index_set.h:1643
ElementIterator end() const
Definition: index_set.h:1084
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:918
static std::size_t memory_consumption()
Definition: index_set.h:931
IntervalIterator begin_intervals() const
Definition: index_set.h:1590
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1014
size_type begin
Definition: index_set.h:879
bool is_contiguous() const
Definition: index_set.h:1816
Threads::Mutex compress_mutex
Definition: index_set.h:986
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:361
size_type pop_front()
Definition: index_set.cc:360
static const char N
IntervalIterator end_intervals() const
Definition: index_set.h:1602
bool is_compressed
Definition: index_set.h:963
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1133
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1384
ElementIterator begin() const
Definition: index_set.h:1075
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:436
bool is_element(const size_type index) const
Definition: index_set.h:1766
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:1324
ElementIterator at(const size_type global_index) const
Definition: index_set.h:1531
const types::global_dof_index invalid_dof_index
Definition: types.h:206
std::ptrdiff_t difference_type
Definition: index_set.h:795
std::ptrdiff_t difference_type
Definition: index_set.h:704
ElementIterator begin() const
Definition: index_set.h:1519
void read(std::istream &in)
Definition: index_set.cc:451
size_type n_elements() const
Definition: index_set.h:1833
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1881
bool is_empty() const
Definition: index_set.h:1825
size_type largest_range
Definition: index_set.h:980
IntervalIterator & operator++()
Definition: index_set.h:1181
size_type nth_index_in_set
Definition: index_set.h:882
static ::ExceptionBase & ExcInternalError()
size_type operator*() const
Definition: index_set.h:1311
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1226