Reference documentation for deal.II version 9.6.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\}}\)
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
23#include <deal.II/base/mutex.h>
25
27
28#include <boost/container/small_vector.hpp>
29
30#include <algorithm>
31#include <vector>
32
33
34#ifdef DEAL_II_WITH_TRILINOS
35# include <Epetra_Map.h>
36# ifdef DEAL_II_TRILINOS_WITH_TPETRA
37# include <Tpetra_Map.hpp>
38# endif
39#endif
40
41#ifdef DEAL_II_WITH_PETSC
42# include <petscis.h>
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
137
138# ifdef DEAL_II_TRILINOS_WITH_TPETRA
142 template <typename NodeType>
143 explicit IndexSet(
144 const Teuchos::RCP<
145 const Tpetra::Map<int, types::signed_global_dof_index, NodeType>> &map);
146# endif // DEAL_II_TRILINOS_WITH_TPETRA
147
151 explicit IndexSet(const Epetra_BlockMap &map);
152#endif // DEAL_II_WITH_TRILINOS
153
158 void
159 clear();
160
167 void
168 set_size(const size_type size);
169
178 size() const;
179
186 void
187 add_range(const size_type begin, const size_type end);
188
197 void
198 add_index(const size_type index);
199
224 template <typename ForwardIterator>
225 void
226 add_indices(const ForwardIterator &begin, const ForwardIterator &end);
227
243 void
244 add_indices(const IndexSet &other, const size_type offset = 0);
245
249 bool
250 is_element(const size_type index) const;
251
256 bool
257 is_contiguous() const;
258
268 bool
269 is_empty() const;
270
280 bool
281 is_ascending_and_one_to_one(const MPI_Comm communicator) const;
282
287 n_elements() const;
288
295 nth_index_in_set(const size_type local_index) const;
296
304 index_within_set(const size_type global_index) const;
305
314 unsigned int
315 n_intervals() const;
316
330
335 void
336 compress() const;
337
352 bool
353 operator==(const IndexSet &is) const;
354
369 bool
370 operator!=(const IndexSet &is) const;
371
379 operator&(const IndexSet &is) const;
380
398 get_view(const size_type begin, const size_type end) const;
399
446 get_view(const IndexSet &mask) const;
447
453 std::vector<IndexSet>
455 const std::vector<types::global_dof_index> &n_indices_per_block) const;
456
462 void
463 subtract_set(const IndexSet &other);
464
482 tensor_product(const IndexSet &other) const;
483
489 pop_back();
490
496 pop_front();
497
506 std::vector<size_type>
507 get_index_vector() const;
508
516 void
517 fill_index_vector(std::vector<size_type> &indices) const;
518
531 template <typename VectorType>
532 void
533 fill_binary_vector(VectorType &vector) const;
534
542 bool
543 is_subset_of(const IndexSet &other) const;
544
549 template <typename StreamType>
550 void
551 print(StreamType &out) const;
552
557 void
558 write(std::ostream &out) const;
559
564 void
565 read(std::istream &in);
566
571 void
572 block_write(std::ostream &out) const;
573
578 void
579 block_read(std::istream &in);
580
581#ifdef DEAL_II_WITH_TRILINOS
610 Epetra_Map
611 make_trilinos_map(const MPI_Comm communicator = MPI_COMM_WORLD,
612 const bool overlapping = false) const;
613
614# ifdef DEAL_II_TRILINOS_WITH_TPETRA
615 template <
616 typename NodeType =
618 Tpetra::Map<int, types::signed_global_dof_index, NodeType>
619 make_tpetra_map(const MPI_Comm communicator = MPI_COMM_WORLD,
620 const bool overlapping = false) const;
621
622 template <
623 typename NodeType =
625 Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index, NodeType>>
626 make_tpetra_map_rcp(const MPI_Comm communicator = MPI_COMM_WORLD,
627 const bool overlapping = false) const;
628# endif
629#endif
630
631#ifdef DEAL_II_WITH_PETSC
632 IS
633 make_petsc_is(const MPI_Comm communicator = MPI_COMM_WORLD) const;
634#endif
635
636
641 std::size_t
642 memory_consumption() const;
643
645 size_type,
646 << "The global index " << arg1
647 << " is not an element of this set.");
648
654 template <class Archive>
655 void
656 serialize(Archive &ar, const unsigned int version);
657
658
670 {
671 public:
676 IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
677
681 explicit IntervalAccessor(const IndexSet *idxset);
682
687 n_elements() const;
688
692 bool
693 is_valid() const;
694
699 begin() const;
700
706 end() const;
707
712 last() const;
713
714 private:
723 operator=(const IntervalAccessor &other);
724
728 bool
729 operator==(const IntervalAccessor &other) const;
733 bool
734 operator<(const IntervalAccessor &other) const;
739 void
740 advance();
745
751
752 friend class IntervalIterator;
753 };
754
760 {
761 public:
766 IntervalIterator(const IndexSet *idxset, const size_type range_idx);
767
771 explicit IntervalIterator(const IndexSet *idxset);
772
777
781 IntervalIterator(const IntervalIterator &other) = default;
782
787 operator=(const IntervalIterator &other) = default;
788
793 operator++();
794
799 operator++(int);
800
804 const IntervalAccessor &
805 operator*() const;
806
810 const IntervalAccessor *
811 operator->() const;
812
816 bool
817 operator==(const IntervalIterator &) const;
818
822 bool
823 operator!=(const IntervalIterator &) const;
824
828 bool
829 operator<(const IntervalIterator &) const;
830
837 int
838 operator-(const IntervalIterator &p) const;
839
845 using iterator_category = std::forward_iterator_tag;
847 using difference_type = std::ptrdiff_t;
850
851 private:
856 };
857
863 {
864 public:
869 ElementIterator(const IndexSet *idxset,
870 const size_type range_idx,
871 const size_type index);
872
876 explicit ElementIterator(const IndexSet *idxset);
877
883 operator*() const;
884
888 bool
889 is_valid() const;
890
895 operator++();
896
901 operator++(int);
902
906 bool
907 operator==(const ElementIterator &) const;
908
912 bool
913 operator!=(const ElementIterator &) const;
914
918 bool
919 operator<(const ElementIterator &) const;
920
929 std::ptrdiff_t
930 operator-(const ElementIterator &p) const;
931
937 using iterator_category = std::forward_iterator_tag;
939 using difference_type = std::ptrdiff_t;
942
943 private:
947 void
948 advance();
949
962 };
963
969 begin() const;
970
986 at(const size_type global_index) const;
987
993 end() const;
994
999 begin_intervals() const;
1000
1006 end_intervals() const;
1007
1012private:
1021 struct Range
1022 {
1025
1027
1036 Range();
1037
1045 Range(const size_type i1, const size_type i2);
1046
1047 friend inline bool
1048 operator<(const Range &range_1, const Range &range_2)
1049 {
1050 return (
1051 (range_1.begin < range_2.begin) ||
1052 ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
1053 }
1054
1055 static bool
1057 {
1058 return x.end < y.end;
1059 }
1060
1061 static bool
1063 {
1064 return (x.nth_index_in_set + (x.end - x.begin) <
1065 y.nth_index_in_set + (y.end - y.begin));
1066 }
1067
1068 friend inline bool
1069 operator==(const Range &range_1, const Range &range_2)
1070 {
1071 return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
1072 }
1073
1074 static std::size_t
1076 {
1077 return sizeof(Range);
1078 }
1079
1085 template <class Archive>
1086 void
1087 serialize(Archive &ar, const unsigned int version);
1088 };
1089
1098 mutable std::vector<Range> ranges;
1099
1108 mutable bool is_compressed;
1109
1115
1126
1132
1136 void
1137 do_compress() const;
1138
1145 bool
1146 is_element_binary_search(const size_type local_index) const;
1147
1154 size_type
1155 nth_index_in_set_binary_search(const size_type local_index) const;
1156
1163 size_type
1164 index_within_set_binary_search(const size_type global_index) const;
1165
1171 void
1172 add_range_lower_bound(const Range &range);
1173
1177 void
1179 boost::container::small_vector<std::pair<size_type, size_type>, 200>
1180 &tmp_ranges,
1181 const bool ranges_are_sorted);
1182};
1183
1184
1185
1203inline IndexSet
1205{
1206 IndexSet is(N);
1207 is.add_range(0, N);
1208 is.compress();
1209 return is;
1210}
1211
1212/* ------------------ inline functions ------------------ */
1213
1214
1215/* IntervalAccessor */
1216
1218 const IndexSet *idxset,
1219 const IndexSet::size_type range_idx)
1220 : index_set(idxset)
1221 , range_idx(range_idx)
1222{
1224 ExcInternalError("Invalid range index"));
1225}
1226
1227
1228
1230 : index_set(idxset)
1231 , range_idx(numbers::invalid_dof_index)
1232{}
1233
1234
1235
1237 const IndexSet::IntervalAccessor &other)
1238 : index_set(other.index_set)
1239 , range_idx(other.range_idx)
1240{
1242 ExcMessage("invalid iterator"));
1243}
1244
1245
1246
1249{
1250 Assert(is_valid(), ExcMessage("invalid iterator"));
1251 return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1252}
1253
1254
1255
1256inline bool
1258{
1259 return index_set != nullptr && range_idx < index_set->n_intervals();
1260}
1261
1262
1263
1266{
1267 Assert(is_valid(), ExcMessage("invalid iterator"));
1268 return {index_set, range_idx, index_set->ranges[range_idx].begin};
1269}
1270
1271
1272
1275{
1276 Assert(is_valid(), ExcMessage("invalid iterator"));
1277
1278 // point to first index in next interval unless we are the last interval.
1279 if (range_idx < index_set->ranges.size() - 1)
1280 return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1281 else
1282 return index_set->end();
1283}
1284
1285
1286
1289{
1290 Assert(is_valid(), ExcMessage("invalid iterator"));
1291
1292 return index_set->ranges[range_idx].end - 1;
1293}
1294
1295
1296
1299{
1300 index_set = other.index_set;
1301 range_idx = other.range_idx;
1302 Assert(range_idx == numbers::invalid_dof_index || is_valid(),
1303 ExcMessage("invalid iterator"));
1304 return *this;
1305}
1306
1307
1308
1309inline bool
1311 const IndexSet::IntervalAccessor &other) const
1312{
1313 Assert(index_set == other.index_set,
1314 ExcMessage(
1315 "Can not compare accessors pointing to different IndexSets"));
1316 return range_idx == other.range_idx;
1317}
1318
1319
1320
1321inline bool
1323 const IndexSet::IntervalAccessor &other) const
1324{
1325 Assert(index_set == other.index_set,
1326 ExcMessage(
1327 "Can not compare accessors pointing to different IndexSets"));
1328 return range_idx < other.range_idx;
1329}
1330
1331
1332
1333inline void
1335{
1336 Assert(
1337 is_valid(),
1338 ExcMessage(
1339 "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1340 ++range_idx;
1341
1342 // set ourselves to invalid if we walk off the end
1343 if (range_idx >= index_set->ranges.size())
1344 range_idx = numbers::invalid_dof_index;
1345}
1346
1347
1348/* IntervalIterator */
1349
1351 const IndexSet *idxset,
1352 const IndexSet::size_type range_idx)
1353 : accessor(idxset, range_idx)
1354{}
1355
1356
1357
1359 : accessor(nullptr)
1360{}
1361
1362
1363
1365 : accessor(idxset)
1366{}
1367
1368
1369
1372{
1373 accessor.advance();
1374 return *this;
1375}
1376
1377
1378
1381{
1382 const IndexSet::IntervalIterator iter = *this;
1383 accessor.advance();
1384 return iter;
1385}
1386
1387
1388
1389inline const IndexSet::IntervalAccessor &
1391{
1392 return accessor;
1393}
1394
1395
1396
1397inline const IndexSet::IntervalAccessor *
1399{
1400 return &accessor;
1401}
1402
1403
1404
1405inline bool
1407 const IndexSet::IntervalIterator &other) const
1408{
1409 return accessor == other.accessor;
1410}
1411
1412
1413
1414inline bool
1416 const IndexSet::IntervalIterator &other) const
1417{
1418 return !(*this == other);
1419}
1420
1421
1422
1423inline bool
1425 const IndexSet::IntervalIterator &other) const
1426{
1427 return accessor < other.accessor;
1428}
1429
1430
1431
1432inline int
1434 const IndexSet::IntervalIterator &other) const
1435{
1436 Assert(accessor.index_set == other.accessor.index_set,
1437 ExcMessage(
1438 "Can not compare iterators belonging to different IndexSets"));
1439
1440 const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index) ?
1441 accessor.index_set->ranges.size() :
1442 accessor.range_idx;
1443 const size_type rhs =
1445 accessor.index_set->ranges.size() :
1446 other.accessor.range_idx;
1447
1448 if (lhs > rhs)
1449 return static_cast<int>(lhs - rhs);
1450 else
1451 return -static_cast<int>(rhs - lhs);
1452}
1453
1454
1455
1456/* ElementIterator */
1457
1459 const IndexSet *idxset,
1460 const IndexSet::size_type range_idx,
1461 const IndexSet::size_type index)
1462 : index_set(idxset)
1463 , range_idx(range_idx)
1464 , idx(index)
1465{
1467 ExcMessage(
1468 "Invalid range index for IndexSet::ElementIterator constructor."));
1469 Assert(
1470 idx >= index_set->ranges[range_idx].begin &&
1473 "Invalid index argument for IndexSet::ElementIterator constructor."));
1474}
1475
1476
1477
1479 : index_set(idxset)
1480 , range_idx(numbers::invalid_dof_index)
1481 , idx(numbers::invalid_dof_index)
1482{}
1483
1484
1485
1486inline bool
1488{
1489 Assert((range_idx == numbers::invalid_dof_index &&
1491 (range_idx < index_set->ranges.size() &&
1492 idx < index_set->ranges[range_idx].end),
1493 ExcInternalError("Invalid ElementIterator state."));
1494
1495 return (range_idx < index_set->ranges.size() &&
1496 idx < index_set->ranges[range_idx].end);
1497}
1498
1499
1500
1503{
1504 Assert(
1505 is_valid(),
1506 ExcMessage(
1507 "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1508 return idx;
1509}
1510
1511
1512
1513inline bool
1515 const IndexSet::ElementIterator &other) const
1516{
1517 Assert(index_set == other.index_set,
1518 ExcMessage(
1519 "Can not compare iterators belonging to different IndexSets"));
1520 return range_idx == other.range_idx && idx == other.idx;
1521}
1522
1523
1524
1525inline void
1527{
1528 Assert(
1529 is_valid(),
1530 ExcMessage(
1531 "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1532 if (idx < index_set->ranges[range_idx].end)
1533 ++idx;
1534 // end of this range?
1535 if (idx == index_set->ranges[range_idx].end)
1536 {
1537 // point to first element in next interval if possible
1538 if (range_idx < index_set->ranges.size() - 1)
1539 {
1540 ++range_idx;
1541 idx = index_set->ranges[range_idx].begin;
1542 }
1543 else
1544 {
1545 // we just fell off the end, set to invalid:
1546 range_idx = numbers::invalid_dof_index;
1548 }
1549 }
1550}
1551
1552
1553
1556{
1557 advance();
1558 return *this;
1559}
1560
1561
1562
1565{
1566 const IndexSet::ElementIterator it = *this;
1567 advance();
1568 return it;
1569}
1570
1571
1572
1573inline bool
1575 const IndexSet::ElementIterator &other) const
1576{
1577 return !(*this == other);
1578}
1579
1580
1581
1582inline bool
1584 const IndexSet::ElementIterator &other) const
1585{
1586 Assert(index_set == other.index_set,
1587 ExcMessage(
1588 "Can not compare iterators belonging to different IndexSets"));
1589 return range_idx < other.range_idx ||
1590 (range_idx == other.range_idx && idx < other.idx);
1591}
1592
1593
1594
1595inline std::ptrdiff_t
1597 const IndexSet::ElementIterator &other) const
1598{
1599 Assert(index_set == other.index_set,
1600 ExcMessage(
1601 "Can not compare iterators belonging to different IndexSets"));
1602 if (*this == other)
1603 return 0;
1604 if (!(*this < other))
1605 return -(other - *this);
1606
1607 // only other can be equal to end() because of the checks above.
1608 Assert(is_valid(), ExcInternalError());
1609
1610 // Note: we now compute how far advance *this in "*this < other" to get other,
1611 // so we need to return -c at the end.
1612
1613 // first finish the current range:
1614 std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1615
1616 // now walk in steps of ranges (need to start one behind our current one):
1617 for (size_type range = range_idx + 1;
1618 range < index_set->ranges.size() && range <= other.range_idx;
1619 ++range)
1620 c += index_set->ranges[range].end - index_set->ranges[range].begin;
1621
1622 Assert(
1623 other.range_idx < index_set->ranges.size() ||
1625 ExcMessage(
1626 "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1627
1628 // We might have walked too far because we went until the end of
1629 // other.range_idx, so walk backwards to other.idx:
1631 c -= index_set->ranges[other.range_idx].end - other.idx;
1632
1633 return -c;
1634}
1635
1636
1637/* Range */
1638
1640 : begin(numbers::invalid_dof_index)
1641 , end(numbers::invalid_dof_index)
1642 , nth_index_in_set(numbers::invalid_dof_index)
1643{}
1644
1645
1646
1648 : begin(i1)
1649 , end(i2)
1650 , nth_index_in_set(numbers::invalid_dof_index)
1651{}
1652
1653
1654
1655/* IndexSet itself */
1656
1658 : is_compressed(true)
1659 , index_space_size(0)
1660 , largest_range(numbers::invalid_unsigned_int)
1661{}
1662
1663
1664
1666 : is_compressed(true)
1667 , index_space_size(size)
1668 , largest_range(numbers::invalid_unsigned_int)
1669{}
1670
1671
1672
1673inline IndexSet::IndexSet(IndexSet &&is) noexcept
1674 : ranges(std::move(is.ranges))
1675 , is_compressed(is.is_compressed)
1676 , index_space_size(is.index_space_size)
1677 , largest_range(is.largest_range)
1678{
1679 is.ranges.clear();
1680 is.is_compressed = true;
1681 is.index_space_size = 0;
1682 is.largest_range = numbers::invalid_unsigned_int;
1683
1684 compress();
1685}
1686
1687
1688
1689inline IndexSet &
1691{
1692 ranges = std::move(is.ranges);
1693 is_compressed = is.is_compressed;
1694 index_space_size = is.index_space_size;
1695 largest_range = is.largest_range;
1696
1697 is.ranges.clear();
1698 is.is_compressed = true;
1699 is.index_space_size = 0;
1700 is.largest_range = numbers::invalid_unsigned_int;
1701
1702 compress();
1703
1704 return *this;
1705}
1706
1707
1708
1711{
1712 compress();
1713 if (ranges.size() > 0)
1714 return {this, 0, ranges[0].begin};
1715 else
1716 return end();
1717}
1718
1719
1720
1723{
1724 compress();
1725 return IndexSet::ElementIterator(this);
1726}
1727
1728
1729
1732{
1733 compress();
1734 if (ranges.size() > 0)
1735 return IndexSet::IntervalIterator(this, 0);
1736 else
1737 return end_intervals();
1738}
1739
1740
1741
1744{
1745 compress();
1746 return IndexSet::IntervalIterator(this);
1747}
1748
1749
1750
1751inline void
1753{
1754 // reset so that there are no indices in the set any more; however,
1755 // as documented, the index set retains its size
1756 ranges.clear();
1757 is_compressed = true;
1759}
1760
1761
1762
1763inline void
1765{
1766 Assert(ranges.empty(),
1767 ExcMessage("This function can only be called if the current "
1768 "object does not yet contain any elements."));
1769 index_space_size = sz;
1770 is_compressed = true;
1771}
1772
1773
1774
1777{
1778 return index_space_size;
1779}
1780
1781
1782
1783inline void
1785{
1786 if (is_compressed == true)
1787 return;
1788
1789 do_compress();
1790}
1791
1792
1793
1794inline void
1796{
1797 add_range(index, index + 1);
1798}
1799
1800
1801
1802inline void
1804{
1811
1812 if (begin != end)
1813 {
1814 // the new index might be larger than the last index present in the
1815 // ranges. Then we can skip the binary search
1816 if (ranges.empty() || begin > ranges.back().end)
1817 ranges.emplace_back(begin, end);
1818 else if (begin == ranges.back().end)
1819 ranges.back().end = end;
1820 else
1822
1823 is_compressed = false;
1824 }
1825}
1826
1827
1828
1829template <typename ForwardIterator>
1830inline void
1831IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1832{
1833 if (begin == end)
1834 return;
1835
1836 // identify ranges in the given iterator range by checking whether some
1837 // indices happen to be consecutive. to avoid quadratic complexity when
1838 // calling add_range many times (as add_range() going into the middle of an
1839 // already existing range must shift entries around), we first collect a
1840 // vector of ranges.
1841 boost::container::small_vector<std::pair<size_type, size_type>, 200>
1842 tmp_ranges;
1843 bool ranges_are_sorted = true;
1844 for (ForwardIterator p = begin; p != end;)
1845 {
1846 // Starting with the current iterator 'p', find an iterator
1847 // 'q' so that the indices pointed to by the iterators in
1848 // the range [p,q) are consecutive. These indices then form
1849 // a range that is contiguous, and that can be added all
1850 // at once.
1851 const size_type begin_index = *p;
1852 size_type end_index = begin_index + 1;
1853
1854 // Start looking at the position after 'p', and keep iterating while
1855 // 'q' points to a duplicate of 'p':
1856 ForwardIterator q = p;
1857 ++q;
1858 while ((q != end) && (*q == *p))
1859 ++q;
1860
1861 // Now we know that 'q' is either past the end, or points to a value
1862 // other than 'p'. If it points to 'end_index', we are still good with
1863 // a contiguous range; then increment the end index of that range, and
1864 // move to the next iterator that is not a duplicate of what
1865 // we were just looking at:
1866 while ((q != end) && (static_cast<size_type>(*q) == end_index))
1867 {
1868 ++q;
1869 while ((q != end) && (static_cast<size_type>(*q) == end_index))
1870 ++q;
1871
1872 ++end_index;
1873 }
1874
1875 // Add this range:
1876 tmp_ranges.emplace_back(begin_index, end_index);
1877
1878 // Then move on to the next element in the input range.
1879 // If the starting index of the next go-around of the for loop is less
1880 // than the end index of the one just identified, then we will have at
1881 // least one pair of ranges that are not sorted, and consequently the
1882 // whole collection of ranges is not sorted.
1883 p = q;
1884 if ((p != end) && (static_cast<size_type>(*p) < end_index))
1885 ranges_are_sorted = false;
1886 }
1887
1888 add_ranges_internal(tmp_ranges, ranges_are_sorted);
1889}
1890
1891
1892
1893inline bool
1895{
1896 if (ranges.empty() == false)
1897 {
1898 compress();
1899
1900 // fast check whether the index is in the largest range
1902 if (index >= ranges[largest_range].begin &&
1903 index < ranges[largest_range].end)
1904 return true;
1905 else if (ranges.size() > 1)
1906 return is_element_binary_search(index);
1907 else
1908 return false;
1909 }
1910 else
1911 return false;
1912}
1913
1914
1915
1916inline bool
1918{
1919 compress();
1920 return (ranges.size() <= 1);
1921}
1922
1923
1924
1925inline bool
1927{
1928 return ranges.empty();
1929}
1930
1931
1932
1935{
1936 // make sure we have non-overlapping ranges
1937 compress();
1938
1939 size_type v = 0;
1940 if (!ranges.empty())
1941 {
1942 const Range &r = ranges.back();
1943 v = r.nth_index_in_set + r.end - r.begin;
1944 }
1945
1946#ifdef DEBUG
1947 size_type s = 0;
1948 for (const auto &range : ranges)
1949 s += (range.end - range.begin);
1950 Assert(s == v, ExcInternalError());
1951#endif
1952
1953 return v;
1954}
1955
1956
1957
1958inline unsigned int
1960{
1961 compress();
1962 return ranges.size();
1963}
1964
1965
1966
1969{
1970 Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1971
1972 compress();
1973 const std::vector<Range>::const_iterator main_range =
1974 ranges.begin() + largest_range;
1975
1976 return main_range->nth_index_in_set;
1977}
1978
1979
1980
1983{
1985
1986 compress();
1987
1988 // first check whether the index is in the largest range
1990 const auto main_range = ranges.begin() + largest_range;
1991 if (n >= main_range->nth_index_in_set &&
1992 n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1993 return main_range->begin + (n - main_range->nth_index_in_set);
1994 else
1996}
1997
1998
1999
2002{
2003 // to make this call thread-safe, compress() must not be called through this
2004 // function
2005 Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
2006 AssertIndexRange(n, size());
2007
2008 // return immediately if the index set is empty
2009 if (is_empty())
2011
2012 // check whether the index is in the largest range. use the result to
2013 // perform a one-sided binary search afterward
2016 return (n - ranges[largest_range].begin) +
2017 ranges[largest_range].nth_index_in_set;
2018 else if (ranges.size() > 1)
2020 else
2022}
2023
2024
2025
2026inline bool
2028{
2029 // If one of the two index sets has size zero, the other one has to
2030 // have size zero as well:
2031 if (size() == 0)
2032 return (is.size() == 0);
2033 if (is.size() == 0)
2034 return (size() == 0);
2035
2036 // Otherwise, they must have the same size (see the documentation):
2037 Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
2038
2039 compress();
2040 is.compress();
2041
2042 return (ranges == is.ranges);
2043}
2044
2045
2046
2047inline bool
2049{
2050 // If one of the two index sets has size zero, the other one has to
2051 // have a non-zero size for inequality:
2052 if (size() == 0)
2053 return (is.size() != 0);
2054 if (is.size() == 0)
2055 return (size() != 0);
2056
2057 // Otherwise, they must have the same size (see the documentation):
2058 Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
2059
2060 compress();
2061 is.compress();
2062
2063 return (ranges != is.ranges);
2064}
2065
2066
2067
2068template <typename Vector>
2069void
2071{
2072 Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
2073
2074 compress();
2075 // first fill all elements of the vector with zeroes.
2076 std::fill(vector.begin(), vector.end(), 0);
2077
2078 // then write ones into the elements whose indices are contained in the
2079 // index set
2080 for (const auto &range : ranges)
2081 for (size_type i = range.begin; i < range.end; ++i)
2082 vector[i] = 1;
2083}
2084
2085
2086
2087template <typename StreamType>
2088inline void
2089IndexSet::print(StreamType &out) const
2090{
2091 compress();
2092 out << '{';
2093 std::vector<Range>::const_iterator p;
2094 for (p = ranges.begin(); p != ranges.end(); ++p)
2095 {
2096 if (p->end - p->begin == 1)
2097 out << p->begin;
2098 else
2099 out << '[' << p->begin << ',' << p->end - 1 << ']';
2100
2101 if (p != --ranges.end())
2102 out << ", ";
2103 }
2104 out << '}' << std::endl;
2105}
2106
2107
2108
2109template <class Archive>
2110inline void
2111IndexSet::Range::serialize(Archive &ar, const unsigned int)
2112{
2114}
2115
2116
2117
2118template <class Archive>
2119inline void
2120IndexSet::serialize(Archive &ar, const unsigned int)
2121{
2123}
2124
2126
2127#endif
size_type operator*() const
Definition index_set.h:1502
bool operator==(const ElementIterator &) const
Definition index_set.h:1514
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition index_set.h:1458
std::ptrdiff_t difference_type
Definition index_set.h:939
ElementIterator & operator++()
Definition index_set.h:1555
bool operator<(const ElementIterator &) const
Definition index_set.h:1583
bool operator!=(const ElementIterator &) const
Definition index_set.h:1574
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition index_set.h:1596
std::forward_iterator_tag iterator_category
Definition index_set.h:937
const IndexSet * index_set
Definition index_set.h:953
size_type last() const
Definition index_set.h:1288
ElementIterator end() const
Definition index_set.h:1274
size_type n_elements() const
Definition index_set.h:1248
bool operator==(const IntervalAccessor &other) const
Definition index_set.h:1310
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition index_set.h:1298
bool operator<(const IntervalAccessor &other) const
Definition index_set.h:1322
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition index_set.h:1217
const IndexSet * index_set
Definition index_set.h:744
ElementIterator begin() const
Definition index_set.h:1265
std::ptrdiff_t difference_type
Definition index_set.h:847
const IntervalAccessor & operator*() const
Definition index_set.h:1390
IntervalIterator(const IntervalIterator &other)=default
bool operator==(const IntervalIterator &) const
Definition index_set.h:1406
int operator-(const IntervalIterator &p) const
Definition index_set.h:1433
bool operator<(const IntervalIterator &) const
Definition index_set.h:1424
std::forward_iterator_tag iterator_category
Definition index_set.h:845
IntervalIterator & operator=(const IntervalIterator &other)=default
IntervalIterator & operator++()
Definition index_set.h:1371
const IntervalAccessor * operator->() const
Definition index_set.h:1398
IntervalAccessor accessor
Definition index_set.h:855
bool operator!=(const IntervalIterator &) const
Definition index_set.h:1415
bool is_subset_of(const IndexSet &other) const
Definition index_set.cc:706
bool is_element_binary_search(const size_type local_index) const
Definition index_set.cc:801
size_type index_within_set_binary_search(const size_type global_index) const
Definition index_set.cc:854
size_type largest_range
Definition index_set.h:1125
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:1917
unsigned int n_intervals() const
Definition index_set.h:1959
Tpetra::Map< int, types::signed_global_dof_index, NodeType > make_tpetra_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition index_set.cc:957
IndexSet(const IndexSet &)=default
IntervalIterator end_intervals() const
Definition index_set.h:1743
void do_compress() const
Definition index_set.cc:139
ElementIterator at(const size_type global_index) const
Definition index_set.cc:876
size_type size() const
Definition index_set.h:1776
void fill_index_vector(std::vector< size_type > &indices) const
Definition index_set.cc:945
bool is_empty() const
Definition index_set.h:1926
size_type index_within_set(const size_type global_index) const
Definition index_set.h:2001
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
Definition index_set.cc:444
size_type n_elements() const
Definition index_set.h:1934
bool operator==(const IndexSet &is) const
Definition index_set.h:2027
void add_range_lower_bound(const Range &range)
Definition index_set.cc:594
bool is_element(const size_type index) const
Definition index_set.h:1894
void serialize(Archive &ar, const unsigned int version)
Definition index_set.h:2120
ElementIterator begin() const
Definition index_set.h:1710
size_type pop_front()
Definition index_set.cc:572
signed int value_type
Definition index_set.h:99
void set_size(const size_type size)
Definition index_set.h:1764
bool operator!=(const IndexSet &is) const
Definition index_set.h:2048
void read(std::istream &in)
Definition index_set.cc:745
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
void clear()
Definition index_set.h:1752
size_type largest_range_starting_index() const
Definition index_set.h:1968
void add_index(const size_type index)
Definition index_set.h:1795
void write(std::ostream &out) const
Definition index_set.cc:730
void block_read(std::istream &in)
Definition index_set.cc:781
bool is_compressed
Definition index_set.h:1108
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:610
IntervalIterator begin_intervals() const
Definition index_set.h:1731
IndexSet & operator=(const IndexSet &)=default
void fill_binary_vector(VectorType &vector) const
std::vector< Range > ranges
Definition index_set.h:1098
void subtract_set(const IndexSet &other)
Definition index_set.cc:473
ElementIterator end() const
Definition index_set.h:1722
Threads::Mutex compress_mutex
Definition index_set.h:1131
IndexSet complete_index_set(const IndexSet::size_type N)
Definition index_set.h:1204
size_type index_space_size
Definition index_set.h:1114
void block_write(std::ostream &out) const
Definition index_set.cc:767
IndexSet get_view(const size_type begin, const size_type end) const
Definition index_set.cc:273
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1803
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1982
std::size_t memory_consumption() const
void print(StreamType &out) const
Definition index_set.h:2089
size_type nth_index_in_set_binary_search(const size_type local_index) const
Definition index_set.cc:837
void compress() const
Definition index_set.h:1784
size_type pop_back()
Definition index_set.cc:554
std::vector< size_type > get_index_vector() const
Definition index_set.cc:926
Teuchos::RCP< Tpetra::Map< int, types::signed_global_dof_index, NodeType > > make_tpetra_map_rcp(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition index_set.cc:967
types::global_dof_index size_type
Definition index_set.h:80
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1831
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:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
static ::ExceptionBase & ExcIndexRangeType(T arg1, T arg2, T arg3)
#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)
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space > NodeType
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:1056
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition index_set.h:1062
friend bool operator==(const Range &range_1, const Range &range_2)
Definition index_set.h:1069
size_type end
Definition index_set.h:1024
size_type nth_index_in_set
Definition index_set.h:1026
static std::size_t memory_consumption()
Definition index_set.h:1075
size_type begin
Definition index_set.h:1023
void serialize(Archive &ar, const unsigned int version)
Definition index_set.h:2111
friend bool operator<(const Range &range_1, const Range &range_2)
Definition index_set.h:1048
void advance(std::tuple< I1, I2 > &t, const unsigned int n)