explicitgridfactory.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_EXPLICITGRIDFACTORY_HH
#define DUNE_MMESH_GRID_EXPLICITGRIDFACTORY_HH
 
// System includes
#include <algorithm>
#include <array>
#include <limits>
#include <map>
#include <memory>
 
// 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>
 
namespace Dune {
 
template <class Grid>
class MMeshExplicitGridFactory : public GridFactoryInterface<Grid> {
  typedef MMeshExplicitGridFactory<Grid> This;
  typedef GridFactoryInterface<Grid> Base;
 
 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;
 
  using Tuple = std::tuple<std::vector<unsigned int>, std::size_t, std::size_t>;
 
  using Base::insertElement;
 
 public:
  static const bool supportsBoundaryIds = true;
  static const bool supportPeriodicity = false;
 
  MMeshExplicitGridFactory() { tr_.tds().clear(); }
 
  void insertElement(const GeometryType& type,
                     const std::vector<unsigned int>& v) {
    insertElement(type, v, 0);
  }
 
  void insertElement(const GeometryType& type,
                     const std::vector<unsigned int>& v,
                     const size_t domainMarker) {
    assert(type == GeometryTypes::simplex(dimension));
    assert(v.size() == dimension + 1);
    auto w = v;
 
    // Create element
    storeElement(w, countElements, domainMarker);
 
    // Increase element count
    countElements++;
 
    // Store vertices to check occurence of element later
    elementVerticesList.push_back(w);
  };
 
 private:
  void storeElement(std::vector<unsigned int>& v, const size_t insertionIndex,
                    const size_t domainMarker) {
    elements_.push_back(std::make_tuple(v, insertionIndex, domainMarker));
  }
 
  template <int d = dimension>
  std::enable_if_t<d == 2, void> createElement(std::vector<unsigned int>& v,
                                               const size_t insertionIndex,
                                               const size_t domainMarker) {
    auto& p0 = vhs_[v[0]]->point();
    auto& p1 = vhs_[v[1]]->point();
    auto& p2 = vhs_[v[2]]->point();
 
    // Check orientation of vertices
    auto orientation = (p1.y() - p0.y()) * (p2.x() - p1.x()) -
                       (p1.x() - p0.x()) * (p2.y() - p1.y());
    // if clockwise, swap two vertices
    if (orientation > 0.0) std::swap(v[0], v[1]);
 
    auto&& v0 = vhs_[v[0]];
    auto&& v1 = vhs_[v[1]];
    auto&& v2 = vhs_[v[2]];
 
    // Create face with vertices v0, v1, v2
    Element_handle face = tr_.tds().create_face(v0, v1, v2);
 
    // Set insertion index
    face->info().insertionIndex = insertionIndex;
    face->info().domainMarker = domainMarker;
 
    // Set this face in vertices v0, v1, v2
    v0->set_face(face);
    v1->set_face(face);
    v2->set_face(face);
 
    // Add facets to neighbor map to obtain connectivity
    addFacetToMap({v[0], v[1]}, face, 2);
    addFacetToMap({v[0], v[2]}, face, 1);
    addFacetToMap({v[1], v[2]}, face, 0);
  }
 
  template <int d = dimension>
  std::enable_if_t<d == 3, void> createElement(
      const std::vector<unsigned int>& v, const size_t insertionIndex,
      const size_t domainMarker) {
    auto&& v0 = vhs_[v[0]];
    auto&& v1 = vhs_[v[1]];
    auto&& v2 = vhs_[v[2]];
    auto&& v3 = vhs_[v[3]];
 
    // Create cell with vertices v0, v1, v2, v3
    Element_handle cell = tr_.tds().create_cell(v0, v1, v2, v3);
 
    // Set insertion index
    cell->info().insertionIndex = insertionIndex;
    cell->info().domainMarker = domainMarker;
 
    // Set this cell in vertices v0, v1, v2, v3
    v0->set_cell(cell);
    v1->set_cell(cell);
    v2->set_cell(cell);
    v3->set_cell(cell);
 
    // Add facets to neighbor map to obtain connectivity
    addFacetToMap({v[0], v[1], v[2]}, cell, 3);
    addFacetToMap({v[0], v[1], v[3]}, cell, 2);
    addFacetToMap({v[0], v[2], v[3]}, cell, 1);
    addFacetToMap({v[1], v[2], v[3]}, cell, 0);
  }
 
