interface.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_INTERFACE_HH
#define DUNE_INTERFACE_HH
 
#if HAVE_MPI
 
#include "remoteindices.hh"
#include <dune/common/enumset.hh>
 
namespace Dune
{
  class InterfaceBuilder
  {
  public:
    class RemoteIndicesStateError : public InvalidStateException
    {};
 
    virtual ~InterfaceBuilder()
    {}
 
  protected:
    InterfaceBuilder()
    {}
 
    template<class R, class T1, class T2, class Op, bool send>
    void buildInterface (const R& remoteIndices,
                         const T1& sourceFlags, const T2& destFlags,
                         Op& functor) const;
  };
 
  class InterfaceInformation
  {
 
  public:
 
    size_t size() const
    {
      return size_;
    }
    std::size_t& operator[](size_t i)
    {
      assert(i<size_);
      return indices_[i];
    }
    std::size_t operator[](size_t i) const
    {
      assert(i<size_);
      return indices_[i];
    }
    void reserve(size_t size)
    {
      indices_ = new std::size_t[size];
      maxSize_ = size;
 
    }
    void free()
    {
      if(indices_)
        delete[] indices_;
      maxSize_ = 0;
      size_=0;
      indices_=0;
    }
    void add(std::size_t index)
    {
      assert(size_<maxSize_);
      indices_[size_++]=index;
    }
 
    InterfaceInformation()
      : size_(0), maxSize_(0), indices_(0)
    {}
 
    virtual ~InterfaceInformation()
    {}
 
    bool operator!=(const InterfaceInformation& o) const
    {
      return !operator==(o);
    }
 
    bool operator==(const InterfaceInformation& o) const
    {
      if(size_!=o.size_)
        return false;
      for(std::size_t i=0; i< size_; ++i)
        if(indices_[i]!=o.indices_[i])
          return false;
      return true;
    }
 
  private:
    size_t size_;
    size_t maxSize_;
    std::size_t* indices_;
  };
 
  class Interface : public InterfaceBuilder
  {
 
  public:
    typedef std::map<int,std::pair<InterfaceInformation,InterfaceInformation> > InformationMap;
 
    template<typename R, typename T1, typename T2>
    void build(const R& remoteIndices, const T1& sourceFlags,
               const T2& destFlags);
 
    void free();
 
    MPI_Comm communicator() const;
 
    const InformationMap& interfaces() const;
 
    Interface(MPI_Comm comm)
      : communicator_(comm), interfaces_()
    {}
 
    Interface()
      : communicator_(MPI_COMM_NULL), interfaces_()
    {}
 
    void print() const;
 
    bool operator!=(const Interface& o) const
    {
      return ! operator==(o);
    }
 
    bool operator==(const Interface& o) const
    {
      if(communicator_!=o.communicator_)
        return false;
      if(interfaces_.size()!=o.interfaces_.size())
        return false;
      typedef InformationMap::const_iterator MIter;
 
      for(MIter m=interfaces_.begin(), om=o.interfaces_.begin();
          m!=interfaces_.end(); ++m, ++om)
      {
        if(om->first!=m->first)
          return false;
        if(om->second.first!=om->second.first)
          return false;
        if(om->second.second!=om->second.second)
          return false;
      }
      return true;
    }
 
    virtual ~Interface();
 
    void strip();
  protected:
 
    InformationMap& interfaces();
 
    MPI_Comm communicator_;
 
  private:
    InformationMap interfaces_;
 
    template<bool send>
    class InformationBuilder
    {
    public:
      InformationBuilder(InformationMap& interfaces)
        : interfaces_(interfaces)
      {}
 
      void reserve(int proc, int size)
      {
        if(send)
          interfaces_[proc].first.reserve(size);
        else
          interfaces_[proc].second.reserve(size);
      }
      void add(int proc, std::size_t local)
      {
        if(send) {
          interfaces_[proc].first.add(local);
        }else{
          interfaces_[proc].second.add(local);
        }
      }
 
    private:
      InformationMap& interfaces_;
    };
  };
 
