globalindex.hh

/** \file
*
*  \brief Implement the functionality that provides a globally unique index for
*         all kinds of entities of a Dune grid. Such functionality is relevant
*         for a number of cases:
*         e.g. to map a degree of freedom associated wit an entity to its
*         location in a global matrix or global vector; e.g. to associate data
*         with entities in an adaptive mesh refinement scheme.
*
*  Nota bene: the functionality to provide a globally unique index for an entity, over all
*             processes in a distributed memory parallel compute environment is at present
*  not available in the Dune framework; therefore, we implement such functionality in a class
*  that is templatized by the GridView Type and the codimension of the entity for which the
*  globaly unique index shall be computed.
*
*  Method: (1) we use the UniqueEntityPartition class that reports if a specific entity belongs
*              to a specific process;
*
*          (2) we compute the number of entities that are owned by the specific process;
*
*          (3) we communicate the index of entities that are owned by the process to processes
*              that also contain these entities but do not own them, so that on a non-owner process
*              we have information on the index of the entity that it got from the owner-process;
*
*  Thus, we can access the same element in a global, distributed matrix that is managed
*  by all processes; a typical case would be using matrix and vector routines from the PETSc
*  or trilinos parallel linear algebra packages for distributed memory parallel computers.
*
*  Copyright by Benedikt Oswald and Patrick Leidenberger, 2002-2010. All rights reserved.
*
*  \author    Benedikt Oswald, Patrick Leidenberger
*
*  \warning   None
*
*  \attention globally unique indices are ONLY provided for entities of the
*             InteriorBorder_Partition type, NOT for the Ghost_Partition type !!!
*/

#ifndef GLOBALUNIQUEINDEX_HH_
#define GLOBALUNIQUEINDEX_HH_

/** \brief Include standard header files. */
#include <vector>
#include <iostream>
#include <fstream>
#include <map>
#include <utility>

/** include base class functionality for the communication interface */
#include <dune/grid/common/datahandleif.hh>

// Include proprietary header files.
#include "uniqueentitypartition.hh"

#include <dune/common/version.hh>

/** include parallel capability */
#if HAVE_MPI
  #include <dune/common/parallel/mpihelper.hh>
#endif


/*********************************************************************************************/
/* calculate globally unique index over all processes in a Dune grid that is distributed     */
/* over all the processes in a distributed memory parallel environment, for all entities     */
/* of a given codim in a given GridView, assuming they all have the same geometry,           */
/* i.e. codim, type                                                                          */
/*********************************************************************************************/
template<class GridView, int CODIM>
class GlobalUniqueIndex
{
private:
  /** define data types */
  typedef typename GridView::Grid Grid;

  typedef typename GridView::Grid::GlobalIdSet GlobalIdSet;
  typedef typename GridView::Grid::GlobalIdSet::IdType IdType;
  typedef typename GridView::Traits::template Codim<0>::Iterator Iterator;
  typedef typename GridView::Traits::template Codim<CODIM>::Entity Entity;


#if HAVE_MPI
  /** define data type related to collective communication */
  typedef typename Dune::template CollectiveCommunication<MPI_Comm> CollectiveCommunication;
#else
  typedef typename Grid::CollectiveCommunication CollectiveCommunication;
#endif

  typedef std::map<IdType,int> MapId2Index;
  typedef std::map<int,int>    IndexMap;

  typedef UniqueEntityPartition<GridView,CODIM> UniqueEntityPartitionType;

private:
  /* A DataHandle class to communicate the global index from the
   * owning to the non-owning entity; the class is based on the MinimumExchange
   * class in the parallelsolver.hh header file.
   */
  template<class GlobalIdSet, class IdType, class Map, typename ValueType>
  class IndexExchange           :       public Dune::CommDataHandleIF<IndexExchange<GlobalIdSet,IdType,Map,ValueType>,ValueType>
  {
  public:
    //! returns true if data for this codim should be communicated
    bool contains (int dim, int codim) const
    {
      return(codim==CODIM);
    }

    //! returns true if size per entity of given dim and codim is a constant
    bool fixedsize (int dim, int codim) const
    {
      return(true);
    }

