grid.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_MMESH_INTERFACE_GRID_HH
#define DUNE_MMESH_INTERFACE_GRID_HH
 
#include <string>
#include <unordered_set>
#include <unordered_map>
#include <memory>
 
#include <dune/common/deprecated.hh>
#include <dune/common/parallel/communication.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/version.hh>
#include <dune/grid/common/capabilities.hh>
#include <dune/grid/common/grid.hh>
#include <dune/mmesh/grid/multiid.hh>
 
// The components of the MMesh interface
#include <dune/mmesh/grid/mmesh.hh>
#include "entity.hh"
#include "entityseed.hh"
#include "connectedcomponent.hh"
#include "dgfparser.hh"
#include "geometry.hh"
#include "gridfactory.hh"
#include "hierarchiciterator.hh"
#include "indexsets.hh"
#include "intersectioniterator.hh"
#include "intersections.hh"
#include "leafiterator.hh"
#include "traits.hh"
 
namespace Dune
{
 
#if HAVE_MPI
  using Comm = MPI_Comm;
#else
  using Comm = No_Comm;
#endif
 
  // MMeshInterfaceGrid family
  template<class MMesh>
  struct MMeshInterfaceGridFamily
  {
 
  public:
 
    typedef GridTraits<
        MMesh::dimension-1, // grid dimension
        MMesh::dimension, // world dimension
        MMeshInterfaceGrid<MMesh>,
        MMeshInterfaceGridGeometry,
        MMeshInterfaceGridEntity,
        MMeshInterfaceGridLeafIterator, // LevelIterator
        MMeshInterfaceGridLeafIntersection,
        MMeshInterfaceGridLeafIntersection, // LevelIntersection
        MMeshInterfaceGridLeafIntersectionIterator,
        MMeshInterfaceGridLeafIntersectionIterator, // LevelIntersectionIterator
        MMeshInterfaceGridHierarchicIterator,
        MMeshInterfaceGridLeafIterator,
        MMeshInterfaceGridLeafIndexSet< const MMeshInterfaceGrid<MMesh> >, // LevelIndexSet
        MMeshInterfaceGridLeafIndexSet< const MMeshInterfaceGrid<MMesh> >,
        MMeshInterfaceGridGlobalIdSet< const MMeshInterfaceGrid<MMesh> >,
        MMeshImpl::MultiId, // GlobalIdSet::IdType,
        MMeshInterfaceGridGlobalIdSet< const MMeshInterfaceGrid<MMesh> >, // LocalIdSet
        MMeshImpl::MultiId, // LocalIdSet::IdType,
        Communication< Comm >,
        DefaultLevelGridViewTraits,
        DefaultLeafGridViewTraits,
        MMeshInterfaceGridEntitySeed
        > Traits;
  };
 
