first commit in ahf_grid

dune/dec/AHFGrid/AHFGrid.hpp

+#pragma once
+
+#include <array>
+#include <cstddef>
+#include <list>
+#include <map>
+#include <set>
+#include <vector>
+
+#include <dune/common/fvector.hh>
+#include <dune/geometry/referenceelements.hh>
+
+#include <dune/dec/common/GeometryType.hpp>
+#include <dune/dec/common/Output.hpp>
+#include <dune/dec/ranges/IteratorRange.hpp>
+#include <dune/dec/utility/CapacityVector.hpp>
+
+namespace Dec
+{
+  template <int dim, int dow>
+  class AHFGrid
+  {
+    using ElementTopology = typename SimplexTopology<dim>::type;
+
+    template <int codim>
+    using SubTopology = typename SimplexTopology<dim - codim>::type;
+
+    struct HalfFacet
+    {
+      static const std::uint32_t invalid = std::uint32_t(-1);
+      static const std::uint8_t num_facets = ElementTopology::numFacets;
+      static const std::uint8_t border = num_facets;
+
+      std::uint32_t eid = invalid;  // element ID
+      std::uint8_t lfid = border;   // local facet ID
+
+      constexpr operator bool() const { return eid != invalid; }
+
+      bool outside() const { return lfid == border; }
+
+      HalfFacet next() const { return {eid, std::uint8_t((lfid + num_facets - 1) % num_facets)}; }
+      HalfFacet prev() const { return {eid, std::uint8_t((lfid + 1)   % num_facets)}; }
+    };
+
+
+    template <class T, class Super = std::vector<std::array<T, HalfFacet::num_facets+1>> >
+    struct HFVector
+        : public Super
+    {
+      HFVector(std::size_t num_elements = 0u)
+        : Super(num_elements)
+      {}
+
+      T const& operator[](HalfFacet const& hf) const { return Super::operator[](hf.eid)[hf.lfid]; }
+      T      & operator[](HalfFacet const& hf)       { return Super::operator[](hf.eid)[hf.lfid]; }
+    };
+
+    template <class T, class Super = std::vector<T> >
+    struct HFBorderVector
+        : public Super
+    {
+      HFBorderVector(std::size_t num_border = 0u)
+        : Super(num_border)
+      {}
+
+      T const& operator[](HalfFacet const& hf) const { return Super::operator[](hf.eid); }
+      T      & operator[](HalfFacet const& hf)       { return Super::operator[](hf.eid); }
+
+      using Super::operator[];
+
+      HalfFacet find(HalfFacet const& hf) const
+      {
+        if (hf.lfid != HalfFacet::border) {
+          error_exit("Invalid argument: expects a border half-facet.");
+          return {HalfFacet::invalid, HalfFacet::border};
+        }
+
+        for (std::uint32_t e_i = 0; e_i < this->size(); ++e_i)
+          if (Super::operator[](e_i).eid == hf.eid)
+            return {e_i, HalfFacet::border};
+
+        error_exit("half-facet not found.");
+        return {HalfFacet::invalid, HalfFacet::border};
+      }
+    };
+
+
+    using size_type = std::size_t;
+    using GlobalCoordinate = Dune::FieldVector<double, dow>;
+
+    using Elements = std::vector<std::vector<size_type>>;
+    using Coordinates = std::vector<GlobalCoordinate>;
+
+    using RefElement = Dune::ReferenceElement<double, dim>;
+    using RefElements = Dune::ReferenceElements<double, dim>;
+
+
+  public:
+
+    template <class E, class C>
+    void init(E const& elements, C const& coordinates);
+
+    // queries
+    // -------
+
+    size_type vertex_of_face (size_type f, std::uint8_t v_i) const { return elements_[f][v_i]; }
+    size_type edge_of_face   (size_type f, std::uint8_t e_i) const { return edge(hf(f,e_i)); }
+    size_type vertex_of_edge (size_type e, std::uint8_t v_i) const
+    {
+      auto const& hf = edge2hf_[e];
+      auto const& refElem = refElement();
+      return elements_[hf.eid][refElem.subEntity(e, dim-1, v_i, dim)];
+    }
+
+    size_type origin (HalfFacet const& hf) const
+    {
+      auto const& refElem = refElement();
+      return elements_[hf.eid][refElem.subEntity(hf.lfid, 1, 0, dim)];
+    }
+    size_type target (HalfFacet const& hf) const
+    {
+      auto const& refElem = refElement();
+      return elements_[hf.eid][refElem.subEntity(hf.lfid, 1, refElem.size(hf.lfid,1,dim)-1, dim)];
+    }
+
+    size_type face   (HalfFacet const& hf) const { return hf.eid; }
+    size_type edge   (HalfFacet const& hf) const { return hf2edge_[hf]; }
+    size_type vertex (HalfFacet const& hf) const { return target(hf); }
+
+    // cyclic_list<size_type>
+    auto edges_of_vertex(size_type v) const
+    {
+      auto hf = hf_of_vertex(v);
+      return ranges::make_cyclic_range(hf, hf,
+        [this](HalfFacet& hf) { hf = this->opposite(this->next(hf)); },
+        [this](HalfFacet const& hf) { return this->edge(hf); });
+    }
+
+    auto faces_of_vertex(size_type v) const
+    {
+      auto hf = hf_of_vertex(v);
+      return ranges::make_cyclic_range(hf, hf,
+        [this](HalfFacet& hf) {
+          auto face_hf = this->opposite(this->next(hf));
+          while (this->outside(face_hf))
+            face_hf = this->opposite(this->next(hf));
+          hf = face_hf;
+        },
+        [this](HalfFacet const& hf) { return this->face(hf); });
+    }
+
+    auto faces_of_edge(size_type e) const
+    {
+      auto hf = hf_of_edge(e);
+      CapacityVector<size_type,2> faces;
+      faces.push_back(this->face(hf));
+
+      auto neighbour = this->opposite(hf);
+      if (!this->outside(neighbour))
+        faces.push_back(neighbour);
+      return faces;
+    }
+
+
+    HalfFacet hf_of_vertex (size_type v) const { return vertex2hf_[v]; } // incident half-facet
+    HalfFacet hf_of_edge   (size_type e) const { return edge2hf_[e]; }
+    HalfFacet hf_of_face   (size_type f) const { return hf(f,0); }
+
+    HalfFacet hf (std::uint32_t f, std::uint8_t e_i) const { return {f,e_i}; }
+
+    HalfFacet opposite (HalfFacet const& hf) const { return hf2opposite_[hf]; }
+
+    HalfFacet next (HalfFacet const& hf) const
+    {
+      return outside(hf) ? hf2next_[hf] : hf.next();
+    }
+
+    HalfFacet prev (HalfFacet const& hf) const
+    {
+      return outside(hf) ? hf2next_.find(hf) : hf.prev();
+    }
+
+    bool outside (HalfFacet const& hf) const { return hf.outside(); }
+    bool outside (size_type v)         const { return vertex2hf_[v].outside(); }
+
+
+    RefElement const& refElement() const
+    {
+      return IsSimplex<ElementTopology>::value ? RefElements::simplex() : RefElements::cube();
+    }
+
+    std::size_t memory() const
+    {
+      std::size_t mem = 0;
+      mem += hf2opposite_.capacity()*(ElementTopology::numFacets+1)*sizeof(HalfFacet) + sizeof(hf2opposite_);
+      mem += hf2next_.capacity()*sizeof(HalfFacet) + sizeof(hf2next_);
+      mem += vertex2hf_.capacity()*sizeof(HalfFacet) + sizeof(vertex2hf_);
+      mem += edge2hf_.capacity()*sizeof(HalfFacet) + sizeof(edge2hf_);
+      mem += hf2edge_.capacity()*(ElementTopology::numFacets+1)*sizeof(size_type) + sizeof(hf2edge_);
+      mem += sizeof(std::uint32_t);
+      // mem += elements_.capacity()*(elements_[0].capacity()*sizeof(elements_[0][0]) + sizeof(elements_[0])) + sizeof(elements_);
+      // mem += coordinates_.capacity()*sizeof(coordinates_[0]) + sizeof(coordinates_);
+
+      return mem;
+    }
+
+  private:
+
+    HFVector<HalfFacet> hf2opposite_;
+    HFBorderVector<HalfFacet> hf2next_; // for border hf only
+
+    std::vector<HalfFacet> vertex2hf_; // half-edge pointing to vertex
+
+    // representation of facets (in 2d: edges)
+    std::vector<HalfFacet> edge2hf_;
+    HFVector<size_type> hf2edge_;
+
+    std::uint32_t numBorder_;
+
+    Elements elements_;
+    Coordinates coordinates_;
+  };
+
+} // end namespace Dec
+
+#include "AHFGrid.impl.hpp"

