Reference documentation for deal.II version GIT relicensing-1062-gc06da148b8 2024-07-15 19:20:02+00:00
\(\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\}}\)
Loading...
Searching...
No Matches
index_set.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2009 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_index_set_h
16#define dealii_index_set_h
17
18#include <deal.II/base/config.h>
19
22#include <deal.II/base/mutex.h>
24
25#include <boost/container/small_vector.hpp>
26
27#include <algorithm>
28#include <vector>
29
30
31#ifdef DEAL_II_WITH_TRILINOS
32# include <Epetra_Map.h>
33# ifdef DEAL_II_TRILINOS_WITH_TPETRA
34# include <Tpetra_Map.hpp>
35# endif
36#endif
37
38#ifdef DEAL_II_WITH_PETSC
39# include <petscis.h>
40#endif
41
43
67{
68public:
69 // forward declarations:
70 class ElementIterator;
71 class IntervalIterator;
72
78
96 using value_type = signed int;
97
98
102 IndexSet();
103
107 explicit IndexSet(const size_type size);
108
112 IndexSet(const IndexSet &) = default;
113
117 IndexSet &
118 operator=(const IndexSet &) = default;
119
124 IndexSet(IndexSet &&is) noexcept;
125
130 IndexSet &
131 operator=(IndexSet &&is) noexcept;
132
133#ifdef DEAL_II_WITH_TRILINOS
134
135# ifdef DEAL_II_TRILINOS_WITH_TPETRA
139 explicit IndexSet(
140 const Teuchos::RCP<const Tpetra::Map<int, types::signed_global_dof_index>>
141 &map);
142# endif // DEAL_II_TRILINOS_WITH_TPETRA
143
147 explicit IndexSet(const Epetra_BlockMap &map);
148#endif // DEAL_II_WITH_TRILINOS
149
154 void
155 clear();
156
163 void
164 set_size(const size_type size);
165
174 size() const;
175
182 void
183 add_range(const size_type begin, const size_type end);
184
193 void
194 add_index(const size_type index);
195
220 template <typename ForwardIterator>
221 void
222 add_indices(const ForwardIterator &begin, const ForwardIterator &end);
223
239 void
240 add_indices(const IndexSet &other, const size_type offset = 0);
241
245 bool
246 is_element(const size_type index) const;
247
252 bool
253 is_contiguous() const;
254
264 bool
265 is_empty() const;
266
276 bool
277 is_ascending_and_one_to_one(const MPI_Comm communicator) const;
278
283 n_elements() const;
284
291 nth_index_in_set(const size_type local_index) const;
292
300 index_within_set(const size_type global_index) const;
301
310 unsigned int
311 n_intervals() const;
312
326
331 void
332 compress() const;
333
348 bool
349 operator==(const IndexSet &is) const;
350
365 bool
366 operator!=(const IndexSet &is) const;
367
375 operator&(const IndexSet &is) const;
376
394 get_view(const size_type begin, const size_type end) const;
395
442 get_view(const IndexSet &mask) const;
443
449 std::vector<IndexSet>
451 const std::vector<types::global_dof_index> &n_indices_per_block) const;
452
458 void
459 subtract_set(const IndexSet &other);
460
478 tensor_product(const IndexSet &other) const;
479
485 pop_back();
486
492 pop_front();
493
502 std::vector<size_type>
503 get_index_vector() const;
504
512 void
513 fill_index_vector(std::vector<size_type> &indices) const;
514
527 template <typename VectorType>
528 void
529 fill_binary_vector(VectorType &vector) const;
530
538 bool
539 is_subset_of(const IndexSet &other) const;
540
545 template <typename StreamType>
546 void
547 print(StreamType &out) const;
548
553 void
554 write(std::ostream &out) const;
555
560 void
561 read(std::istream &in);
562
567 void
568 block_write(std::ostream &out) const;
569
574 void
575 block_read(std::istream &in);
576
577#ifdef DEAL_II_WITH_TRILINOS
606 Epetra_Map
607 make_trilinos_map(const MPI_Comm communicator = MPI_COMM_WORLD,
608 const bool overlapping = false) const;
609
610# ifdef DEAL_II_TRILINOS_WITH_TPETRA
611 Tpetra::Map<int, types::signed_global_dof_index>
612 make_tpetra_map(const MPI_Comm communicator = MPI_COMM_WORLD,
613 const bool overlapping = false) const;
614
615 Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index>>
616 make_tpetra_map_rcp(const MPI_Comm communicator = MPI_COMM_WORLD,
617 const bool overlapping = false) const;
618# endif
619#endif
620
621#ifdef DEAL_II_WITH_PETSC
622 IS
623 make_petsc_is(const MPI_Comm communicator = MPI_COMM_WORLD) const;
624#endif
625
626
631 std::size_t
632 memory_consumption() const;
633
635 size_type,
636 << "The global index " << arg1
637 << " is not an element of this set.");
638
644 template <class Archive>
645 void
646 serialize(Archive &ar, const unsigned int version);
647
648
660 {
661 public:
666 IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
667
671 explicit IntervalAccessor(const IndexSet *idxset);
672
677 n_elements() const;
678
682 bool
683 is_valid() const;
684
689 begin() const;
690
696 end() const;
697
702 last() const;
703
704 private:
713 operator=(const IntervalAccessor &other);
714
718 bool
719 operator==(const IntervalAccessor &other) const;
723 bool
724 operator<(const IntervalAccessor &other) const;
729 void
730 advance();
735
741
742 friend class IntervalIterator;
743 };
744
750 {
751 public:
756 IntervalIterator(const IndexSet *idxset, const size_type range_idx);
757
761 explicit IntervalIterator(const IndexSet *idxset);
762
767
771 IntervalIterator(const IntervalIterator &other) = default;
772
777 operator=(const IntervalIterator &other) = default;
778
783 operator++();
784
789 operator++(int);
790
794 const IntervalAccessor &
795 operator*() const;
796
800 const IntervalAccessor *
801 operator->() const;
802
806 bool
807 operator==(const IntervalIterator &) const;
808
812 bool
813 operator!=(const IntervalIterator &) const;
814
818 bool
819 operator<(const IntervalIterator &) const;
820
827 int
828 operator-(const IntervalIterator &p) const;
829
835 using iterator_category = std::forward_iterator_tag;
837 using difference_type = std::ptrdiff_t;
840
841 private:
846 };
847
853 {
854 public:
859 ElementIterator(const IndexSet *idxset,
860 const size_type range_idx,
861 const size_type index);
862
866 explicit ElementIterator(const IndexSet *idxset);
867
873 operator*() const;
874
878 bool
879 is_valid() const;
880
885 operator++();
886
891 operator++(int);
892
896 bool
897 operator==(const ElementIterator &) const;
898
902 bool
903 operator!=(const ElementIterator &) const;
904
908 bool
909 operator<(const ElementIterator &) const;
910
919 std::ptrdiff_t
920 operator-(const ElementIterator &p) const;
921
927 using iterator_category = std::forward_iterator_tag;
929 using difference_type = std::ptrdiff_t;
932
933 private:
937 void
938 advance();
939
952 };
953
959 begin() const;
960
976 at(const size_type global_index) const;
977
983 end() const;
984
989 begin_intervals() const;
990
996 end_intervals() const;
997
1002private:
1011 struct Range
1012 {
1015
1017
1026 Range();
1027
1035 Range(const size_type i1, const size_type i2);
1036
1037 friend inline bool
1038 operator<(const Range &range_1, const Range &range_2)
1039 {
1040 return (
1041 (range_1.begin < range_2.begin) ||
1042 ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
1043 }
1044
1045 static bool
1047 {
1048 return x.end < y.end;
1049 }
1050
1051 static bool
1053 {
1054 return (x.nth_index_in_set + (x.end - x.begin) <
1055 y.nth_index_in_set + (y.end - y.begin));
1056 }
1057
1058 friend inline bool
1059 operator==(const Range &range_1, const Range &range_2)
1060 {
1061 return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
1062 }
1063
1064 static std::size_t
1066 {
1067 return sizeof(Range);
1068 }
1069
1075 template <class Archive>
1076 void
1077 serialize(Archive &ar, const unsigned int version);
1078 };
1079
1088 mutable std::vector<Range> ranges;
1089
1098 mutable bool is_compressed;
1099
1105
1116
1122
1126 void
1127 do_compress() const;
1128
1135 bool
1136 is_element_binary_search(const size_type local_index) const;
1137
1144 size_type
1145 nth_index_in_set_binary_search(const size_type local_index) const;
1146
1153 size_type
1154 index_within_set_binary_search(const size_type global_index) const;
1155
1161 void
1162 add_range_lower_bound(const Range &range);
1163
1167 void
1169 boost::container::small_vector<std::pair<size_type, size_type>, 200>
1170 &tmp_ranges,
1171 const bool ranges_are_sorted);
1172};
1173
1174
1175
1193inline IndexSet
1195{
1196 IndexSet is(N);
1197 is.add_range(0, N);
1198 is.compress();
1199 return is;
1200}
1201
1202/* ------------------ inline functions ------------------ */
1203
1204
1205/* IntervalAccessor */
1206
1208 const IndexSet *idxset,
1209 const IndexSet::size_type range_idx)
1210 : index_set(idxset)
1211 , range_idx(range_idx)
1212{
1213 Assert(range_idx < idxset->n_intervals(),
1214 ExcInternalError("Invalid range index"));
1215}
1216
1217
1218
1220 : index_set(idxset)
1221 , range_idx(numbers::invalid_dof_index)
1222{}
1223
1224
1225
1227 const IndexSet::IntervalAccessor &other)
1228 : index_set(other.index_set)
1229 , range_idx(other.range_idx)
1230{
1232 ExcMessage("invalid iterator"));
1233}
1234
1235
1236
1239{
1240 Assert(is_valid(), ExcMessage("invalid iterator"));
1241 return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1242}
1243
1244
1245
1246inline bool
1248{
1249 return index_set != nullptr && range_idx < index_set->n_intervals();
1250}
1251
1252
1253
1256{
1257 Assert(is_valid(), ExcMessage("invalid iterator"));
1258 return {index_set, range_idx, index_set->ranges[range_idx].begin};
1259}
1260
1261
1262
1265{
1266 Assert(is_valid(), ExcMessage("invalid iterator"));
1267
1268 // point to first index in next interval unless we are the last interval.
1269 if (range_idx < index_set->ranges.size() - 1)
1270 return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1271 else
1272 return index_set->end();
1273}
1274
1275
1276
1279{
1280 Assert(is_valid(), ExcMessage("invalid iterator"));
1281
1282 return index_set->ranges[range_idx].end - 1;
1283}
1284
1285
1286
1289{
1290 index_set = other.index_set;
1291 range_idx = other.range_idx;
1292 Assert(range_idx == numbers::invalid_dof_index || is_valid(),
1293 ExcMessage("invalid iterator"));
1294 return *this;
1295}
1296
1297
1298
1299inline bool
1301 const IndexSet::IntervalAccessor &other) const
1302{
1303 Assert(index_set == other.index_set,
1304 ExcMessage(
1305 "Can not compare accessors pointing to different IndexSets"));
1306 return range_idx == other.range_idx;
1307}
1308
1309
1310
1311inline bool
1313 const IndexSet::IntervalAccessor &other) const
1314{
1315 Assert(index_set == other.index_set,
1316 ExcMessage(
1317 "Can not compare accessors pointing to different IndexSets"));
1318 return range_idx < other.range_idx;
1319}
1320
1321
1322
1323inline void
1325{
1326 Assert(
1327 is_valid(),
1328 ExcMessage(
1329 "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1330 ++range_idx;
1331
1332 // set ourselves to invalid if we walk off the end
1333 if (range_idx >= index_set->ranges.size())
1334 range_idx = numbers::invalid_dof_index;
1335}
1336
1337
1338/* IntervalIterator */
1339
1341 const IndexSet *idxset,
1342 const IndexSet::size_type range_idx)
1343 : accessor(idxset, range_idx)
1344{}
1345
1346
1347
1349 : accessor(nullptr)
1350{}
1351
1352
1353
1355 : accessor(idxset)
1356{}
1357
1358
1359
1362{
1363 accessor.advance();
1364 return *this;
1365}
1366
1367
1368
1371{
1372 const IndexSet::IntervalIterator iter = *this;
1373 accessor.advance();
1374 return iter;
1375}
1376
1377
1378
1379inline const IndexSet::IntervalAccessor &
1381{
1382 return accessor;
1383}
1384
1385
1386
1387inline const IndexSet::IntervalAccessor *
1389{
1390 return &accessor;
1391}
1392
1393
1394
1395inline bool
1397 const IndexSet::IntervalIterator &other) const
1398{
1399 return accessor == other.accessor;
1400}
1401
1402
1403
1404inline bool
1406 const IndexSet::IntervalIterator &other) const
1407{
1408 return !(*this == other);
1409}
1410
1411
1412
1413inline bool
1415 const IndexSet::IntervalIterator &other) const
1416{
1417 return accessor < other.accessor;
1418}
1419
1420
1421
1422inline int
1424 const IndexSet::IntervalIterator &other) const
1425{
1426 Assert(accessor.index_set == other.accessor.index_set,
1427 ExcMessage(
1428 "Can not compare iterators belonging to different IndexSets"));
1429
1430 const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index) ?
1431 accessor.index_set->ranges.size() :
1432 accessor.range_idx;
1433 const size_type rhs =
1435 accessor.index_set->ranges.size() :
1436 other.accessor.range_idx;
1437
1438 if (lhs > rhs)
1439 return static_cast<int>(lhs - rhs);
1440 else
1441 return -static_cast<int>(rhs - lhs);
1442}
1443
1444
1445
1446/* ElementIterator */
1447
1449 const IndexSet *idxset,
1450 const IndexSet::size_type range_idx,
1451 const IndexSet::size_type index)
1452 : index_set(idxset)
1453 , range_idx(range_idx)
1454 , idx(index)
1455{
1456 Assert(range_idx < index_set->ranges.size(),
1457 ExcMessage(
1458 "Invalid range index for IndexSet::ElementIterator constructor."));
1459 Assert(
1460 idx >= index_set->ranges[range_idx].begin &&
1461 idx < index_set->ranges[range_idx].end,
1463 "Invalid index argument for IndexSet::ElementIterator constructor."));
1464}
1465
1466
1467
1469 : index_set(idxset)
1470 , range_idx(numbers::invalid_dof_index)
1471 , idx(numbers::invalid_dof_index)
1472{}
1473
1474
1475
1476inline bool
1478{
1479 Assert((range_idx == numbers::invalid_dof_index &&
1481 (range_idx < index_set->ranges.size() &&
1482 idx < index_set->ranges[range_idx].end),
1483 ExcInternalError("Invalid ElementIterator state."));
1484
1485 return (range_idx < index_set->ranges.size() &&
1486 idx < index_set->ranges[range_idx].end);
1487}
1488
1489
1490
1493{
1494 Assert(
1495 is_valid(),
1496 ExcMessage(
1497 "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1498 return idx;
1499}
1500
1501
1502
1503inline bool
1505 const IndexSet::ElementIterator &other) const
1506{
1507 Assert(index_set == other.index_set,
1508 ExcMessage(
1509 "Can not compare iterators belonging to different IndexSets"));
1510 return range_idx == other.range_idx && idx == other.idx;
1511}
1512
1513
1514
1515inline void
1517{
1518 Assert(
1519 is_valid(),
1520 ExcMessage(
1521 "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1522 if (idx < index_set->ranges[range_idx].end)
1523 ++idx;
1524 // end of this range?
1525 if (idx == index_set->ranges[range_idx].end)
1526 {
1527 // point to first element in next interval if possible
1528 if (range_idx < index_set->ranges.size() - 1)
1529 {
1530 ++range_idx;
1531 idx = index_set->ranges[range_idx].begin;
1532 }
1533 else
1534 {
1535 // we just fell off the end, set to invalid:
1536 range_idx = numbers::invalid_dof_index;
1538 }
1539 }
1540}
1541
1542
1543
1546{
1547 advance();
1548 return *this;
1549}
1550
1551
1552
1555{
1556 const IndexSet::ElementIterator it = *this;
1557 advance();
1558 return it;
1559}
1560
1561
1562
1563inline bool
1565 const IndexSet::ElementIterator &other) const
1566{
1567 return !(*this == other);
1568}
1569
1570
1571
1572inline bool
1574 const IndexSet::ElementIterator &other) const
1575{
1576 Assert(index_set == other.index_set,
1577 ExcMessage(
1578 "Can not compare iterators belonging to different IndexSets"));
1579 return range_idx < other.range_idx ||
1580 (range_idx == other.range_idx && idx < other.idx);
1581}
1582
1583
1584
1585inline std::ptrdiff_t
1587 const IndexSet::ElementIterator &other) const
1588{
1589 Assert(index_set == other.index_set,
1590 ExcMessage(
1591 "Can not compare iterators belonging to different IndexSets"));
1592 if (*this == other)
1593 return 0;
1594 if (!(*this < other))
1595 return -(other - *this);
1596
1597 // only other can be equal to end() because of the checks above.
1598 Assert(is_valid(), ExcInternalError());
1599
1600 // Note: we now compute how far advance *this in "*this < other" to get other,
1601 // so we need to return -c at the end.
1602
1603 // first finish the current range:
1604 std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1605
1606 // now walk in steps of ranges (need to start one behind our current one):
1607 for (size_type range = range_idx + 1;
1608 range < index_set->ranges.size() && range <= other.range_idx;
1609 ++range)
1610 c += index_set->ranges[range].end - index_set->ranges[range].begin;
1611
1612 Assert(
1613 other.range_idx < index_set->ranges.size() ||
1615 ExcMessage(
1616 "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1617
1618 // We might have walked too far because we went until the end of
1619 // other.range_idx, so walk backwards to other.idx:
1621 c -= index_set->ranges[other.range_idx].end - other.idx;
1622
1623 return -c;
1624}
1625
1626
1627/* Range */
1628
1630 : begin(numbers::invalid_dof_index)
1631 , end(numbers::invalid_dof_index)
1632 , nth_index_in_set(numbers::invalid_dof_index)
1633{}
1634
1635
1636
1638 : begin(i1)
1639 , end(i2)
1640 , nth_index_in_set(numbers::invalid_dof_index)
1641{}
1642
1643
1644
1645/* IndexSet itself */
1646
1648 : is_compressed(true)
1649 , index_space_size(0)
1650 , largest_range(numbers::invalid_unsigned_int)
1651{}
1652
1653
1654
1656 : is_compressed(true)
1657 , index_space_size(size)
1658 , largest_range(numbers::invalid_unsigned_int)
1659{}
1660
1661
1662
1663inline IndexSet::IndexSet(IndexSet &&is) noexcept
1664 : ranges(std::move(is.ranges))
1665 , is_compressed(is.is_compressed)
1666 , index_space_size(is.index_space_size)
1667 , largest_range(is.largest_range)
1668{
1669 is.ranges.clear();
1670 is.is_compressed = true;
1671 is.index_space_size = 0;
1672 is.largest_range = numbers::invalid_unsigned_int;
1673
1674 compress();
1675}
1676
1677
1678
1679inline IndexSet &
1681{
1682 ranges = std::move(is.ranges);
1683 is_compressed = is.is_compressed;
1684 index_space_size = is.index_space_size;
1685 largest_range = is.largest_range;
1686
1687 is.ranges.clear();
1688 is.is_compressed = true;
1689 is.index_space_size = 0;
1690 is.largest_range = numbers::invalid_unsigned_int;
1691
1692 compress();
1693
1694 return *this;
1695}
1696
1697
1698
1701{
1702 compress();
1703 if (ranges.size() > 0)
1704 return {this, 0, ranges[0].begin};
1705 else
1706 return end();
1707}
1708
1709
1710
1713{
1714 compress();
1715 return IndexSet::ElementIterator(this);
1716}
1717
1718
1719
1722{
1723 compress();
1724 if (ranges.size() > 0)
1725 return IndexSet::IntervalIterator(this, 0);
1726 else
1727 return end_intervals();
1728}
1729
1730
1731
1734{
1735 compress();
1736 return IndexSet::IntervalIterator(this);
1737}
1738
1739
1740
1741inline void
1743{
1744 // reset so that there are no indices in the set any more; however,
1745 // as documented, the index set retains its size
1746 ranges.clear();
1747 is_compressed = true;
1749}
1750
1751
1752
1753inline void
1755{
1756 Assert(ranges.empty(),
1757 ExcMessage("This function can only be called if the current "
1758 "object does not yet contain any elements."));
1759 index_space_size = sz;
1760 is_compressed = true;
1761}
1762
1763
1764
1767{
1768 return index_space_size;
1769}
1770
1771
1772
1773inline void
1775{
1776 if (is_compressed == true)
1777 return;
1778
1779 do_compress();
1780}
1781
1782
1783
1784inline void
1786{
1787 add_range(index, index + 1);
1788}
1789
1790
1791
1792inline void
1794{
1797 ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1799 ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1801
1802 if (begin != end)
1803 {
1804 // the new index might be larger than the last index present in the
1805 // ranges. Then we can skip the binary search
1806 if (ranges.empty() || begin > ranges.back().end)
1807 ranges.emplace_back(begin, end);
1808 else if (begin == ranges.back().end)
1809 ranges.back().end = end;
1810 else
1812
1813 is_compressed = false;
1814 }
1815}
1816
1817
1818
1819template <typename ForwardIterator>
1820inline void
1821IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1822{
1823 if (begin == end)
1824 return;
1825
1826 // identify ranges in the given iterator range by checking whether some
1827 // indices happen to be consecutive. to avoid quadratic complexity when
1828 // calling add_range many times (as add_range() going into the middle of an
1829 // already existing range must shift entries around), we first collect a
1830 // vector of ranges.
1831 boost::container::small_vector<std::pair<size_type, size_type>, 200>
1832 tmp_ranges;
1833 bool ranges_are_sorted = true;
1834 for (ForwardIterator p = begin; p != end;)
1835 {
1836 // Starting with the current iterator 'p', find an iterator
1837 // 'q' so that the indices pointed to by the iterators in
1838 // the range [p,q) are consecutive. These indices then form
1839 // a range that is contiguous, and that can be added all
1840 // at once.
1841 const size_type begin_index = *p;
1842 size_type end_index = begin_index + 1;
1843
1844 // Start looking at the position after 'p', and keep iterating while
1845 // 'q' points to a duplicate of 'p':
1846 ForwardIterator q = p;
1847 ++q;
1848 while ((q != end) && (*q == *p))
1849 ++q;
1850
1851 // Now we know that 'q' is either past the end, or points to a value
1852 // other than 'p'. If it points to 'end_index', we are still good with
1853 // a contiguous range; then increment the end index of that range, and
1854 // move to the next iterator that is not a duplicate of what
1855 // we were just looking at:
1856 while ((q != end) && (static_cast<size_type>(*q) == end_index))
1857 {
1858 ++q;
1859 while ((q != end) && (static_cast<size_type>(*q) == end_index))
1860 ++q;
1861
1862 ++end_index;
1863 }
1864
1865 // Add this range:
1866 tmp_ranges.emplace_back(begin_index, end_index);
1867
1868 // Then move on to the next element in the input range.
1869 // If the starting index of the next go-around of the for loop is less
1870 // than the end index of the one just identified, then we will have at
1871 // least one pair of ranges that are not sorted, and consequently the
1872 // whole collection of ranges is not sorted.
1873 p = q;
1874 if ((p != end) && (static_cast<size_type>(*p) < end_index))
1875 ranges_are_sorted = false;
1876 }
1877
1878 add_ranges_internal(tmp_ranges, ranges_are_sorted);
1879}
1880
1881
1882
1883inline bool
1885{
1886 if (ranges.empty() == false)
1887 {
1888 compress();
1889
1890 // fast check whether the index is in the largest range
1892 if (index >= ranges[largest_range].begin &&
1893 index < ranges[largest_range].end)
1894 return true;
1895 else if (ranges.size() > 1)
1896 return is_element_binary_search(index);
1897 else
1898 return false;
1899 }
1900 else
1901 return false;
1902}
1903
1904
1905
1906inline bool
1908{
1909 compress();
1910 return (ranges.size() <= 1);
1911}
1912
1913
1914
1915inline bool
1917{
1918 return ranges.empty();
1919}
1920
1921
1922
1925{
1926 // make sure we have non-overlapping ranges
1927 compress();
1928
1929 size_type v = 0;
1930 if (!ranges.empty())
1931 {
1932 const Range &r = ranges.back();
1933 v = r.nth_index_in_set + r.end - r.begin;
1934 }
1935
1936#ifdef DEBUG
1937 size_type s = 0;
1938 for (const auto &range : ranges)
1939 s += (range.end - range.begin);
1940 Assert(s == v, ExcInternalError());
1941#endif
1942
1943 return v;
1944}
1945
1946
1947
1948inline unsigned int
1950{
1951 compress();
1952 return ranges.size();
1953}
1954
1955
1956
1959{
1960 Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1961
1962 compress();
1963 const std::vector<Range>::const_iterator main_range =
1964 ranges.begin() + largest_range;
1965
1966 return main_range->nth_index_in_set;
1967}
1968
1969
1970
1973{
1975
1976 compress();
1977
1978 // first check whether the index is in the largest range
1980 const auto main_range = ranges.begin() + largest_range;
1981 if (n >= main_range->nth_index_in_set &&
1982 n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1983 return main_range->begin + (n - main_range->nth_index_in_set);
1984 else
1986}
1987
1988
1989
1992{
1993 // to make this call thread-safe, compress() must not be called through this
1994 // function
1995 Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1996 AssertIndexRange(n, size());
1997
1998 // return immediately if the index set is empty
1999 if (is_empty())
2001
2002 // check whether the index is in the largest range. use the result to
2003 // perform a one-sided binary search afterward
2006 return (n - ranges[largest_range].begin) +
2007 ranges[largest_range].nth_index_in_set;
2008 else if (ranges.size() > 1)
2010 else
2012}
2013
2014
2015
2016inline bool
2018{
2019 // If one of the two index sets has size zero, the other one has to
2020 // have size zero as well:
2021 if (size() == 0)
2022 return (is.size() == 0);
2023 if (is.size() == 0)
2024 return (size() == 0);
2025
2026 // Otherwise, they must have the same size (see the documentation):
2027 Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
2028
2029 compress();
2030 is.compress();
2031
2032 return (ranges == is.ranges);
2033}
2034
2035
2036
2037inline bool
2039{
2040 // If one of the two index sets has size zero, the other one has to
2041 // have a non-zero size for inequality:
2042 if (size() == 0)
2043 return (is.size() != 0);
2044 if (is.size() == 0)
2045 return (size() != 0);
2046
2047 // Otherwise, they must have the same size (see the documentation):
2048 Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
2049
2050 compress();
2051 is.compress();
2052
2053 return (ranges != is.ranges);
2054}
2055
2056
2057
2058template <typename Vector>
2059void
2061{
2062 Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
2063
2064 compress();
2065 // first fill all elements of the vector with zeroes.
2066 std::fill(vector.begin(), vector.end(), 0);
2067
2068 // then write ones into the elements whose indices are contained in the
2069 // index set
2070 for (const auto &range : ranges)
2071 for (size_type i = range.begin; i < range.end; ++i)
2072 vector[i] = 1;
2073}
2074
2075
2076
2077template <typename StreamType>
2078inline void
2079IndexSet::print(StreamType &out) const
2080{
2081 compress();
2082 out << '{';
2083 std::vector<Range>::const_iterator p;
2084 for (p = ranges.begin(); p != ranges.end(); ++p)
2085 {
2086 if (p->end - p->begin == 1)
2087 out << p->begin;
2088 else
2089 out << '[' << p->begin << ',' << p->end - 1 << ']';
2090
2091 if (p != --ranges.end())
2092 out << ", ";
2093 }
2094 out << '}' << std::endl;
2095}
2096
2097
2098
2099template <class Archive>
2100inline void
2101IndexSet::Range::serialize(Archive &ar, const unsigned int)
2102{
2104}
2105
2106
2107
2108template <class Archive>
2109inline void
2110IndexSet::serialize(Archive &ar, const unsigned int)
2111{
2113}
2114
2116
2117#endif
size_type operator*() const
Definition index_set.h:1492
bool operator==(const ElementIterator &) const
Definition index_set.h:1504
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition index_set.h:1448
std::ptrdiff_t difference_type
Definition index_set.h:929
ElementIterator & operator++()
Definition index_set.h:1545
bool operator<(const ElementIterator &) const
Definition index_set.h:1573
bool operator!=(const ElementIterator &) const
Definition index_set.h:1564
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition index_set.h:1586
std::forward_iterator_tag iterator_category
Definition index_set.h:927
const IndexSet * index_set
Definition index_set.h:943
size_type last() const
Definition index_set.h:1278
ElementIterator end() const
Definition index_set.h:1264
size_type n_elements() const
Definition index_set.h:1238
bool operator==(const IntervalAccessor &other) const
Definition index_set.h:1300
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition index_set.h:1288
bool operator<(const IntervalAccessor &other) const
Definition index_set.h:1312
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition index_set.h:1207
const IndexSet * index_set
Definition index_set.h:734
ElementIterator begin() const
Definition index_set.h:1255
std::ptrdiff_t difference_type
Definition index_set.h:837
const IntervalAccessor & operator*() const
Definition index_set.h:1380
IntervalIterator(const IntervalIterator &other)=default
bool operator==(const IntervalIterator &) const
Definition index_set.h:1396
int operator-(const IntervalIterator &p) const
Definition index_set.h:1423
bool operator<(const IntervalIterator &) const
Definition index_set.h:1414
std::forward_iterator_tag iterator_category
Definition index_set.h:835
IntervalIterator & operator=(const IntervalIterator &other)=default
IntervalIterator & operator++()
Definition index_set.h:1361
const IntervalAccessor * operator->() const
Definition index_set.h:1388
IntervalAccessor accessor
Definition index_set.h:845
bool operator!=(const IntervalIterator &) const
Definition index_set.h:1405
bool is_subset_of(const IndexSet &other) const
Definition index_set.cc:704
bool is_element_binary_search(const size_type local_index) const
Definition index_set.cc:799
size_type index_within_set_binary_search(const size_type global_index) const
Definition index_set.cc:852
size_type largest_range
Definition index_set.h:1115
IS make_petsc_is(const MPI_Comm communicator=MPI_COMM_WORLD) const
bool is_ascending_and_one_to_one(const MPI_Comm communicator) const
bool is_contiguous() const
Definition index_set.h:1907
unsigned int n_intervals() const
Definition index_set.h:1949
IndexSet(const IndexSet &)=default
IntervalIterator end_intervals() const
Definition index_set.h:1733
void do_compress() const
Definition index_set.cc:137
ElementIterator at(const size_type global_index) const
Definition index_set.cc:874
size_type size() const
Definition index_set.h:1766
void fill_index_vector(std::vector< size_type > &indices) const
Definition index_set.cc:943
bool is_empty() const
Definition index_set.h:1916
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1991
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
Definition index_set.cc:442
size_type n_elements() const
Definition index_set.h:1924
bool operator==(const IndexSet &is) const
Definition index_set.h:2017
void add_range_lower_bound(const Range &range)
Definition index_set.cc:592
bool is_element(const size_type index) const
Definition index_set.h:1884
void serialize(Archive &ar, const unsigned int version)
Definition index_set.h:2110
ElementIterator begin() const
Definition index_set.h:1700
size_type pop_front()
Definition index_set.cc:570
signed int value_type
Definition index_set.h:96
void set_size(const size_type size)
Definition index_set.h:1754
bool operator!=(const IndexSet &is) const
Definition index_set.h:2038
void read(std::istream &in)
Definition index_set.cc:743
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
void clear()
Definition index_set.h:1742
size_type largest_range_starting_index() const
Definition index_set.h:1958
void add_index(const size_type index)
Definition index_set.h:1785
void write(std::ostream &out) const
Definition index_set.cc:728
void block_read(std::istream &in)
Definition index_set.cc:779
bool is_compressed
Definition index_set.h:1098
void add_ranges_internal(boost::container::small_vector< std::pair< size_type, size_type >, 200 > &tmp_ranges, const bool ranges_are_sorted)
Definition index_set.cc:608
IntervalIterator begin_intervals() const
Definition index_set.h:1721
Tpetra::Map< int, types::signed_global_dof_index > make_tpetra_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition index_set.cc:954
IndexSet & operator=(const IndexSet &)=default
void fill_binary_vector(VectorType &vector) const
std::vector< Range > ranges
Definition index_set.h:1088
void subtract_set(const IndexSet &other)
Definition index_set.cc:471
ElementIterator end() const
Definition index_set.h:1712
Threads::Mutex compress_mutex
Definition index_set.h:1121
IndexSet complete_index_set(const IndexSet::size_type N)
Definition index_set.h:1194
size_type index_space_size
Definition index_set.h:1104
void block_write(std::ostream &out) const
Definition index_set.cc:765
IndexSet get_view(const size_type begin, const size_type end) const
Definition index_set.cc:271
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1793
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1972
std::size_t memory_consumption() const
void print(StreamType &out) const
Definition index_set.h:2079
size_type nth_index_in_set_binary_search(const size_type local_index) const
Definition index_set.cc:835
void compress() const
Definition index_set.h:1774
Teuchos::RCP< Tpetra::Map< int, types::signed_global_dof_index > > make_tpetra_map_rcp(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition index_set.cc:963
size_type pop_back()
Definition index_set.cc:552
std::vector< size_type > get_index_vector() const
Definition index_set.cc:924
types::global_dof_index size_type
Definition index_set.h:77
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1821
IndexSet operator&(const IndexSet &is) const
virtual size_type size() const override
iterator end()
iterator begin()
#define DEAL_II_DEPRECATED
Definition config.h:207
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
static ::ExceptionBase & ExcIndexNotPresent(size_type arg1)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
static ::ExceptionBase & ExcMessage(std::string arg1)
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::global_dof_index invalid_dof_index
Definition types.h:252
unsigned int global_dof_index
Definition types.h:81
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition index_set.h:1046
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition index_set.h:1052
friend bool operator==(const Range &range_1, const Range &range_2)
Definition index_set.h:1059
size_type end
Definition index_set.h:1014
size_type nth_index_in_set
Definition index_set.h:1016
static std::size_t memory_consumption()
Definition index_set.h:1065
size_type begin
Definition index_set.h:1013
void serialize(Archive &ar, const unsigned int version)
Definition index_set.h:2101
friend bool operator<(const Range &range_1, const Range &range_2)
Definition index_set.h:1038
void advance(std::tuple< I1, I2 > &t, const unsigned int n)