  //**********************************************************************
  //
  // --MMeshInterfaceGrid class
  //
  //************************************************************************
  template <class MMesh>
  class MMeshInterfaceGrid
#ifndef DOXYGEN_SHOULD_SKIP_THIS
   : public GridDefaultImplementation< MMesh::dimension-1, MMesh::dimension,
                                       typename MMesh::FieldType,
                                       MMeshInterfaceGridFamily<MMesh> >
#endif /* DOXYGEN_SHOULD_SKIP_THIS */
  {
  public:
    static constexpr int dimension = MMesh::dimension-1;
    static constexpr int dimensionworld = MMesh::dimension;
 
    using FieldType = typename MMesh::FieldType;
    using IdType = typename MMesh::IdType;
 
    using BoundarySegments = std::unordered_map< IdType, std::size_t >;
 
    //**********************************************************
    // The Interface Methods
    //**********************************************************
 
    typedef MMeshInterfaceGridFamily<MMesh> GridFamily;
    typedef typename GridFamily::Traits Traits;
 
    using GridImp = typename GridFamily::Traits::Grid;
 
    typedef std::unique_ptr< GridImp > GridPtrType;
 
    using HostGridType = typename MMesh::HostGridType;
 
    using MMeshType = MMesh;
 
    typedef typename Traits::template Codim<0>::LeafIterator LeafIterator;
 
    typedef FieldType ctype;
 
    typedef Dune::FieldVector<ctype, dimensionworld> GlobalCoordinate;
 
    template<int cd>
    using MMeshInterfaceEntity = typename HostGridEntityChooser_<HostGridType, dimensionworld, cd+1>::type;
 
    using VertexHandle = MMeshInterfaceEntity<dimension>;
 
    using EdgeHandle = MMeshInterfaceEntity<dimension-1>;
 
    using ElementOutput = std::list<MMeshInterfaceEntity<0>>;
 
    using Entity = typename Traits::template Codim<0>::Entity;
 
    using Vertex = typename Traits::template Codim<dimension>::Entity;
 
    using ConnectedComponent = MMeshInterfaceConnectedComponent<const GridImp>;
 
    using RemeshingIndicator = RatioIndicator<GridImp>;
 
    explicit MMeshInterfaceGrid(MMesh* mMesh, BoundarySegments boundarySegments = {})
     : mMesh_(mMesh),
       boundarySegments_(boundarySegments),
       #ifdef HAVE_MPI
       comm_( MPIHelper::getCommunicator() )
      #endif
    {
      leafIndexSet_ = std::make_unique<MMeshInterfaceGridLeafIndexSet<const GridImp>>( this );
      globalIdSet_ = std::make_unique<MMeshInterfaceGridGlobalIdSet<const GridImp>>( this );
    }
 
    int maxLevel() const {
      return 0;
    }
 
    int size (int level, int codim) const {
      // we only have one level
      assert(level == 0);
      return size(codim);
    }
 
    size_t numBoundarySegments () const {
      return localBoundarySegments_.size();
    }
 
    const BoundarySegments& boundarySegments() const
    {
      return localBoundarySegments_;
    }
 
    void addBoundarySegment ( const std::vector< std::size_t >& ids, std::size_t bndSegIdx )
    {
      boundarySegments_[ ids ] = bndSegIdx;
    }
 
    void setBoundarySegments( const BoundarySegments boundarySegments )
    {
      boundarySegments_ = boundarySegments;
    }
 
    int size (int codim) const {
      return leafIndexSet().size(codim);
    }
 
 
    int size (int level, GeometryType type) const {
      // we only have one level
      assert(level == 0);
      return size(type);
    }
 
 
    int size (GeometryType type) const
    {
      return leafIndexSet().size(type);
    }
 
    const typename Traits::GlobalIdSet& globalIdSet() const {
      return *globalIdSet_;
    }
 
 
    const typename Traits::LocalIdSet& localIdSet() const {
      return *globalIdSet_;
    }
 
 
    const MMeshInterfaceGridLeafIndexSet<const GridImp>& levelIndexSet(int level) const
    {
      if (level != 0)
        DUNE_THROW(GridError, "levelIndexSet of nonexisting level " << level << " requested!");
      return *leafIndexSet_;
    }
 
    const MMeshInterfaceGridLeafIndexSet<const GridImp>& leafIndexSet() const
    {
      return *leafIndexSet_;
    }
 
    template < class EntitySeed >
    typename Traits::template Codim<EntitySeed::codimension>::Entity
    entity(const EntitySeed& seed) const
    {
      typedef MMeshInterfaceGridEntity<
        EntitySeed::codimension,
        dimension,
        const typename Traits::Grid
        > EntityImp;
 
      auto hostEntity = seed.impl().hostEntity();
      assert( hostEntity != decltype(hostEntity)() );
      return EntityImp(this, hostEntity);
    }
 
    typename Traits::template Codim<dimension>::Entity entity(const MMeshInterfaceEntity<dimension>& hostEntity) const
    {
      return entity( typename Traits::template Codim<dimension>::EntitySeed( hostEntity ) );
    }
 
    typename Traits::template Codim<0>::Entity entity(const MMeshInterfaceEntity<0>& hostEntity) const
    {
      return entity( typename Traits::template Codim<0>::EntitySeed( hostEntity ) );
    }
 
    template<int dimw = dimension+1>
    std::enable_if_t<dimw == 3, typename Traits::template Codim<1>::Entity> entity(const MMeshInterfaceEntity<1>& hostEntity) const
    {
      return entity( typename Traits::template Codim<1>::EntitySeed( hostEntity ) );
    }
 
    template<int codim>
    typename Traits::template Codim<codim>::LevelIterator lbegin (int level) const {
      return MMeshInterfaceGridLeafIterator<codim,All_Partition, const GridImp>(this);
    }
 
 
    template<int codim>
    typename Traits::template Codim<codim>::LevelIterator lend (int level) const {
      return MMeshInterfaceGridLeafIterator<codim,All_Partition, const GridImp>(this, true);
    }
 
 
    template<int codim, PartitionIteratorType PiType>
    typename Traits::template Codim<codim>::template Partition<PiType>::LevelIterator lbegin (int level) const {
      return MMeshInterfaceGridLeafIterator<codim,PiType, const GridImp>(this);
    }
 
 
    template<int codim, PartitionIteratorType PiType>
    typename Traits::template Codim<codim>::template Partition<PiType>::LevelIterator lend (int level) const {
      return MMeshInterfaceGridLeafIterator<codim,PiType, const GridImp>(this, true);
    }
 
 
    template<int codim>
    typename Traits::template Codim<codim>::LeafIterator leafbegin() const {
      return MMeshInterfaceGridLeafIterator<codim,All_Partition, const GridImp>(this);
    }
 
 
    template<int codim>
    typename Traits::template Codim<codim>::LeafIterator leafend() const {
      return MMeshInterfaceGridLeafIterator<codim,All_Partition, const GridImp>(this, true);
    }
 
 
    template<int codim, PartitionIteratorType PiType>
    typename Traits::template Codim<codim>::template Partition<PiType>::LeafIterator leafbegin() const {
      return MMeshInterfaceGridLeafIterator<codim,PiType, const GridImp>(this);
    }
 
 
    template<int codim, PartitionIteratorType PiType>
    typename Traits::template Codim<codim>::template Partition<PiType>::LeafIterator leafend() const {
      return MMeshInterfaceGridLeafIterator<codim,PiType, const GridImp>(this, true);
    }
 
    void globalRefine(int steps = 1)
    {
      for(int i = 0; i < steps; ++i)
      {
        // mark all elements
        for (const auto& element : elements(this->leafGridView()))
          mark( 1, element );
 
        preAdapt();
        adapt();
        postAdapt();
      }
    }
 
 
    bool mark(int refCount, const typename Traits::template Codim<0>::Entity & e) const
    {
      if ( refCount != 0 )
      {
        mark_[ globalIdSet_->id( e ) ] = refCount;
        return true;
      }
      return false;
    }
 
    int getMark(const typename Traits::template Codim<0>::Entity & e) const
    {
      auto entry = mark_.find( globalIdSet_->id( e ) );
      if ( entry != mark_.end() )
        return entry->second;
      else
        return 0;
    }
 
    bool preAdapt()
    {
      return mark_.size() > 0;
    }
 
    bool adapt()
    {
      return mMesh_->adapt();
    }
 
    template< class GridImp, class DataHandle >
    bool adapt ( AdaptDataHandleInterface< GridImp, DataHandle > &handle )
    {
      preAdapt();
      sequence_ += 1;
 
      for (const auto& ielement : elements( this->leafGridView() ))
        if (ielement.isNew())
        {
          bool initialize = true;
          const auto &connComp = ielement.impl().connectedComponent();
 
          for ( const auto& old : connComp.children() )
          {
            const Entity& father = old;
 
            if (father.geometry().volume() > ielement.geometry().volume())
            {
              ielement.impl().bindFather( father );
              handle.prolongLocal( father, ielement, true );
            }
            else
            {
              father.impl().bindFather( ielement );
              handle.restrictLocal( ielement, father, initialize );
            }
 
            initialize = false;
          }
        }
 
      postAdapt();
      return true;
    }
 
    void postAdapt()
    {
      mark_.clear();
      childrenConnectedComponentMap_.clear();
    }
 
  public:
    unsigned int overlapSize(int codim) const {
      return 0;
    }
 
 
    unsigned int ghostSize(int codim) const
    {
      std::size_t size = 0;
 
      if (codim == 0)
        for (const auto& e : elements(this->leafGridView()))
          if (e.partitionType() == GhostEntity)
            size++;
 
      if (codim == dimension)
        for (const auto& v : vertices(this->leafGridView()))
          if (v.partitionType() == GhostEntity)
            size++;
 
      return size;
    }
 
 
    unsigned int overlapSize(int level, int codim) const {
      return overlapSize(codim);
    }
 
 
    unsigned int ghostSize(int level, int codim) const {
      return ghostSize(codim);
    }
 
 
    void loadBalance() {
      mMesh_->loadBalance();
    };
 
    template<class T>
    bool loadBalance( const T& t ) { return false; };
 
    void loadBalance(int strategy, int minlevel, int depth, int maxlevel, int minelement){
      DUNE_THROW(NotImplemented, "MMeshInterfaceGrid::loadBalance()");
    }
 
 
    const Communication< Comm >& comm () const
    {
      return comm_;
    }
 
    template< class Data, class InterfaceType, class CommunicationDirection >
    void communicate (
      Data &dataHandle,
      InterfaceType interface,
      CommunicationDirection direction,
      int level = 0 ) const
    {
      if (comm().size() <= 1)
        return;
 
#if HAVE_MPI
      if( (interface == InteriorBorder_All_Interface) || (interface == All_All_Interface) )
      {
        MMeshCommunication<GridImp, MMeshType> communication( getMMesh().partitionHelper() );
        const auto& gv = this->leafGridView();
 
        switch( direction )
        {
        case ForwardCommunication:
          communication( gv.template begin< 0, Interior_Partition >(), gv.template end< 0, Interior_Partition >(),
                         gv.template begin< 0, Ghost_Partition >(), gv.template end< 0, Ghost_Partition >(),
                         dataHandle, InteriorEntity, GhostEntity,
                         (interface == All_All_Interface) );
          break;
 
        case BackwardCommunication:
          communication( gv.template begin< 0, Ghost_Partition >(), gv.template end< 0, Ghost_Partition >(),
                         gv.template begin< 0, Interior_Partition >(), gv.template end< 0, Interior_Partition >(),
                         dataHandle, GhostEntity, InteriorEntity,
                         (interface == All_All_Interface) );
          break;
        }
      }
      else
        DUNE_THROW( NotImplemented, "Communication on interface type " << interface << " not implemented." );
#else
        DUNE_THROW( NotImplemented, "MPI not found!" );
#endif //HAVE_MPI
    }
 
    bool isInterface( const MMeshInterfaceEntity<0>& segment ) const
    {
      int count = mMesh_->interfaceSegments().count( getVertexIds_( segment ) );
      assert( count <= 1 );
      return ( count > 0 );
    }
 
    template< int d = dimension >
    std::enable_if_t< d == 2, bool > isInterface( const MMeshInterfaceEntity<1>& edge ) const
    {
      auto circulator = mMesh_->getHostGrid().incident_facets( edge );
      for ( std::size_t i = 0; i < CGAL::circulator_size( circulator ); ++i, ++circulator )
        if ( isInterface( *circulator ) )
          return true;
      return false;
    }
 
    bool isInterface( const VertexHandle& vertex ) const
    {
      return vertex->info().isInterface;
    }
 
    bool isInterface( const typename Traits::LeafIntersection& intersection ) const
    {
      return isInterface( intersection.impl().getHostIntersection() );
    }
 
    template< class Entity >
    bool isInterface( const Entity& entity ) const
    {
      return isInterface( entity.impl().hostEntity() );
    }
 
    template< class Entity >
    std::size_t domainMarker( const Entity& entity ) const
    {
      assert( isInterface( entity ) );
      return mMesh_->interfaceSegments()[ getVertexIds_( entity.impl().hostEntity() ) ];
    }
 
    void markAsRefined( const std::vector< std::vector< std::size_t > >& children, const ConnectedComponent connectedComponent )
    {
      for( const auto& child : children )
        childrenConnectedComponentMap_.insert( { child, connectedComponent } );
    }
 
    bool hasConnectedComponent( const Entity& entity ) const
    {
      int count = childrenConnectedComponentMap_.count( getVertexIds_( entity.impl().hostEntity() ) );
      assert( count <= 1 );
      return (count == 1 );
    }
 
    const ConnectedComponent& getConnectedComponent( const Entity& entity ) const
    {
      auto it = childrenConnectedComponentMap_.find( getVertexIds_( entity.impl().hostEntity() ) );
      assert( it != childrenConnectedComponentMap_.end() );
      return it->second;
    }
 
    ConnectedComponent& getConnectedComponent( const Entity& entity )
    {
      auto it = childrenConnectedComponentMap_.find( getVertexIds_( entity.impl().hostEntity() ) );
      assert( it != childrenConnectedComponentMap_.end() );
      return it->second;
    }
 
    // Return if MMeshInterfaceEntity can be mirrored
    bool canBeMirrored(const MMeshInterfaceEntity<0>& hostEntity) const
    {
      using HostGridEntity = typename GridImp::MMeshType::template HostGridEntity<0>;
 
      if (hostEntity.first->neighbor(hostEntity.second) == HostGridEntity())
        return false;
 
      if constexpr (dimensionworld == 2)
        return hostEntity.first->dimension() >= 1;
 
      return true;
    }
 
    template <int d = dimensionworld>
    std::enable_if_t< d == 2, const MMeshInterfaceEntity<0> >
    mirrorHostEntity(const MMeshInterfaceEntity<0>& other) const
    {
      return mMesh_->getHostGrid().mirror_edge( other );
    }
 
    template <int d = dimensionworld>
    std::enable_if_t< d == 3, const MMeshInterfaceEntity<0> >
    mirrorHostEntity(const MMeshInterfaceEntity<0>& other) const
    {
      return mMesh_->getHostGrid().mirror_facet( other );
    }
 
    // **********************************************************
    // End of Interface Methods
    // **********************************************************
 
  private:
    Communication<Comm> comm_;
 
    static inline auto getVertexIds_( const MMeshInterfaceEntity<0>& entity )
    {
      std::vector<std::size_t> ids( dimensionworld );
      for( int i = 0; i < dimensionworld; ++i )
        ids[ i ] = entity.first->vertex((entity.second+i+1)%(dimensionworld+1))->info().id;
      std::sort(ids.begin(), ids.end());
      return ids;
    }
 
  public:
    void setIds()
    {
      globalIdSet_->update(this);
    }
 
    void setIndices()
    {
      leafIndexSet_->update(this);
 
      localBoundarySegments_.clear();
      std::size_t count = 0;
      for (const auto& e : elements(this->leafGridView()))
        for (const auto& is : intersections(this->leafGridView(), e))
          if (is.boundary())
          {
            auto iid = globalIdSet_->id( entity( is.impl().getHostIntersection() ) );
            localBoundarySegments_.insert( std::make_pair(iid, count++) );
          }
    }
 
    const MMesh& getMMesh() const
    {
      return *mMesh_;
    }
 
    const HostGridType& getHostGrid() const
    {
      return mMesh_->getHostGrid();
    }
 
    int sequence() const
    {
      return sequence_;
    }
 
    const auto& partitionHelper() const
    {
      return mMesh_->partitionHelper();
    }
 
  private:
    std::unique_ptr<MMeshInterfaceGridLeafIndexSet<const GridImp>> leafIndexSet_;
 
    std::unique_ptr<MMeshInterfaceGridGlobalIdSet<const GridImp>> globalIdSet_;
 
    std::unordered_map< std::vector< std::size_t >, ConnectedComponent, HashUIntVector > childrenConnectedComponentMap_;
 
  protected:
    MMesh* mMesh_;
 
    BoundarySegments boundarySegments_, localBoundarySegments_;
 
    mutable std::unordered_map< MMeshImpl::MultiId, int > mark_;
 
    int sequence_ = 0;
 
  }; // end Class MMeshInterfaceGrid
 
 
  // DGFGridInfo for MMeshInterfaceGrid
  // ----------------------------------
 
