subdomaingrid.hh

#ifndef DUNE_MULTIDOMAINGRID_SUBDOMAINGRID_HH
#define DUNE_MULTIDOMAINGRID_SUBDOMAINGRID_HH
 
#include <memory>
#include <string>
 
#include <dune/common/exceptions.hh>
 
#include <dune/grid/multidomaingrid/hostgridaccessor.hh>
#include <dune/grid/multidomaingrid/subdomainset.hh>
 
#include <dune/grid/multidomaingrid/subdomaingrid/geometry.hh>
#include <dune/grid/multidomaingrid/subdomaingrid/localgeometry.hh>
#include <dune/grid/multidomaingrid/subdomaingrid/entity.hh>
#include <dune/grid/multidomaingrid/subdomaingrid/iterator.hh>
#include <dune/grid/multidomaingrid/subdomaingrid/hierarchiciterator.hh>
#include <dune/grid/multidomaingrid/subdomaingrid/intersection.hh>
#include <dune/grid/multidomaingrid/subdomaingrid/intersectioniterator.hh>
#include <dune/grid/multidomaingrid/subdomaingrid/idsets.hh>
#include <dune/grid/multidomaingrid/subdomaingrid/indexsets.hh>
#include <dune/grid/multidomaingrid/subdomaingrid/gridview.hh>
 
 
namespace Dune {
 
namespace mdgrid {
 
 
// forward declaration in correct namespace
template<typename,typename>
class MultiDomainGrid;
 
namespace subdomain {
 
template<typename MDGrid>
class SubDomainGrid;
 
template<typename MDGrid>
struct SubDomainGridFamily {
 
private:
 
  static const int dim  = MDGrid::dimension;
  static const int dimw = MDGrid::dimensionworld;
 
public:
 
  struct Traits
  {
    using Grid = SubDomainGrid<MDGrid>;
 
 
    using LeafIntersection = Dune::Intersection<
      const Grid,
      IntersectionWrapper<
        const Grid,
        IndexSetWrapper<
          const Grid,
          typename MDGrid::LeafIndexSetImp
          >,
        typename MDGrid::LeafGridView::Intersection
        >
      >;
 
    using LevelIntersection = Dune::Intersection<
      const Grid,
      IntersectionWrapper<
        const Grid,
        IndexSetWrapper<
          const Grid,
          typename MDGrid::LevelIndexSetImp
          >,
        typename MDGrid::LevelGridView::Intersection
        >
      >;
 
    using LeafIntersectionIterator = Dune::IntersectionIterator<
      const Grid,
      IntersectionIteratorWrapper<
        const Grid,
        IndexSetWrapper<
          const Grid,
          typename MDGrid::LeafIndexSetImp
          >,
        typename MDGrid::LeafGridView::IntersectionIterator
        >,
      IntersectionWrapper<
        const Grid,
        IndexSetWrapper<
          const Grid,
          typename MDGrid::LeafIndexSetImp
          >,
        typename MDGrid::LeafGridView::Intersection
        >
      >;
 
    using LevelIntersectionIterator = Dune::IntersectionIterator<
      const Grid,
      IntersectionIteratorWrapper<
        const Grid,
        IndexSetWrapper<
          const Grid,
          typename MDGrid::LevelIndexSetImp
          >,
        typename MDGrid::LevelGridView::IntersectionIterator
        >,
      IntersectionWrapper<
        const Grid,
        IndexSetWrapper<
          const Grid,
          typename MDGrid::LevelIndexSetImp
          >,
        typename MDGrid::LevelGridView::Intersection
        >
      >;
 
    using HierarchicIterator = Dune::EntityIterator< 0, const Grid, HierarchicIteratorWrapper< const Grid > >;
 
    using LevelGridView = Dune::GridView<LevelGridViewTraits<const Grid>>;
    using LeafGridView = Dune::GridView<LeafGridViewTraits<const Grid>>;
 
    template <int cd>
    struct Codim
    {
 
      using Geometry      = Dune::Geometry<dim-cd, dimw, const Grid, GeometryWrapper>;
      using LocalGeometry = Dune::Geometry<dim-cd, dim, const Grid, LocalGeometryWrapper>;
 
      using Entity        = Dune::Entity<cd, dim, const Grid, EntityWrapper>;
 
