multidomainmcmgmapper.hh

#ifndef DUNE_MULTIDOMAINGRID_MULTDIDOMAINMCMGMAPPER_HH
#define DUNE_MULTIDOMAINGRID_MULTDIDOMAINMCMGMAPPER_HH
 
#include <iostream>
#include <map>
#include <numeric>
 
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/geometry/referenceelements.hh>
 
namespace Dune {
 
namespace mdgrid {
 
namespace {
 
  template<typename GV, bool max_subomain_index_is_static>
  struct MCMGMapperStorageProvider
  {
 
  protected:
 
    typedef typename GV::IndexSet::IndexType IndexType;
 
    void update(const GV& gv) {}
 
    std::array<std::vector<IndexType>,GV::Grid::maxSubDomainIndex()> _offsets; // provide a map with all geometry types
 
  };
 
  template<typename GV>
  struct MCMGMapperStorageProvider<GV,false>
  {
 
  protected:
 
    typedef typename GV::IndexSet::IndexType IndexType;
 
    void update(const GV& gv) {
      _offsets.resize(gv.grid().maxSubDomainIndex());
    }
 
    std::vector<std::vector<IndexType> > _offsets; // provide a map with all geometry types
 
  };
 
 
}
 
template <typename GV>
class MultiDomainMCMGMapper : public MultipleCodimMultipleGeomTypeMapper<GV>,
                              public MCMGMapperStorageProvider<GV,GV::Grid::maxSubDomainIndexIsStatic()>
 
{
 
  typedef MultipleCodimMultipleGeomTypeMapper<GV> Base;
 
  typedef MCMGMapperStorageProvider<GV,GV::Grid::maxSubDomainIndexIsStatic()> StorageProvider;
 
  using StorageProvider::_offsets;
 
public:
 
  typedef typename GV::IndexSet::IndexType IndexType;
  typedef typename GV::Grid::SubDomainIndex SubDomainIndex;
 
  MultiDomainMCMGMapper (const GV& gv, const MCMGLayout& layout)
    : Base(gv,layout)
  {
    update(gv);
  }
 
  MultiDomainMCMGMapper (const GV& gv)
    : Base(gv)
  {
    update(gv);
  }
 
  using Base::gridView;
 
  template<class EntityType>
  int map (SubDomainIndex subDomain, const EntityType& e) const
  {
    return gridView().indexSet().index(subDomain,e) + _offsets[subDomain][GlobalGeometryTypeIndex::index(e.type())];
  }
 
  int map (SubDomainIndex subDomain, const typename GV::template Codim<0>::Entity& e, int i, unsigned int codim) const
  {
    GeometryType gt =
        ReferenceElements<double, GV::dimension>::general(e.type()).type(i,
                                                                         codim);
    return gridView().indexSet().subIndex(subDomain, e, i, codim) +
           _offsets[subDomain][GlobalGeometryTypeIndex::index(gt)];
  }
 
  int size (SubDomainIndex subDomain) const
  {
    return _offsets[subDomain].back();
  }
 
  template<class EntityType>
  bool contains (SubDomainIndex subDomain, const EntityType& e, IndexType& result) const
  {
    if (!gridView().indexSet().contains(subDomain, e) ||
        !Base::layout()(e.type(), GV::dimension)) {
      result = 0;
      return false;
    }
    result = map(subDomain,e);
    return true;
  }
 
  template<int cc> // this is now the subentity's codim
  bool contains (SubDomainIndex subDomain, const typename GV::template Codim<0>::Entity& e, int i, IndexType& result) const
  {
    GeometryType gt = ReferenceElements<double,GV::dimension>::general(e.type()).type(i,cc);
    // if the entity is contained in the subdomain, all of its subentities are contained as well
    if (!gridView().indexSet().contains(subDomain, e) ||
        !Base::layout()(gt, GV::dimension)) {
      result = 0;
      return false;
    }
    result = gridView().indexSet().subIndex(subDomain, e, i, cc) +
             _offsets[subDomain][GlobalGeometryTypeIndex::index(gt)];
    return true;
  }
 
 
  void update(const GV& gv)
  {
    static_cast<Base*>(this)->update(gv);
    StorageProvider::update(gv);
    for (SubDomainIndex subDomain = 0;
         subDomain < gridView().grid().maxSubDomainIndex(); ++subDomain) {
      std::vector<IndexType>& offsets = _offsets[subDomain];
      offsets.resize(GlobalGeometryTypeIndex::size(GV::dimension) + 1);
      std::fill(offsets.begin(),offsets.end(),0);
 
      // Compute offsets for the different geometry types.
      // Note that mapper becomes invalid when the grid is modified.
      for (int cc = 0; cc <= GV::dimension; ++cc)
        {
          for (auto gt : gridView().indexSet().types(subDomain,cc))
            offsets[GlobalGeometryTypeIndex::index(gt) + 1] =
                gridView().indexSet().size(subDomain, gt);
          // convert sizes to offset
          // last entry stores total size
          std::partial_sum(offsets.begin(),offsets.end(),offsets.begin());
        }
    }
  }
};
 
} // namespace mdgrid
 
} // namespace Dune
 
#endif // DUNE_MULTIDOMAINGRID_MULTDIDOMAINMCMGMAPPER_HH