Reference documentation for deal.II version Git 2375602 2018-06-23 12:24:17 -0500
index_set.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2018 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 #endif
35 
36 #if defined(DEAL_II_WITH_MPI) || defined(DEAL_II_WITH_PETSC)
37 # include <mpi.h>
38 #else
39 typedef int MPI_Comm;
40 # ifndef MPI_COMM_WORLD
41 # define MPI_COMM_WORLD 0
42 # endif
43 #endif
44 
45 DEAL_II_NAMESPACE_OPEN
46 
71 class IndexSet
72 {
73 public:
74  // forward declarations:
75  class ElementIterator;
77 
83 
101  typedef signed int value_type;
102 
103 
107  IndexSet();
108 
112  explicit IndexSet(const size_type size);
113 
117  IndexSet(const IndexSet &) = default;
118 
122  IndexSet &
123  operator=(const IndexSet &) = default;
124 
129  IndexSet(IndexSet &&is) noexcept;
130 
135  IndexSet &
136  operator=(IndexSet &&is) noexcept;
137 
138 #ifdef DEAL_II_WITH_TRILINOS
139 
142  explicit IndexSet(const Epetra_Map &map);
143 #endif
144 
149  void
150  clear();
151 
158  void
159  set_size(const size_type size);
160 
168  size_type
169  size() const;
170 
177  void
178  add_range(const size_type begin, const size_type end);
179 
183  void
184  add_index(const size_type index);
185 
195  template <typename ForwardIterator>
196  void
197  add_indices(const ForwardIterator &begin, const ForwardIterator &end);
198 
214  void
215  add_indices(const IndexSet &other, const unsigned int offset = 0);
216 
220  bool
221  is_element(const size_type index) const;
222 
227  bool
228  is_contiguous() const;
229 
234  bool
235  is_empty() const;
236 
246  bool
247  is_ascending_and_one_to_one(const MPI_Comm &communicator) const;
248 
252  size_type
253  n_elements() const;
254 
260  size_type
261  nth_index_in_set(const unsigned int local_index) const;
262 
269  size_type
270  index_within_set(const size_type global_index) const;
271 
280  unsigned int
281  n_intervals() const;
282 
287  unsigned int
289 
294  void
295  compress() const;
296 
302  bool
303  operator==(const IndexSet &is) const;
304 
310  bool
311  operator!=(const IndexSet &is) const;
312 
319  IndexSet operator&(const IndexSet &is) const;
320 
334  IndexSet
335  get_view(const size_type begin, const size_type end) const;
336 
342  void
343  subtract_set(const IndexSet &other);
344 
349  size_type
350  pop_back();
351 
356  size_type
357  pop_front();
358 
362  void
363  fill_index_vector(std::vector<size_type> &indices) const;
364 
377  template <typename VectorType>
378  void
379  fill_binary_vector(VectorType &vector) const;
380 
385  template <class StreamType>
386  void
387  print(StreamType &out) const;
388 
393  void
394  write(std::ostream &out) const;
395 
400  void
401  read(std::istream &in);
402 
407  void
408  block_write(std::ostream &out) const;
409 
414  void
415  block_read(std::istream &in);
416 
417 #ifdef DEAL_II_WITH_TRILINOS
418 
446  Epetra_Map
447  make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
448  const bool overlapping = false) const;
449 #endif
450 
451 
456  std::size_t
457  memory_consumption() const;
458 
460  size_type,
461  << "The global index " << arg1
462  << " is not an element of this set.");
463 
468  template <class Archive>
469  void
470  serialize(Archive &ar, const unsigned int version);
471 
472 
484  {
485  public:
490  IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
491 
495  explicit IntervalAccessor(const IndexSet *idxset);
496 
500  size_type
501  n_elements() const;
502 
506  bool
507  is_valid() const;
508 
513  begin() const;
514 
520  end() const;
521 
525  size_type
526  last() const;
527 
528  private:
532  IntervalAccessor(const IntervalAccessor &other);
537  operator=(const IntervalAccessor &other);
538 
542  bool
543  operator==(const IntervalAccessor &other) const;
547  bool
548  operator<(const IntervalAccessor &other) const;
553  void
554  advance();
559 
564  size_type range_idx;
565 
566  friend class IntervalIterator;
567  };
568 
574  {
575  public:
580  IntervalIterator(const IndexSet *idxset, const size_type range_idx);
581 
585  explicit IntervalIterator(const IndexSet *idxset);
586 
591 
595  IntervalIterator(const IntervalIterator &other);
596 
601  operator=(const IntervalIterator &other);
602 
607  operator++();
608 
613  operator++(int);
614 
618  const IntervalAccessor &operator*() const;
619 
623  const IntervalAccessor *operator->() const;
624 
628  bool
629  operator==(const IntervalIterator &) const;
630 
634  bool
635  operator!=(const IntervalIterator &) const;
636 
640  bool
641  operator<(const IntervalIterator &) const;
642 
649  int
650  operator-(const IntervalIterator &p) const;
651 
657  typedef std::forward_iterator_tag iterator_category;
658  typedef IntervalAccessor value_type;
659  typedef std::ptrdiff_t difference_type;
660  typedef IntervalAccessor * pointer;
661  typedef IntervalAccessor & reference;
662 
663  private:
668  };
669 
675  {
676  public:
681  ElementIterator(const IndexSet *idxset,
682  const size_type range_idx,
683  const size_type index);
684 
688  explicit ElementIterator(const IndexSet *idxset);
689 
694  size_type operator*() const;
695 
699  bool
700  is_valid() const;
701 
706  operator++();
707 
712  operator++(int);
713 
717  bool
718  operator==(const ElementIterator &) const;
719 
723  bool
724  operator!=(const ElementIterator &) const;
725 
729  bool
730  operator<(const ElementIterator &) const;
731 
740  std::ptrdiff_t
741  operator-(const ElementIterator &p) const;
742 
748  typedef std::forward_iterator_tag iterator_category;
749  typedef size_type value_type;
750  typedef std::ptrdiff_t difference_type;
751  typedef size_type * pointer;
752  typedef size_type & reference;
753 
754  private:
758  void
759  advance();
760 
768  size_type range_idx;
772  size_type idx;
773  };
774 
780  begin() const;
781 
797  at(const size_type global_index) const;
798 
804  end() const;
805 
810  begin_intervals() const;
811 
817  end_intervals() const;
818 
823 private:
832  struct Range
833  {
834  size_type begin;
835  size_type end;
836 
837  size_type nth_index_in_set;
838 
847  Range();
848 
856  Range(const size_type i1, const size_type i2);
857 
858  friend inline bool
859  operator<(const Range &range_1, const Range &range_2)
860  {
861  return (
862  (range_1.begin < range_2.begin) ||
863  ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
864  }
865 
866  static bool
867  end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
868  {
869  return x.end < y.end;
870  }
871 
872  static bool
873  nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
874  {
875  return (x.nth_index_in_set + (x.end - x.begin) <
876  y.nth_index_in_set + (y.end - y.begin));
877  }
878 
879  friend inline bool
880  operator==(const Range &range_1, const Range &range_2)
881  {
882  return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
883  }
884 
885  static std::size_t
886  memory_consumption()
887  {
888  return sizeof(Range);
889  }
890 
895  template <class Archive>
896  void
897  serialize(Archive &ar, const unsigned int version);
898  };
899 
908  mutable std::vector<Range> ranges;
909 
918  mutable bool is_compressed;
919 
924  size_type index_space_size;
925 
935  mutable size_type largest_range;
936 
942 
946  void
947  do_compress() const;
948 };
949 
950 
968 inline IndexSet
969 complete_index_set(const unsigned int N)
970 {
971  IndexSet is(N);
972  is.add_range(0, N);
973  is.compress();
974  return is;
975 }
976 
977 /* ------------------ inline functions ------------------ */
978 
979 
980 /* IntervalAccessor */
981 
983  const IndexSet * idxset,
984  const IndexSet::size_type range_idx)
985  : index_set(idxset)
986  , range_idx(range_idx)
987 {
988  Assert(range_idx < idxset->n_intervals(),
989  ExcInternalError("Invalid range index"));
990 }
991 
992 
993 
995  : index_set(idxset)
996  , range_idx(numbers::invalid_dof_index)
997 {}
998 
999 
1000 
1002  const IndexSet::IntervalAccessor &other)
1003  : index_set(other.index_set)
1004  , range_idx(other.range_idx)
1005 {
1007  ExcMessage("invalid iterator"));
1008 }
1009 
1010 
1011 
1012 inline IndexSet::size_type
1014 {
1015  Assert(is_valid(), ExcMessage("invalid iterator"));
1016  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1017 }
1018 
1019 
1020 
1021 inline bool
1023 {
1024  return index_set != nullptr && range_idx < index_set->n_intervals();
1025 }
1026 
1027 
1028 
1031 {
1032  Assert(is_valid(), ExcMessage("invalid iterator"));
1033  return IndexSet::ElementIterator(index_set,
1034  range_idx,
1035  index_set->ranges[range_idx].begin);
1036 }
1037 
1038 
1039 
1042 {
1043  Assert(is_valid(), ExcMessage("invalid iterator"));
1044 
1045  // point to first index in next interval unless we are the last interval.
1046  if (range_idx < index_set->ranges.size() - 1)
1047  return IndexSet::ElementIterator(index_set,
1048  range_idx + 1,
1049  index_set->ranges[range_idx + 1].begin);
1050  else
1051  return index_set->end();
1052 }
1053 
1054 
1055 
1056 inline IndexSet::size_type
1058 {
1059  Assert(is_valid(), ExcMessage("invalid iterator"));
1060 
1061  return index_set->ranges[range_idx].end - 1;
1062 }
1063 
1064 
1065 
1068 {
1069  index_set = other.index_set;
1070  range_idx = other.range_idx;
1071  Assert(range_idx == numbers::invalid_dof_index || is_valid(),
1072  ExcMessage("invalid iterator"));
1073  return *this;
1074 }
1075 
1076 
1077 
1078 inline bool
1081 {
1082  Assert(index_set == other.index_set,
1083  ExcMessage(
1084  "Can not compare accessors pointing to different IndexSets"));
1085  return range_idx == other.range_idx;
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 void
1104 {
1105  Assert(
1106  is_valid(),
1107  ExcMessage(
1108  "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1109  ++range_idx;
1110 
1111  // set ourselves to invalid if we walk off the end
1112  if (range_idx >= index_set->ranges.size())
1113  range_idx = numbers::invalid_dof_index;
1114 }
1115 
1116 
1117 /* IntervalIterator */
1118 
1120  const IndexSet * idxset,
1121  const IndexSet::size_type range_idx)
1122  : accessor(idxset, range_idx)
1123 {}
1124 
1125 
1126 
1128  : accessor(nullptr)
1129 {}
1130 
1131 
1132 
1134  : accessor(idxset)
1135 {}
1136 
1137 
1138 
1140  const IndexSet::IntervalIterator &other)
1141  : accessor(other.accessor)
1142 {}
1143 
1144 
1145 
1148 {
1149  accessor = other.accessor;
1150  return *this;
1151 }
1152 
1153 
1154 
1157 {
1158  accessor.advance();
1159  return *this;
1160 }
1161 
1162 
1163 
1166 {
1167  const IndexSet::IntervalIterator iter = *this;
1168  accessor.advance();
1169  return iter;
1170 }
1171 
1172 
1173 
1175  operator*() const
1176 {
1177  return accessor;
1178 }
1179 
1180 
1181 
1184 {
1185  return &accessor;
1186 }
1187 
1188 
1189 
1190 inline bool
1193 {
1194  return accessor == other.accessor;
1195 }
1196 
1197 
1198 
1199 inline bool
1202 {
1203  return !(*this == other);
1204 }
1205 
1206 
1207 
1208 inline bool
1211 {
1212  return accessor < other.accessor;
1213 }
1214 
1215 
1216 
1217 inline int
1220 {
1221  Assert(accessor.index_set == other.accessor.index_set,
1222  ExcMessage(
1223  "Can not compare iterators belonging to different IndexSets"));
1224 
1225  const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index) ?
1226  accessor.index_set->ranges.size() :
1227  accessor.range_idx;
1228  const size_type rhs =
1230  accessor.index_set->ranges.size() :
1231  other.accessor.range_idx;
1232 
1233  if (lhs > rhs)
1234  return static_cast<int>(lhs - rhs);
1235  else
1236  return -static_cast<int>(rhs - lhs);
1237 }
1238 
1239 
1240 
1241 /* ElementIterator */
1242 
1244  const IndexSet * idxset,
1245  const IndexSet::size_type range_idx,
1246  const IndexSet::size_type index)
1247  : index_set(idxset)
1248  , range_idx(range_idx)
1249  , idx(index)
1250 {
1251  Assert(range_idx < index_set->ranges.size(),
1252  ExcMessage(
1253  "Invalid range index for IndexSet::ElementIterator constructor."));
1254  Assert(
1255  idx >= index_set->ranges[range_idx].begin &&
1256  idx < index_set->ranges[range_idx].end,
1258  "Invalid index argument for IndexSet::ElementIterator constructor."));
1259 }
1260 
1261 
1262 
1264  : index_set(idxset)
1265  , range_idx(numbers::invalid_dof_index)
1266  , idx(numbers::invalid_dof_index)
1267 {}
1268 
1269 
1270 
1271 inline bool
1273 {
1274  Assert((range_idx == numbers::invalid_dof_index &&
1275  idx == numbers::invalid_dof_index) ||
1276  (range_idx < index_set->ranges.size() &&
1277  idx < index_set->ranges[range_idx].end),
1278  ExcInternalError("Invalid ElementIterator state."));
1279 
1280  return (range_idx < index_set->ranges.size() &&
1281  idx < index_set->ranges[range_idx].end);
1282 }
1283 
1284 
1285 
1287 {
1288  Assert(
1289  is_valid(),
1290  ExcMessage(
1291  "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1292  return idx;
1293 }
1294 
1295 
1296 
1297 inline bool
1300 {
1301  Assert(index_set == other.index_set,
1302  ExcMessage(
1303  "Can not compare iterators belonging to different IndexSets"));
1304  return range_idx == other.range_idx && idx == other.idx;
1305 }
1306 
1307 
1308 
1309 inline void
1311 {
1312  Assert(
1313  is_valid(),
1314  ExcMessage(
1315  "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1316  if (idx < index_set->ranges[range_idx].end)
1317  ++idx;
1318  // end of this range?
1319  if (idx == index_set->ranges[range_idx].end)
1320  {
1321  // point to first element in next interval if possible
1322  if (range_idx < index_set->ranges.size() - 1)
1323  {
1324  ++range_idx;
1325  idx = index_set->ranges[range_idx].begin;
1326  }
1327  else
1328  {
1329  // we just fell off the end, set to invalid:
1330  range_idx = numbers::invalid_dof_index;
1332  }
1333  }
1334 }
1335 
1336 
1337 
1340 {
1341  advance();
1342  return *this;
1343 }
1344 
1345 
1346 
1349 {
1350  const IndexSet::ElementIterator it = *this;
1351  advance();
1352  return it;
1353 }
1354 
1355 
1356 
1357 inline bool
1360 {
1361  return !(*this == other);
1362 }
1363 
1364 
1365 
1366 inline bool
1369 {
1370  Assert(index_set == other.index_set,
1371  ExcMessage(
1372  "Can not compare iterators belonging to different IndexSets"));
1373  return range_idx < other.range_idx ||
1374  (range_idx == other.range_idx && idx < other.idx);
1375 }
1376 
1377 
1378 
1379 inline std::ptrdiff_t
1382 {
1383  Assert(index_set == other.index_set,
1384  ExcMessage(
1385  "Can not compare iterators belonging to different IndexSets"));
1386  if (*this == other)
1387  return 0;
1388  if (!(*this < other))
1389  return -(other - *this);
1390 
1391  // only other can be equal to end() because of the checks above.
1392  Assert(is_valid(), ExcInternalError());
1393 
1394  // Note: we now compute how far advance *this in "*this < other" to get other,
1395  // so we need to return -c at the end.
1396 
1397  // first finish the current range:
1398  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1399 
1400  // now walk in steps of ranges (need to start one behind our current one):
1401  for (size_type range = range_idx + 1;
1402  range < index_set->ranges.size() && range <= other.range_idx;
1403  ++range)
1404  c += index_set->ranges[range].end - index_set->ranges[range].begin;
1405 
1406  Assert(
1407  other.range_idx < index_set->ranges.size() ||
1409  ExcMessage(
1410  "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1411 
1412  // We might have walked too far because we went until the end of
1413  // other.range_idx, so walk backwards to other.idx:
1415  c -= index_set->ranges[other.range_idx].end - other.idx;
1416 
1417  return -c;
1418 }
1419 
1420 
1421 /* Range */
1422 
1424  : begin(numbers::invalid_dof_index)
1425  , end(numbers::invalid_dof_index)
1426  , nth_index_in_set(numbers::invalid_dof_index)
1427 {}
1428 
1429 
1430 
1431 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1432  : begin(i1)
1433  , end(i2)
1434  , nth_index_in_set(numbers::invalid_dof_index)
1435 {}
1436 
1437 
1438 
1439 /* IndexSet itself */
1440 
1442  : is_compressed(true)
1443  , index_space_size(0)
1444  , largest_range(numbers::invalid_unsigned_int)
1445 {}
1446 
1447 
1448 
1449 inline IndexSet::IndexSet(const size_type size)
1450  : is_compressed(true)
1451  , index_space_size(size)
1452  , largest_range(numbers::invalid_unsigned_int)
1453 {}
1454 
1455 
1456 
1457 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1458  : ranges(std::move(is.ranges))
1459  , is_compressed(is.is_compressed)
1460  , index_space_size(is.index_space_size)
1461  , largest_range(is.largest_range)
1462 {
1463  is.ranges.clear();
1464  is.is_compressed = true;
1465  is.index_space_size = 0;
1466  is.largest_range = numbers::invalid_unsigned_int;
1467 
1468  compress();
1469 }
1470 
1471 
1472 
1473 inline IndexSet &
1475 {
1476  ranges = std::move(is.ranges);
1477  is_compressed = is.is_compressed;
1478  index_space_size = is.index_space_size;
1479  largest_range = is.largest_range;
1480 
1481  is.ranges.clear();
1482  is.is_compressed = true;
1483  is.index_space_size = 0;
1484  is.largest_range = numbers::invalid_unsigned_int;
1485 
1486  compress();
1487 
1488  return *this;
1489 }
1490 
1491 
1492 
1495 {
1496  compress();
1497  if (ranges.size() > 0)
1498  return IndexSet::ElementIterator(this, 0, ranges[0].begin);
1499  else
1500  return end();
1501 }
1502 
1503 
1504 
1506 IndexSet::at(const size_type global_index) const
1507 {
1508  compress();
1509  Assert(global_index < size(),
1510  ExcIndexRangeType<size_type>(global_index, 0, size()));
1511 
1512  if (ranges.empty())
1513  return end();
1514 
1515  std::vector<Range>::const_iterator main_range =
1516  ranges.begin() + largest_range;
1517 
1518  Range r(global_index, global_index + 1);
1519  // This optimization makes the bounds for lower_bound smaller by checking
1520  // the largest range first.
1521  std::vector<Range>::const_iterator range_begin, range_end;
1522  if (global_index < main_range->begin)
1523  {
1524  range_begin = ranges.begin();
1525  range_end = main_range;
1526  }
1527  else
1528  {
1529  range_begin = main_range;
1530  range_end = ranges.end();
1531  }
1532 
1533  // This will give us the first range p=[a,b[ with b>=global_index using
1534  // a binary search
1535  const std::vector<Range>::const_iterator p =
1536  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1537 
1538  // We couldn't find a range, which means we have no range that contains
1539  // global_index and also no range behind it, meaning we need to return end().
1540  if (p == ranges.end())
1541  return end();
1542 
1543  // Finally, we can have two cases: Either global_index is not in [a,b[,
1544  // which means we need to return an iterator to a because global_index, ...,
1545  // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
1546  // [a,b[ and we will return an iterator pointing directly at global_index
1547  // (else branch).
1548  if (global_index < p->begin)
1549  return IndexSet::ElementIterator(this, p - ranges.begin(), p->begin);
1550  else
1551  return IndexSet::ElementIterator(this, p - ranges.begin(), global_index);
1552 }
1553 
1554 
1555 
1558 {
1559  compress();
1560  return IndexSet::ElementIterator(this);
1561 }
1562 
1563 
1564 
1567 {
1568  compress();
1569  if (ranges.size() > 0)
1570  return IndexSet::IntervalIterator(this, 0);
1571  else
1572  return end_intervals();
1573 }
1574 
1575 
1576 
1579 {
1580  compress();
1581  return IndexSet::IntervalIterator(this);
1582 }
1583 
1584 
1585 
1586 inline void
1588 {
1589  // reset so that there are no indices in the set any more; however,
1590  // as documented, the index set retains its size
1591  ranges.clear();
1592  is_compressed = true;
1594 }
1595 
1596 
1597 
1598 inline void
1599 IndexSet::set_size(const size_type sz)
1600 {
1601  Assert(ranges.empty(),
1602  ExcMessage("This function can only be called if the current "
1603  "object does not yet contain any elements."));
1604  index_space_size = sz;
1605  is_compressed = true;
1606 }
1607 
1608 
1609 
1610 inline IndexSet::size_type
1612 {
1613  return index_space_size;
1614 }
1615 
1616 
1617 
1618 inline void
1620 {
1621  if (is_compressed == true)
1622  return;
1623 
1624  do_compress();
1625 }
1626 
1627 
1628 
1629 inline void
1630 IndexSet::add_index(const size_type index)
1631 {
1632  Assert(index < index_space_size,
1633  ExcIndexRangeType<size_type>(index, 0, index_space_size));
1634 
1635  const Range new_range(index, index + 1);
1636  if (ranges.size() == 0 || index > ranges.back().end)
1637  ranges.push_back(new_range);
1638  else if (index == ranges.back().end)
1639  ranges.back().end++;
1640  else
1641  ranges.insert(Utilities::lower_bound(ranges.begin(),
1642  ranges.end(),
1643  new_range),
1644  new_range);
1645  is_compressed = false;
1646 }
1647 
1648 
1649 
1650 template <typename ForwardIterator>
1651 inline void
1652 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1653 {
1654  // insert each element of the range. if some of them happen to be
1655  // consecutive, merge them to a range
1656  for (ForwardIterator p = begin; p != end;)
1657  {
1658  const size_type begin_index = *p;
1659  size_type end_index = begin_index + 1;
1660  ForwardIterator q = p;
1661  ++q;
1662  while ((q != end) && (*q == end_index))
1663  {
1664  ++end_index;
1665  ++q;
1666  }
1667 
1668  add_range(begin_index, end_index);
1669  p = q;
1670  }
1671 }
1672 
1673 
1674 
1675 inline bool
1676 IndexSet::is_element(const size_type index) const
1677 {
1678  if (ranges.empty() == false)
1679  {
1680  compress();
1681 
1682  // fast check whether the index is in the largest range
1684  if (index >= ranges[largest_range].begin &&
1685  index < ranges[largest_range].end)
1686  return true;
1687 
1688  // get the element after which we would have to insert a range that
1689  // consists of all elements from this element to the end of the index
1690  // range plus one. after this call we know that if p!=end() then
1691  // p->begin<=index unless there is no such range at all
1692  //
1693  // if the searched for element is an element of this range, then we're
1694  // done. otherwise, the element can't be in one of the following ranges
1695  // because otherwise p would be a different iterator
1696  //
1697  // since we already know the position relative to the largest range (we
1698  // called compress!), we can perform the binary search on ranges with
1699  // lower/higher number compared to the largest range
1700  std::vector<Range>::const_iterator p = std::upper_bound(
1701  ranges.begin() +
1702  (index < ranges[largest_range].begin ? 0 : largest_range + 1),
1703  index < ranges[largest_range].begin ? ranges.begin() + largest_range :
1704  ranges.end(),
1705  Range(index, size() + 1));
1706 
1707  if (p == ranges.begin())
1708  return ((index >= p->begin) && (index < p->end));
1709 
1710  Assert((p == ranges.end()) || (p->begin > index), ExcInternalError());
1711 
1712  // now move to that previous range
1713  --p;
1714  Assert(p->begin <= index, ExcInternalError());
1715 
1716  return (p->end > index);
1717  }
1718 
1719  // didn't find this index, so it's not in the set
1720  return false;
1721 }
1722 
1723 
1724 
1725 inline bool
1727 {
1728  compress();
1729  return (ranges.size() <= 1);
1730 }
1731 
1732 
1733 
1734 inline bool
1736 {
1737  return ranges.empty();
1738 }
1739 
1740 
1741 
1742 inline IndexSet::size_type
1744 {
1745  // make sure we have non-overlapping ranges
1746  compress();
1747 
1748  size_type v = 0;
1749  if (!ranges.empty())
1750  {
1751  Range &r = ranges.back();
1752  v = r.nth_index_in_set + r.end - r.begin;
1753  }
1754 
1755 #ifdef DEBUG
1756  size_type s = 0;
1757  for (std::vector<Range>::iterator range = ranges.begin();
1758  range != ranges.end();
1759  ++range)
1760  s += (range->end - range->begin);
1761  Assert(s == v, ExcInternalError());
1762 #endif
1763 
1764  return v;
1765 }
1766 
1767 
1768 
1769 inline unsigned int
1771 {
1772  compress();
1773  return ranges.size();
1774 }
1775 
1776 
1777 
1778 inline unsigned int
1780 {
1781  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1782 
1783  compress();
1784  const std::vector<Range>::const_iterator main_range =
1785  ranges.begin() + largest_range;
1786 
1787  return main_range->nth_index_in_set;
1788 }
1789 
1790 
1791 
1792 inline IndexSet::size_type
1793 IndexSet::nth_index_in_set(const unsigned int n) const
1794 {
1795  Assert(n < n_elements(), ExcIndexRangeType<size_type>(n, 0, n_elements()));
1796 
1797  compress();
1798 
1799  // first check whether the index is in the largest range
1801  std::vector<Range>::const_iterator main_range =
1802  ranges.begin() + largest_range;
1803  if (n >= main_range->nth_index_in_set &&
1804  n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1805  return main_range->begin + (n - main_range->nth_index_in_set);
1806 
1807  // find out which chunk the local index n belongs to by using a binary
1808  // search. the comparator is based on the end of the ranges. Use the
1809  // position relative to main_range to subdivide the ranges
1810  Range r(n, n + 1);
1811  r.nth_index_in_set = n;
1812  std::vector<Range>::const_iterator range_begin, range_end;
1813  if (n < main_range->nth_index_in_set)
1814  {
1815  range_begin = ranges.begin();
1816  range_end = main_range;
1817  }
1818  else
1819  {
1820  range_begin = main_range + 1;
1821  range_end = ranges.end();
1822  }
1823 
1824  const std::vector<Range>::const_iterator p =
1825  Utilities::lower_bound(range_begin, range_end, r, Range::nth_index_compare);
1826 
1827  Assert(p != ranges.end(), ExcInternalError());
1828  return p->begin + (n - p->nth_index_in_set);
1829 }
1830 
1831 
1832 
1833 inline IndexSet::size_type
1834 IndexSet::index_within_set(const size_type n) const
1835 {
1836  // to make this call thread-safe, compress() must not be called through this
1837  // function
1838  Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1839  Assert(n < size(), ExcIndexRangeType<size_type>(n, 0, size()));
1840 
1841  // return immediately if the index set is empty
1842  if (is_empty())
1844 
1845  // check whether the index is in the largest range. use the result to
1846  // perform a one-sided binary search afterward
1848  std::vector<Range>::const_iterator main_range =
1849  ranges.begin() + largest_range;
1850  if (n >= main_range->begin && n < main_range->end)
1851  return (n - main_range->begin) + main_range->nth_index_in_set;
1852 
1853  Range r(n, n);
1854  std::vector<Range>::const_iterator range_begin, range_end;
1855  if (n < main_range->begin)
1856  {
1857  range_begin = ranges.begin();
1858  range_end = main_range;
1859  }
1860  else
1861  {
1862  range_begin = main_range + 1;
1863  range_end = ranges.end();
1864  }
1865 
1866  std::vector<Range>::const_iterator p =
1867  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1868 
1869  // if n is not in this set
1870  if (p == range_end || p->end == n || p->begin > n)
1872 
1873  Assert(p != ranges.end(), ExcInternalError());
1874  Assert(p->begin <= n, ExcInternalError());
1875  Assert(n < p->end, ExcInternalError());
1876  return (n - p->begin) + p->nth_index_in_set;
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 inline bool
1896 {
1897  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1898 
1899  compress();
1900  is.compress();
1901 
1902  return ranges != is.ranges;
1903 }
1904 
1905 
1906 
1907 template <typename Vector>
1908 void
1909 IndexSet::fill_binary_vector(Vector &vector) const
1910 {
1911  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
1912 
1913  compress();
1914  // first fill all elements of the vector with zeroes.
1915  std::fill(vector.begin(), vector.end(), 0);
1916 
1917  // then write ones into the elements whose indices are contained in the
1918  // index set
1919  for (std::vector<Range>::iterator it = ranges.begin(); it != ranges.end();
1920  ++it)
1921  for (size_type i = it->begin; i < it->end; ++i)
1922  vector[i] = 1;
1923 }
1924 
1925 
1926 
1927 template <class StreamType>
1928 inline void
1929 IndexSet::print(StreamType &out) const
1930 {
1931  compress();
1932  out << "{";
1933  std::vector<Range>::const_iterator p;
1934  for (p = ranges.begin(); p != ranges.end(); ++p)
1935  {
1936  if (p->end - p->begin == 1)
1937  out << p->begin;
1938  else
1939  out << "[" << p->begin << "," << p->end - 1 << "]";
1940 
1941  if (p != --ranges.end())
1942  out << ", ";
1943  }
1944  out << "}" << std::endl;
1945 }
1946 
1947 
1948 
1949 template <class Archive>
1950 inline void
1951 IndexSet::Range::serialize(Archive &ar, const unsigned int)
1952 {
1953  ar &begin &end &nth_index_in_set;
1954 }
1955 
1956 
1957 
1958 template <class Archive>
1959 inline void
1960 IndexSet::serialize(Archive &ar, const unsigned int)
1961 {
1963 }
1964 
1965 DEAL_II_NAMESPACE_CLOSE
1966 
1967 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:912
bool operator<(const ElementIterator &) const
Definition: index_set.h:1368
static const unsigned int invalid_unsigned_int
Definition: types.h:173
void add_index(const size_type index)
Definition: index_set.h:1630
IntervalIterator begin_intervals() const
Definition: index_set.h:1566
size_type n_elements() const
Definition: index_set.h:1013
const IndexSet * index_set
Definition: index_set.h:764
IndexSet & operator=(const IndexSet &)=default
bool is_valid() const
Definition: index_set.h:1272
void block_read(std::istream &in)
Definition: index_set.cc:483
signed int value_type
Definition: index_set.h:101
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:593
size_type n_elements() const
Definition: index_set.h:1743
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1652
void fill_binary_vector(VectorType &vector) const
bool operator==(const IndexSet &is) const
Definition: index_set.h:1882
STL namespace.
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1192
bool is_empty() const
Definition: index_set.h:1735
size_type size() const
Definition: index_set.h:1611
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition: index_set.h:1243
LinearAlgebra::distributed::Vector< Number > Vector
void print(StreamType &out) const
Definition: index_set.h:1929
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1080
size_type index_space_size
Definition: index_set.h:924
ElementIterator begin() const
Definition: index_set.h:1030
void block_write(std::ostream &out) const
Definition: index_set.cc:469
bool is_contiguous() const
Definition: index_set.h:1726
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1951
unsigned int n_intervals() const
Definition: index_set.h:1770
IndexSet complete_index_set(const unsigned int N)
Definition: index_set.h:969
void clear()
Definition: index_set.h:1587
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1960
static::ExceptionBase & ExcMessage(std::string arg1)
void write(std::ostream &out) const
Definition: index_set.cc:432
size_type operator*() const
Definition: index_set.h:1286
IntervalIterator & operator=(const IntervalIterator &other)
Definition: index_set.h:1147
std::vector< Range > ranges
Definition: index_set.h:908
static::ExceptionBase & ExcIndexNotPresent(size_type arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:408
unsigned int global_dof_index
Definition: types.h:89
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1895
void subtract_set(const IndexSet &other)
Definition: index_set.cc:265
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:503
#define Assert(cond, exc)
Definition: exceptions.h:1227
unsigned int largest_range_starting_index() const
Definition: index_set.h:1779
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const IntervalAccessor & operator*() const
Definition: index_set.h:1175
IntervalAccessor accessor
Definition: index_set.h:667
IndexSet operator&(const IndexSet &is) const
Definition: index_set.cc:189
ElementIterator at(const size_type global_index) const
Definition: index_set.h:1506
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1067
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:236
void compress() const
Definition: index_set.h:1619
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1201
const IndexSet * index_set
Definition: index_set.h:558
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:523
size_type pop_back()
Definition: index_set.cc:337
ElementIterator begin() const
Definition: index_set.h:1494
void add_range(const size_type begin, const size_type end)
Definition: index_set.cc:96
void set_size(const size_type size)
Definition: index_set.h:1599
std::forward_iterator_tag iterator_category
Definition: index_set.h:748
ElementIterator & operator++()
Definition: index_set.h:1339
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:982
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1210
ElementIterator end() const
Definition: index_set.h:1041
Threads::Mutex compress_mutex
Definition: index_set.h:941
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1092
void do_compress() const
Definition: index_set.cc:125
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1359
size_type pop_front()
Definition: index_set.cc:355
bool is_compressed
Definition: index_set.h:918
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1834
bool is_element(const size_type index) const
Definition: index_set.h:1676
std::size_t memory_consumption() const
Definition: index_set.cc:656
ElementIterator end() const
Definition: index_set.h:1557
const types::global_dof_index invalid_dof_index
Definition: types.h:188
size_type nth_index_in_set(const unsigned int local_index) const
Definition: index_set.h:1793
bool operator==(const ElementIterator &) const
Definition: index_set.h:1299
size_type last() const
Definition: index_set.h:1057
types::global_dof_index size_type
Definition: index_set.h:76
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1219
void read(std::istream &in)
Definition: index_set.cc:447
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1381
size_type largest_range
Definition: index_set.h:935
const IntervalAccessor * operator->() const
Definition: index_set.h:1183
IntervalIterator & operator++()
Definition: index_set.h:1156
std::forward_iterator_tag iterator_category
Definition: index_set.h:657
IntervalIterator end_intervals() const
Definition: index_set.h:1578
static::ExceptionBase & ExcInternalError()