      using EntitySeed    = EntitySeedWrapper<typename MDGrid::HostGrid::template Codim<cd>::EntitySeed>;
 
      template <PartitionIteratorType pitype>
      struct Partition
      {
 
        using LevelIterator = Dune::EntityIterator<
          cd,
          const Grid,
          IteratorWrapper<
            typename Traits::LevelGridView,
            typename MDGrid::LevelGridView::template Codim<cd>::template Partition<pitype>::Iterator,
            cd,
            pitype,
            const Grid
            >
          >;
 
        using LeafIterator = Dune::EntityIterator<
          cd,
          const Grid,
          IteratorWrapper<
            typename Traits::LeafGridView,
            typename MDGrid::LeafGridView::template Codim<cd>::template Partition<pitype>::Iterator,
            cd,
            pitype,
            const Grid
            >
          >;
      };
 
      using LeafIterator  = typename Partition< All_Partition >::LeafIterator;
      using LevelIterator = typename Partition< All_Partition >::LevelIterator;
 
    private:
      friend class Dune::Entity<cd, dim, const Grid, EntityWrapper>;
 
    };
 
    using LevelIndexSet = IndexSetWrapper<const Grid, typename MDGrid::LevelIndexSetImp>;
    using LeafIndexSet  = IndexSetWrapper<const Grid, typename MDGrid::LeafIndexSetImp>;
 
    using GlobalIdSet = IdSet<
      const Grid,
      IdSetWrapper<
        const Grid,
        typename MDGrid::HostGrid::Traits::GlobalIdSet
        >,
      typename MDGrid::HostGrid::Traits::GlobalIdSet::IdType
      >;
 
    using LocalIdSet = IdSet<
      const Grid,
      IdSetWrapper<
        const Grid,
        typename MDGrid::HostGrid::Traits::LocalIdSet
        >,
      typename MDGrid::HostGrid::Traits::LocalIdSet::IdType
      >;
 