    /*! how many objects of type DataType have to be sent for a given entity
     *
     *  Note: Only the sender side needs to know this size. */
    template<class EntityType>
    size_t size (EntityType& e) const
    {
      return(1);
    }

    /*! pack data from user to message buffer */
    template<class MessageBuffer, class EntityType>
    void gather (MessageBuffer& buff, const EntityType& e) const
    {
      IdType id=globalidset_.id(e);

      buff.write((*mapid2entity_.find(id)).second);
    }

    /*! unpack data from message buffer to user

      n is the number of objects sent by the sender
    */
    template<class MessageBuffer, class EntityType>
    void scatter (MessageBuffer& buff, const EntityType& entity, size_t n)
    {
      ValueType x;
      buff.read(x);

      /** only if the incoming index is a valid one,
	  i.e. if it is greater than zero, will it be
	  inserted as the global index; it is made
	  sure in the upper class, i.e. GlobalUniqueIndex,
	  that non-owning processes use -1 to mark an entity
	  that they do not own.
      */
      if(x >= 0)
	{
	  IdType id = globalidset_.id(entity);
	  mapid2entity_.erase(id);
	  mapid2entity_.insert(std::make_pair(id,x));

	  const int lindex = indexSet_.index(entity);
	  localGlobalMap_[lindex] = x;
	  globalLocalMap_[x]      = lindex;
	}
    }

    //! constructor
    IndexExchange (const GlobalIdSet& globalidset, Map& mapid2entity, int &rank,
		   const typename GridView::IndexSet& localIndexSet, IndexMap& localGlobal, IndexMap& globalLocal)
      : globalidset_(globalidset),
	mapid2entity_(mapid2entity),
	rank_(rank),
	indexSet_(localIndexSet),
	localGlobalMap_(localGlobal),
	globalLocalMap_(globalLocal)
    {
    }

  private:
    const GlobalIdSet& globalidset_;
    Map& mapid2entity_;
    int& rank_;

