implicitgridfactory.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_MMESH_GRID_IMPLICITGRIDFACTORY_HH
#define DUNE_MMESH_GRID_IMPLICITGRIDFACTORY_HH
 
// System includes
#include <algorithm>
#include <array>
#include <limits>
#include <map>
#include <memory>
#include <unordered_map>
#include <unordered_set>
 
// Dune includes
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/gridfactory.hh>
#include <dune/grid/common/boundarysegment.hh>
 
// MMesh includes
#include <dune/mmesh/grid/mmesh.hh>
#include <dune/mmesh/grid/common.hh>
#include <dune/mmesh/grid/pointfieldvector.hh>
 
namespace Dune
{
 
  template< class Grid >
  class MMeshImplicitGridFactory
    : public GridFactoryInterface< Grid >
  {
    typedef MMeshImplicitGridFactory< Grid > This;
 
  public:
    typedef typename Grid::ctype ctype;
    typedef typename Grid::HostGridType HostGrid;
 
    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 typename Grid::IdType IdType;
 
    typedef std::unordered_map< IdType, std::size_t > BoundarySegments;
    typedef std::unordered_map< std::size_t, std::size_t > BoundaryIds;
 
    typedef std::unordered_map< IdType, std::size_t > InterfaceSegments;
 
    template< int codim >
    struct Codim
    {
      typedef typename Grid::template Codim< codim >::Entity Entity;
    };
 
  private:
    typedef typename HostGrid::Point Point;
    typedef typename HostGrid::Vertex_handle Vertex_handle;
    typedef typename Grid::template HostGridEntity<0> Element_handle;
    typedef typename Grid::template HostGridEntity<1> Face_handle;
 
  public:
    static const bool supportsBoundaryIds = true;
    static const bool supportPeriodicity = false;
 
    MMeshImplicitGridFactory ()
    {}
 
    void insertElement ( const GeometryType &type,
                         const std::vector< unsigned int > &v )
    {
      Element_handle fh;
      if( isElement( v, fh ) )
        fh->info().insertionIndex = countElements;
 
      // Increase element count in each case
      countElements++;
    };
 
    void insertElement ( const GeometryType &type,
                         const std::vector< unsigned int > &v,
                         const size_t domainMarker )
    {
      insertElement( type, v );
    }
 
    template< int d = dimension >
    std::enable_if_t< d == 2, bool >
    isElement( const std::vector< unsigned int >& v, Element_handle& fh ) const
    {
        assert( v.size() == dimension+1 );
        return tr_.is_face( vhs_[v[0]], vhs_[v[1]], vhs_[v[2]], fh );
    }
 
    template< int d = dimension >
    std::enable_if_t< d == 3, bool >
    isElement( const std::vector< unsigned int >& v, Element_handle& fh ) const
    {
        assert( v.size() == dimension+1 );
        return tr_.is_cell( vhs_[v[0]], vhs_[v[1]], vhs_[v[2]], vhs_[v[3]], fh );
    }
 
    void insertBoundarySegment(const std::vector<unsigned int>& vertices)
    {
      assert( vertices.size() == dimension );
 
      std::vector< std::size_t > sorted_vertices;
      for ( const auto& v : vertices )
        sorted_vertices.push_back( vhs_[v]->info().id );
      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<Dune::BoundarySegment<dimension,dimension>>& boundarySegment)
    {
      DUNE_THROW( NotImplemented, "insertBoundarySegments with Dune::BoundarySegment" );
    }
 
    void insertInterfaceBoundarySegment ( const std::vector< unsigned int >& vertices )
    {
      assert( vertices.size() == dimension-1 );
 
      std::vector< std::size_t > sorted_vertices;
      for ( const auto& v : vertices )
        sorted_vertices.push_back( vhs_[v]->info().id );
      std::sort(sorted_vertices.begin(), sorted_vertices.end());
 
      if( boundarySegments_.find( sorted_vertices ) != boundarySegments_.end() )
          DUNE_THROW( GridError, "A boundary segment was inserted twice." );
 
      interfaceBoundarySegments_.insert( std::make_pair( sorted_vertices, countInterfaceBoundarySegments++ ) );
    }
 
    void insertVertex ( const WorldVector &pos )
    {
      // Insert vertex
      Vertex_handle vh = tr_.insert( makePoint( pos ) );
 
      // Store insertion index
      vh->info().id = countVertices;
      vh->info().idWasSet = true;
 
      // Increase vertex counter
      countVertices++;
 
      // Store vertex handle for later use
      vhs_.push_back( vh );
    }
 
    void insertInterface ( const std::vector< unsigned int > &vertices, const std::size_t marker = 1 )
    {
      assert( vertices.size() == dimension );
 
      std::vector< std::size_t > sorted_vertices;
      for( const auto& v : vertices )
        sorted_vertices.push_back( vhs_[v]->info().id );
      std::sort(sorted_vertices.begin(), sorted_vertices.end());
      interfaceSegments_.insert( std::make_pair( sorted_vertices, marker ) );
    }
 
    unsigned int insertionIndex ( const typename Codim<0>::Entity &entity ) const
    {
      return entity.impl().hostEntity()->info().insertionIndex;
    }
 
    unsigned int insertionIndex ( const typename Codim< dimension >::Entity &entity ) const
    {
      return entity.impl().hostEntity()->info().id;
    }
 
    unsigned int insertionIndex ( const typename Grid::LeafIntersection &intersection ) const
    {
      return intersection.impl().boundaryId();
    }
 
    const BoundarySegments& boundarySegments() const
    {
      return boundarySegments_;
    }
 
    const BoundaryIds& boundaryIds() const
    {
      return boundaryIds_;
    }
 
    void addBoundaryId( std::size_t boundarySegmentIndex, std::size_t boundaryId )
    {
      boundaryIds_.insert( std::make_pair( boundarySegmentIndex, boundaryId ) );
    }
 
    std::unique_ptr<Grid> createGrid ()
    {
      return std::make_unique<Grid> (
        std::move(tr_),
        std::move(boundarySegments_),
        std::move(interfaceBoundarySegments_),
        std::move(boundaryIds_),
        std::move(interfaceSegments_)
      );
    }
 
  private:
    template< int dim = dimension >
    std::enable_if_t< dim == 2, Point >
    makePoint( const WorldVector& v ) const
    {
      return Point ( v[ 0 ], v[ 1 ] );
    }
 
    template< int dim = dimension >
    std::enable_if_t< dim == 3, Point >
    makePoint( const WorldVector& v ) const
    {
      return Point ( v[ 0 ], v[ 1 ], v[ 2 ] );
    }
 
    HostGrid tr_;
    std::vector< Vertex_handle > vhs_;
    BoundarySegments boundarySegments_, interfaceBoundarySegments_;
    BoundaryIds boundaryIds_;
    InterfaceSegments interfaceSegments_;
    std::size_t countVertices = 0, countElements = 0;
    std::size_t countBoundarySegments = 0, countInterfaceBoundarySegments = 0;
  };
 
} // end namespace Dune
 
#endif