Reference documentation for deal.II version 9.4.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
index_set.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2009 - 2022 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
24
25#include <algorithm>
26#include <vector>
27
28
29#ifdef DEAL_II_WITH_TRILINOS
30# include <Epetra_Map.h>
31# ifdef DEAL_II_TRILINOS_WITH_TPETRA
32# include <Tpetra_Map.hpp>
33# endif
34#endif
35
36#if defined(DEAL_II_WITH_MPI) || defined(DEAL_II_WITH_PETSC)
37# include <mpi.h>
38#else
39using MPI_Comm = int;
40# ifndef MPI_COMM_WORLD
41# define MPI_COMM_WORLD 0
42# endif
43#endif
44
46
70{
71public:
72 // forward declarations:
73 class ElementIterator;
74 class IntervalIterator;
75
81
99 using value_type = signed int;
100
101
105 IndexSet();
106
110 explicit IndexSet(const size_type size);
111
115 IndexSet(const IndexSet &) = default;
116
120 IndexSet &
121 operator=(const IndexSet &) = default;
122
127 IndexSet(IndexSet &&is) noexcept;
128
133 IndexSet &
134 operator=(IndexSet &&is) noexcept;
135
136#ifdef DEAL_II_WITH_TRILINOS
140 explicit IndexSet(const Epetra_BlockMap &map);
141#endif
142
147 void
148 clear();
149
156 void
157 set_size(const size_type size);
158
167 size() const;
168
175 void
176 add_range(const size_type begin, const size_type end);
177
181 void
182 add_index(const size_type index);
183
193 template <typename ForwardIterator>
194 void
195 add_indices(const ForwardIterator &begin, const ForwardIterator &end);
196
212 void
213 add_indices(const IndexSet &other, const size_type offset = 0);
214
218 bool
219 is_element(const size_type index) const;
220
225 bool
226 is_contiguous() const;
227
232 bool
233 is_empty() const;
234
244 bool
245 is_ascending_and_one_to_one(const MPI_Comm &communicator) const;
246
251 n_elements() const;
252
259 nth_index_in_set(const size_type local_index) const;
260
269
278 unsigned int
279 n_intervals() const;
280
294
299 void
300 compress() const;
301
307 bool
308 operator==(const IndexSet &is) const;
309
315 bool
316 operator!=(const IndexSet &is) const;
317
325 operator&(const IndexSet &is) const;
326
341 get_view(const size_type begin, const size_type end) const;
342
348 std::vector<IndexSet>
350 const std::vector<types::global_dof_index> &n_indices_per_block) const;
351
357 void
358 subtract_set(const IndexSet &other);
359
377 tensor_product(const IndexSet &other) const;
378
384 pop_back();
385
391 pop_front();
392
396 void
397 fill_index_vector(std::vector<size_type> &indices) const;
398
411 template <typename VectorType>
412 void
413 fill_binary_vector(VectorType &vector) const;
414
419 template <class StreamType>
420 void
421 print(StreamType &out) const;
422
427 void
428 write(std::ostream &out) const;
429
434 void
435 read(std::istream &in);
436
441 void
442 block_write(std::ostream &out) const;
443
448 void
449 block_read(std::istream &in);
450
451#ifdef DEAL_II_WITH_TRILINOS
480 Epetra_Map
481 make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
482 const bool overlapping = false) const;
483
484# ifdef DEAL_II_TRILINOS_WITH_TPETRA
485 Tpetra::Map<int, types::global_dof_index>
486 make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
487 const bool overlapping = false) const;
488# endif
489#endif
490
491
496 std::size_t
497 memory_consumption() const;
498
500 size_type,
501 << "The global index " << arg1
502 << " is not an element of this set.");
503
509 template <class Archive>
510 void
511 serialize(Archive &ar, const unsigned int version);
512
513
525 {
526 public:
531 IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
532
536 explicit IntervalAccessor(const IndexSet *idxset);
537
542 n_elements() const;
543
547 bool
548 is_valid() const;
549
554 begin() const;
555
561 end() const;
562
567 last() const;
568
569 private:
578 operator=(const IntervalAccessor &other);
579
583 bool
584 operator==(const IntervalAccessor &other) const;
588 bool
589 operator<(const IntervalAccessor &other) const;
594 void
595 advance();
600
606
607 friend class IntervalIterator;
608 };
609
615 {
616 public:
621 IntervalIterator(const IndexSet *idxset, const size_type range_idx);
622
626 explicit IntervalIterator(const IndexSet *idxset);
627
632
636 IntervalIterator(const IntervalIterator &other) = default;
637
642 operator=(const IntervalIterator &other) = default;
643
648 operator++();
649
654 operator++(int);
655
659 const IntervalAccessor &
660 operator*() const;
661
665 const IntervalAccessor *
666 operator->() const;
667
671 bool
672 operator==(const IntervalIterator &) const;
673
677 bool
678 operator!=(const IntervalIterator &) const;
679
683 bool
684 operator<(const IntervalIterator &) const;
685
692 int
693 operator-(const IntervalIterator &p) const;
694
700 using iterator_category = std::forward_iterator_tag;
702 using difference_type = std::ptrdiff_t;
705
706 private:
711 };
712
718 {
719 public:
724 ElementIterator(const IndexSet *idxset,
725 const size_type range_idx,
726 const size_type index);
727
731 explicit ElementIterator(const IndexSet *idxset);
732
738 operator*() const;
739
743 bool
744 is_valid() const;
745
750 operator++();
751
756 operator++(int);
757
761 bool
762 operator==(const ElementIterator &) const;
763
767 bool
768 operator!=(const ElementIterator &) const;
769
773 bool
774 operator<(const ElementIterator &) const;
775
784 std::ptrdiff_t
785 operator-(const ElementIterator &p) const;
786
792 using iterator_category = std::forward_iterator_tag;
794 using difference_type = std::ptrdiff_t;
797
798 private:
802 void
803 advance();
804
817 };
818
824 begin() const;
825
841 at(const size_type global_index) const;
842
848 end() const;
849
854 begin_intervals() const;
855
861 end_intervals() const;
862
867private:
876 struct Range
877 {
880
882
891 Range();
892
900 Range(const size_type i1, const size_type i2);
901
902 friend inline bool
903 operator<(const Range &range_1, const Range &range_2)
904 {
905 return (
906 (range_1.begin < range_2.begin) ||
907 ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
908 }
909
910 static bool
912 {
913 return x.end < y.end;
914 }
915
916 static bool
918 {
919 return (x.nth_index_in_set + (x.end - x.begin) <
920 y.nth_index_in_set + (y.end - y.begin));
921 }
922
923 friend inline bool
924 operator==(const Range &range_1, const Range &range_2)
925 {
926 return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
927 }
928
929 static std::size_t
931 {
932 return sizeof(Range);
933 }
934
940 template <class Archive>
941 void
942 serialize(Archive &ar, const unsigned int version);
943 };
944
953 mutable std::vector<Range> ranges;
954
963 mutable bool is_compressed;
964
970
981
987
991 void
992 do_compress() const;
993};
994
995
1013inline IndexSet
1015{
1016 IndexSet is(N);
1017 is.add_range(0, N);
1018 is.compress();
1019 return is;
1020}
1021
1022/* ------------------ inline functions ------------------ */
1023
1024
1025/* IntervalAccessor */
1026
1028 const IndexSet * idxset,
1029 const IndexSet::size_type range_idx)
1030 : index_set(idxset)
1031 , range_idx(range_idx)
1032{
1033 Assert(range_idx < idxset->n_intervals(),
1034 ExcInternalError("Invalid range index"));
1035}
1036
1037
1038
1040 : index_set(idxset)
1041 , range_idx(numbers::invalid_dof_index)
1042{}
1043
1044
1045
1047 const IndexSet::IntervalAccessor &other)
1048 : index_set(other.index_set)
1049 , range_idx(other.range_idx)
1050{
1052 ExcMessage("invalid iterator"));
1053}
1054
1055
1056
1059{
1060 Assert(is_valid(), ExcMessage("invalid iterator"));
1061 return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1062}
1063
1064
1065
1066inline bool
1068{
1069 return index_set != nullptr && range_idx < index_set->n_intervals();
1070}
1071
1072
1073
1076{
1077 Assert(is_valid(), ExcMessage("invalid iterator"));
1078 return {index_set, range_idx, index_set->ranges[range_idx].begin};
1079}
1080
1081
1082
1085{
1086 Assert(is_valid(), ExcMessage("invalid iterator"));
1087
1088 // point to first index in next interval unless we are the last interval.
1089 if (range_idx < index_set->ranges.size() - 1)
1090 return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1091 else
1092 return index_set->end();
1093}
1094
1095
1096
1099{
1100 Assert(is_valid(), ExcMessage("invalid iterator"));
1101
1102 return index_set->ranges[range_idx].end - 1;
1103}
1104
1105
1106
1109{
1110 index_set = other.index_set;
1111 range_idx = other.range_idx;
1112 Assert(range_idx == numbers::invalid_dof_index || is_valid(),
1113 ExcMessage("invalid iterator"));
1114 return *this;
1115}
1116
1117
1118
1119inline bool
1121 const IndexSet::IntervalAccessor &other) const
1122{
1123 Assert(index_set == other.index_set,
1124 ExcMessage(
1125 "Can not compare accessors pointing to different IndexSets"));
1126 return range_idx == other.range_idx;
1127}
1128
1129
1130
1131inline bool
1133 const IndexSet::IntervalAccessor &other) const
1134{
1135 Assert(index_set == other.index_set,
1136 ExcMessage(
1137 "Can not compare accessors pointing to different IndexSets"));
1138 return range_idx < other.range_idx;
1139}
1140
1141
1142
1143inline void
1145{
1146 Assert(
1147 is_valid(),
1148 ExcMessage(
1149 "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1150 ++range_idx;
1151
1152 // set ourselves to invalid if we walk off the end
1153 if (range_idx >= index_set->ranges.size())
1154 range_idx = numbers::invalid_dof_index;
1155}
1156
1157
1158/* IntervalIterator */
1159
1161 const IndexSet * idxset,
1162 const IndexSet::size_type range_idx)
1163 : accessor(idxset, range_idx)
1164{}
1165
1166
1167
1169 : accessor(nullptr)
1170{}
1171
1172
1173
1175 : accessor(idxset)
1176{}
1177
1178
1179
1182{
1183 accessor.advance();
1184 return *this;
1185}
1186
1187
1188
1191{
1192 const IndexSet::IntervalIterator iter = *this;
1193 accessor.advance();
1194 return iter;
1195}
1196
1197
1198
1199inline const IndexSet::IntervalAccessor &
1201{
1202 return accessor;
1203}
1204
1205
1206
1207inline const IndexSet::IntervalAccessor *
1209{
1210 return &accessor;
1211}
1212
1213
1214
1215inline bool
1217 const IndexSet::IntervalIterator &other) const
1218{
1219 return accessor == other.accessor;
1220}
1221
1222
1223
1224inline bool
1226 const IndexSet::IntervalIterator &other) const
1227{
1228 return !(*this == other);
1229}
1230
1231
1232
1233inline bool
1235 const IndexSet::IntervalIterator &other) const
1236{
1237 return accessor < other.accessor;
1238}
1239
1240
1241
1242inline int
1244 const IndexSet::IntervalIterator &other) const
1245{
1246 Assert(accessor.index_set == other.accessor.index_set,
1247 ExcMessage(
1248 "Can not compare iterators belonging to different IndexSets"));
1249
1250 const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index) ?
1251 accessor.index_set->ranges.size() :
1252 accessor.range_idx;
1253 const size_type rhs =
1255 accessor.index_set->ranges.size() :
1256 other.accessor.range_idx;
1257
1258 if (lhs > rhs)
1259 return static_cast<int>(lhs - rhs);
1260 else
1261 return -static_cast<int>(rhs - lhs);
1262}
1263
1264
1265
1266/* ElementIterator */
1267
1269 const IndexSet * idxset,
1270 const IndexSet::size_type range_idx,
1271 const IndexSet::size_type index)
1272 : index_set(idxset)
1273 , range_idx(range_idx)
1274 , idx(index)
1275{
1276 Assert(range_idx < index_set->ranges.size(),
1277 ExcMessage(
1278 "Invalid range index for IndexSet::ElementIterator constructor."));
1279 Assert(
1280 idx >= index_set->ranges[range_idx].begin &&
1281 idx < index_set->ranges[range_idx].end,
1283 "Invalid index argument for IndexSet::ElementIterator constructor."));
1284}
1285
1286
1287
1289 : index_set(idxset)
1290 , range_idx(numbers::invalid_dof_index)
1292{}
1293
1294
1295
1296inline bool
1298{
1299 Assert((range_idx == numbers::invalid_dof_index &&
1301 (range_idx < index_set->ranges.size() &&
1302 idx < index_set->ranges[range_idx].end),
1303 ExcInternalError("Invalid ElementIterator state."));
1304
1305 return (range_idx < index_set->ranges.size() &&
1306 idx < index_set->ranges[range_idx].end);
1307}
1308
1309
1310
1313{
1314 Assert(
1315 is_valid(),
1316 ExcMessage(
1317 "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1318 return idx;
1319}
1320
1321
1322
1323inline bool
1325 const IndexSet::ElementIterator &other) const
1326{
1327 Assert(index_set == other.index_set,
1328 ExcMessage(
1329 "Can not compare iterators belonging to different IndexSets"));
1330 return range_idx == other.range_idx && idx == other.idx;
1331}
1332
1333
1334
1335inline void
1337{
1338 Assert(
1339 is_valid(),
1340 ExcMessage(
1341 "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1342 if (idx < index_set->ranges[range_idx].end)
1343 ++idx;
1344 // end of this range?
1345 if (idx == index_set->ranges[range_idx].end)
1346 {
1347 // point to first element in next interval if possible
1348 if (range_idx < index_set->ranges.size() - 1)
1349 {
1350 ++range_idx;
1351 idx = index_set->ranges[range_idx].begin;
1352 }
1353 else
1354 {
1355 // we just fell off the end, set to invalid:
1356 range_idx = numbers::invalid_dof_index;
1358 }
1359 }
1360}
1361
1362
1363
1366{
1367 advance();
1368 return *this;
1369}
1370
1371
1372
1375{
1376 const IndexSet::ElementIterator it = *this;
1377 advance();
1378 return it;
1379}
1380
1381
1382
1383inline bool
1385 const IndexSet::ElementIterator &other) const
1386{
1387 return !(*this == other);
1388}
1389
1390
1391
1392inline bool
1394 const IndexSet::ElementIterator &other) const
1395{
1396 Assert(index_set == other.index_set,
1397 ExcMessage(
1398 "Can not compare iterators belonging to different IndexSets"));
1399 return range_idx < other.range_idx ||
1400 (range_idx == other.range_idx && idx < other.idx);
1401}
1402
1403
1404
1405inline std::ptrdiff_t
1407 const IndexSet::ElementIterator &other) const
1408{
1409 Assert(index_set == other.index_set,
1410 ExcMessage(
1411 "Can not compare iterators belonging to different IndexSets"));
1412 if (*this == other)
1413 return 0;
1414 if (!(*this < other))
1415 return -(other - *this);
1416
1417 // only other can be equal to end() because of the checks above.
1418 Assert(is_valid(), ExcInternalError());
1419
1420 // Note: we now compute how far advance *this in "*this < other" to get other,
1421 // so we need to return -c at the end.
1422
1423 // first finish the current range:
1424 std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1425
1426 // now walk in steps of ranges (need to start one behind our current one):
1427 for (size_type range = range_idx + 1;
1428 range < index_set->ranges.size() && range <= other.range_idx;
1429 ++range)
1430 c += index_set->ranges[range].end - index_set->ranges[range].begin;
1431
1432 Assert(
1433 other.range_idx < index_set->ranges.size() ||
1435 ExcMessage(
1436 "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1437
1438 // We might have walked too far because we went until the end of
1439 // other.range_idx, so walk backwards to other.idx:
1441 c -= index_set->ranges[other.range_idx].end - other.idx;
1442
1443 return -c;
1444}
1445
1446
1447/* Range */
1448
1453{}
1454
1455
1456
1458 : begin(i1)
1459 , end(i2)
1461{}
1462
1463
1464
1465/* IndexSet itself */
1466
1468 : is_compressed(true)
1469 , index_space_size(0)
1471{}
1472
1473
1474
1476 : is_compressed(true)
1477 , index_space_size(size)
1478 , largest_range(numbers::invalid_unsigned_int)
1479{}
1480
1481
1482
1483inline IndexSet::IndexSet(IndexSet &&is) noexcept
1484 : ranges(std::move(is.ranges))
1485 , is_compressed(is.is_compressed)
1486 , index_space_size(is.index_space_size)
1487 , largest_range(is.largest_range)
1488{
1489 is.ranges.clear();
1490 is.is_compressed = true;
1491 is.index_space_size = 0;
1492 is.largest_range = numbers::invalid_unsigned_int;
1493
1494 compress();
1495}
1496
1497
1498
1499inline IndexSet &
1501{
1502 ranges = std::move(is.ranges);
1503 is_compressed = is.is_compressed;
1504 index_space_size = is.index_space_size;
1505 largest_range = is.largest_range;
1506
1507 is.ranges.clear();
1508 is.is_compressed = true;
1509 is.index_space_size = 0;
1510 is.largest_range = numbers::invalid_unsigned_int;
1511
1512 compress();
1513
1514 return *this;
1515}
1516
1517
1518
1521{
1522 compress();
1523 if (ranges.size() > 0)
1524 return {this, 0, ranges[0].begin};
1525 else
1526 return end();
1527}
1528
1529
1530
1533{
1534 compress();
1536
1537 if (ranges.empty())
1538 return end();
1539
1540 std::vector<Range>::const_iterator main_range =
1541 ranges.begin() + largest_range;
1542
1544 // This optimization makes the bounds for lower_bound smaller by checking
1545 // the largest range first.
1546 std::vector<Range>::const_iterator range_begin, range_end;
1547 if (global_index < main_range->begin)
1548 {
1549 range_begin = ranges.begin();
1550 range_end = main_range;
1551 }
1552 else
1553 {
1554 range_begin = main_range;
1555 range_end = ranges.end();
1556 }
1557
1558 // This will give us the first range p=[a,b[ with b>=global_index using
1559 // a binary search
1560 const std::vector<Range>::const_iterator p =
1561 Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1562
1563 // We couldn't find a range, which means we have no range that contains
1564 // global_index and also no range behind it, meaning we need to return end().
1565 if (p == ranges.end())
1566 return end();
1567
1568 // Finally, we can have two cases: Either global_index is not in [a,b[,
1569 // which means we need to return an iterator to a because global_index, ...,
1570 // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
1571 // [a,b[ and we will return an iterator pointing directly at global_index
1572 // (else branch).
1573 if (global_index < p->begin)
1574 return {this, static_cast<size_type>(p - ranges.begin()), p->begin};
1575 else
1576 return {this, static_cast<size_type>(p - ranges.begin()), global_index};
1577}
1578
1579
1580
1583{
1584 compress();
1585 return IndexSet::ElementIterator(this);
1586}
1587
1588
1589
1592{
1593 compress();
1594 if (ranges.size() > 0)
1595 return IndexSet::IntervalIterator(this, 0);
1596 else
1597 return end_intervals();
1598}
1599
1600
1601
1604{
1605 compress();
1606 return IndexSet::IntervalIterator(this);
1607}
1608
1609
1610
1611inline void
1613{
1614 // reset so that there are no indices in the set any more; however,
1615 // as documented, the index set retains its size
1616 ranges.clear();
1617 is_compressed = true;
1619}
1620
1621
1622
1623inline void
1625{
1626 Assert(ranges.empty(),
1627 ExcMessage("This function can only be called if the current "
1628 "object does not yet contain any elements."));
1629 index_space_size = sz;
1630 is_compressed = true;
1631}
1632
1633
1634
1637{
1638 return index_space_size;
1639}
1640
1641
1642
1643inline void
1645{
1646 if (is_compressed == true)
1647 return;
1648
1649 do_compress();
1650}
1651
1652
1653
1654inline void
1656{
1658
1659 const Range new_range(index, index + 1);
1660 if (ranges.size() == 0 || index > ranges.back().end)
1661 ranges.push_back(new_range);
1662 else if (index == ranges.back().end)
1663 ranges.back().end++;
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
1674inline void
1676{
1679 ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1681 ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1683
1684 if (begin != end)
1685 {
1686 const Range new_range(begin, end);
1687
1688 // the new index might be larger than the last index present in the
1689 // ranges. Then we can skip the binary search
1690 if (ranges.size() == 0 || begin > ranges.back().end)
1691 ranges.push_back(new_range);
1692 else
1693 ranges.insert(Utilities::lower_bound(ranges.begin(),
1694 ranges.end(),
1695 new_range),
1696 new_range);
1697 is_compressed = false;
1698 }
1699}
1700
1701
1702
1703template <typename ForwardIterator>
1704inline void
1705IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1706{
1707 if (begin == end)
1708 return;
1709
1710 // identify ranges in the given iterator range by checking whether some
1711 // indices happen to be consecutive. to avoid quadratic complexity when
1712 // calling add_range many times (as add_range() going into the middle of an
1713 // already existing range must shift entries around), we first collect a
1714 // vector of ranges.
1715 std::vector<std::pair<size_type, size_type>> tmp_ranges;
1716 bool ranges_are_sorted = true;
1717 for (ForwardIterator p = begin; p != end;)
1718 {
1719 const size_type begin_index = *p;
1720 size_type end_index = begin_index + 1;
1721 ForwardIterator q = p;
1722 ++q;
1723 while ((q != end) && (static_cast<size_type>(*q) == end_index))
1724 {
1725 ++end_index;
1726 ++q;
1727 }
1728
1729 tmp_ranges.emplace_back(begin_index, end_index);
1730 p = q;
1731
1732 // if the starting index of the next go-around of the for loop is less
1733 // than the end index of the one just identified, then we will have at
1734 // least one pair of ranges that are not sorted, and consequently the
1735 // whole collection of ranges is not sorted.
1736 if (p != end && static_cast<size_type>(*p) < end_index)
1737 ranges_are_sorted = false;
1738 }
1739
1740 if (!ranges_are_sorted)
1741 std::sort(tmp_ranges.begin(), tmp_ranges.end());
1742
1743 // if we have many ranges, we first construct a temporary index set (where
1744 // we add ranges in a consecutive way, so fast), otherwise, we work with
1745 // add_range(). the number 9 is chosen heuristically given the fact that
1746 // there are typically up to 8 independent ranges when adding the degrees of
1747 // freedom on a 3D cell or 9 when adding degrees of freedom of faces. if
1748 // doing cell-by-cell additions, we want to avoid repeated calls to
1749 // IndexSet::compress() which gets called upon merging two index sets, so we
1750 // want to be in the other branch then.
1751 if (tmp_ranges.size() > 9)
1752 {
1753 IndexSet tmp_set(size());
1754 tmp_set.ranges.reserve(tmp_ranges.size());
1755 for (const auto &i : tmp_ranges)
1756 tmp_set.add_range(i.first, i.second);
1757 this->add_indices(tmp_set);
1758 }
1759 else
1760 for (const auto &i : tmp_ranges)
1761 add_range(i.first, i.second);
1762}
1763
1764
1765
1766inline bool
1768{
1769 if (ranges.empty() == false)
1770 {
1771 compress();
1772
1773 // fast check whether the index is in the largest range
1775 if (index >= ranges[largest_range].begin &&
1776 index < ranges[largest_range].end)
1777 return true;
1778
1779 // get the element after which we would have to insert a range that
1780 // consists of all elements from this element to the end of the index
1781 // range plus one. after this call we know that if p!=end() then
1782 // p->begin<=index unless there is no such range at all
1783 //
1784 // if the searched for element is an element of this range, then we're
1785 // done. otherwise, the element can't be in one of the following ranges
1786 // because otherwise p would be a different iterator
1787 //
1788 // since we already know the position relative to the largest range (we
1789 // called compress!), we can perform the binary search on ranges with
1790 // lower/higher number compared to the largest range
1791 std::vector<Range>::const_iterator p = std::upper_bound(
1792 ranges.begin() +
1793 (index < ranges[largest_range].begin ? 0 : largest_range + 1),
1794 index < ranges[largest_range].begin ? ranges.begin() + largest_range :
1795 ranges.end(),
1796 Range(index, size() + 1));
1797
1798 if (p == ranges.begin())
1799 return ((index >= p->begin) && (index < p->end));
1800
1801 Assert((p == ranges.end()) || (p->begin > index), ExcInternalError());
1802
1803 // now move to that previous range
1804 --p;
1805 Assert(p->begin <= index, ExcInternalError());
1806
1807 return (p->end > index);
1808 }
1809
1810 // didn't find this index, so it's not in the set
1811 return false;
1812}
1813
1814
1815
1816inline bool
1818{
1819 compress();
1820 return (ranges.size() <= 1);
1821}
1822
1823
1824
1825inline bool
1827{
1828 return ranges.empty();
1829}
1830
1831
1832
1835{
1836 // make sure we have non-overlapping ranges
1837 compress();
1838
1839 size_type v = 0;
1840 if (!ranges.empty())
1841 {
1842 Range &r = ranges.back();
1843 v = r.nth_index_in_set + r.end - r.begin;
1844 }
1845
1846#ifdef DEBUG
1847 size_type s = 0;
1848 for (const auto &range : ranges)
1849 s += (range.end - range.begin);
1850 Assert(s == v, ExcInternalError());
1851#endif
1852
1853 return v;
1854}
1855
1856
1857
1858inline unsigned int
1860{
1861 compress();
1862 return ranges.size();
1863}
1864
1865
1866
1869{
1870 Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1871
1872 compress();
1873 const std::vector<Range>::const_iterator main_range =
1874 ranges.begin() + largest_range;
1875
1876 return main_range->nth_index_in_set;
1877}
1878
1879
1880
1883{
1885
1886 compress();
1887
1888 // first check whether the index is in the largest range
1890 std::vector<Range>::const_iterator main_range =
1891 ranges.begin() + largest_range;
1892 if (n >= main_range->nth_index_in_set &&
1893 n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1894 return main_range->begin + (n - main_range->nth_index_in_set);
1895
1896 // find out which chunk the local index n belongs to by using a binary
1897 // search. the comparator is based on the end of the ranges. Use the
1898 // position relative to main_range to subdivide the ranges
1899 Range r(n, n + 1);
1900 r.nth_index_in_set = n;
1901 std::vector<Range>::const_iterator range_begin, range_end;
1902 if (n < main_range->nth_index_in_set)
1903 {
1904 range_begin = ranges.begin();
1905 range_end = main_range;
1906 }
1907 else
1908 {
1909 range_begin = main_range + 1;
1910 range_end = ranges.end();
1911 }
1912
1913 const std::vector<Range>::const_iterator p =
1914 Utilities::lower_bound(range_begin, range_end, r, Range::nth_index_compare);
1915
1916 Assert(p != ranges.end(), ExcInternalError());
1917 return p->begin + (n - p->nth_index_in_set);
1918}
1919
1920
1921
1924{
1925 // to make this call thread-safe, compress() must not be called through this
1926 // function
1927 Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1928 AssertIndexRange(n, size());
1929
1930 // return immediately if the index set is empty
1931 if (is_empty())
1933
1934 // check whether the index is in the largest range. use the result to
1935 // perform a one-sided binary search afterward
1937 std::vector<Range>::const_iterator main_range =
1938 ranges.begin() + largest_range;
1939 if (n >= main_range->begin && n < main_range->end)
1940 return (n - main_range->begin) + main_range->nth_index_in_set;
1941
1942 Range r(n, n);
1943 std::vector<Range>::const_iterator range_begin, range_end;
1944 if (n < main_range->begin)
1945 {
1946 range_begin = ranges.begin();
1947 range_end = main_range;
1948 }
1949 else
1950 {
1951 range_begin = main_range + 1;
1952 range_end = ranges.end();
1953 }
1954
1955 std::vector<Range>::const_iterator p =
1956 Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1957
1958 // if n is not in this set
1959 if (p == range_end || p->end == n || p->begin > n)
1961
1962 Assert(p != ranges.end(), ExcInternalError());
1963 Assert(p->begin <= n, ExcInternalError());
1964 Assert(n < p->end, ExcInternalError());
1965 return (n - p->begin) + p->nth_index_in_set;
1966}
1967
1968
1969
1970inline bool
1972{
1973 Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1974
1975 compress();
1976 is.compress();
1977
1978 return ranges == is.ranges;
1979}
1980
1981
1982
1983inline bool
1985{
1986 Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
1987
1988 compress();
1989 is.compress();
1990
1991 return ranges != is.ranges;
1992}
1993
1994
1995
1996template <typename Vector>
1997void
1999{
2000 Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
2001
2002 compress();
2003 // first fill all elements of the vector with zeroes.
2004 std::fill(vector.begin(), vector.end(), 0);
2005
2006 // then write ones into the elements whose indices are contained in the
2007 // index set
2008 for (const auto &range : ranges)
2009 for (size_type i = range.begin; i < range.end; ++i)
2010 vector[i] = 1;
2011}
2012
2013
2014
2015template <class StreamType>
2016inline void
2017IndexSet::print(StreamType &out) const
2018{
2019 compress();
2020 out << '{';
2021 std::vector<Range>::const_iterator p;
2022 for (p = ranges.begin(); p != ranges.end(); ++p)
2023 {
2024 if (p->end - p->begin == 1)
2025 out << p->begin;
2026 else
2027 out << '[' << p->begin << ',' << p->end - 1 << ']';
2028
2029 if (p != --ranges.end())
2030 out << ", ";
2031 }
2032 out << '}' << std::endl;
2033}
2034
2035
2036
2037template <class Archive>
2038inline void
2039IndexSet::Range::serialize(Archive &ar, const unsigned int)
2040{
2042}
2043
2044
2045
2046template <class Archive>
2047inline void
2048IndexSet::serialize(Archive &ar, const unsigned int)
2049{
2051}
2052
2054
2055#endif
size_type operator*() const
Definition: index_set.h:1312
bool operator==(const ElementIterator &) const
Definition: index_set.h:1324
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition: index_set.h:1268
std::ptrdiff_t difference_type
Definition: index_set.h:794
ElementIterator & operator++()
Definition: index_set.h:1365
bool operator<(const ElementIterator &) const
Definition: index_set.h:1393
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1384
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1406
std::forward_iterator_tag iterator_category
Definition: index_set.h:792
const IndexSet * index_set
Definition: index_set.h:808
size_type last() const
Definition: index_set.h:1098
ElementIterator end() const
Definition: index_set.h:1084
size_type n_elements() const
Definition: index_set.h:1058
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1120
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:1108
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1132
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:1027
const IndexSet * index_set
Definition: index_set.h:599
ElementIterator begin() const
Definition: index_set.h:1075
std::ptrdiff_t difference_type
Definition: index_set.h:702
const IntervalAccessor & operator*() const
Definition: index_set.h:1200
IntervalIterator(const IntervalIterator &other)=default
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1216
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1243
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1234
std::forward_iterator_tag iterator_category
Definition: index_set.h:700
IntervalIterator & operator=(const IntervalIterator &other)=default
IntervalIterator & operator++()
Definition: index_set.h:1181
const IntervalAccessor * operator->() const
Definition: index_set.h:1208
IntervalAccessor accessor
Definition: index_set.h:710
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1225
size_type largest_range
Definition: index_set.h:980
bool is_contiguous() const
Definition: index_set.h:1817
unsigned int n_intervals() const
Definition: index_set.h:1859
IndexSet(const IndexSet &)=default
IntervalIterator end_intervals() const
Definition: index_set.h:1603
void do_compress() const
Definition: index_set.cc:98
ElementIterator at(const size_type global_index) const
Definition: index_set.h:1532
size_type size() const
Definition: index_set.h:1636
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:512
bool is_empty() const
Definition: index_set.h:1826
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1923
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
Definition: index_set.cc:239
size_type n_elements() const
Definition: index_set.h:1834
bool operator==(const IndexSet &is) const
Definition: index_set.h:1971
bool is_element(const size_type index) const
Definition: index_set.h:1767
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2048
ElementIterator begin() const
Definition: index_set.h:1520
size_type pop_front()
Definition: index_set.cc:365
signed int value_type
Definition: index_set.h:99
void set_size(const size_type size)
Definition: index_set.h:1624
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1984
void read(std::istream &in)
Definition: index_set.cc:456
void clear()
Definition: index_set.h:1612
size_type largest_range_starting_index() const
Definition: index_set.h:1868
Tpetra::Map< int, types::global_dof_index > make_tpetra_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:532
IndexSet tensor_product(const IndexSet &other) const
Definition: index_set.cc:335
void add_index(const size_type index)
Definition: index_set.h:1655
void write(std::ostream &out) const
Definition: index_set.cc:441
void block_read(std::istream &in)
Definition: index_set.cc:492
bool is_compressed
Definition: index_set.h:963
IntervalIterator begin_intervals() const
Definition: index_set.h:1591
IndexSet & operator=(const IndexSet &)=default
void fill_binary_vector(VectorType &vector) const
std::vector< Range > ranges
Definition: index_set.h:953
void subtract_set(const IndexSet &other)
Definition: index_set.cc:266
ElementIterator end() const
Definition: index_set.h:1582
Threads::Mutex compress_mutex
Definition: index_set.h:986
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1014
size_type index_space_size
Definition: index_set.h:969
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:601
void block_write(std::ostream &out) const
Definition: index_set.cc:478
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:212
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1675
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1882
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:671
std::size_t memory_consumption() const
Definition: index_set.cc:737
void print(StreamType &out) const
Definition: index_set.h:2017
void compress() const
Definition: index_set.h:1644
size_type pop_back()
Definition: index_set.cc:347
types::global_dof_index size_type
Definition: index_set.h:80
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1705
IndexSet operator&(const IndexSet &is) const
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcIndexNotPresent(size_type arg1)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcMessage(std::string arg1)
iterator end()
size_type size() const
iterator begin()
static const char N
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
std::string compress(const std::string &input)
Definition: utilities.cc:392
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1153
static const unsigned int invalid_unsigned_int
Definition: types.h:201
const types::global_dof_index invalid_dof_index
Definition: types.h:216
unsigned int global_dof_index
Definition: types.h:76
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:911
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition: index_set.h:917
friend bool operator==(const Range &range_1, const Range &range_2)
Definition: index_set.h:924
size_type end
Definition: index_set.h:879
size_type nth_index_in_set
Definition: index_set.h:881
static std::size_t memory_consumption()
Definition: index_set.h:930
size_type begin
Definition: index_set.h:878
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:2039
friend bool operator<(const Range &range_1, const Range &range_2)
Definition: index_set.h:903
void advance(std::tuple< I1, I2 > &t, const unsigned int n)