Reference documentation for deal.II version Git 05e4468 2017-09-21 10:18:23 +0200
index_set.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2017 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_index_set_h
17 #define dealii_index_set_h
18 
19 #include <deal.II/base/config.h>
20 #include <deal.II/base/utilities.h>
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/thread_management.h>
23 #include <boost/serialization/vector.hpp>
24 #include <vector>
25 #include <algorithm>
26 #include <iterator>
27 
28 
29 #ifdef DEAL_II_WITH_TRILINOS
30 # include <Epetra_Map.h>
31 #endif
32 
33 #if defined(DEAL_II_WITH_MPI) || defined(DEAL_II_WITH_PETSC)
34 #include <mpi.h>
35 #else
36 typedef int MPI_Comm;
37 # ifndef MPI_COMM_WORLD
38 # define MPI_COMM_WORLD 0
39 # endif
40 #endif
41 
42 DEAL_II_NAMESPACE_OPEN
43 
68 class IndexSet
69 {
70 public:
71  // forward declarations:
72  class ElementIterator;
74 
80 
98  typedef signed int value_type;
99 
100 
104  IndexSet ();
105 
109  explicit IndexSet (const size_type size);
110 
114  IndexSet (const IndexSet &) = default;
115 
119  IndexSet &operator= (const IndexSet &) = default;
120 
125  IndexSet (IndexSet &&is);
126 
132 
133 #ifdef DEAL_II_WITH_TRILINOS
134 
137  explicit IndexSet(const Epetra_Map &map);
138 #endif
139 
144  void clear ();
145 
152  void set_size (const size_type size);
153 
161  size_type size () const;
162 
169  void add_range (const size_type begin,
170  const size_type end);
171 
175  void add_index (const size_type index);
176 
186  template <typename ForwardIterator>
187  void add_indices (const ForwardIterator &begin,
188  const ForwardIterator &end);
189 
205  void add_indices(const IndexSet &other,
206  const unsigned int offset = 0);
207 
211  bool is_element (const size_type index) const;
212 
217  bool is_contiguous () const;
218 
223  bool is_empty () const;
224 
234  bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const;
235 
239  size_type n_elements () const;
240 
246  size_type nth_index_in_set (const unsigned int local_index) const;
247 
254  size_type index_within_set (const size_type global_index) const;
255 
264  unsigned int n_intervals () const;
265 
270  unsigned int largest_range_starting_index() const;
271 
276  void compress () const;
277 
283  bool operator == (const IndexSet &is) const;
284 
290  bool operator != (const IndexSet &is) const;
291 
298  IndexSet operator & (const IndexSet &is) const;
299 
313  IndexSet get_view (const size_type begin,
314  const size_type end) const;
315 
321  void subtract_set (const IndexSet &other);
322 
327  size_type pop_back ();
328 
333  size_type pop_front ();
334 
338  void fill_index_vector(std::vector<size_type> &indices) const;
339 
352  template <typename VectorType>
353  void fill_binary_vector (VectorType &vector) const;
354 
359  template <class StreamType>
360  void print(StreamType &out) const;
361 
366  void write(std::ostream &out) const;
367 
372  void read(std::istream &in);
373 
378  void block_write(std::ostream &out) const;
379 
384  void block_read(std::istream &in);
385 
386 #ifdef DEAL_II_WITH_TRILINOS
387 
415  Epetra_Map make_trilinos_map (const MPI_Comm &communicator = MPI_COMM_WORLD,
416  const bool overlapping = false) const;
417 #endif
418 
419 
424  std::size_t memory_consumption () const;
425 
427  << "The global index " << arg1
428  << " is not an element of this set.");
429 
434  template <class Archive>
435  void serialize (Archive &ar, const unsigned int version);
436 
437 
449  {
450  public:
455  IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
456 
460  explicit IntervalAccessor(const IndexSet *idxset);
461 
465  size_type n_elements() const;
466 
470  bool is_valid() const;
471 
475  ElementIterator begin() const;
476 
481  ElementIterator end() const;
482 
486  size_type last() const;
487 
488  private:
492  IntervalAccessor(const IntervalAccessor &other);
497 
501  bool operator == (const IntervalAccessor &other) const;
505  bool operator < (const IntervalAccessor &other) const;
510  void advance ();
515 
520  size_type range_idx;
521 
522  friend class IntervalIterator;
523  };
524 
529  class IntervalIterator : public std::iterator<std::forward_iterator_tag,IntervalAccessor>
530  {
531  public:
536  IntervalIterator(const IndexSet *idxset, const size_type range_idx);
537 
541  explicit IntervalIterator(const IndexSet *idxset);
542 
547 
551  IntervalIterator(const IntervalIterator &other);
552 
557 
562 
567 
571  const IntervalAccessor &operator* () const;
572 
576  const IntervalAccessor *operator-> () const;
577 
581  bool operator == (const IntervalIterator &) const;
582 
586  bool operator != (const IntervalIterator &) const;
587 
591  bool operator < (const IntervalIterator &) const;
592 
599  int operator - (const IntervalIterator &p) const;
600 
601  private:
606  };
607 
612  class ElementIterator : public std::iterator<std::forward_iterator_tag,size_type>
613  {
614  public:
619  ElementIterator(const IndexSet *idxset,
620  const size_type range_idx,
621  const size_type index);
622 
626  explicit ElementIterator(const IndexSet *idxset);
627 
632  size_type operator* () const;
633 
637  bool is_valid () const;
638 
643 
648 
652  bool operator == (const ElementIterator &) const;
653 
657  bool operator != (const ElementIterator &) const;
658 
662  bool operator < (const ElementIterator &) const;
663 
672  std::ptrdiff_t operator - (const ElementIterator &p) const;
673 
674  private:
678  void advance ();
679 
687  size_type range_idx;
691  size_type idx;
692  };
693 
698  ElementIterator begin() const;
699 
714  ElementIterator at(const size_type global_index) const;
715 
720  ElementIterator end() const;
721 
726 
732 
737 private:
746  struct Range
747  {
748  size_type begin;
749  size_type end;
750 
751  size_type nth_index_in_set;
752 
761  Range ();
762 
770  Range (const size_type i1,
771  const size_type i2);
772 
773  friend
774  inline bool operator< (const Range &range_1,
775  const Range &range_2)
776  {
777  return ((range_1.begin < range_2.begin)
778  ||
779  ((range_1.begin == range_2.begin)
780  &&
781  (range_1.end < range_2.end)));
782  }
783 
784  static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
785  {
786  return x.end < y.end;
787  }
788 
789  static bool nth_index_compare (const IndexSet::Range &x,
790  const IndexSet::Range &y)
791  {
792  return (x.nth_index_in_set+(x.end-x.begin) <
793  y.nth_index_in_set+(y.end-y.begin));
794  }
795 
796  friend
797  inline bool operator== (const Range &range_1,
798  const Range &range_2)
799  {
800  return ((range_1.begin == range_2.begin)
801  &&
802  (range_1.end == range_2.end));
803  }
804 
805  static std::size_t memory_consumption ()
806  {
807  return sizeof(Range);
808  }
809 
814  template <class Archive>
815  void serialize (Archive &ar, const unsigned int version);
816  };
817 
826  mutable std::vector<Range> ranges;
827 
836  mutable bool is_compressed;
837 
842  size_type index_space_size;
843 
853  mutable size_type largest_range;
854 
860 
864  void do_compress() const;
865 };
866 
867 
885 inline
886 IndexSet complete_index_set (const unsigned int N)
887 {
888  IndexSet is (N);
889  is.add_range(0, N);
890  is.compress();
891  return is;
892 }
893 
894 /* ------------------ inline functions ------------------ */
895 
896 
897 /* IntervalAccessor */
898 
899 inline
902  const IndexSet::size_type range_idx)
903  :
904  index_set(idxset),
905  range_idx(range_idx)
906 {
907  Assert(range_idx < idxset->n_intervals(), ExcInternalError("Invalid range index"));
908 }
909 
910 
911 
912 inline
915  :
916  index_set(idxset),
917  range_idx(numbers::invalid_dof_index)
918 {}
919 
920 
921 
922 inline
925  :
926  index_set (other.index_set),
927  range_idx(other.range_idx)
928 {
930  ExcMessage("invalid iterator"));
931 }
932 
933 
934 
935 inline
938 {
939  Assert(is_valid(), ExcMessage("invalid iterator"));
940  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
941 }
942 
943 
944 
945 inline
946 bool
948 {
949  return index_set != nullptr && range_idx < index_set->n_intervals();
950 }
951 
952 
953 
954 inline
957 {
958  Assert(is_valid(), ExcMessage("invalid iterator"));
959  return IndexSet::ElementIterator(index_set, range_idx,
960  index_set->ranges[range_idx].begin);
961 }
962 
963 
964 
965 inline
968 {
969  Assert(is_valid(), ExcMessage("invalid iterator"));
970 
971  // point to first index in next interval unless we are the last interval.
972  if (range_idx < index_set->ranges.size()-1)
973  return IndexSet::ElementIterator(index_set, range_idx+1, index_set->ranges[range_idx+1].begin);
974  else
975  return index_set->end();
976 }
977 
978 
979 
980 inline
983 {
984  Assert(is_valid(), ExcMessage("invalid iterator"));
985 
986  return index_set->ranges[range_idx].end-1;
987 }
988 
989 
990 
991 inline
994 {
995  index_set = other.index_set;
996  range_idx = other.range_idx;
997  Assert( range_idx == numbers::invalid_dof_index || is_valid(),
998  ExcMessage("invalid iterator"));
999  return *this;
1000 }
1001 
1002 
1003 
1004 inline
1005 bool
1007 {
1008  Assert (index_set == other.index_set,
1009  ExcMessage("Can not compare accessors pointing to different IndexSets"));
1010  return range_idx == other.range_idx;
1011 }
1012 
1013 
1014 
1015 inline
1016 bool
1018 {
1019  Assert (index_set == other.index_set,
1020  ExcMessage("Can not compare accessors pointing to different IndexSets"));
1021  return range_idx < other.range_idx;
1022 }
1023 
1024 
1025 
1026 inline
1028 {
1029  Assert (is_valid(),
1030  ExcMessage("Impossible to advance an IndexSet::IntervalIterator that is invalid"));
1031  ++range_idx;
1032 
1033  // set ourselves to invalid if we walk off the end
1034  if (range_idx>=index_set->ranges.size())
1035  range_idx = numbers::invalid_dof_index;
1036 }
1037 
1038 
1039 /* IntervalIterator */
1040 
1041 inline
1043  const IndexSet::size_type range_idx)
1044  :
1045  accessor(idxset, range_idx)
1046 {}
1047 
1048 
1049 
1050 inline
1052  :
1053  accessor(nullptr)
1054 {}
1055 
1056 
1057 
1058 inline
1060  :
1061  accessor(idxset)
1062 {}
1063 
1064 
1065 
1066 inline
1068  :
1069  accessor(other.accessor)
1070 {}
1071 
1072 
1073 
1074 inline
1077 {
1078  accessor = other.accessor;
1079  return *this;
1080 }
1081 
1082 
1083 
1084 inline
1087 {
1088  accessor.advance();
1089  return *this;
1090 }
1091 
1092 
1093 
1094 inline
1097 {
1098  const IndexSet::IntervalIterator iter = *this;
1099  accessor.advance ();
1100  return iter;
1101 }
1102 
1103 
1104 
1105 inline
1108 {
1109  return accessor;
1110 }
1111 
1112 
1113 
1114 inline
1117 {
1118  return &accessor;
1119 }
1120 
1121 
1122 
1123 inline
1124 bool
1126 {
1127  return accessor == other.accessor;
1128 }
1129 
1130 
1131 
1132 inline
1133 bool
1135 {
1136  return !(*this == other);
1137 }
1138 
1139 
1140 
1141 inline
1142 bool
1144 {
1145  return accessor < other.accessor;
1146 }
1147 
1148 
1149 
1150 inline
1151 int
1153 {
1154  Assert (accessor.index_set == other.accessor.index_set,
1155  ExcMessage("Can not compare iterators belonging to different IndexSets"));
1156 
1157  const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index)
1158  ?
1159  accessor.index_set->ranges.size()
1160  :
1161  accessor.range_idx;
1163  ?
1164  accessor.index_set->ranges.size()
1165  :
1166  other.accessor.range_idx;
1167 
1168  if (lhs > rhs)
1169  return static_cast<int>(lhs - rhs);
1170  else
1171  return -static_cast<int>(rhs - lhs);
1172 }
1173 
1174 
1175 
1176 /* ElementIterator */
1177 
1178 inline
1181  const IndexSet::size_type range_idx,
1182  const IndexSet::size_type index)
1183  :
1184  index_set(idxset),
1185  range_idx(range_idx),
1186  idx(index)
1187 {
1188  Assert(range_idx < index_set->ranges.size(),
1189  ExcMessage("Invalid range index for IndexSet::ElementIterator constructor."));
1190  Assert(idx >= index_set->ranges[range_idx].begin
1191  &&
1192  idx < index_set->ranges[range_idx].end,
1193  ExcInternalError("Invalid index argument for IndexSet::ElementIterator constructor."));
1194 }
1195 
1196 
1197 
1198 inline
1200  :
1201  index_set(idxset),
1202  range_idx(numbers::invalid_dof_index),
1203  idx(numbers::invalid_dof_index)
1204 {}
1205 
1206 
1207 
1208 inline
1209 bool
1211 {
1213  ||
1214  (range_idx < index_set->ranges.size() && idx<index_set->ranges[range_idx].end),
1215  ExcInternalError("Invalid ElementIterator state."));
1216 
1217  return (range_idx < index_set->ranges.size() &&
1218  idx<index_set->ranges[range_idx].end);
1219 }
1220 
1221 
1222 
1223 
1224 inline
1227 {
1228  Assert (is_valid(),
1229  ExcMessage("Impossible to dereference an IndexSet::ElementIterator that is invalid"));
1230  return idx;
1231 }
1232 
1233 
1234 
1235 inline
1236 bool
1238 {
1239  Assert (index_set == other.index_set,
1240  ExcMessage("Can not compare iterators belonging to different IndexSets"));
1241  return range_idx == other.range_idx && idx==other.idx;
1242 }
1243 
1244 
1245 
1246 inline
1247 void
1249 {
1250  Assert (is_valid(),
1251  ExcMessage("Impossible to advance an IndexSet::ElementIterator that is invalid"));
1252  if (idx < index_set->ranges[range_idx].end)
1253  ++idx;
1254  // end of this range?
1255  if (idx == index_set->ranges[range_idx].end)
1256  {
1257  // point to first element in next interval if possible
1258  if (range_idx < index_set->ranges.size()-1)
1259  {
1260  ++range_idx;
1261  idx = index_set->ranges[range_idx].begin;
1262  }
1263  else
1264  {
1265  // we just fell off the end, set to invalid:
1266  range_idx = numbers::invalid_dof_index;
1268  }
1269  }
1270 }
1271 
1272 
1273 
1274 inline
1277 {
1278  advance();
1279  return *this;
1280 }
1281 
1282 
1283 
1284 inline
1287 {
1288  const IndexSet::ElementIterator it = *this;
1289  advance();
1290  return it;
1291 }
1292 
1293 
1294 
1295 inline
1297 {
1298  return !(*this == other);
1299 }
1300 
1301 
1302 
1303 inline
1305 {
1306  Assert (index_set == other.index_set,
1307  ExcMessage("Can not compare iterators belonging to different IndexSets"));
1308  return range_idx < other.range_idx || (range_idx == other.range_idx && idx<other.idx);
1309 }
1310 
1311 
1312 
1313 inline
1314 std::ptrdiff_t
1316 {
1317  Assert (index_set == other.index_set,
1318  ExcMessage("Can not compare iterators belonging to different IndexSets"));
1319  if (*this == other)
1320  return 0;
1321  if (!(*this < other))
1322  return -(other-*this);
1323 
1324  // only other can be equal to end() because of the checks above.
1325  Assert (is_valid(), ExcInternalError());
1326 
1327  // Note: we now compute how far advance *this in "*this < other" to get other, so we need to return -c at the end.
1328 
1329  // first finish the current range:
1330  std::ptrdiff_t c = index_set->ranges[range_idx].end-idx;
1331 
1332  // now walk in steps of ranges (need to start one behind our current one):
1333  for (size_type range=range_idx+1; range<index_set->ranges.size() && range<=other.range_idx; ++range)
1334  c += index_set->ranges[range].end-index_set->ranges[range].begin;
1335 
1336  Assert(other.range_idx < index_set->ranges.size() || other.range_idx == numbers::invalid_dof_index,
1337  ExcMessage("Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
1338 
1339  // We might have walked too far because we went until the end of other.range_idx, so walk backwards to other.idx:
1341  c -= index_set->ranges[other.range_idx].end - other.idx;
1342 
1343  return -c;
1344 }
1345 
1346 
1347 /* Range */
1348 
1349 inline
1351  :
1352  begin(numbers::invalid_dof_index),
1353  end(numbers::invalid_dof_index),
1354  nth_index_in_set(numbers::invalid_dof_index)
1355 {}
1356 
1357 
1358 
1359 inline
1361  const size_type i2)
1362  :
1363  begin(i1),
1364  end(i2),
1365  nth_index_in_set(numbers::invalid_dof_index)
1366 {}
1367 
1368 
1369 
1370 /* IndexSet itself */
1371 
1372 inline
1374  :
1375  is_compressed (true),
1376  index_space_size (0),
1377  largest_range (numbers::invalid_unsigned_int)
1378 {}
1379 
1380 
1381 
1382 inline
1384  :
1385  is_compressed (true),
1386  index_space_size (size),
1387  largest_range (numbers::invalid_unsigned_int)
1388 {}
1389 
1390 
1391 
1392 inline
1394  :
1395  ranges (std::move(is.ranges)),
1396  is_compressed (is.is_compressed),
1397  index_space_size (is.index_space_size),
1398  largest_range (is.largest_range)
1399 {
1400  is.ranges.clear ();
1401  is.is_compressed = true;
1402  is.index_space_size = 0;
1403  is.largest_range = numbers::invalid_unsigned_int;
1404 
1405  compress ();
1406 }
1407 
1408 
1409 
1410 inline
1412 {
1413  ranges = std::move (is.ranges);
1414  is_compressed = is.is_compressed;
1415  index_space_size = is.index_space_size;
1416  largest_range = is.largest_range;
1417 
1418  is.ranges.clear ();
1419  is.is_compressed = true;
1420  is.index_space_size = 0;
1421  is.largest_range = numbers::invalid_unsigned_int;
1422 
1423  compress ();
1424 
1425  return *this;
1426 }
1427 
1428 
1429 
1430 inline
1432 {
1433  compress();
1434  if (ranges.size()>0)
1435  return IndexSet::ElementIterator(this, 0, ranges[0].begin);
1436  else
1437  return end();
1438 }
1439 
1440 
1441 
1442 inline
1444 {
1445  compress();
1446  Assert (global_index < size(),
1447  ExcIndexRangeType<size_type> (global_index, 0, size()));
1448 
1449  if (ranges.empty())
1450  return end();
1451 
1452  std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
1453 
1454  Range r (global_index, global_index+1);
1455  // This optimization makes the bounds for lower_bound smaller by checking
1456  // the largest range first.
1457  std::vector<Range>::const_iterator range_begin, range_end;
1458  if (global_index<main_range->begin)
1459  {
1460  range_begin = ranges.begin();
1461  range_end = main_range;
1462  }
1463  else
1464  {
1465  range_begin = main_range;
1466  range_end = ranges.end();
1467  }
1468 
1469  // This will give us the first range p=[a,b[ with b>=global_index using
1470  // a binary search
1471  const std::vector<Range>::const_iterator
1472  p = Utilities::lower_bound(range_begin, range_end, r, Range::end_compare);
1473 
1474  // We couldn't find a range, which means we have no range that contains
1475  // global_index and also no range behind it, meaning we need to return end().
1476  if (p == ranges.end())
1477  return end();
1478 
1479  // Finally, we can have two cases: Either global_index is not in [a,b[,
1480  // which means we need to return an iterator to a because global_index, ...,
1481  // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in
1482  // [a,b[ and we will return an iterator pointing directly at global_index
1483  // (else branch).
1484  if (global_index < p->begin)
1485  return IndexSet::ElementIterator(this, p-ranges.begin(), p->begin);
1486  else
1487  return IndexSet::ElementIterator(this, p-ranges.begin(), global_index);
1488 }
1489 
1490 
1491 
1492 inline
1495 {
1496  compress();
1497  return IndexSet::ElementIterator(this);
1498 }
1499 
1500 
1501 
1502 inline
1505 {
1506  compress();
1507  if (ranges.size()>0)
1508  return IndexSet::IntervalIterator(this, 0);
1509  else
1510  return end_intervals();
1511 }
1512 
1513 
1514 
1515 inline
1518 {
1519  compress();
1520  return IndexSet::IntervalIterator(this);
1521 }
1522 
1523 
1524 
1525 inline
1526 void
1528 {
1529  // reset so that there are no indices in the set any more; however,
1530  // as documented, the index set retains its size
1531  ranges.clear ();
1532  is_compressed = true;
1534 }
1535 
1536 
1537 
1538 inline
1539 void
1541 {
1542  Assert (ranges.empty(),
1543  ExcMessage ("This function can only be called if the current "
1544  "object does not yet contain any elements."));
1545  index_space_size = sz;
1546  is_compressed = true;
1547 }
1548 
1549 
1550 
1551 inline
1554 {
1555  return index_space_size;
1556 }
1557 
1558 
1559 
1560 inline
1561 void
1563 {
1564  if (is_compressed == true)
1565  return;
1566 
1567  do_compress();
1568 }
1569 
1570 
1571 
1572 inline
1573 void
1575 {
1576  Assert (index < index_space_size,
1577  ExcIndexRangeType<size_type> (index, 0, index_space_size));
1578 
1579  const Range new_range(index, index+1);
1580  if (ranges.size() == 0 || index > ranges.back().end)
1581  ranges.push_back(new_range);
1582  else if (index == ranges.back().end)
1583  ranges.back().end++;
1584  else
1585  ranges.insert (Utilities::lower_bound (ranges.begin(),
1586  ranges.end(),
1587  new_range),
1588  new_range);
1589  is_compressed = false;
1590 }
1591 
1592 
1593 
1594 template <typename ForwardIterator>
1595 inline
1596 void
1597 IndexSet::add_indices (const ForwardIterator &begin,
1598  const ForwardIterator &end)
1599 {
1600  // insert each element of the range. if some of them happen to be
1601  // consecutive, merge them to a range
1602  for (ForwardIterator p=begin; p!=end;)
1603  {
1604  const size_type begin_index = *p;
1605  size_type end_index = begin_index + 1;
1606  ForwardIterator q = p;
1607  ++q;
1608  while ((q != end) && (*q == end_index))
1609  {
1610  ++end_index;
1611  ++q;
1612  }
1613 
1614  add_range (begin_index, end_index);
1615  p = q;
1616  }
1617 }
1618 
1619 
1620 
1621 inline
1622 bool
1623 IndexSet::is_element (const size_type index) const
1624 {
1625  if (ranges.empty() == false)
1626  {
1627  compress ();
1628 
1629  // fast check whether the index is in the largest range
1631  if (index >= ranges[largest_range].begin &&
1632  index < ranges[largest_range].end)
1633  return true;
1634 
1635  // get the element after which we would have to insert a range that
1636  // consists of all elements from this element to the end of the index
1637  // range plus one. after this call we know that if p!=end() then
1638  // p->begin<=index unless there is no such range at all
1639  //
1640  // if the searched for element is an element of this range, then we're
1641  // done. otherwise, the element can't be in one of the following ranges
1642  // because otherwise p would be a different iterator
1643  //
1644  // since we already know the position relative to the largest range (we
1645  // called compress!), we can perform the binary search on ranges with
1646  // lower/higher number compared to the largest range
1647  std::vector<Range>::const_iterator
1648  p = std::upper_bound (ranges.begin() + (index<ranges[largest_range].begin?
1649  0 : largest_range+1),
1650  index<ranges[largest_range].begin ?
1651  ranges.begin() + largest_range:
1652  ranges.end(),
1653  Range (index, size()+1));
1654 
1655  if (p == ranges.begin())
1656  return ((index >= p->begin) && (index < p->end));
1657 
1658  Assert ((p == ranges.end()) || (p->begin > index),
1659  ExcInternalError());
1660 
1661  // now move to that previous range
1662  --p;
1663  Assert (p->begin <= index, ExcInternalError());
1664 
1665  return (p->end > index);
1666  }
1667 
1668  // didn't find this index, so it's not in the set
1669  return false;
1670 }
1671 
1672 
1673 
1674 inline
1675 bool
1677 {
1678  compress ();
1679  return (ranges.size() <= 1);
1680 }
1681 
1682 
1683 
1684 inline
1685 bool
1687 {
1688  return ranges.empty();
1689 }
1690 
1691 
1692 
1693 inline
1696 {
1697  // make sure we have non-overlapping ranges
1698  compress ();
1699 
1700  size_type v = 0;
1701  if (!ranges.empty())
1702  {
1703  Range &r = ranges.back();
1704  v = r.nth_index_in_set + r.end - r.begin;
1705  }
1706 
1707 #ifdef DEBUG
1708  size_type s = 0;
1709  for (std::vector<Range>::iterator range = ranges.begin();
1710  range != ranges.end();
1711  ++range)
1712  s += (range->end - range->begin);
1713  Assert(s==v, ExcInternalError());
1714 #endif
1715 
1716  return v;
1717 }
1718 
1719 
1720 
1721 inline
1722 unsigned int
1724 {
1725  compress ();
1726  return ranges.size();
1727 }
1728 
1729 
1730 
1731 inline
1732 unsigned int
1734 {
1735  Assert(ranges.empty()==false, ExcMessage("IndexSet cannot be empty."));
1736 
1737  compress();
1738  const std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
1739 
1740  return main_range->nth_index_in_set;
1741 }
1742 
1743 
1744 
1745 inline
1747 IndexSet::nth_index_in_set (const unsigned int n) const
1748 {
1749  Assert (n < n_elements(), ExcIndexRangeType<size_type> (n, 0, n_elements()));
1750 
1751  compress ();
1752 
1753  // first check whether the index is in the largest range
1755  std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
1756  if (n>=main_range->nth_index_in_set &&
1757  n<main_range->nth_index_in_set+(main_range->end-main_range->begin))
1758  return main_range->begin + (n-main_range->nth_index_in_set);
1759 
1760  // find out which chunk the local index n belongs to by using a binary
1761  // search. the comparator is based on the end of the ranges. Use the
1762  // position relative to main_range to subdivide the ranges
1763  Range r (n,n+1);
1764  r.nth_index_in_set = n;
1765  std::vector<Range>::const_iterator range_begin, range_end;
1766  if (n<main_range->nth_index_in_set)
1767  {
1768  range_begin = ranges.begin();
1769  range_end = main_range;
1770  }
1771  else
1772  {
1773  range_begin = main_range + 1;
1774  range_end = ranges.end();
1775  }
1776 
1777  const std::vector<Range>::const_iterator
1778  p = Utilities::lower_bound(range_begin, range_end, r,
1779  Range::nth_index_compare);
1780 
1781  Assert (p != ranges.end(), ExcInternalError());
1782  return p->begin + (n-p->nth_index_in_set);
1783 }
1784 
1785 
1786 
1787 inline
1790 {
1791  // to make this call thread-safe, compress() must not be called through this
1792  // function
1793  Assert (is_compressed == true, ExcMessage ("IndexSet must be compressed."));
1794  Assert (n < size(), ExcIndexRangeType<size_type> (n, 0, size()));
1795 
1796  // return immediately if the index set is empty
1797  if (is_empty())
1799 
1800  // check whether the index is in the largest range. use the result to
1801  // perform a one-sided binary search afterward
1803  std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
1804  if (n >= main_range->begin && n < main_range->end)
1805  return (n-main_range->begin) + main_range->nth_index_in_set;
1806 
1807  Range r(n, n);
1808  std::vector<Range>::const_iterator range_begin, range_end;
1809  if (n<main_range->begin)
1810  {
1811  range_begin = ranges.begin();
1812  range_end = main_range;
1813  }
1814  else
1815  {
1816  range_begin = main_range + 1;
1817  range_end = ranges.end();
1818  }
1819 
1820  std::vector<Range>::const_iterator
1821  p = Utilities::lower_bound(range_begin, range_end, r,
1822  Range::end_compare);
1823 
1824  // if n is not in this set
1825  if (p==range_end || p->end == n || p->begin > n)
1827 
1828  Assert(p!=ranges.end(), ExcInternalError());
1829  Assert(p->begin<=n, ExcInternalError());
1830  Assert(n<p->end, ExcInternalError());
1831  return (n-p->begin) + p->nth_index_in_set;
1832 }
1833 
1834 
1835 
1836 inline
1837 bool
1839 {
1840  Assert (size() == is.size(),
1841  ExcDimensionMismatch (size(), is.size()));
1842 
1843  compress ();
1844  is.compress ();
1845 
1846  return ranges == is.ranges;
1847 }
1848 
1849 
1850 
1851 inline
1852 bool
1854 {
1855  Assert (size() == is.size(),
1856  ExcDimensionMismatch (size(), is.size()));
1857 
1858  compress ();
1859  is.compress ();
1860 
1861  return ranges != is.ranges;
1862 }
1863 
1864 
1865 
1866 template <typename Vector>
1867 void
1868 IndexSet::fill_binary_vector (Vector &vector) const
1869 {
1870  Assert (vector.size() == size(),
1871  ExcDimensionMismatch (vector.size(), size()));
1872 
1873  compress();
1874  // first fill all elements of the vector with zeroes.
1875  std::fill (vector.begin(), vector.end(), 0);
1876 
1877  // then write ones into the elements whose indices are contained in the
1878  // index set
1879  for (std::vector<Range>::iterator it = ranges.begin();
1880  it != ranges.end();
1881  ++it)
1882  for (size_type i=it->begin; i<it->end; ++i)
1883  vector[i] = 1;
1884 }
1885 
1886 
1887 
1888 template <class StreamType>
1889 inline
1890 void
1891 IndexSet::print (StreamType &out) const
1892 {
1893  compress();
1894  out << "{";
1895  std::vector<Range>::const_iterator p;
1896  for (p = ranges.begin(); p != ranges.end(); ++p)
1897  {
1898  if (p->end-p->begin==1)
1899  out << p->begin;
1900  else
1901  out << "[" << p->begin << "," << p->end-1 << "]";
1902 
1903  if (p != --ranges.end())
1904  out << ", ";
1905  }
1906  out << "}" << std::endl;
1907 }
1908 
1909 
1910 
1911 template <class Archive>
1912 inline
1913 void
1914 IndexSet::Range::serialize (Archive &ar, const unsigned int)
1915 {
1916  ar &begin &end &nth_index_in_set;
1917 }
1918 
1919 
1920 
1921 template <class Archive>
1922 inline
1923 void
1924 IndexSet::serialize (Archive &ar, const unsigned int)
1925 {
1927 }
1928 
1929 DEAL_II_NAMESPACE_CLOSE
1930 
1931 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:626
bool operator<(const ElementIterator &) const
Definition: index_set.h:1304
static const unsigned int invalid_unsigned_int
Definition: types.h:173
void add_index(const size_type index)
Definition: index_set.h:1574
IntervalIterator begin_intervals() const
Definition: index_set.h:1504
size_type n_elements() const
Definition: index_set.h:937
const IndexSet * index_set
Definition: index_set.h:683
IndexSet & operator=(const IndexSet &)=default
bool is_valid() const
Definition: index_set.h:1210
void block_read(std::istream &in)
Definition: index_set.cc:507
signed int value_type
Definition: index_set.h:98
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:618
size_type n_elements() const
Definition: index_set.h:1695
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1597
void fill_binary_vector(VectorType &vector) const
bool operator==(const IndexSet &is) const
Definition: index_set.h:1838
STL namespace.
bool operator==(const IntervalIterator &) const
Definition: index_set.h:1125
bool is_empty() const
Definition: index_set.h:1686
size_type size() const
Definition: index_set.h:1553
ElementIterator(const IndexSet *idxset, const size_type range_idx, const size_type index)
Definition: index_set.h:1180
iterator end()
void print(StreamType &out) const
Definition: index_set.h:1891
bool operator==(const IntervalAccessor &other) const
Definition: index_set.h:1006
size_type index_space_size
Definition: index_set.h:842
ElementIterator begin() const
Definition: index_set.h:956
void block_write(std::ostream &out) const
Definition: index_set.cc:492
bool is_contiguous() const
Definition: index_set.h:1676
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1914
unsigned int n_intervals() const
Definition: index_set.h:1723
IndexSet complete_index_set(const unsigned int N)
Definition: index_set.h:886
void clear()
Definition: index_set.h:1527
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1924
static::ExceptionBase & ExcMessage(std::string arg1)
void write(std::ostream &out) const
Definition: index_set.cc:455
size_type operator*() const
Definition: index_set.h:1226
IntervalIterator & operator=(const IntervalIterator &other)
Definition: index_set.h:1076
std::vector< Range > ranges
Definition: index_set.h:826
static::ExceptionBase & ExcIndexNotPresent(size_type arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:593
unsigned int global_dof_index
Definition: types.h:88
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1853
void subtract_set(const IndexSet &other)
Definition: index_set.cc:288
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:526
#define Assert(cond, exc)
Definition: exceptions.h:337
unsigned int largest_range_starting_index() const
Definition: index_set.h:1733
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const IntervalAccessor & operator*() const
Definition: index_set.h:1107
IntervalAccessor accessor
Definition: index_set.h:605
IndexSet operator&(const IndexSet &is) const
Definition: index_set.cc:199
std::size_t size() const
ElementIterator at(const size_type global_index) const
Definition: index_set.h:1443
IntervalAccessor & operator=(const IntervalAccessor &other)
Definition: index_set.h:993
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:254
void compress() const
Definition: index_set.h:1562
bool operator!=(const IntervalIterator &) const
Definition: index_set.h:1134
iterator begin()
const IndexSet * index_set
Definition: index_set.h:514
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:549
size_type pop_back()
Definition: index_set.cc:362
ElementIterator begin() const
Definition: index_set.h:1431
void add_range(const size_type begin, const size_type end)
Definition: index_set.cc:98
void set_size(const size_type size)
Definition: index_set.h:1540
ElementIterator & operator++()
Definition: index_set.h:1276
IntervalAccessor(const IndexSet *idxset, const size_type range_idx)
Definition: index_set.h:901
bool operator<(const IntervalIterator &) const
Definition: index_set.h:1143
ElementIterator end() const
Definition: index_set.h:967
Threads::Mutex compress_mutex
Definition: index_set.h:859
bool operator<(const IntervalAccessor &other) const
Definition: index_set.h:1017
void do_compress() const
Definition: index_set.cc:130
bool operator!=(const ElementIterator &) const
Definition: index_set.h:1296
size_type pop_front()
Definition: index_set.cc:379
bool is_compressed
Definition: index_set.h:836
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1789
bool is_element(const size_type index) const
Definition: index_set.h:1623
std::size_t memory_consumption() const
Definition: index_set.cc:688
ElementIterator end() const
Definition: index_set.h:1494
bool is_valid() const
Definition: index_set.h:947
const types::global_dof_index invalid_dof_index
Definition: types.h:187
size_type nth_index_in_set(const unsigned int local_index) const
Definition: index_set.h:1747
bool operator==(const ElementIterator &) const
Definition: index_set.h:1237
size_type last() const
Definition: index_set.h:982
types::global_dof_index size_type
Definition: index_set.h:73
int operator-(const IntervalIterator &p) const
Definition: index_set.h:1152
void read(std::istream &in)
Definition: index_set.cc:470
std::ptrdiff_t operator-(const ElementIterator &p) const
Definition: index_set.h:1315
size_type largest_range
Definition: index_set.h:853
const IntervalAccessor * operator->() const
Definition: index_set.h:1116
IntervalIterator & operator++()
Definition: index_set.h:1086
IntervalIterator end_intervals() const
Definition: index_set.h:1517
static::ExceptionBase & ExcInternalError()