Reference documentation for deal.II version Git 9155a679b0 2020-01-17 08:18:08 -0500
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
371  IndexSet
372  tensor_product(const IndexSet &other) const;
373 
378  size_type
379  pop_back();
380 
385  size_type
386  pop_front();
387 
391  void
392  fill_index_vector(std::vector<size_type> &indices) const;
393 
406  template <typename VectorType>
407  void
408  fill_binary_vector(VectorType &vector) const;
409 
414  template <class StreamType>
415  void
416  print(StreamType &out) const;
417 
422  void
423  write(std::ostream &out) const;
424 
429  void
430  read(std::istream &in);
431 
436  void
437  block_write(std::ostream &out) const;
438 
443  void
444  block_read(std::istream &in);
445 
446 #ifdef DEAL_II_WITH_TRILINOS
447 
475  Epetra_Map
476  make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
477  const bool overlapping = false) const;
478 
479 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
480  Tpetra::Map<int, types::global_dof_index>
481  make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
482  const bool overlapping = false) const;
483 # endif
484 #endif
485 
486 
491  std::size_t
492  memory_consumption() const;
493 
495  size_type,
496  << "The global index " << arg1
497  << " is not an element of this set.");
498 
503  template <class Archive>
504  void
505  serialize(Archive &ar, const unsigned int version);
506 
507 
519  {
520  public:
525  IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
526 
530  explicit IntervalAccessor(const IndexSet *idxset);
531 
535  size_type
536  n_elements() const;
537 
541  bool
542  is_valid() const;
543 
548  begin() const;
549 
555  end() const;
556 
560  size_type
561  last() const;
562 
563  private:
567  IntervalAccessor(const IntervalAccessor &other);
572  operator=(const IntervalAccessor &other);
573 
577  bool
578  operator==(const IntervalAccessor &other) const;
582  bool
583  operator<(const IntervalAccessor &other) const;
588  void
589  advance();
594 
600 
601  friend class IntervalIterator;
602  };
603 
609  {
610  public:
615  IntervalIterator(const IndexSet *idxset, const size_type range_idx);
616 
620  explicit IntervalIterator(const IndexSet *idxset);
621 
626 
630  IntervalIterator(const IntervalIterator &other) = default;
631 
636  operator=(const IntervalIterator &other) = default;
637 
642  operator++();
643 
648  operator++(int);
649 
653  const IntervalAccessor &operator*() const;
654 
658  const IntervalAccessor *operator->() const;
659 
663  bool
664  operator==(const IntervalIterator &) const;
665 
669  bool
670  operator!=(const IntervalIterator &) const;
671 
675  bool
676  operator<(const IntervalIterator &) const;
677 
684  int
685  operator-(const IntervalIterator &p) const;
686 
692  using iterator_category = std::forward_iterator_tag;
694  using difference_type = std::ptrdiff_t;
695  using pointer = IntervalAccessor *;
696  using reference = IntervalAccessor &;
697 
698  private:
703  };
704 
710  {
711  public:
716  ElementIterator(const IndexSet *idxset,
717  const size_type range_idx,
718  const size_type index);
719 
723  explicit ElementIterator(const IndexSet *idxset);
724 
729  size_type operator*() const;
730 
734  bool
735  is_valid() const;
736 
741  operator++();
742 
747  operator++(int);
748 
752  bool
753  operator==(const ElementIterator &) const;
754 
758  bool
759  operator!=(const ElementIterator &) const;
760 
764  bool
765  operator<(const ElementIterator &) const;
766 
775  std::ptrdiff_t
776  operator-(const ElementIterator &p) const;
777 
783  using iterator_category = std::forward_iterator_tag;
784  using value_type = size_type;
785  using difference_type = std::ptrdiff_t;
786  using pointer = size_type *;
787  using reference = size_type &;
788 
789  private:
793  void
794  advance();
795 
808  };
809 
815  begin() const;
816 
832  at(const size_type global_index) const;
833 
839  end() const;
840 
845  begin_intervals() const;
846 
852  end_intervals() const;
853 
858 private:
867  struct Range
868  {
870  size_type end;
871 
873 
882  Range();
883 
891  Range(const size_type i1, const size_type i2);
892 
893  friend inline bool
894  operator<(const Range &range_1, const Range &range_2)
895  {
896  return (
897  (range_1.begin < range_2.begin) ||
898  ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
899  }
900 
901  static bool
902  end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
903  {
904  return x.end < y.end;
905  }
906 
907  static bool
908  nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
909  {
910  return (x.nth_index_in_set + (x.end - x.begin) <
911  y.nth_index_in_set + (y.end - y.begin));
912  }
913 
914  friend inline bool
915  operator==(const Range &range_1, const Range &range_2)
916  {
917  return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
918  }
919 
920  static std::size_t
922  {
923  return sizeof(Range);
924  }
925 
930  template <class Archive>
931  void
932  serialize(Archive &ar, const unsigned int version);
933  };
934 
943  mutable std::vector<Range> ranges;
944 
953  mutable bool is_compressed;
954 
960 
971 
977 
981  void
982  do_compress() const;
983 };
984 
985 
1003 inline IndexSet
1005 {
1006  IndexSet is(N);
1007  is.add_range(0, N);
1008  is.compress();
1009  return is;
1010 }
1011 
1012 /* ------------------ inline functions ------------------ */
1013 
1014 
1015 /* IntervalAccessor */
1016 
1018  const IndexSet * idxset,
1020  : index_set(idxset)
1021  , range_idx(range_idx)
1022 {
1023  Assert(range_idx < idxset->n_intervals(),
1024  ExcInternalError("Invalid range index"));
1025 }
1026 
1027 
1028 
1030  : index_set(idxset)
1031  , range_idx(numbers::invalid_dof_index)
1032 {}
1033 
1034 
1035 
1037  const IndexSet::IntervalAccessor &other)
1038  : index_set(other.index_set)
1039  , range_idx(other.range_idx)
1040 {
1042  ExcMessage("invalid iterator"));
1043 }
1044 
1045 
1046 
1047 inline IndexSet::size_type
1049 {
1050  Assert(is_valid(), ExcMessage("invalid iterator"));
1051  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1052 }
1053 
1054 
1055 
1056 inline bool
1058 {
1059  return index_set != nullptr && range_idx < index_set->n_intervals();
1060 }
1061 
1062 
1063 
1066 {
1067  Assert(is_valid(), ExcMessage("invalid iterator"));
1068  return {index_set, range_idx, index_set->ranges[range_idx].begin};
1069 }
1070 
1071 
1072 
1075 {
1076  Assert(is_valid(), ExcMessage("invalid iterator"));
1077 
1078  // point to first index in next interval unless we are the last interval.
1079  if (range_idx < index_set->ranges.size() - 1)
1080  return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1081  else
1082  return index_set->end();
1083 }
1084 
1085 
1086 
1087 inline IndexSet::size_type
1089 {
1090  Assert(is_valid(), ExcMessage("invalid iterator"));
1091 
1092  return index_set->ranges[range_idx].end - 1;
1093 }
1094 
1095 
1096 
1099 {
1100  index_set = other.index_set;
1101  range_idx = other.range_idx;
1103  ExcMessage("invalid iterator"));
1104  return *this;
1105 }
1106 
1107 
1108 
1109 inline bool
1112 {
1113  Assert(index_set == other.index_set,
1114  ExcMessage(
1115  "Can not compare accessors pointing to different IndexSets"));
1116  return range_idx == other.range_idx;
1117 }
1118 
1119 
1120 
1121 inline bool
1124 {
1125  Assert(index_set == other.index_set,
1126  ExcMessage(
1127  "Can not compare accessors pointing to different IndexSets"));
1128  return range_idx < other.range_idx;
1129 }
1130 
1131 
1132 
1133 inline void
1135 {
1136  Assert(
1137  is_valid(),
1138  ExcMessage(
1139  "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1140  ++range_idx;
1141 
1142  // set ourselves to invalid if we walk off the end
1143  if (range_idx >= index_set->ranges.size())
1145 }
1146 
1147 
1148 /* IntervalIterator */
1149 
1151  const IndexSet * idxset,
1153  : accessor(idxset, range_idx)
1154 {}
1155 
1156 
1157 
1159  : accessor(nullptr)
1160 {}
1161 
1162 
1163 
1165  : accessor(idxset)
1166 {}
1167 
1168 
1169 
1172 {
1173  accessor.advance();
1174  return *this;
1175 }
1176 
1177 
1178 
1181 {
1182  const IndexSet::IntervalIterator iter = *this;
1183  accessor.advance();
1184  return iter;
1185 }
1186 
1187 
1188 
1190  operator*() const
1191 {
1192  return accessor;
1193 }
1194 
1195 
1196 
1199 {
1200  return &accessor;
1201 }
1202 
1203 
1204 
1205 inline bool
1208 {
1209  return accessor == other.accessor;
1210 }
1211 
1212 
1213 
1214 inline bool
1217 {
1218  return !(*this == other);
1219 }
1220 
1221 
1222 
1223 inline bool
1226 {
1227  return accessor < other.accessor;
1228 }
1229 
1230 
1231 
1232 inline int
1235 {
1237  ExcMessage(
1238  "Can not compare iterators belonging to different IndexSets"));
1239 
1241  accessor.index_set->ranges.size() :
1243  const size_type rhs =
1245  accessor.index_set->ranges.size() :
1246  other.accessor.range_idx;
1247 
1248  if (lhs > rhs)
1249  return static_cast<int>(lhs - rhs);
1250  else
1251  return -static_cast<int>(rhs - lhs);
1252 }
1253 
1254 
1255 
1256 /* ElementIterator */
1257 
1259  const IndexSet * idxset,
1260  const IndexSet::size_type range_idx,
1261  const IndexSet::size_type index)
1262  : index_set(idxset)
1263  , range_idx(range_idx)
1264  , idx(index)
1265 {
1266  Assert(range_idx < index_set->ranges.size(),
1267  ExcMessage(
1268  "Invalid range index for IndexSet::ElementIterator constructor."));
1269  Assert(
1270  idx >= index_set->ranges[range_idx].begin &&
1271  idx < index_set->ranges[range_idx].end,
1273  "Invalid index argument for IndexSet::ElementIterator constructor."));
1274 }
1275 
1276 
1277 
1279  : index_set(idxset)
1280  , range_idx(numbers::invalid_dof_index)
1281  , idx(numbers::invalid_dof_index)
1282 {}
1283 
1284 
1285 
1286 inline bool
1288 {
1291  (range_idx < index_set->ranges.size() &&
1292  idx < index_set->ranges[range_idx].end),
1293  ExcInternalError("Invalid ElementIterator state."));
1294 
1295  return (range_idx < index_set->ranges.size() &&
1296  idx < index_set->ranges[range_idx].end);
1297 }
1298 
1299 
1300 
1302 {
1303  Assert(
1304  is_valid(),
1305  ExcMessage(
1306  "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1307  return idx;
1308 }
1309 
1310 
1311 
1312 inline bool
1315 {
1316  Assert(index_set == other.index_set,
1317  ExcMessage(
1318  "Can not compare iterators belonging to different IndexSets"));
1319  return range_idx == other.range_idx && idx == other.idx;
1320 }
1321 
1322 
1323 
1324 inline void
1326 {
1327  Assert(
1328  is_valid(),
1329  ExcMessage(
1330  "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1331  if (idx < index_set->ranges[range_idx].end)
1332  ++idx;
1333  // end of this range?
1334  if (idx == index_set->ranges[range_idx].end)
1335  {
1336  // point to first element in next interval if possible
1337  if (range_idx < index_set->ranges.size() - 1)
1338  {
1339  ++range_idx;
1340  idx = index_set->ranges[range_idx].begin;
1341  }
1342  else
1343  {
1344  // we just fell off the end, set to invalid:
1347  }
1348  }
1349 }
1350 
1351 
1352 
1355 {
1356  advance();
1357  return *this;
1358 }
1359 
1360 
1361 
1364 {
1365  const IndexSet::ElementIterator it = *this;
1366  advance();
1367  return it;
1368 }
1369 
1370 
1371 
1372 inline bool
1375 {
1376  return !(*this == other);
1377 }
1378 
1379 
1380 
1381 inline bool
1384 {
1385  Assert(index_set == other.index_set,
1386  ExcMessage(
1387  "Can not compare iterators belonging to different IndexSets"));
1388  return range_idx < other.range_idx ||
1389  (range_idx == other.range_idx && idx < other.idx);
1390 }
1391 
1392 
1393 
1394 inline std::ptrdiff_t
1397 {
1398  Assert(index_set == other.index_set,
1399  ExcMessage(
1400  "Can not compare iterators belonging to different IndexSets"));
1401  if (*this == other)
1402  return 0;
1403  if (!(*this < other))
1404  return -(other - *this);
1405 
1406  // only other can be equal to end() because of the checks above.
1408 
1409  // Note: we now compute how far advance *this in "*this < other" to get other,
1410  // so we need to return -c at the end.
1411 
1412  // first finish the current range:
1413  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1414 
1415  // now walk in steps of ranges (need to start one behind our current one):
1416  for (size_type range = range_idx + 1;
1417  range < index_set->ranges.size() && range <= other.range_idx;
1418  ++range)
1419  c += index_set->ranges[range].end - index_set->ranges[range].begin;
1420 
1421  Assert(
1422  other.range_idx < index_set->ranges.size() ||
1424  ExcMessage(
1425  "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1426 
1427  // We might have walked too far because we went until the end of
1428  // other.range_idx, so walk backwards to other.idx:
1430  c -= index_set->ranges[other.range_idx].end - other.idx;
1431 
1432  return -c;
1433 }
1434 
1435 
1436 /* Range */
1437 
1439  : begin(numbers::invalid_dof_index)
1440  , end(numbers::invalid_dof_index)
1441  , nth_index_in_set(numbers::invalid_dof_index)
1442 {}
1443 
1444 
1445 
1446 inline IndexSet::Range::Range(const size_type i1, const size_type i2)
1447  : begin(i1)
1448  , end(i2)
1449  , nth_index_in_set(numbers::invalid_dof_index)
1450 {}
1451 
1452 
1453 
1454 /* IndexSet itself */
1455 
1457  : is_compressed(true)
1458  , index_space_size(0)
1459  , largest_range(numbers::invalid_unsigned_int)
1460 {}
1461 
1462 
1463 
1465  : is_compressed(true)
1466  , index_space_size(size)
1467  , largest_range(numbers::invalid_unsigned_int)
1468 {}
1469 
1470 
1471 
1472 inline IndexSet::IndexSet(IndexSet &&is) noexcept
1473  : ranges(std::move(is.ranges))
1474  , is_compressed(is.is_compressed)
1475  , index_space_size(is.index_space_size)
1476  , largest_range(is.largest_range)
1477 {
1478  is.ranges.clear();
1479  is.is_compressed = true;
1480  is.index_space_size = 0;
1481  is.largest_range = numbers::invalid_unsigned_int;
1482 
1483  compress();
1484 }
1485 
1486 
1487 
1488 inline IndexSet &
1490 {
1491  ranges = std::move(is.ranges);
1492  is_compressed = is.is_compressed;
1493  index_space_size = is.index_space_size;
1494  largest_range = is.largest_range;
1495 
1496  is.ranges.clear();
1497  is.is_compressed = true;
1498  is.index_space_size = 0;
1499  is.largest_range = numbers::invalid_unsigned_int;
1500 
1501  compress();
1502 
1503  return *this;
1504 }
1505 
1506 
1507 
1510 {
1511  compress();
1512  if (ranges.size() > 0)
1513  return {this, 0, ranges[0].begin};
1514  else
1515  return end();
1516 }
1517 
1518 
1519 
1521 IndexSet::at(const size_type global_index) const
1522 {
1523  compress();
1524  AssertIndexRange(global_index, size());
1525 
1526  if (ranges.empty())
1527  return end();
1528 
1529  std::vector<Range>::const_iterator main_range =
1530  ranges.begin() + largest_range;
1531 
1532  Range r(global_index, global_index + 1);
1533  // This optimization makes the bounds for lower_bound smaller by checking
1534  // the largest range first.
1535  std::vector<Range>::const_iterator range_begin, range_end;
1536  if (global_index < main_range->begin)
1537  {
1538  range_begin = ranges.begin();
1539  range_end = main_range;
1540  }
1541  else
1542  {
1543  range_begin = main_range;
1544  range_end = ranges.end();
1545  }
1546 
1547  // This will give us the first range p=[a,b[ with b>=global_index using
1548  // a binary search
1549  const std::vector<Range>::const_iterator p =
1550  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1551 
1552  // We couldn't find a range, which means we have no range that contains
1553  // global_index and also no range behind it, meaning we need to return end().
1554  if (p == ranges.end())
1555  return end();
1556 
1557  // Finally, we can have two cases: Either global_index is not in [a,b[,
1558  // which means we need to return an iterator to a because global_index, ...,
1559  // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
1560  // [a,b[ and we will return an iterator pointing directly at global_index
1561  // (else branch).
1562  if (global_index < p->begin)
1563  return {this, static_cast<size_type>(p - ranges.begin()), p->begin};
1564  else
1565  return {this, static_cast<size_type>(p - ranges.begin()), global_index};
1566 }
1567 
1568 
1569 
1572 {
1573  compress();
1574  return IndexSet::ElementIterator(this);
1575 }
1576 
1577 
1578 
1581 {
1582  compress();
1583  if (ranges.size() > 0)
1584  return IndexSet::IntervalIterator(this, 0);
1585  else
1586  return end_intervals();
1587 }
1588 
1589 
1590 
1593 {
1594  compress();
1595  return IndexSet::IntervalIterator(this);
1596 }
1597 
1598 
1599 
1600 inline void
1602 {
1603  // reset so that there are no indices in the set any more; however,
1604  // as documented, the index set retains its size
1605  ranges.clear();
1606  is_compressed = true;
1608 }
1609 
1610 
1611 
1612 inline void
1614 {
1615  Assert(ranges.empty(),
1616  ExcMessage("This function can only be called if the current "
1617  "object does not yet contain any elements."));
1618  index_space_size = sz;
1619  is_compressed = true;
1620 }
1621 
1622 
1623 
1624 inline IndexSet::size_type
1626 {
1627  return index_space_size;
1628 }
1629 
1630 
1631 
1632 inline void
1634 {
1635  if (is_compressed == true)
1636  return;
1637 
1638  do_compress();
1639 }
1640 
1641 
1642 
1643 inline void
1645 {
1647 
1648  const Range new_range(index, index + 1);
1649  if (ranges.size() == 0 || index > ranges.back().end)
1650  ranges.push_back(new_range);
1651  else if (index == ranges.back().end)
1652  ranges.back().end++;
1653  else
1654  ranges.insert(Utilities::lower_bound(ranges.begin(),
1655  ranges.end(),
1656  new_range),
1657  new_range);
1658  is_compressed = false;
1659 }
1660 
1661 
1662 
1663 inline void
1665 {
1666  Assert((begin < index_space_size) ||
1667  ((begin == index_space_size) && (end == index_space_size)),
1668  ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1669  Assert(end <= index_space_size,
1670  ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1671  AssertIndexRange(begin, end + 1);
1672 
1673  if (begin != end)
1674  {
1675  const Range new_range(begin, end);
1676 
1677  // the new index might be larger than the last index present in the
1678  // ranges. Then we can skip the binary search
1679  if (ranges.size() == 0 || begin > ranges.back().end)
1680  ranges.push_back(new_range);
1681  else
1682  ranges.insert(Utilities::lower_bound(ranges.begin(),
1683  ranges.end(),
1684  new_range),
1685  new_range);
1686  is_compressed = false;
1687  }
1688 }
1689 
1690 
1691 
1692 template <typename ForwardIterator>
1693 inline void
1694 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1695 {
1696  if (begin == end)
1697  return;
1698 
1699  // identify ranges in the given iterator range by checking whether some
1700  // indices happen to be consecutive. to avoid quadratic complexity when
1701  // calling add_range many times (as add_range() going into the middle of an
1702  // already existing range must shift entries around), we first collect a
1703  // vector of ranges.
1704  std::vector<std::pair<size_type, size_type>> tmp_ranges;
1705  bool ranges_are_sorted = true;
1706  for (ForwardIterator p = begin; p != end;)
1707  {
1708  const size_type begin_index = *p;
1709  size_type end_index = begin_index + 1;
1710  ForwardIterator q = p;
1711  ++q;
1712  while ((q != end) && (*q == end_index))
1713  {
1714  ++end_index;
1715  ++q;
1716  }
1717 
1718  tmp_ranges.emplace_back(begin_index, end_index);
1719  p = q;
1720 
1721  // if the starting index of the next go-around of the for loop is less
1722  // than the end index of the one just identified, then we will have at
1723  // least one pair of ranges that are not sorted, and consequently the
1724  // whole collection of ranges is not sorted.
1725  if (p != end && *p < end_index)
1726  ranges_are_sorted = false;
1727  }
1728 
1729  if (!ranges_are_sorted)
1730  std::sort(tmp_ranges.begin(), tmp_ranges.end());
1731 
1732  // if we have many ranges, we first construct a temporary index set (where
1733  // we add ranges in a consecutive way, so fast), otherwise, we work with
1734  // add_range(). the number 9 is chosen heuristically given the fact that
1735  // there are typically up to 8 independent ranges when adding the degrees of
1736  // freedom on a 3D cell or 9 when adding degrees of freedom of faces. if
1737  // doing cell-by-cell additions, we want to avoid repeated calls to
1738  // IndexSet::compress() which gets called upon merging two index sets, so we
1739  // want to be in the other branch then.
1740  if (tmp_ranges.size() > 9)
1741  {
1742  IndexSet tmp_set(size());
1743  tmp_set.ranges.reserve(tmp_ranges.size());
1744  for (const auto &i : tmp_ranges)
1745  tmp_set.add_range(i.first, i.second);
1746  this->add_indices(tmp_set);
1747  }
1748  else
1749  for (const auto &i : tmp_ranges)
1750  add_range(i.first, i.second);
1751 }
1752 
1753 
1754 
1755 inline bool
1757 {
1758  if (ranges.empty() == false)
1759  {
1760  compress();
1761 
1762  // fast check whether the index is in the largest range
1764  if (index >= ranges[largest_range].begin &&
1765  index < ranges[largest_range].end)
1766  return true;
1767 
1768  // get the element after which we would have to insert a range that
1769  // consists of all elements from this element to the end of the index
1770  // range plus one. after this call we know that if p!=end() then
1771  // p->begin<=index unless there is no such range at all
1772  //
1773  // if the searched for element is an element of this range, then we're
1774  // done. otherwise, the element can't be in one of the following ranges
1775  // because otherwise p would be a different iterator
1776  //
1777  // since we already know the position relative to the largest range (we
1778  // called compress!), we can perform the binary search on ranges with
1779  // lower/higher number compared to the largest range
1780  std::vector<Range>::const_iterator p = std::upper_bound(
1781  ranges.begin() +
1782  (index < ranges[largest_range].begin ? 0 : largest_range + 1),
1783  index < ranges[largest_range].begin ? ranges.begin() + largest_range :
1784  ranges.end(),
1785  Range(index, size() + 1));
1786 
1787  if (p == ranges.begin())
1788  return ((index >= p->begin) && (index < p->end));
1789 
1790  Assert((p == ranges.end()) || (p->begin > index), ExcInternalError());
1791 
1792  // now move to that previous range
1793  --p;
1794  Assert(p->begin <= index, ExcInternalError());
1795 
1796  return (p->end > index);
1797  }
1798 
1799  // didn't find this index, so it's not in the set
1800  return false;
1801 }
1802 
1803 
1804 
1805 inline bool
1807 {
1808  compress();
1809  return (ranges.size() <= 1);
1810 }
1811 
1812 
1813 
1814 inline bool
1816 {
1817  return ranges.empty();
1818 }
1819 
1820 
1821 
1822 inline IndexSet::size_type
1824 {
1825  // make sure we have non-overlapping ranges
1826  compress();
1827 
1828  size_type v = 0;
1829  if (!ranges.empty())
1830  {
1831  Range &r = ranges.back();
1832  v = r.nth_index_in_set + r.end - r.begin;
1833  }
1834 
1835 #ifdef DEBUG
1836  size_type s = 0;
1837  for (const auto &range : ranges)
1838  s += (range.end - range.begin);
1839  Assert(s == v, ExcInternalError());
1840 #endif
1841 
1842  return v;
1843 }
1844 
1845 
1846 
1847 inline unsigned int
1849 {
1850  compress();
1851  return ranges.size();
1852 }
1853 
1854 
1855 
1856 inline IndexSet::size_type
1858 {
1859  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1860 
1861  compress();
1862  const std::vector<Range>::const_iterator main_range =
1863  ranges.begin() + largest_range;
1864 
1865  return main_range->nth_index_in_set;
1866 }
1867 
1868 
1869 
1870 inline IndexSet::size_type
1872 {
1874 
1875  compress();
1876 
1877  // first check whether the index is in the largest range
1879  std::vector<Range>::const_iterator main_range =
1880  ranges.begin() + largest_range;
1881  if (n >= main_range->nth_index_in_set &&
1882  n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1883  return main_range->begin + (n - main_range->nth_index_in_set);
1884 
1885  // find out which chunk the local index n belongs to by using a binary
1886  // search. the comparator is based on the end of the ranges. Use the
1887  // position relative to main_range to subdivide the ranges
1888  Range r(n, n + 1);
1889  r.nth_index_in_set = n;
1890  std::vector<Range>::const_iterator range_begin, range_end;
1891  if (n < main_range->nth_index_in_set)
1892  {
1893  range_begin = ranges.begin();
1894  range_end = main_range;
1895  }
1896  else
1897  {
1898  range_begin = main_range + 1;
1899  range_end = ranges.end();
1900  }
1901 
1902  const std::vector<Range>::const_iterator p =
1903  Utilities::lower_bound(range_begin, range_end, r, Range::nth_index_compare);
1904 
1905  Assert(p != ranges.end(), ExcInternalError());
1906  return p->begin + (n - p->nth_index_in_set);
1907 }
1908 
1909 
1910 
1911 inline IndexSet::size_type
1913 {
1914  // to make this call thread-safe, compress() must not be called through this
1915  // function
1916  Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1917  AssertIndexRange(n, size());
1918 
1919  // return immediately if the index set is empty
1920  if (is_empty())
1922 
1923  // check whether the index is in the largest range. use the result to
1924  // perform a one-sided binary search afterward
1926  std::vector<Range>::const_iterator main_range =
1927  ranges.begin() + largest_range;
1928  if (n >= main_range->begin && n < main_range->end)
1929  return (n - main_range->begin) + main_range->nth_index_in_set;
1930 
1931  Range r(n, n);
1932  std::vector<Range>::const_iterator range_begin, range_end;
1933  if (n < main_range->begin)
1934  {
1935  range_begin = ranges.begin();
1936  range_end = main_range;
1937  }
1938  else
1939  {
1940  range_begin = main_range + 1;
1941  range_end = ranges.end();
1942  }
1943 
1944  std::vector<Range>::const_iterator p =
1945  Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1946 
1947  // if n is not in this set
1948  if (p == range_end || p->end == n || p->begin > n)
1950 
1951  Assert(p != ranges.end(), ExcInternalError());
1952  Assert(p->begin <= n, ExcInternalError());
1953  Assert(n < p->end, ExcInternalError());
1954  return (n - p->begin) + p->nth_index_in_set;
1955 }
1956 
1957 
1958 
1959 inline bool
1961 {
1962  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1963 
1964  compress();
1965  is.compress();
1966 
1967  return ranges == is.ranges;
1968 }
1969 
1970 
1971 
1972 inline bool
1974 {
1975  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1976 
1977  compress();
1978  is.compress();
1979 
1980  return ranges != is.ranges;
1981 }
1982 
1983 
1984 
1985 template <typename Vector>
1986 void
1987 IndexSet::fill_binary_vector(Vector &vector) const
1988 {
1989  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
1990 
1991  compress();
1992  // first fill all elements of the vector with zeroes.
1993  std::fill(vector.begin(), vector.end(), 0);
1994 
1995  // then write ones into the elements whose indices are contained in the
1996  // index set
1997  for (const auto &range : ranges)
1998  for (size_type i = range.begin; i < range.end; ++i)
1999  vector[i] = 1;
2000 }
2001 
2002 
2003 
2004 template <class StreamType>
2005 inline void
2006 IndexSet::print(StreamType &out) const
2007 {
2008  compress();
2009  out << "{";
2010  std::vector<Range>::const_iterator p;
2011  for (p = ranges.begin(); p != ranges.end(); ++p)
2012  {
2013  if (p->end - p->begin == 1)
2014  out << p->begin;
2015  else
2016  out << "[" << p->begin << "," << p->end - 1 << "]";
2017 
2018  if (p != --ranges.end())
2019  out << ", ";
2020  }
2021  out << "}" << std::endl;
2022 }
2023 
2024 
2025 
2026 template <class Archive>
2027 inline void
2028 IndexSet::Range::serialize(Archive &ar, const unsigned int)
2029 {
2030  ar &begin &end &nth_index_in_set;
2031 }
2032 
2033 
2034 
2035 template <class Archive>
2036 inline void
2037 IndexSet::serialize(Archive &ar, const unsigned int)
2038 {
2040 }
2041 
2042 DEAL_II_NAMESPACE_CLOSE
2043 
2044 #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:1571
void add_index(const size_type index)
Definition: index_set.h:1644
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:1198
const IntervalAccessor & operator*() const
Definition: index_set.h:1190
const IndexSet * index_set
Definition: index_set.h:799
IndexSet & operator=(const IndexSet &)=default
void block_read(std::istream &in)
Definition: index_set.cc:467
size_type largest_range_starting_index() const
Definition: index_set.h:1857
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1396
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1694
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:1258
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1225
std::forward_iterator_tag iterator_category
Definition: index_set.h:692
std::size_t memory_consumption() const
Definition: index_set.cc:709
LinearAlgebra::distributed::Vector< Number > Vector
size_type n_elements() const
Definition: index_set.h:1048
size_type index_space_size
Definition: index_set.h:959
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:646
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:576
size_type size() const
Definition: index_set.h:1625
void print(StreamType &out) const
Definition: index_set.h:2006
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2028
void clear()
Definition: index_set.h:1601
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2037
static ::ExceptionBase & ExcMessage(std::string arg1)
void block_write(std::ostream &out) const
Definition: index_set.cc:453
std::vector< Range > ranges
Definition: index_set.h:943
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:783
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
bool operator<(const ElementIterator &) const
Definition: index_set.h:1383
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1912
bool operator==(const IndexSet &is) const
Definition: index_set.h:1960
signed int value_type
Definition: index_set.h:104
IntervalAccessor accessor
Definition: index_set.h:702
void do_compress() const
Definition: index_set.cc:98
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1098
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1207
size_type last() const
Definition: index_set.h:1088
const IndexSet * index_set
Definition: index_set.h:593
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1111
unsigned int n_intervals() const
Definition: index_set.h:1848
IndexSet tensor_product(const IndexSet &other) const
Definition: index_set.cc:310
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1234
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1973
size_type pop_back()
Definition: index_set.cc:322
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:487
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1664
IndexSet operator &(const IndexSet &is) const
void set_size(const size_type size)
Definition: index_set.h:1613
ElementIterator & operator++()
Definition: index_set.h:1354
unsigned int global_dof_index
Definition: types.h:89
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:1017
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:209
void compress() const
Definition: index_set.h:1633
ElementIterator end() const
Definition: index_set.h:1074
IntervalIterator begin_intervals() const
Definition: index_set.h:1580
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1004
bool is_contiguous() const
Definition: index_set.h:1806
Threads::Mutex compress_mutex
Definition: index_set.h:976
size_type pop_front()
Definition: index_set.cc:340
IntervalIterator end_intervals() const
Definition: index_set.h:1592
bool is_compressed
Definition: index_set.h:953
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1123
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1374
ElementIterator begin() const
Definition: index_set.h:1065
void write(std::ostream &out) const
Definition: index_set.cc:416
bool is_element(const size_type index) const
Definition: index_set.h:1756
bool operator==(const ElementIterator &) const
Definition: index_set.h:1314
ElementIterator at(const size_type global_index) const
Definition: index_set.h:1521
const types::global_dof_index invalid_dof_index
Definition: types.h:202
ElementIterator begin() const
Definition: index_set.h:1509
void read(std::istream &in)
Definition: index_set.cc:431
size_type n_elements() const
Definition: index_set.h:1823
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1871
bool is_empty() const
Definition: index_set.h:1815
size_type largest_range
Definition: index_set.h:970
IntervalIterator & operator++()
Definition: index_set.h:1171
static ::ExceptionBase & ExcInternalError()
size_type operator*() const
Definition: index_set.h:1301
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1216