globalp2mapper.hh

#ifndef DUNE_GFE_PARALLEL_GLOBALP2MAPPER_HH
#define DUNE_GFE_PARALLEL_GLOBALP2MAPPER_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 Dune header files
#include <dune/common/version.hh>

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

namespace Dune {

  template <class GridView>
  class GlobalP2Mapper
  {
  public:

    /** \brief The integer number type used for indices */
    typedef typename GridView::IndexSet::IndexType Index;

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

    GlobalP2Mapper(const GridView& gridView)
    : p2Mapper_(gridView)
    {
      static_assert(GridView::dimension==2, "Only implemented for two-dimensional grids");

      GlobalIndexSet<GridView> globalVertexIndex(gridView,2);
      GlobalIndexSet<GridView> globalEdgeIndex(gridView,1);
      GlobalIndexSet<GridView> globalElementIndex(gridView,0);

      // total number of degrees of freedom
      size_ = globalVertexIndex.size(2) + globalEdgeIndex.size(1) + globalElementIndex.size(0);

      // Determine
      for (auto it = gridView.template begin<0>(); it != gridView.template end<0>(); ++it)
      {
        // Loop over all vertices
#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
        for (size_t i=0; i<it->subEntities(2); i++)
#else
        for (size_t i=0; i<it->template count<2>(); i++)
#endif
        {
          //int localIndex  = globalVertexIndex.localIndex (*it->template subEntity<2>(i));
          int localIndex  = p2Mapper_.map(*it, i, 2);
          int globalIndex = globalVertexIndex.index(*it->template subEntity<2>(i));

          localGlobalMap_[localIndex]  = globalIndex;
          globalLocalMap_[globalIndex] = localIndex;
        }

        // Loop over all edges
#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
        for (size_t i=0; i<it->subEntities(1); i++)
#else
        for (size_t i=0; i<it->template count<1>(); i++)
#endif
        {
          //int localIndex  = globalEdgeIndex.localIndex (*it->template subEntity<1>(i)) + gridView.size(2);
          int localIndex  = p2Mapper_.map(*it, i, 1);
          int globalIndex = globalEdgeIndex.index(*it->template subEntity<1>(i)) + globalVertexIndex.size(2);

          localGlobalMap_[localIndex]  = globalIndex;
          globalLocalMap_[globalIndex] = localIndex;
        }

        // One element degree of freedom for quadrilaterals
        if (not it->type().isQuadrilateral())
          DUNE_THROW(Dune::NotImplemented, "for non-quadrilaterals");

        if (it->type().isQuadrilateral())
        {
          //int localIndex  = globalEdgeIndex.localIndex (*it->template subEntity<1>(i)) + gridView.size(2);
          int localIndex  = p2Mapper_.map(*it, 0, 0);
          int globalIndex = globalElementIndex.index(*it->template subEntity<0>(0))
                            + globalEdgeIndex.size(1)
                            + globalVertexIndex.size(2);

          localGlobalMap_[localIndex]  = globalIndex;
          globalLocalMap_[globalIndex] = localIndex;
        }

      }

    }

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

    template <class Entity>
    Index subIndex(const Entity& entity, uint i, uint codim) const
    {
      int localIndex = p2Mapper_.map(entity, i, codim);
      return localGlobalMap_.find(localIndex)->second;
    }

    template <class Entity>
    bool contains(const Entity& entity, uint i, uint codim, Index& result) const
    {
      Index localIndex;
      if (not p2Mapper_.contains(entity, i, codim,localIndex))
        return false;
      result = localGlobalMap_.find(localIndex)->second;
      return true;
    }

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

    unsigned int size() const
    {
      return size_;
    }

    P2BasisMapper<GridView> p2Mapper_;

    IndexMap localGlobalMap_;
    IndexMap globalLocalMap_;

    size_t nOwnedLocalEntity_;
    size_t size_;

  };

}
#endif /* GLOBALUNIQUEINDEX_HH_ */