  template<class R, class T1, class T2, class Op, bool send>
  void InterfaceBuilder::buildInterface(const R& remoteIndices, const T1& sourceFlags, const T2& destFlags, Op& interfaceInformation) const
  {
 
    if(!remoteIndices.isSynced())
      DUNE_THROW(RemoteIndicesStateError,"RemoteIndices is not in sync with the index set. Call RemoteIndices::rebuild first!");
    // Allocate the memory for the data type construction.
    typedef R RemoteIndices;
    typedef typename RemoteIndices::RemoteIndexMap::const_iterator const_iterator;
 
    const const_iterator end=remoteIndices.end();
 
    int rank;
 
    MPI_Comm_rank(remoteIndices.communicator(), &rank);
 
    // Allocate memory for the type construction.
    for(const_iterator process=remoteIndices.begin(); process != end; ++process) {
      // Messure the number of indices send to the remote process first
      int size=0;
      typedef typename RemoteIndices::RemoteIndexList::const_iterator RemoteIterator;
      const RemoteIterator remoteEnd = send ? process->second.first->end() :
                                       process->second.second->end();
      RemoteIterator remote = send ? process->second.first->begin() : process->second.second->begin();
 
      while(remote!=remoteEnd) {
        if( send ?  destFlags.contains(remote->attribute()) :
            sourceFlags.contains(remote->attribute())) {
 
          // do we send the index?
          if( send ? sourceFlags.contains(remote->localIndexPair().local().attribute()) :
              destFlags.contains(remote->localIndexPair().local().attribute()))
            ++size;
        }
        ++remote;
      }
      interfaceInformation.reserve(process->first, size);
    }
 
    // compare the local and remote indices and set up the types
 
    for(const_iterator process=remoteIndices.begin(); process != end; ++process) {
      typedef typename RemoteIndices::RemoteIndexList::const_iterator RemoteIterator;
      const RemoteIterator remoteEnd = send ? process->second.first->end() :
                                       process->second.second->end();
      RemoteIterator remote = send ? process->second.first->begin() : process->second.second->begin();
 
      while(remote!=remoteEnd) {
        if( send ?  destFlags.contains(remote->attribute()) :
            sourceFlags.contains(remote->attribute())) {
          // do we send the index?
          if( send ? sourceFlags.contains(remote->localIndexPair().local().attribute()) :
              destFlags.contains(remote->localIndexPair().local().attribute()))
            interfaceInformation.add(process->first,remote->localIndexPair().local().local());
        }
        ++remote;
      }
    }
  }
 
  inline MPI_Comm Interface::communicator() const
  {
    return communicator_;
 
  }
 
 
  inline const std::map<int,std::pair<InterfaceInformation,InterfaceInformation> >& Interface::interfaces() const
  {
    return interfaces_;
  }
 
  inline std::map<int,std::pair<InterfaceInformation,InterfaceInformation> >& Interface::interfaces()
  {
    return interfaces_;
  }
 
  inline void Interface::print() const
  {
    typedef InformationMap::const_iterator const_iterator;
    const const_iterator end=interfaces_.end();
    int rank;
    MPI_Comm_rank(communicator(), &rank);
 
    for(const_iterator infoPair=interfaces_.begin(); infoPair!=end; ++infoPair) {
      {
        std::cout<<rank<<": send for process "<<infoPair->first<<": ";
        const InterfaceInformation& info(infoPair->second.first);
        for(size_t i=0; i < info.size(); i++)
          std::cout<<info[i]<<" ";
        std::cout<<std::endl;
      } {
 
        std::cout<<rank<<": receive for process "<<infoPair->first<<": ";
        const InterfaceInformation& info(infoPair->second.second);
        for(size_t i=0; i < info.size(); i++)
          std::cout<<info[i]<<" ";
        std::cout<<std::endl;
      }
 
    }
  }
 
  template<typename R, typename T1, typename T2>
  inline void Interface::build(const R& remoteIndices, const T1& sourceFlags,
                               const T2& destFlags)
  {
    communicator_=remoteIndices.communicator();
 
    assert(interfaces_.empty());
 
    // Build the send interface
    InformationBuilder<true> sendInformation(interfaces_);
    this->template buildInterface<R,T1,T2,InformationBuilder<true>,true>(remoteIndices, sourceFlags,
                                                                         destFlags, sendInformation);
 
    // Build the receive interface
    InformationBuilder<false> recvInformation(interfaces_);
    this->template buildInterface<R,T1,T2,InformationBuilder<false>,false>(remoteIndices,sourceFlags,
                                                                           destFlags, recvInformation);
    strip();
  }
  inline void Interface::strip()
  {
    typedef InformationMap::iterator const_iterator;
    for(const_iterator interfacePair = interfaces_.begin(); interfacePair != interfaces_.end();)
      if(interfacePair->second.first.size()==0 && interfacePair->second.second.size()==0) {
        interfacePair->second.first.free();
        interfacePair->second.second.free();
        const_iterator toerase=interfacePair++;
        interfaces_.erase(toerase);
      }else
        ++interfacePair;
  }
 
  inline void Interface::free()
  {
    typedef InformationMap::iterator iterator;
    typedef InformationMap::const_iterator const_iterator;
    const const_iterator end = interfaces_.end();
    for(iterator interfacePair = interfaces_.begin(); interfacePair != end; ++interfacePair) {
      interfacePair->second.first.free();
      interfacePair->second.second.free();
    }
    interfaces_.clear();
  }
 
  inline Interface::~Interface()
  {
    free();
  }
}
 
namespace std
{
  inline ostream& operator<<(ostream& os, const Dune::Interface& interface)
  {
    typedef Dune::Interface::InformationMap InfoMap;
    typedef InfoMap::const_iterator Iter;
    for(Iter i=interface.interfaces().begin(), end = interface.interfaces().end();
        i!=end; ++i)
    {
      os<<i->first<<": [ source=[";
      for(std::size_t j=0; j < i->second.first.size(); ++j)
        os<<i->second.first[j]<<" ";
      os<<"] size="<<i->second.first.size()<<", target=[";
      for(std::size_t j=0; j < i->second.second.size(); ++j)
        os<<i->second.second[j]<<" ";
      os<<"] size="<<i->second.second.size()<<"\n";
    }
    return os;
  }
} // end namespace std
#endif // HAVE_MPI
 
#endif