indexset.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_INDEXSET_HH
#define DUNE_INDEXSET_HH
 
#include <algorithm>
#include <dune/common/arraylist.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/unused.hh>
#include <iostream>
 
#include "localindex.hh"
 
#include <stdint.h> // for uint32_t
 
namespace Dune
{
  // forward declarations
 
  template<class TG, class TL>
  class IndexPair;
 
  template<class TG, class TL>
  std::ostream& operator<<(std::ostream& os, const IndexPair<TG,TL>& pair);
 
  template<class TG, class TL>
  bool operator==(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
 
  template<class TG, class TL>
  bool operator!=(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
 
  template<class TG, class TL>
  bool operator<(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
 
  template<class TG, class TL>
  bool operator>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
 
  template<class TG, class TL>
  bool operator<=(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
 
  template<class TG, class TL>
  bool operator >=(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
 
  template<class TG, class TL>
  bool operator==(const IndexPair<TG,TL>&, const TG&);
 
  template<class TG, class TL>
  bool operator!=(const IndexPair<TG,TL>&, const TG&);
 
  template<class TG, class TL>
  bool operator<(const IndexPair<TG,TL>&, const TG&);
 
  template<class TG, class TL>
  bool operator>(const IndexPair<TG,TL>&, const TG&);
 
  template<class TG, class TL>
  bool operator<=(const IndexPair<TG,TL>&, const TG&);
 
  template<class TG, class TL>
  bool operator >=(const IndexPair<TG,TL>&, const TG&);
 
  template<typename T>
  struct MPITraits;
 
  template<class TG, class TL>
  class IndexPair
  {
    friend std::ostream& operator<<<>(std::ostream&, const IndexPair<TG,TL>&);
    friend bool operator==<>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
    friend bool operator!=<>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
    friend bool operator< <>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
    friend bool operator><>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
    friend bool operator<=<>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
    friend bool operator>=<>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
    friend bool operator==<>(const IndexPair<TG,TL>&, const TG &);
    friend bool operator!=<>(const IndexPair<TG,TL>&, const TG &);
    friend bool operator< <>(const IndexPair<TG,TL>&, const TG &);
    friend bool operator> <>(const IndexPair<TG,TL>&, const TG &);
    friend bool operator<=<>(const IndexPair<TG,TL>&, const TG &);
    friend bool operator>=<>(const IndexPair<TG,TL>&, const TG &);
    friend struct MPITraits<IndexPair<TG,TL> >;
 
  public:
    typedef TG GlobalIndex;
 
    typedef TL LocalIndex;
 
    IndexPair(const GlobalIndex& global, const LocalIndex& local);
 
    IndexPair();
    IndexPair(const GlobalIndex& global);
 
    inline const GlobalIndex& global() const;
 
    inline LocalIndex& local();
 
    inline const LocalIndex& local() const;
 
    inline void setLocal(int index);
  private:
    GlobalIndex global_;
    LocalIndex local_;
  };
 
  enum ParallelIndexSetState
  {
    GROUND,
    RESIZE
  };
 
  class InvalidIndexSetState : public InvalidStateException {};
 
  // Forward declaration
  template<class I> class GlobalLookupIndexSet;
 
  template<typename TG, typename TL, int N=100>
  class ParallelIndexSet
  {
    friend class GlobalLookupIndexSet<ParallelIndexSet<TG,TL,N> >;
 
  public:
    typedef TG GlobalIndex;
 
    typedef TL LocalIndex;
 
    typedef Dune::IndexPair<GlobalIndex,LocalIndex> IndexPair;
 
    enum {
      arraySize= (N>0) ? N : 1
    };
 
    class iterator :
      public ArrayList<IndexPair,N>::iterator
    {
      typedef typename ArrayList<IndexPair,N>::iterator
      Father;
      friend class ParallelIndexSet<GlobalIndex,LocalIndex,N>;
    public:
      iterator(ParallelIndexSet<TG,TL,N>& indexSet, const Father& father)
        : Father(father), indexSet_(&indexSet)
      {}
 
      iterator(const iterator& other)
        : Father(other), indexSet_(other.indexSet_)
      {}
 
      iterator& operator==(const iterator& other)
      {
        Father::operator==(other);
        indexSet_ = other.indexSet_;
      }
 
    private:
      inline void markAsDeleted() const throw(InvalidIndexSetState)
      {
#ifndef NDEBUG
        if(indexSet_->state_ != RESIZE)
          DUNE_THROW(InvalidIndexSetState, "Indices can only be removed "
                     <<"while in RESIZE state!");
#endif
        Father::operator*().local().setState(DELETED);
      }
 
      ParallelIndexSet<TG,TL,N>* indexSet_;
 
    };
 
 
 
    typedef typename
    ArrayList<IndexPair,N>::const_iterator
    const_iterator;
 
    ParallelIndexSet();
 
    inline const ParallelIndexSetState& state()
    {
      return state_;
    }
 
    void beginResize() throw(InvalidIndexSetState);
 
    inline void add(const GlobalIndex& global) throw(InvalidIndexSetState);
 
    inline void add(const GlobalIndex& global, const LocalIndex& local)
    throw(InvalidIndexSetState);
 
    inline void markAsDeleted(const iterator& position)
    throw(InvalidIndexSetState);
 
    void endResize() throw(InvalidIndexSetState);
 
    inline IndexPair&
    operator[](const GlobalIndex& global);
 
    inline IndexPair&
    at(const GlobalIndex& global);
 
    inline const IndexPair&
    operator[](const GlobalIndex& global) const;
 
    inline const IndexPair&
    at(const GlobalIndex& global) const;
 
    inline iterator begin();
 
    inline iterator end();
 
    inline const_iterator begin() const;
 
    inline const_iterator end() const;
 
    inline void renumberLocal();
 
    inline int seqNo() const;
 
    inline size_t size() const;
 
  private:
    ArrayList<IndexPair,N> localIndices_;
    ArrayList<IndexPair,N> newIndices_;
    ParallelIndexSetState state_;
    int seqNo_;
    bool deletedEntries_;
    inline void merge();
  };
 
 
  template<class TG, class TL, int N>
  std::ostream& operator<<(std::ostream& os, const ParallelIndexSet<TG,TL,N>& indexSet);
 
  template<class I>
  class GlobalLookupIndexSet
  {
  public:
    typedef I ParallelIndexSet;
 
    typedef typename ParallelIndexSet::LocalIndex LocalIndex;
 
    typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
 
    typedef typename ParallelIndexSet::const_iterator const_iterator;
 
    typedef Dune::IndexPair<typename I::GlobalIndex, typename I::LocalIndex> IndexPair;
 
    GlobalLookupIndexSet(const ParallelIndexSet& indexset, std::size_t size);
 
    GlobalLookupIndexSet(const ParallelIndexSet& indexset);
 
    ~GlobalLookupIndexSet();
 
    inline const IndexPair&
    operator[](const GlobalIndex& global) const;
 
    inline const IndexPair*
    pair(const std::size_t& local) const;
 
    inline const_iterator begin() const;
 
    inline const_iterator end() const;
 
    inline int seqNo() const;
 
    inline size_t size() const;
  private:
    const ParallelIndexSet& indexSet_;
 
    std::size_t size_;
 
    std::vector<const IndexPair*> indices_;
 
  };
 
 
  template<typename T>
  struct LocalIndexComparator
  {
    static bool compare(const T& t1, const T& t2){
      DUNE_UNUSED_PARAMETER(t1);
      DUNE_UNUSED_PARAMETER(t2);
      return false;
    }
  };
 
  template<class TG, class TL>
  struct IndexSetSortFunctor
  {
    bool operator()(const IndexPair<TG,TL>& i1, const IndexPair<TG,TL>& i2)
    {
      return i1.global()<i2.global() || (i1.global()==i2.global() &&
                                         LocalIndexComparator<TL>::compare(i1.local(),
                                                                           i2.local()));
    }
  };
 
 
 
  template<class TG, class TL>
  inline std::ostream& operator<<(std::ostream& os, const IndexPair<TG,TL>& pair)
  {
    os<<"{global="<<pair.global_<<", local="<<pair.local_<<"}";
    return os;
  }
 
  template<class TG, class TL, int N>
  inline std::ostream& operator<<(std::ostream& os, const ParallelIndexSet<TG,TL,N>& indexSet)
  {
    typedef typename ParallelIndexSet<TG,TL,N>::const_iterator Iterator;
    Iterator end = indexSet.end();
    os<<"{";
    for(Iterator index = indexSet.begin(); index != end; ++index)
      os<<*index<<" ";
    os<<"}";
    return os;
 
  }
 
  template<class TG, class TL>
  inline bool operator==(const IndexPair<TG,TL>& a, const IndexPair<TG,TL>& b)
  {
    return a.global_==b.global_;
  }
 
  template<class TG, class TL>
  inline bool operator!=(const IndexPair<TG,TL>& a, const IndexPair<TG,TL>& b)
  {
    return a.global_!=b.global_;
  }
 
  template<class TG, class TL>
  inline bool operator<(const IndexPair<TG,TL>& a, const IndexPair<TG,TL>& b)
  {
    return a.global_<b.global_;
  }
 
  template<class TG, class TL>
  inline bool operator>(const IndexPair<TG,TL>& a, const IndexPair<TG,TL>& b)
  {
    return a.global_>b.global_;
  }
 
  template<class TG, class TL>
  inline bool operator<=(const IndexPair<TG,TL>& a, const IndexPair<TG,TL>& b)
  {
    return a.global_<=b.global_;
  }
 
  template<class TG, class TL>
  inline bool operator >=(const IndexPair<TG,TL>& a, const IndexPair<TG,TL>& b)
  {
    return a.global_>=b.global_;
  }
 
  template<class TG, class TL>
  inline bool operator==(const IndexPair<TG,TL>& a, const TG& b)
  {
    return a.global_==b;
  }
 
  template<class TG, class TL>
  inline bool operator!=(const IndexPair<TG,TL>& a, const TG& b)
  {
    return a.global_!=b;
  }
 
  template<class TG, class TL>
  inline bool operator<(const IndexPair<TG,TL>& a, const TG& b)
  {
    return a.global_<b;
  }
 
  template<class TG, class TL>
  inline bool operator>(const IndexPair<TG,TL>& a, const TG& b)
  {
    return a.global_>b;
  }
 
  template<class TG, class TL>
  inline bool operator<=(const IndexPair<TG,TL>& a, const TG& b)
  {
    return a.global_<=b;
  }
 
  template<class TG, class TL>
  inline bool operator >=(const IndexPair<TG,TL>& a, const TG& b)
  {
    return a.global_>=b;
  }
 
#ifndef DOXYGEN
 
  template<class TG, class TL>
  IndexPair<TG,TL>::IndexPair(const TG& global, const TL& local)
    : global_(global), local_(local){}
 
  template<class TG, class TL>
  IndexPair<TG,TL>::IndexPair(const TG& global)
    : global_(global), local_(){}
 
  template<class TG, class TL>
  IndexPair<TG,TL>::IndexPair()
    : global_(), local_(){}
 
  template<class TG, class TL>
  inline const TG& IndexPair<TG,TL>::global() const {
    return global_;
  }
 
  template<class TG, class TL>
  inline TL& IndexPair<TG,TL>::local() {
    return local_;
  }
 
  template<class TG, class TL>
  inline const TL& IndexPair<TG,TL>::local() const {
    return local_;
  }
 
  template<class TG, class TL>
  inline void IndexPair<TG,TL>::setLocal(int local){
    local_=local;
  }
 
  template<class TG, class TL, int N>
  ParallelIndexSet<TG,TL,N>::ParallelIndexSet()
    : state_(GROUND), seqNo_(0)
  {}
 
  template<class TG, class TL, int N>
  void ParallelIndexSet<TG,TL,N>::beginResize() throw(InvalidIndexSetState)
  {
 
    // Checks in unproductive code
#ifndef NDEBUG
    if(state_!=GROUND)
      DUNE_THROW(InvalidIndexSetState,
                 "IndexSet has to be in GROUND state, when "
                 << "beginResize() is called!");
#endif
 
    state_ = RESIZE;
    deletedEntries_ = false;
  }
 
  template<class TG, class TL, int N>
  inline void ParallelIndexSet<TG,TL,N>::add(const GlobalIndex& global)
  throw(InvalidIndexSetState)
  {
    // Checks in unproductive code
#ifndef NDEBUG
    if(state_ != RESIZE)
      DUNE_THROW(InvalidIndexSetState, "Indices can only be added "
                 <<"while in RESIZE state!");
#endif
    newIndices_.push_back(IndexPair(global));
  }
 
  template<class TG, class TL, int N>
  inline void ParallelIndexSet<TG,TL,N>::add(const TG& global, const TL& local)
  throw(InvalidIndexSetState)
  {
    // Checks in unproductive code
#ifndef NDEBUG
    if(state_ != RESIZE)
      DUNE_THROW(InvalidIndexSetState, "Indices can only be added "
                 <<"while in RESIZE state!");
#endif
    newIndices_.push_back(IndexPair(global,local));
  }
 
  template<class TG, class TL, int N>
  inline void ParallelIndexSet<TG,TL,N>::markAsDeleted(const iterator& global)
  throw(InvalidIndexSetState){
    // Checks in unproductive code
#ifndef NDEBUG
    if(state_ != RESIZE)
      DUNE_THROW(InvalidIndexSetState, "Indices can only be removed "
                 <<"while in RESIZE state!");
#endif
    deletedEntries_ = true;
 
    global.markAsDeleted();
  }
 
  template<class TG, class TL, int N>
  void ParallelIndexSet<TG,TL,N>::endResize() throw(InvalidIndexSetState){
    // Checks in unproductive code
#ifndef NDEBUG
    if(state_ != RESIZE)
      DUNE_THROW(InvalidIndexSetState, "endResize called while not "
                 <<"in RESIZE state!");
#endif
 
    std::sort(newIndices_.begin(), newIndices_.end(), IndexSetSortFunctor<TG,TL>());
    merge();
    seqNo_++;
    state_ = GROUND;
  }
 
 
  template<class TG, class TL, int N>
  inline void ParallelIndexSet<TG,TL,N>::merge(){
    if(localIndices_.size()==0)
    {
      localIndices_=newIndices_;
      newIndices_.clear();
    }
    else if(newIndices_.size()>0 || deletedEntries_)
    {
      ArrayList<IndexPair,N> tempPairs;
      typedef typename ArrayList<IndexPair,N>::iterator iterator;
      typedef typename ArrayList<IndexPair,N>::const_iterator const_iterator;
 
      iterator old=localIndices_.begin();
      iterator added=newIndices_.begin();
      const const_iterator endold=localIndices_.end();
      const const_iterator endadded=newIndices_.end();
 
      while(old != endold && added!= endadded)
      {
        if(old->local().state()==DELETED) {
          old.eraseToHere();
        }
        else
        {
          if(old->global() < added->global() ||
             (old->global() == added->global()
              && LocalIndexComparator<TL>::compare(old->local(),added->local())))
          {
            tempPairs.push_back(*old);
            old.eraseToHere();
            continue;
          }else
          {
            tempPairs.push_back(*added);
            added.eraseToHere();
          }
        }
      }
 
      while(old != endold)
      {
        if(old->local().state()!=DELETED) {
          tempPairs.push_back(*old);
        }
        old.eraseToHere();
      }
 
      while(added!= endadded)
      {
        tempPairs.push_back(*added);
        added.eraseToHere();
      }
      localIndices_ = tempPairs;
    }
  }
 
 
  template<class TG, class TL, int N>
  inline const IndexPair<TG,TL>&
  ParallelIndexSet<TG,TL,N>::at(const TG& global) const
  {
    // perform a binary search
    int low=0, high=localIndices_.size()-1, probe=-1;
 
    while(low<high)
    {
      probe = (high + low) / 2;
      if(global <= localIndices_[probe].global())
        high = probe;
      else
        low = probe+1;
    }
 
    if(probe==-1)
      DUNE_THROW(RangeError, "No entries!");
 
    if( localIndices_[low].global() != global)
      DUNE_THROW(RangeError, "Could not find entry of "<<global);
    else
      return localIndices_[low];
  }
 
  template<class TG, class TL, int N>
  inline const IndexPair<TG,TL>&
  ParallelIndexSet<TG,TL,N>::operator[](const TG& global) const
  {
    // perform a binary search
    int low=0, high=localIndices_.size()-1, probe=-1;
 
    while(low<high)
    {
      probe = (high + low) / 2;
      if(global <= localIndices_[probe].global())
        high = probe;
      else
        low = probe+1;
    }
 
    return localIndices_[low];
  }
  template<class TG, class TL, int N>
  inline IndexPair<TG,TL>& ParallelIndexSet<TG,TL,N>::at(const TG& global)
  {
    // perform a binary search
    int low=0, high=localIndices_.size()-1, probe=-1;
 
    while(low<high)
    {
      probe = (high + low) / 2;
      if(localIndices_[probe].global() >= global)
        high = probe;
      else
        low = probe+1;
    }
 
    if(probe==-1)
      DUNE_THROW(RangeError, "No entries!");
 
    if( localIndices_[low].global() != global)
      DUNE_THROW(RangeError, "Could not find entry of "<<global);
    else
      return localIndices_[low];
  }
 
  template<class TG, class TL, int N>
  inline IndexPair<TG,TL>& ParallelIndexSet<TG,TL,N>::operator[](const TG& global)
  {
    // perform a binary search
    int low=0, high=localIndices_.size()-1, probe=-1;
 
    while(low<high)
    {
      probe = (high + low) / 2;
      if(localIndices_[probe].global() >= global)
        high = probe;
      else
        low = probe+1;
    }
 
    return localIndices_[low];
  }
  template<class TG, class TL, int N>
  inline typename ParallelIndexSet<TG,TL,N>::iterator
  ParallelIndexSet<TG,TL,N>::begin()
  {
    return iterator(*this, localIndices_.begin());
  }
 
 
  template<class TG, class TL, int N>
  inline typename ParallelIndexSet<TG,TL,N>::iterator
  ParallelIndexSet<TG,TL,N>::end()
  {
    return iterator(*this,localIndices_.end());
  }
 
  template<class TG, class TL, int N>
  inline typename ParallelIndexSet<TG,TL,N>::const_iterator
  ParallelIndexSet<TG,TL,N>::begin() const
  {
    return localIndices_.begin();
  }
 
 
  template<class TG, class TL, int N>
  inline typename ParallelIndexSet<TG,TL,N>::const_iterator
  ParallelIndexSet<TG,TL,N>::end() const
  {
    return localIndices_.end();
  }
 
  template<class TG, class TL, int N>
  void ParallelIndexSet<TG,TL,N>::renumberLocal(){
#ifndef NDEBUG
    if(state_==RESIZE)
      DUNE_THROW(InvalidIndexSetState, "IndexSet has to be in "
                 <<"GROUND state for renumberLocal()");
#endif
 
    typedef typename ArrayList<IndexPair,N>::iterator iterator;
    const const_iterator end_ = end();
    uint32_t index=0;
 
    for(iterator pair=begin(); pair!=end_; index++, ++pair)
      pair->local()=index;
  }
 
  template<class TG, class TL, int N>
  inline int ParallelIndexSet<TG,TL,N>::seqNo() const
  {
    return seqNo_;
  }
 
  template<class TG, class TL, int N>
  inline size_t ParallelIndexSet<TG,TL,N>::size() const
  {
    return localIndices_.size();
  }
 
  template<class I>
  GlobalLookupIndexSet<I>::GlobalLookupIndexSet(const I& indexset,
                                                std::size_t size)
    : indexSet_(indexset), size_(size),
      indices_(size_, static_cast<const IndexPair*>(0))
  {
    const_iterator end_ = indexSet_.end();
 
    for(const_iterator pair = indexSet_.begin(); pair!=end_; ++pair) {
      assert(pair->local()<size_);
      indices_[pair->local()] = &(*pair);
    }
  }
 
  template<class I>
  GlobalLookupIndexSet<I>::GlobalLookupIndexSet(const I& indexset)
    : indexSet_(indexset), size_(0)
  {
    const_iterator end_ = indexSet_.end();
    for(const_iterator pair = indexSet_.begin(); pair!=end_; ++pair)
      size_=std::max(size_,static_cast<std::size_t>(pair->local()));
 
    indices_.resize(++size_,  0);
 
    for(const_iterator pair = indexSet_.begin(); pair!=end_; ++pair)
      indices_[pair->local()] = &(*pair);
  }
 
  template<class I>
  GlobalLookupIndexSet<I>::~GlobalLookupIndexSet()
  {}
 
  template<class I>
  inline const IndexPair<typename I::GlobalIndex, typename I::LocalIndex>*
  GlobalLookupIndexSet<I>::pair(const std::size_t& local) const
  {
    return indices_[local];
  }
 
  template<class I>
  inline const IndexPair<typename I::GlobalIndex, typename I::LocalIndex>&
  GlobalLookupIndexSet<I>::operator[](const GlobalIndex& global) const
  {
    return indexSet_[global];
  }
 
  template<class I>
  typename I::const_iterator GlobalLookupIndexSet<I>::begin() const
  {
    return indexSet_.begin();
  }
 
  template<class I>
  typename I::const_iterator GlobalLookupIndexSet<I>::end() const
  {
    return indexSet_.end();
  }
 
  template<class I>
  inline size_t GlobalLookupIndexSet<I>::size() const
  {
    return size_;
  }
 
  template<class I>
  inline int GlobalLookupIndexSet<I>::seqNo() const
  {
    return indexSet_.seqNo();
  }
 
  template<typename TG, typename TL, int N, typename TG1, typename TL1, int N1>
  bool operator==(const ParallelIndexSet<TG,TL,N>& idxset,
                  const ParallelIndexSet<TG1,TL1,N1>& idxset1)
  {
    if(idxset.size()!=idxset1.size())
      return false;
    typedef typename ParallelIndexSet<TG,TL,N>::const_iterator Iter;
    typedef typename ParallelIndexSet<TG1,TL1,N1>::const_iterator Iter1;
    Iter iter=idxset.begin();
    for(Iter1 iter1=idxset1.begin(); iter1 != idxset1.end(); ++iter, ++iter1) {
      if(iter1->global()!=iter->global())
        return false;
      typedef typename ParallelIndexSet<TG,TL,N>::LocalIndex PI;
      const PI& pi=iter->local(), pi1=iter1->local();
 
      if(pi!=pi1)
        return false;
    }
    return true;
  }
 
  template<typename TG, typename TL, int N, typename TG1, typename TL1, int N1>
  bool operator!=(const ParallelIndexSet<TG,TL,N>& idxset,
                  const ParallelIndexSet<TG1,TL1,N1>& idxset1)
  {
    return !(idxset==idxset1);
  }
 
 
#endif // DOXYGEN
 
}
#endif