dune/dec/AHFGrid/AHFGrid.impl.hpp

+namespace Dec {
+
+template <int dim, int dow>
+  template <class E, class C>
+void AHFGrid<dim,dow>::init(E const& elements, C const& coordinates)
+{
+  elements_ = elements;
+  coordinates_ = coordinates;
+  std::size_t num_vertices = coordinates.size();
+
+  RefElement const& refElem = refElement();
+
+  // vertex of facet with largest ID
+  auto largest_vertex = [&elem=elements_, &refElem](std::uint32_t e_i, std::uint8_t f_i)
+  {
+    auto const& e = elem[e_i];
+    size_type v = e[refElem.subEntity(f_i, 1, 0, dim)];
+    int i0 = 0;
+    for (int i = 1; i < refElem.size(f_i,1,dim); ++i) {
+      int v_i = refElem.subEntity(f_i, 1, i, dim);
+      if (e[v_i] > v) {
+        v = e[v_i];
+        i0 = i;
+      }
+    }
+    return std::make_pair(v, i0);
+  };
+
+
+  { // steps 1. and 2.
+
+    // a mapping from each vertex to its incident half-facets in which the vertex has largest ID
+    std::vector<std::list<HalfFacet>> v2hfs(num_vertices);
+    // a mapping from each vertex to its adjacent vertices in each of the above incident half-facets
+    std::vector<std::map<HalfFacet,std::set<size_type>>> v2adj(num_vertices);
+
+    // 1. fill v2hfs and v2adj
+    for (std::uint32_t e_i = 0; e_i < elements_.size(); ++e_i) {
+      for (std::uint8_t f_i = 0; f_i < refElem.size(1); ++f_i) {
+        // auto const [v, i0] = largest_vertex(e_i, f_i);
+        size_type v = 0; int i0 = 0;
+        std::tie(v, i0) = largest_vertex(e_i, f_i);
+
+        auto hf = HalfFacet{e_i, f_i};
+        v2hfs[v].push_back(hf);
+
+        // set of adjacent vertices of v in facet (simplex only!)
+        auto& us = v2adj[v][hf];
+        for (int i = 0; i < refElem.size(f_i,1,dim); ++i)
+          if (i != i0)
+            us.insert(refElem.subEntity(f_i, 1, i, dim));
+      }
+    }
+
+    // 2. fill hf2opposite_
+    hf2opposite_.resize(elements_.size());
+    std::uint32_t border_id = 0;
+    for (std::uint32_t e_i = 0; e_i < elements_.size(); ++e_i) {
+      for (std::uint8_t f_i = 0; f_i < refElem.size(1); ++f_i) {
+        auto hf = HalfFacet{e_i, f_i};
+        if (!hf2opposite_[hf]) {
+          size_type v = 0; int i0 = 0;
+          std::tie(v, i0) = largest_vertex(e_i, f_i);
+
+          auto const& us = v2adj[v][hf];
+
+          std::vector<HalfFacet> hfs; hfs.reserve(2);
+          for (auto const& hf_ : v2hfs[v]) {
+            if (v2adj[v][hf_] == us)
+              hfs.push_back(hf_);
+          }
+
+          if (hfs.size() > 1) {
+            // form a cyclic mapping
+            for (std::size_t i = 0; i < hfs.size(); ++i)
+              hf2opposite_[hfs[i]] = hfs[(i+1)%hfs.size()];
+          } else {
+            // assign boundary half-edges
+            auto b = hf2opposite_[hfs[0]] = HalfFacet{border_id++, HalfFacet::border};
+            hf2opposite_[b] = hfs[0];
+          }
+        }
+      }
+    }
+
+    numBorder_ = border_id;
+  }
+
+  // 3. fill vertex2hf_
+  vertex2hf_.resize(num_vertices);
+  for (std::uint32_t e_i = 0; e_i < elements_.size(); ++e_i) {
+    auto const& e = elements_[e_i];
+    for (std::uint8_t f_i = 0; f_i < refElem.size(1); ++f_i) {
+      int v_i = refElem.subEntity(f_i, 1, refElem.size(f_i,1,dim)-1, dim); // target vertex
+      size_type v = e[v_i];
+      if (!vertex2hf_[v])
+        vertex2hf_[v] = HalfFacet{e_i, f_i};
+    }
+  }
+
+  { // 4. fill edge2hf_
+
+    // create unique ID for each facet
+    using Facet = SortedArray<size_type, SubTopology<1>::numCorners>;
+    std::map<Facet, std::pair<size_type, HalfFacet>> f2hf;
+    size_type facet_id = 0;
+    for (std::uint32_t e_i = 0; e_i < elements_.size(); ++e_i) {
+      auto const& e = elements_[e_i];
+      for (std::uint8_t f_i = 0; f_i < refElem.size(1); ++f_i) {
+        std::array<size_type, SubTopology<1>::numCorners> facet_unsorted;
+        for (int i = 0; i < refElem.size(f_i,1,dim); ++i)
+          facet_unsorted[i] = e[refElem.subEntity(f_i, 1, i, dim)];
+        auto facet = Facet(facet_unsorted);
+
+        auto hf = HalfFacet{e_i, f_i};
+        auto it = f2hf.insert(std::make_pair(facet, std::make_pair(facet_id, hf)));
+        if (it.second)
+          facet_id++;
+      }
+    }
+
+    edge2hf_.resize(f2hf.size());
+    // for (auto&& [key, facet] : f2hf) {
+    for (auto const& facet : f2hf) {
+      size_type id = facet.second.first;
+      edge2hf_[id] = facet.second.second;
+    }
+  }
+
+  // 5. fill hf2edge_
+  hf2edge_.resize(hf2opposite_.size());
+  for (size_type e_i = 0; e_i < edge2hf_.size(); ++e_i) {
+    auto const& hf = edge2hf_[e_i];
+    hf2edge_[hf] = e_i;
+    hf2edge_[hf2opposite_[hf]] = e_i;
+  }
+
+  // 6. fill next of border half-facets
+  hf2next_.resize(numBorder_);
+  for (std::uint32_t b_i = 0; b_i < hf2next_.size(); ++b_i) {
+    if (!hf2next_[b_i]) {
+      auto hf = opposite(HalfFacet{b_i, HalfFacet::border}).prev();
+      while (!opposite(hf).outside())
+        hf = opposite(hf).prev();
+
+      hf2next_[b_i] = opposite(hf);
+    }
+  }
+}
+
+#if 0
+
+template <int dim, int dow>
+  template <class E, class C>
+void AHFGrid<dim,dow>::init(E const& elements, C const& coordinates)
+{
+  elements_ = elements;
+  coordinates_ = coordinates;
+  std::size_t num_vertices = coordinates.size();
+
+  RefElement const& refElem = refElement();
+
+  hf2opposite_.resize(elements_.size());
+  hf2facet_.resize(elements.size());
+
+  // 1. create unique ID for each facet
+  using Facet = SortedArray<size_type, SubTopology<1>::numCorners>;
+  std::map<Facet, std::pair<size_type, HalfFacet>> f2hf;
+
+  size_type facet_id = 0;
+  for (std::uint32_t e_i = 0; e_i < elements_.size(); ++e_i) {
+    auto const& e = elements_[e_i];
+    for (std::uint8_t f_i = 0; f_i < refElem.size(1); ++f_i) {
+      std::array<size_type, SubTopology<1>::numCorners> facet_unsorted;
+      for (int i = 0; i < refElem.size(f_i,1,dim); ++i)
+        facet_unsorted[i] = e[refElem.subEntity(f_i, 1, i, dim)];
+      auto facet = Facet(facet_unsorted);
+
+      auto hf = HalfFacet{e_i, f_i};
+      auto it = f2hf.insert(std::make_pair(facet, std::make_pair(facet_id, hf)));
+      if (it.second)
+        facet_id++;
+      else
+        hf2opposite_[hf] = it.first->second.second;
+
+      hf2edge_[hf] = it.first->second.first;
+
+      size_type v = facet_unsorted.back();
+      if (!vertex2hf_[v])
+        vertex2hf_[v] = hf;
+    }
+  }
+
+  std::uint32_t border_id = 0;
+
+  // 2. assign hf to facet
+  edge2hf_.resize(f2hf.size());
+  for (auto const& facet : f2hf) {
+    size_type id = facet.second.first;
+    auto const& hf = facet.second.second;
+    edge2hf_[id] = hf;
+
+    auto& opp = hf2opposite_[hf];
+    if (!hf2opposite_[hf]) {
+      // assign boundary half-edges
+      opp = HalfFacet{border_id++, HalfFacet::border};
+      hf2opposite_[opp] = hf;
+    }
+  }
+
+  numBorder_ = border_id;
+
+  // 3. fill next of border half-facets
+  hf2next_.resize(numBorder_);
+  for (std::uint32_t b_i = 0; b_i < hf2next_.size(); ++b_i) {
+    if (!hf2next_[b_i]) {
+      auto hf = opposite(HalfFacet{b_i, HalfFacet::border}).prev();
+      while (!opposite(hf).outside())
+        hf = opposite(hf).prev();
+
+      hf2next_[b_i] = opposite(hf);
+    }
+  }
+}
+
+#endif
+
+} // end namespace Dec

