adaptive_mesh_storage.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
#include "kernel/util/string.hpp"
#include "kernel/util/tiny_algebra.hpp"
#include <iostream>
#include <kernel/shape.hpp>
#include <kernel/util/assertion.hpp>
#include <kernel/util/slotmap.hpp>
 
#include <array>
#include <cstdint>
#include <utility>
 
namespace FEAT::Geometry::Intern
{
 
  struct Layer
  {
    Index idx;
  };
 
  inline bool operator==(Layer lhs, Layer rhs)
  {
    return lhs.idx == rhs.idx;
  }
 
  inline bool operator!=(Layer lhs, Layer rhs)
  {
    return !(lhs == rhs);
  }
 
  inline bool operator<(Layer lhs, Layer rhs)
  {
    return lhs.idx < rhs.idx;
  }
 
  inline bool operator>(Layer lhs, Layer rhs)
  {
    return lhs.idx > rhs.idx;
  }
 
  inline bool operator<=(Layer lhs, Layer rhs)
  {
    return lhs.idx <= rhs.idx;
  }
 
  inline bool operator>=(Layer lhs, Layer rhs)
  {
    return lhs.idx >= rhs.idx;
  }
 
  template<int n>
  struct ElementKey
  {
    ElementKey() = default;
    ElementKey(std::uint32_t idx, std::uint32_t gen) : index(idx), generation(gen)
    {
    }
 
    ElementKey(const ElementKey& other) = default;
    ElementKey(ElementKey&& other) = default;
 
    ElementKey& operator=(const ElementKey& other) = default;
    ElementKey& operator=(ElementKey&& other) = default;
 
    std::uint32_t index = 0;
    std::uint32_t generation = 0;
 
    Layer layer = {0};
    bool is_permanent = false;
  };
 
  template<int dim_>
  inline bool operator==(const ElementKey<dim_>& lhs, const ElementKey<dim_>& rhs)
  {
    return lhs.index == rhs.index && lhs.generation == rhs.generation && lhs.layer == rhs.layer &&
           lhs.is_permanent == rhs.is_permanent;
  }
 
  template<int dim_>
  inline bool operator!=(const ElementKey<dim_>& lhs, const ElementKey<dim_>& rhs)
  {
    return !(lhs == rhs);
  }
 
  template<int dim_>
  inline std::ostream& operator<<(std::ostream& stream, const ElementKey<dim_> key)
  {
    stream
      << "ElementKey<" << stringify(dim_)
      << ">{ index = " << stringify(key.index)
      << ", generation = " << stringify(key.generation)
      << ", layer = Layer {" << stringify(key.layer.idx)
      << "}, is_permanent = " << stringify(key.is_permanent) << " }";
    return stream;
  }
 
  using VertexKey = ElementKey<0>;
  using EdgeKey = ElementKey<1>;
  using FaceKey = ElementKey<2>;
  using CellKey = ElementKey<3>;
 
  template<int n>
  struct OrientedElement
  {
    OrientedElement() = default;
    OrientedElement(int o, ElementKey<n> k) : orientation(o), key(k)
    {
    }
 
    int orientation = 0;
    ElementKey<n> key = ElementKey<n>();
  };
 
  using OrientedEdge = OrientedElement<1>;
  using OrientedFace = OrientedElement<2>;
 
  template<int dim_>
  inline bool operator==(const OrientedElement<dim_>& lhs, const OrientedElement<dim_>& rhs)
  {
    return lhs.key == rhs.key && lhs.orientation == rhs.orientation;
  }
 
  template<int dim_>
  inline bool operator!=(const OrientedElement<dim_>& lhs, const OrientedElement<dim_>& rhs)
  {
    return !(lhs == rhs);
  }
 
  template<typename TemplateSet_, typename Shape_, int dim_ = Shape_::dimension>
  struct ElementChildren : public ElementChildren<TemplateSet_, Shape_, dim_ - 1>
  {
    static constexpr int max_children = TemplateSet_::template max_children<Shape_::dimension, dim_>();
 
    std::array<Intern::ElementKey<dim_>, max_children> children;
 
    template<int query_dim_>
    std::array<ElementKey<query_dim_>, TemplateSet_::template max_children<Shape_::dimension, query_dim_>()>& by_dim()
    {
      static_assert(query_dim_ <= Shape_::dimension);
 
      return ElementChildren<TemplateSet_, Shape_, query_dim_>::children;
    }
 
    template<int query_dim_>
    const std::array<ElementKey<query_dim_>, TemplateSet_::template max_children<Shape_::dimension, query_dim_>()>& by_dim() const
    {
      static_assert(query_dim_ <= Shape_::dimension);
 
      return ElementChildren<TemplateSet_, Shape_, query_dim_>::children;
    }
  };
 
  template<typename TemplateSet_, typename Shape_>
  struct ElementChildren<TemplateSet_, Shape_, 0>
  {
    static constexpr int max_children = TemplateSet_::template max_children<Shape_::dimension, 0>();
 
    std::array<Intern::ElementKey<0>, max_children> children;
 
    template<int query_dim_>
    std::array<ElementKey<query_dim_>, TemplateSet_::template max_children<Shape_::dimension, query_dim_>()>& by_dim()
    {
      static_assert(query_dim_ == 0);
 
      return children;
    }
 
