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/boundarysegment.hh>
#include <dune/grid/common/gridfactory.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