indexsets.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_MMESH_GRID_INDEXSETS_HH
#define DUNE_MMESH_GRID_INDEXSETS_HH
 
// Dune includes
#include <dune/grid/common/indexidset.hh>
#include <dune/mmesh/grid/multiid.hh>
 
namespace Dune {
 
template <class GridImp>
class MMeshLeafIndexSet : public IndexSet<GridImp, MMeshLeafIndexSet<GridImp>> {
  typedef typename std::remove_const<GridImp>::type::HostGridType HostGrid;
 
  template <int codim>
  using HostGridEntity = typename GridImp::template HostGridEntity<codim>;
 
 public:
  typedef std::size_t IndexType;
  typedef const std::vector<GeometryType> Types;
 
  typedef std::unordered_map<MMeshImpl::MultiId, std::size_t> CodimIndexMap;
 
  /*
   * We use the remove_const to extract the Type from the mutable class,
   * because the const class is not instantiated yet.
   */
  enum { dim = std::remove_const<GridImp>::type::dimension };
 
  MMeshLeafIndexSet(const GridImp* grid) : grid_(grid) {}
 
  MMeshLeafIndexSet(const MMeshLeafIndexSet& leafIndexSet)
      : grid_(leafIndexSet.grid_), sizeOfCodim_(leafIndexSet.sizeOfCodim_) {}
 
  template <int codim>
  std::enable_if_t<codim == 0 || codim == dim, IndexType> index(
      const Entity<codim, dim, GridImp, MMeshEntity>& e) const {
    auto hostEntity = e.impl().hostEntity();
    IndexType index = hostEntity->info().index;
    assert(index <= size(codim));
    return index;
  }
 
  template <int codim>
  std::enable_if_t<codim == 1, IndexType> index(
      const Entity<codim, dim, GridImp, MMeshEntity>& e) const {
    try {
      return codimIndexMap_[0].at(grid_->globalIdSet().id(e));
    } catch (...) {
      DUNE_THROW(InvalidStateException, "Id of codim 1 entity not found.");
    }
  }
 
  template <int codim>
  std::enable_if_t<codim == 2 && dim == 3, IndexType> index(
      const Entity<codim, dim, GridImp, MMeshEntity>& e) const {
    try {
      return codimIndexMap_[1].at(grid_->globalIdSet().id(e));
    } catch (...) {
      DUNE_THROW(InvalidStateException, "Id of codim 2 entity not found.");
    }
  }
 
  template <class Entity>
  IndexType subIndex(const Entity& e, int i, int codim) {
    return subIndex<Entity::codimension>(e, i, codim);
  }
 
  template <int cc>
  std::enable_if_t<cc == dim, IndexType> subIndex(
      const typename std::remove_const<GridImp>::type::Traits::template Codim<
          cc>::Entity& e,
      int i, int codim) const {
    const HostGridEntity<dim> hostEntity = e.impl().hostEntity();
    return hostEntity->info().index;
  };
 
  template <int cc>
  std::enable_if_t<cc == 0, IndexType> subIndex(
      const typename std::remove_const<GridImp>::type::Traits::template Codim<
          cc>::Entity& e,
      int i, int codim) const {
    assert(codim >= 0 && codim <= dim);
 
    if (codim == 0) return index(e);
    if (codim == dim)
      return e.impl()
          .template subEntity<dim>(i)
          .impl()
          .hostEntity()
          ->info()
          .index;
    else if (codim == 1)
      return codimIndexMap_[0].at(
          grid_->globalIdSet().id(e.impl().template subEntity<1>(i)));
    else if (codim == 2)
      return codimIndexMap_[1].at(
          grid_->globalIdSet().id(e.impl().template subEntity<2>(i)));
    else
      DUNE_THROW(InvalidStateException,
                 "subIndex() was called for codim " << codim);
 
    return 0;
  }
 
