entity.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_GRID_ENTITY_HH
#define DUNE_MMESH_GRID_ENTITY_HH
 
// Dune includes
#include <dune/grid/common/grid.hh>
 
namespace Dune
{
  // Forward declarations
  template<int codim, int dim, class GridImp>
  class MMeshEntity;
}
 
 
// MMesh includes
#include <dune/mmesh/grid/connectedcomponent.hh>
#include <dune/mmesh/grid/polygoncutting.hh>
#include <dune/mmesh/misc/partitionhelper.hh>
 
namespace Dune
{
  template<class GridImp>
  class MMeshLeafIntersectionIterator;
 
  template<class GridImp>
  class MMeshHierarchicIterator;
 
  // External forward declarations
  template< class Grid >
  struct HostGridAccess;
 
 
  //**********************************************************************
  //
  // --MMeshEntity
  // --Entity
  //
  template<int codim, int dim, class GridImp>
  class MMeshEntity
#ifndef DOXYGEN_SHOULD_SKIP_THIS
   : public EntityDefaultImplementation <codim,dim,GridImp,MMeshEntity>
#endif /* DOXYGEN_SHOULD_SKIP_THIS */
  {
    template <class GridImp_>
    friend class MMeshLeafIndexSet;
 
    template <class GridImp_>
    friend class MMeshLocalIdSet;
 
    template <class GridImp_>
    friend class MMeshGlobalIdSet;
 
    friend struct HostGridAccess< typename std::remove_const< GridImp >::type >;
 
    using IdType = MMeshImpl::MultiId;
 
  private:
    typedef typename GridImp::ctype ctype;
 
    // The equivalent entity in the host grid
    typedef typename GridImp::template HostGridEntity<codim> HostGridEntity;
 
    typedef typename GridImp::template Codim<codim>::Entity Entity;
 
  public:
    typedef typename GridImp::template Codim<codim>::Geometry Geometry;
 
    typedef typename GridImp::template Codim<codim>::EntitySeed EntitySeed;
 
    MMeshEntity()
      : mMesh_(nullptr)
    {}
 
    MMeshEntity(const GridImp* mMesh, const HostGridEntity& hostEntity)
      : hostEntity_(hostEntity)
      , mMesh_(mMesh)
    {}
 
    MMeshEntity(const GridImp* mMesh, HostGridEntity&& hostEntity)
      : hostEntity_(std::move(hostEntity))
      , mMesh_(mMesh)
    {}
 
    MMeshEntity(const MMeshEntity& original)
      : hostEntity_(original.hostEntity_)
      , mMesh_(original.mMesh_)
    {}
 
    MMeshEntity(MMeshEntity&& original)
      : hostEntity_(std::move(original.hostEntity_))
      , mMesh_(original.mMesh_)
    {}
 
    MMeshEntity& operator=(const MMeshEntity& original)
    {
      if (this != &original)
      {
        mMesh_ = original.mMesh_;
        hostEntity_ = original.hostEntity_;
      }
      return *this;
    }
 
    MMeshEntity& operator=(MMeshEntity&& original)
    {
      if (this != &original)
      {
        mMesh_ = original.mMesh_;
        hostEntity_ = std::move(original.hostEntity_);
      }
      return *this;
    }
 
    // Comparator for vertices
    template <int cc = codim>
    std::enable_if_t< cc == dim, bool >
    equals(const MMeshEntity& other) const
    {
      return hostEntity_ == other.hostEntity_;
    }
 
    // Comparator for edges
    template <int cc = codim>
    std::enable_if_t< cc == 1 && dim == 2, bool >
    equals(const MMeshEntity& other) const
    {
      return (hostEntity_ == other.hostEntity_)
       || ( mMesh_->getHostGrid().mirror_edge( hostEntity_ ) == other.hostEntity_ );
    }
 
    // Comparator for facets
    template <int cc = codim>
    std::enable_if_t< cc == 1 && dim == 3, bool >
    equals(const MMeshEntity& other) const
    {
      return (hostEntity_ == other.hostEntity_)
       || ( hostEntity_ == mMesh_->getHostGrid().mirror_facet( other.hostEntity_ ) );
    }
 
    // Comparator for edges in 3d
    template <int cc = codim>
    std::enable_if_t< cc == 2 && dim == 3, bool >
    equals(const MMeshEntity& other) const
    {
      const auto& v0 = hostEntity_.first->vertex( hostEntity_.second );
      const auto& v1 = hostEntity_.first->vertex( hostEntity_.third );
      const auto& w0 = other.hostEntity_.first->vertex( other.hostEntity_.second );
      const auto& w1 = other.hostEntity_.first->vertex( other.hostEntity_.third );
      return (v0 == w0 && v1 == w1) || (v0 == w1 && v1 == w0);
    }
 
    bool hasFather () const {
      return false;
    }
 
    EntitySeed seed () const
    {
      return EntitySeed( hostEntity_ );
    }
 
    int level () const {
      // we only have one level
      return 0;
    }
 
    PartitionType partitionType () const
    {
      const auto& ph = mMesh_->partitionHelper();
      if constexpr (codim == 0 || codim == dim)
        return ph.partitionType( hostEntity_->info().partition );
      else
        return ph.partitionType( grid().entity(hostEntity_) );
    }
 
    unsigned int subEntities (unsigned int cc) const
    {
      // we have a simplex grid
      int n = dim-codim+1;
      int k = dim-cc+1;
 
      // binomial: n over k
      int binomial=1;
      for (int i=n-k+1; i<=n; i++)
        binomial *= i;
      for (long i=2; i<=k; i++)
        binomial /= i;
 
      return binomial;
    }
 
    template <int cc>
    std::enable_if_t< codim == 1 && cc == dim, typename GridImp::template Codim<dim>::Entity >
    subEntity (unsigned int i) const
    {
      assert( i < subEntities( cc ) );
      const auto& cell = hostEntity_.first;
      auto facetIdx = MMeshImpl::cgalFacetToDuneFacet<dim, HostGridEntity>( hostEntity_ );
      const auto i0 = cgalIndex( MMeshImpl::ref<dim>().subEntity(facetIdx, 1, i, dim) );
 
      return MMeshEntity<cc, dim, GridImp> (
        mMesh_,
        typename GridImp::template HostGridEntity<dim> (
          cell->vertex( i0 )
        )
      );
    }
 
    template <int cc>
    std::enable_if_t< codim == 2 && cc == 3, typename GridImp::template Codim<3>::Entity >
    subEntity (unsigned int i) const
    {
      assert( i < subEntities( cc ) );
      const auto& cell = hostEntity_.first;
 
      auto edgeIdx = MMeshImpl::cgalEdgeToDuneEdge<3, HostGridEntity>( hostEntity_ );
 
      const auto i0 = cgalIndex( MMeshImpl::ref<dim>().subEntity(edgeIdx, 2, i, 3) );
 
      return MMeshEntity<cc, dim, GridImp> (
        mMesh_,
        typename GridImp::template HostGridEntity<dim> (
          cell->vertex( i0 )
        )
      );
    }
 
    template <int cc>
    std::enable_if_t< codim == dim && cc == dim, typename GridImp::template Codim<dim>::Entity >
    subEntity (unsigned int i) const
    {
      assert( i < subEntities( cc ) );
      return MMeshEntity<cc, dim, GridImp> (
        mMesh_,
        hostEntity_
      );
    }
 
    template <bool enable = true>
    std::enable_if_t< codim == dim && enable, MMeshIncidentIterator<GridImp> >
    incidentBegin () const {
      using Impl = typename MMeshIncidentIterator<GridImp>::Implementation;
      return MMeshIncidentIterator<GridImp>( Impl( mMesh_, hostEntity_) );
    }
 
    template <bool enable = true>
    std::enable_if_t< codim == dim && enable, MMeshIncidentIterator<GridImp> >
    incidentEnd () const {
      using Impl = typename MMeshIncidentIterator<GridImp>::Implementation;
      return MMeshIncidentIterator<GridImp>( Impl( mMesh_, hostEntity_, true ) );
    }
 
    template <bool enable = true>
    std::enable_if_t< codim == dim-1 && enable, MMeshEdgeIncidentIterator<GridImp> >
    incidentBegin () const {
      using Impl = typename MMeshEdgeIncidentIterator<GridImp>::Implementation;
      return MMeshEdgeIncidentIterator<GridImp>( Impl( mMesh_, hostEntity_) );
    }
 
    template <bool enable = true>
    std::enable_if_t< codim == dim-1 && enable, MMeshEdgeIncidentIterator<GridImp> >
    incidentEnd () const {
      using Impl = typename MMeshEdgeIncidentIterator<GridImp>::Implementation;
      return MMeshEdgeIncidentIterator<GridImp>( Impl( mMesh_, hostEntity_, true ) );
    }
 
    template <bool enable = true>
    std::enable_if_t< codim == dim && enable, MMeshIncidentFacetsIterator<GridImp> >
    incidentFacetsBegin () const {
      using Impl = typename MMeshIncidentFacetsIterator<GridImp>::Implementation;
      return MMeshIncidentFacetsIterator<GridImp>( Impl( mMesh_, hostEntity_) );
    }
 
    template <bool enable = true>
    std::enable_if_t< codim == dim && enable, MMeshIncidentFacetsIterator<GridImp> >
    incidentFacetsEnd () const {
      using Impl = typename MMeshIncidentFacetsIterator<GridImp>::Implementation;
      return MMeshIncidentFacetsIterator<GridImp>( Impl( mMesh_, hostEntity_, true ) );
    }
 
    template <bool enable = true>
    std::enable_if_t< codim == dim && enable, MMeshIncidentVerticesIterator<GridImp>>
    incidentVerticesBegin ( bool includeInfinite ) const {
      using Impl = typename MMeshIncidentVerticesIterator<GridImp>::Implementation;
      return MMeshIncidentVerticesIterator<GridImp>( Impl( mMesh_, hostEntity_, includeInfinite) );
    }
 
    template <bool enable = true>
    std::enable_if_t< codim == dim && enable, MMeshIncidentVerticesIterator<GridImp>>
    incidentVerticesEnd ( bool includeInfinite ) const {
      using Impl = typename MMeshIncidentVerticesIterator<GridImp>::Implementation;
      return MMeshIncidentVerticesIterator<GridImp>( Impl( mMesh_, hostEntity_, includeInfinite, true ) );
    }
 
    template <int cd = codim>
    std::enable_if_t< cd == dim, std::size_t >
    insertionLevel() const
    {
      if ( codim == dim )
        return hostEntity_->info().insertionLevel;
    }
 
    template <int cd = codim>
    std::enable_if_t< cd != dim, std::size_t >
    insertionLevel() const
    {
      std::size_t insertionLevel = 0;
      for( std::size_t i = 0; i < subEntities(dim); ++i )
        insertionLevel = std::max( insertionLevel, this->template subEntity<dim>(i).impl().insertionLevel() );
      return insertionLevel;
    }
 
    bool isInterface() const
    {
      if constexpr( codim == dim )
        return hostEntity_->info().isInterface;
      else
        DUNE_THROW( NotImplemented, "isInterface for codim != dim" );
    }
 
    int boundaryFlag() const
    {
      if constexpr( codim == dim )
        return hostEntity_->info().boundaryFlag;
      else
        DUNE_THROW( NotImplemented, "boudnaryFlag for codim != dim" );
    }
 
    Geometry geometry () const
    {
      return Geometry( hostEntity_ );
    }
 
    GeometryType type () const
    {
      return GeometryTypes::simplex(dim-codim);
    }
 
    const HostGridEntity& hostEntity () const
    {
      return hostEntity_;
    }
 
    HostGridEntity& hostEntity ()
    {
      return hostEntity_;
    }
 
    const GridImp& grid () const
    {
      return *mMesh_;
    }
 
    IdType id() const
    {
      // cache id
      if (id_ == IdType())
      {
        typename IdType::VT idlist( dim+1-codim );
        for( std::size_t i = 0; i < this->subEntities(dim); ++i )
          idlist[i] = this->template subEntity<dim>(i).impl().hostEntity()->info().id;
        std::sort( idlist.begin(), idlist.end() );
        id_ = IdType( idlist );
      }
 
      return id_;
    }
 
  private:
 
    const auto cgalIndex ( const std::size_t& i ) const
    {
      if constexpr( codim != dim )
        return hostEntity_.first->info().cgalIndex[i];
      else
      {
        DUNE_THROW( NotImplemented, "cgalIndex(i) for codim != dim");
        return 0;
      }
    }
 
 
    mutable IdType id_;
 
    HostGridEntity hostEntity_;
 
    const GridImp* mMesh_;
  };
 
 
  //***********************
  //
  //  --MMeshEntity
  //
  //***********************
  template<int dim, class GridImp>
  class MMeshEntity<0,dim,GridImp>
#ifndef DOXYGEN_SHOULD_SKIP_THIS
   : public EntityDefaultImplementation<0,dim,GridImp, MMeshEntity>
#endif /* DOXYGEN_SHOULD_SKIP_THIS */
  {
    friend struct HostGridAccess< typename std::remove_const< GridImp >::type >;
 
    typedef Entity< 0, dim, GridImp, MMeshEntity > EntityType;
 
  public:
    // equivalent entity in the host grid
    typedef typename GridImp::template HostGridEntity<0> HostGridEntity;
 
    typedef typename GridImp::template Codim<0>::Geometry Geometry;
 
    typedef typename GridImp::template Codim<0>::LocalGeometry LocalGeometry;
 
    // type of global coordinate
    typedef typename Geometry::GlobalCoordinate GlobalCoordinate;
 
    // ctype
    typedef typename GridImp::ctype ctype;
 
    // Grid
    typedef GridImp Grid;
 
    typedef MMeshLeafIntersectionIterator<GridImp> LeafIntersectionIterator;
 
    typedef MMeshHierarchicIterator<GridImp> HierarchicIterator;
 
    typedef typename GridImp::template Codim<0>::EntitySeed EntitySeed;
 
    using ConnectedComponent = MMeshConnectedComponent< GridImp >;
 
    using IdType = MMeshImpl::MultiId;
 
    using VertexStorage = std::array<GlobalCoordinate, dim+1>;
 
    MMeshEntity()
      : isLeaf_(true), mMesh_(nullptr)
    {}
 
    MMeshEntity(const GridImp* mMesh, const HostGridEntity& hostEntity)
      : hostEntity_(hostEntity), isLeaf_(true)
      , mMesh_(mMesh)
    {}
 
    MMeshEntity(const GridImp* mMesh, HostGridEntity&& hostEntity)
      : hostEntity_(std::move(hostEntity)), isLeaf_(true)
      , mMesh_(mMesh)
    {}
 
    MMeshEntity(const GridImp* mMesh, const HostGridEntity& hostEntity, const IdType& id)
      : hostEntity_(hostEntity)
      , id_(id), isLeaf_(false)
      , mMesh_(mMesh)
    {}
 
    MMeshEntity(const GridImp* mMesh, const VertexStorage& vertex)
      : id_( /*caching id*/ IdType({ std::size_t(-4), std::size_t(-3), std::size_t(-2) }) ), isLeaf_(false), mMesh_(mMesh), vertex_(vertex)
    {}
 
    MMeshEntity(const MMeshEntity& original)
      : hostEntity_(original.hostEntity_), id_(original.id_), isLeaf_(original.isLeaf_)
      , mMesh_(original.mMesh_), vertex_(original.vertex_)
    {}
 
    MMeshEntity(MMeshEntity&& original)
      : hostEntity_(std::move(original.hostEntity_)), id_(original.id_), isLeaf_(original.isLeaf_)
      , mMesh_(original.mMesh_), vertex_(original.vertex_)
    {}
 
    MMeshEntity& operator=(const MMeshEntity& original)
    {
      if (this != &original)
      {
        mMesh_ = original.mMesh_;
        hostEntity_ = original.hostEntity_;
        isLeaf_ = original.isLeaf_;
        id_ = original.id_;
        vertex_ = original.vertex_;
      }
      return *this;
    }
 
    MMeshEntity& operator=(MMeshEntity&& original)
    {
      if (this != &original)
      {
        mMesh_ = original.mMesh_;
        hostEntity_ = std::move(original.hostEntity_);
        isLeaf_ = original.isLeaf_;
        id_ = original.id_;
        vertex_ = original.vertex_;
      }
      return *this;
    }
 
    bool equals(const MMeshEntity& other) const
    {
      return hostEntity_ == other.hostEntity_;
    }
 
    bool operator==(const MMeshEntity& other) const
    {
      return this->equals(other);
    }
 
    bool operator<(const MMeshEntity& other) const
    {
      return hostEntity_ < other.hostEntity_;
    }
 
    MMeshEntity father () const
    {
      DUNE_THROW( InvalidStateException, "MMesh entities do no have a father, but a connectedComponent instead!" );
      return *this;
    }
 
    bool hasFather () const
    {
      return false;
    }
 
    /* \brief Returns the connected component
     */
    const ConnectedComponent& connectedComponent () const
    {
      assert( isNew() == true );
      return mMesh_->getConnectedComponent(*this);
    }
 
    const bool isNew () const
    {
      return hostEntity_->info().isNew;
    }
 
    void setIsNew ( bool isNew ) const
    {
      hostEntity_->info().isNew = isNew;
    }
 
    const bool mightVanish () const
    {
      return hostEntity_->info().mightVanish;
    }
 
    void setWillVanish ( bool mightVanish ) const
    {
      hostEntity_->info().mightVanish = mightVanish;
    }
 
    void mark ( int refCount ) const
    {
      hostEntity_->info().mark = refCount;
    }
 
    int getMark () const
    {
      return hostEntity_->info().mark;
    }
 
    EntitySeed seed () const
    {
      return EntitySeed(hostEntity_);
    }
 
    int level () const
    {
      // we only have one level
      return 0;
    }
 
    PartitionType partitionType () const
    {
      // for caching entities
      if (hostEntity_ == decltype(hostEntity_)())
        return InteriorEntity;
 
      const auto& ph = mMesh_->partitionHelper();
      return ph.partitionType( hostEntity_->info().partition );
    }
 
    Geometry geometry () const
    {
      if (isLeaf_)
        return Geometry( hostEntity_ );
      else
        return Geometry( this->vertex_ );
    }
 
    void bindFather( const EntityType& father )
    {
      father_ = &father;
    }
 
    LocalGeometry geometryInFather() const
    {
      if constexpr( dim != 2 )
        DUNE_THROW(NotImplemented, "geometryInFather() for dim != 2");
      else
      {
        assert( father_ != nullptr );
 
        auto thisPoints = this->vertex_;
 
        if( isLeaf_ )
          for ( int i = 0; i < 3; ++i )
            thisPoints[i] = geometry().corner(i);
 
        std::array< GlobalCoordinate, 3 > local;
        for ( int i = 0; i < 3; ++i )
          local[i] = father_->impl().geometry().local( thisPoints[i] );
 
        return LocalGeometry( local );
      }
    }
 
    unsigned int subEntities (unsigned int cc) const
    {
      return MMeshImpl::ref<dim>().size(cc);
    }
 
    template <int cc>
    std::enable_if_t< cc == 0, typename GridImp::template Codim<cc>::Entity >
    subEntity (unsigned int i) const {
      assert( i < subEntities( cc ) );
      return *this;
    }
 
    template <int cc>
    std::enable_if_t< cc == dim, typename GridImp::template Codim<cc>::Entity >
    subEntity (unsigned int i) const {
      assert( i < subEntities( cc ) );
      if (hostEntity_ != HostGridEntity())
      {
        const auto i0 = cgalIndex( i );
        return MMeshEntity<cc, dim, GridImp>( mMesh_, hostEntity_->vertex( i0 ) );
      }
      else {
        auto vh = mMesh_->getHostGrid().infinite_vertex();
        // TODO should distinguish for different i
        return MMeshEntity<cc, dim, GridImp>( mMesh_, vh );
      }
    }
 
    template <int cc>
    std::enable_if_t< cc == 1, typename GridImp::template Codim<cc>::Entity >
    subEntity (unsigned int i) const {
      assert( i < subEntities( cc ) );
      auto second = MMeshImpl::duneFacetToCgalSecond<dim>( i, cgalIndex() );
      return MMeshEntity<cc, dim, GridImp> ( mMesh_, typename GridImp::template HostGridEntity<1> ( hostEntity_, second ) );
    }
 
    template <int cc>
    std::enable_if_t< cc == 2 && dim == 3, typename GridImp::template Codim<cc>::Entity >
    subEntity (unsigned int i) const {
      assert( i < subEntities( cc ) );
 
      const auto i0 = cgalIndex( MMeshImpl::ref<dim>().subEntity(i, 2, 0, 3) );
      const auto i1 = cgalIndex( MMeshImpl::ref<dim>().subEntity(i, 2, 1, 3) );
 
      return MMeshEntity<cc, dim, GridImp> ( mMesh_, typename GridImp::template HostGridEntity<2> ( hostEntity_, i0, i1 ) );
    }
 
    MMeshLeafIntersectionIterator<GridImp> ileafbegin () const {
      return MMeshLeafIntersectionIterator<GridImp>(
        mMesh_,
        hostEntity_);
    }
 
    MMeshLeafIntersectionIterator<GridImp> ileafend () const {
      return MMeshLeafIntersectionIterator<GridImp>(
         mMesh_,
         hostEntity_,
         true);
    }
 
    MMeshLeafIntersectionIterator<GridImp> ilevelbegin () const {
      return ileafbegin();
    }
 
    MMeshLeafIntersectionIterator<GridImp> ilevelend () const {
      return ileafend();
    }
 
    MMeshHierarchicIterator<GridImp> hbegin (int maxlevel) const {
      return MMeshHierarchicIterator<GridImp>(
      mMesh_,
      *this,
      maxlevel);
    }
 
    MMeshHierarchicIterator<GridImp> hend (int maxlevel) const {
      return MMeshHierarchicIterator<GridImp>(
        mMesh_,
         *this,
         maxlevel,
         true);
    }
 
    bool isLeaf() const {
      return isLeaf_ && !isNew();
    }
 
    bool wasRefined () const
    {
      return false;
    }
 
    bool mightBeCoarsened () const
    {
      return false;
    }
 
    GeometryType type () const
    {
      return GeometryTypes::simplex(dim);
    }
 
    std::size_t domainMarker() const
    {
      return hostEntity_->info().domainMarker;
    }
 
    const HostGridEntity& hostEntity () const
    {
      return hostEntity_;
    }
 
    HostGridEntity& hostEntity ()
    {
      return hostEntity_;
    }
 
    const GridImp& grid () const
    {
      return *mMesh_;
    }
 
    IdType id() const
    {
      // cache id
      if (id_ == IdType())
      {
        typename IdType::VT idlist( dim+1 );
        for( std::size_t i = 0; i < this->subEntities(dim); ++i )
          idlist[i] = hostEntity_->vertex(i)->info().id;
        std::sort( idlist.begin(), idlist.end() );
        id_ = IdType( idlist );
      }
 
      return id_;
    }
 
  private:
 
    const auto& cgalIndex ( const std::size_t& i ) const
    {
      return hostEntity_->info().cgalIndex[i];
    }
 
    const auto& cgalIndex () const
    {
      return hostEntity_->info().cgalIndex;
    }
 
    HostGridEntity hostEntity_;
 
    mutable IdType id_;
 
    bool isLeaf_;
 
    const GridImp* mMesh_;
 
  protected:
    VertexStorage vertex_;
 
    const EntityType* father_;
 
  }; // end of MMeshEntity codim = 0
 
} // namespace Dune
 
#include <dune/mmesh/grid/cachingentity.hh>
 
#endif