Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
index_set.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2019 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 
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/thread_management.h>
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 
48 DEAL_II_NAMESPACE_OPEN
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 unsigned int 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 unsigned int local_index) const;
265 
272  size_type
273  index_within_set(const size_type global_index) const;
274 
283  unsigned int
284  n_intervals() const;
285 
290  unsigned int
292 
297  void
298  compress() const;
299 
305  bool
306  operator==(const IndexSet &is) const;
307 
313  bool
314  operator!=(const IndexSet &is) const;
315 
322  IndexSet operator&(const IndexSet &is) const;
323 
337  IndexSet
338  get_view(const size_type begin, const size_type end) const;
339 
345  void
346  subtract_set(const IndexSet &other);
347 
352  size_type
353  pop_back();
354 
359  size_type
360  pop_front();
361 
365  void
366  fill_index_vector(std::vector<size_type> &indices) const;
367 
380  template <typename VectorType>
381  void
382  fill_binary_vector(VectorType &vector) const;
383 
388  template <class StreamType>
389  void
390  print(StreamType &out) const;
391 
396  void
397  write(std::ostream &out) const;
398 
403  void
404  read(std::istream &in);
405 
410  void
411  block_write(std::ostream &out) const;
412 
417  void
418  block_read(std::istream &in);
419 
420 #ifdef DEAL_II_WITH_TRILINOS
421 
449  Epetra_Map
450  make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
451  const bool overlapping = false) const;
452 
453 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
454  Tpetra::Map<int, types::global_dof_index>
455  make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
456  const bool overlapping = false) const;
457 # endif
458 #endif
459 
460 
465  std::size_t
466  memory_consumption() const;
467 
469  size_type,
470  << "The global index " << arg1
471  << " is not an element of this set.");
472 
477  template <class Archive>
478  void
479  serialize(Archive &ar, const unsigned int version);
480 
481 
493  {
494  public:
499  IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
500 
504  explicit IntervalAccessor(const IndexSet *idxset);
505 
509  size_type
510  n_elements() const;
511 
515  bool
516  is_valid() const;
517 
522  begin() const;
523 
529  end() const;
530 
534  size_type
535  last() const;
536 
537  private:
541  IntervalAccessor(const IntervalAccessor &other);
546  operator=(const IntervalAccessor &other);
547 
551  bool
552  operator==(const IntervalAccessor &other) const;
556  bool
557  operator<(const IntervalAccessor &other) const;
562  void
563  advance();
568 
574 
575  friend class IntervalIterator;
576  };
577 
583  {
584  public:
589  IntervalIterator(const IndexSet *idxset, const size_type range_idx);
590 
594  explicit IntervalIterator(const IndexSet *idxset);
595 
600 
604  IntervalIterator(const IntervalIterator &other) = default;
605 
610  operator=(const IntervalIterator &other) = default;
611 
616  operator++();
617 
622  operator++(int);
623 
627  const IntervalAccessor &operator*() const;
628 
632  const IntervalAccessor *operator->() const;
633 
637  bool
638  operator==(const IntervalIterator &) const;
639 
643  bool
644  operator!=(const IntervalIterator &) const;
645 
649  bool
650  operator<(const IntervalIterator &) const;
651 
658  int
659  operator-(const IntervalIterator &p) const;
660 
666  using iterator_category = std::forward_iterator_tag;
668  using difference_type = std::ptrdiff_t;
669  using pointer = IntervalAccessor *;
670  using reference = IntervalAccessor &;
671 
672  private:
677  };
678 
684  {
685  public:
690  ElementIterator(const IndexSet *idxset,
691  const size_type range_idx,
692  const size_type index);
693 
697  explicit ElementIterator(const IndexSet *idxset);
698 
703  size_type operator*() const;
704 
708  bool
709  is_valid() const;
710 
715  operator++();
716 
721  operator++(int);
722 
726  bool
727  operator==(const ElementIterator &) const;
728 
732  bool
733  operator!=(const ElementIterator &) const;
734 
738  bool
739  operator<(const ElementIterator &) const;
740 
749  std::ptrdiff_t
750  operator-(const ElementIterator &p) const;
751 
757  using iterator_category = std::forward_iterator_tag;
758  using value_type = size_type;
759  using difference_type = std::ptrdiff_t;
760  using pointer = size_type *;
761  using reference = size_type &;
762 
763  private:
767  void
768  advance();
769 
782  };
783 
789  begin() const;
790 
806  at(const size_type global_index) const;
807 
813  end() const;
814 
819  begin_intervals() const;
820 
826  end_intervals() const;
827 
832 private:
841  struct Range
842  {
843  size_type begin;
844  size_type end;
845 
846  size_type nth_index_in_set;
847 
856  Range();
857 
865  Range(const size_type i1, const size_type i2);
866 
867  friend inline bool
868  operator<(const Range &range_1, const Range &range_2)
869  {
870  return (
871  (range_1.begin < range_2.begin) ||
872  ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
873  }
874 
875  static bool
876  end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
877  {
878  return x.end < y.end;
879  }
880 
881  static bool
882  nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
883  {
884  return (x.nth_index_in_set + (x.end - x.begin) <
885  y.nth_index_in_set + (y.end - y.begin));
886  }
887 
888  friend inline bool
889  operator==(const Range &range_1, const Range &range_2)
890  {
891  return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
892  }
893 
894  static std::size_t
895  memory_consumption()
896  {
897  return sizeof(Range);
898  }
899 
904  template <class Archive>
905  void
906  serialize(Archive &ar, const unsigned int version);
907  };
908 
917  mutable std::vector<Range> ranges;
918 
927  mutable bool is_compressed;
928 
934 
945 
951 
955  void
956  do_compress() const;
957 };
958 
959 
977 inline IndexSet
978 complete_index_set(const unsigned int N)
979 {
980  IndexSet is(N);
981  is.add_range(0, N);
982  is.compress();
983  return is;
984 }
985 
986 /* ------------------ inline functions ------------------ */
987 
988 
989 /* IntervalAccessor */
990 
992  const IndexSet * idxset,
993  const IndexSet::size_type range_idx)
994  : index_set(idxset)
995  , range_idx(range_idx)
996 {
997  Assert(range_idx < idxset->n_intervals(),
998  ExcInternalError("Invalid range index"));
999 }
1000 
1001 
1002 
1004  : index_set(idxset)
1005  , range_idx(numbers::invalid_dof_index)
1006 {}
1007 
1008 
1009 
1011  const IndexSet::IntervalAccessor &other)
1012  : index_set(other.index_set)
1013  , range_idx(other.range_idx)
1014 {
1016  ExcMessage("invalid iterator"));
1017 }
1018 
1019 
1020 
1021 inline IndexSet::size_type
1023 {
1024  Assert(is_valid(), ExcMessage("invalid iterator"));
1025  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1026 }
1027 
1028 
1029 
1030 inline bool
1032 {
1033  return index_set != nullptr && range_idx < index_set->n_intervals();
1034 }
1035 
1036 
1037 
1040 {
1041  Assert(is_valid(), ExcMessage("invalid iterator"));
1042  return {index_set, range_idx, index_set->ranges[range_idx].begin};
1043 }
1044 
1045 
1046 
1049 {
1050  Assert(is_valid(), ExcMessage("invalid iterator"));
1051 
1052  // point to first index in next interval unless we are the last interval.
1053  if (range_idx < index_set->ranges.size() - 1)
1054  return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1055  else
1056  return index_set->end();
1057 }
1058 
1059 
1060 
1061 inline IndexSet::size_type
1063 {
1064  Assert(is_valid(), ExcMessage("invalid iterator"));
1065 
1066  return index_set->ranges[range_idx].end - 1;
1067 }
1068 
1069 
1070 
1073 {
1074  index_set = other.index_set;
1075  range_idx = other.range_idx;
1076  Assert(range_idx == numbers::invalid_dof_index || is_valid(),
1077  ExcMessage("invalid iterator"));
1078  return *this;
1079 }
1080 
1081 
1082 
1083 inline bool
1086 {
1087  Assert(index_set == other.index_set,
1088  ExcMessage(
1089  "Can not compare accessors pointing to different IndexSets"));
1090  return range_idx == other.range_idx;
1091 }
1092 
1093 
1094 
1095 inline bool
1098 {
1099  Assert(index_set == other.index_set,
1100  ExcMessage(
1101  "Can not compare accessors pointing to different IndexSets"));
1102  return range_idx < other.range_idx;
1103 }
1104 
1105 
1106 
1107 inline void
1109 {
1110  Assert(
1111  is_valid(),
1112  ExcMessage(
1113  "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1114  ++range_idx;
1115 
1116  // set ourselves to invalid if we walk off the end
1117  if (range_idx >= index_set->ranges.size())
1118  range_idx = numbers::invalid_dof_index;
1119 }
1120 
1121 
1122 /* IntervalIterator */
1123 
1125  const IndexSet * idxset,
1126  const IndexSet::size_type range_idx)
1127  : accessor(idxset, range_idx)
1128 {}
1129 
1130 
1131 
1133  : accessor(nullptr)
1134 {}
1135 
1136 
1137 
1139  : accessor(idxset)
1140 {}
1141 
1142 
1143 
1146 {
1147  accessor.advance();
1148  return *this;
1149 }
1150 
1151 
1152 
1155 {
1156  const IndexSet::IntervalIterator iter = *this;
1157  accessor.advance();
1158  return iter;
1159 }
1160 
1161 
1162 
1164  operator*() const
1165 {
1166  return accessor;
1167 }
1168 
1169 
1170 
1173 {
1174  return &accessor;
1175 }
1176 
1177 
1178 
1179 inline bool
1182 {
1183  return accessor == other.accessor;
1184 }
1185 
1186 
1187 
1188 inline bool
1191 {
1192  return !(*this == other);
1193 }
1194 
1195 
1196 
1197 inline bool
1200 {
1201  return accessor < other.accessor;
1202 }
1203 
1204 
1205 
1206 inline int
1209 {
1210  Assert(accessor.index_set == other.accessor.index_set,
1211  ExcMessage(
1212  "Can not compare iterators belonging to different IndexSets"));
1213 
1214  const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index) ?
1215  accessor.index_set->ranges.size() :
1216  accessor.range_idx;
1217  const size_type rhs =
1219  accessor.index_set->ranges.size() :
1220  other.accessor.range_idx;
1221 
1222  if (lhs > rhs)
1223  return static_cast<int>(lhs - rhs);
1224  else
1225  return -static_cast<int>(rhs - lhs);
1226 }
1227 
1228 
1229 
1230 /* ElementIterator */
1231 
1233  const IndexSet * idxset,
1234  const IndexSet::size_type range_idx,
1235  const IndexSet::size_type index)
1236  : index_set(idxset)
1237  , range_idx(range_idx)
1238  , idx(index)
1239 {
1240  Assert(range_idx < index_set->ranges.size(),
1241  ExcMessage(
1242  "Invalid range index for IndexSet::ElementIterator constructor."));
1243  Assert(
1244  idx >= index_set->ranges[range_idx].begin &&
1245  idx < index_set->ranges[range_idx].end,
1247  "Invalid index argument for IndexSet::ElementIterator constructor."));
1248 }
1249 
1250 
1251 
1253  : index_set(idxset)
1254  , range_idx(numbers::invalid_dof_index)
1255  , idx(numbers::invalid_dof_index)
1256 {}
1257 
1258 
1259 
1260 inline bool
1262 {
1263  Assert((range_idx == numbers::invalid_dof_index &&
1264  idx == numbers::invalid_dof_index) ||
1265  (range_idx < index_set->ranges.size() &&
1266  idx < index_set->ranges[range_idx].end),
1267  ExcInternalError("Invalid ElementIterator state."));
1268 
1269  return (range_idx < index_set->ranges.size() &&
1270  idx < index_set->ranges[range_idx].end);
1271 }
1272 
1273 
1274 
1276 {
1277  Assert(
1278  is_valid(),
1279  ExcMessage(
1280  "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1281  return idx;
1282 }
1283 
1284 
1285 
1286 inline bool
1289 {
1290  Assert(index_set == other.index_set,
1291  ExcMessage(
1292  "Can not compare iterators belonging to different IndexSets"));
1293  return range_idx == other.range_idx && idx == other.idx;
1294 }
1295 
1296 
1297 
1298 inline void
1300 {
1301  Assert(
1302  is_valid(),
1303  ExcMessage(
1304  "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1305  if (idx < index_set->ranges[range_idx].end)
1306  ++idx;
1307  // end of this range?
1308  if (idx == index_set->ranges[range_idx].end)
1309  {
1310  // point to first element in next interval if possible
1311  if (range_idx < index_set->ranges.size() - 1)
1312  {
1313  ++range_idx;
1314  idx = index_set->ranges[range_idx].begin;
1315  }
1316  else
1317  {
1318  // we just fell off the end, set to invalid:
1319  range_idx = numbers::invalid_dof_index;
1321  }
1322  }
1323 }
1324 
1325 
1326 
1329 {
1330  advance();
1331  return *this;
1332 }
1333 
1334 
1335 
1338 {
1339  const IndexSet::ElementIterator it = *this;
1340  advance();
1341  return it;
1342 }
1343 
1344 
1345 
1346 inline bool
1349 {
1350  return !(*this == other);
1351 }
1352 
1353 
1354 
1355 inline bool
1358 {
1359  Assert(index_set == other.index_set,
1360  ExcMessage(
1361  "Can not compare iterators belonging to different IndexSets"));
1362  return range_idx < other.range_idx ||
1363  (range_idx == other.range_idx && idx < other.idx);
1364 }
1365 
1366 
1367 
1368 inline std::ptrdiff_t
1371 {
1372  Assert(index_set == other.index_set,
1373  ExcMessage(
1374  "Can not compare iterators belonging to different IndexSets"));
1375  if (*this == other)
1376  return 0;
1377  if (!(*this < other))
1378  return -(other - *this);
1379 
1380  // only other can be equal to end() because of the checks above.
1381  Assert(is_valid(), ExcInternalError());
1382 
1383  // Note: we now compute how far advance *this in "*this < other" to get other,
1384  // so we need to return -c at the end.
1385 
1386  // first finish the current range:
1387  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1388 
1389  // now walk in steps of ranges (need to start one behind our current one):
1390  for (size_type range = range_idx + 1;
1391  range < index_set->ranges.size() && range <= other.range_idx;
1392  ++range)
1393  c += index_set->ranges[range].end - index_set->ranges[range].begin;
1394 
1395  Assert(
1396  other.range_idx < index_set->ranges.size() ||
1398  ExcMessage(
1399  "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1400 
1401  // We might have walked too far because we went until the end of
1402  // other.range_idx, so walk backwards to other.idx:
1404  c -= index_set->ranges[other.range_idx].end - other.idx;
1405 
1406  return -c;
1407 }
1408 
1409 
1410 /* Range */
1411 
1413  : begin(numbers::invalid_dof_index)
1414  , end(numbers::invalid_dof_index)
1415  , nth_index_in_set(numbers::invalid_dof_index)
1416 {}
1417 
1418 
1419 
1420 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1421  : begin(i1)
1422  , end(i2)
1423  , nth_index_in_set(numbers::invalid_dof_index)
1424 {}
1425 
1426 
1427 
1428 /* IndexSet itself */
1429 
1431  : is_compressed(true)
1432  , index_space_size(0)
1433  , largest_range(numbers::invalid_unsigned_int)
1434 {}
1435 
1436 
1437 
1438 inline IndexSet::IndexSet(const size_type size)
1439  : is_compressed(true)
1440  , index_space_size(size)
1441  , largest_range(numbers::invalid_unsigned_int)
1442 {}
1443 
1444 
1445 
1446 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1447  : ranges(std::move(is.ranges))
1448  , is_compressed(is.is_compressed)
1449  , index_space_size(is.index_space_size)
1450  , largest_range(is.largest_range)
1451 {
1452  is.ranges.clear();
1453  is.is_compressed = true;
1454  is.index_space_size = 0;
1455  is.largest_range = numbers::invalid_unsigned_int;
1456 
1457  compress();
1458 }
1459 
1460 
1461 
1462 inline IndexSet &
1464 {
1465  ranges = std::move(is.ranges);
1466  is_compressed = is.is_compressed;
1467  index_space_size = is.index_space_size;
1468  largest_range = is.largest_range;
1469 
1470  is.ranges.clear();
1471  is.is_compressed = true;
1472  is.index_space_size = 0;
1473  is.largest_range = numbers::invalid_unsigned_int;
1474 
1475  compress();
1476 
1477  return *this;
1478 }
1479 
1480 
1481 
1484 {
1485  compress();
1486  if (ranges.size() > 0)
1487  return {this, 0, ranges[0].begin};
1488  else
1489  return end();
1490 }
1491 
1492 
1493 
1495 IndexSet::at(const size_type global_index) const
1496 {
1497  compress();
1498  Assert(global_index < size(),
1499  ExcIndexRangeType<size_type>(global_index, 0, size()));
1500 
1501  if (ranges.empty())
1502  return end();
1503 
1504  std::vector<Range>::const_iterator main_range =
1505  ranges.begin() + largest_range;
1506 
1507  Range r(global_index, global_index + 1);
1508  // This optimization makes the bounds for lower_bound smaller by checking
1509  // the largest range first.
1510  std::vector<Range>::const_iterator range_begin, range_end;
1511  if (global_index < main_range->begin)
1512  {
1513  range_begin = ranges.begin();
1514  range_end = main_range;
1515  }
1516  else
1517  {
1518  range_begin = main_range;
1519  range_end = ranges.end();
1520  }
1521 
1522  // This will give us the first range p=[a,b[ with b>=global_index using
1523  // a binary search
1524  const std::vector<Range>::const_iterator p =
1525  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1526 
1527  // We couldn't find a range, which means we have no range that contains
1528  // global_index and also no range behind it, meaning we need to return end().
1529  if (p == ranges.end())
1530  return end();
1531 
1532  // Finally, we can have two cases: Either global_index is not in [a,b[,
1533  // which means we need to return an iterator to a because global_index, ...,
1534  // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
1535  // [a,b[ and we will return an iterator pointing directly at global_index
1536  // (else branch).
1537  if (global_index < p->begin)
1538  return {this, static_cast<size_type>(p - ranges.begin()), p->begin};
1539  else
1540  return {this, static_cast<size_type>(p - ranges.begin()), global_index};
1541 }
1542 
1543 
1544 
1547 {
1548  compress();
1549  return IndexSet::ElementIterator(this);
1550 }
1551 
1552 
1553 
1556 {
1557  compress();
1558  if (ranges.size() > 0)
1559  return IndexSet::IntervalIterator(this, 0);
1560  else
1561  return end_intervals();
1562 }
1563 
1564 
1565 
1568 {
1569  compress();
1570  return IndexSet::IntervalIterator(this);
1571 }
1572 
1573 
1574 
1575 inline void
1577 {
1578  // reset so that there are no indices in the set any more; however,
1579  // as documented, the index set retains its size
1580  ranges.clear();
1581  is_compressed = true;
1583 }
1584 
1585 
1586 
1587 inline void
1589 {
1590  Assert(ranges.empty(),
1591  ExcMessage("This function can only be called if the current "
1592  "object does not yet contain any elements."));
1593  index_space_size = sz;
1594  is_compressed = true;
1595 }
1596 
1597 
1598 
1599 inline IndexSet::size_type
1601 {
1602  return index_space_size;
1603 }
1604 
1605 
1606 
1607 inline void
1609 {
1610  if (is_compressed == true)
1611  return;
1612 
1613  do_compress();
1614 }
1615 
1616 
1617 
1618 inline void
1620 {
1621  Assert(index < index_space_size,
1622  ExcIndexRangeType<size_type>(index, 0, index_space_size));
1623 
1624  const Range new_range(index, index + 1);
1625  if (ranges.size() == 0 || index > ranges.back().end)
1626  ranges.push_back(new_range);
1627  else if (index == ranges.back().end)
1628  ranges.back().end++;
1629  else
1630  ranges.insert(Utilities::lower_bound(ranges.begin(),
1631  ranges.end(),
1632  new_range),
1633  new_range);
1634  is_compressed = false;
1635 }
1636 
1637 
1638 
1639 template <typename ForwardIterator>
1640 inline void
1641 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1642 {
1643  // insert each element of the range. if some of them happen to be
1644  // consecutive, merge them to a range
1645  for (ForwardIterator p = begin; p != end;)
1646  {
1647  const size_type begin_index = *p;
1648  size_type end_index = begin_index + 1;
1649  ForwardIterator q = p;
1650  ++q;
1651  while ((q != end) && (*q == end_index))
1652  {
1653  ++end_index;
1654  ++q;
1655  }
1656 
1657  add_range(begin_index, end_index);
1658  p = q;
1659  }
1660 }
1661 
1662 
1663 
1664 inline bool
1666 {
1667  if (ranges.empty() == false)
1668  {
1669  compress();
1670 
1671  // fast check whether the index is in the largest range
1673  if (index >= ranges[largest_range].begin &&
1674  index < ranges[largest_range].end)
1675  return true;
1676 
1677  // get the element after which we would have to insert a range that
1678  // consists of all elements from this element to the end of the index
1679  // range plus one. after this call we know that if p!=end() then
1680  // p->begin<=index unless there is no such range at all
1681  //
1682  // if the searched for element is an element of this range, then we're
1683  // done. otherwise, the element can't be in one of the following ranges
1684  // because otherwise p would be a different iterator
1685  //
1686  // since we already know the position relative to the largest range (we
1687  // called compress!), we can perform the binary search on ranges with
1688  // lower/higher number compared to the largest range
1689  std::vector<Range>::const_iterator p = std::upper_bound(
1690  ranges.begin() +
1691  (index < ranges[largest_range].begin ? 0 : largest_range + 1),
1692  index < ranges[largest_range].begin ? ranges.begin() + largest_range :
1693  ranges.end(),
1694  Range(index, size() + 1));
1695 
1696  if (p == ranges.begin())
1697  return ((index >= p->begin) && (index < p->end));
1698 
1699  Assert((p == ranges.end()) || (p->begin > index), ExcInternalError());
1700 
1701  // now move to that previous range
1702  --p;
1703  Assert(p->begin <= index, ExcInternalError());
1704 
1705  return (p->end > index);
1706  }
1707 
1708  // didn't find this index, so it's not in the set
1709  return false;
1710 }
1711 
1712 
1713 
1714 inline bool
1716 {
1717  compress();
1718  return (ranges.size() <= 1);
1719 }
1720 
1721 
1722 
1723 inline bool
1725 {
1726  return ranges.empty();
1727 }
1728 
1729 
1730 
1731 inline IndexSet::size_type
1733 {
1734  // make sure we have non-overlapping ranges
1735  compress();
1736 
1737  size_type v = 0;
1738  if (!ranges.empty())
1739  {
1740  Range &r = ranges.back();
1741  v = r.nth_index_in_set + r.end - r.begin;
1742  }
1743 
1744 #ifdef DEBUG
1745  size_type s = 0;
1746  for (const auto &range : ranges)
1747  s += (range.end - range.begin);
1748  Assert(s == v, ExcInternalError());
1749 #endif
1750 
1751  return v;
1752 }
1753 
1754 
1755 
1756 inline unsigned int
1758 {
1759  compress();
1760  return ranges.size();
1761 }
1762 
1763 
1764 
1765 inline unsigned int
1767 {
1768  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1769 
1770  compress();
1771  const std::vector<Range>::const_iterator main_range =
1772  ranges.begin() + largest_range;
1773 
1774  return main_range->nth_index_in_set;
1775 }
1776 
1777 
1778 
1779 inline IndexSet::size_type
1780 IndexSet::nth_index_in_set(const unsigned int n) const
1781 {
1782  Assert(n < n_elements(), ExcIndexRangeType<size_type>(n, 0, n_elements()));
1783 
1784  compress();
1785 
1786  // first check whether the index is in the largest range
1788  std::vector<Range>::const_iterator main_range =
1789  ranges.begin() + largest_range;
1790  if (n >= main_range->nth_index_in_set &&
1791  n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1792  return main_range->begin + (n - main_range->nth_index_in_set);
1793 
1794  // find out which chunk the local index n belongs to by using a binary
1795  // search. the comparator is based on the end of the ranges. Use the
1796  // position relative to main_range to subdivide the ranges
1797  Range r(n, n + 1);
1798  r.nth_index_in_set = n;
1799  std::vector<Range>::const_iterator range_begin, range_end;
1800  if (n < main_range->nth_index_in_set)
1801  {
1802  range_begin = ranges.begin();
1803  range_end = main_range;
1804  }
1805  else
1806  {
1807  range_begin = main_range + 1;
1808  range_end = ranges.end();
1809  }
1810 
1811  const std::vector<Range>::const_iterator p =
1812  Utilities::lower_bound(range_begin, range_end, r, Range::nth_index_compare);
1813 
1814  Assert(p != ranges.end(), ExcInternalError());
1815  return p->begin + (n - p->nth_index_in_set);
1816 }
1817 
1818 
1819 
1820 inline IndexSet::size_type
1822 {
1823  // to make this call thread-safe, compress() must not be called through this
1824  // function
1825  Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1826  Assert(n < size(), ExcIndexRangeType<size_type>(n, 0, size()));
1827 
1828  // return immediately if the index set is empty
1829  if (is_empty())
1831 
1832  // check whether the index is in the largest range. use the result to
1833  // perform a one-sided binary search afterward
1835  std::vector<Range>::const_iterator main_range =
1836  ranges.begin() + largest_range;
1837  if (n >= main_range->begin && n < main_range->end)
1838  return (n - main_range->begin) + main_range->nth_index_in_set;
1839 
1840  Range r(n, n);
1841  std::vector<Range>::const_iterator range_begin, range_end;
1842  if (n < main_range->begin)
1843  {
1844  range_begin = ranges.begin();
1845  range_end = main_range;
1846  }
1847  else
1848  {
1849  range_begin = main_range + 1;
1850  range_end = ranges.end();
1851  }
1852 
1853  std::vector<Range>::const_iterator p =
1854  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1855 
1856  // if n is not in this set
1857  if (p == range_end || p->end == n || p->begin > n)
1859 
1860  Assert(p != ranges.end(), ExcInternalError());
1861  Assert(p->begin <= n, ExcInternalError());
1862  Assert(n < p->end, ExcInternalError());
1863  return (n - p->begin) + p->nth_index_in_set;
1864 }
1865 
1866 
1867 
1868 inline bool
1870 {
1871  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1872 
1873  compress();
1874  is.compress();
1875 
1876  return ranges == is.ranges;
1877 }
1878 
1879 
1880 
1881 inline bool
1883 {
1884  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1885 
1886  compress();
1887  is.compress();
1888 
1889  return ranges != is.ranges;
1890 }
1891 
1892 
1893 
1894 template <typename Vector>
1895 void
1896 IndexSet::fill_binary_vector(Vector &vector) const
1897 {
1898  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
1899 
1900  compress();
1901  // first fill all elements of the vector with zeroes.
1902  std::fill(vector.begin(), vector.end(), 0);
1903 
1904  // then write ones into the elements whose indices are contained in the
1905  // index set
1906  for (const auto &range : ranges)
1907  for (size_type i = range.begin; i < range.end; ++i)
1908  vector[i] = 1;
1909 }
1910 
1911 
1912 
1913 template <class StreamType>
1914 inline void
1915 IndexSet::print(StreamType &out) const
1916 {
1917  compress();
1918  out << "{";
1919  std::vector<Range>::const_iterator p;
1920  for (p = ranges.begin(); p != ranges.end(); ++p)
1921  {
1922  if (p->end - p->begin == 1)
1923  out << p->begin;
1924  else
1925  out << "[" << p->begin << "," << p->end - 1 << "]";
1926 
1927  if (p != --ranges.end())
1928  out << ", ";
1929  }
1930  out << "}" << std::endl;
1931 }
1932 
1933 
1934 
1935 template <class Archive>
1936 inline void
1937 IndexSet::Range::serialize(Archive &ar, const unsigned int)
1938 {
1939  ar &begin &end &nth_index_in_set;
1940 }
1941 
1942 
1943 
1944 template <class Archive>
1945 inline void
1946 IndexSet::serialize(Archive &ar, const unsigned int)
1947 {
1949 }
1950 
1951 DEAL_II_NAMESPACE_CLOSE
1952 
1953 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1067
void fill_binary_vector(VectorType &vector) const
static const unsigned int invalid_unsigned_int
Definition: types.h:173
ElementIterator end() const
Definition: index_set.h:1546
void add_index(const size_type index)
Definition: index_set.h:1619
const IntervalAccessor * operator->() const
Definition: index_set.h:1172
const IntervalAccessor & operator*() const
Definition: index_set.h:1164
size_type nth_index_in_set(const unsigned int local_index) const
Definition: index_set.h:1780
const IndexSet * index_set
Definition: index_set.h:773
IndexSet & operator=(const IndexSet &)=default
void block_read(std::istream &in)
Definition: index_set.cc:485
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1370
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1641
unsigned int largest_range_starting_index() const
Definition: index_set.h:1766
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:1232
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1199
std::forward_iterator_tag iterator_category
Definition: index_set.h:666
std::size_t memory_consumption() const
Definition: index_set.cc:727
LinearAlgebra::distributed::Vector< Number > Vector
size_type n_elements() const
Definition: index_set.h:1022
size_type index_space_size
Definition: index_set.h:933
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:664
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:594
size_type size() const
Definition: index_set.h:1600
void print(StreamType &out) const
Definition: index_set.h:1915
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1937
IndexSet complete_index_set(const unsigned int N)
Definition: index_set.h:978
void clear()
Definition: index_set.h:1576
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1946
static ::ExceptionBase & ExcMessage(std::string arg1)
void block_write(std::ostream &out) const
Definition: index_set.cc:471
std::vector< Range > ranges
Definition: index_set.h:917
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:267
#define Assert(cond, exc)
Definition: exceptions.h:1407
std::forward_iterator_tag iterator_category
Definition: index_set.h:757
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
bool operator<(const ElementIterator &) const
Definition: index_set.h:1357
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1821
bool operator==(const IndexSet &is) const
Definition: index_set.h:1869
signed int value_type
Definition: index_set.h:104
IntervalAccessor accessor
Definition: index_set.h:676
void do_compress() const
Definition: index_set.cc:127
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1072
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1181
size_type last() const
Definition: index_set.h:1062
const IndexSet * index_set
Definition: index_set.h:567
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1085
unsigned int n_intervals() const
Definition: index_set.h:1757
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1208
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1882
size_type pop_back()
Definition: index_set.cc:339
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:505
void add_range(const size_type begin, const size_type end)
Definition: index_set.cc:98
IndexSet operator &(const IndexSet &is) const
IntervalIterator & operator=(const IntervalIterator &other)=default
void set_size(const size_type size)
Definition: index_set.h:1588
ElementIterator & operator++()
Definition: index_set.h:1328
unsigned int global_dof_index
Definition: types.h:89
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:991
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:238
void compress() const
Definition: index_set.h:1608
ElementIterator end() const
Definition: index_set.h:1048
IntervalIterator begin_intervals() const
Definition: index_set.h:1555
bool is_contiguous() const
Definition: index_set.h:1715
Threads::Mutex compress_mutex
Definition: index_set.h:950
size_type pop_front()
Definition: index_set.cc:357
IntervalIterator end_intervals() const
Definition: index_set.h:1567
bool is_compressed
Definition: index_set.h:927
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1097
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1348
ElementIterator begin() const
Definition: index_set.h:1039
void write(std::ostream &out) const
Definition: index_set.cc:434
bool is_element(const size_type index) const
Definition: index_set.h:1665
bool operator==(const ElementIterator &) const
Definition: index_set.h:1288
ElementIterator at(const size_type global_index) const
Definition: index_set.h:1495
const types::global_dof_index invalid_dof_index
Definition: types.h:188
ElementIterator begin() const
Definition: index_set.h:1483
void read(std::istream &in)
Definition: index_set.cc:449
size_type n_elements() const
Definition: index_set.h:1732
bool is_empty() const
Definition: index_set.h:1724
size_type largest_range
Definition: index_set.h:944
IntervalIterator & operator++()
Definition: index_set.h:1145
static ::ExceptionBase & ExcInternalError()
size_type operator*() const
Definition: index_set.h:1275
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1190