  template <int cc>
  std::enable_if_t<cc != 0 && cc != dim, IndexType> subIndex(
      const typename std::remove_const<GridImp>::type::Traits::template Codim<
          cc>::Entity& e,
      int i, int codim) const {
    DUNE_THROW(NotImplemented, "SubIndices for codim != 0 || codim != dim.");
  };
 
  std::size_t size(GeometryType type) const {
    if (type == GeometryTypes::vertex)
      return size(dim);
    else if (type == GeometryTypes::line)
      return size(dim - 1);
    else if (type == GeometryTypes::triangle)
      return size(dim - 2);
    else if (type == GeometryTypes::tetrahedron)
      return size(dim - 3);
    else
      return 0;
  }
 
  std::size_t size(int codim) const {
    assert((0 <= codim) && (codim <= dim));
    return sizeOfCodim_[codim];
  }
 
  const Types geomTypes(int codim) const { return types(codim); }
 
  Types types(int codim) const {
    switch (dim - codim) {
      case 0:
        return {{GeometryTypes::vertex}};
      case 1:
        return {{GeometryTypes::line}};
      case 2:
        return {{GeometryTypes::triangle}};
      case 3:
        return {{GeometryTypes::tetrahedron}};
      default:
        DUNE_THROW(InvalidStateException,
                   "Codim is not within 0 <= codim <= dim.");
    }
  }
 
  template <class EntityType, int d = dim>
  std::enable_if_t<d == 2 && EntityType::codimension == 0, bool> contains(
      const EntityType& e) const {
    const auto hostEntity = e.impl().hostEntity();
    return grid_->getHostGrid().is_face(
        hostEntity->vertex(0), hostEntity->vertex(1), hostEntity->vertex(2));
  }
 
  template <class EntityType, int d = dim>
  std::enable_if_t<d == 2 && EntityType::codimension == 1, bool> contains(
      const EntityType& e) const {
    const auto hostEntity = e.impl().hostEntity();
    return grid_->getHostGrid().is_edge(
        hostEntity.first->vertex((hostEntity.second + 1) % 3),
        hostEntity.first->vertex((hostEntity.second + 2) % 3));
  }
 
  template <class EntityType, int d = dim>
  std::enable_if_t<d == 2 && EntityType::codimension == 2, bool> contains(
      const EntityType& e) const {
    const auto hostEntity = e.impl().hostEntity();
    return grid_->getHostGrid().tds().is_vertex(hostEntity);
  }
 
  template <class EntityType, int d = dim>
  std::enable_if_t<d == 3 && EntityType::codimension == 0, bool> contains(
      const EntityType& e) const {
    const auto hostEntity = e.impl().hostEntity();
    return grid_->getHostGrid().is_cell(hostEntity);
  }
 
  template <class EntityType, int d = dim>
  std::enable_if_t<d == 3 && EntityType::codimension == 1, bool> contains(
      const EntityType& e) const {
    const auto hostEntity = e.impl().hostEntity();
    return grid_->getHostGrid().tds().is_facet(hostEntity.first,
                                               hostEntity.second);
  }
 
  template <class EntityType, int d = dim>
  std::enable_if_t<d == 3 && EntityType::codimension == 2, bool> contains(
      const EntityType& e) const {
    const auto hostEntity = e.impl().hostEntity();
    return grid_->getHostGrid().tds().is_edge(
        hostEntity.first, hostEntity.second, hostEntity.third);
  }
 
  template <class EntityType, int d = dim>
  std::enable_if_t<d == 3 && EntityType::codimension == 3, bool> contains(
      const EntityType& e) const {
    const auto hostEntity = e.impl().hostEntity();
    return grid_->getHostGrid().tds().is_vertex(hostEntity);
  }
 