  void addFacetToMap(const std::vector<unsigned int>& v,
                     const Element_handle& element, const int fi) {
    assert(v.size() == dimension);
 
    // Make a set of the vertex indices
    std::set<unsigned int> facetIndices(v.begin(), v.end());
 
    // Try to insert this set into the neighborMap
    auto entry = neighborMap.insert(
        {facetIndices, std::pair<Element_handle, int>(element, fi)});
 
    // If facetIndices was already in neighborMap, connect the neighbors and
    // remove the entry
    if (!entry.second) {
      auto facet = entry.first->second;
      Element_handle neighbor = facet.first;
      int ni = facet.second;
 
      element->set_neighbor(fi, neighbor);
      neighbor->set_neighbor(ni, element);
 
      neighborMap.erase(facetIndices);
    }
  }
 
 public:
  template <int d = dimension>
  std::enable_if_t<d == 2, bool> isElement(
      const std::vector<unsigned int>& v) const {
    assert(v.size() == dimension + 1);
    return tr_.is_face(vhs_[v[0]], vhs_[v[1]], vhs_[v[2]]);
  }
 
  template <int d = dimension>
  std::enable_if_t<d == 3, bool> isElement(
      const std::vector<unsigned int>& v) const {
    assert(v.size() == dimension + 1);
    Element_handle cell;
    return tr_.is_cell(vhs_[v[0]], vhs_[v[1]], vhs_[v[2]], vhs_[v[3]], cell);
  }
 
  virtual 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<BoundarySegment>& 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
    auto vh = tr_.tds().create_vertex();
    vh->set_point(makePoint(pos));
 
    // Store insertion index
    vh->info().id = countVertices;
    vh->info().idWasSet = true;
 
    // Increase vertex counter
    countVertices++;
 
    // Store the inserted 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() {
    // Sort elements by x-coordinate of center for partitioning
    std::sort(elements_.begin(), elements_.end(),
              [this](const auto& a, const auto& b) {
                const auto& va = std::get<0>(a);
                const auto& vb = std::get<0>(b);
 
                double xa = 0.0;
                double xb = 0.0;
                for (int i = 0; i < dimensionworld + 1; ++i) {
                  xa += vhs_[va[i]]->point().x();
                  xb += vhs_[vb[i]]->point().x();
                }
                return xa < xb;
              });
 
    for (auto& t : elements_)
      createElement(std::get<0>(t), std::get<1>(t), std::get<2>(t));
 
    // Create the infinite cells (neighbors of boundary cells)
    createInfiniteVertex();
 
    // Create the infinite cells (neighbors of boundary cells)
    createInfiniteCells();
 
    // Remove interfaceSegments_ from boundarySegments_
    for (const auto& interfaceSeg : interfaceSegments_)
      boundarySegments_.erase(interfaceSeg.first);
 
    // Mark interface vertices as isInterface
    for (const auto& interfaceSeg : interfaceSegments_)
      for (const auto& v : interfaceSeg.first.vt())
        vhs_[v]->info().isInterface = true;
 
    // Check if all inserted elements really exist in the triangulation
    checkOccurenceOfAllElements();
 
    // Return pointer to grid
    return std::make_unique<Grid>(std::move(tr_), std::move(boundarySegments_),
                                  std::move(interfaceBoundarySegments_),
                                  std::move(boundaryIds_),
                                  std::move(interfaceSegments_));
  }
 
  const std::vector<Vertex_handle>& vertexHandles() const { return vhs_; }
 
 private:
  void createInfiniteVertex() {
    if (countElements > 0)
      tr_.tds().set_dimension(dimension);
    else
      tr_.tds().set_dimension(0);
 
    Vertex_handle infinite = tr_.tds().create_vertex();
    infinite->info().id = std::size_t(-1);
    tr_.set_infinite_vertex(infinite);
  }
 
