index_set.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2009 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_index_set_h
#define dealii_index_set_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/memory_space.h>
#include <deal.II/base/mpi_stub.h>
#include <deal.II/base/mutex.h>
#include <deal.II/base/trilinos_utilities.h>
 
#include <deal.II/lac/trilinos_tpetra_types.h>
 
#include <boost/container/small_vector.hpp>
 
#include <algorithm>
#include <vector>
 
 
#ifdef DEAL_II_WITH_TRILINOS
#  include <Epetra_Map.h>
#  ifdef DEAL_II_TRILINOS_WITH_TPETRA
#    include <Tpetra_Map.hpp>
#  endif
#endif
 
#ifdef DEAL_II_WITH_PETSC
#  include <petscis.h>
#endif
 
DEAL_II_NAMESPACE_OPEN
 
class IndexSet
{
public:
  // forward declarations:
  class ElementIterator;
  class IntervalIterator;
 
  using size_type = types::global_dof_index;
 
  using value_type = signed int;
 
 
  IndexSet();
 
  explicit IndexSet(const size_type size);
 
  IndexSet(const IndexSet &) = default;
 
  IndexSet &
  operator=(const IndexSet &) = default;
 
  IndexSet(IndexSet &&is) noexcept;
 
  IndexSet &
  operator=(IndexSet &&is) noexcept;
 
#ifdef DEAL_II_WITH_TRILINOS
 
#  ifdef DEAL_II_TRILINOS_WITH_TPETRA
  template <typename NodeType>
  explicit IndexSet(
    const Teuchos::RCP<
      const Tpetra::Map<int, types::signed_global_dof_index, NodeType>> &map);
#  endif // DEAL_II_TRILINOS_WITH_TPETRA
 
  explicit IndexSet(const Epetra_BlockMap &map);
#endif // DEAL_II_WITH_TRILINOS
 
  void
  clear();
 
  void
  set_size(const size_type size);
 
  size_type
  size() const;
 
  void
  add_range(const size_type begin, const size_type end);
 
  void
  add_index(const size_type index);
 
  template <typename ForwardIterator>
  void
  add_indices(const ForwardIterator &begin, const ForwardIterator &end);
 
  void
  add_indices(const IndexSet &other, const size_type offset = 0);
 
  bool
  is_element(const size_type index) const;
 
  bool
  is_contiguous() const;
 
  bool
  is_empty() const;
 
  bool
  is_ascending_and_one_to_one(const MPI_Comm communicator) const;
 
  size_type
  n_elements() const;
 
  size_type
  nth_index_in_set(const size_type local_index) const;
 
  size_type
  index_within_set(const size_type global_index) const;
 
  unsigned int
  n_intervals() const;
 
  size_type
  largest_range_starting_index() const;
 
  void
  compress() const;
 
  bool
  operator==(const IndexSet &is) const;
 
  bool
  operator!=(const IndexSet &is) const;
 
  IndexSet
  operator&(const IndexSet &is) const;
 
  IndexSet
  get_view(const size_type begin, const size_type end) const;
 
  IndexSet
  get_view(const IndexSet &mask) const;
 
  std::vector<IndexSet>
  split_by_block(
    const std::vector<types::global_dof_index> &n_indices_per_block) const;
 
  void
  subtract_set(const IndexSet &other);
 
  IndexSet
  tensor_product(const IndexSet &other) const;
 
  size_type
  pop_back();
 
  size_type
  pop_front();
 
  std::vector<size_type>
  get_index_vector() const;
 
  DEAL_II_DEPRECATED
  void
  fill_index_vector(std::vector<size_type> &indices) const;
 
  template <typename VectorType>
  void
  fill_binary_vector(VectorType &vector) const;
 
  bool
  is_subset_of(const IndexSet &other) const;
 
  template <typename StreamType>
  void
  print(StreamType &out) const;
 
  void
  write(std::ostream &out) const;
 
  void
  read(std::istream &in);
 
  void
  block_write(std::ostream &out) const;
 
