globalindexset.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

#ifndef DUNE_GRID_UTILITY_GLOBALINDEXSET_HH
#define DUNE_GRID_UTILITY_GLOBALINDEXSET_HH

#include <vector>
#include <iostream>
#include <fstream>
#include <memory>
#include <map>
#include <utility>
#include <algorithm>

#include <dune/grid/common/gridenums.hh>
#include <dune/grid/common/datahandleif.hh>

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

namespace Dune
{

  template<class GridView>
  class GlobalIndexSet
  {
  public:
    typedef int Index;

    template <class Entity, int Codim>
    struct SubPartitionTypeProvider
    {
      static PartitionType get(const Entity& entity, int codim, int i)
      {
        if (codim==Codim)
          return entity.template subEntity<Codim>(i).partitionType();
        else
          return SubPartitionTypeProvider<Entity,Codim-1>::get(entity, codim, i);
      }
    };

    template <class Entity>
    struct SubPartitionTypeProvider<Entity,0>
    {
      static PartitionType get(const Entity& entity, int /*codim*/, int i)
      {
        return entity.template subEntity<0>(i).partitionType();
      }
    };

  private:
    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 Grid::Communication Communication;

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

    /*********************************************************************************************/
    /* calculate unique partitioning for all entities of a given codim in a given GridView,      */
    /* assuming they all have the same geometry, i.e. codim, type                                */
    /*********************************************************************************************/
    class UniqueEntityPartition
    {
    private:
      /* A DataHandle class to calculate the minimum of a std::vector which is accompanied by an index set */
      template<class IS, class V> // mapper type and vector type
      class MinimumExchange
      : public Dune::CommDataHandleIF<MinimumExchange<IS,V>,typename V::value_type>
      {
      public:
        typedef typename V::value_type DataType;

        bool contains (int /*dim*/, unsigned int codim) const
        {
          return codim==indexSetCodim_;
        }

        bool fixedSize (int /*dim*/, int /*codim*/) const
        {
          return true ;
        }

        template<class EntityType>
        std::size_t size (EntityType&) const
        {
          return 1 ;
        }

        template<class MessageBuffer, class EntityType>
        void gather (MessageBuffer& buff, const EntityType& e) const
        {
          buff.write(v_[indexset_.index(e)]);
        }

        template<class MessageBuffer, class EntityType>
        void scatter (MessageBuffer& buff, const EntityType& e, std::size_t)
        {
          DataType x;
          buff.read(x);
          if (x>=0) // other is -1 means, he does not want it
            v_[indexset_.index(e)] = std::min(x,v_[indexset_.index(e)]);
        }

        MinimumExchange (const IS& indexset, V& v, unsigned int indexSetCodim)
        : indexset_(indexset),
          v_(v),
          indexSetCodim_(indexSetCodim)
        {}

      private:
        const IS& indexset_;
        V& v_;
        unsigned int indexSetCodim_;
      };

    public:
      UniqueEntityPartition (const GridView& gridview, unsigned int codim)
      : assignment_(gridview.size(codim))
      {
        typedef typename GridView::IndexSet IndexSet;

        // assign own rank to entities that I might have
        for (auto it = gridview.template begin<0>(); it!=gridview.template end<0>(); ++it)
          for (unsigned int i=0; i<it->subEntities(codim); i++)
          {
            // Evil hack: I need to call subEntity, which needs the entity codimension as a static parameter.
            // However, we only have it as a run-time parameter.
            PartitionType subPartitionType = SubPartitionTypeProvider<typename GridView::template Codim<0>::Entity, GridView::dimension>::get(*it,codim,i);

            assignment_[gridview.indexSet().subIndex(*it,i,codim)]
              = ( subPartitionType==Dune::InteriorEntity or subPartitionType==Dune::BorderEntity )
              ? gridview.comm().rank()  // set to own rank
              : - 1;   // it is a ghost entity, I will not possibly own it.
          }

        MinimumExchange<IndexSet,std::vector<Index> > dh(gridview.indexSet(),assignment_,codim);

        gridview.communicate(dh,Dune::All_All_Interface,Dune::ForwardCommunication);
      }

      int owner(std::size_t i)
      {
        return assignment_[i];
      }

      std::size_t numOwners(int rank) const
      {
        return std::count(assignment_.begin(), assignment_.end(), rank);
      }