    template<int query_dim_>
    const std::array<ElementKey<query_dim_>, TemplateSet_::template max_children<Shape_::dimension, query_dim_>()>& by_dim() const
    {
      static_assert(query_dim_ == 0);
 
      return children;
    }
  };
 
  template<typename Shape_, int dim_ = Shape_::dimension - 1>
  struct ElementTopology : public ElementTopology<Shape_, dim_ - 1>
  {
    static_assert(dim_ >= 0);
 
    static constexpr int num_entities = Shape::FaceTraits<Shape_, dim_>::count;
 
    std::array<OrientedElement<dim_>, num_entities> entities;
 
    template<int key_dim_>
    OrientedElement<key_dim_> key_by_dim(int idx) const
    {
      static_assert(key_dim_ < Shape_::dimension);
      ASSERT(idx < (ElementTopology<Shape_, key_dim_>::num_entities));
 
      return ElementTopology<Shape_, key_dim_>::entities[std::size_t(idx)];
    }
 
    template<int query_dim_>
    std::array<OrientedElement<query_dim_>, Shape::FaceTraits<Shape_, query_dim_>::count>& by_dim()
    {
      static_assert(query_dim_ < Shape_::dimension);
 
      return ElementTopology<Shape_, query_dim_>::entities;
    }
 
    template<int query_dim_>
    const std::array<OrientedElement<query_dim_>, Shape::FaceTraits<Shape_, query_dim_>::count>& by_dim() const
    {
      static_assert(query_dim_ < Shape_::dimension);
 
      return ElementTopology<Shape_, query_dim_>::entities;
    }
 
    template<int query_dim_>
    constexpr int size()
    {
      static_assert(query_dim_ <= Shape_::dimension);
 
      return Shape::FaceTraits<Shape_, query_dim_>::count;
    }
  };
 
  template<typename Shape_>
  struct ElementTopology<Shape_, 0>
  {
    static constexpr int num_entities = Shape::FaceTraits<Shape_, 0>::count;
 
    std::array<OrientedElement<0>, num_entities> entities;
 
    template<int key_dim_>
    OrientedElement<key_dim_> key_by_dim(int idx) const
    {
      static_assert(key_dim_ == 0);
      ASSERT(idx < num_entities);
 
      return entities[std::size_t(idx)];
    }
 
    template<int query_dim_>
    std::array<OrientedElement<query_dim_>, Shape::FaceTraits<Shape_, query_dim_>::count>& by_dim()
    {
      static_assert(query_dim_ == 0);
 
      return entities;
    }
 
    template<int query_dim_>
    const std::array<OrientedElement<query_dim_>, Shape::FaceTraits<Shape_, query_dim_>::count>& by_dim() const
    {
      static_assert(query_dim_ == 0);
 
      return entities;
    }
 
    template<int query_dim_>
    constexpr int size()
    {
      static_assert(query_dim_ == 0);
 
      return Shape::FaceTraits<Shape_, query_dim_>::count;
    }
  };
 
  template<typename TemplateSet_, typename Shape_>
  struct AdaptiveElement
  {
    static_assert(Shape_::dimension > 0);
 
    using TemplateSet = TemplateSet_;
 
    using RefinementType = typename TemplateSet::template RefinementTypeByDim<Shape_::dimension>;
 
    AdaptiveElement(RefinementType t, Layer l, const ElementTopology<Shape_> topo) :
      type(t),
      layer(l),
      topology(std::move(topo))
    {
    }
 
    RefinementType type;
    Index index = 0;
    Layer layer = {0};
 
    // OPTIMIZATION: We are allocating the space for the maximum number of
    // children of any template of the template set. We can instead store
    // children out-of-band in Util::SecondaryMaps with two advantages:
    // * we do not need to allocate space for children of zero type elements
    // * improved cache-behaviour when iterating over mesh elements, since most
    // algorithms do not need the parent-child relations
    ElementChildren<TemplateSet_, Shape_> children;
    ElementTopology<Shape_> topology;
 
    static constexpr int num_neighbors = Shape::FaceTraits<Shape_, Shape_::dimension - 1>::count;
    std::array<Intern::ElementKey<Shape_::dimension>, num_neighbors> neighbors;
  };
 
  template<typename VertexType>
  struct AdaptiveVertex
  {
    AdaptiveVertex(VertexType v, Layer l) : vertex(v), layer(l), last_changed(l)
    {
    }
 
    VertexType vertex;
    Index index = 0;
    Layer layer = {0};
 
    // NOTE: last_changed is (so far) exclusively used for the
    // BPXJacobiPreconditioner. The mesh should probably not directly know that this
    // property is needed. Maybe introduce something like the PMP property system
    // that allows dynamically adding data to any mesh entities.
 
    Layer last_changed;
  };
 
  template<typename TemplateSet_, typename MeshShape_, typename VertexType_, int element_dim_ = MeshShape_::dimension>
  class MeshLayer : public MeshLayer<TemplateSet_, MeshShape_, VertexType_, element_dim_ - 1>
  {
  protected:
    using BaseClass = MeshLayer<TemplateSet_, MeshShape_, VertexType_, element_dim_>;
 
    using ShapeType = typename Shape::FaceTraits<MeshShape_, element_dim_>::ShapeType;
 
    using AdaptiveElementType = AdaptiveElement<TemplateSet_, ShapeType>;
    using ElementKeyType = ElementKey<element_dim_>;
 