  void
  block_read(std::istream &in);
 
#ifdef DEAL_II_WITH_TRILINOS
  Epetra_Map
  make_trilinos_map(const MPI_Comm communicator = MPI_COMM_WORLD,
                    const bool     overlapping  = false) const;
 
#  ifdef DEAL_II_TRILINOS_WITH_TPETRA
  template <
    typename NodeType =
      LinearAlgebra::TpetraWrappers::TpetraTypes::NodeType<MemorySpace::Host>>
  Tpetra::Map<int, types::signed_global_dof_index, NodeType>
  make_tpetra_map(const MPI_Comm communicator = MPI_COMM_WORLD,
                  const bool     overlapping  = false) const;
 
  template <
    typename NodeType =
      LinearAlgebra::TpetraWrappers::TpetraTypes::NodeType<MemorySpace::Host>>
  Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index, NodeType>>
  make_tpetra_map_rcp(const MPI_Comm communicator = MPI_COMM_WORLD,
                      const bool     overlapping  = false) const;
#  endif
#endif
 
#ifdef DEAL_II_WITH_PETSC
  IS
  make_petsc_is(const MPI_Comm communicator = MPI_COMM_WORLD) const;
#endif
 
 
  std::size_t
  memory_consumption() const;
 
  DeclException1(ExcIndexNotPresent,
                 size_type,
                 << "The global index " << arg1
                 << " is not an element of this set.");
 
  template <class Archive>
  void
  serialize(Archive &ar, const unsigned int version);
 
 
  class IntervalAccessor
  {
  public:
    IntervalAccessor(const IndexSet *idxset, const size_type range_idx);
 
    explicit IntervalAccessor(const IndexSet *idxset);
 
    size_type
    n_elements() const;
 
    bool
    is_valid() const;
 
    ElementIterator
    begin() const;
 
    ElementIterator
    end() const;
 
    size_type
    last() const;
 
  private:
    IntervalAccessor(const IntervalAccessor &other);
    IntervalAccessor &
    operator=(const IntervalAccessor &other);
 
    bool
    operator==(const IntervalAccessor &other) const;
    bool
    operator<(const IntervalAccessor &other) const;
    void
    advance();
    const IndexSet *index_set;
 
    size_type range_idx;
 
    friend class IntervalIterator;
  };
 
  class IntervalIterator
  {
  public:
    IntervalIterator(const IndexSet *idxset, const size_type range_idx);
 
    explicit IntervalIterator(const IndexSet *idxset);
 
    IntervalIterator();
 
    IntervalIterator(const IntervalIterator &other) = default;
 
    IntervalIterator &
    operator=(const IntervalIterator &other) = default;
 
    IntervalIterator &
    operator++();
 
    IntervalIterator
    operator++(int);
 
    const IntervalAccessor &
    operator*() const;
 
    const IntervalAccessor *
    operator->() const;
 
    bool
    operator==(const IntervalIterator &) const;
 
    bool
    operator!=(const IntervalIterator &) const;
 
    bool
    operator<(const IntervalIterator &) const;
 
    int
    operator-(const IntervalIterator &p) const;
 
    using iterator_category = std::forward_iterator_tag;
    using value_type        = IntervalAccessor;
    using difference_type   = std::ptrdiff_t;
    using pointer           = IntervalAccessor *;
    using reference         = IntervalAccessor &;
 
  private:
    IntervalAccessor accessor;
  };
 
  class ElementIterator
  {
  public:
    ElementIterator(const IndexSet *idxset,
                    const size_type range_idx,
                    const size_type index);
 
    explicit ElementIterator(const IndexSet *idxset);
 
    size_type
    operator*() const;
 
    bool
    is_valid() const;
 
    ElementIterator &
    operator++();
 
    ElementIterator
    operator++(int);
 
    bool
    operator==(const ElementIterator &) const;
 
    bool
    operator!=(const ElementIterator &) const;
 
    bool
    operator<(const ElementIterator &) const;
 
    std::ptrdiff_t
    operator-(const ElementIterator &p) const;
 
