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 <dune/functions/functionspacebases/lagrangebasis.hh>

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);

      auto localView = p2Mapper_.localView();
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
      auto localIndexSet = p2Mapper_.localIndexSet();
#endif

      // Determine
      for (const auto& element : elements(gridView))
      {
        localView.bind(element);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
        localIndexSet.bind(localView);
#endif

        // Loop over all local degrees of freedom
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
        for (size_t i=0; i<localIndexSet.size(); i++)
#else
        for (size_t i=0; i<localView.size(); i++)
#endif
        {
          int codim = localView.tree().finiteElement().localCoefficients().localKey(i).codim();
          int entity   = localView.tree().finiteElement().localCoefficients().localKey(i).subEntity();

#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
          int localIndex  = localIndexSet.index(i);
#else
          auto localIndex  = localView.index(i);
#endif
          int globalIndex;
          switch (codim)
          {
            case 2:  // vertex dofs
              globalIndex = globalVertexIndex.index(element.template subEntity<2>(entity));
              break;

            case 1:  // edge dofs
              globalIndex = globalEdgeIndex.index(element.template subEntity<1>(entity)) + globalVertexIndex.size(2);
              break;

            case 0:  // element dofs
              globalIndex = globalElementIndex.index(element.template subEntity<0>(entity))
                            + globalEdgeIndex.size(1)
                            + globalVertexIndex.size(2);
              break;

            default:
              DUNE_THROW(Dune::Exception, "Impossible codimension!");
          }

          localGlobalMap_[localIndex]  = globalIndex;
        }

      }

    }

    /** \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 subEntity, uint codim, Index& result) const
    {
      auto localView = p2Mapper_.localView();
      localView.bind(entity);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
      auto localIndexSet = p2Mapper_.localIndexSet();
      localIndexSet.bind(localView);
#endif

      Index localIndex;
      bool dofFound = false;
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
      for (size_t i=0; i<localIndexSet.size(); i++)
#else
      for (size_t i=0; i<localView.size(); i++)
#endif
      {
        if (localView.tree().finiteElement().localCoefficients().localKey(i).subEntity() == subEntity
          and localView.tree().finiteElement().localCoefficients().localKey(i).codim() == codim)
        {
          dofFound = true;
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
          localIndex = localIndexSet.index(i);
#else
          localIndex = localView.index(i);
#endif
          break;
        }
      }

      if (not dofFound)
        return false;

      result = localGlobalMap_.find(localIndex)->second;
      return true;
    }

    unsigned int size() const
    {
      return size_;
    }

    Functions::LagrangeBasis<GridView,2> p2Mapper_;

    IndexMap localGlobalMap_;

    size_t nOwnedLocalEntity_;
    size_t size_;

  };

}
#endif   // DUNE_GFE_PARALLEL_GLOBALP2MAPPER_HH