dune/dec/common/GeometryType.hpp

+#pragma once
+
+#include <string>
+#include <type_traits>
+
+namespace Dec
+{
+  enum TopologyConstruction { pyramidConstruction = 0, prismConstruction = 1 };
+
+  // Basic Topology Types
+  // --------------------
+
+  struct Point
+  {
+    static const unsigned int dimension = 0;
+    static const unsigned int numCorners = 1;
+    static const unsigned int numFacets = 0;
+
+    static const unsigned int id = 0;
+
+    static std::string name () { return "p"; }
+  };
+
+
+  template< class BaseTopology >
+  struct Prism
+  {
+    static const unsigned int dimension = BaseTopology::dimension + 1;
+    static const unsigned int numCorners = 2 * BaseTopology::numCorners;
+    static const unsigned int numFacets = BaseTopology::numFacets + 2;
+
+    static const unsigned int id = BaseTopology::id | ((unsigned int)prismConstruction << (dimension-1));
+
+    static std::string name () { return BaseTopology::name() + "l"; }
+  };
+
+
+  template< class BaseTopology >
+  struct Pyramid
+  {
+    static const unsigned int dimension = BaseTopology::dimension + 1;
+    static const unsigned int numCorners = BaseTopology::numCorners + 1;
+    static const unsigned int numFacets = dimension == 1 ? 2 : BaseTopology::numFacets + 1;
+
+    static const unsigned int id = BaseTopology::id | ((unsigned int)pyramidConstruction << (dimension-1));
+
+    static std::string name () { return BaseTopology::name() + "o"; }
+  };
+
+
+
+  // Properties of Topologies
+  // ------------------------
+
+  template< class Topology >
+  struct IsSimplex
+      : public std::integral_constant< bool, (Topology::id >> 1) == 0 >
+  {};
+
+  template< class Topology >
+  struct IsCube
+      : public std::integral_constant< bool,  (Topology::id | 1) == (1 << Topology::dimension) - 1 >
+  {};
+
+
+
+  // SimplexTopology
+  // ---------------
+
+  template< unsigned int dim >
+  struct SimplexTopology
+  {
+    using type = Pyramid< typename SimplexTopology< dim-1 >::type >;
+  };
+
+  template<>
+  struct SimplexTopology< 0 >
+  {
+    using type = Point;
+  };
+
+
+
+  // CubeTopology
+  // ------------
+
+  template< unsigned int dim >
+  struct CubeTopology
+  {
+    using type = Prism< typename CubeTopology< dim-1 >::type >;
+  };
+
+  template<>
+  struct CubeTopology< 0 >
+  {
+    using type = Point;
+  };
+
+
+
+  // PyramidTopology
+  // ---------------
+
+  template< unsigned int dim >
+  struct PyramidTopology
+  {
+    using type = Pyramid< typename CubeTopology< dim-1 >::type >;
+  };
+
+
+
+  // PrismTopology
+  // -------------
+
+  template< unsigned int dim >
+  struct PrismTopology
+  {
+    using type = Prism< typename SimplexTopology< dim-1 >::type >;
+  };
+
+} // end namespace Dec

dune/dec/halfedgegrid/HalfEdgeGrid.hpp

+#pragma once
+
+#include "HalfEdge.hpp"
+
+namespace Dec
+{
+  template <int dim, int dow>
+  class HalfEdgeGrid
+  {
+  public: // initialization of data structures:
+
+    /// Create half-edge structure based on the given Dune `gridView`.
+    template <class BaseGridView>
+    void init(BaseGridView const& gridView);
+
+
+  private: // implementation details:
+
+    // add a chain of boundary half-edges, starting with a new half-edge with index `he`
+    // at the opposite of `next(he_start)`, `he_last` is the last inserted element of the chain.
+    IndexType add_boundary_chain(IndexType he, IndexType he_start, IndexType he_last);
+
+
+  public:
+
+    /// Returns the half-edge with index `he`
+    HalfEdge const& halfEdge(IndexType he) const
+    {
+      assert(