    using iterator_category = std::forward_iterator_tag;
    using value_type        = size_type;
    using difference_type   = std::ptrdiff_t;
    using pointer           = size_type *;
    using reference         = size_type &;
 
  private:
    void
    advance();
 
    const IndexSet *index_set;
    size_type range_idx;
    size_type idx;
  };
 
  ElementIterator
  begin() const;
 
  ElementIterator
  at(const size_type global_index) const;
 
  ElementIterator
  end() const;
 
  IntervalIterator
  begin_intervals() const;
 
  IntervalIterator
  end_intervals() const;
 
private:
  struct Range
  {
    size_type begin;
    size_type end;
 
    size_type nth_index_in_set;
 
    Range();
 
    Range(const size_type i1, const size_type i2);
 
    friend inline bool
    operator<(const Range &range_1, const Range &range_2)
    {
      return (
        (range_1.begin < range_2.begin) ||
        ((range_1.begin == range_2.begin) && (range_1.end < range_2.end)));
    }
 
    static bool
    end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
    {
      return x.end < y.end;
    }
 
    static bool
    nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
    {
      return (x.nth_index_in_set + (x.end - x.begin) <
              y.nth_index_in_set + (y.end - y.begin));
    }
 
    friend inline bool
    operator==(const Range &range_1, const Range &range_2)
    {
      return ((range_1.begin == range_2.begin) && (range_1.end == range_2.end));
    }
 
    static std::size_t
    memory_consumption()
    {
      return sizeof(Range);
    }
 
    template <class Archive>
    void
    serialize(Archive &ar, const unsigned int version);
  };
 
  mutable std::vector<Range> ranges;
 
  mutable bool is_compressed;
 
  size_type index_space_size;
 
  mutable size_type largest_range;
 
  mutable Threads::Mutex compress_mutex;
 
  void
  do_compress() const;
 
  bool
  is_element_binary_search(const size_type local_index) const;
 
  size_type
  nth_index_in_set_binary_search(const size_type local_index) const;
 
  size_type
  index_within_set_binary_search(const size_type global_index) const;
 
  void
  add_range_lower_bound(const Range &range);
 
  void
  add_ranges_internal(
    boost::container::small_vector<std::pair<size_type, size_type>, 200>
              &tmp_ranges,
    const bool ranges_are_sorted);
};
 
 
 
inline IndexSet
complete_index_set(const IndexSet::size_type N)
{
  IndexSet is(N);
  is.add_range(0, N);
  is.compress();
  return is;
}
 
/* ------------------ inline functions ------------------ */
 
 
/* IntervalAccessor */
 
inline IndexSet::IntervalAccessor::IntervalAccessor(
  const IndexSet           *idxset,
  const IndexSet::size_type range_idx)
  : index_set(idxset)
  , range_idx(range_idx)
{
  Assert(range_idx < idxset->n_intervals(),
         ExcInternalError("Invalid range index"));
}
 
 
 
inline IndexSet::IntervalAccessor::IntervalAccessor(const IndexSet *idxset)
  : index_set(idxset)
  , range_idx(numbers::invalid_dof_index)
{}
 
 
 
inline IndexSet::IntervalAccessor::IntervalAccessor(
  const IndexSet::IntervalAccessor &other)
  : index_set(other.index_set)
  , range_idx(other.range_idx)
{
  Assert(range_idx == numbers::invalid_dof_index || is_valid(),
         ExcMessage("invalid iterator"));
}
 
 
 
inline IndexSet::size_type
IndexSet::IntervalAccessor::n_elements() const
{
  Assert(is_valid(), ExcMessage("invalid iterator"));
  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
}
 
 
 
inline bool
IndexSet::IntervalAccessor::is_valid() const
{
  return index_set != nullptr && range_idx < index_set->n_intervals();
}
 
 
 
inline IndexSet::ElementIterator
IndexSet::IntervalAccessor::begin() const
{
  Assert(is_valid(), ExcMessage("invalid iterator"));
  return {index_set, range_idx, index_set->ranges[range_idx].begin};
}
 
 
 