  template <int d = dimension>
  std::enable_if_t<d == 2, void> createInfiniteCells() {
    // Iterate over all unique facets
    for (const auto& entry : neighborMap) {
      auto facet = entry.second;
      Element_handle face = facet.first;
      int fi = facet.second;
 
      // Remove real boundary segments from interfaceSegments_
      std::vector<std::size_t> vertices;
      vertices.push_back(face->vertex((fi + 2) % 3)->info().id);
      vertices.push_back(face->vertex((fi + 1) % 3)->info().id);
      std::sort(vertices.begin(), vertices.end());
 
      auto it = boundarySegments_.find(vertices);
      if (it != boundarySegments_.end()) interfaceSegments_.erase(vertices);
 
      // Create infinite face with correct orientation
      Element_handle iface = tr_.tds().create_face(face->vertex((fi + 2) % 3),
                                                   face->vertex((fi + 1) % 3),
                                                   tr_.infinite_vertex());
 
      tr_.infinite_vertex()->set_face(iface);
 
      face->set_neighbor(fi, iface);
      iface->set_neighbor(2, face);
 
      // Map infinite neighbors
      for (int i = 0; i < 2; ++i) {
        std::set<std::size_t> vertexIndex({iface->vertex(i)->info().id});
        auto entry = infiniteNeighborMap.insert(
            {vertexIndex, std::pair<Element_handle, int>(iface, (i + 1) % 2)});
 
        // If vertexIndex was already in map, connect neighbors and remove entry
        if (!entry.second) {
          auto pair = entry.first->second;
          Element_handle neighbor = pair.first;
          int ni = pair.second;
 
          iface->set_neighbor((i + 1) % 2, neighbor);
          neighbor->set_neighbor(ni, iface);
 
          infiniteNeighborMap.erase(vertexIndex);
        }
      }
    }
 
    // Assert that all boundary facets are mapped
    assert(infiniteNeighborMap.size() == 0);
  }
 
  template <int d = dimension>
  std::enable_if_t<d == 3, void> createInfiniteCells() {
    // Iterate over all unique facets
    for (const auto& entry : neighborMap) {
      auto facet = entry.second;
      Element_handle cell = facet.first;
      int fi = facet.second;
 
      // Remove real boundary segments from interfaceSegments_
      std::vector<std::size_t> vertices;
      vertices.push_back(
          cell->vertex((fi % 2 == 1) ? (fi + 2) & 3 : (fi + 1) & 3)->info().id);
      vertices.push_back(
          cell->vertex((fi % 2 == 1) ? (fi + 1) & 3 : (fi + 2) & 3)->info().id);
      vertices.push_back(cell->vertex((fi + 3) & 3)->info().id);
      std::sort(vertices.begin(), vertices.end());
 
      auto it = boundarySegments_.find(vertices);
      if (it != boundarySegments_.end()) interfaceSegments_.erase(vertices);
 
      // Create infinite cell with correct orientation
      Element_handle icell = tr_.tds().create_cell(
          cell->vertex((fi % 2 == 1) ? (fi + 2) & 3 : (fi + 1) & 3),
          cell->vertex((fi % 2 == 1) ? (fi + 1) & 3 : (fi + 2) & 3),
          cell->vertex((fi + 3) & 3), tr_.infinite_vertex());
 
      tr_.infinite_vertex()->set_cell(icell);
 
      cell->set_neighbor(fi, icell);
      icell->set_neighbor(3, cell);
 
      // Map infinite neighbors
      for (int i = 0; i < 3; ++i) {
        std::set<std::size_t> edgeIndices(
            {icell->vertex(i)->info().id,
             icell->vertex((i + 1) % 3)->info().id});
        auto entry = infiniteNeighborMap.insert(
            {edgeIndices, std::pair<Element_handle, int>(icell, (i + 2) % 3)});
 
        // If edgeIndices was already in map, connect neighbors and remove entry
        if (!entry.second) {
          auto pair = entry.first->second;
          Element_handle neighbor = pair.first;
          int ni = pair.second;
 
          icell->set_neighbor((i + 2) % 3, neighbor);
          neighbor->set_neighbor(ni, icell);
 
          infiniteNeighborMap.erase(edgeIndices);
        }
      }
    }
 
    // Assert that all boundary facets are mapped
    assert(infiniteNeighborMap.size() == 0);
  }
 
  void checkOccurenceOfAllElements() const {
    for (const auto& v : elementVerticesList)
      if (!isElement(v))
        DUNE_THROW(InvalidStateException,
                   "Inserted element was not found in CGAL triangulation.");
  }
 
  HostGrid tr_;
  std::vector<Vertex_handle> vhs_;
  std::vector<Tuple> elements_;
  BoundarySegments boundarySegments_, interfaceBoundarySegments_;
  BoundaryIds boundaryIds_;
  InterfaceSegments interfaceSegments_;
  std::vector<std::vector<unsigned int> > elementVerticesList;
  std::map<std::set<unsigned int>, std::pair<Element_handle, int> > neighborMap;
  std::map<std::set<std::size_t>, std::pair<Element_handle, int> >
      infiniteNeighborMap;
  std::size_t countElements = 0, countVertices = 0;
  std::size_t countBoundarySegments = 0, countInterfaceBoundarySegments = 0;
};
 
}  // end namespace Dune
 
#endif