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_INTERFACE_INTERSECTIONS_HH
#define DUNE_MMESH_INTERFACE_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 MMeshInterfaceGridLeafIntersection
  {
    friend class MMeshInterfaceGridLeafIntersectionIterator<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 MMeshInterfaceEntity<0> MMeshInterfaceEntity;
    typedef typename GridImp::template MMeshInterfaceEntity<1> HostLeafIntersection;
 
    typedef typename GridImp::MMeshType::template Codim<0>::Entity MMeshEntity;
 
    using LocalIndexMap = std::unordered_map< std::size_t, std::size_t >;
    typedef typename GridImp::MMeshType::template Codim<dimworld>::Entity MMeshVertex;
 
  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<0>::Entity Entity;
    typedef FieldVector<ctype, dimensionworld> NormalVector;
 
    MMeshInterfaceGridLeafIntersection()
    {}
 
    MMeshInterfaceGridLeafIntersection(const GridImp* grid,
                          const MMeshInterfaceEntity& hostEntity,
                          const std::size_t index,
                          const std::size_t nbIdx)
      : grid_(grid),
        interfaceEntity_(hostEntity),
        index_(index),
        nbIdx_(nbIdx),
        cgalIndex_( MMeshInterfaceImpl::computeCGALIndices<MMeshInterfaceEntity, dim>( interfaceEntity_ ) )
    {
      const auto& indexSet = grid_->leafIndexSet();
 
      std::array< std::size_t, dimension > ids;
 
      if constexpr (dim == 1)
      {
        ids[0] = indexSet.vertexIndexMap().at( interfaceEntity_.first->vertex( cgalIndex_[index_] )->info().id );
      }
      else // dim == 2
      {
        ids[0] = indexSet.vertexIndexMap().at( interfaceEntity_.first->vertex( cgalIndex_[index_==2 ? 1 : 0] )->info().id );
        ids[1] = indexSet.vertexIndexMap().at( interfaceEntity_.first->vertex( cgalIndex_[index_==0 ? 1 : 2] )->info().id );
      }
      std::sort(ids.begin(), ids.end());
      try {
        localIndexMap_ = indexSet.indexMap().at( ids );
      } catch (std::exception &e) {
        DUNE_THROW(InvalidStateException, e.what());
      }
      assert( 0 <= nbIdx_ && ( nbIdx_ < numOutside() || boundary() ) );
    }
 
    bool equals(const MMeshInterfaceGridLeafIntersection& other) const
    {
      return interfaceEntity_ == other.interfaceEntity_
          && index_ == other.index_
          && nbIdx_ == other.nbIdx_;
    }
 
    Entity inside() const {
      return MMeshInterfaceGridEntity<0,dimension,GridImp>(grid_, interfaceEntity_);
    }
 
    std::size_t numOutside() const {
      return localIndexMap_.size() - 1;
    }
 
    Entity outside() const
    {
      assert( neighbor() );
 
      const auto& vertexIndexMap = grid_->leafIndexSet().vertexIndexMap();
 
      std::size_t myLastVertexIndex = vertexIndexMap.at(
        interfaceEntity_.first->vertex(cgalIndex_[dim-index_])->info().id
      );
 
      // get the i-th index map entry
      auto it = localIndexMap_.begin();
      // exclude myself
      while ( it->first == myLastVertexIndex )
        ++it;
 
      for( std::size_t count = 0; count < nbIdx_; count++ )
      {
        ++it;
 
        // exclude myself
        if ( it->first == myLastVertexIndex )
          count--;
      }
 
      // obtain neighbor by searching in incident elements of intersection vertex
      std::size_t v0idx = index_;
      if constexpr (dim == 2)
        if (index_ == 1)
          v0idx = 0;
 
      const auto& vertex0 = interfaceEntity_.first->vertex( cgalIndex_[v0idx] );
 
      std::array< std::size_t, dimensionworld > vIdx;
      vIdx[0] = vertexIndexMap.at( vertex0->info().id );
      vIdx[1] = it->first;
 
      if constexpr ( dimensionworld == 3 )
        try {
          std::size_t v1idx = (index_ == 1) ? 2 : 1;
          vIdx[2] = vertexIndexMap.at( interfaceEntity_.first->vertex( cgalIndex_[v1idx] )->info().id );
        } catch (std::exception &e) {
          DUNE_THROW(InvalidStateException, e.what());
        }
      std::sort(vIdx.begin(), vIdx.end());
 
      // search in incident facets for the right entity
      for( const auto& facet : incidentFacets( MMeshVertex {{ &grid_->getMMesh(), vertex0 }} ) )
      {
        std::array< std::size_t, dimensionworld > tmpIdx;
        bool notInterface = false;
        for( std::size_t i = 0; i < facet.subEntities(dimensionworld); ++i )
        {
          std::size_t idx = facet.impl().template subEntity<dimensionworld>(i).impl().hostEntity()->info().id;
          auto it = vertexIndexMap.find( idx );
          if( it != vertexIndexMap.end() )
            tmpIdx[i] = it->second;
          else
            notInterface = true;
        }
        if( notInterface )
          continue;
        std::sort(tmpIdx.begin(), tmpIdx.end());
 
        if ( vIdx == tmpIdx )
          return MMeshInterfaceGridEntity<0,dimension,GridImp>(grid_, facet.impl().hostEntity());
      }
      DUNE_THROW( InvalidStateException, "Neighbor could not be determined!" );
    }
 
    bool boundary () const {
      return localIndexMap_.size() == 1;
    }
 
    NormalVector centerUnitOuterNormal () const {
      return unitOuterNormal( FieldVector<ctype, dimension-1> ( 0 ) );
    }
 
    bool neighbor() const {
      return !boundary();
    }
 
    size_t boundarySegmentIndex() const
    {
      auto iid = grid_->globalIdSet().id( grid_->entity( getHostIntersection() ) );
      auto it = grid_->boundarySegments().find( iid );
      if( it == grid_->boundarySegments().end() )
        return 0;
 
      return it->second;
    }
 
    std::size_t boundaryId() const
    {
      return 0; // TODO
    }
 
    bool conforming () const {
      // we are always conforming
      return true;
    }
 
    GeometryType type () const {
      return GeometryTypes::simplex(dimension-1);
    }
 
    LocalGeometry geometryInInside () const
    {
      return LocalGeometry( indexInInside() );
    }
 
    LocalGeometry geometryInOutside () const
    {
      assert( neighbor() );
      return LocalGeometry( indexInOutside() );
    }
 
    template< int d = dimension >
    typename std::enable_if_t< d == 1, Geometry >
    geometry () const
    {
      return Geometry( interfaceEntity_.first->vertex( cgalIndex_[index_] ) );
    }
 
    template< int d = dimension >
    typename std::enable_if_t< d == 2, Geometry >
    geometry () const
    {
      int vIdx0 = cgalIndex_[index_==2 ? 1 : 0];
      int vIdx1 = cgalIndex_[index_==0 ? 1 : 2];
 
      return Geometry( {{ interfaceEntity_.first, vIdx0, vIdx1 }} );
    }
 
    int indexInInside () const
    {
      return index_;
    }
 
    int indexInOutside () const
    {
      const auto& subEn = inside().template subEntity<1>( index_ );
 
      const auto neighbor = outside();
      for ( int i = 0; i < dimensionworld; ++i )
        if ( neighbor.template subEntity<1>( i ) == subEn )
          return i;
 
      DUNE_THROW( InvalidStateException, "indexInOutside() could not be determined!" );
    }
 
    FieldVector<ctype, GridImp::dimensionworld> outerNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const {
      return integrationOuterNormal( local );
    }
 
    template< int d = dimension >
    typename std::enable_if_t< d == 1, NormalVector >
    integrationOuterNormal (const FieldVector<ctype, dimension-1>& local) const
    {
      auto face = interfaceEntity_.first;
 
      const auto& p1 = face->vertex( cgalIndex_[ index_ ] )->point();
      const auto& p2 = face->vertex( cgalIndex_[1-index_] )->point();
 
      NormalVector n ( { p1.x() - p2.x(), p1.y() - p2.y() } );
      n /= n.two_norm();
      return n;
    }
 
    template< int d = dimension >
    typename std::enable_if_t< d == 2, NormalVector >
    integrationOuterNormal (const FieldVector<ctype, dimension-1>& local) const
    {
      const auto& cell = interfaceEntity_.first;
 
      int vIdx0 = cgalIndex_[index_==2 ? 1 : 0];
      int vIdx1 = cgalIndex_[index_==0 ? 1 : 2];
      int vIdxLast = cgalIndex_[2-index_];
 
      auto a = makeFieldVector( cell->vertex( vIdx0 )->point() );
      auto b = makeFieldVector( cell->vertex( vIdx1 )->point() );
      auto v = makeFieldVector( cell->vertex( vIdxLast )->point() );
 
      // return vector that is orthogonal to edge a-b and triangle a-b-v
      NormalVector ab = b - a;
      NormalVector vb = b - v;
 
      NormalVector elementNormal;
      elementNormal[0] = ab[1]*vb[2] - ab[2]*vb[1];
      elementNormal[1] = ab[2]*vb[0] - ab[0]*vb[2];
      elementNormal[2] = ab[0]*vb[1] - ab[1]*vb[0];
 
      NormalVector outerNormal;
      outerNormal[0] = ab[1]*elementNormal[2] - ab[2]*elementNormal[1];
      outerNormal[1] = ab[2]*elementNormal[0] - ab[0]*elementNormal[2];
      outerNormal[2] = ab[0]*elementNormal[1] - ab[1]*elementNormal[0];
 
      if( outerNormal * vb < 0.0 )
        outerNormal *= -1.0;
 
      return outerNormal;
    }
 
    FieldVector<ctype, GridImp::dimensionworld> unitOuterNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const {
      NormalVector n = integrationOuterNormal( local );
      n /= n.two_norm();
      return n;
    }
 
    const auto getHostVertex() const
    {
      return interfaceEntity_.first->vertex((interfaceEntity_.second+index_+1)%(dimensionworld+1));
    }
 
    const MMeshInterfaceEntity& getHostIntersection() const
    {
      return interfaceEntity_;
    }
 
  private:
    const GridImp* grid_;
    MMeshInterfaceEntity interfaceEntity_;
    std::size_t index_;
    std::size_t nbIdx_;
    LocalIndexMap localIndexMap_;
    std::array< std::size_t, dim+1 > cgalIndex_;
  };
 
}  // namespace Dune
 
#endif