    Util::SlotMap<AdaptiveElementType, ElementKeyType> _permanent_entities;
    Util::SlotMap<AdaptiveElementType, ElementKeyType> _transient_entities;
 
  public:
    template<typename Shape_>
    ElementKey<Shape_::dimension> insert_permanent(const AdaptiveElement<TemplateSet_, Shape_>& element)
    {
      auto key = MeshLayer<TemplateSet_, MeshShape_, VertexType_, Shape_::dimension>::_permanent_entities.insert(element);
      key.is_permanent = true;
      return key;
    }
 
    template<typename Shape_>
    ElementKey<Shape_::dimension> insert_transient(const AdaptiveElement<TemplateSet_, Shape_>& element)
    {
      auto key = MeshLayer<TemplateSet_, MeshShape_, VertexType_, Shape_::dimension>::_transient_entities.insert(element);
      key.is_permanent = false;
      return key;
    }
 
    template<int dim_>
    void erase(ElementKey<dim_> key)
    {
      if constexpr(dim_ == 0)
      {
        MeshLayer<TemplateSet_, MeshShape_, VertexType_, 0>::_vertices.erase(key);
      }
      else
      {
        if(key.is_permanent)
        {
          MeshLayer<TemplateSet_, MeshShape_, VertexType_, dim_>::_permanent_entities.erase(key);
        }
        else
        {
          MeshLayer<TemplateSet_, MeshShape_, VertexType_, dim_>::_transient_entities.erase(key);
        }
      }
    }
 
    template<int dim_>
    Index num_permanent_entities() const
    {
      if constexpr(dim_ == 0)
      {
        return MeshLayer<TemplateSet_, MeshShape_, VertexType_, 0>::_vertices.size();
      }
      else
      {
        return MeshLayer<TemplateSet_, MeshShape_, VertexType_, dim_>::_permanent_entities.size();
      }
    }
 
    template<int dim_>
    Index num_transient_entities() const
    {
      if constexpr(dim_ == 0)
      {
        return 0;
      }
      else
      {
        return MeshLayer<TemplateSet_, MeshShape_, VertexType_, dim_>::_transient_entities.size();
      }
    }
 
    template<int dim_>
    Util::SlotMap<AdaptiveElement<TemplateSet_, typename Shape::FaceTraits<MeshShape_, dim_>::ShapeType>, ElementKey<dim_>>&
    permanent_entities()
    {
      return MeshLayer<TemplateSet_, MeshShape_, VertexType_, dim_>::_permanent_entities;
    }
 
    template<int dim_>
    const Util::SlotMap<AdaptiveElement<TemplateSet_, typename Shape::FaceTraits<MeshShape_, dim_>::ShapeType>, ElementKey<dim_>>&
    permanent_entities() const
    {
      return MeshLayer<TemplateSet_, MeshShape_, VertexType_, dim_>::_permanent_entities;
    }
 
    template<int dim_>
    Util::SlotMap<AdaptiveElement<TemplateSet_, typename Shape::FaceTraits<MeshShape_, dim_>::ShapeType>, ElementKey<dim_>>&
    transient_entities()
    {
      return MeshLayer<TemplateSet_, MeshShape_, VertexType_, dim_>::_transient_entities;
    }
 
    template<int dim_>
    const Util::SlotMap<AdaptiveElement<TemplateSet_, typename Shape::FaceTraits<MeshShape_, dim_>::ShapeType>, ElementKey<dim_>>&
    transient_entities() const
    {
      return MeshLayer<TemplateSet_, MeshShape_, VertexType_, dim_>::_transient_entities;
    }
 
    std::size_t bytes() const
    {
      return _permanent_entities.bytes() + _transient_entities.bytes() + BaseClass::bytes();
    }
 
    template<int dim_>
    AdaptiveElement<TemplateSet_, typename Shape::FaceTraits<MeshShape_, dim_>::ShapeType>& get(ElementKey<dim_> key)
    {
      static_assert(dim_ >= 1);
      if(key.is_permanent)
      {
        return MeshLayer<TemplateSet_, MeshShape_, VertexType_, dim_>::_permanent_entities[key];
      }
      else
      {
        return MeshLayer<TemplateSet_, MeshShape_, VertexType_, dim_>::_transient_entities[key];
      }
    }
 
    template<int dim_>
    const AdaptiveElement<TemplateSet_, typename Shape::FaceTraits<MeshShape_, dim_>::ShapeType>& get(ElementKey<dim_> key) const
    {
      static_assert(dim_ >= 1);
      if(key.is_permanent)
      {
        return MeshLayer<TemplateSet_, MeshShape_, VertexType_, dim_>::_permanent_entities[key];
      }
      else
      {
        return MeshLayer<TemplateSet_, MeshShape_, VertexType_, dim_>::_transient_entities[key];
      }
    }
  };
 
  template<typename TemplateSet_, typename MeshShape_, typename VertexType_>
  class MeshLayer<TemplateSet_, MeshShape_, VertexType_, 0>
  {
  protected:
    Util::SlotMap<AdaptiveVertex<VertexType_>, ElementKey<0>> _vertices;
 
  public:
    ElementKey<0> insert_vertex(const AdaptiveVertex<VertexType_>& vertex)
    {
      return _vertices.insert(vertex);
    }
 
    void erase(ElementKey<0> key)
    {
      _vertices.erase(key);
    }
 
    Index num_vertices() const
    {
      return _vertices.size();
    }
 
