gridfactory.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_MMESH_INTERFACE_GRIDFACTORY_HH
#define DUNE_MMESH_INTERFACE_GRIDFACTORY_HH
 
// Dune includes
#include <dune/grid/common/gridfactory.hh>
#include <dune/grid/common/boundarysegment.hh>
 
namespace Dune
{
 
  template< class MMeshImp >
  class GridFactory< MMeshInterfaceGrid<MMeshImp> >
    : public GridFactoryInterface< MMeshInterfaceGrid<MMeshImp> >
  {
    typedef GridFactory< MMeshInterfaceGrid<MMeshImp> > This;
 
  public:
    typedef MMeshInterfaceGrid<MMeshImp> Grid;
 
    typedef typename Grid::ctype ctype;
    typedef typename Grid::HostGridType HostGrid;
    typedef typename HostGrid::Vertex_handle VertexHandle;
 
    typedef MMeshImp MMesh;
 
    static const int dimension = Grid::dimension;
    static const int dimensionworld = Grid::dimensionworld;
 
    typedef FieldVector< ctype, dimensionworld > WorldVector;
    typedef FieldMatrix< ctype, dimensionworld, dimensionworld > WorldMatrix;
 
    typedef Dune::BoundarySegment< dimension, dimensionworld > BoundarySegment;
    typedef std::unordered_map< std::vector< std::size_t >, std::size_t, HashUIntVector > BoundarySegments;
    typedef std::map< std::vector< std::size_t >, unsigned int > InsertionIndexMap;
 
    typedef std::map< std::size_t, std::size_t > VertexIdMap;
 
    template< int codim >
    struct Codim
    {
      typedef typename Grid::template Codim< codim >::Entity Entity;
    };
 
  public:
    static const bool supportsBoundaryIds = true;
    static const bool supportPeriodicity = false;
 
    GridFactory ( const std::shared_ptr<MMesh> mMesh )
     : mMesh_( mMesh )
    {}
 
    GridFactory ()
    {
      DUNE_THROW(NotImplemented, "GridFactory() for MMeshInterfaceGrid");
    }
 
    void insertElement ( const GeometryType &type,
                         const std::vector< unsigned int > &vertices )
    {
      // mark vertices as interface segment in host mmesh
      std::vector< std::size_t > ids;
      for ( const auto& v : vertices )
        ids.push_back( vertexIdMap_.at( v ) );
 
      std::sort(ids.begin(), ids.end());
 
      (mMesh_->interfaceSegments()).insert( std::make_pair(ids, 1) );
 
      insertionIndexMap_.insert( { ids, countElements++ } );
    };
 
    virtual void insertBoundarySegment ( const std::vector< unsigned int >& vertices )
    {
      std::vector< std::size_t > sorted_vertices;
      for ( const auto& v : vertices )
        sorted_vertices.push_back( vertexIdMap_.at( v ) );
      std::sort(sorted_vertices.begin(), sorted_vertices.end());
 
      if( boundarySegments_.find( sorted_vertices ) != boundarySegments_.end() )
        DUNE_THROW( GridError, "A boundary segment was inserted twice." );
 
      boundarySegments_.insert( std::make_pair( sorted_vertices, countBoundarySegments++ ) );
    }
 
    void insertBoundarySegment ( const std::vector< unsigned int >& vertices,
                                 const std::shared_ptr< BoundarySegment >& boundarySegment )
    {
      DUNE_THROW( NotImplemented, "insertBoundarySegments with Dune::BoundarySegment" );
    }
 
    void insertVertex ( const WorldVector &pos )
    {
      // get the vertex handle from the host mmesh
      VertexHandle vh = mMesh_->getHostGrid().insert( makePoint( pos ) );
 
      vertexIdMap_.insert( { countVertices, vh->info().id } );
      vh->info().isInterface = true;
 
      countVertices++;
    }
 
    void addVertexHandle ( const VertexHandle &vh )
    {
      vertexIdMap_.insert( { countVertices, vh->info().id } );
      vh->info().isInterface = true;
 
      countVertices++;
    }
 
    unsigned int insertionIndex ( const typename Codim<0>::Entity &entity ) const
    {
      std::vector< std::size_t > ids;
      for( std::size_t i = 0; i < entity.subEntities(dimension); ++i )
        ids.push_back( entity.template subEntity<dimension>(i).impl().hostEntity()->info().id );
      std::sort(ids.begin(), ids.end());
      auto it = insertionIndexMap_.find( ids );
      if( it != insertionIndexMap_.end() )
        return it->second;
      else
        return 0; // should not happen
    }
 
    unsigned int insertionIndex ( const typename Codim< dimension >::Entity &entity ) const
    {
      std::size_t index = mMesh_->interfaceGrid().leafIndexSet().index( entity );
      assert( index < std::numeric_limits<unsigned int>::max() );
      return index;
    }
 
    typename Grid::GridPtrType createGrid ()
    {
      DUNE_THROW( InvalidStateException, "The interface grid cannot be created, get the pointer by getGrid()!" );
      return nullptr;
    }
 
    auto getGrid ()
    {
      mMesh_->interfaceGridPtr()->setIds();
      mMesh_->interfaceGridPtr()->setIndices();
      mMesh_->interfaceGridPtr()->setBoundarySegments( boundarySegments_ );
 
      // Return pointer to grid
      return mMesh_->interfaceGridPtr();
    }
 
  private:
    std::shared_ptr<MMesh> mMesh_;
    BoundarySegments boundarySegments_;
    std::size_t countBoundarySegments = 0;
    VertexIdMap vertexIdMap_;
    InsertionIndexMap insertionIndexMap_;
    unsigned int countElements = 0;
    std::size_t countVertices = 0;
  };
 
} // end namespace Dune
 
#endif