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;
/**
- * 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;
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
}
/**
- * 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;
};
-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)
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;
-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
{
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),
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;
}
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;
+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;
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)
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)
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
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());
}
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());
}
&&
(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))
||
(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
-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
+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
{
-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();
projects / dealii.git / commitdiff