    Util::SlotMap<AdaptiveVertex<VertexType_>, ElementKey<0>>& vertices()
    {
      return _vertices;
    }
 
    const Util::SlotMap<AdaptiveVertex<VertexType_>, ElementKey<0>>& vertices() const
    {
      return _vertices;
    }
 
    std::size_t bytes() const
    {
      return _vertices.bytes();
    }
 
    AdaptiveVertex<VertexType_>& get_vertex(ElementKey<0> key)
    {
      return _vertices[key];
    }
 
    const AdaptiveVertex<VertexType_>& get_vertex(ElementKey<0> key) const
    {
      return _vertices[key];
    }
  };
 
  template<typename MeshShape_, typename TemplateSet_, typename VertexType_>
  class AdaptiveMeshStorage
  {
  public:
    using VertexType = VertexType_;
 
    using CoordType = typename VertexType_::ValueType;
 
    using TemplateSet = TemplateSet_;
 
    template<int dim_>
    using ElementRefByDim = ElementKey<dim_>;
 
    template<int dim_>
    using OrientedElementRefByDim = OrientedElement<dim_>;
 
    template<int dim_>
    using ElementTopologyByDim = ElementTopology<typename Shape::FaceTraits<MeshShape_, dim_>::ShapeType>;
 
    template<int dim_>
    using ElementChildrenByDim = ElementChildren<TemplateSet_, typename Shape::FaceTraits<MeshShape_, dim_>::ShapeType>;
 
    template<int dim_>
    using AdaptiveElementByDim = AdaptiveElement<TemplateSet_, typename Shape::FaceTraits<MeshShape_, dim_>::ShapeType>;
 
    using AdaptiveVertexType = AdaptiveVertex<VertexType>;
 
    enum class IterationOrder
    {
      PreOrder,
      PostOrder
    };
 
  private:
    using MeshLayerType = MeshLayer<TemplateSet, MeshShape_, VertexType>;
 
    static constexpr int num_coords = VertexType::n;
 
    static const constexpr ElementKey<MeshShape_::dimension> neighbor_sentinel = ElementKey<MeshShape_::dimension>();
 
    std::vector<MeshLayerType> _layers;
 
    bool _is_indexed = false;
 
    bool _has_neighbors = false;
 
  public:
    template<int dim>
    typename TemplateSet::template RefinementTypeByDim<dim> type(ElementKey<dim> key)
    {
      return (*this)[key].type;
    }
 
    template<int dim_>
    std::array<Index, 4> erase(ElementKey<dim_> key)
    {
      ASSERTM(key.layer.idx < _layers.size(), "Trying to erase adaptive mesh element of non-existing layer");
 
      // Erasing an entity produces a gap in that layers indices. The mesh is
      // thus no longer correctly indexed.
      _is_indexed = false;
      // Erasing an entity (most-likely) removed a neighbor of an element. Neighbors need to be re-determined.
      _has_neighbors = false;
 
      std::array<Index, 4> num_erased = {};
      if constexpr(dim_ > 0)
      {
        // NOTE: Different mesh layers are stored separately. Deleting a child
        // does thus not invalidate the references to the current entity or its
        // children.
        auto& entity = (*this)[key];
        num_erased =
          _erase_children<typename Shape::FaceTraits<MeshShape_, dim_>::ShapeType, dim_>(entity.type, entity.children);
      }
 
      // Actually delete current entity
      _layers[key.layer.idx].erase(key);
      num_erased[dim_] += 1;
 
      return num_erased;
    }
 
    VertexKey insert(const AdaptiveVertexType& vertex)
    {
      // Adding a new vertex to a layer shifts the indices of all vertices on
      // further layers. These vertices stored indices are thus wrong, and need
      // to be recomputed.
      _is_indexed = false;
 
      // Add additional mesh layers, if required
      Layer layer = vertex.layer;
      _add_mesh_layers(layer.idx);
 
      // Add vertex
      VertexKey key = _layers[layer.idx].insert_vertex(vertex);
      key.layer = layer;
      return key;
    }
 
    VertexKey insert(VertexType vertex, Layer layer)
    {
      return insert(AdaptiveVertex(vertex, layer));
    }
 
    template<typename Shape_>
    ElementKey<Shape_::dimension> insert(const AdaptiveElement<TemplateSet, Shape_>& element)
    {
      // Adding an entity might cause the indices of entities on further layers
      // to be shifted. These entites stored indices are thus wrong and need to
      // be recomputed.
      _is_indexed = false;
      // Neighbors for this new entity need to be determined.
      _has_neighbors = false;
 
      // Add additional mesh layers, if required
      Layer layer = element.layer;
      _add_mesh_layers(layer.idx);
 
      ElementKey<Shape_::dimension> key;
      bool is_permanent = element.type.is_zero_refinement();
 
      if(is_permanent)
      {
        key = _layers[layer.idx].insert_permanent(element);
      }
      else
      {
        key = _layers[layer.idx].insert_transient(element);
      }
 
      // Add metadata to key
      key.is_permanent = is_permanent;
      key.layer = layer;
 
      return key;
    }
 
    Index num_layers() const
    {
      return _layers.size();
    }
 