    private:
      std::vector<int> assignment_;
    };

  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.
     */
    class IndexExchange
    : public Dune::CommDataHandleIF<IndexExchange,Index>
    {
    public:
      bool contains (int /*dim*/, unsigned int codim) const
      {
        return codim==indexSetCodim_;
      }

      bool fixedSize (int /*dim*/, int /*codim*/) const
      {
        return true;
      }

      template<class EntityType>
      std::size_t size (EntityType&) const
      {
        return 1;
      }

      template<class MessageBuffer, class EntityType>
      void gather (MessageBuffer& buff, const EntityType& e) const
      {
        IdType id=globalidset_.id(e);

        if (indexSetCodim_==0)
          buff.write(mapid2entity_[id]);
        else
          buff.write((*mapid2entity_.find(id)).second);
      }

      template<class MessageBuffer, class EntityType>
      void scatter (MessageBuffer& buff, const EntityType& entity, std::size_t)
      {
        Index x;
        buff.read(x);

        if(x >= 0) {
          const IdType id = globalidset_.id(entity);

          if (indexSetCodim_==0)
            mapid2entity_[id] = x;
          else
          {
            mapid2entity_.erase(id);
            mapid2entity_.insert(std::make_pair(id,x));

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

      IndexExchange (const GlobalIdSet& globalidset, MapId2Index& mapid2entity,
                     const typename GridView::IndexSet& localIndexSet, IndexMap& localGlobal,
                     unsigned int indexSetCodim)
      : globalidset_(globalidset),
      mapid2entity_(mapid2entity),
      indexSet_(localIndexSet),
      localGlobalMap_(localGlobal),
      indexSetCodim_(indexSetCodim)
      {}

    private:
      const GlobalIdSet& globalidset_;
      MapId2Index& mapid2entity_;

      const typename GridView::IndexSet& indexSet_;
      IndexMap& localGlobalMap_;
      unsigned int indexSetCodim_;
    };

  public:
    GlobalIndexSet(const GridView& gridview, int codim)
    : gridview_(gridview),
      codim_(codim)
    {
      int rank = gridview.comm().rank();
      int size = gridview.comm().size();

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

      std::unique_ptr<UniqueEntityPartition> uniqueEntityPartition;
      if (codim_!=0)
        uniqueEntityPartition = std::make_unique<UniqueEntityPartition>(gridview,codim_);

      int nLocalEntity = (codim_==0)
                    ? std::distance(gridview.template begin<0, Dune::Interior_Partition>(), gridview.template end<0, Dune::Interior_Partition>())
                    : uniqueEntityPartition->numOwners(rank);

      // 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. */
      nGlobalEntity_ = gridview.comm().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 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);

      gridview_.comm().template allgather<int>(&nLocalEntity, 1, offset.data());

      int myoffset = 0;
      for (int i=1; i<rank+1; i++)
        myoffset += offset[i-1];

      /*  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
      // initialize map that stores an entity's global index via it's globally unique id as key
      globalIndex_.clear();

      const GlobalIdSet& globalIdSet = gridview_.grid().globalIdSet();      
      Index globalcontrib = 0;      
      if (codim_==0)  // This case is simpler
      {
        for (Iterator iter = gridview_.template begin<0>(); iter!=gridview_.template end<0>(); ++iter)
        {
          const IdType id = globalIdSet.id(*iter);      
          if (iter->partitionType() == Dune::InteriorEntity)
          {
            const Index gindex = myoffset + globalcontrib;    
            globalIndex_[id] = gindex;                      
            globalcontrib++;                                
          }

          else
          {
            globalIndex_[id] = -1;     
          }
        }
      }
      else  // if (codim==0) else
      {
      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)
      {
        for (std::size_t i=0; i<iter->subEntities(codim_); i++)
        {
          IdType id=globalIdSet.subId(*iter,i,codim_);

          Index idx = gridview_.indexSet().subIndex(*iter,i,codim_);

          if (!firstTime[idx] )
            continue;

          firstTime[idx] = false;

          if (uniqueEntityPartition->owner(idx) == rank)  
          {
            const Index gindex = myoffset + globalcontrib;    
            globalIndex_.insert(std::make_pair(id,gindex)); 
            const Index lindex = idx;
            localGlobalMap_[lindex] = gindex;

            globalcontrib++;                                
          }
          else 
          {
            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 dataHandle(globalIdSet,globalIndex_,indexSet,localGlobalMap_,codim_);
      gridview_.communicate(dataHandle, Dune::All_All_Interface, Dune::ForwardCommunication);
    }

    template <class Entity>
    Index index(const Entity& entity) const
    {
      if (codim_==0)
      {
        const GlobalIdSet& globalIdSet = gridview_.grid().globalIdSet(); 
        const IdType id = globalIdSet.id(entity);                        
        const Index gindex = globalIndex_.find(id)->second;                
        return gindex;
      }
      else
        return localGlobalMap_.find(gridview_.indexSet().index(entity))->second;
    }

    template <class Entity>
    Index subIndex(const Entity& entity, unsigned int i, unsigned int codim) const
    {
      if (codim_==0)
      {
        const GlobalIdSet& globalIdSet = gridview_.grid().globalIdSet(); 
        const IdType id = globalIdSet.subId(entity,i,codim);                        
        const Index gindex = globalIndex_.find(id)->second;                
        return gindex;
      }
      else
        return localGlobalMap_.find(gridview_.indexSet().subIndex(entity,i,codim))->second;
    }

    unsigned int size(unsigned int codim) const
    {
      return (codim_==codim) ? nGlobalEntity_ : 0;
    }

  protected:
    const GridView gridview_;

    unsigned int codim_;

    int nGlobalEntity_;

    IndexMap localGlobalMap_;

    MapId2Index globalIndex_;
  };

}  // namespace Dune

#endif /* DUNE_GRID_UTILITY_GLOBALINDEXSET_HH */