]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix quadratic complexity in number of ranges in two methods. Do not directly insert...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 16 Jan 2014 12:40:23 +0000 (12:40 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 16 Jan 2014 12:40:23 +0000 (12:40 +0000)
git-svn-id: https://svn.dealii.org/trunk@32224 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/base/index_set.h
deal.II/source/base/index_set.cc
deal.II/source/multigrid/mg_tools.cc

index 4c21c7596b20894ca312770d5482e0e5ea32cd27..54e3ae91720864008a30a92793d1e5c4d2b4f3d6 100644 (file)
@@ -85,6 +85,13 @@ inconvenience this causes.
 
 <ol>
 
+  <li> 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.
+  <br>
+  (Martin Kronbichler, 2014/01/15)
+  </li>
+
   <li> 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.
index b5b439d9901d9ffe870737e9c2e8ce2e7be1f7c6..166a642a235894a9cee92231f3fcf4576c51bf7a 100644 (file)
@@ -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
-   * <code>[begin,end)</code>.
+   * Add a whole set of indices described by dereferencing every element of
+   * the the iterator range <code>[begin,end)</code>.
    */
   template <typename ForwardIterator>
   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 <code>[0,size())</code> 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
+   * <code>[0,size())</code> 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
-   * <tt>[begin, end)</tt> and returns
-   * the intersection of the current
-   * index set with the interval, shifted
-   * to the range <tt>[0,
-   * end-begin)</tt>.
+   * This command takes an interval <tt>[begin, end)</tt> and returns the
+   * intersection of the current index set with the interval, shifted to the
+   * range <tt>[0, end-begin)</tt>.
    *
-   * In other words, the result of this operation is the
-   * intersection of the set represented by the current object
-   * and the interval <tt>[begin, end)</tt>, as seen
-   * <i>within the interval <tt>[begin, end)</tt></i> by
-   * shifting the result of the intersection operation to
-   * the left by <tt>begin</tt>. This corresponds to the notion
-   * of a <i>view</i>: The interval <tt>[begin, end)</tt> is
-   * a <i>window</i> 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 <tt>[begin,
+   * end)</tt>, as seen <i>within the interval <tt>[begin, end)</tt></i> by
+   * shifting the result of the intersection operation to the left by
+   * <tt>begin</tt>. This corresponds to the notion of a <i>view</i>: The
+   * interval <tt>[begin, end)</tt> is a <i>window</i> 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<types::global_dof_index> &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
-   * <code>vector.operator[]</code>. 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@<bool@>,
-   * std::vector@<int@>, and
-   * std::vector@<double@>.
+   * The given argument is filled with integer values zero and one, using
+   * <code>vector.operator[]</code>. 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@<bool@>, std::vector@<int@>,
+   * and std::vector@<double@>.
    */
   template <typename Vector>
   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 <class STREAM>
   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
-   * <code>[begin,end)</code>.
+   * A type that denotes the half open index range <code>[begin,end)</code>.
    *
-   * 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 <class Archive>
     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<Range> 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
-   * <tt>O(1)</tt> access to the range with
-   * most elements, while general access
-   * costs <tt>O(log(n_ranges))</tt>. 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 <tt>O(1)</tt> access to the range with most elements, while general
+   * access costs <tt>O(log(n_ranges))</tt>. 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<types::global_dof_index> (begin, 0, index_space_size));
-  Assert (end <= index_space_size,
-          ExcIndexRangeType<types::global_dof_index> (end, 0, index_space_size+1));
-  Assert (begin <= end,
-          ExcIndexRangeType<types::global_dof_index> (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<Range>::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<Range>::const_iterator
       p = std::upper_bound (ranges.begin() + (index<ranges[largest_range].begin?
                                               0 : largest_range+1),
@@ -860,16 +649,14 @@ IndexSet::is_element (const types::global_dof_index index) const
       Assert ((p == ranges.end()) || (p->begin > 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<types::global_dof_index> (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<Range>::const_iterator main_range=ranges.begin()+largest_range;
   if (n>=main_range->nth_index_in_set &&
       n<main_range->nth_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<Range>::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<types::global_dof_index> (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<Range>::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<Range>::iterator it = ranges.begin();
        it != ranges.end();
        ++it)
index ce4b5bcfb716dfd71bc2f57b696fc20017fc77df..17d1c555b74f511da05be32ef9eb0fc108b43e9c 100644 (file)
 
 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<types::global_dof_index> (begin, 0, index_space_size));
+  Assert (end <= index_space_size,
+          ExcIndexRangeType<types::global_dof_index> (end, 0, index_space_size+1));
+  Assert (begin <= end,
+          ExcIndexRangeType<types::global_dof_index> (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<Range>::iterator
-       i = ranges.begin();
+  std::vector<Range>::iterator store = ranges.begin();
+  for (std::vector<Range>::iterator i = ranges.begin();
        i != ranges.end(); )
     {
       std::vector<Range>::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<Range> 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<Range>::iterator
-       i = ranges.begin();
-       i != ranges.end();
+  for (std::vector<Range>::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<Range> temp_list;
+
+  std::vector<Range>::iterator own_it = ranges.begin();
+  std::vector<Range>::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<Range>::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<Range>::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<Range>::const_iterator r1 = ranges.begin(),
+                                     r2 = other.ranges.begin();
+
+  std::vector<Range> 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<Range> temp_list;
-
-  std::vector<Range>::iterator own_it = ranges.begin();
-  std::vector<Range>::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<Range>::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<Range>::iterator it = temp_list.begin();
-       it != temp_list.end();
-       ++it)
-    add_range(it->begin, it->end);
-
-  compress();
-}
-
-
 void IndexSet::fill_index_vector(std::vector<types::global_dof_index> &indices) const
 {
   compress();
index 7321c64d2e4f6fb68ad3afc84a75cc45fc40904a..65c253e033b0598a97afc77f8055b13fd26f55bf 100644 (file)
@@ -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<mg_dof_handler.get_tria().n_global_levels(); ++l)
-      {
-        interface_dofs[l].clear();
-        boundary_interface_dofs[l].clear();
-      }
+    std::vector<std::vector<types::global_dof_index> >
+      tmp_interface_dofs(interface_dofs.size());
+    std::vector<std::vector<types::global_dof_index> >
+      tmp_boundary_interface_dofs(interface_dofs.size());
 
     const FiniteElement<dim,spacedim> &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<GeometryInfo<dim>::faces_per_cell; ++face_nr)
-            if (cell->at_boundary(face_nr))
-              for (unsigned int j=0; j<dofs_per_face; ++j)
+        if (has_coarser_neighbor == false)
+          continue;
+
+        for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
+          if (cell->at_boundary(face_nr))
+            for (unsigned int j=0; j<dofs_per_face; ++j)
 //            if (cell_dofs[fe.face_to_cell_index(j,face_nr)] == true) //is this necessary?
-                boundary_cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
+              boundary_cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
 
 
         const unsigned int level = cell->level();
@@ -1647,12 +1648,31 @@ namespace MGTools
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           {
             if (cell_dofs[i])
-              interface_dofs[level].add_index(local_dof_indices[i]);
+              tmp_interface_dofs[level].push_back(local_dof_indices[i]);
 
             if (boundary_cell_dofs[i])
-              boundary_interface_dofs[level].add_index(local_dof_indices[i]);
+              tmp_boundary_interface_dofs[level].push_back(local_dof_indices[i]);
           }
       }
+
+
+    for (unsigned int l=0; l<mg_dof_handler.get_tria().n_global_levels(); ++l)
+      {
+        interface_dofs[l].clear();
+        std::sort(tmp_interface_dofs[l].begin(), tmp_interface_dofs[l].end());
+        interface_dofs[l].add_indices(tmp_interface_dofs[l].begin(),
+                                      std::unique(tmp_interface_dofs[l].begin(),
+                                                  tmp_interface_dofs[l].end()));
+        interface_dofs[l].compress();
+        boundary_interface_dofs[l].clear();
+        std::sort(tmp_boundary_interface_dofs[l].begin(),
+                  tmp_boundary_interface_dofs[l].end());
+        boundary_interface_dofs[l].add_indices(tmp_boundary_interface_dofs[l].begin(),
+                                               std::unique(tmp_boundary_interface_dofs[l].begin(),
+                                                           tmp_boundary_interface_dofs[l].end()));
+        boundary_interface_dofs[l].compress();
+      }
+
   }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.