    template<typename VisitorType_, int dim_>
    void walk_subtree(
      ElementKey<dim_> root_key,
      VisitorType_& visitor,
      Layer max_depth,
      IterationOrder order = IterationOrder::PreOrder) const
    {
      auto& root = (*this)[root_key];
 
      // Abort if max_depth is reached
      if(root.layer > max_depth)
      {
        return;
      }
 
      // Visit Root
      if(order == IterationOrder::PreOrder)
      {
        visitor(root_key);
      }
 
      if constexpr(dim_ > 0)
      {
        _walk_subtree_children<typename Shape::FaceTraits<MeshShape_, dim_>::ShapeType>(root, visitor, max_depth);
      }
 
      // Visit Root
      if(order == IterationOrder::PostOrder)
      {
        visitor(root_key);
      }
    }
 
    template<typename VisitorType_>
    void walk_tree(VisitorType_& visitor, Layer max_depth, IterationOrder order = IterationOrder::PreOrder) const
    {
      // There is nothing to iterate over if the mesh has no layers
      if(_layers.size() == 0)
      {
        return;
      }
 
      const MeshLayerType& root_layer = _layers[0];
 
      if constexpr(MeshShape_::dimension >= 0)
      {
        for(auto entry : root_layer.vertices())
        {
          walk_subtree(entry.key, visitor, max_depth, order);
        }
      }
 
      if constexpr(MeshShape_::dimension >= 1)
      {
        for(auto entry : root_layer.template permanent_entities<1>())
        {
          walk_subtree(entry.key, visitor, max_depth, order);
        }
        for(auto entry : root_layer.template transient_entities<1>())
        {
          walk_subtree(entry.key, visitor, max_depth, order);
        }
      }
 
      if constexpr(MeshShape_::dimension >= 2)
      {
        for(auto entry : root_layer.template permanent_entities<2>())
        {
          walk_subtree(entry.key, visitor, max_depth, order);
        }
        for(auto entry : root_layer.template transient_entities<2>())
        {
          walk_subtree(entry.key, visitor, max_depth, order);
        }
      }
 
      if constexpr(MeshShape_::dimension >= 3)
      {
        for(auto entry : root_layer.template permanent_entities<3>())
        {
          walk_subtree(entry.key, visitor, max_depth, order);
        }
        for(auto entry : root_layer.template transient_entities<3>())
        {
          walk_subtree(entry.key, visitor, max_depth, order);
        }
      }
    }
 
    Index num_entities(Layer layer, int dim) const
    {
      if(dim == 0)
      {
        return _num_entities<0>(layer);
      }
      if(dim == 1)
      {
        return _num_entities<1>(layer);
      }
      if(dim == 2)
      {
        return _num_entities<2>(layer);
      }
      if(dim == 3)
      {
        return _num_entities<3>(layer);
      }
      XABORTM("unsupported dimension in MeshStorage::num_entities");
    }
 
    template<int dim_>
    Index num_total_entities() const
    {
      Index result = 0;
      for(auto& layer : _layers)
      {
        result += layer.template num_permanent_entities<dim_>();
        result += layer.template num_transient_entities<dim_>();
      }
 
      return result;
    }
 