  template<class MMesh>
  struct DGFGridInfo< MMeshInterfaceGrid<MMesh> >
  {
    static int refineStepsForHalf ()
    {
      return MMesh::dimension-1;
    }
 
    static double refineWeight ()
    {
      return -1;
    }
  };
 
  // Capabilites of MMeshInterfaceGrid
  // ---------------------------------
 
  namespace Capabilities
  {
    template<class MMesh>
    struct hasSingleGeometryType<MMeshInterfaceGrid<MMesh>>
    {
      static const bool v = true;
      static const unsigned int topologyId = Dune::GeometryType::simplex;
    };
 
    template<class MMesh, int codim>
    struct hasEntity<MMeshInterfaceGrid<MMesh>, codim>
    {
      static const bool v = (codim >= 0 || codim <= MMesh::dimension-1);
    };
 
    template<class MMesh, int codim>
    struct hasEntityIterator<MMeshInterfaceGrid<MMesh>, codim>
    {
      static const bool v = (codim >= 0 || codim <= MMesh::dimension-1);
    };
 
    template<class MMesh>
    struct isLevelwiseConforming<MMeshInterfaceGrid<MMesh>>
    {
      static const bool v = true;
    };
 
    template<class MMesh>
    struct isLeafwiseConforming<MMeshInterfaceGrid<MMesh>>
    {
      static const bool v = true;
    };
  } // end namespace Capabilities
 
} // namespace Dune
 
#endif // DUNE_MMESH_INTERFACE_GRID_HH