  template <int d = dim>
  std::enable_if_t<d == 2, void> update(const GridImp* grid) {
    grid_ = grid;
    const auto& hostgrid = grid_->getHostGrid();
 
    // Store face indices within face infos
    std::size_t elementCount = 0;
    for (const auto& element : elements(grid_->leafGridView(), Partitions::all))
      element.impl().hostEntity()->info().index = elementCount++;
 
    // Store vertex indices within vertex infos
    std::size_t vertexCount = 0;
    for (const auto& vertex : vertices(grid_->leafGridView(), Partitions::all))
      vertex.impl().hostEntity()->info().index = vertexCount++;
 
    // Store the finite edge indices in a map
    codimIndexMap_[0].clear();
    std::size_t edgeCount = 0;
    for (const auto& edge : edges(grid_->leafGridView(), Partitions::all))
      codimIndexMap_[0][grid_->globalIdSet().id(edge)] = edgeCount++;
 
    // Cache sizes since it is expensive to compute them
    sizeOfCodim_[0] = elementCount;
    sizeOfCodim_[1] = edgeCount;
    sizeOfCodim_[2] = vertexCount;
  }
 
  template <int d = dim>
  std::enable_if_t<d == 3, void> update(const GridImp* grid) {
    grid_ = grid;
    const auto& hostgrid = grid_->getHostGrid();
 
    // Store cell indices within cell infos
    std::size_t elementCount = 0;
    for (const auto& element : elements(grid_->leafGridView(), Partitions::all))
      element.impl().hostEntity()->info().index = elementCount++;
 
    // Store vertex indices within vertex infos
    std::size_t vertexCount = 0;
    for (const auto& vertex : vertices(grid_->leafGridView(), Partitions::all))
      vertex.impl().hostEntity()->info().index = vertexCount++;
 
    // Store the finite facet indices in a map
    codimIndexMap_[0].clear();
    std::size_t facetCount = 0;
    for (const auto& facet : facets(grid_->leafGridView(), Partitions::all))
      codimIndexMap_[0][grid_->globalIdSet().id(facet)] = facetCount++;
 
    // Store the finite edge indices in a map
    codimIndexMap_[1].clear();
    std::size_t edgeCount = 0;
    for (const auto& edge : edges(grid_->leafGridView(), Partitions::all))
      codimIndexMap_[1][grid_->globalIdSet().id(edge)] = edgeCount++;
 
    // Cache sizes since it is expensive to compute them
    sizeOfCodim_[0] = elementCount;
    sizeOfCodim_[1] = facetCount;
    sizeOfCodim_[2] = edgeCount;
    sizeOfCodim_[3] = vertexCount;
  }
 
