intersections.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_INTERSECTIONS_HH
#define DUNE_MMESH_GRID_INTERSECTIONS_HH
 
#include <unordered_set>
 
// Dune MMesh includes
#include <dune/mmesh/grid/entity.hh>
#include <dune/mmesh/grid/leafiterator.hh>
#include <dune/mmesh/grid/pointfieldvector.hh>
#include <dune/mmesh/misc/twistutility.hh>
 
// local includes
 
namespace Dune
{
 
  // External forward declarations
  template< class Grid >
  struct HostGridAccess;
 
  template<class GridImp>
  class MMeshLeafIntersection
  {
    friend class MMeshLeafIntersectionIterator<GridImp>;
 
    friend struct HostGridAccess< typename std::remove_const< GridImp >::type >;
 
    enum {dim=GridImp::dimension};
 
    enum {dimworld=GridImp::dimensionworld};
 
    // The type used to store coordinates
    typedef typename GridImp::ctype ctype;
 
    typedef typename GridImp::template HostGridEntity<0> HostGridEntity;
    typedef typename GridImp::template HostGridEntity<1> HostLeafIntersection;
 
  public:
    enum {dimension=GridImp::dimension};
    enum {dimensionworld=GridImp::dimensionworld};
    typedef typename GridImp::template Codim<1>::Geometry Geometry;
    typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
    typedef typename GridImp::template Codim<1>::LocalGeometry::Implementation LocalGeometryImpl;
    typedef typename GridImp::template Codim<0>::Entity Entity;
    typedef FieldVector<ctype, dimworld> NormalVector;
 
    MMeshLeafIntersection()
    {}
 
    MMeshLeafIntersection(const GridImp* mMesh,
                          const HostGridEntity& hostEntity,
                          const int index)
      : mMesh_(mMesh)
      , hostIntersection_(hostEntity, index)
    {}
 
    MMeshLeafIntersection(const GridImp* mMesh,
                          HostLeafIntersection&& hostIntersection)
      : mMesh_(mMesh)
      , hostIntersection_(std::move(hostIntersection))
    {}
 
    bool equals(const MMeshLeafIntersection& other) const
    {
      return hostIntersection_ == other.hostIntersection_;
    }
 
    Entity inside() const {
      return MMeshEntity<0,dim,GridImp>(mMesh_, hostIntersection_.first);
    }
 
    Entity outside() const {
      assert( neighbor() );
      auto neighborHostEntity = (hostIntersection_.first)->neighbor(hostIntersection_.second);
      return MMeshEntity<0,dim,GridImp>(mMesh_, neighborHostEntity);
    }
 
    bool boundary () const {
      auto neighborHostEntity = (hostIntersection_.first)->neighbor(hostIntersection_.second);
      return mMesh_->getHostGrid().is_infinite(neighborHostEntity);
    }
 
    NormalVector centerUnitOuterNormal () const {
      return unitOuterNormal( FieldVector<ctype, dim-1> ( 0 ) );
    }
 
    bool neighbor() const {
      return !boundary();
    }
 
    std::size_t boundarySegmentIndex() const
    {
      HostGridEntity cell  = hostIntersection_.first;
      const auto& facetIdx = hostIntersection_.second;
 
      std::vector< std::size_t > vertices;
      for( std::size_t i = 0; i < dim; ++i )
        vertices.push_back( cell->vertex( (facetIdx+i+1)%(dim+1) )->info().id );
      std::sort(vertices.begin(), vertices.end());
 
      auto it = mMesh_->boundarySegments().find( vertices );
      if( it == mMesh_->boundarySegments().end() )
        return 0; // default to 0
 
      return it->second;
    }
 
    std::size_t boundaryId() const
    {
      auto it = mMesh_->boundaryIds().find( boundarySegmentIndex() );
      if( it == mMesh_->boundaryIds().end() )
        return 0; // default
 
      return it->second;
    }
 
    bool conforming () const {
      // we are always conforming
      return true;
    }
 
    GeometryType type () const {
      return GeometryTypes::simplex(dim-1);
    }
 
    LocalGeometry geometryInInside () const
    {
      LocalGeometryImpl impl( indexInInside() );
      return LocalGeometry( impl );
    }
 
    LocalGeometry geometryInOutside () const
    {
      assert( neighbor() );
      LocalGeometryImpl impl( indexInOutside() );
      return LocalGeometry( impl );
    }
 
    Geometry geometry () const
    {
      return Geometry( hostIntersection_ );
    }
 
    int indexInInside () const {
      return MMeshImpl::cgalFacetToDuneFacet<dim, HostLeafIntersection>( hostIntersection_ );
    }
 
    int indexInOutside () const {
      const auto& neighbor = hostIntersection_.first->neighbor( hostIntersection_.second );
      const auto& second = mMesh_->getHostGrid().mirror_index( hostIntersection_.first, hostIntersection_.second );
      HostLeafIntersection facetFromOutside ( { neighbor, second } );
      return MMeshImpl::cgalFacetToDuneFacet<dim, HostLeafIntersection>( facetFromOutside );
    }
 
    FieldVector<ctype, GridImp::dimensionworld> outerNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const {
      return integrationOuterNormal( local );
    }
 
    template< int d = dim >
    typename std::enable_if_t< d == 2, NormalVector >
    integrationOuterNormal (const FieldVector<ctype, dim-1>& local) const
    {
      HostGridEntity face = hostIntersection_.first;
      const auto& edgeIdx = hostIntersection_.second;
 
      const auto& p1 = face->vertex( face->cw ( edgeIdx ) )->point();
      const auto& p2 = face->vertex( face->ccw( edgeIdx ) )->point();
 
      return NormalVector ( { p1.y() - p2.y(), - p1.x() + p2.x() } );
    }
 
    template< int d = dim >
    typename std::enable_if_t< d == 3, NormalVector >
    integrationOuterNormal (const FieldVector<ctype, dim-1>& local) const
    {
      HostGridEntity cell  = hostIntersection_.first;
      const auto& facetIdx = hostIntersection_.second;
 
      const auto& p1 = cell->vertex( (facetIdx+1)%4 )->point();
      const auto& p2 = cell->vertex( (facetIdx+2)%4 )->point();
      const auto& p3 = cell->vertex( (facetIdx+3)%4 )->point();
 
      return makeFieldVector(
        ( facetIdx % 2 == 0) ?
          CGAL::normal(p1, p2, p3) :
          CGAL::normal(p1, p3, p2)
      );
    }
 
    FieldVector<ctype, GridImp::dimensionworld> unitOuterNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const {
      NormalVector n = integrationOuterNormal( local );
      n /= n.two_norm();
      return n;
    }
 
    const HostLeafIntersection& getHostIntersection() const
    {
      return hostIntersection_;
    }
 
  private:
    const GridImp* mMesh_;
    HostLeafIntersection hostIntersection_;
  };
 
}  // namespace Dune
 
#endif