    void fill_neighbors()
    {
      static constexpr int shape_dim = MeshShape_::dimension;
      static constexpr int num_facets = Shape::FaceTraits<MeshShape_, shape_dim - 1>::count;
 
      // OPTIMIZATION: This is the ConformalMesh::fill_neighbors() algorithm
      // ported to adaptive mesh keys. This information can be retrieved more
      // cheaply from the refinement templates themselves. Add this capability to
      // the template query (maybe optionally, with a boolean indicating support
      // on the template set), support it in the AdaptiveMesh code, and don't
      // invalidate neighbors when adding elements here if the template provide
      // neighbor information.
      // This method can then be left as a fallback.
 
      if(_has_neighbors)
      {
        // Neighbors have already been determined since the last adaptation of
        // the mesh. No work to be done.
        return;
      }
 
      std::cout << "Recalculating neighbors\n";
 
      // In the following we are calling mesh entities of dimension shape_dim - 1 facets.
      // We are using the fact that each facet is shared by at most 2 elements.
      // We are goint to first determine which elements are adjacent to each
      // facet, by iterating over the elements and adding their keys to the
      // facet_neighbors of each of their facets. Then we can iterate over the
      // elements a second time and determine their neighbors as the element
      // registered for the facet, that is not itself.
 
      // Maps of facet neighbors
      std::vector<std::vector<Util::SecondaryMap<std::array<ElementKey<shape_dim>, 2>, ElementKey<shape_dim - 1>>>>
        permanent_facet_neighbors(num_layers());
 
      std::vector<std::vector<Util::SecondaryMap<std::array<ElementKey<shape_dim>, 2>, ElementKey<shape_dim - 1>>>>
        transient_facet_neighbors(num_layers());
 
      for(Index i(0); i < num_layers(); i++)
      {
        permanent_facet_neighbors[i].resize(num_layers());
        transient_facet_neighbors[i].resize(num_layers());
      }
 
      // Register element with their adjacent facets. Fills the facet_neighbors map.
      auto determine_facet_neighbors = [&](Layer layer, auto slotmap_entry)
      {
        // HACK: We have to recreate the "full" element key here.
        // The keys returned by the SlotMap iterators do not contain the
        // layer and permanency information we have added to the usual
        // slotmap data.
        // The proper way to this is to wrap the iterators in the MeshLayer
        // class and add this information there. As this is currently the
        // only place where we directly iterate over all elements of a layer
        // via the iterators, I am leaving that boilerplate out for now.
        ElementKey<shape_dim> element_key = slotmap_entry.key;
        element_key.layer = layer;
        element_key.is_permanent = slotmap_entry.value.type.is_zero_refinement();
 
        // Iterate over all facets
        for(int facet_idx = 0; facet_idx < num_facets; facet_idx++)
        {
          // Retrieve key of facet
          auto key = slotmap_entry.value.topology.template key_by_dim<shape_dim - 1>(facet_idx).key;
 
          // Init secondary map, if this is the first time we see this facet
          if(key.is_permanent && !permanent_facet_neighbors[layer.idx][key.layer.idx].contains_key(key))
          {
            permanent_facet_neighbors[layer.idx][key.layer.idx].insert(key, {neighbor_sentinel, neighbor_sentinel});
          }
 
          if(!key.is_permanent && !transient_facet_neighbors[layer.idx][key.layer.idx].contains_key(key))
          {
            transient_facet_neighbors[layer.idx][key.layer.idx].insert(key, {neighbor_sentinel, neighbor_sentinel});
          }
 
          auto& neighbors = key.is_permanent ? permanent_facet_neighbors[layer.idx][key.layer.idx][key] : transient_facet_neighbors[layer.idx][key.layer.idx][key];
 
          if(neighbors[0] == neighbor_sentinel)
          {
            neighbors[0] = element_key;
          }
          else if(neighbors[1] == neighbor_sentinel)
          {
            neighbors[1] = element_key;
          }
          else
          {
            XABORTM("Facet has more than two neighbors!");
          }
        }
      };
 
      // Determine neighbors by looking at registered facet neighbors.
      auto determine_neighbors = [&](Layer layer, auto slotmap_entry)
      {
        auto& element = slotmap_entry.value;
 
        // Iterate over all facets
        for(int facet_idx = 0; facet_idx < num_facets; facet_idx++)
        {
          // Retrieve key of facet
          auto key = slotmap_entry.value.topology.template key_by_dim<shape_dim - 1>(facet_idx).key;
 
          // HACK: We have to recreate the "full" element key here.
          // The keys returned by the SlotMap iterators do not contain the
          // layer and permanency information we have added to the usual
          // slotmap data.
          // The proper way to this is to wrap the iterators in the MeshLayer
          // class and add this information there. As this is currently the
          // only place where we directly iterate over all elements of a layer
          // via the iterators, I am leaving that boilerplate out for now.
          ElementKey<shape_dim> element_key = slotmap_entry.key;
          element_key.layer = layer;
          element_key.is_permanent = slotmap_entry.value.type.is_zero_refinement();
 
          // Neighbor is the entity that is not the current entity
          auto& candidates = key.is_permanent ? permanent_facet_neighbors[layer.idx][key.layer.idx][key] : transient_facet_neighbors[layer.idx][key.layer.idx][key];
          if(candidates[0] == element_key)
          {
            element.neighbors[std::size_t(facet_idx)] = candidates[1];
          }
          else if(candidates[1] == element_key)
          {
            element.neighbors[std::size_t(facet_idx)] = candidates[0];
          }
        }
      };
 
      // Find elements adjacent to all facets
      for(Index i(0); i < _layers.size(); i++)
      {
        for(auto entry : _layers[i].template permanent_entities<shape_dim>())
        {
          determine_facet_neighbors(Layer{i}, entry);
        }
        for(auto entry : _layers[i].template transient_entities<shape_dim>())
        {
          determine_facet_neighbors(Layer{i}, entry);
        }
      }
 
      for(Index i(0); i < _layers.size(); i++)
      {
        for(auto entry : _layers[i].template permanent_entities<shape_dim>())
        {
          determine_neighbors(Layer{i}, entry);
        }
        for(auto entry : _layers[i].template transient_entities<shape_dim>())
        {
          determine_neighbors(Layer{i}, entry);
        }
      }
 
      _has_neighbors = true;
    }
 
    void transform(const VertexType& origin, const VertexType& angles, const VertexType& offset)
    {
      // create rotation matrix
      Tiny::Matrix<CoordType, num_coords, num_coords> rot;
      _aux_rot_mat(rot, angles);
 
      // transform all vertices
      for(auto& layer : _layers)
      {
        for(auto entry : layer.vertices())
        {
          entry.value.set_mat_vec_mul(rot, entry.value - origin) += offset;
        }
      }
    }
 
    void reindex()
    {
      // NOTE: The index-schema works as follows. On each layer of the mesh,
      // first the permanent entities get assigned indices, then the transient
      // entities get assigned indices. Indices are assigned in order,
      // starting on each layer with the index one beyond the last index of the
      // previous layers permanent indices. In other words, the transient
      // indices of each layer get reused by the next layer.  This is fine,
      // because transient entities are replaced with their children on the
      // next layer and thus no longer require an index.
      //
      // Visually the indeces are thus staggered like this:
      //
      // L0: |-- Perm. Indices 0 --|---- Trans. Indices 0 ----|
      // L1: |-- Perm. Indices 0 --|-- Perm. Indices 1 --|- Trans. Indices 1 -|
      // L2: |-- Perm. Indices 0 --|-- Perm. Indices 1 --|---- Perm. Indices 2 ----|-- Trans. Indices 2 --|
 
      Index next_vertex_index = 0;
      Index next_edge_index = 0;
      Index next_face_index = 0;
      Index next_cell_index = 0;
 
      for(auto& layer : _layers)
      {
        if constexpr(MeshShape_::dimension >= 0)
        {
          // There are only permanent vertices. We can just index them without
          // having to worry about transient indices.
          for(auto entry : layer.vertices())
          {
            entry.value.index = next_vertex_index++;
          }
        }
 
        if constexpr(MeshShape_::dimension >= 1)
        {
          _index_layer<1>(layer, next_edge_index);
        }
 
        if constexpr(MeshShape_::dimension >= 2)
        {
          _index_layer<2>(layer, next_face_index);
        }
 
        if constexpr(MeshShape_::dimension >= 3)
        {
          _index_layer<3>(layer, next_cell_index);
        }
      }
 
      // Mesh is freshly indexed, by definition
      _is_indexed = true;
    }
 