inline IndexSet::ElementIterator
IndexSet::IntervalAccessor::end() const
{
  Assert(is_valid(), ExcMessage("invalid iterator"));
 
  // point to first index in next interval unless we are the last interval.
  if (range_idx < index_set->ranges.size() - 1)
    return {index_set, range_idx + 1, index_set->ranges[range_idx + 1].begin};
  else
    return index_set->end();
}
 
 
 
inline IndexSet::size_type
IndexSet::IntervalAccessor::last() const
{
  Assert(is_valid(), ExcMessage("invalid iterator"));
 
  return index_set->ranges[range_idx].end - 1;
}
 
 
 
inline IndexSet::IntervalAccessor &
IndexSet::IntervalAccessor::operator=(const IndexSet::IntervalAccessor &other)
{
  index_set = other.index_set;
  range_idx = other.range_idx;
  Assert(range_idx == numbers::invalid_dof_index || is_valid(),
         ExcMessage("invalid iterator"));
  return *this;
}
 
 
 
inline bool
IndexSet::IntervalAccessor::operator==(
  const IndexSet::IntervalAccessor &other) const
{
  Assert(index_set == other.index_set,
         ExcMessage(
           "Can not compare accessors pointing to different IndexSets"));
  return range_idx == other.range_idx;
}
 
 
 
inline bool
IndexSet::IntervalAccessor::operator<(
  const IndexSet::IntervalAccessor &other) const
{
  Assert(index_set == other.index_set,
         ExcMessage(
           "Can not compare accessors pointing to different IndexSets"));
  return range_idx < other.range_idx;
}
 
 
 
inline void
IndexSet::IntervalAccessor::advance()
{
  Assert(
    is_valid(),
    ExcMessage(
      "Impossible to advance an IndexSet::IntervalIterator that is invalid"));
  ++range_idx;
 
  // set ourselves to invalid if we walk off the end
  if (range_idx >= index_set->ranges.size())
    range_idx = numbers::invalid_dof_index;
}
 
 
/* IntervalIterator */
 
inline IndexSet::IntervalIterator::IntervalIterator(
  const IndexSet           *idxset,
  const IndexSet::size_type range_idx)
  : accessor(idxset, range_idx)
{}
 
 
 
inline IndexSet::IntervalIterator::IntervalIterator()
  : accessor(nullptr)
{}
 
 
 
inline IndexSet::IntervalIterator::IntervalIterator(const IndexSet *idxset)
  : accessor(idxset)
{}
 
 
 
inline IndexSet::IntervalIterator &
IndexSet::IntervalIterator::operator++()
{
  accessor.advance();
  return *this;
}
 
 
 
inline IndexSet::IntervalIterator
IndexSet::IntervalIterator::operator++(int)
{
  const IndexSet::IntervalIterator iter = *this;
  accessor.advance();
  return iter;
}
 
 
 
inline const IndexSet::IntervalAccessor &
IndexSet::IntervalIterator::operator*() const
{
  return accessor;
}
 
 
 
inline const IndexSet::IntervalAccessor *
IndexSet::IntervalIterator::operator->() const
{
  return &accessor;
}
 
 
 
inline bool
IndexSet::IntervalIterator::operator==(
  const IndexSet::IntervalIterator &other) const
{
  return accessor == other.accessor;
}
 
 
 
inline bool
IndexSet::IntervalIterator::operator!=(
  const IndexSet::IntervalIterator &other) const
{
  return !(*this == other);
}
 
 
 
inline bool
IndexSet::IntervalIterator::operator<(
  const IndexSet::IntervalIterator &other) const
{
  return accessor < other.accessor;
}
 
 
 
