Reference documentation for deal.II version GIT relicensing-35-gbe64a602e6 2024-03-01 21:20:01+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 Teuchos::RCP<const Tpetra::Map<int, types::signed_global_dof_index>> map);
141# endif // DEAL_II_TRILINOS_WITH_TPETRA
142
146 explicit IndexSet(const Epetra_BlockMap &map);
147#endif // DEAL_II_WITH_TRILINOS
148
153 void
154 clear();
155
162 void
163 set_size(const size_type size);
164
173 size() const;
174
181 void
182 add_range(const size_type begin, const size_type end);
183
192 void
193 add_index(const size_type index);
194
219 template <typename ForwardIterator>
220 void
221 add_indices(const ForwardIterator &begin, const ForwardIterator &end);
222
238 void
239 add_indices(const IndexSet &other, const size_type offset = 0);
240
244 bool
245 is_element(const size_type index) const;
246
251 bool
252 is_contiguous() const;
253
263 bool
264 is_empty() const;
265
275 bool
276 is_ascending_and_one_to_one(const MPI_Comm communicator) const;
277
282 n_elements() const;
283
290 nth_index_in_set(const size_type local_index) const;
291
299 index_within_set(const size_type global_index) const;
300
309 unsigned int
310 n_intervals() const;
311
325
330 void
331 compress() const;
332
347 bool
348 operator==(const IndexSet &is) const;
349
364 bool
365 operator!=(const IndexSet &is) const;
366
374 operator&(const IndexSet &is) const;
375
393 get_view(const size_type begin, const size_type end) const;
394
441 get_view(const IndexSet &mask) const;
442
448 std::vector<IndexSet>
450 const std::vector<types::global_dof_index> &n_indices_per_block) const;
451
457 void
458 subtract_set(const IndexSet &other);
459
477 tensor_product(const IndexSet &other) const;
478
484 pop_back();
485
491 pop_front();
492
501 std::vector<size_type>
502 get_index_vector() const;
503
511 void
512 fill_index_vector(std::vector<size_type> &indices) const;
513
526 template <typename VectorType>
527 void
528 fill_binary_vector(VectorType &vector) const;
529
537 bool
538 is_subset_of(const IndexSet &other) const;
539
544 template <typename StreamType>
545 void
546 print(StreamType &out) const;
547
552 void
553 write(std::ostream &out) const;
554
559 void
560 read(std::istream &in);
561
566 void
567 block_write(std::ostream &out) const;
568
573 void
574 block_read(std::istream &in);
575
576#ifdef DEAL_II_WITH_TRILINOS
605 Epetra_Map
606 make_trilinos_map(const MPI_Comm communicator = MPI_COMM_WORLD,
607 const bool overlapping = false) const;
608
609# ifdef DEAL_II_TRILINOS_WITH_TPETRA
610 Tpetra::Map<int, types::signed_global_dof_index>
611 make_tpetra_map(const MPI_Comm communicator = MPI_COMM_WORLD,
612 const bool overlapping = false) const;
613
614 Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index>>
615 make_tpetra_map_rcp(const MPI_Comm communicator = MPI_COMM_WORLD,
616 const bool overlapping = false) const;
617# endif
618#endif
619
620#ifdef DEAL_II_WITH_PETSC
621 IS
622 make_petsc_is(const MPI_Comm communicator = MPI_COMM_WORLD) const;
623#endif
624
625
630 std::size_t
631 memory_consumption() const;
632
634 size_type,
635 << "The global index " << arg1
636 << " is not an element of this set.");
637
643 template <class Archive>
644 void
645 serialize(Archive &ar, const unsigned int version);
646
647
659 {
660 public:
665 IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
666
670 explicit IntervalAccessor(const IndexSet *idxset);
671
676 n_elements() const;
677
681 bool
682 is_valid() const;
683
688 begin() const;
689
695 end() const;
696
701 last() const;
702
703 private:
712 operator=(const IntervalAccessor &other);
713
717 bool
718 operator==(const IntervalAccessor &other) const;
722 bool
723 operator<(const IntervalAccessor &other) const;
728 void
729 advance();
734
740
741 friend class IntervalIterator;
742 };
743
749 {
750 public:
755 IntervalIterator(const IndexSet *idxset, const size_type range_idx);
756
760 explicit IntervalIterator(const IndexSet *idxset);
761
766
770 IntervalIterator(const IntervalIterator &other) = default;
771
776 operator=(const IntervalIterator &other) = default;
777
782 operator++();
783
788 operator++(int);
789
793 const IntervalAccessor &
794 operator*() const;
795
799 const IntervalAccessor *
800 operator->() const;
801
805 bool
806 operator==(const IntervalIterator &) const;
807
811 bool
812 operator!=(const IntervalIterator &) const;
813
817 bool
818 operator<(const IntervalIterator &) const;
819
826 int
827 operator-(const IntervalIterator &p) const;
828
834 using iterator_category = std::forward_iterator_tag;
836 using difference_type = std::ptrdiff_t;
839
840 private:
845 };
846
852 {
853 public:
858 ElementIterator(const IndexSet *idxset,
859 const size_type range_idx,
860 const size_type index);
861
865 explicit ElementIterator(const IndexSet *idxset);
866
872 operator*() const;
873
877 bool
878 is_valid() const;
879
884 operator++();
885
890 operator++(int);
891
895 bool
896 operator==(const ElementIterator &) const;
897
901 bool
902 operator!=(const ElementIterator &) const;
903
907 bool
908 operator<(const ElementIterator &) const;
909
918 std::ptrdiff_t
919 operator-(const ElementIterator &p) const;
920
926 using iterator_category = std::forward_iterator_tag;
928 using difference_type = std::ptrdiff_t;
931
932 private:
936 void
937 advance();
938
951 };
952
958 begin() const;
959
975 at(const size_type global_index) const;
976
982 end() const;
983
988 begin_intervals() const;
989
995 end_intervals() const;
996
1001private:
1010 struct Range
1011 {
1014
1016
1025 Range();
1026
1034 Range(const size_type i1, const size_type i2);
1035
1036 friend inline bool
1037 operator<(const Range &range_1, const Range &range_2)
1038 {
1039 return (
1040 (range_1.begin < range_2.begin) ||
1041 ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
1042 }
1043
1044 static bool
1046 {
1047 return x.end < y.end;
1048 }
1049
1050 static bool
1052 {
1053 return (x.nth_index_in_set + (x.end - x.begin) <
1054 y.nth_index_in_set + (y.end - y.begin));
1055 }
1056
1057 friend inline bool
1058 operator==(const Range &range_1, const Range &range_2)
1059 {
1060 return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
1061 }
1062
1063 static std::size_t
1065 {
1066 return sizeof(Range);
1067 }
1068
1074 template <class Archive>
1075 void
1076 serialize(Archive &ar, const unsigned int version);
1077 };
1078
1087 mutable std::vector<Range> ranges;
1088
1097 mutable bool is_compressed;
1098
1104
1115
1121
1125 void
1126 do_compress() const;
1127
1134 bool
1135 is_element_binary_search(const size_type local_index) const;
1136
1143 size_type
1144 nth_index_in_set_binary_search(const size_type local_index) const;
1145
1152 size_type
1153 index_within_set_binary_search(const size_type global_index) const;
1154
1160 void
1161 add_range_lower_bound(const Range &range);
1162
1166 void
1168 boost::container::small_vector<std::pair<size_type, size_type>, 200>
1169 &tmp_ranges,
1170 const bool ranges_are_sorted);
1171};
1172
1173
1174
1192inline IndexSet
1194{
1195 IndexSet is(N);
1196 is.add_range(0, N);
1197 is.compress();
1198 return is;
1199}
1200
1201/* ------------------ inline functions ------------------ */
1202
1203
1204/* IntervalAccessor */
1205
1207 const IndexSet *idxset,
1208 const IndexSet::size_type range_idx)
1209 : index_set(idxset)
1210 , range_idx(range_idx)
1211{
1212 Assert(range_idx < idxset->n_intervals(),
1213 ExcInternalError("Invalid range index"));
1214}
1215
1216
1217
1219 : index_set(idxset)
1220 , range_idx(numbers::invalid_dof_index)
1221{}
1222
1223
1224
1226 const IndexSet::IntervalAccessor &other)
1227 : index_set(other.index_set)
1228 , range_idx(other.range_idx)
1229{
1231 ExcMessage("invalid iterator"));
1232}
1233
1234
1235
1238{
1239 Assert(is_valid(), ExcMessage("invalid iterator"));
1240 return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
1241}
1242
1243
1244
1245inline bool
1247{
1248 return index_set != nullptr && range_idx < index_set->n_intervals();
1249}
1250
1251
1252
1255{
1256 Assert(is_valid(), ExcMessage("invalid iterator"));
1257 return {index_set, range_idx, index_set->ranges[range_idx].begin};
1258}
1259
1260
1261
1264{
1265 Assert(is_valid(), ExcMessage("invalid iterator"));
1266
1267 // point to first index in next interval unless we are the last interval.
1268 if (range_idx < index_set->ranges.size() - 1)
1269 return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
1270 else
1271 return index_set->end();
1272}
1273
1274
1275
1278{
1279 Assert(is_valid(), ExcMessage("invalid iterator"));
1280
1281 return index_set->ranges[range_idx].end - 1;
1282}
1283
1284
1285
1288{
1289 index_set = other.index_set;
1290 range_idx = other.range_idx;
1291 Assert(range_idx == numbers::invalid_dof_index || is_valid(),
1292 ExcMessage("invalid iterator"));
1293 return *this;
1294}
1295
1296
1297
1298inline bool
1300 const IndexSet::IntervalAccessor &other) const
1301{
1302 Assert(index_set == other.index_set,
1303 ExcMessage(
1304 "Can not compare accessors pointing to different IndexSets"));
1305 return range_idx == other.range_idx;
1306}
1307
1308
1309
1310inline bool
1312 const IndexSet::IntervalAccessor &other) const
1313{
1314 Assert(index_set == other.index_set,
1315 ExcMessage(
1316 "Can not compare accessors pointing to different IndexSets"));
1317 return range_idx < other.range_idx;
1318}
1319
1320
1321
1322inline void
1324{
1325 Assert(
1326 is_valid(),
1327 ExcMessage(
1328 "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1329 ++range_idx;
1330
1331 // set ourselves to invalid if we walk off the end
1332 if (range_idx >= index_set->ranges.size())
1333 range_idx = numbers::invalid_dof_index;
1334}
1335
1336
1337/* IntervalIterator */
1338
1340 const IndexSet *idxset,
1341 const IndexSet::size_type range_idx)
1342 : accessor(idxset, range_idx)
1343{}
1344
1345
1346
1348 : accessor(nullptr)
1349{}
1350
1351
1352
1354 : accessor(idxset)
1355{}
1356
1357
1358
1361{
1362 accessor.advance();
1363 return *this;
1364}
1365
1366
1367
1370{
1371 const IndexSet::IntervalIterator iter = *this;
1372 accessor.advance();
1373 return iter;
1374}
1375
1376
1377
1378inline const IndexSet::IntervalAccessor &
1380{
1381 return accessor;
1382}
1383
1384
1385
1386inline const IndexSet::IntervalAccessor *
1388{
1389 return &accessor;
1390}
1391
1392
1393
1394inline bool
1396 const IndexSet::IntervalIterator &other) const
1397{
1398 return accessor == other.accessor;
1399}
1400
1401
1402
1403inline bool
1405 const IndexSet::IntervalIterator &other) const
1406{
1407 return !(*this == other);
1408}
1409
1410
1411
1412inline bool
1414 const IndexSet::IntervalIterator &other) const
1415{
1416 return accessor < other.accessor;
1417}
1418
1419
1420
1421inline int
1423 const IndexSet::IntervalIterator &other) const
1424{
1425 Assert(accessor.index_set == other.accessor.index_set,
1426 ExcMessage(
1427 "Can not compare iterators belonging to different IndexSets"));
1428
1429 const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index) ?
1430 accessor.index_set->ranges.size() :
1431 accessor.range_idx;
1432 const size_type rhs =
1434 accessor.index_set->ranges.size() :
1435 other.accessor.range_idx;
1436
1437 if (lhs > rhs)
1438 return static_cast<int>(lhs - rhs);
1439 else
1440 return -static_cast<int>(rhs - lhs);
1441}
1442
1443
1444
1445/* ElementIterator */
1446
1448 const IndexSet *idxset,
1449 const IndexSet::size_type range_idx,
1450 const IndexSet::size_type index)
1451 : index_set(idxset)
1452 , range_idx(range_idx)
1453 , idx(index)
1454{
1455 Assert(range_idx < index_set->ranges.size(),
1456 ExcMessage(
1457 "Invalid range index for IndexSet::ElementIterator constructor."));
1458 Assert(
1459 idx >= index_set->ranges[range_idx].begin &&
1460 idx < index_set->ranges[range_idx].end,
1462 "Invalid index argument for IndexSet::ElementIterator constructor."));
1463}
1464
1465
1466
1468 : index_set(idxset)
1469 , range_idx(numbers::invalid_dof_index)
1470 , idx(numbers::invalid_dof_index)
1471{}
1472
1473
1474
1475inline bool
1477{
1478 Assert((range_idx == numbers::invalid_dof_index &&
1480 (range_idx < index_set->ranges.size() &&
1481 idx < index_set->ranges[range_idx].end),
1482 ExcInternalError("Invalid ElementIterator state."));
1483
1484 return (range_idx < index_set->ranges.size() &&
1485 idx < index_set->ranges[range_idx].end);
1486}
1487
1488
1489
1492{
1493 Assert(
1494 is_valid(),
1495 ExcMessage(
1496 "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1497 return idx;
1498}
1499
1500
1501
1502inline bool
1504 const IndexSet::ElementIterator &other) const
1505{
1506 Assert(index_set == other.index_set,
1507 ExcMessage(
1508 "Can not compare iterators belonging to different IndexSets"));
1509 return range_idx == other.range_idx && idx == other.idx;
1510}
1511
1512
1513
1514inline void
1516{
1517 Assert(
1518 is_valid(),
1519 ExcMessage(
1520 "Impossible to advance an IndexSet::ElementIterator that is invalid"));
1521 if (idx < index_set->ranges[range_idx].end)
1522 ++idx;
1523 // end of this range?
1524 if (idx == index_set->ranges[range_idx].end)
1525 {
1526 // point to first element in next interval if possible
1527 if (range_idx < index_set->ranges.size() - 1)
1528 {
1529 ++range_idx;
1530 idx = index_set->ranges[range_idx].begin;
1531 }
1532 else
1533 {
1534 // we just fell off the end, set to invalid:
1535 range_idx = numbers::invalid_dof_index;
1537 }
1538 }
1539}
1540
1541
1542
1545{
1546 advance();
1547 return *this;
1548}
1549
1550
1551
1554{
1555 const IndexSet::ElementIterator it = *this;
1556 advance();
1557 return it;
1558}
1559
1560
1561
1562inline bool
1564 const IndexSet::ElementIterator &other) const
1565{
1566 return !(*this == other);
1567}
1568
1569
1570
1571inline bool
1573 const IndexSet::ElementIterator &other) const
1574{
1575 Assert(index_set == other.index_set,
1576 ExcMessage(
1577 "Can not compare iterators belonging to different IndexSets"));
1578 return range_idx < other.range_idx ||
1579 (range_idx == other.range_idx && idx < other.idx);
1580}
1581
1582
1583
1584inline std::ptrdiff_t
1586 const IndexSet::ElementIterator &other) const
1587{
1588 Assert(index_set == other.index_set,
1589 ExcMessage(
1590 "Can not compare iterators belonging to different IndexSets"));
1591 if (*this == other)
1592 return 0;
1593 if (!(*this < other))
1594 return -(other - *this);
1595
1596 // only other can be equal to end() because of the checks above.
1597 Assert(is_valid(), ExcInternalError());
1598
1599 // Note: we now compute how far advance *this in "*this < other" to get other,
1600 // so we need to return -c at the end.
1601
1602 // first finish the current range:
1603 std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
1604
1605 // now walk in steps of ranges (need to start one behind our current one):
1606 for (size_type range = range_idx + 1;
1607 range < index_set->ranges.size() && range <= other.range_idx;
1608 ++range)
1609 c += index_set->ranges[range].end - index_set->ranges[range].begin;
1610
1611 Assert(
1612 other.range_idx < index_set->ranges.size() ||
1614 ExcMessage(
1615 "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1616
1617 // We might have walked too far because we went until the end of
1618 // other.range_idx, so walk backwards to other.idx:
1620 c -= index_set->ranges[other.range_idx].end - other.idx;
1621
1622 return -c;
1623}
1624
1625
1626/* Range */
1627
1629 : begin(numbers::invalid_dof_index)
1630 , end(numbers::invalid_dof_index)
1631 , nth_index_in_set(numbers::invalid_dof_index)
1632{}
1633
1634
1635
1637 : begin(i1)
1638 , end(i2)
1639 , nth_index_in_set(numbers::invalid_dof_index)
1640{}
1641
1642
1643
1644/* IndexSet itself */
1645
1647 : is_compressed(true)
1648 , index_space_size(0)
1649 , largest_range(numbers::invalid_unsigned_int)
1650{}
1651
1652
1653
1655 : is_compressed(true)
1656 , index_space_size(size)
1657 , largest_range(numbers::invalid_unsigned_int)
1658{}
1659
1660
1661
1662inline IndexSet::IndexSet(IndexSet &&is) noexcept
1663 : ranges(std::move(is.ranges))
1664 , is_compressed(is.is_compressed)
1665 , index_space_size(is.index_space_size)
1666 , largest_range(is.largest_range)
1667{
1668 is.ranges.clear();
1669 is.is_compressed = true;
1670 is.index_space_size = 0;
1671 is.largest_range = numbers::invalid_unsigned_int;
1672
1673 compress();
1674}
1675
1676
1677
1678inline IndexSet &
1680{
1681 ranges = std::move(is.ranges);
1682 is_compressed = is.is_compressed;
1683 index_space_size = is.index_space_size;
1684 largest_range = is.largest_range;
1685
1686 is.ranges.clear();
1687 is.is_compressed = true;
1688 is.index_space_size = 0;
1689 is.largest_range = numbers::invalid_unsigned_int;
1690
1691 compress();
1692
1693 return *this;
1694}
1695
1696
1697
1700{
1701 compress();
1702 if (ranges.size() > 0)
1703 return {this, 0, ranges[0].begin};
1704 else
1705 return end();
1706}
1707
1708
1709
1712{
1713 compress();
1714 return IndexSet::ElementIterator(this);
1715}
1716
1717
1718
1721{
1722 compress();
1723 if (ranges.size() > 0)
1724 return IndexSet::IntervalIterator(this, 0);
1725 else
1726 return end_intervals();
1727}
1728
1729
1730
1733{
1734 compress();
1735 return IndexSet::IntervalIterator(this);
1736}
1737
1738
1739
1740inline void
1742{
1743 // reset so that there are no indices in the set any more; however,
1744 // as documented, the index set retains its size
1745 ranges.clear();
1746 is_compressed = true;
1748}
1749
1750
1751
1752inline void
1754{
1755 Assert(ranges.empty(),
1756 ExcMessage("This function can only be called if the current "
1757 "object does not yet contain any elements."));
1758 index_space_size = sz;
1759 is_compressed = true;
1760}
1761
1762
1763
1766{
1767 return index_space_size;
1768}
1769
1770
1771
1772inline void
1774{
1775 if (is_compressed == true)
1776 return;
1777
1778 do_compress();
1779}
1780
1781
1782
1783inline void
1785{
1786 add_range(index, index + 1);
1787}
1788
1789
1790
1791inline void
1793{
1796 ExcIndexRangeType<size_type>(begin, 0, index_space_size));
1798 ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
1800
1801 if (begin != end)
1802 {
1803 // the new index might be larger than the last index present in the
1804 // ranges. Then we can skip the binary search
1805 if (ranges.empty() || begin > ranges.back().end)
1806 ranges.emplace_back(begin, end);
1807 else if (begin == ranges.back().end)
1808 ranges.back().end = end;
1809 else
1811
1812 is_compressed = false;
1813 }
1814}
1815
1816
1817
1818template <typename ForwardIterator>
1819inline void
1820IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
1821{
1822 if (begin == end)
1823 return;
1824
1825 // identify ranges in the given iterator range by checking whether some
1826 // indices happen to be consecutive. to avoid quadratic complexity when
1827 // calling add_range many times (as add_range() going into the middle of an
1828 // already existing range must shift entries around), we first collect a
1829 // vector of ranges.
1830 boost::container::small_vector<std::pair<size_type, size_type>, 200>
1831 tmp_ranges;
1832 bool ranges_are_sorted = true;
1833 for (ForwardIterator p = begin; p != end;)
1834 {
1835 // Starting with the current iterator 'p', find an iterator
1836 // 'q' so that the indices pointed to by the iterators in
1837 // the range [p,q) are consecutive. These indices then form
1838 // a range that is contiguous, and that can be added all
1839 // at once.
1840 const size_type begin_index = *p;
1841 size_type end_index = begin_index + 1;
1842
1843 // Start looking at the position after 'p', and keep iterating while
1844 // 'q' points to a duplicate of 'p':
1845 ForwardIterator q = p;
1846 ++q;
1847 while ((q != end) && (*q == *p))
1848 ++q;
1849
1850 // Now we know that 'q' is either past the end, or points to a value
1851 // other than 'p'. If it points to 'end_index', we are still good with
1852 // a contiguous range; then increment the end index of that range, and
1853 // move to the next iterator that is not a duplicate of what
1854 // we were just looking at:
1855 while ((q != end) && (static_cast<size_type>(*q) == end_index))
1856 {
1857 ++q;
1858 while ((q != end) && (static_cast<size_type>(*q) == end_index))
1859 ++q;
1860
1861 ++end_index;
1862 }
1863
1864 // Add this range:
1865 tmp_ranges.emplace_back(begin_index, end_index);
1866
1867 // Then move on to the next element in the input range.
1868 // If the starting index of the next go-around of the for loop is less
1869 // than the end index of the one just identified, then we will have at
1870 // least one pair of ranges that are not sorted, and consequently the
1871 // whole collection of ranges is not sorted.
1872 p = q;
1873 if ((p != end) && (static_cast<size_type>(*p) < end_index))
1874 ranges_are_sorted = false;
1875 }
1876
1877 add_ranges_internal(tmp_ranges, ranges_are_sorted);
1878}
1879
1880
1881
1882inline bool
1884{
1885 if (ranges.empty() == false)
1886 {
1887 compress();
1888
1889 // fast check whether the index is in the largest range
1891 if (index >= ranges[largest_range].begin &&
1892 index < ranges[largest_range].end)
1893 return true;
1894 else if (ranges.size() > 1)
1895 return is_element_binary_search(index);
1896 else
1897 return false;
1898 }
1899 else
1900 return false;
1901}
1902
1903
1904
1905inline bool
1907{
1908 compress();
1909 return (ranges.size() <= 1);
1910}
1911
1912
1913
1914inline bool
1916{
1917 return ranges.empty();
1918}
1919
1920
1921
1924{
1925 // make sure we have non-overlapping ranges
1926 compress();
1927
1928 size_type v = 0;
1929 if (!ranges.empty())
1930 {
1931 const Range &r = ranges.back();
1932 v = r.nth_index_in_set + r.end - r.begin;
1933 }
1934
1935#ifdef DEBUG
1936 size_type s = 0;
1937 for (const auto &range : ranges)
1938 s += (range.end - range.begin);
1939 Assert(s == v, ExcInternalError());
1940#endif
1941
1942 return v;
1943}
1944
1945
1946
1947inline unsigned int
1949{
1950 compress();
1951 return ranges.size();
1952}
1953
1954
1955
1958{
1959 Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
1960
1961 compress();
1962 const std::vector<Range>::const_iterator main_range =
1963 ranges.begin() + largest_range;
1964
1965 return main_range->nth_index_in_set;
1966}
1967
1968
1969
1972{
1974
1975 compress();
1976
1977 // first check whether the index is in the largest range
1979 const auto main_range = ranges.begin() + largest_range;
1980 if (n >= main_range->nth_index_in_set &&
1981 n < main_range->nth_index_in_set + (main_range->end - main_range->begin))
1982 return main_range->begin + (n - main_range->nth_index_in_set);
1983 else
1985}
1986
1987
1988
1991{
1992 // to make this call thread-safe, compress() must not be called through this
1993 // function
1994 Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
1995 AssertIndexRange(n, size());
1996
1997 // return immediately if the index set is empty
1998 if (is_empty())
2000
2001 // check whether the index is in the largest range. use the result to
2002 // perform a one-sided binary search afterward
2005 return (n - ranges[largest_range].begin) +
2006 ranges[largest_range].nth_index_in_set;
2007 else if (ranges.size() > 1)
2009 else
2011}
2012
2013
2014
2015inline bool
2017{
2018 // If one of the two index sets has size zero, the other one has to
2019 // have size zero as well:
2020 if (size() == 0)
2021 return (is.size() == 0);
2022 if (is.size() == 0)
2023 return (size() == 0);
2024
2025 // Otherwise, they must have the same size (see the documentation):
2026 Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
2027
2028 compress();
2029 is.compress();
2030
2031 return (ranges == is.ranges);
2032}
2033
2034
2035
2036inline bool
2038{
2039 // If one of the two index sets has size zero, the other one has to
2040 // have a non-zero size for inequality:
2041 if (size() == 0)
2042 return (is.size() != 0);
2043 if (is.size() == 0)
2044 return (size() != 0);
2045
2046 // Otherwise, they must have the same size (see the documentation):
2047 Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
2048
2049 compress();
2050 is.compress();
2051
2052 return (ranges != is.ranges);
2053}
2054
2055
2056
2057template <typename Vector>
2058void
2060{
2061 Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
2062
2063 compress();
2064 // first fill all elements of the vector with zeroes.
2065 std::fill(vector.begin(), vector.end(), 0);
2066
2067 // then write ones into the elements whose indices are contained in the
2068 // index set
2069 for (const auto &range : ranges)
2070 for (size_type i = range.begin; i < range.end; ++i)
2071 vector[i] = 1;
2072}
2073
2074
2075
2076template <typename StreamType>
2077inline void
2078IndexSet::print(StreamType &out) const
2079{
2080 compress();
2081 out << '{';
2082 std::vector<Range>::const_iterator p;
2083 for (p = ranges.begin(); p != ranges.end(); ++p)
2084 {
2085 if (p->end - p->begin == 1)
2086 out << p->begin;
2087 else
2088 out << '[' << p->begin << ',' << p->end - 1 << ']';
2089
2090 if (p != --ranges.end())
2091 out << ", ";
2092 }
2093 out << '}' << std::endl;
2094}
2095
2096
2097
2098template <class Archive>
2099inline void
2100IndexSet::Range::serialize(Archive &ar, const unsigned int)
2101{
2103}
2104
2105
2106
2107template <class Archive>
2108inline void
2109IndexSet::serialize(Archive &ar, const unsigned int)
2110{
2112}
2113
2115
2116#endif
size_type operator*() const
Definition index_set.h:1491
bool operator==(const ElementIterator &) const
Definition index_set.h:1503
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition index_set.h:1447
std::ptrdiff_t difference_type
Definition index_set.h:928
ElementIterator & operator++()
Definition index_set.h:1544
bool operator<(const ElementIterator &) const
Definition index_set.h:1572
bool operator!=(const ElementIterator &) const
Definition index_set.h:1563
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition index_set.h:1585
std::forward_iterator_tag iterator_category
Definition index_set.h:926
const IndexSet * index_set
Definition index_set.h:942
size_type last() const
Definition index_set.h:1277
ElementIterator end() const
Definition index_set.h:1263
size_type n_elements() const
Definition index_set.h:1237
bool operator==(const IntervalAccessor &other) const
Definition index_set.h:1299
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition index_set.h:1287
bool operator<(const IntervalAccessor &other) const
Definition index_set.h:1311
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition index_set.h:1206
const IndexSet * index_set
Definition index_set.h:733
ElementIterator begin() const
Definition index_set.h:1254
std::ptrdiff_t difference_type
Definition index_set.h:836
const IntervalAccessor & operator*() const
Definition index_set.h:1379
IntervalIterator(const IntervalIterator &other)=default
bool operator==(const IntervalIterator &) const
Definition index_set.h:1395
int operator-(const IntervalIterator &p) const
Definition index_set.h:1422
bool operator<(const IntervalIterator &) const
Definition index_set.h:1413
std::forward_iterator_tag iterator_category
Definition index_set.h:834
IntervalIterator & operator=(const IntervalIterator &other)=default
IntervalIterator & operator++()
Definition index_set.h:1360
const IntervalAccessor * operator->() const
Definition index_set.h:1387
IntervalAccessor accessor
Definition index_set.h:844
bool operator!=(const IntervalIterator &) const
Definition index_set.h:1404
bool is_subset_of(const IndexSet &other) const
Definition index_set.cc:703
bool is_element_binary_search(const size_type local_index) const
Definition index_set.cc:798
size_type index_within_set_binary_search(const size_type global_index) const
Definition index_set.cc:851
size_type largest_range
Definition index_set.h:1114
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:1906
unsigned int n_intervals() const
Definition index_set.h:1948
IndexSet(const IndexSet &)=default
IntervalIterator end_intervals() const
Definition index_set.h:1732
void do_compress() const
Definition index_set.cc:136
ElementIterator at(const size_type global_index) const
Definition index_set.cc:873
size_type size() const
Definition index_set.h:1765
void fill_index_vector(std::vector< size_type > &indices) const
Definition index_set.cc:942
bool is_empty() const
Definition index_set.h:1915
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1990
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
Definition index_set.cc:441
size_type n_elements() const
Definition index_set.h:1923
bool operator==(const IndexSet &is) const
Definition index_set.h:2016
void add_range_lower_bound(const Range &range)
Definition index_set.cc:591
bool is_element(const size_type index) const
Definition index_set.h:1883
void serialize(Archive &ar, const unsigned int version)
Definition index_set.h:2109
ElementIterator begin() const
Definition index_set.h:1699
size_type pop_front()
Definition index_set.cc:569
signed int value_type
Definition index_set.h:96
void set_size(const size_type size)
Definition index_set.h:1753
bool operator!=(const IndexSet &is) const
Definition index_set.h:2037
void read(std::istream &in)
Definition index_set.cc:742
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
void clear()
Definition index_set.h:1741
size_type largest_range_starting_index() const
Definition index_set.h:1957
IndexSet tensor_product(const IndexSet &other) const
Definition index_set.cc:539
void add_index(const size_type index)
Definition index_set.h:1784
void write(std::ostream &out) const
Definition index_set.cc:727
void block_read(std::istream &in)
Definition index_set.cc:778
bool is_compressed
Definition index_set.h:1097
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:607
IntervalIterator begin_intervals() const
Definition index_set.h:1720
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:953
IndexSet & operator=(const IndexSet &)=default
void fill_binary_vector(VectorType &vector) const
std::vector< Range > ranges
Definition index_set.h:1087
void subtract_set(const IndexSet &other)
Definition index_set.cc:470
ElementIterator end() const
Definition index_set.h:1711
Threads::Mutex compress_mutex
Definition index_set.h:1120
IndexSet complete_index_set(const IndexSet::size_type N)
Definition index_set.h:1193
size_type index_space_size
Definition index_set.h:1103
void block_write(std::ostream &out) const
Definition index_set.cc:764
IndexSet get_view(const size_type begin, const size_type end) const
Definition index_set.cc:270
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1792
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1971
std::size_t memory_consumption() const
void print(StreamType &out) const
Definition index_set.h:2078
size_type nth_index_in_set_binary_search(const size_type local_index) const
Definition index_set.cc:834
void compress() const
Definition index_set.h:1773
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:962
size_type pop_back()
Definition index_set.cc:551
std::vector< size_type > get_index_vector() const
Definition index_set.cc:923
types::global_dof_index size_type
Definition index_set.h:77
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1820
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:1045
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
Definition index_set.h:1051
friend bool operator==(const Range &range_1, const Range &range_2)
Definition index_set.h:1058
size_type end
Definition index_set.h:1013
size_type nth_index_in_set
Definition index_set.h:1015
static std::size_t memory_consumption()
Definition index_set.h:1064
size_type begin
Definition index_set.h:1012
void serialize(Archive &ar, const unsigned int version)
Definition index_set.h:2100
friend bool operator<(const Range &range_1, const Range &range_2)
Definition index_set.h:1037
void advance(std::tuple< I1, I2 > &t, const unsigned int n)