    using Communication = typename MDGrid::Communication;
  };
 
};
 
 
template<typename MDGrid>
class SubDomainGrid :
    public GridDefaultImplementation<MDGrid::dimension,
                                     MDGrid::dimensionworld,
                                     typename MDGrid::ctype,
                                     SubDomainGridFamily<MDGrid> > {
 
  template<typename,typename>
  friend class ::Dune::mdgrid::MultiDomainGrid;
 
  template<int codim, int dim, typename GridImp>
  friend class EntityWrapperBase;
 
  template<int codim, int dim, typename GridImp>
  friend class EntityWrapper;
 
  template<typename, typename, int codim, PartitionIteratorType pitype, typename GridImp>
  friend class IteratorWrapper;
 
  template<typename GridImp>
  friend class HierarchicIteratorWrapper;
 
  template<int mydim, int coorddim, typename GridImp>
  friend class GeometryWrapper;
 
  template<int mydim, int coorddim, typename GridImp>
  friend class LocalGeometryWrapper;
 
  template<typename GridImp, typename WrappedIndexSet>
  friend class IndexSetWrapper;
 
  template<typename GridImp, typename WrappedIdSet>
  friend class IdSetWrapper;
 
  template<typename GridImp>
  friend struct ::Dune::mdgrid::detail::HostGridAccessor;
 
  template<typename,typename,typename>
  friend class IntersectionIteratorWrapper;
 
  template<typename,typename,typename>
  friend class IntersectionWrapper;
 
  template<typename>
  friend class LevelGridView;
 
  template<typename>
  friend class LeafGridView;
 
  typedef GridDefaultImplementation<MDGrid::dimension,
                                    MDGrid::dimensionworld,
                                    typename MDGrid::ctype,
                                    SubDomainGridFamily<MDGrid>
                                    > BaseT;
 
  using GridImp = SubDomainGrid<MDGrid>;
 
public:
 
  using MultiDomainGrid = MDGrid;
 
  using HostGrid = typename MDGrid::HostGrid;
 
private:
 
  typedef IndexSetWrapper<const SubDomainGrid<MDGrid>, typename MDGrid::LevelIndexSetImp> LevelIndexSetImp;
 
  typedef IndexSetWrapper<const SubDomainGrid<MDGrid>, typename MDGrid::LeafIndexSetImp> LeafIndexSetImp;
 
  typedef IdSetWrapper<const SubDomainGrid<MDGrid>, typename HostGrid::Traits::GlobalIdSet> GlobalIdSetImp;
 
  typedef IdSetWrapper<const SubDomainGrid<MDGrid>, typename HostGrid::Traits::LocalIdSet> LocalIdSetImp;
 
public:
 
  typedef SubDomainGridFamily<MDGrid> GridFamily;
  typedef typename GridFamily::Traits Traits;
 
  typedef typename MDGrid::ctype ctype;
 
  typedef typename MDGrid::SubDomainIndex SubDomainIndex;
 
  enum IntersectionType { neighbor, foreign, boundary, processor };
 
  using BaseT::dimension;
  using BaseT::dimensionworld;
 
 
  template<typename EntitySeed>
  typename Traits::template Codim<EntitySeed::codimension>::Entity
  entity(const EntitySeed& entitySeed) const
  {
    return
      EntityWrapper<EntitySeed::codimension,dimension,const GridImp>(
        this,
        typename MDGrid::template Codim<EntitySeed::codimension>::Entity(
          _grid.entity(entitySeed)
        )
      );
  }
 
  int maxLevel() const {
    return _grid.maxLevel();
  }
 
  template<int codim>
  typename Traits::template Codim<codim>::LevelIterator lbegin(int level) const {
    return IteratorWrapper<
      typename Traits::LevelGridView,
      typename MultiDomainGrid::LevelGridView::template Codim<codim>::template Partition<All_Partition>::Iterator,
      codim,
      All_Partition,
      const GridImp>(
                     this,
                     &this->levelGridView(level).indexSet(),
                     _grid.levelGridView(level).template begin<codim>(),
                     _grid.levelGridView(level).template end<codim>()
                     );
  }
 
  template<int codim>
  typename Traits::template Codim<codim>::LevelIterator lend(int level) const {
    return IteratorWrapper<
      typename Traits::LevelGridView,
      typename MultiDomainGrid::LevelGridView::template Codim<codim>::template Partition<All_Partition>::Iterator,
      codim,
      All_Partition,
      const GridImp>(
                     this,
                     &this->levelGridView(level).indexSet(),
                     _grid.levelGridView(level).template end<codim>(),
                     _grid.levelGridView(level).template end<codim>()
                     );
  }
 
  template<int codim, PartitionIteratorType pitype>
  typename Traits::template Codim<codim>::template Partition<pitype>::LevelIterator lbegin(int level) const {
    return IteratorWrapper<
      typename Traits::LevelGridView,
      typename MultiDomainGrid::LevelGridView::template Codim<codim>::template Partition<pitype>::Iterator,
      codim,
      pitype,
      const GridImp>(
                     this,
                     &this->levelGridView(level).indexSet(),
                     _grid.levelGridView(level).template begin<codim,pitype>(),
                     _grid.levelGridView(level).template end<codim,pitype>()
                     );
  }
 
  template<int codim, PartitionIteratorType pitype>
  typename Traits::template Codim<codim>::template Partition<pitype>::LevelIterator lend(int level) const {
    return IteratorWrapper<
      typename Traits::LevelGridView,
      typename MultiDomainGrid::LevelGridView::template Codim<codim>::template Partition<pitype>::Iterator,
      codim,
      pitype,
      const GridImp>(
                     this,
                     &this->levelGridView(level).indexSet(),
                     _grid.levelGridView(level).template end<codim,pitype>(),
                     _grid.levelGridView(level).template end<codim,pitype>()
                     );
  }
 
  template<int codim>
  typename Traits::template Codim<codim>::LeafIterator leafbegin() const {
    return IteratorWrapper<
      typename Traits::LeafGridView,
      typename MultiDomainGrid::LeafGridView::template Codim<codim>::template Partition<All_Partition>::Iterator,
      codim,
      All_Partition,
      const GridImp>(
                     this,
                     &this->leafGridView().indexSet(),
                     _grid.leafGridView().template begin<codim>(),
                     _grid.leafGridView().template end<codim>()
                     );
  }
 
  template<int codim>
  typename Traits::template Codim<codim>::LeafIterator leafend() const {
    return IteratorWrapper<
      typename Traits::LeafGridView,
      typename MultiDomainGrid::LeafGridView::template Codim<codim>::template Partition<All_Partition>::Iterator,
      codim,
      All_Partition,
      const GridImp>(
                     this,
                     &this->leafGridView().indexSet(),
                     _grid.leafGridView().template end<codim>(),
                     _grid.leafGridView().template end<codim>()
                     );
  }
 
  template<int codim, PartitionIteratorType pitype>
  typename Traits::template Codim<codim>::template Partition<pitype>::LeafIterator leafbegin() const {
    return IteratorWrapper<
      typename Traits::LeafGridView,
      typename MultiDomainGrid::LeafGridView::template Codim<codim>::template Partition<pitype>::Iterator,
      codim,
      pitype,
      const GridImp>(
                     this,
                     &this->leafGridView().indexSet(),
                     _grid.leafGridView().template begin<codim,pitype>(),
                     _grid.leafGridView().template end<codim,pitype>()
                     );
  }
 
  template<int codim, PartitionIteratorType pitype>
  typename Traits::template Codim<codim>::template Partition<pitype>::LeafIterator leafend() const {
    return IteratorWrapper<
      typename Traits::LeafGridView,
      typename MultiDomainGrid::LeafGridView::template Codim<codim>::template Partition<pitype>::Iterator,
      codim,
      pitype,
      const GridImp>(
                     this,
                     &this->leafGridView().indexSet(),
                     _grid.leafGridView().template end<codim,pitype>(),
                     _grid.leafGridView().template end<codim,pitype>()
                     );
  }
 
  int size(int level, int codim) const {
    assert(level <= maxLevel());
    return _levelIndexSets[level]->size(codim);
  }
 
  int size(int codim) const {
    return _leafIndexSet.size(codim);
  }
 
  int size(int level, GeometryType type) const {
    assert(level <= maxLevel());
    return _levelIndexSets[level]->size(type);
  }
 
  int size(GeometryType type) const {
    // TODO: check this (perhaps return sum over levelindexsets?)
    return _leafIndexSet.size(type);
  }
 
  const typename Traits::GlobalIdSet& globalIdSet() const {
    return _globalIdSet;
  }
 
  const typename Traits::LocalIdSet& localIdSet() const {
    return _localIdSet;
  }
 
  const typename Traits::LevelIndexSet& levelIndexSet(int level) const {
    if (!_grid.supportLevelIndexSets()) {
      DUNE_THROW(GridError,"level index set support not enabled for this grid");
    }
    assert(level <= maxLevel());
    assert(_levelIndexSets[level]);
    return *_levelIndexSets[level];
  }
 
  const typename Traits::LeafIndexSet& leafIndexSet() const {
    return _leafIndexSet;
  }
 
 
  void globalRefine(int refCount) {
    DUNE_THROW(NotImplemented,"grid modification only allowed on the MultiDomainGrid");
  }
 
 
  bool mark(int refCount, const typename Traits::template Codim<0>::Entity& e) {
    return _grid.mark(refCount,multiDomainEntity(e));
  }
 
 
  int getMark(const typename Traits::template Codim<0>::Entity& e) {
    return _grid.getMark(multiDomainEntity(e));
  }
 
 
  bool preAdapt() {
    DUNE_THROW(NotImplemented,"grid modification only allowed on the MultiDomainGrid");
  }
 
 
  bool adapt() {
    DUNE_THROW(NotImplemented,"grid modification only allowed on the MultiDomainGrid");
  }
 
 
  void postAdapt() {
    DUNE_THROW(NotImplemented,"grid modification only allowed on the MultiDomainGrid");
  }
 
  int overlapSize(int level, int codim) const {
    return _grid.overlapSize(level,codim);
  }
 
  int overlapSize(int codim) const {
    return _grid.overlapSize(codim);
  }
 
  int ghostSize(int level, int codim) const {
    return _grid.ghostSize(level,codim);
  }
 
  int ghostSize(int codim) const {
    return _grid.ghostSize(codim);
  }
 
  const typename Traits::Communication& comm() const {
    return _grid.comm();
  }
 
  template<typename DataHandleImp, typename DataTypeImp>
  void communicate (CommDataHandleIF<DataHandleImp,DataTypeImp> &data,
                    InterfaceType iftype,
                    CommunicationDirection dir,
                    int level) const
  {
    DataHandleWrapper<CommDataHandleIF<DataHandleImp,DataTypeImp> > datahandle(data,*this);
    _grid.hostGrid().levelGridView(level).communicate(datahandle,iftype,dir);
  }
 
  template<typename DataHandleImp, typename DataTypeImp>
  void communicate (CommDataHandleIF<DataHandleImp,DataTypeImp> &data,
                    InterfaceType iftype,
                    CommunicationDirection dir) const
  {
    DataHandleWrapper<CommDataHandleIF<DataHandleImp,DataTypeImp> > datahandle(data,*this);
    _grid.hostGrid().leafGridView().communicate(datahandle,iftype,dir);
  }
 
  size_t numBoundarySegments() const
  {
    return _grid.numBoundarySegments();
  }
 
  const MDGrid& multiDomainGrid() const {
    return _grid;
  }
 
  SubDomainIndex domain() const {
    return _subDomain;
  }
 
  void update() {
    if (_grid.supportLevelIndexSets()) {
      while (_levelIndexSets.size() <= static_cast<std::size_t>(maxLevel())) {
        _levelIndexSets.push_back(std::make_shared<LevelIndexSetImp>(*this,_grid.levelIndexSet(_levelIndexSets.size())));
      }
    }
  }
 
  bool operator==(const SubDomainGrid& rhs) const {
    return (&_grid == &rhs._grid && _subDomain == rhs._subDomain);
  }
 
 
  template <
      class EntityType,
      std::enable_if_t<std::is_same_v<typename MDGrid::Traits::template Codim<EntityType::codimension>::Entity,EntityType>, int> = 0>
  typename Traits::template Codim<EntityType::codimension>::Entity subDomainEntity(const EntityType& mdEntity) const {
    return EntityWrapper<EntityType::codimension,dimension,const GridImp>(this,mdEntity);
  }
 
  template <
      class EntityType,
      std::enable_if_t<std::is_same_v<typename MDGrid::Traits::template Codim<EntityType::codimension>::Entity,EntityType>, int> = 0>
  static const EntityType&
  multiDomainEntity(const EntityType &e) {
    return e;
  }
 
  template <
      class EntityType,
      std::enable_if_t<!std::is_same_v<typename MDGrid::Traits::template Codim<EntityType::codimension>::Entity,EntityType>,int> = 0>
  static const typename MDGrid::Traits::template Codim<EntityType::codimension>::Entity &
  multiDomainEntity(const EntityType &e) {
    return e.impl().multiDomainEntity();
  }
 
  template<typename EntityType>
  static const typename MDGrid::template HostEntity<EntityType>::type& hostEntity(const EntityType& e) {
    return e.impl().hostEntity();
  }
 
  static const auto& multiDomainIntersection(const typename Traits::LeafIntersection& is) {
    return is.impl().multiDomainIntersection();
  }
 
  static const auto& multiDomainIntersection(const typename Traits::LevelIntersection& is) {
    return is.impl().multiDomainIntersection();
  }
 
  typename Traits::LeafIntersectionIterator subDomainIntersectionIterator(const typename MDGrid::LeafSubDomainInterfaceIterator it) const {
    assert(_subDomain == it->subDomain1() || _subDomain == it->subDomain2());
    if (_subDomain == it->subDomain1())
      return IntersectionIteratorWrapper<
        const GridImp,
        typename GridImp::LeafGridView::IndexSet,
        typename MDGrid::LeafGridView::IntersectionIterator
        >(
          &this->leafGridView().indexSet(),
          it->firstMultiDomainIntersectionIterator()
          );
    else
      return IntersectionIteratorWrapper<
        const GridImp,
        typename GridImp::LeafGridView::IndexSet,
        typename MDGrid::LeafGridView::IntersectionIterator
        >(
          &this->leafGridView().indexSet(),
          it->secondMultiDomainIntersectionIterator()
          );
  }
 
  typename Traits::LevelIntersectionIterator subDomainIntersectionIterator(const typename MDGrid::LevelSubDomainInterfaceIterator it) const {
    assert(_subDomain == it->subDomain1() || _subDomain == it->subDomain2());
    if (_subDomain == it->subDomain1())
      return IntersectionIteratorWrapper<
        const GridImp,
        typename GridImp::LevelGridView::IndexSet,
        typename MDGrid::LevelGridView::IntersectionIterator
        >(
          &this->levelGridView(it->firstMultiDomainIntersectionIterator()->inside().level()).indexSet(),
          it->firstMultiDomainIntersectionIterator()
          );
    else
      return IntersectionIteratorWrapper<
        const GridImp,
        typename GridImp::LevelGridView::IndexSet,
        typename MDGrid::LevelGridView::IntersectionIterator
        >(
          &this->levelGridView(it->secondMultiDomainIntersectionIterator()->inside().level()).indexSet(),
          it->secondMultiDomainIntersectionIterator()
          );
  }
 
  template<typename Intersection>
  IntersectionType intersectionType(const Intersection& intersection) const {
    return intersection.impl().intersectionType();
  }
 
private:
 
  MDGrid& _grid;
  SubDomainIndex _subDomain;
  GlobalIdSetImp _globalIdSet;
  LocalIdSetImp _localIdSet;
  LeafIndexSetImp _leafIndexSet;
  std::vector<std::shared_ptr<LevelIndexSetImp> > _levelIndexSets;
 
  SubDomainGrid(MDGrid& grid, SubDomainIndex subDomain) :
    _grid(grid),
    _subDomain(subDomain),
    _globalIdSet(*this,grid.hostGrid().globalIdSet()),
    _localIdSet(*this,grid.hostGrid().localIdSet()),
    _leafIndexSet(*this,grid.leafIndexSet())
  {
    update();
  }
 
  template<typename EntityType>
  bool containsMultiDomainEntity(const EntityType& e) const {
    if (_grid.supportLevelIndexSets())
      return levelIndexSet(e.level()).containsMultiDomainEntity(e);
    else
      return leafIndexSet().containsMultiDomainEntity(e);
  }
 
  template<typename EntityType>
  bool containsHostEntity(const EntityType& e) const {
    if (_grid.supportLevelIndexSets())
      return levelIndexSet(e.level()).containsHostEntity(e);
    else
      return leafIndexSet().containsHostEntity(e);
  }
 
  SubDomainGrid(const SubDomainGrid& rv);
  SubDomainGrid& operator=(const SubDomainGrid& rv);
 
 
  template<typename Impl>
  struct DataHandleWrapper
    : public Dune::CommDataHandleIF<DataHandleWrapper<Impl>,
                                    typename Impl::DataType
                                    >
  {
 
    bool contains(int dim, int codim) const
    {
      return _impl.contains(dim,codim); // TODO: check if codim supported
    }
 
    bool fixedSize(int dim, int codim) const
    {
      //_impl.fixedSize(dim,codim); // TODO: warning if true?
      return false;
    }
 
    template<typename Entity>
    std::size_t size(const Entity& e) const
    {
      if (_grid.containsHostEntity(e))
        return _impl.size(_grid.subDomainEntity(_grid._grid.wrapHostEntity(e)));
      else
        return 0;
    }
 
    template<typename MessageBufferImp, typename Entity>
    void gather(MessageBufferImp& buf, const Entity& e) const
    {
      if (_grid.containsHostEntity(e))
        _impl.gather(buf,_grid.subDomainEntity(_grid._grid.wrapHostEntity(e)));
    }
 
    template<typename MessageBufferImp, typename Entity>
    void scatter(MessageBufferImp& buf, const Entity& e, std::size_t n)
    {
      if (_grid.containsHostEntity(e))
        _impl.scatter(buf,_grid.subDomainEntity(_grid._grid.wrapHostEntity(e)),n);
    }
 
    DataHandleWrapper(Impl& impl, const SubDomainGrid<MDGrid>& grid)
      : _impl(impl)
      , _grid(grid)
    {}
 
    Impl& _impl;
    const SubDomainGrid<MDGrid>& _grid;
 
  };
 
};
 
} // namespace subdomain
 
} // namespace mdgrid
 
} // namespace Dune
 
#endif // DUNE_MULTIDOMAINGRID_MULTIDOMAINGRID_HH