inline int
IndexSet::IntervalIterator::operator-(
  const IndexSet::IntervalIterator &other) const
{
  Assert(accessor.index_set == other.accessor.index_set,
         ExcMessage(
           "Can not compare iterators belonging to different IndexSets"));
 
  const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index) ?
                          accessor.index_set->ranges.size() :
                          accessor.range_idx;
  const size_type rhs =
    (other.accessor.range_idx == numbers::invalid_dof_index) ?
      accessor.index_set->ranges.size() :
      other.accessor.range_idx;
 
  if (lhs > rhs)
    return static_cast<int>(lhs - rhs);
  else
    return -static_cast<int>(rhs - lhs);
}
 
 
 
/* ElementIterator */
 
inline IndexSet::ElementIterator::ElementIterator(
  const IndexSet           *idxset,
  const IndexSet::size_type range_idx,
  const IndexSet::size_type index)
  : index_set(idxset)
  , range_idx(range_idx)
  , idx(index)
{
  Assert(range_idx < index_set->ranges.size(),
         ExcMessage(
           "Invalid range index for IndexSet::ElementIterator constructor."));
  Assert(
    idx >= index_set->ranges[range_idx].begin &&
      idx < index_set->ranges[range_idx].end,
    ExcInternalError(
      "Invalid index argument for IndexSet::ElementIterator constructor."));
}
 
 
 
inline IndexSet::ElementIterator::ElementIterator(const IndexSet *idxset)
  : index_set(idxset)
  , range_idx(numbers::invalid_dof_index)
  , idx(numbers::invalid_dof_index)
{}
 
 
 
inline bool
IndexSet::ElementIterator::is_valid() const
{
  Assert((range_idx == numbers::invalid_dof_index &&
          idx == numbers::invalid_dof_index) ||
           (range_idx < index_set->ranges.size() &&
            idx < index_set->ranges[range_idx].end),
         ExcInternalError("Invalid ElementIterator state."));
 
  return (range_idx < index_set->ranges.size() &&
          idx < index_set->ranges[range_idx].end);
}
 
 
 
inline IndexSet::size_type
IndexSet::ElementIterator::operator*() const
{
  Assert(
    is_valid(),
    ExcMessage(
      "Impossible to dereference an IndexSet::ElementIterator that is invalid"));
  return idx;
}
 
 
 
inline bool
IndexSet::ElementIterator::operator==(
  const IndexSet::ElementIterator &other) const
{
  Assert(index_set == other.index_set,
         ExcMessage(
           "Can not compare iterators belonging to different IndexSets"));
  return range_idx == other.range_idx && idx == other.idx;
}
 
 
 
inline void
IndexSet::ElementIterator::advance()
{
  Assert(
    is_valid(),
    ExcMessage(
      "Impossible to advance an IndexSet::ElementIterator that is invalid"));
  if (idx < index_set->ranges[range_idx].end)
    ++idx;
  // end of this range?
  if (idx == index_set->ranges[range_idx].end)
    {
      // point to first element in next interval if possible
      if (range_idx < index_set->ranges.size() - 1)
        {
          ++range_idx;
          idx = index_set->ranges[range_idx].begin;
        }
      else
        {
          // we just fell off the end, set to invalid:
          range_idx = numbers::invalid_dof_index;
          idx       = numbers::invalid_dof_index;
        }
    }
}
 
 
 
inline IndexSet::ElementIterator &
IndexSet::ElementIterator::operator++()
{
  advance();
  return *this;
}
 
 
 
inline IndexSet::ElementIterator
IndexSet::ElementIterator::operator++(int)
{
  const IndexSet::ElementIterator it = *this;
  advance();
  return it;
}
 
 
 
inline bool
IndexSet::ElementIterator::operator!=(
  const IndexSet::ElementIterator &other) const
{
  return !(*this == other);
}
 
 
 
inline bool
IndexSet::ElementIterator::operator<(
  const IndexSet::ElementIterator &other) const
{
  Assert(index_set == other.index_set,
         ExcMessage(
           "Can not compare iterators belonging to different IndexSets"));
  return range_idx < other.range_idx ||
         (range_idx == other.range_idx && idx < other.idx);
}
 
 
 