  GridImp* grid_;
  std::array<std::size_t, dim + 1> sizeOfCodim_;
  std::array<CodimIndexMap, dim - 1> codimIndexMap_;
};
 
template <class GridImp>
class MMeshGlobalIdSet
    : public IdSet<GridImp, MMeshGlobalIdSet<GridImp>, MMeshImpl::MultiId> {
  typedef typename std::remove_const<GridImp>::type::HostGridType HostGrid;
 
  /*
   * We use the remove_const to extract the Type from the mutable class,
   * because the const class is not instantiated yet.
   */
  enum { dim = std::remove_const<GridImp>::type::dimension };
 
  template <int codim>
  using HostGridEntity = typename GridImp::template HostGridEntity<codim>;
 
 public:
  using IdType = MMeshImpl::MultiId;
 
  MMeshGlobalIdSet(const GridImp* g) : grid_(g), nextVertexId_(0) { init(); }
 
  void init() {
    const auto& hostgrid = grid_->getHostGrid();
 
    // Determine nextVertexId_
    for (auto vh = hostgrid.finite_vertices_begin();
         vh != hostgrid.finite_vertices_end(); ++vh)
      if (vh->info().idWasSet && vh->info().id >= nextVertexId_)
        nextVertexId_ = vh->info().id + 1;
  }
 
  /*
     We use the remove_const to extract the Type from the mutable class,
     because the const class is not instantiated yet.
   */
  template <int cd>
  std::enable_if_t<cd == dim, IdType> id(
      const typename std::remove_const<GridImp>::type::Traits::template Codim<
          cd>::Entity& e) const {
    return IdType(e.impl().hostEntity()->info().id);
  }
 
  template <int cd>
  std::enable_if_t<cd != dim, IdType> id(
      const typename std::remove_const<GridImp>::type::Traits::template Codim<
          cd>::Entity& e) const {
    return IdType(e.impl().id());
  }
 
  template <int d = dim>
  std::enable_if_t<d == 2, IdType> subId(
      const typename std::remove_const<GridImp>::type::Traits::template Codim<
          0>::Entity& e,
      int i, int codim) const {
    assert(0 <= codim && codim <= dim);
    IdType dummyId({std::size_t(-4), std::size_t(-3), std::size_t(-2)});
    switch (codim) {
      case 0:
        return e.impl().id();
      case 1:
        if (e.impl().id() != dummyId) {
          auto id0 = e.impl().id().vt()[i < 2 ? 0 : 1];
          auto id1 = e.impl().id().vt()[i == 0 ? 1 : 2];
          return IdType({std::min(id0, id1), std::max(id0, id1)});
        } else {
          std::size_t id0 = (i < 2 ? -4 : -3);
          std::size_t id1 = (i == 0 ? -3 : -2);
          return IdType({id0, id1});
        }
      case 2:
        if (e.impl().id() != dummyId)
          return e.impl().id().vt()[i];
        else
          return IdType(std::size_t(-4 + i));
    };
    return IdType();
  }
 
  template <int d = dim>
  std::enable_if_t<d == 3, IdType> subId(
      const typename std::remove_const<GridImp>::type::Traits::template Codim<
          0>::Entity& e,
      int i, int codim) const {
    assert(0 <= codim && codim <= dim);
    switch (codim) {
      case 0:
        return id<0>(e.impl().template subEntity<0>(i));
      case 1:
        return id<1>(e.impl().template subEntity<1>(i));
      case 2:
        return id<2>(e.impl().template subEntity<2>(i));
      case 3:
        return id<3>(e.impl().template subEntity<3>(i));
    };
    return IdType();
  }
 
  template <int d = dim>
  std::enable_if_t<d == 2, void> update(const GridImp* grid) {
    grid_ = grid;
    const auto& hostgrid = grid_->getHostGrid();
 
    // Store vertex ids within vertex infos
    for (auto vh = hostgrid.finite_vertices_begin();
         vh != hostgrid.finite_vertices_end(); ++vh)
      if (!vh->info().idWasSet) {
        vh->info().id = nextVertexId_++;
        vh->info().idWasSet = true;
      }
 
    // Compute mapping DUNE vertex index to CGAL vertex index
    for (auto fh = hostgrid.all_faces_begin(); fh != hostgrid.all_faces_end();
         ++fh)
      fh->info().cgalIndex = MMeshImpl::computeCGALIndices<decltype(fh), 2>(fh);
  }
 
  template <int d = dim>
  std::enable_if_t<d == 3, void> update(const GridImp* grid) {
    grid_ = grid;
    const auto& hostgrid = grid_->getHostGrid();
 
    // Store vertex ids within vertex infos
    for (auto vh = hostgrid.finite_vertices_begin();
         vh != hostgrid.finite_vertices_end(); ++vh)
      if (!vh->info().idWasSet) {
        vh->info().id = nextVertexId_++;
        vh->info().idWasSet = true;
      }
 
    // Compute mapping DUNE vertex index to CGAL vertex index
    for (auto ch = hostgrid.all_cells_begin(); ch != hostgrid.all_cells_end();
         ++ch)
      ch->info().cgalIndex = MMeshImpl::computeCGALIndices<decltype(ch), 3>(ch);
  }
 
  std::size_t setNextId(HostGridEntity<dim> vh) const {
    assert(!vh->info().idWasSet);
    vh->info().id = nextVertexId_++;
    vh->info().idWasSet = true;
    return vh->info().id;
  }
 
  GridImp* grid_;
  mutable std::size_t nextVertexId_;
};
 
}  // end namespace Dune
 
#endif