00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__index_set_h
00013 #define __deal2__index_set_h
00014
00015 #include <deal.II/base/config.h>
00016 #include <deal.II/base/utilities.h>
00017 #include <deal.II/base/exceptions.h>
00018
00019 #include <vector>
00020 #include <algorithm>
00021
00022 #ifdef DEAL_II_USE_TRILINOS
00023 # include <Epetra_Map.h>
00024 #endif
00025
00026 #if defined(DEAL_II_COMPILER_SUPPORTS_MPI) || defined(DEAL_II_USE_PETSC)
00027 #include <mpi.h>
00028 #else
00029 typedef int MPI_Comm;
00030 #define MPI_COMM_WORLD 0
00031 #endif
00032
00033 DEAL_II_NAMESPACE_OPEN
00034
00054 class IndexSet
00055 {
00056 public:
00060 IndexSet ();
00061
00067 IndexSet (const unsigned int size);
00068
00078 void set_size (const unsigned int size);
00079
00090 unsigned int size () const;
00091
00098 void add_range (const unsigned int begin,
00099 const unsigned int end);
00100
00105 void add_index (const unsigned int index);
00106
00114 template <typename ForwardIterator>
00115 void add_indices (const ForwardIterator &begin,
00116 const ForwardIterator &end);
00117
00123 void add_indices(const IndexSet & other);
00124
00130 bool is_element (const unsigned int index) const;
00131
00139 bool is_contiguous () const;
00140
00145 unsigned int n_elements () const;
00146
00154 unsigned int nth_index_in_set (const unsigned int local_index) const;
00155
00166 unsigned int index_within_set (const unsigned int global_index) const;
00167
00184 unsigned int n_intervals () const;
00185
00194 void compress () const;
00195
00205 bool operator == (const IndexSet &is) const;
00206
00216 bool operator != (const IndexSet &is) const;
00217
00229 IndexSet operator & (const IndexSet &is) const;
00230
00239 IndexSet get_view (const unsigned int begin,
00240 const unsigned int end) const;
00241
00242
00250 void subtract_set (const IndexSet & other);
00251
00252
00257 void fill_index_vector(std::vector<unsigned int> & indices) const;
00258
00259
00265 template <class STREAM>
00266 void print(STREAM &out) const;
00267
00273 void write(std::ostream & out) const;
00274
00281 void read(std::istream & in);
00282
00289 void block_write(std::ostream & out) const;
00290
00297 void block_read(std::istream & in);
00298
00299 #ifdef DEAL_II_USE_TRILINOS
00300
00357 Epetra_Map make_trilinos_map (const MPI_Comm &communicator = MPI_COMM_WORLD,
00358 const bool overlapping = false) const;
00359 #endif
00360
00361
00367 std::size_t memory_consumption () const;
00368
00369 DeclException1 (ExcIndexNotPresent, int,
00370 << "The global index " << arg1
00371 << " is not an element of this set.");
00372
00377 template <class Archive>
00378 void serialize (Archive & ar, const unsigned int version);
00379
00380 private:
00396 struct Range
00397 {
00398 unsigned int begin;
00399 unsigned int end;
00400
00401 unsigned int nth_index_in_set;
00402
00412 Range ();
00413
00420 Range (const unsigned int i1,
00421 const unsigned int i2);
00422
00423 friend
00424 inline bool operator< (const Range &range_1,
00425 const Range &range_2)
00426 {
00427 return ((range_1.begin < range_2.begin)
00428 ||
00429 ((range_1.begin == range_2.begin)
00430 &&
00431 (range_1.end < range_2.end)));
00432 }
00433
00434 static bool end_compare(const IndexSet::Range & x, const IndexSet::Range & y)
00435 {
00436 return x.end < y.end;
00437 }
00438
00439 static bool nth_index_compare (const IndexSet::Range & x,
00440 const IndexSet::Range & y)
00441 {
00442 return (x.nth_index_in_set+(x.end-x.begin) <
00443 y.nth_index_in_set+(y.end-y.begin));
00444 }
00445
00446 friend
00447 inline bool operator== (const Range &range_1,
00448 const Range &range_2)
00449 {
00450 return ((range_1.begin == range_2.begin)
00451 ||
00452 (range_1.begin == range_2.begin));
00453 }
00454
00455 std::size_t memory_consumption () const
00456 {
00457 return sizeof(Range);
00458 }
00459
00464 template <class Archive>
00465 void serialize (Archive & ar, const unsigned int version);
00466 };
00467
00482 mutable std::vector<Range> ranges;
00483
00497 mutable bool is_compressed;
00498
00505 unsigned int index_space_size;
00506
00522 mutable unsigned int largest_range;
00523
00528 void do_compress() const;
00529 };
00530
00531
00532
00533
00534 inline
00535 IndexSet::Range::Range ()
00536 :
00537 begin(numbers::invalid_unsigned_int),
00538 end(numbers::invalid_unsigned_int)
00539 {}
00540
00541
00542 inline
00543 IndexSet::Range::Range (const unsigned int i1,
00544 const unsigned int i2)
00545 :
00546 begin(i1),
00547 end(i2)
00548 {}
00549
00550
00551 inline
00552 IndexSet::IndexSet ()
00553 :
00554 is_compressed (true),
00555 index_space_size (0),
00556 largest_range (deal_II_numbers::invalid_unsigned_int)
00557 {}
00558
00559
00560
00561 inline
00562 IndexSet::IndexSet (const unsigned int size)
00563 :
00564 is_compressed (true),
00565 index_space_size (size),
00566 largest_range (deal_II_numbers::invalid_unsigned_int)
00567 {}
00568
00569
00570
00571 inline
00572 void
00573 IndexSet::set_size (const unsigned int sz)
00574 {
00575 Assert (ranges.empty(),
00576 ExcMessage ("This function can only be called if the current "
00577 "object does not yet contain any elements."));
00578 index_space_size = sz;
00579 is_compressed = true;
00580 }
00581
00582
00583
00584 inline
00585 unsigned int
00586 IndexSet::size () const
00587 {
00588 return index_space_size;
00589 }
00590
00591
00592
00593 inline
00594 void
00595 IndexSet::compress () const
00596 {
00597 if (is_compressed == true)
00598 return;
00599
00600 do_compress();
00601 }
00602
00603
00604
00605 inline
00606 void
00607 IndexSet::add_range (const unsigned int begin,
00608 const unsigned int end)
00609 {
00610 Assert ((begin < index_space_size)
00611 ||
00612 ((begin == index_space_size) && (end == index_space_size)),
00613 ExcIndexRange (begin, 0, index_space_size));
00614 Assert (end <= index_space_size,
00615 ExcIndexRange (end, 0, index_space_size+1));
00616 Assert (begin <= end,
00617 ExcIndexRange (begin, 0, end));
00618
00619 if (begin != end)
00620 {
00621 const Range new_range(begin,end);
00622
00623
00624
00625
00626 if (ranges.size() == 0 || begin > ranges.back().end)
00627 ranges.push_back(new_range);
00628 else
00629 ranges.insert (Utilities::lower_bound (ranges.begin(),
00630 ranges.end(),
00631 new_range),
00632 new_range);
00633 is_compressed = false;
00634 }
00635 }
00636
00637
00638
00639 inline
00640 void
00641 IndexSet::add_index (const unsigned int index)
00642 {
00643 Assert (index < index_space_size,
00644 ExcIndexRange (index, 0, index_space_size));
00645
00646 const Range new_range(index, index+1);
00647 if (ranges.size() == 0 || index > ranges.back().end)
00648 ranges.push_back(new_range);
00649 else if (index == ranges.back().end)
00650 ranges.back().end++;
00651 else
00652 ranges.insert (Utilities::lower_bound (ranges.begin(),
00653 ranges.end(),
00654 new_range),
00655 new_range);
00656 is_compressed = false;
00657 }
00658
00659
00660
00661 template <typename ForwardIterator>
00662 inline
00663 void
00664 IndexSet::add_indices (const ForwardIterator &begin,
00665 const ForwardIterator &end)
00666 {
00667
00668
00669
00670
00671 for (ForwardIterator p=begin; p!=end;)
00672 {
00673 const unsigned int begin_index = *p;
00674 unsigned int end_index = begin_index + 1;
00675 ForwardIterator q = p;
00676 ++q;
00677 while ((q != end) && (*q == end_index))
00678 {
00679 ++end_index;
00680 ++q;
00681 }
00682
00683 add_range (begin_index, end_index);
00684 p = q;
00685 }
00686 }
00687
00688
00689
00690 inline
00691 void
00692 IndexSet::add_indices(const IndexSet & other)
00693 {
00694 if (this == &other)
00695 return;
00696
00697 for (std::vector<Range>::iterator range = other.ranges.begin();
00698 range != other.ranges.end();
00699 ++range)
00700 {
00701 add_range(range->begin, range->end);
00702 }
00703
00704 compress();
00705 }
00706
00707
00708
00709 inline
00710 bool
00711 IndexSet::is_element (const unsigned int index) const
00712 {
00713 if (ranges.empty() == false)
00714 {
00715 compress ();
00716
00717
00718
00719 Assert (largest_range < ranges.size(), ExcInternalError());
00720 if (index >= ranges[largest_range].begin &&
00721 index < ranges[largest_range].end)
00722 return true;
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749 std::vector<Range>::const_iterator
00750 p = std::upper_bound (ranges.begin() + (index<ranges[largest_range].begin?
00751 0 : largest_range+1),
00752 index<ranges[largest_range].begin ?
00753 ranges.begin() + largest_range:
00754 ranges.end(),
00755 Range (index, size()+1));
00756
00757 if (p == ranges.begin())
00758 return ((index >= p->begin) && (index < p->end));
00759
00760 Assert ((p == ranges.end()) || (p->begin > index),
00761 ExcInternalError());
00762
00763
00764
00765 --p;
00766 Assert (p->begin <= index, ExcInternalError());
00767
00768 return (p->end > index);
00769 }
00770
00771
00772
00773 return false;
00774 }
00775
00776
00777
00778 inline
00779 bool
00780 IndexSet::is_contiguous () const
00781 {
00782 compress ();
00783 return (ranges.size() <= 1);
00784 }
00785
00786
00787
00788 inline
00789 unsigned int
00790 IndexSet::n_elements () const
00791 {
00792
00793
00794 compress ();
00795
00796 unsigned int v = 0;
00797 if (!ranges.empty())
00798 {
00799 Range & r = ranges.back();
00800 v = r.nth_index_in_set + r.end - r.begin;
00801 }
00802
00803 #ifdef DEBUG
00804 unsigned int s = 0;
00805 for (std::vector<Range>::iterator range = ranges.begin();
00806 range != ranges.end();
00807 ++range)
00808 s += (range->end - range->begin);
00809 Assert(s==v, ExcInternalError());
00810 #endif
00811
00812 return v;
00813 }
00814
00815
00816
00817 inline
00818 unsigned int
00819 IndexSet::nth_index_in_set (const unsigned int n) const
00820 {
00821
00822
00823 Assert (is_compressed == true, ExcMessage ("IndexSet must be compressed."));
00824 Assert (n < n_elements(), ExcIndexRange (n, 0, n_elements()));
00825
00826
00827
00828 Assert (largest_range < ranges.size(), ExcInternalError());
00829 std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
00830 if (n>=main_range->nth_index_in_set &&
00831 n<main_range->nth_index_in_set+(main_range->end-main_range->begin))
00832 return main_range->begin + (n-main_range->nth_index_in_set);
00833
00834
00835
00836
00837
00838
00839 Range r (n,n+1);
00840 r.nth_index_in_set = n;
00841 std::vector<Range>::const_iterator range_begin, range_end;
00842 if (n<main_range->nth_index_in_set)
00843 {
00844 range_begin = ranges.begin();
00845 range_end = main_range;
00846 }
00847 else
00848 {
00849 range_begin = main_range + 1;
00850 range_end = ranges.end();
00851 }
00852
00853 std::vector<Range>::const_iterator
00854 p = Utilities::lower_bound(range_begin, range_end, r,
00855 Range::nth_index_compare);
00856
00857 if (p != ranges.end())
00858 return p->begin + (n-p->nth_index_in_set);
00859 else
00860 {
00861 Assert (false, ExcInternalError());
00862 return numbers::invalid_unsigned_int;
00863 }
00864 }
00865
00866
00867
00868 inline
00869 unsigned int
00870 IndexSet::index_within_set (const unsigned int n) const
00871 {
00872
00873
00874 Assert (is_compressed == true, ExcMessage ("IndexSet must be compressed."));
00875 Assert (is_element(n) == true, ExcIndexNotPresent (n));
00876 Assert (n < size(), ExcIndexRange (n, 0, size()));
00877
00878
00879
00880
00881 Assert (largest_range < ranges.size(), ExcInternalError());
00882 std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
00883 if (n >= main_range->begin && n < main_range->end)
00884 return (n-main_range->begin) + main_range->nth_index_in_set;
00885
00886 Range r(n, n);
00887 std::vector<Range>::const_iterator range_begin, range_end;
00888 if (n<main_range->begin)
00889 {
00890 range_begin = ranges.begin();
00891 range_end = main_range;
00892 }
00893 else
00894 {
00895 range_begin = main_range + 1;
00896 range_end = ranges.end();
00897 }
00898
00899 std::vector<Range>::const_iterator
00900 p = Utilities::lower_bound(range_begin, range_end, r,
00901 Range::end_compare);
00902
00903 Assert(p!=ranges.end(), ExcInternalError());
00904 Assert(p->begin<=n, ExcInternalError());
00905 Assert(n<p->end, ExcInternalError());
00906 return (n-p->begin) + p->nth_index_in_set;
00907 }
00908
00909
00910
00911 inline
00912 bool
00913 IndexSet::operator == (const IndexSet &is) const
00914 {
00915 Assert (size() == is.size(),
00916 ExcDimensionMismatch (size(), is.size()));
00917
00918 compress ();
00919 is.compress ();
00920
00921 return ranges == is.ranges;
00922 }
00923
00924
00925
00926 inline
00927 bool
00928 IndexSet::operator != (const IndexSet &is) const
00929 {
00930 Assert (size() == is.size(),
00931 ExcDimensionMismatch (size(), is.size()));
00932
00933 compress ();
00934 is.compress ();
00935
00936 return ranges != is.ranges;
00937 }
00938
00939 template <class STREAM>
00940 inline
00941 void
00942 IndexSet::print (STREAM &out) const
00943 {
00944 compress();
00945 out << "{";
00946 std::vector<Range>::const_iterator p;
00947 for (p = ranges.begin(); p != ranges.end(); ++p)
00948 {
00949 if (p->end-p->begin==1)
00950 out << p->begin;
00951 else
00952 out << "[" << p->begin << "," << p->end-1 << "]";
00953
00954 if (p !=--ranges.end())
00955 out << ", ";
00956 }
00957 out << "}" << std::endl;
00958 }
00959
00960
00961
00962 template <class Archive>
00963 inline
00964 void
00965 IndexSet::Range::serialize (Archive & ar, const unsigned int)
00966 {
00967 ar & begin & end & nth_index_in_set;
00968 }
00969
00970
00971
00972 template <class Archive>
00973 inline
00974 void
00975 IndexSet::serialize (Archive & ar, const unsigned int)
00976 {
00977 ar & ranges & is_compressed & index_space_size & largest_range;
00978 }
00979
00980 DEAL_II_NAMESPACE_CLOSE
00981
00982 #endif
00983
00984