inline std::ptrdiff_t
IndexSet::ElementIterator::operator-(
  const IndexSet::ElementIterator &other) const
{
  Assert(index_set == other.index_set,
         ExcMessage(
           "Can not compare iterators belonging to different IndexSets"));
  if (*this == other)
    return 0;
  if (!(*this < other))
    return -(other - *this);
 
  // only other can be equal to end() because of the checks above.
  Assert(is_valid(), ExcInternalError());
 
  // Note: we now compute how far advance *this in "*this < other" to get other,
  // so we need to return -c at the end.
 
  // first finish the current range:
  std::ptrdiff_t c = index_set->ranges[range_idx].end - idx;
 
  // now walk in steps of ranges (need to start one behind our current one):
  for (size_type range = range_idx + 1;
       range < index_set->ranges.size() && range <= other.range_idx;
       ++range)
    c += index_set->ranges[range].end - index_set->ranges[range].begin;
 
  Assert(
    other.range_idx < index_set->ranges.size() ||
      other.range_idx == numbers::invalid_dof_index,
    ExcMessage(
      "Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));
 
  // We might have walked too far because we went until the end of
  // other.range_idx, so walk backwards to other.idx:
  if (other.range_idx != numbers::invalid_dof_index)
    c -= index_set->ranges[other.range_idx].end - other.idx;
 
  return -c;
}
 
 
/* Range */
 
inline IndexSet::Range::Range()
  : begin(numbers::invalid_dof_index)
  , end(numbers::invalid_dof_index)
  , nth_index_in_set(numbers::invalid_dof_index)
{}
 
 
 
inline IndexSet::Range::Range(const size_type i1, const size_type i2)
  : begin(i1)
  , end(i2)
  , nth_index_in_set(numbers::invalid_dof_index)
{}
 
 
 
/* IndexSet itself */
 
inline IndexSet::IndexSet()
  : is_compressed(true)
  , index_space_size(0)
  , largest_range(numbers::invalid_unsigned_int)
{}
 
 
 
inline IndexSet::IndexSet(const size_type size)
  : is_compressed(true)
  , index_space_size(size)
  , largest_range(numbers::invalid_unsigned_int)
{}
 
 
 
inline IndexSet::IndexSet(IndexSet &&is) noexcept
  : ranges(std::move(is.ranges))
  , is_compressed(is.is_compressed)
  , index_space_size(is.index_space_size)
  , largest_range(is.largest_range)
{
  is.ranges.clear();
  is.is_compressed    = true;
  is.index_space_size = 0;
  is.largest_range    = numbers::invalid_unsigned_int;
 
  compress();
}
 
 
 
inline IndexSet &
IndexSet::operator=(IndexSet &&is) noexcept
{
  ranges           = std::move(is.ranges);
  is_compressed    = is.is_compressed;
  index_space_size = is.index_space_size;
  largest_range    = is.largest_range;
 
  is.ranges.clear();
  is.is_compressed    = true;
  is.index_space_size = 0;
  is.largest_range    = numbers::invalid_unsigned_int;
 
  compress();
 
  return *this;
}
 
 
 
inline IndexSet::ElementIterator
IndexSet::begin() const
{
  compress();
  if (ranges.size() > 0)
    return {this, 0, ranges[0].begin};
  else
    return end();
}
 
 
 
inline IndexSet::ElementIterator
IndexSet::end() const
{
  compress();
  return IndexSet::ElementIterator(this);
}
 
 
 
inline IndexSet::IntervalIterator
IndexSet::begin_intervals() const
{
  compress();
  if (ranges.size() > 0)
    return IndexSet::IntervalIterator(this, 0);
  else
    return end_intervals();
}
 
 
 
inline IndexSet::IntervalIterator
IndexSet::end_intervals() const
{
  compress();
  return IndexSet::IntervalIterator(this);
}
 
 
 
inline void
IndexSet::clear()
{
  // reset so that there are no indices in the set any more; however,
  // as documented, the index set retains its size
  ranges.clear();
  is_compressed = true;
  largest_range = numbers::invalid_unsigned_int;
}
 
 
 
inline void
IndexSet::set_size(const size_type sz)
{
  Assert(ranges.empty(),
         ExcMessage("This function can only be called if the current "
                    "object does not yet contain any elements."));
  index_space_size = sz;
  is_compressed    = true;
}
 
 
 
inline IndexSet::size_type
IndexSet::size() const
{
  return index_space_size;
}
 
 
 
inline void
IndexSet::compress() const
{
  if (is_compressed == true)
    return;
 
  do_compress();
}
 
 
 
inline void
IndexSet::add_index(const size_type index)
{
  add_range(index, index + 1);
}
 
 
 
inline void
IndexSet::add_range(const size_type begin, const size_type end)
{
  Assert((begin < index_space_size) ||
           ((begin == index_space_size) && (end == index_space_size)),
         ExcIndexRangeType<size_type>(begin, 0, index_space_size));
  Assert(end <= index_space_size,
         ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
  AssertIndexRange(begin, end + 1);
 
  if (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.empty() || begin > ranges.back().end)
        ranges.emplace_back(begin, end);
      else if (begin == ranges.back().end)
        ranges.back().end = end;
      else
        add_range_lower_bound(Range(begin, end));
 
      is_compressed = false;
    }
}
 
 
 
template <typename ForwardIterator>
inline void
IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
{
  if (begin == end)
    return;
 
  // identify ranges in the given iterator range by checking whether some
  // indices happen to be consecutive. to avoid quadratic complexity when
  // calling add_range many times (as add_range() going into the middle of an
  // already existing range must shift entries around), we first collect a
  // vector of ranges.
  boost::container::small_vector<std::pair<size_type, size_type>, 200>
       tmp_ranges;
  bool ranges_are_sorted = true;
  for (ForwardIterator p = begin; p != end;)
    {
      // Starting with the current iterator 'p', find an iterator
      // 'q' so that the indices pointed to by the iterators in
      // the range [p,q) are consecutive. These indices then form
      // a range that is contiguous, and that can be added all
      // at once.
      const size_type begin_index = *p;
      size_type       end_index   = begin_index + 1;
 
      // Start looking at the position after 'p', and keep iterating while
      // 'q' points to a duplicate of 'p':
      ForwardIterator q = p;
      ++q;
      while ((q != end) && (*q == *p))
        ++q;
 
      // Now we know that 'q' is either past the end, or points to a value
      // other than 'p'. If it points to 'end_index', we are still good with
      // a contiguous range; then increment the end index of that range, and
      // move to the next iterator that is not a duplicate of what
      // we were just looking at:
      while ((q != end) && (static_cast<size_type>(*q) == end_index))
        {
          ++q;
          while ((q != end) && (static_cast<size_type>(*q) == end_index))
            ++q;
 
          ++end_index;
        }
 
      // Add this range:
      tmp_ranges.emplace_back(begin_index, end_index);
 
      // Then move on to the next element in the input range.
      // If the starting index of the next go-around of the for loop is less
      // than the end index of the one just identified, then we will have at
      // least one pair of ranges that are not sorted, and consequently the
      // whole collection of ranges is not sorted.
      p = q;
      if ((p != end) && (static_cast<size_type>(*p) < end_index))
        ranges_are_sorted = false;
    }
 
  add_ranges_internal(tmp_ranges, ranges_are_sorted);
}
 
 
 
inline bool
IndexSet::is_element(const size_type index) const
{
  if (ranges.empty() == false)
    {
      compress();
 
      // 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;
      else if (ranges.size() > 1)
        return is_element_binary_search(index);
      else
        return false;
    }
  else
    return false;
}
 
 
 
inline bool
IndexSet::is_contiguous() const
{
  compress();
  return (ranges.size() <= 1);
}
 
 
 
inline bool
IndexSet::is_empty() const
{
  return ranges.empty();
}
 
 
 
inline IndexSet::size_type
IndexSet::n_elements() const
{
  // make sure we have non-overlapping ranges
  compress();
 
  size_type v = 0;
  if (!ranges.empty())
    {
      const Range &r = ranges.back();
      v              = r.nth_index_in_set + r.end - r.begin;
    }
 
#ifdef DEBUG
  size_type s = 0;
  for (const auto &range : ranges)
    s += (range.end - range.begin);
  Assert(s == v, ExcInternalError());
#endif
 
  return v;
}
 
 
 
inline unsigned int
IndexSet::n_intervals() const
{
  compress();
  return ranges.size();
}
 
 
 
inline IndexSet::size_type
IndexSet::largest_range_starting_index() const
{
  Assert(ranges.empty() == false, ExcMessage("IndexSet cannot be empty."));
 
  compress();
  const std::vector<Range>::const_iterator main_range =
    ranges.begin() + largest_range;
 
  return main_range->nth_index_in_set;
}
 
 
 
inline IndexSet::size_type
IndexSet::nth_index_in_set(const size_type n) const
{
  AssertIndexRange(n, n_elements());
 
  compress();
 
  // first check whether the index is in the largest range
  Assert(largest_range < ranges.size(), ExcInternalError());
  const auto 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);
  else
    return nth_index_in_set_binary_search(n);
}
 
 
 
inline IndexSet::size_type
IndexSet::index_within_set(const size_type n) const
{
  // to make this call thread-safe, compress() must not be called through this
  // function
  Assert(is_compressed == true, ExcMessage("IndexSet must be compressed."));
  AssertIndexRange(n, size());
 
  // return immediately if the index set is empty
  if (is_empty())
    return numbers::invalid_dof_index;
 
  // 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());
  if (n >= ranges[largest_range].begin && n < ranges[largest_range].end)
    return (n - ranges[largest_range].begin) +
           ranges[largest_range].nth_index_in_set;
  else if (ranges.size() > 1)
    return index_within_set_binary_search(n);
  else
    return numbers::invalid_dof_index;
}
 
 
 
inline bool
IndexSet::operator==(const IndexSet &is) const
{
  // If one of the two index sets has size zero, the other one has to
  // have size zero as well:
  if (size() == 0)
    return (is.size() == 0);
  if (is.size() == 0)
    return (size() == 0);
 
  // Otherwise, they must have the same size (see the documentation):
  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
 
  compress();
  is.compress();
 
  return (ranges == is.ranges);
}
 
 
 
inline bool
IndexSet::operator!=(const IndexSet &is) const
{
  // If one of the two index sets has size zero, the other one has to
  // have a non-zero size for inequality:
  if (size() == 0)
    return (is.size() != 0);
  if (is.size() == 0)
    return (size() != 0);
 
  // Otherwise, they must have the same size (see the documentation):
  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
 
  compress();
  is.compress();
 
  return (ranges != is.ranges);
}
 
 
 
template <typename Vector>
void
IndexSet::fill_binary_vector(Vector &vector) const
{
  Assert(vector.size() == size(), ExcDimensionMismatch(vector.size(), size()));
 
  compress();
  // 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
  for (const auto &range : ranges)
    for (size_type i = range.begin; i < range.end; ++i)
      vector[i] = 1;
}
 
 
 
template <typename StreamType>
inline void
IndexSet::print(StreamType &out) const
{
  compress();
  out << '{';
  std::vector<Range>::const_iterator p;
  for (p = ranges.begin(); p != ranges.end(); ++p)
    {
      if (p->end - p->begin == 1)
        out << p->begin;
      else
        out << '[' << p->begin << ',' << p->end - 1 << ']';
 
      if (p != --ranges.end())
        out << ", ";
    }
  out << '}' << std::endl;
}
 
 
 
template <class Archive>
inline void
IndexSet::Range::serialize(Archive &ar, const unsigned int)
{
  ar &begin &end &nth_index_in_set;
}
 
 
 
template <class Archive>
inline void
IndexSet::serialize(Archive &ar, const unsigned int)
{
  ar &ranges &is_compressed &index_space_size &largest_range;
}
 
DEAL_II_NAMESPACE_CLOSE
 
#endif