Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
index_set.h
1 // ---------------------------------------------------------------------
2 // @f$Id: index_set.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 2009 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__index_set_h
18 #define __deal2__index_set_h
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/utilities.h>
23 #include <vector>
24 #include <algorithm>
25 
26 #ifdef DEAL_II_WITH_TRILINOS
27 # include <Epetra_Map.h>
28 #endif
29 
30 #if defined(DEAL_II_WITH_MPI) || defined(DEAL_II_WITH_PETSC)
31 #include <mpi.h>
32 #else
33 typedef int MPI_Comm;
34 # ifndef MPI_COMM_WORLD
35 # define MPI_COMM_WORLD 0
36 # endif
37 #endif
38 
40 
60 class IndexSet
61 {
62 public:
66  IndexSet ();
67 
73  explicit IndexSet (const types::global_dof_index size);
74 
80  void clear ();
81 
93  void set_size (const types::global_dof_index size);
94 
105  types::global_dof_index size () const;
106 
113  void add_range (const types::global_dof_index begin,
114  const types::global_dof_index end);
115 
120  void add_index (const types::global_dof_index index);
121 
129  template <typename ForwardIterator>
130  void add_indices (const ForwardIterator &begin,
131  const ForwardIterator &end);
132 
152  void add_indices(const IndexSet &other,
153  const unsigned int offset = 0);
154 
160  bool is_element (const types::global_dof_index index) const;
161 
169  bool is_contiguous () const;
170 
176 
184  types::global_dof_index nth_index_in_set (const unsigned int local_index) const;
185 
197 
214  unsigned int n_intervals () const;
215 
224  void compress () const;
225 
235  bool operator == (const IndexSet &is) const;
236 
246  bool operator != (const IndexSet &is) const;
247 
259  IndexSet operator & (const IndexSet &is) const;
260 
280  const types::global_dof_index end) const;
281 
282 
290  void subtract_set (const IndexSet &other);
291 
292 
297  void fill_index_vector(std::vector<types::global_dof_index> &indices) const;
298 
320  template <typename Vector>
321  void fill_binary_vector (Vector &vector) const;
322 
328  template <class STREAM>
329  void print(STREAM &out) const;
330 
336  void write(std::ostream &out) const;
337 
344  void read(std::istream &in);
345 
352  void block_write(std::ostream &out) const;
353 
360  void block_read(std::istream &in);
361 
362 #ifdef DEAL_II_WITH_TRILINOS
363 
420  Epetra_Map make_trilinos_map (const MPI_Comm &communicator = MPI_COMM_WORLD,
421  const bool overlapping = false) const;
422 #endif
423 
424 
430  std::size_t memory_consumption () const;
431 
432  DeclException1 (ExcIndexNotPresent, types::global_dof_index,
433  << "The global index " << arg1
434  << " is not an element of this set.");
435 
440  template <class Archive>
441  void serialize (Archive &ar, const unsigned int version);
442 
443 private:
459  struct Range
460  {
463 
464  types::global_dof_index nth_index_in_set;
465 
475  Range ();
476 
483  Range (const types::global_dof_index i1,
484  const types::global_dof_index i2);
485 
486  friend
487  inline bool operator< (const Range &range_1,
488  const Range &range_2)
489  {
490  return ((range_1.begin < range_2.begin)
491  ||
492  ((range_1.begin == range_2.begin)
493  &&
494  (range_1.end < range_2.end)));
495  }
496 
497  static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
498  {
499  return x.end < y.end;
500  }
501 
502  static bool nth_index_compare (const IndexSet::Range &x,
503  const IndexSet::Range &y)
504  {
505  return (x.nth_index_in_set+(x.end-x.begin) <
506  y.nth_index_in_set+(y.end-y.begin));
507  }
508 
509  friend
510  inline bool operator== (const Range &range_1,
511  const Range &range_2)
512  {
513  return ((range_1.begin == range_2.begin)
514  &&
515  (range_1.end == range_2.end));
516  }
517 
518  std::size_t memory_consumption () const
519  {
520  return sizeof(Range);
521  }
522 
527  template <class Archive>
528  void serialize (Archive &ar, const unsigned int version);
529  };
530 
545  mutable std::vector<Range> ranges;
546 
560  mutable bool is_compressed;
561 
569 
586 
591  void do_compress() const;
592 };
593 
594 
613 inline
614 IndexSet complete_index_set (const unsigned int N)
615 {
616  IndexSet is (N);
617  is.add_range(0, N);
618  return is;
619 }
620 
621 /* ------------------ inline functions ------------------ */
622 
623 inline
625  :
626  begin(numbers::invalid_dof_index),
627  end(numbers::invalid_dof_index)
628 {}
629 
630 
631 inline
633  const types::global_dof_index i2)
634  :
635  begin(i1),
636  end(i2)
637 {}
638 
639 
640 inline
642  :
643  is_compressed (true),
644  index_space_size (0),
645  largest_range (deal_II_numbers::invalid_unsigned_int)
646 {}
647 
648 
649 
650 inline
652  :
653  is_compressed (true),
654  index_space_size (size),
655  largest_range (deal_II_numbers::invalid_unsigned_int)
656 {}
657 
658 
659 
660 inline
661 void
663 {
664  ranges.clear ();
665  largest_range = 0;
666  is_compressed = true;
667 }
668 
669 
670 inline
671 void
673 {
674  Assert (ranges.empty(),
675  ExcMessage ("This function can only be called if the current "
676  "object does not yet contain any elements."));
677  index_space_size = sz;
678  is_compressed = true;
679 }
680 
681 
682 
683 inline
686 {
687  return index_space_size;
688 }
689 
690 
691 
692 inline
693 void
695 {
696  if (is_compressed == true)
697  return;
698 
699  do_compress();
700 }
701 
702 
703 
704 inline
705 void
707  const types::global_dof_index end)
708 {
709  Assert ((begin < index_space_size)
710  ||
711  ((begin == index_space_size) && (end == index_space_size)),
712  ExcIndexRangeType<types::global_dof_index> (begin, 0, index_space_size));
713  Assert (end <= index_space_size,
714  ExcIndexRangeType<types::global_dof_index> (end, 0, index_space_size+1));
715  Assert (begin <= end,
716  ExcIndexRangeType<types::global_dof_index> (begin, 0, end));
717 
718  if (begin != end)
719  {
720  const Range new_range(begin,end);
721 
722  // the new index might be larger than the last
723  // index present in the ranges. Then we can
724  // skip the binary search
725  if (ranges.size() == 0 || begin > ranges.back().end)
726  ranges.push_back(new_range);
727  else
728  ranges.insert (Utilities::lower_bound (ranges.begin(),
729  ranges.end(),
730  new_range),
731  new_range);
732  is_compressed = false;
733  }
734 }
735 
736 
737 
738 inline
739 void
741 {
742  Assert (index < index_space_size,
743  ExcIndexRangeType<types::global_dof_index> (index, 0, index_space_size));
744 
745  const Range new_range(index, index+1);
746  if (ranges.size() == 0 || index > ranges.back().end)
747  ranges.push_back(new_range);
748  else if (index == ranges.back().end)
749  ranges.back().end++;
750  else
751  ranges.insert (Utilities::lower_bound (ranges.begin(),
752  ranges.end(),
753  new_range),
754  new_range);
755  is_compressed = false;
756 }
757 
758 
759 
760 template <typename ForwardIterator>
761 inline
762 void
763 IndexSet::add_indices (const ForwardIterator &begin,
764  const ForwardIterator &end)
765 {
766  // insert each element of the
767  // range. if some of them happen to
768  // be consecutive, merge them to a
769  // range
770  for (ForwardIterator p=begin; p!=end;)
771  {
772  const types::global_dof_index begin_index = *p;
773  types::global_dof_index end_index = begin_index + 1;
774  ForwardIterator q = p;
775  ++q;
776  while ((q != end) && (*q == end_index))
777  {
778  ++end_index;
779  ++q;
780  }
781 
782  add_range (begin_index, end_index);
783  p = q;
784  }
785 }
786 
787 
788 
789 inline
790 void
792  const unsigned int offset)
793 {
794  if ((this == &other) && (offset == 0))
795  return;
796 
797  for (std::vector<Range>::iterator range = other.ranges.begin();
798  range != other.ranges.end();
799  ++range)
800  {
801  add_range(range->begin+offset, range->end+offset);
802  }
803 
804  compress();
805 }
806 
807 
808 
809 inline
810 bool
812 {
813  if (ranges.empty() == false)
814  {
815  compress ();
816 
817  // fast check whether the index is in the
818  // largest range
820  if (index >= ranges[largest_range].begin &&
821  index < ranges[largest_range].end)
822  return true;
823 
824  // get the element after which
825  // we would have to insert a
826  // range that consists of all
827  // elements from this element
828  // to the end of the index
829  // range plus one. after this
830  // call we know that if
831  // p!=end() then
832  // p->begin<=index unless there
833  // is no such range at all
834  //
835  // if the searched for element
836  // is an element of this range,
837  // then we're done. otherwise,
838  // the element can't be in one
839  // of the following ranges
840  // because otherwise p would be
841  // a different iterator
842  //
843  // since we already know the position
844  // relative to the largest range (we
845  // called compress!), we can perform
846  // the binary search on ranges with
847  // lower/higher number compared to the
848  // largest range
849  std::vector<Range>::const_iterator
850  p = std::upper_bound (ranges.begin() + (index<ranges[largest_range].begin?
851  0 : largest_range+1),
852  index<ranges[largest_range].begin ?
853  ranges.begin() + largest_range:
854  ranges.end(),
855  Range (index, size()+1));
856 
857  if (p == ranges.begin())
858  return ((index >= p->begin) && (index < p->end));
859 
860  Assert ((p == ranges.end()) || (p->begin > index),
861  ExcInternalError());
862 
863  // now move to that previous
864  // range
865  --p;
866  Assert (p->begin <= index, ExcInternalError());
867 
868  return (p->end > index);
869  }
870 
871  // didn't find this index, so it's
872  // not in the set
873  return false;
874 }
875 
876 
877 
878 inline
879 bool
881 {
882  compress ();
883  return (ranges.size() <= 1);
884 }
885 
886 
887 
888 inline
891 {
892  // make sure we have
893  // non-overlapping ranges
894  compress ();
895 
897  if (!ranges.empty())
898  {
899  Range &r = ranges.back();
900  v = r.nth_index_in_set + r.end - r.begin;
901  }
902 
903 #ifdef DEBUG
905  for (std::vector<Range>::iterator range = ranges.begin();
906  range != ranges.end();
907  ++range)
908  s += (range->end - range->begin);
909  Assert(s==v, ExcInternalError());
910 #endif
911 
912  return v;
913 }
914 
915 
916 
917 inline
919 IndexSet::nth_index_in_set (const unsigned int n) const
920 {
921  // to make this call thread-safe, compress()
922  // must not be called through this function
923  Assert (is_compressed == true, ExcMessage ("IndexSet must be compressed."));
924  Assert (n < n_elements(), ExcIndexRangeType<types::global_dof_index> (n, 0, n_elements()));
925 
926  // first check whether the index is in the
927  // largest range
929  std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
930  if (n>=main_range->nth_index_in_set &&
931  n<main_range->nth_index_in_set+(main_range->end-main_range->begin))
932  return main_range->begin + (n-main_range->nth_index_in_set);
933 
934  // find out which chunk the local index n
935  // belongs to by using a binary search. the
936  // comparator is based on the end of the
937  // ranges. Use the position relative to main_range to
938  // subdivide the ranges
939  Range r (n,n+1);
940  r.nth_index_in_set = n;
941  std::vector<Range>::const_iterator range_begin, range_end;
942  if (n<main_range->nth_index_in_set)
943  {
944  range_begin = ranges.begin();
945  range_end = main_range;
946  }
947  else
948  {
949  range_begin = main_range + 1;
950  range_end = ranges.end();
951  }
952 
953  std::vector<Range>::const_iterator
954  p = Utilities::lower_bound(range_begin, range_end, r,
955  Range::nth_index_compare);
956 
957  if (p != ranges.end())
958  return p->begin + (n-p->nth_index_in_set);
959  else
960  {
961  Assert (false, ExcInternalError());
963  }
964 }
965 
966 
967 
968 inline
971 {
972  // to make this call thread-safe, compress()
973  // must not be called through this function
974  Assert (is_compressed == true, ExcMessage ("IndexSet must be compressed."));
975  Assert (is_element(n) == true, ExcIndexNotPresent (n));
976  Assert (n < size(), ExcIndexRangeType<types::global_dof_index> (n, 0, size()));
977 
978  // check whether the index is in the largest
979  // range. use the result to perform a
980  // one-sided binary search afterward
982  std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
983  if (n >= main_range->begin && n < main_range->end)
984  return (n-main_range->begin) + main_range->nth_index_in_set;
985 
986  Range r(n, n);
987  std::vector<Range>::const_iterator range_begin, range_end;
988  if (n<main_range->begin)
989  {
990  range_begin = ranges.begin();
991  range_end = main_range;
992  }
993  else
994  {
995  range_begin = main_range + 1;
996  range_end = ranges.end();
997  }
998 
999  std::vector<Range>::const_iterator
1000  p = Utilities::lower_bound(range_begin, range_end, r,
1001  Range::end_compare);
1002 
1003  Assert(p!=ranges.end(), ExcInternalError());
1004  Assert(p->begin<=n, ExcInternalError());
1005  Assert(n<p->end, ExcInternalError());
1006  return (n-p->begin) + p->nth_index_in_set;
1007 }
1008 
1009 
1010 
1011 inline
1012 bool
1014 {
1015  Assert (size() == is.size(),
1016  ExcDimensionMismatch (size(), is.size()));
1017 
1018  compress ();
1019  is.compress ();
1020 
1021  return ranges == is.ranges;
1022 }
1023 
1024 
1025 
1026 inline
1027 bool
1029 {
1030  Assert (size() == is.size(),
1031  ExcDimensionMismatch (size(), is.size()));
1032 
1033  compress ();
1034  is.compress ();
1035 
1036  return ranges != is.ranges;
1037 }
1038 
1039 
1040 
1041 template <typename Vector>
1042 void
1044 {
1045  Assert (vector.size() == size(),
1046  ExcDimensionMismatch (vector.size(), size()));
1047 
1048  compress();
1049  // first fill all elements of the vector
1050  // with zeroes.
1051  std::fill (vector.begin(), vector.end(), 0);
1052 
1053  // then write ones into the elements whose
1054  // indices are contained in the index set
1055  for (std::vector<Range>::iterator it = ranges.begin();
1056  it != ranges.end();
1057  ++it)
1058  for (types::global_dof_index i=it->begin; i<it->end; ++i)
1059  vector[i] = 1;
1060 }
1061 
1062 
1063 
1064 template <class STREAM>
1065 inline
1066 void
1067 IndexSet::print (STREAM &out) const
1068 {
1069  compress();
1070  out << "{";
1071  std::vector<Range>::const_iterator p;
1072  for (p = ranges.begin(); p != ranges.end(); ++p)
1073  {
1074  if (p->end-p->begin==1)
1075  out << p->begin;
1076  else
1077  out << "[" << p->begin << "," << p->end-1 << "]";
1078 
1079  if (p !=--ranges.end())
1080  out << ", ";
1081  }
1082  out << "}" << std::endl;
1083 }
1084 
1085 
1086 
1087 template <class Archive>
1088 inline
1089 void
1090 IndexSet::Range::serialize (Archive &ar, const unsigned int)
1091 {
1092  ar &begin &end &nth_index_in_set;
1093 }
1094 
1095 
1096 
1097 template <class Archive>
1098 inline
1099 void
1100 IndexSet::serialize (Archive &ar, const unsigned int)
1101 {
1103 }
1104 
1105 DEAL_II_NAMESPACE_CLOSE
1106 
1107 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:731
types::global_dof_index index_within_set(const types::global_dof_index global_index) const
Definition: index_set.h:970
types::global_dof_index size() const
Definition: index_set.h:685
void block_read(std::istream &in)
::ExceptionBase & ExcMessage(std::string arg1)
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:763
bool operator==(const IndexSet &is) const
Definition: index_set.h:1013
types::global_dof_index largest_range
Definition: index_set.h:585
void add_index(const types::global_dof_index index)
Definition: index_set.h:740
void fill_binary_vector(Vector &vector) const
Definition: index_set.h:1043
void block_write(std::ostream &out) const
bool is_contiguous() const
Definition: index_set.h:880
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1090
unsigned int n_intervals() const
IndexSet complete_index_set(const unsigned int N)
Definition: index_set.h:614
void clear()
Definition: index_set.h:662
void serialize(Archive &ar, const unsigned int version)
Definition: index_set.h:1100
void write(std::ostream &out) const
std::vector< Range > ranges
Definition: index_set.h:545
unsigned int global_dof_index
Definition: types.h:100
bool operator!=(const IndexSet &is) const
Definition: index_set.h:1028
void subtract_set(const IndexSet &other)
#define Assert(cond, exc)
Definition: exceptions.h:299
types::global_dof_index nth_index_in_set(const unsigned int local_index) const
Definition: index_set.h:919
types::global_dof_index index_space_size
Definition: index_set.h:568
IndexSet operator&(const IndexSet &is) const
void set_size(const types::global_dof_index size)
Definition: index_set.h:672
void add_range(const types::global_dof_index begin, const types::global_dof_index end)
Definition: index_set.h:706
iterator begin()
void compress() const
Definition: index_set.h:694
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
iterator end()
IndexSet()
Definition: index_set.h:641
IndexSet get_view(const types::global_dof_index begin, const types::global_dof_index end) const
void fill_index_vector(std::vector< types::global_dof_index > &indices) const
void do_compress() const
bool is_compressed
Definition: index_set.h:560
std::size_t memory_consumption() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const types::global_dof_index invalid_dof_index
Definition: types.h:217
::ExceptionBase & ExcInternalError()
void read(std::istream &in)
void print(STREAM &out) const
Definition: index_set.h:1067
bool is_element(const types::global_dof_index index) const
Definition: index_set.h:811
std::size_t size() const
types::global_dof_index n_elements() const
Definition: index_set.h:890