From f98484f7e540ce59b32212dc9a4e7dc6d69975cd Mon Sep 17 00:00:00 2001 From: kronbichler Date: Thu, 16 Jan 2014 12:40:23 +0000 Subject: [PATCH] Fix quadratic complexity in number of ranges in two methods. Do not directly insert into an IndexSet with MGTools::extract_inner_interface_dofs. git-svn-id: https://svn.dealii.org/trunk@32224 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 7 + deal.II/include/deal.II/base/index_set.h | 557 +++++++---------------- deal.II/source/base/index_set.cc | 323 +++++++------ deal.II/source/multigrid/mg_tools.cc | 44 +- 4 files changed, 398 insertions(+), 533 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 4c21c7596b..54e3ae9172 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -85,6 +85,13 @@ inconvenience this causes.
    +
  1. Fixed: The methods IndexSet::do_compress() and + IndexSet::add_indices(IndexSet&) had quadratic complexity in the number of + ranges. The algorithms have been changed into linear complexity ones. +
    + (Martin Kronbichler, 2014/01/15) +
    2. +
  2. Fixed: There were several bugs in functions like FEValues::get_function_values() where the code did not properly handle the case of FE_Nothing. This is now fixed. diff --git a/deal.II/include/deal.II/base/index_set.h b/deal.II/include/deal.II/base/index_set.h index b5b439d990..166a642a23 100644 --- a/deal.II/include/deal.II/base/index_set.h +++ b/deal.II/include/deal.II/base/index_set.h @@ -66,356 +66,243 @@ public: IndexSet (); /** - * Constructor that also sets the - * overall size of the index - * range. + * Constructor that also sets the overall size of the index range. */ explicit IndexSet (const types::global_dof_index size); /** - * Remove all indices from this - * index set. The index set retains - * its size, however. + * Remove all indices from this index set. The index set retains its size, + * however. */ void clear (); /** - * Set the maximal size of the - * indices upon which this object - * operates. + * Set the maximal size of the indices upon which this object operates. * - * This function can only be - * called if the index set does - * not yet contain any elements. - * This can be achieved by calling - * clear(), for example. + * This function can only be called if the index set does not yet contain + * any elements. This can be achieved by calling clear(), for example. */ void set_size (const types::global_dof_index size); /** - * Return the size of the index - * space of which this index set - * is a subset of. + * Return the size of the index space of which this index set is a subset + * of. * - * Note that the result is not equal to - * the number of indices within this - * set. The latter information is - * returned by n_elements(). + * Note that the result is not equal to the number of indices within this + * set. The latter information is returned by n_elements(). */ types::global_dof_index size () const; /** - * Add the half-open range - * $[\text{begin},\text{end})$ to - * the set of indices represented - * by this class. + * Add the half-open range $[\text{begin},\text{end})$ to the set of indices + * represented by this class. */ void add_range (const types::global_dof_index begin, const types::global_dof_index end); /** - * Add an individual index to the - * set of indices. + * Add an individual index to the set of indices. */ void add_index (const types::global_dof_index index); /** - * Add a whole set of indices - * described by dereferencing - * every element of the the - * iterator range - * [begin,end). + * Add a whole set of indices described by dereferencing every element of + * the the iterator range [begin,end). */ template void add_indices (const ForwardIterator &begin, const ForwardIterator &end); /** - * Add the given IndexSet @p other to the - * current one, constructing the union of - * *this and @p other. + * Add the given IndexSet @p other to the current one, constructing the + * union of *this and @p other. * - * If the @p offset argument is nonzero, then every - * index in @p other is shifted by @p offset before being - * added to the current index set. This allows to construct, - * for example, one index set from several others that are - * supposed to represent index sets corresponding to - * different ranges (e.g., when constructing the set of - * nonzero entries of a block vector from the sets of nonzero - * elements of the individual blocks of a vector). + * If the @p offset argument is nonzero, then every index in @p other is + * shifted by @p offset before being added to the current index set. This + * allows to construct, for example, one index set from several others that + * are supposed to represent index sets corresponding to different ranges + * (e.g., when constructing the set of nonzero entries of a block vector + * from the sets of nonzero elements of the individual blocks of a vector). * - * This function will generate an exception if any of the - * (possibly shifted) indices of the @p other index set - * lie outside the range [0,size()) represented - * by the current object. + * This function will generate an exception if any of the (possibly shifted) + * indices of the @p other index set lie outside the range + * [0,size()) represented by the current object. */ void add_indices(const IndexSet &other, const unsigned int offset = 0); /** - * Return whether the specified - * index is an element of the - * index set. + * Return whether the specified index is an element of the index set. */ bool is_element (const types::global_dof_index index) const; /** - * Return whether the index set - * stored by this object defines - * a contiguous range. This is - * true also if no indices are - * stored at all. + * Return whether the index set stored by this object defines a contiguous + * range. This is true also if no indices are stored at all. */ bool is_contiguous () const; /** - * Return the number of elements - * stored in this index set. + * Return the number of elements stored in this index set. */ types::global_dof_index n_elements () const; /** - * Return the global index of the local - * index with number @p local_index - * stored in this index set. @p - * local_index obviously needs to be less - * than n_elements(). + * Return the global index of the local index with number @p local_index + * stored in this index set. @p local_index obviously needs to be less than + * n_elements(). */ types::global_dof_index nth_index_in_set (const unsigned int local_index) const; /** - * Return the how-manyth element of this - * set (counted in ascending order) @p - * global_index is. @p global_index needs - * to be less than the size(). This - * function throws an exception if the - * index @p global_index is not actually - * a member of this index set, i.e. if - * is_element(global_index) is false. + * Return the how-manyth element of this set (counted in ascending order) @p + * global_index is. @p global_index needs to be less than the size(). This + * function throws an exception if the index @p global_index is not actually + * a member of this index set, i.e. if is_element(global_index) is false. */ types::global_dof_index index_within_set (const types::global_dof_index global_index) const; /** - * Each index set can be - * represented as the union of a - * number of contiguous intervals - * of indices, where if necessary - * intervals may only consist of - * individual elements to - * represent isolated members of - * the index set. + * Each index set can be represented as the union of a number of contiguous + * intervals of indices, where if necessary intervals may only consist of + * individual elements to represent isolated members of the index set. * - * This function returns the - * minimal number of such - * intervals that are needed to - * represent the index set under - * consideration. + * This function returns the minimal number of such intervals that are + * needed to represent the index set under consideration. */ unsigned int n_intervals () const; /** - * Compress the internal - * representation by merging - * individual elements with - * contiguous ranges, etc. This - * function does not have any - * external effect. + * Compress the internal representation by merging individual elements with + * contiguous ranges, etc. This function does not have any external effect. */ void compress () const; /** - * Comparison for equality of - * index sets. This operation is - * only allowed if the size of - * the two sets is the same - * (though of course they do not - * have to have the same number - * of indices). + * Comparison for equality of index sets. This operation is only allowed if + * the size of the two sets is the same (though of course they do not have + * to have the same number of indices). */ bool operator == (const IndexSet &is) const; /** - * Comparison for inequality of - * index sets. This operation is - * only allowed if the size of - * the two sets is the same - * (though of course they do not - * have to have the same number - * of indices). + * Comparison for inequality of index sets. This operation is only allowed + * if the size of the two sets is the same (though of course they do not + * have to have the same number of indices). */ bool operator != (const IndexSet &is) const; /** - * Return the intersection of the - * current index set and the - * argument given, i.e. a set of - * indices that are elements of - * both index sets. The two index - * sets must have the same size - * (though of course they do not - * have to have the same number - * of indices). + * Return the intersection of the current index set and the argument given, + * i.e. a set of indices that are elements of both index sets. The two index + * sets must have the same size (though of course they do not have to have + * the same number of indices). */ IndexSet operator & (const IndexSet &is) const; /** - * This command takes an interval - * [begin, end) and returns - * the intersection of the current - * index set with the interval, shifted - * to the range [0, - * end-begin). + * This command takes an interval [begin, end) and returns the + * intersection of the current index set with the interval, shifted to the + * range [0, end-begin). * - * In other words, the result of this operation is the - * intersection of the set represented by the current object - * and the interval [begin, end), as seen - * within the interval [begin, end) by - * shifting the result of the intersection operation to - * the left by begin. This corresponds to the notion - * of a view: The interval [begin, end) is - * a window through which we see the set represented - * by the current object. + * In other words, the result of this operation is the intersection of the + * set represented by the current object and the interval [begin, + * end), as seen within the interval [begin, end) by + * shifting the result of the intersection operation to the left by + * begin. This corresponds to the notion of a view: The + * interval [begin, end) is a window through which we see + * the set represented by the current object. */ IndexSet get_view (const types::global_dof_index begin, const types::global_dof_index end) const; /** - * Removes all elements contained in @p - * other from this set. In other words, - * if $x$ is the current object and $o$ - * the argument, then we compute $x + * Removes all elements contained in @p other from this set. In other words, + * if $x$ is the current object and $o$ the argument, then we compute $x * \leftarrow x \backslash o$. */ void subtract_set (const IndexSet &other); /** - * Fills the given vector with all - * indices contained in this IndexSet. + * Fills the given vector with all indices contained in this IndexSet. */ void fill_index_vector(std::vector &indices) const; /** - * Fill the given vector with either - * zero or one elements, providing - * a binary representation of this - * index set. The given vector is - * assumed to already have the correct - * size. + * Fill the given vector with either zero or one elements, providing a + * binary representation of this index set. The given vector is assumed to + * already have the correct size. * - * The given argument is filled with - * integer values zero and one, using - * vector.operator[]. Thus, - * any object that has such an operator - * can be used as long as it allows - * conversion of integers zero and one to - * elements of the vector. Specifically, - * this is the case for classes Vector, - * BlockVector, but also - * std::vector@, - * std::vector@, and - * std::vector@. + * The given argument is filled with integer values zero and one, using + * vector.operator[]. Thus, any object that has such an + * operator can be used as long as it allows conversion of integers zero and + * one to elements of the vector. Specifically, this is the case for classes + * Vector, BlockVector, but also std::vector@, std::vector@, + * and std::vector@. */ template void fill_binary_vector (Vector &vector) const; /** - * Outputs a text representation of this - * IndexSet to the given stream. Used for - * testing. + * Outputs a text representation of this IndexSet to the given stream. Used + * for testing. */ template void print(STREAM &out) const; /** - * Writes the IndexSet into a text based - * file format, that can be read in again - * using the read() function. + * Writes the IndexSet into a text based file format, that can be read in + * again using the read() function. */ void write(std::ostream &out) const; /** - * Constructs the IndexSet from a text - * based representation given by the - * stream @param in written by the - * write() function. + * Constructs the IndexSet from a text based representation given by the + * stream @param in written by the write() function. */ void read(std::istream &in); /** - * Writes the IndexSet into a binary, - * compact representation, that can be - * read in again using the block_read() - * function. + * Writes the IndexSet into a binary, compact representation, that can be + * read in again using the block_read() function. */ void block_write(std::ostream &out) const; /** - * Constructs the IndexSet from a binary - * representation given by the stream - * @param in written by the write_block() - * function. + * Constructs the IndexSet from a binary representation given by the stream + * @param in written by the write_block() function. */ void block_read(std::istream &in); #ifdef DEAL_II_WITH_TRILINOS /** - * Given an MPI communicator, - * create a Trilinos map object - * that represents a distribution - * of vector elements or matrix - * rows in which we will locally - * store those elements or rows - * for which we store the index - * in the current index set, and - * all the other elements/rows - * elsewhere on one of the other + * Given an MPI communicator, create a Trilinos map object that represents a + * distribution of vector elements or matrix rows in which we will locally + * store those elements or rows for which we store the index in the current + * index set, and all the other elements/rows elsewhere on one of the other * MPI processes. * - * The last argument only plays a - * role if the communicator is a - * parallel one, distributing - * computations across multiple - * processors. In that case, if - * the last argument is false, - * then it is assumed that the - * index sets this function is - * called on on all processors - * are mutually exclusive but - * together enumerate each index - * exactly once. In other words, - * if you call this function on - * two processors, then the index - * sets this function is called - * with must together have all - * possible indices from zero to - * size()-1, and no index must - * appear in both index - * sets. This corresponds, for - * example, to the case where we - * want to split the elements of - * vectors into unique subsets to - * be stored on different - * processors -- no element - * should be owned by more than - * one processor, but each - * element must be owned by one. + * The last argument only plays a role if the communicator is a parallel + * one, distributing computations across multiple processors. In that case, + * if the last argument is false, then it is assumed that the index sets + * this function is called on on all processors are mutually exclusive but + * together enumerate each index exactly once. In other words, if you call + * this function on two processors, then the index sets this function is + * called with must together have all possible indices from zero to + * size()-1, and no index must appear in both index sets. This corresponds, + * for example, to the case where we want to split the elements of vectors + * into unique subsets to be stored on different processors -- no element + * should be owned by more than one processor, but each element must be + * owned by one. * - * On the other hand, if the - * second argument is true, then - * the index sets can be - * overlapping, though they still - * need to contain each index - * exactly once on all processors - * taken together. This is a - * useful operation if we want to - * create vectors that not only - * contain the locally owned - * indices, but for example also - * the elements that correspond - * to degrees of freedom located - * on ghost cells. + * On the other hand, if the second argument is true, then the index sets + * can be overlapping, though they still need to contain each index exactly + * once on all processors taken together. This is a useful operation if we + * want to create vectors that not only contain the locally owned indices, + * but for example also the elements that correspond to degrees of freedom + * located on ghost cells. */ Epetra_Map make_trilinos_map (const MPI_Comm &communicator = MPI_COMM_WORLD, const bool overlapping = false) const; @@ -423,8 +310,7 @@ public: /** - * Determine an estimate for the memory - * consumption (in bytes) of this + * Determine an estimate for the memory consumption (in bytes) of this * object. */ std::size_t memory_consumption () const; @@ -442,18 +328,11 @@ public: private: /** - * A type that denotes the half - * open index range - * [begin,end). + * A type that denotes the half open index range [begin,end). * - * The nth_index_in_set denotes - * the how many-th index within - * this IndexSet the first - * element of the current range - * is. This information is only - * accurate if - * IndexSet::compress() has been - * called after the last + * The nth_index_in_set denotes the how many-th index within this IndexSet + * the first element of the current range is. This information is only + * accurate if IndexSet::compress() has been called after the last * insertion. */ struct Range @@ -521,72 +400,52 @@ private: } /** - * Write or read the data of this object to or - * from a stream for the purpose of serialization + * Write or read the data of this object to or from a stream for the + * purpose of serialization */ template void serialize (Archive &ar, const unsigned int version); }; /** - * A set of contiguous ranges of - * indices that make up (part of) - * this index set. This variable - * is always kept sorted. + * A set of contiguous ranges of indices that make up (part of) this index + * set. This variable is always kept sorted. * - * The variable is marked - * "mutable" so that it can be - * changed by compress(), though - * this of course doesn't change - * anything about the external - * representation of this index - * set. + * The variable is marked "mutable" so that it can be changed by compress(), + * though this of course doesn't change anything about the external + * representation of this index set. */ mutable std::vector ranges; /** - * True if compress() has been - * called after the last change - * in the set of indices. + * True if compress() has been called after the last change in the set of + * indices. * - * The variable is marked - * "mutable" so that it can be - * changed by compress(), though - * this of course doesn't change - * anything about the external - * representation of this index - * set. + * The variable is marked "mutable" so that it can be changed by compress(), + * though this of course doesn't change anything about the external + * representation of this index set. */ mutable bool is_compressed; /** - * The overall size of the index - * range. Elements of this index - * set have to have a smaller - * number than this value. + * The overall size of the index range. Elements of this index set have to + * have a smaller number than this value. */ types::global_dof_index index_space_size; /** - * This integer caches the index of the - * largest range in @p ranges. This gives - * O(1) access to the range with - * most elements, while general access - * costs O(log(n_ranges)). The - * largest range is needed for the - * methods @p is_element(), @p - * index_within_set(), @p - * nth_index_in_set. In many - * applications, the largest range - * contains most elements (the locally - * owned range), whereas there are only a - * few other elements (ghosts). + * This integer caches the index of the largest range in @p ranges. This + * gives O(1) access to the range with most elements, while general + * access costs O(log(n_ranges)). The largest range is needed for + * the methods @p is_element(), @p index_within_set(), @p + * nth_index_in_set. In many applications, the largest range contains most + * elements (the locally owned range), whereas there are only a few other + * elements (ghosts). */ mutable types::global_dof_index largest_range; /** - * Actually perform the compress() - * operation. + * Actually perform the compress() operation. */ void do_compress() const; }; @@ -701,40 +560,6 @@ IndexSet::compress () const -inline -void -IndexSet::add_range (const types::global_dof_index begin, - const types::global_dof_index end) -{ - Assert ((begin < index_space_size) - || - ((begin == index_space_size) && (end == index_space_size)), - ExcIndexRangeType (begin, 0, index_space_size)); - Assert (end <= index_space_size, - ExcIndexRangeType (end, 0, index_space_size+1)); - Assert (begin <= end, - ExcIndexRangeType (begin, 0, end)); - - if (begin != end) - { - const Range new_range(begin,end); - - // the new index might be larger than the last - // index present in the ranges. Then we can - // skip the binary search - if (ranges.size() == 0 || begin > ranges.back().end) - ranges.push_back(new_range); - else - ranges.insert (Utilities::lower_bound (ranges.begin(), - ranges.end(), - new_range), - new_range); - is_compressed = false; - } -} - - - inline void IndexSet::add_index (const types::global_dof_index index) @@ -763,10 +588,8 @@ void IndexSet::add_indices (const ForwardIterator &begin, const ForwardIterator &end) { - // insert each element of the - // range. if some of them happen to - // be consecutive, merge them to a - // range + // insert each element of the range. if some of them happen to be + // consecutive, merge them to a range for (ForwardIterator p=begin; p!=end;) { const types::global_dof_index begin_index = *p; @@ -786,26 +609,6 @@ IndexSet::add_indices (const ForwardIterator &begin, -inline -void -IndexSet::add_indices(const IndexSet &other, - const unsigned int offset) -{ - if ((this == &other) && (offset == 0)) - return; - - for (std::vector::iterator range = other.ranges.begin(); - range != other.ranges.end(); - ++range) - { - add_range(range->begin+offset, range->end+offset); - } - - compress(); -} - - - inline bool IndexSet::is_element (const types::global_dof_index index) const @@ -814,38 +617,24 @@ IndexSet::is_element (const types::global_dof_index index) const { compress (); - // fast check whether the index is in the - // largest range + // fast check whether the index is in the largest range Assert (largest_range < ranges.size(), ExcInternalError()); if (index >= ranges[largest_range].begin && index < ranges[largest_range].end) return true; - // get the element after which - // we would have to insert a - // range that consists of all - // elements from this element - // to the end of the index - // range plus one. after this - // call we know that if - // p!=end() then - // p->begin<=index unless there - // is no such range at all + // get the element after which we would have to insert a range that + // consists of all elements from this element to the end of the index + // range plus one. after this call we know that if p!=end() then + // p->begin<=index unless there is no such range at all // - // if the searched for element - // is an element of this range, - // then we're done. otherwise, - // the element can't be in one - // of the following ranges - // because otherwise p would be - // a different iterator + // if the searched for element is an element of this range, then we're + // done. otherwise, the element can't be in one of the following ranges + // because otherwise p would be a different iterator // - // since we already know the position - // relative to the largest range (we - // called compress!), we can perform - // the binary search on ranges with - // lower/higher number compared to the - // largest range + // since we already know the position relative to the largest range (we + // called compress!), we can perform the binary search on ranges with + // lower/higher number compared to the largest range std::vector::const_iterator p = std::upper_bound (ranges.begin() + (indexbegin > index), ExcInternalError()); - // now move to that previous - // range + // now move to that previous range --p; Assert (p->begin <= index, ExcInternalError()); return (p->end > index); } - // didn't find this index, so it's - // not in the set + // didn't find this index, so it's not in the set return false; } @@ -889,8 +676,7 @@ inline types::global_dof_index IndexSet::n_elements () const { - // make sure we have - // non-overlapping ranges + // make sure we have non-overlapping ranges compress (); types::global_dof_index v = 0; @@ -914,28 +700,35 @@ IndexSet::n_elements () const +inline +unsigned int +IndexSet::n_intervals () const +{ + compress (); + return ranges.size(); +} + + + inline types::global_dof_index IndexSet::nth_index_in_set (const unsigned int n) const { - // to make this call thread-safe, compress() - // must not be called through this function + // to make this call thread-safe, compress() must not be called through this + // function Assert (is_compressed == true, ExcMessage ("IndexSet must be compressed.")); Assert (n < n_elements(), ExcIndexRangeType (n, 0, n_elements())); - // first check whether the index is in the - // largest range + // first check whether the index is in the largest range Assert (largest_range < ranges.size(), ExcInternalError()); std::vector::const_iterator main_range=ranges.begin()+largest_range; if (n>=main_range->nth_index_in_set && nnth_index_in_set+(main_range->end-main_range->begin)) return main_range->begin + (n-main_range->nth_index_in_set); - // find out which chunk the local index n - // belongs to by using a binary search. the - // comparator is based on the end of the - // ranges. Use the position relative to main_range to - // subdivide the ranges + // find out which chunk the local index n belongs to by using a binary + // search. the comparator is based on the end of the ranges. Use the + // position relative to main_range to subdivide the ranges Range r (n,n+1); r.nth_index_in_set = n; std::vector::const_iterator range_begin, range_end; @@ -969,15 +762,14 @@ inline types::global_dof_index IndexSet::index_within_set (const types::global_dof_index n) const { - // to make this call thread-safe, compress() - // must not be called through this function + // to make this call thread-safe, compress() must not be called through this + // function Assert (is_compressed == true, ExcMessage ("IndexSet must be compressed.")); Assert (is_element(n) == true, ExcIndexNotPresent (n)); Assert (n < size(), ExcIndexRangeType (n, 0, size())); - // check whether the index is in the largest - // range. use the result to perform a - // one-sided binary search afterward + // check whether the index is in the largest range. use the result to + // perform a one-sided binary search afterward Assert (largest_range < ranges.size(), ExcInternalError()); std::vector::const_iterator main_range=ranges.begin()+largest_range; if (n >= main_range->begin && n < main_range->end) @@ -1046,12 +838,11 @@ IndexSet::fill_binary_vector (Vector &vector) const ExcDimensionMismatch (vector.size(), size())); compress(); - // first fill all elements of the vector - // with zeroes. + // first fill all elements of the vector with zeroes. std::fill (vector.begin(), vector.end(), 0); - // then write ones into the elements whose - // indices are contained in the index set + // then write ones into the elements whose indices are contained in the + // index set for (std::vector::iterator it = ranges.begin(); it != ranges.end(); ++it) diff --git a/deal.II/source/base/index_set.cc b/deal.II/source/base/index_set.cc index ce4b5bcfb7..17d1c555b7 100644 --- a/deal.II/source/base/index_set.cc +++ b/deal.II/source/base/index_set.cc @@ -28,16 +28,48 @@ DEAL_II_NAMESPACE_OPEN + +void +IndexSet::add_range (const types::global_dof_index begin, + const types::global_dof_index end) +{ + Assert ((begin < index_space_size) + || + ((begin == index_space_size) && (end == index_space_size)), + ExcIndexRangeType (begin, 0, index_space_size)); + Assert (end <= index_space_size, + ExcIndexRangeType (end, 0, index_space_size+1)); + Assert (begin <= end, + ExcIndexRangeType (begin, 0, end)); + + if (begin != end) + { + const Range new_range(begin,end); + + // the new index might be larger than the last index present in the + // ranges. Then we can skip the binary search + if (ranges.size() == 0 || begin > ranges.back().end) + ranges.push_back(new_range); + else + ranges.insert (Utilities::lower_bound (ranges.begin(), + ranges.end(), + new_range), + new_range); + is_compressed = false; + } +} + + + void IndexSet::do_compress () const { - // see if any of the - // contiguous ranges can be - // merged. since they are sorted by - // their first index, determining + // see if any of the contiguous ranges can be merged. do not use + // std::vector::erase in-place as it is quadratic in the number of + // ranges. since the ranges are sorted by their first index, determining // overlap isn't all that hard - for (std::vector::iterator - i = ranges.begin(); + std::vector::iterator store = ranges.begin(); + for (std::vector::iterator i = ranges.begin(); i != ranges.end(); ) { std::vector::iterator @@ -47,37 +79,29 @@ IndexSet::do_compress () const types::global_dof_index first_index = i->begin; types::global_dof_index last_index = i->end; - // see if we can merge any of - // the following ranges - bool can_merge = false; + // see if we can merge any of the following ranges while (next != ranges.end() && (next->begin <= last_index)) { last_index = std::max (last_index, next->end); ++next; - can_merge = true; } + i = next; - if (can_merge == true) - { - // delete the old ranges - // and insert the new range - // in place of the previous - // one - *i = Range(first_index, last_index); - i = ranges.erase (i+1, next); - } - else - ++i; + // store the new range in the slot we last occupied + *store = Range(first_index, last_index); + ++store; + } + // use a compact array with exactly the right amount of storage + if (store != ranges.end()) + { + std::vector new_ranges(ranges.begin(), store); + ranges.swap(new_ranges); } - - // now compute indices within set and the - // range with most elements + // now compute indices within set and the range with most elements types::global_dof_index next_index = 0, largest_range_size = 0; - for (std::vector::iterator - i = ranges.begin(); - i != ranges.end(); + for (std::vector::iterator i = ranges.begin(); i != ranges.end(); ++i) { Assert(i->begin < i->end, ExcInternalError()); @@ -92,11 +116,8 @@ IndexSet::do_compress () const } is_compressed = true; - // check that next_index is - // correct. needs to be after the - // previous statement because we - // otherwise will get into an - // endless loop + // check that next_index is correct. needs to be after the previous + // statement because we otherwise will get into an endless loop Assert (next_index == n_elements(), ExcInternalError()); } @@ -119,18 +140,15 @@ IndexSet::operator & (const IndexSet &is) const && (r2 != is.ranges.end())) { - // if r1 and r2 do not overlap - // at all, then move the - // pointer that sits to the - // left of the other up by one + // if r1 and r2 do not overlap at all, then move the pointer that sits + // to the left of the other up by one if (r1->end <= r2->begin) ++r1; else if (r2->end <= r1->begin) ++r2; else { - // the ranges must overlap - // somehow + // the ranges must overlap somehow Assert (((r1->begin <= r2->begin) && (r1->end > r2->begin)) || @@ -138,21 +156,15 @@ IndexSet::operator & (const IndexSet &is) const (r2->end > r1->begin)), ExcInternalError()); - // add the overlapping - // range to the result + // add the overlapping range to the result result.add_range (std::max (r1->begin, r2->begin), std::min (r1->end, r2->end)); - // now move that iterator - // that ends earlier one - // up. note that it has to - // be this one because a - // subsequent range may - // still have a chance of - // overlapping with the - // range that ends later + // now move that iterator that ends earlier one up. note that it has + // to be this one because a subsequent range may still have a chance + // of overlapping with the range that ends later if (r1->end <= r2->end) ++r1; else @@ -166,15 +178,6 @@ IndexSet::operator & (const IndexSet &is) const -unsigned int -IndexSet::n_intervals () const -{ - compress (); - return ranges.size(); -} - - - IndexSet IndexSet::get_view (const types::global_dof_index begin, const types::global_dof_index end) const @@ -209,6 +212,131 @@ IndexSet::get_view (const types::global_dof_index begin, +void +IndexSet::subtract_set (const IndexSet &other) +{ + compress(); + other.compress(); + is_compressed = false; + + + // we save new ranges to be added to our IndexSet in an temporary list and + // add all of them in one go at the end. This is necessary because a growing + // ranges vector invalidates iterators. + std::list temp_list; + + std::vector::iterator own_it = ranges.begin(); + std::vector::iterator other_it = other.ranges.begin(); + + while (own_it != ranges.end() && other_it != other.ranges.end()) + { + //advance own iterator until we get an overlap + if (own_it->end <= other_it->begin) + { + ++own_it; + continue; + } + //we are done with other_it, so advance + if (own_it->begin >= other_it->end) + { + ++other_it; + continue; + } + + //Now own_it and other_it overlap. First save the part of own_it that + //is before other_it (if not empty). + if (own_it->begin < other_it->begin) + { + Range r(own_it->begin, other_it->begin); + r.nth_index_in_set = 0; //fix warning of unused variable + temp_list.push_back(r); + } + // change own_it to the sub range behind other_it. Do not delete own_it + // in any case. As removal would invalidate iterators, we just shrink + // the range to an empty one. + own_it->begin = other_it->end; + if (own_it->begin > own_it->end) + { + own_it->begin = own_it->end; + ++own_it; + } + + // continue without advancing iterators, the right one will be advanced + // next. + } + + // Now delete all empty ranges we might + // have created. + for (std::vector::iterator it = ranges.begin(); + it != ranges.end(); ) + { + if (it->begin >= it->end) + it = ranges.erase(it); + else + ++it; + } + + // done, now add the temporary ranges + for (std::list::iterator it = temp_list.begin(); + it != temp_list.end(); + ++it) + add_range(it->begin, it->end); + + compress(); +} + + + +void +IndexSet::add_indices(const IndexSet &other, + const unsigned int offset) +{ + if ((this == &other) && (offset == 0)) + return; + + compress(); + other.compress(); + + std::vector::const_iterator r1 = ranges.begin(), + r2 = other.ranges.begin(); + + std::vector new_ranges; + // just get the start and end of the ranges right in this method, everything + // else will be done in compress() + while (r1 != ranges.end() || r2 != other.ranges.end()) + { + // the two ranges do not overlap or we are at the end of one of the + // ranges + if (r2 == other.ranges.end() || + (r1 != ranges.end() && r1->end < (r2->begin+offset))) + { + new_ranges.push_back(*r1); + ++r1; + } + else if (r1 == ranges.end() || (r2->end+offset) < r1->begin) + { + new_ranges.push_back(Range(r2->begin+offset,r2->end+offset)); + ++r2; + } + else + { + // ok, we do overlap, so just take the combination of the current + // range (do not bother to merge with subsequent ranges) + Range next(std::min(r1->begin, r2->begin+offset), + std::max(r1->end, r2->end+offset)); + new_ranges.push_back(next); + ++r1; + ++r2; + } + } + ranges.swap(new_ranges); + + is_compressed = false; + compress(); +} + + + void IndexSet::write(std::ostream &out) const { @@ -277,87 +405,6 @@ IndexSet::block_read(std::istream &in) -void -IndexSet::subtract_set (const IndexSet &other) -{ - compress(); - other.compress(); - is_compressed = false; - - - // we save new ranges to be added to our - // IndexSet in an temporary list and add - // all of them in one go at the end. This - // is necessary because a growing ranges - // vector invalidates iterators. - std::list temp_list; - - std::vector::iterator own_it = ranges.begin(); - std::vector::iterator other_it = other.ranges.begin(); - - while (own_it != ranges.end() && other_it != other.ranges.end()) - { - //advance own iterator until we get an - //overlap - if (own_it->end <= other_it->begin) - { - ++own_it; - continue; - } - //we are done with other_it, so advance - if (own_it->begin >= other_it->end) - { - ++other_it; - continue; - } - - //Now own_it and other_it overlap. - //First save the part of own_it that is - //before other_it (if not empty). - if (own_it->begin < other_it->begin) - { - Range r(own_it->begin, other_it->begin); - r.nth_index_in_set = 0; //fix warning of unused variable - temp_list.push_back(r); - } - // change own_it to the sub range - // behind other_it. Do not delete - // own_it in any case. As removal would - // invalidate iterators, we just shrink - // the range to an empty one. - own_it->begin = other_it->end; - if (own_it->begin > own_it->end) - { - own_it->begin = own_it->end; - ++own_it; - } - - // continue without advancing - // iterators, the right one will be - // advanced next. - } - - // Now delete all empty ranges we might - // have created. - for (std::vector::iterator it = ranges.begin(); - it != ranges.end(); ) - { - if (it->begin >= it->end) - it = ranges.erase(it); - else - ++it; - } - - // done, now add the temporary ranges - for (std::list::iterator it = temp_list.begin(); - it != temp_list.end(); - ++it) - add_range(it->begin, it->end); - - compress(); -} - - void IndexSet::fill_index_vector(std::vector &indices) const { compress(); diff --git a/deal.II/source/multigrid/mg_tools.cc b/deal.II/source/multigrid/mg_tools.cc index 7321c64d2e..65c253e033 100644 --- a/deal.II/source/multigrid/mg_tools.cc +++ b/deal.II/source/multigrid/mg_tools.cc @@ -1581,11 +1581,10 @@ namespace MGTools ExcDimensionMismatch (boundary_interface_dofs.size(), mg_dof_handler.get_tria().n_global_levels())); - for (unsigned int l=0; l > + tmp_interface_dofs(interface_dofs.size()); + std::vector > + tmp_boundary_interface_dofs(interface_dofs.size()); const FiniteElement &fe = mg_dof_handler.get_fe(); @@ -1633,12 +1632,14 @@ namespace MGTools } } - if (has_coarser_neighbor == true) - for (unsigned int face_nr=0; face_nr::faces_per_cell; ++face_nr) - if (cell->at_boundary(face_nr)) - for (unsigned int j=0; j::faces_per_cell; ++face_nr) + if (cell->at_boundary(face_nr)) + for (unsigned int j=0; jlevel(); @@ -1647,12 +1648,31 @@ namespace MGTools for (unsigned int i=0; i