Reference documentation for deal.II version Git 41b7f7ec9f 2019-11-14 19:01:26 -0600
\(\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 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
273  index_within_set(const size_type global_index) const;
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 
352  void
353  subtract_set(const IndexSet &other);
354 
359  size_type
360  pop_back();
361 
366  size_type
367  pop_front();
368 
372  void
373  fill_index_vector(std::vector<size_type> &indices) const;
374 
387  template <typename VectorType>
388  void
389  fill_binary_vector(VectorType &vector) const;
390 
395  template <class StreamType>
396  void
397  print(StreamType &out) const;
398 
403  void
404  write(std::ostream &out) const;
405 
410  void
411  read(std::istream &in);
412 
417  void
418  block_write(std::ostream &out) const;
419 
424  void
425  block_read(std::istream &in);
426 
427 #ifdef DEAL_II_WITH_TRILINOS
428 
456  Epetra_Map
457  make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
458  const bool overlapping = false) const;
459 
460 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
461  Tpetra::Map<int, types::global_dof_index>
462  make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
463  const bool overlapping = false) const;
464 # endif
465 #endif
466 
467 
472  std::size_t
473  memory_consumption() const;
474 
476  size_type,
477  << "The global index " << arg1
478  << " is not an element of this set.");
479 
484  template <class Archive>
485  void
486  serialize(Archive &ar, const unsigned int version);
487 
488 
500  {
501  public:
506  IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
507 
511  explicit IntervalAccessor(const IndexSet *idxset);
512 
516  size_type
517  n_elements() const;
518 
522  bool
523  is_valid() const;
524 
529  begin() const;
530 
536  end() const;
537 
541  size_type
542  last() const;
543 
544  private:
548  IntervalAccessor(const IntervalAccessor &other);
553  operator=(const IntervalAccessor &other);
554 
558  bool
559  operator==(const IntervalAccessor &other) const;
563  bool
564  operator<(const IntervalAccessor &other) const;
569  void
570  advance();
575 
581 
582  friend class IntervalIterator;
583  };
584 
590  {
591  public:
596  IntervalIterator(const IndexSet *idxset, const size_type range_idx);
597 
601  explicit IntervalIterator(const IndexSet *idxset);
602 
607 
611  IntervalIterator(const IntervalIterator &other) = default;
612 
617  operator=(const IntervalIterator &other) = default;
618 
623  operator++();
624 
629  operator++(int);
630 
634  const IntervalAccessor &operator*() const;
635 
639  const IntervalAccessor *operator->() const;
640 
644  bool
645  operator==(const IntervalIterator &) const;
646 
650  bool
651  operator!=(const IntervalIterator &) const;
652 
656  bool
657  operator<(const IntervalIterator &) const;
658 
665  int
666  operator-(const IntervalIterator &p) const;
667 
673  using iterator_category = std::forward_iterator_tag;
675  using difference_type = std::ptrdiff_t;
676  using pointer = IntervalAccessor *;
677  using reference = IntervalAccessor &;
678 
679  private:
684  };
685 
691  {
692  public:
697  ElementIterator(const IndexSet *idxset,
698  const size_type range_idx,
699  const size_type index);
700 
704  explicit ElementIterator(const IndexSet *idxset);
705 
710  size_type operator*() const;
711 
715  bool
716  is_valid() const;
717 
722  operator++();
723 
728  operator++(int);
729 
733  bool
734  operator==(const ElementIterator &) const;
735 
739  bool
740  operator!=(const ElementIterator &) const;
741 
745  bool
746  operator<(const ElementIterator &) const;
747 
756  std::ptrdiff_t
757  operator-(const ElementIterator &p) const;
758 
764  using iterator_category = std::forward_iterator_tag;
765  using value_type = size_type;
766  using difference_type = std::ptrdiff_t;
767  using pointer = size_type *;
768  using reference = size_type &;
769 
770  private:
774  void
775  advance();
776 
789  };
790 
796  begin() const;
797 
813  at(const size_type global_index) const;
814 
820  end() const;
821 
826  begin_intervals() const;
827 
833  end_intervals() const;
834 
839 private:
848  struct Range
849  {
851  size_type end;
852 
854 
863  Range();
864 
872  Range(const size_type i1, const size_type i2);
873 
874  friend inline bool
875  operator<(const Range &range_1, const Range &range_2)
876  {
877  return (
878  (range_1.begin < range_2.begin) ||
879  ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
880  }
881 
882  static bool
883  end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
884  {
885  return x.end < y.end;
886  }
887 
888  static bool
889  nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
890  {
891  return (x.nth_index_in_set + (x.end - x.begin) <
892  y.nth_index_in_set + (y.end - y.begin));
893  }
894 
895  friend inline bool
896  operator==(const Range &range_1, const Range &range_2)
897  {
898  return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
899  }
900 
901  static std::size_t
903  {
904  return sizeof(Range);
905  }
906 
911  template <class Archive>
912  void
913  serialize(Archive &ar, const unsigned int version);
914  };
915 
924  mutable std::vector<Range> ranges;
925 
934  mutable bool is_compressed;
935 
941 
952 
958 
962  void
963  do_compress() const;
964 };
965 
966 
984 inline IndexSet
986 {
987  IndexSet is(N);
988  is.add_range(0, N);
989  is.compress();
990  return is;
991 }
992 
993 /* ------------------ inline functions ------------------ */
994 
995 
996 /* IntervalAccessor */
997 
999  const IndexSet * idxset,
1001  : index_set(idxset)
1002  , range_idx(range_idx)
1003 {
1004  Assert(range_idx < idxset->n_intervals(),
1005  ExcInternalError("Invalid range index"));
1006 }
1007 
1008 
1009 
1011  : index_set(idxset)
1012  , range_idx(numbers::invalid_dof_index)
1013 {}
1014 
1015 
1016 
1018  const IndexSet::IntervalAccessor &other)
1019  : index_set(other.index_set)
1020  , range_idx(other.range_idx)
1021 {
1023  ExcMessage("invalid iterator"));
1024 }
1025 
1026 
1027 
1028 inline IndexSet::size_type
1030 {
1031  Assert(is_valid(), ExcMessage("invalid iterator"));
1032  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1033 }
1034 
1035 
1036 
1037 inline bool
1039 {
1040  return index_set != nullptr && range_idx < index_set->n_intervals();
1041 }
1042 
1043 
1044 
1047 {
1048  Assert(is_valid(), ExcMessage("invalid iterator"));
1049  return {index_set, range_idx, index_set->ranges[range_idx].begin};
1050 }
1051 
1052 
1053 
1056 {
1057  Assert(is_valid(), ExcMessage("invalid iterator"));
1058 
1059  // point to first index in next interval unless we are the last interval.
1060  if (range_idx < index_set->ranges.size() - 1)
1061  return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1062  else
1063  return index_set->end();
1064 }
1065 
1066 
1067 
1068 inline IndexSet::size_type
1070 {
1071  Assert(is_valid(), ExcMessage("invalid iterator"));
1072 
1073  return index_set->ranges[range_idx].end - 1;
1074 }
1075 
1076 
1077 
1080 {
1081  index_set = other.index_set;
1082  range_idx = other.range_idx;
1084  ExcMessage("invalid iterator"));
1085  return *this;
1086 }
1087 
1088 
1089 
1090 inline bool
1093 {
1094  Assert(index_set == other.index_set,
1095  ExcMessage(
1096  "Can not compare accessors pointing to different IndexSets"));
1097  return range_idx == other.range_idx;
1098 }
1099 
1100 
1101 
1102 inline bool
1105 {
1106  Assert(index_set == other.index_set,
1107  ExcMessage(
1108  "Can not compare accessors pointing to different IndexSets"));
1109  return range_idx < other.range_idx;
1110 }
1111 
1112 
1113 
1114 inline void
1116 {
1117  Assert(
1118  is_valid(),
1119  ExcMessage(
1120  "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1121  ++range_idx;
1122 
1123  // set ourselves to invalid if we walk off the end
1124  if (range_idx >= index_set->ranges.size())
1126 }
1127 
1128 
1129 /* IntervalIterator */
1130 
1132  const IndexSet * idxset,
1134  : accessor(idxset, range_idx)
1135 {}
1136 
1137 
1138 
1140  : accessor(nullptr)
1141 {}
1142 
1143 
1144 
1146  : accessor(idxset)
1147 {}
1148 
1149 
1150 
1153 {
1154  accessor.advance();
1155  return *this;
1156 }
1157 
1158 
1159 
1162 {
1163  const IndexSet::IntervalIterator iter = *this;
1164  accessor.advance();
1165  return iter;
1166 }
1167 
1168 
1169 
1171  operator*() const
1172 {
1173  return accessor;
1174 }
1175 
1176 
1177 
1180 {
1181  return &accessor;
1182 }
1183 
1184 
1185 
1186 inline bool
1189 {
1190  return accessor == other.accessor;
1191 }
1192 
1193 
1194 
1195 inline bool
1198 {
1199  return !(*this == other);
1200 }
1201 
1202 
1203 
1204 inline bool
1207 {
1208  return accessor < other.accessor;
1209 }
1210 
1211 
1212 
1213 inline int
1216 {
1218  ExcMessage(
1219  "Can not compare iterators belonging to different IndexSets"));
1220 
1222  accessor.index_set->ranges.size() :
1224  const size_type rhs =
1226  accessor.index_set->ranges.size() :
1227  other.accessor.range_idx;
1228 
1229  if (lhs > rhs)
1230  return static_cast<int>(lhs - rhs);
1231  else
1232  return -static_cast<int>(rhs - lhs);
1233 }
1234 
1235 
1236 
1237 /* ElementIterator */
1238 
1240  const IndexSet * idxset,
1241  const IndexSet::size_type range_idx,
1242  const IndexSet::size_type index)
1243  : index_set(idxset)
1244  , range_idx(range_idx)
1245  , idx(index)
1246 {
1247  Assert(range_idx < index_set->ranges.size(),
1248  ExcMessage(
1249  "Invalid range index for IndexSet::ElementIterator constructor."));
1250  Assert(
1251  idx >= index_set->ranges[range_idx].begin &&
1252  idx < index_set->ranges[range_idx].end,
1254  "Invalid index argument for IndexSet::ElementIterator constructor."));
1255 }
1256 
1257 
1258 
1260  : index_set(idxset)
1261  , range_idx(numbers::invalid_dof_index)
1262  , idx(numbers::invalid_dof_index)
1263 {}
1264 
1265 
1266 
1267 inline bool
1269 {
1272  (range_idx < index_set->ranges.size() &&
1273  idx < index_set->ranges[range_idx].end),
1274  ExcInternalError("Invalid ElementIterator state."));
1275 
1276  return (range_idx < index_set->ranges.size() &&
1277  idx < index_set->ranges[range_idx].end);
1278 }
1279 
1280 
1281 
1283 {
1284  Assert(
1285  is_valid(),
1286  ExcMessage(
1287  "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1288  return idx;
1289 }
1290 
1291 
1292 
1293 inline bool
1296 {
1297  Assert(index_set == other.index_set,
1298  ExcMessage(
1299  "Can not compare iterators belonging to different IndexSets"));
1300  return range_idx == other.range_idx && idx == other.idx;
1301 }
1302 
1303 
1304 
1305 inline void
1307 {
1308  Assert(
1309  is_valid(),
1310  ExcMessage(
1311  "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1312  if (idx < index_set->ranges[range_idx].end)
1313  ++idx;
1314  // end of this range?
1315  if (idx == index_set->ranges[range_idx].end)
1316  {
1317  // point to first element in next interval if possible
1318  if (range_idx < index_set->ranges.size() - 1)
1319  {
1320  ++range_idx;
1321  idx = index_set->ranges[range_idx].begin;
1322  }
1323  else
1324  {
1325  // we just fell off the end, set to invalid:
1328  }
1329  }
1330 }
1331 
1332 
1333 
1336 {
1337  advance();
1338  return *this;
1339 }
1340 
1341 
1342 
1345 {
1346  const IndexSet::ElementIterator it = *this;
1347  advance();
1348  return it;
1349 }
1350 
1351 
1352 
1353 inline bool
1356 {
1357  return !(*this == other);
1358 }
1359 
1360 
1361 
1362 inline bool
1365 {
1366  Assert(index_set == other.index_set,
1367  ExcMessage(
1368  "Can not compare iterators belonging to different IndexSets"));
1369  return range_idx < other.range_idx ||
1370  (range_idx == other.range_idx && idx < other.idx);
1371 }
1372 
1373 
1374 
1375 inline std::ptrdiff_t
1378 {
1379  Assert(index_set == other.index_set,
1380  ExcMessage(
1381  "Can not compare iterators belonging to different IndexSets"));
1382  if (*this == other)
1383  return 0;
1384  if (!(*this < other))
1385  return -(other - *this);
1386 
1387  // only other can be equal to end() because of the checks above.
1389 
1390  // Note: we now compute how far advance *this in "*this < other" to get other,
1391  // so we need to return -c at the end.
1392 
1393  // first finish the current range:
1394  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1395 
1396  // now walk in steps of ranges (need to start one behind our current one):
1397  for (size_type range = range_idx + 1;
1398  range < index_set->ranges.size() && range <= other.range_idx;
1399  ++range)
1400  c += index_set->ranges[range].end - index_set->ranges[range].begin;
1401 
1402  Assert(
1403  other.range_idx < index_set->ranges.size() ||
1405  ExcMessage(
1406  "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1407 
1408  // We might have walked too far because we went until the end of
1409  // other.range_idx, so walk backwards to other.idx:
1411  c -= index_set->ranges[other.range_idx].end - other.idx;
1412 
1413  return -c;
1414 }
1415 
1416 
1417 /* Range */
1418 
1420  : begin(numbers::invalid_dof_index)
1421  , end(numbers::invalid_dof_index)
1422  , nth_index_in_set(numbers::invalid_dof_index)
1423 {}
1424 
1425 
1426 
1427 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1428  : begin(i1)
1429  , end(i2)
1430  , nth_index_in_set(numbers::invalid_dof_index)
1431 {}
1432 
1433 
1434 
1435 /* IndexSet itself */
1436 
1438  : is_compressed(true)
1439  , index_space_size(0)
1440  , largest_range(numbers::invalid_unsigned_int)
1441 {}
1442 
1443 
1444 
1446  : is_compressed(true)
1447  , index_space_size(size)
1448  , largest_range(numbers::invalid_unsigned_int)
1449 {}
1450 
1451 
1452 
1453 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1454  : ranges(std::move(is.ranges))
1455  , is_compressed(is.is_compressed)
1456  , index_space_size(is.index_space_size)
1457  , largest_range(is.largest_range)
1458 {
1459  is.ranges.clear();
1460  is.is_compressed = true;
1461  is.index_space_size = 0;
1462  is.largest_range = numbers::invalid_unsigned_int;
1463 
1464  compress();
1465 }
1466 
1467 
1468 
1469 inline IndexSet &
1471 {
1472  ranges = std::move(is.ranges);
1473  is_compressed = is.is_compressed;
1474  index_space_size = is.index_space_size;
1475  largest_range = is.largest_range;
1476 
1477  is.ranges.clear();
1478  is.is_compressed = true;
1479  is.index_space_size = 0;
1480  is.largest_range = numbers::invalid_unsigned_int;
1481 
1482  compress();
1483 
1484  return *this;
1485 }
1486 
1487 
1488 
1491 {
1492  compress();
1493  if (ranges.size() > 0)
1494  return {this, 0, ranges[0].begin};
1495  else
1496  return end();
1497 }
1498 
1499 
1500 
1502 IndexSet::at(const size_type global_index) const
1503 {
1504  compress();
1505  Assert(global_index < size(),
1506  ExcIndexRangeType<size_type>(global_index, 0, size()));
1507 
1508  if (ranges.empty())
1509  return end();
1510 
1511  std::vector<Range>::const_iterator main_range =
1512  ranges.begin() + largest_range;
1513 
1514  Range r(global_index, global_index + 1);
1515  // This optimization makes the bounds for lower_bound smaller by checking
1516  // the largest range first.
1517  std::vector<Range>::const_iterator range_begin, range_end;
1518  if (global_index < main_range->begin)
1519  {
1520  range_begin = ranges.begin();
1521  range_end = main_range;
1522  }
1523  else
1524  {
1525  range_begin = main_range;
1526  range_end = ranges.end();
1527  }
1528 
1529  // This will give us the first range p=[a,b[ with b>=global_index using
1530  // a binary search
1531  const std::vector<Range>::const_iterator p =
1532  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1533 
1534  // We couldn't find a range, which means we have no range that contains
1535  // global_index and also no range behind it, meaning we need to return end().
1536  if (p == ranges.end())
1537  return end();
1538 
1539  // Finally, we can have two cases: Either global_index is not in [a,b[,
1540  // which means we need to return an iterator to a because global_index, ...,
1541  // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
1542  // [a,b[ and we will return an iterator pointing directly at global_index
1543  // (else branch).
1544  if (global_index < p->begin)
1545  return {this, static_cast<size_type>(p - ranges.begin()), p->begin};
1546  else
1547  return {this, static_cast<size_type>(p - ranges.begin()), global_index};
1548 }
1549 
1550 
1551 
1554 {
1555  compress();
1556  return IndexSet::ElementIterator(this);
1557 }
1558 
1559 
1560 
1563 {
1564  compress();
1565  if (ranges.size() > 0)
1566  return IndexSet::IntervalIterator(this, 0);
1567  else
1568  return end_intervals();
1569 }
1570 
1571 
1572 
1575 {
1576  compress();
1577  return IndexSet::IntervalIterator(this);
1578 }
1579 
1580 
1581 
1582 inline void
1584 {
1585  // reset so that there are no indices in the set any more; however,
1586  // as documented, the index set retains its size
1587  ranges.clear();
1588  is_compressed = true;
1590 }
1591 
1592 
1593 
1594 inline void
1596 {
1597  Assert(ranges.empty(),
1598  ExcMessage("This function can only be called if the current "
1599  "object does not yet contain any elements."));
1600  index_space_size = sz;
1601  is_compressed = true;
1602 }
1603 
1604 
1605 
1606 inline IndexSet::size_type
1608 {
1609  return index_space_size;
1610 }
1611 
1612 
1613 
1614 inline void
1616 {
1617  if (is_compressed == true)
1618  return;
1619 
1620  do_compress();
1621 }
1622 
1623 
1624 
1625 inline void
1627 {
1628  Assert(index < index_space_size,
1629  ExcIndexRangeType<size_type>(index, 0, index_space_size));
1630 
1631  const Range new_range(index, index + 1);
1632  if (ranges.size() == 0 || index > ranges.back().end)
1633  ranges.push_back(new_range);
1634  else if (index == ranges.back().end)
1635  ranges.back().end++;
1636  else
1637  ranges.insert(Utilities::lower_bound(ranges.begin(),
1638  ranges.end(),
1639  new_range),
1640  new_range);
1641  is_compressed = false;
1642 }
1643 
1644 
1645 
1646 inline void
1648 {
1649  Assert((begin < index_space_size) ||
1650  ((begin == index_space_size) && (end == index_space_size)),
1651  ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1652  Assert(end <= index_space_size,
1653  ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1654  Assert(begin <= end, ExcIndexRangeType<size_type>(begin, 0, end));
1655 
1656  if (begin != end)
1657  {
1658  const Range new_range(begin, end);
1659 
1660  // the new index might be larger than the last index present in the
1661  // ranges. Then we can skip the binary search
1662  if (ranges.size() == 0 || begin > ranges.back().end)
1663  ranges.push_back(new_range);
1664  else
1665  ranges.insert(Utilities::lower_bound(ranges.begin(),
1666  ranges.end(),
1667  new_range),
1668  new_range);
1669  is_compressed = false;
1670  }
1671 }
1672 
1673 
1674 
1675 template <typename ForwardIterator>
1676 inline void
1677 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1678 {
1679  if (begin == end)
1680  return;
1681 
1682  // identify ranges in the given iterator range by checking whether some
1683  // indices happen to be consecutive. to avoid quadratic complexity when
1684  // calling add_range many times (as add_range() going into the middle of an
1685  // already existing range must shift entries around), we first collect a
1686  // vector of ranges.
1687  std::vector<std::pair<size_type, size_type>> tmp_ranges;
1688  bool ranges_are_sorted = true;
1689  for (ForwardIterator p = begin; p != end;)
1690  {
1691  const size_type begin_index = *p;
1692  size_type end_index = begin_index + 1;
1693  ForwardIterator q = p;
1694  ++q;
1695  while ((q != end) && (*q == end_index))
1696  {
1697  ++end_index;
1698  ++q;
1699  }
1700 
1701  tmp_ranges.emplace_back(begin_index, end_index);
1702  p = q;
1703 
1704  // if the starting index of the next go-around of the for loop is less
1705  // than the end index of the one just identified, then we will have at
1706  // least one pair of ranges that are not sorted, and consequently the
1707  // whole collection of ranges is not sorted.
1708  if (p != end && *p < end_index)
1709  ranges_are_sorted = false;
1710  }
1711 
1712  if (!ranges_are_sorted)
1713  std::sort(tmp_ranges.begin(), tmp_ranges.end());
1714 
1715  // if we have many ranges, we first construct a temporary index set (where
1716  // we add ranges in a consecutive way, so fast), otherwise, we work with
1717  // add_range(). the number 9 is chosen heuristically given the fact that
1718  // there are typically up to 8 independent ranges when adding the degrees of
1719  // freedom on a 3D cell or 9 when adding degrees of freedom of faces. if
1720  // doing cell-by-cell additions, we want to avoid repeated calls to
1721  // IndexSet::compress() which gets called upon merging two index sets, so we
1722  // want to be in the other branch then.
1723  if (tmp_ranges.size() > 9)
1724  {
1725  IndexSet tmp_set(size());
1726  tmp_set.ranges.reserve(tmp_ranges.size());
1727  for (const auto &i : tmp_ranges)
1728  tmp_set.add_range(i.first, i.second);
1729  this->add_indices(tmp_set);
1730  }
1731  else
1732  for (const auto &i : tmp_ranges)
1733  add_range(i.first, i.second);
1734 }
1735 
1736 
1737 
1738 inline bool
1740 {
1741  if (ranges.empty() == false)
1742  {
1743  compress();
1744 
1745  // fast check whether the index is in the largest range
1747  if (index >= ranges[largest_range].begin &&
1748  index < ranges[largest_range].end)
1749  return true;
1750 
1751  // get the element after which we would have to insert a range that
1752  // consists of all elements from this element to the end of the index
1753  // range plus one. after this call we know that if p!=end() then
1754  // p->begin<=index unless there is no such range at all
1755  //
1756  // if the searched for element is an element of this range, then we're
1757  // done. otherwise, the element can't be in one of the following ranges
1758  // because otherwise p would be a different iterator
1759  //
1760  // since we already know the position relative to the largest range (we
1761  // called compress!), we can perform the binary search on ranges with
1762  // lower/higher number compared to the largest range
1763  std::vector<Range>::const_iterator p = std::upper_bound(
1764  ranges.begin() +
1765  (index < ranges[largest_range].begin ? 0 : largest_range + 1),
1766  index < ranges[largest_range].begin ? ranges.begin() + largest_range :
1767  ranges.end(),
1768  Range(index, size() + 1));
1769 
1770  if (p == ranges.begin())
1771  return ((index >= p->begin) && (index < p->end));
1772 
1773  Assert((p == ranges.end()) || (p->begin > index), ExcInternalError());
1774 
1775  // now move to that previous range
1776  --p;
1777  Assert(p->begin <= index, ExcInternalError());
1778 
1779  return (p->end > index);
1780  }
1781 
1782  // didn't find this index, so it's not in the set
1783  return false;
1784 }
1785 
1786 
1787 
1788 inline bool
1790 {
1791  compress();
1792  return (ranges.size() <= 1);
1793 }
1794 
1795 
1796 
1797 inline bool
1799 {
1800  return ranges.empty();
1801 }
1802 
1803 
1804 
1805 inline IndexSet::size_type
1807 {
1808  // make sure we have non-overlapping ranges
1809  compress();
1810 
1811  size_type v = 0;
1812  if (!ranges.empty())
1813  {
1814  Range &r = ranges.back();
1815  v = r.nth_index_in_set + r.end - r.begin;
1816  }
1817 
1818 #ifdef DEBUG
1819  size_type s = 0;
1820  for (const auto &range : ranges)
1821  s += (range.end - range.begin);
1822  Assert(s == v, ExcInternalError());
1823 #endif
1824 
1825  return v;
1826 }
1827 
1828 
1829 
1830 inline unsigned int
1832 {
1833  compress();
1834  return ranges.size();
1835 }
1836 
1837 
1838 
1839 inline IndexSet::size_type
1841 {
1842  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1843 
1844  compress();
1845  const std::vector<Range>::const_iterator main_range =
1846  ranges.begin() + largest_range;
1847 
1848  return main_range->nth_index_in_set;
1849 }
1850 
1851 
1852 
1853 inline IndexSet::size_type
1855 {
1856  Assert(n < n_elements(), ExcIndexRangeType<size_type>(n, 0, n_elements()));
1857 
1858  compress();
1859 
1860  // first check whether the index is in the largest range
1862  std::vector<Range>::const_iterator main_range =
1863  ranges.begin() + largest_range;
1864  if (n >= main_range->nth_index_in_set &&
1865  n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1866  return main_range->begin + (n - main_range->nth_index_in_set);
1867 
1868  // find out which chunk the local index n belongs to by using a binary
1869  // search. the comparator is based on the end of the ranges. Use the
1870  // position relative to main_range to subdivide the ranges
1871  Range r(n, n + 1);
1872  r.nth_index_in_set = n;
1873  std::vector<Range>::const_iterator range_begin, range_end;
1874  if (n < main_range->nth_index_in_set)
1875  {
1876  range_begin = ranges.begin();
1877  range_end = main_range;
1878  }
1879  else
1880  {
1881  range_begin = main_range + 1;
1882  range_end = ranges.end();
1883  }
1884 
1885  const std::vector<Range>::const_iterator p =
1886  Utilities::lower_bound(range_begin, range_end, r, Range::nth_index_compare);
1887 
1888  Assert(p != ranges.end(), ExcInternalError());
1889  return p->begin + (n - p->nth_index_in_set);
1890 }
1891 
1892 
1893 
1894 inline IndexSet::size_type
1896 {
1897  // to make this call thread-safe, compress() must not be called through this
1898  // function
1899  Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1900  Assert(n < size(), ExcIndexRangeType<size_type>(n, 0, size()));
1901 
1902  // return immediately if the index set is empty
1903  if (is_empty())
1905 
1906  // check whether the index is in the largest range. use the result to
1907  // perform a one-sided binary search afterward
1909  std::vector<Range>::const_iterator main_range =
1910  ranges.begin() + largest_range;
1911  if (n >= main_range->begin && n < main_range->end)
1912  return (n - main_range->begin) + main_range->nth_index_in_set;
1913 
1914  Range r(n, n);
1915  std::vector<Range>::const_iterator range_begin, range_end;
1916  if (n < main_range->begin)
1917  {
1918  range_begin = ranges.begin();
1919  range_end = main_range;
1920  }
1921  else
1922  {
1923  range_begin = main_range + 1;
1924  range_end = ranges.end();
1925  }
1926 
1927  std::vector<Range>::const_iterator p =
1928  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1929 
1930  // if n is not in this set
1931  if (p == range_end || p->end == n || p->begin > n)
1933 
1934  Assert(p != ranges.end(), ExcInternalError());
1935  Assert(p->begin <= n, ExcInternalError());
1936  Assert(n < p->end, ExcInternalError());
1937  return (n - p->begin) + p->nth_index_in_set;
1938 }
1939 
1940 
1941 
1942 inline bool
1944 {
1945  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1946 
1947  compress();
1948  is.compress();
1949 
1950  return ranges == is.ranges;
1951 }
1952 
1953 
1954 
1955 inline bool
1957 {
1958  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1959 
1960  compress();
1961  is.compress();
1962 
1963  return ranges != is.ranges;
1964 }
1965 
1966 
1967 
1968 template <typename Vector>
1969 void
1970 IndexSet::fill_binary_vector(Vector &vector) const
1971 {
1972  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
1973 
1974  compress();
1975  // first fill all elements of the vector with zeroes.
1976  std::fill(vector.begin(), vector.end(), 0);
1977 
1978  // then write ones into the elements whose indices are contained in the
1979  // index set
1980  for (const auto &range : ranges)
1981  for (size_type i = range.begin; i < range.end; ++i)
1982  vector[i] = 1;
1983 }
1984 
1985 
1986 
1987 template <class StreamType>
1988 inline void
1989 IndexSet::print(StreamType &out) const
1990 {
1991  compress();
1992  out << "{";
1993  std::vector<Range>::const_iterator p;
1994  for (p = ranges.begin(); p != ranges.end(); ++p)
1995  {
1996  if (p->end - p->begin == 1)
1997  out << p->begin;
1998  else
1999  out << "[" << p->begin << "," << p->end - 1 << "]";
2000 
2001  if (p != --ranges.end())
2002  out << ", ";
2003  }
2004  out << "}" << std::endl;
2005 }
2006 
2007 
2008 
2009 template <class Archive>
2010 inline void
2011 IndexSet::Range::serialize(Archive &ar, const unsigned int)
2012 {
2013  ar &begin &end &nth_index_in_set;
2014 }
2015 
2016 
2017 
2018 template <class Archive>
2019 inline void
2020 IndexSet::serialize(Archive &ar, const unsigned int)
2021 {
2023 }
2024 
2025 DEAL_II_NAMESPACE_CLOSE
2026 
2027 #endif
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1079
void fill_binary_vector(VectorType &vector) const
static const unsigned int invalid_unsigned_int
Definition: types.h:187
ElementIterator end() const
Definition: index_set.h:1553
void add_index(const size_type index)
Definition: index_set.h:1626
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
const IntervalAccessor * operator->() const
Definition: index_set.h:1179
const IntervalAccessor & operator*() const
Definition: index_set.h:1171
const IndexSet * index_set
Definition: index_set.h:780
IndexSet & operator=(const IndexSet &)=default
void block_read(std::istream &in)
Definition: index_set.cc:456
size_type largest_range_starting_index() const
Definition: index_set.h:1840
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1377
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1677
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:1239
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1206
std::forward_iterator_tag iterator_category
Definition: index_set.h:673
std::size_t memory_consumption() const
Definition: index_set.cc:698
LinearAlgebra::distributed::Vector< Number > Vector
size_type n_elements() const
Definition: index_set.h:1029
size_type index_space_size
Definition: index_set.h:940
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:635
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:565
size_type size() const
Definition: index_set.h:1607
void print(StreamType &out) const
Definition: index_set.h:1989
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2011
void clear()
Definition: index_set.h:1583
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2020
static ::ExceptionBase & ExcMessage(std::string arg1)
void block_write(std::ostream &out) const
Definition: index_set.cc:442
std::vector< Range > ranges
Definition: index_set.h:924
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:238
#define Assert(cond, exc)
Definition: exceptions.h:1411
std::forward_iterator_tag iterator_category
Definition: index_set.h:764
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
bool operator<(const ElementIterator &) const
Definition: index_set.h:1364
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1895
bool operator==(const IndexSet &is) const
Definition: index_set.h:1943
signed int value_type
Definition: index_set.h:104
IntervalAccessor accessor
Definition: index_set.h:683
void do_compress() const
Definition: index_set.cc:98
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1079
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1188
size_type last() const
Definition: index_set.h:1069
const IndexSet * index_set
Definition: index_set.h:574
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1092
unsigned int n_intervals() const
Definition: index_set.h:1831
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1215
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1956
size_type pop_back()
Definition: index_set.cc:310
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:476
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1647
IndexSet operator &(const IndexSet &is) const
void set_size(const size_type size)
Definition: index_set.h:1595
ElementIterator & operator++()
Definition: index_set.h:1335
unsigned int global_dof_index
Definition: types.h:89
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:998
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:209
void compress() const
Definition: index_set.h:1615
ElementIterator end() const
Definition: index_set.h:1055
IntervalIterator begin_intervals() const
Definition: index_set.h:1562
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:985
bool is_contiguous() const
Definition: index_set.h:1789
Threads::Mutex compress_mutex
Definition: index_set.h:957
size_type pop_front()
Definition: index_set.cc:328
IntervalIterator end_intervals() const
Definition: index_set.h:1574
bool is_compressed
Definition: index_set.h:934
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1104
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1355
ElementIterator begin() const
Definition: index_set.h:1046
void write(std::ostream &out) const
Definition: index_set.cc:405
bool is_element(const size_type index) const
Definition: index_set.h:1739
bool operator==(const ElementIterator &) const
Definition: index_set.h:1295
ElementIterator at(const size_type global_index) const
Definition: index_set.h:1502
const types::global_dof_index invalid_dof_index
Definition: types.h:202
ElementIterator begin() const
Definition: index_set.h:1490
void read(std::istream &in)
Definition: index_set.cc:420
size_type n_elements() const
Definition: index_set.h:1806
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1854
bool is_empty() const
Definition: index_set.h:1798
size_type largest_range
Definition: index_set.h:951
IntervalIterator & operator++()
Definition: index_set.h:1152
static ::ExceptionBase & ExcInternalError()
size_type operator*() const
Definition: index_set.h:1282
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1197