    const typename GridView::IndexSet& indexSet_;
    IndexMap& localGlobalMap_;
    IndexMap& globalLocalMap_;
  };

public:
  /**********************************************************************************************************************/
  /* \brief class constructor - when the class is instantiated by passing a const reference to a GridView object,       */
  /*                      we calculate the complete set of global unique indices so that we can then                    */
  /*  later query the global index, by directly passing the entity in question: then the respective global              */
  /*  index is returned.                                                                                                */
  /**********************************************************************************************************************/
  GlobalUniqueIndex(const GridView& gridview)
    :       gridview_(gridview),
	    uniqueEntityPartition_(gridview)
  {
    int rank = gridview.comm().rank();
    int size = gridview.comm().size();

    nLocalEntity_=0;
    nGlobalEntity_=0;

    nLocalEntity_ = uniqueEntityPartition_.numOwners();


    /** compute the global, non-redundant number of entities, i.e. the number of entities in the set
     *  without double, aka. redundant entities, on the interprocessor boundary via global reduce. */

    const CollectiveCommunication& collective = gridview_.comm();
    nGlobalEntity_ = collective.template sum<int>(nLocalEntity_);


    /** communicate the number of locally owned entities to all other processes so that the respective offset
     *  can be calculated on the respective processor; we use the Dune mpi collective communication facility
     *  for this; first, we gather the number of locally owned entities on the root process and, second, we
     *  broadcast the array to all processes where the respective offset can be calculated. */

    std::vector<int> offset(size);
    std::fill(offset.begin(), offset.end(), 0);

    /** gather number of locally owned entities on root process */
    collective.template gather<int>(&nLocalEntity_,offset.data(),1,0);

    /** broadcast the array containing the number of locally owned entities to all processes */
    collective.template broadcast<int>(offset.data(),size,0);

    indexOffset_.clear();
    indexOffset_.resize(size,0);

    for(unsigned int ii=0; ii<indexOffset_.size(); ++ii)
	for(unsigned int jj=0; jj < ii; ++jj)
	    indexOffset_[ii] += offset[jj];

    /** compute globally unique index over all processes; the idea of the algorithm is as follows: if
     *  an entity is owned by the process, it is assigned an index that is the addition of the offset
     *  specific for this process and a consecutively incremented counter; if the entity is not owned
     *  by the process, it is assigned -1, which signals that this specific entity will get its global
     *  unique index through communication afterwards;
     *
     *  thus, the calculation of the globally unique index is divided into 2 stages:
     *
     *  (1) we calculate the global index independently;
     *
     *  (2) we achieve parallel adjustment by communicating the index
     *      from the owning entity to the non-owning entity.
     *
     */

    /** 1st stage of global index calculation: calculate global index for owned entities */
     /** intialize map that stores an entity's global index via it's globally unique id as key */
    MapId2Index globalIndex;

    const typename GridView::IndexSet& indexSet = gridview.indexSet();

    const GlobalIdSet &globalIdSet=gridview_.grid().globalIdSet();  /** retrieve globally unique Id set */
    int myoffset=indexOffset_[rank];

    int globalcontrib=0;    /** initialize contribution for the global index */

    std::vector<bool> firstTime(gridview_.size(CODIM));
    std::fill(firstTime.begin(), firstTime.end(), true);

    for(Iterator iter = gridview_.template begin<0>();iter!=gridview_.template end<0>(); ++iter)
    {
#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
      for (size_t i=0; i<iter->subEntities(CODIM); i++)
#else
      for (size_t i=0; i<iter->template count<CODIM>(); i++)
#endif
      {
        IdType id=globalIdSet.subId(*iter,i,CODIM);                 /** retrieve the entity's id */

        int idx = gridview_.indexSet().subIndex(*iter,i,CODIM);

        if (! firstTime[idx] )
          continue;

        firstTime[idx] = false;

        if(uniqueEntityPartition_.owner(idx) == true)  /** if the entity is owned by the process, go ahead with computing the global index */
        {
          const int gindex = myoffset + globalcontrib;    /** compute global index */
          globalIndex.insert(std::make_pair(id,gindex)); /** insert pair (key, value) into the map */

          const int lindex = idx;
          localGlobalMap_[lindex] = gindex;
          globalLocalMap_[gindex] = lindex;

          globalcontrib++;                                /** increment contribution to global index */
        }
        else /** if entity is not owned, insert -1 to signal not yet calculated global index */
        {
          globalIndex.insert(std::make_pair(id,-1));
        }
      }

    }

    /** 2nd stage of global index calculation: communicate global index for non-owned entities */

    // Create the data handle and communicate.
    IndexExchange<GlobalIdSet,IdType,MapId2Index,int> dataHandle(globalIdSet,globalIndex,rank,indexSet,localGlobalMap_,globalLocalMap_);
    gridview_.communicate(dataHandle, Dune::All_All_Interface, Dune::ForwardCommunication);
  }


  /** \brief Given a local index, retrieve its index globally unique over all processes. */
  int globalIndex(const int& localIndex) const {
    return localGlobalMap_.find(localIndex)->second;
  }

  int localIndex(const int& globalIndex) const {
    return globalLocalMap_.find(globalIndex)->second;
  }

  int globalIndex(const typename GridView::template Codim<CODIM>::Entity& entity) const
  {
    return localGlobalMap_.find(gridview_.indexSet().index(entity))->second;
  }

  int localIndex(const typename GridView::template Codim<CODIM>::Entity& entity) const {
    return gridview_.indexSet().index(entity);
  }

  inline unsigned int nGlobalEntity() const
  {
    return(nGlobalEntity_);
  }

  inline unsigned int nOwnedLocalEntity() const
  {
    return(nLocalEntity_);
  }

  const GridView& getGridView() const {
    return gridview_;
  }

protected:
  /** store data members */
  const GridView gridview_;                                                       /** store a const reference to a gridview */
  UniqueEntityPartitionType uniqueEntityPartition_;
  int nLocalEntity_;                                                      /** store number of entities that are owned by the local */
  int nGlobalEntity_;                                             /** store global number of entities, i.e. number of entities without rendundant entities on interprocessor boundaries */
  std::vector<int> indexOffset_;                          /** store offset of entity index on every process */

  IndexMap localGlobalMap_;
  IndexMap globalLocalMap_;
};

#endif /* GLOBALUNIQUEINDEX_HH_ */