    // Index-based Accessors
 
    const AdaptiveVertexType& get_vertex_by_index(Layer layer, Index idx) const
    {
      ASSERTM(layer.idx < _layers.size(), "Trying to retrieve vertex of non-existing mesh layer!");
 
      // NOTE: We require the storage to be indexed, despite not making use of
      // the index member of any mesh element. This is so that, given a
      // ElementKey k pointing to a vertex v with index i on layer l, it is
      // impossible that get_index(k) != i or get_vertex_by_index(l, i) != v
      ASSERT(_is_indexed);
 
      // OPTIMIZATION: This is a O(num_layers) search for the right index
      // range.  We can optimize this slightly by precomputing the index range
      // for each layer, eliminating the subtraction of the current idx.  We
      // can then optimize this to O(log num_layers) by building a search tree.
      // And finally we can optimize the average case by caching the
      // index-range of the last access in a thread-local variable, allowing
      // linear iteration in O(1).
 
      // Search through permanent edges of previous and current layer
      for(Index layer_idx = 0; layer_idx <= layer.idx; layer_idx++)
      {
        const MeshLayerType& current_layer = _layers[layer_idx];
 
        if(idx < current_layer.num_vertices())
        {
          return current_layer.vertices().data()[idx];
        }
        idx -= current_layer.num_vertices();
      }
 
      XABORTM("Tried retrieving non-existent vertex index from AdaptiveMeshStorage!");
    }
 
    AdaptiveVertexType& get_vertex_by_index(Layer layer, Index idx)
    {
      // Reuse const accessor by casting this to const and then removing const from the return value
      // SAFE: Casting to const is always fine
      // SAFE: Removing const is fine because this method was called on a non const value to begin with
      return const_cast<AdaptiveVertexType&>(std::as_const(*this).get_vertex_by_index(layer, idx));
    }
 
    template<int dim_>
    const AdaptiveElementByDim<dim_>& get_by_index(Layer layer, Index idx) const
    {
      ASSERTM(layer.idx < _layers.size(), "Trying to retrieve entity of non-existing mesh layer!");
 
      // NOTE: We require the storage to be indexed, despite not making use of
      // the index member of any mesh element. This is so that, given a
      // ElementKey k pointing to a mesh element e with index i on layer l, it is
      // impossible that get_index(k) != i or get_by_index(l, i) != e
      ASSERT(_is_indexed);
 
      // Search through permanent edges of previous and current layer
      for(Index layer_idx = 0; layer_idx <= layer.idx; layer_idx++)
      {
        const MeshLayerType& current_layer = _layers[layer_idx];
 
        if(idx < current_layer.template num_permanent_entities<dim_>())
        {
          return current_layer.template permanent_entities<dim_>().data()[idx];
        }
        idx -= current_layer.template num_permanent_entities<dim_>();
      }
 
      const MeshLayerType& current_layer = _layers[layer.idx];
      if(idx < current_layer.template num_transient_entities<dim_>())
      {
        return current_layer.template transient_entities<dim_>().data()[idx];
      }
 
      XABORTM("Trying to retrieve non-exisiting element from AdaptiveMeshStorage");
    }
 
    template<int dim_>
    AdaptiveElementByDim<dim_>& get_by_index(Layer layer, Index idx)
    {
      // Reuse const accessor by casting this to const and then removing const from the return value
      // SAFE: Casting to const is always fine
      // SAFE: Removing const is fine because this method was called on a non const value to begin with
      return const_cast<AdaptiveElementByDim<dim_>&>(std::as_const(*this).template get_by_index<dim_>(layer, idx));
    }
 
    Index get_neighbor(Layer layer, Index element_idx, Index neighbor_idx) const
    {
      static constexpr int shape_dim = MeshShape_::dimension;
      ASSERT(_has_neighbors);
      ASSERT(_is_indexed);
 
      auto& element = get_by_index<shape_dim>(layer, element_idx);
      if(element.neighbors[neighbor_idx] != neighbor_sentinel)
      {
        return get_index(element.neighbors[neighbor_idx]);
      }
      return ~Index(0);
    }
 
    // Key-based Access Operators
 
    template<int dim_>
    Index get_index(ElementKey<dim_> key) const
    {
      ASSERT(_is_indexed);
      return (*this)[key].index;
    }
 
    template<int dim_>
    Index get_index(OrientedElement<dim_> oriented_key) const
    {
      ASSERT(_is_indexed);
      return (*this)[oriented_key.key].index;
    }
 
    AdaptiveVertexType& operator[](VertexKey key)
    {
      ASSERT(key.layer.idx < _layers.size());
      return _layers[key.layer.idx].get_vertex(key);
    }
 
