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