    const AdaptiveVertexType& operator[](VertexKey key) const
    {
      ASSERT(key.layer.idx < _layers.size());
      return _layers[key.layer.idx].get_vertex(key);
    }
 
    template<int dim_>
    AdaptiveElementByDim<dim_>& operator[](ElementKey<dim_> key)
    {
      ASSERT(key.layer.idx < _layers.size());
      return _layers[key.layer.idx].get(key);
    }
 
    template<int dim_>
    const AdaptiveElementByDim<dim_>& operator[](ElementKey<dim_> key) const
    {
      ASSERT(key.layer.idx < _layers.size());
      return _layers[key.layer.idx].get(key);
    }
 
    AdaptiveVertexType& operator[](OrientedElement<0> ref)
    {
      return (*this)[ref.key];
    }
 
    const AdaptiveVertexType& operator[](OrientedElement<0> ref) const
    {
      return (*this)[ref.key];
    }
 
    template<int dim_>
    AdaptiveElementByDim<dim_>& operator[](OrientedElement<dim_> ref)
    {
      return (*this)[ref.key];
    }
 
    template<int dim_>
    const AdaptiveElementByDim<dim_>& operator[](OrientedElement<dim_> ref) const
    {
      return (*this)[ref.key];
    }
 
    std::size_t bytes() const
    {
      std::size_t total = 0;
      for(auto& layer : _layers)
      {
        total += layer.bytes();
      }
      return total;
    }
 
  private:
    template<int dim_>
    Index _num_entities(Layer layer) const
    {
      if constexpr(dim_ > MeshShape_::dimension)
      {
        return 0;
      }
      else
      {
        ASSERT(layer.idx < _layers.size());
 
        Index result = 0;
        for(Index layer_idx = 0; layer_idx <= layer.idx; layer_idx++)
        {
          result += _layers[layer_idx].template num_permanent_entities<dim_>();
        }
        result += _layers[layer.idx].template num_transient_entities<dim_>();
 
        return result;
      }
    }
 
    template<typename Shape_, int dim_ = Shape_::dimension>
    std::array<Index, 4> _erase_children(
        typename TemplateSet::template RefinementTypeByDim<Shape_::dimension> type,
        ElementChildren<TemplateSet, Shape_>& children)
    {
      if constexpr(dim_ >= 0)
      {
        auto result = _erase_children<Shape_, dim_ - 1>(type, children);
 
        auto& array = children.template by_dim<dim_>();
        for(Index i(0); i < TemplateSet::template num_children<Shape_::dimension, dim_>(type); i++)
        {
          std::array<Index, 4> num_erased = erase(array[i]);
 
          result[0] += num_erased[0];
          result[1] += num_erased[1];
          result[2] += num_erased[2];
          result[3] += num_erased[3];
        }
 
        return result;
      }
      else
      {
        return {};
      }
    }
 
    template<typename Shape_, int dim_ = Shape_::dimension, typename VisitorType>
    void _walk_subtree_children(const AdaptiveElement<TemplateSet, Shape_>& root, VisitorType& visitor, Layer max_depth) const
    {
      if constexpr(dim_ >= 0)
      {
        _walk_subtree_children<Shape_, dim_ - 1>(root, visitor, max_depth);
 
        const auto& children = root.children.template by_dim<dim_>();
        for(Index i = 0; i < TemplateSet::template num_children<Shape_::dimension, dim_>(root.type); i++)
        {
          walk_subtree(children[i], visitor, max_depth);
        }
      }
    }
 
    void _add_mesh_layers(Index new_top_layer)
    {
      while(_layers.size() <= new_top_layer)
      {
        _layers.push_back(MeshLayerType{});
      }
    }
 
    template<int dim_>
    void _index_layer(MeshLayerType& layer, Index& next_permanent_index)
    {
      // First assign indices to all permanent mesh elements
      for(auto entry : layer.template permanent_entities<dim_>())
      {
        entry.value.index = next_permanent_index++;
      }
 
      // Creat a copy of the current next index
      Index next_transient_index = next_permanent_index;
 
      // Assign indices to all transient mesh elements, using the copied index
      // This causes these indices to be reused on the next layer, where the
      // transient elements have been replaced with their children and no longer
      // exist
      for(auto entry : layer.template transient_entities<dim_>())
      {
        entry.value.index = next_transient_index++;
      }
    }
 
    template<int sm_, int sn_, int sv_>
    static void _aux_rot_mat(
      Tiny::Matrix<CoordType, 1, 1, sm_, sn_>& r,
      const Tiny::Vector<CoordType, 1, sv_>& /*unsused*/)
    {
      r.set_identity();
    }
 
    template<int sm_, int sn_, int sv_>
    static void _aux_rot_mat(
      Tiny::Matrix<CoordType, 2, 2, sm_, sn_>& r,
      const Tiny::Vector<CoordType, 2, sv_>& a)
    {
      r.set_rotation_2d(a[0]);
    }
 
    template<int sm_, int sn_, int sv_>
    static void _aux_rot_mat(
      Tiny::Matrix<CoordType, 3, 3, sm_, sn_>& r,
      const Tiny::Vector<CoordType, 3, sv_>& a)
    {
      r.set_rotation_3d(a[0], a[1], a[2]);
    }
  };
} // namespace FEAT::Geometry::Intern