adaptive_mesh.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
 
// includes, FEAT
#include "kernel/geometry/intern/adaptive_refinement_utils.hpp"
#include <kernel/geometry/templates/two_refinement_data.hpp>
#include <kernel/geometry/intern/congruency_mapping.hpp>
#include <kernel/geometry/intern/congruency_sampler.hpp>
#include <kernel/geometry/intern/face_index_mapping.hpp>
#include <kernel/geometry/templates/schneiders_data.hpp>
#include <kernel/geometry/templates/sun_zhao_ma_data.hpp>
#include <kernel/geometry/templates/sun_zhao_ma_expansion_data.hpp>
#include <kernel/shape.hpp>
#include <kernel/base_header.hpp>
#include <kernel/geometry/conformal_mesh.hpp>
#include <kernel/geometry/intern/adaptive_mesh_algorithms.hpp>
#include <kernel/geometry/intern/adaptive_mesh_storage.hpp>
#include <kernel/geometry/intern/refinement_field.hpp>
#include <kernel/geometry/mesh_part.hpp>
#include <kernel/geometry/subdivision_levels.hpp>
#include <kernel/geometry/template_interface.hpp>
#include <kernel/util/assertion.hpp>
#include <kernel/geometry/template_sets.hpp>
#include <kernel/geometry/template_builder.hpp>
 
// includes, system
#include <array>
#include <limits>
#include <set>
#include <unordered_map>
 
namespace FEAT::Geometry
{
  using Layer = Intern::Layer;
 
  // Possible template sets
  using SchneidersTemplates = StandardTemplateSet<SchneidersData>;
  using SunZhaoMaTemplates = IsolatedPointTemplateSet<SunZhaoMaData>;
  using SunZhaoMaExpansionTemplates = IsolatedPointTemplateSet<SunZhaoMaExpansionData>;
  using TwoRefinementTemplates = StandardTemplateSet<TwoRefinementData>;
 
  namespace Intern
  {
    template<int dim_>
    struct MeshRoots : public MeshRoots<dim_ - 1>
    {
      std::unordered_map<Index, Intern::ElementKey<dim_>> roots;
 
      template<int query_dim_>
      std::unordered_map<Index, Intern::ElementKey<query_dim_>>& by_dim()
      {
        static_assert(query_dim_ <= dim_);
        return MeshRoots<query_dim_>::roots;
      }
 
      template<int query_dim_>
      const std::unordered_map<Index, Intern::ElementKey<query_dim_>>& by_dim() const
      {
        static_assert(query_dim_ <= dim_);
        return MeshRoots<query_dim_>::roots;
      }
    };
 
    template<>
    struct MeshRoots<0>
    {
      std::unordered_map<Index, Intern::ElementKey<0>> roots;
 
      template<int query_dim_>
      std::unordered_map<Index, Intern::ElementKey<0>>& by_dim()
      {
        static_assert(query_dim_ == 0);
        return roots;
      }
 
      template<int query_dim_>
      const std::unordered_map<Index, Intern::ElementKey<0>>& by_dim() const
      {
        static_assert(query_dim_ == 0);
        return roots;
      }
    };
  }
 
  struct AdaptionStats
  {
    std::array<Index, 4> added;
 
    std::array<Index, 4> removed;
 
    std::array<Index, 4> kept;
 
    void reset()
    {
      added.fill(0);
      removed.fill(0);
      kept.fill(0);
    }
 
    template<int dim_>
    void added_element(Index count = 1)
    {
      added[dim_] += count;
    }
 
    template<int dim_>
    void removed_element(Index count = 1)
    {
      removed[dim_] += count;
    }
 
    void removed_elements(const std::array<Index, 4>& counts)
    {
      for(Index i = 0; i < 4; i++)
      {
        removed.at(i) += counts.at(i);
      }
    }
 
    template<int dim_>
    void kept_element(Index count = 1)
    {
      kept[dim_] += count;
    }
  };
 
  template<
    typename TemplateSet_,
    typename Shape_,
    int num_coords_ = Shape_::dimension,
    typename Coord_ = Real>
  class AdaptiveMesh
  {
    static_assert(TemplateSet_::template is_shape_compatible<Shape_>());
 
  public:
    using ShapeType = Shape_;
 
    static constexpr int shape_dim = ShapeType::dimension;
 
    static constexpr int world_dim = ShapeType::dimension;
 
    using CoordType = Coord_;
 
    using FoundationMeshType = ConformalMesh<Shape_, num_coords_, Coord_>;
 
    using VertexType = typename FoundationMeshType::VertexType;
 
    using TemplateSet = TemplateSet_;
 
    enum class ImportBehaviour : std::uint8_t
    {
      RequiredOnly, 
      All,          
    };
 
  private:
    using MeshStorage = Intern::AdaptiveMeshStorage<Shape_, TemplateSet, VertexType>;
 
    // The following definitions are helpers for the code below
 
    using VertexKey = typename MeshStorage::template ElementRefByDim<0>;
    using EdgeKey = typename MeshStorage::template ElementRefByDim<1>;
    using FaceKey = typename MeshStorage::template ElementRefByDim<2>;
    using CellKey = typename MeshStorage::template ElementRefByDim<3>;
 
    using OrientedEdge = typename MeshStorage::template OrientedElementRefByDim<1>;
    using OrientedFace = typename MeshStorage::template OrientedElementRefByDim<2>;
 
    using AdaptiveVertex = typename MeshStorage::AdaptiveVertexType;
 
    template<int dim_>
    using LevelTuple =
      Intern::RefinementFieldTuple<
        typename TemplateSet::VertexMarkerType,
        Shape::FaceTraits<typename Shape::FaceTraits<Shape_, dim_>::ShapeType, 0>::count>;
 
    template<int dim_>
    using ElementRef = typename MeshStorage::template ElementRefByDim<dim_>;
 
    template<int dim_>
    using OrientedElementRef = typename MeshStorage::template OrientedElementRefByDim<dim_>;
 
    template<int dim_>
    using ElementTopology = typename MeshStorage::template ElementTopologyByDim<dim_>;
 
    template<int dim_>
    using ElementChildren = typename MeshStorage::template ElementChildrenByDim<dim_>;
 
    template<int dim_>
    using AdaptiveElement = typename MeshStorage::template AdaptiveElementByDim<dim_>;
 
    using MeshRoots = Intern::MeshRoots<shape_dim>;
 
    using VertexMarkerType = typename TemplateSet::VertexMarkerType;
 
    // Algorithms that interact with the foundation mesh are written to recurse
    // over the dimension, similar to other algorithms in kernel/geometry. Some
    // of these algorithms require access to internals of the AdaptiveMesh.
 
    template<typename AdaptiveMeshType_, typename TargetMeshType_, typename AlgShape_>
    friend struct Intern::MeshPartProjector;
 
    template<typename AdaptiveMeshType_, int topology_dim_, int collection_dim_>
    friend struct Intern::FoundationTopologyCollector;
 
    const FoundationMeshType& _foundation_mesh;
 
    Intern::MeshRoots<shape_dim> _roots = {};
 
    MeshStorage _storage = {};
 
    AdaptionStats _stats = {};
 
  public:
    explicit AdaptiveMesh(const FoundationMeshType& mesh) : _foundation_mesh(mesh) {}
 
    AdaptiveMesh(AdaptiveMesh& other) = delete;
    AdaptiveMesh(AdaptiveMesh&& other) = delete;
 
    AdaptiveMesh& operator=(AdaptiveMesh& other) = delete;
    AdaptiveMesh& operator=(AdaptiveMesh&& other) = delete;
 
    AdaptionStats
    adapt(SubdivisionLevels& foundation_levels, ImportBehaviour import_behaviour = ImportBehaviour::RequiredOnly)
    {
      // We expect every vertex to be marked with a subdivision level
      ASSERTM(foundation_levels.size() == _foundation_mesh.get_num_vertices(), "Invalid SDL size");
 
      // Algorithm types
      using EntityCollector = Intern::EntityCollector<ShapeType, TemplateSet, FoundationMeshType>;
 
      // Reset the adaption stats to zero
      _stats.reset();
 
      // Change foundation_levels to match template set requirements.
      Intern::RefinementField<typename TemplateSet::VertexMarkerType> r_field =
        TemplateSet::make_refinement_field(_foundation_mesh, foundation_levels);
 
      // We collect sets of all entities that are required to be added to layer
      // 0 of the adaptive mesh. This includes in all cases all elements with
      // non-zero refinement type and all entities making up those elements. Note
      // that sub-entities are added irrespective of their own refinement type,
      // as they are always required for producing a complete refinement tree.
      Intern::MeshIndexSet set;
 
      // NOTE: Building the sets of all required entities is nice for
      // determining which entities can be safely removed from the mesh, as we
      // can just remove all elements not in the set. Otherwise we would need to
      // determine which entities have both a zero refinement type and are not a
      // sub-entity of an element with non-zero refinement type. We can change
      // this if the hashing becomes a performance bottleneck or the memory-usage
      // causes problems.
 
      // Actually collect those sets of entities
      bool import_all = import_behaviour == ImportBehaviour::All;
      EntityCollector::collect(set, _foundation_mesh, r_field, import_all);
 
      // Create, adapt, or replace refinement trees to match new subdivision levels.
      _adapt_roots(set, r_field);
 
      // Delete all refinement trees of entiites that were previously refined, but are no longer refined.
      _garbage_collect(_roots, set);
 
      // Re-index the storage to allow for index-based access
      _storage.reindex();
 
      // Return a copy of the stats for this adapatation.
      return _stats;
    }
 
    template<typename VertexMarker_>
    static FoundationMeshType create_refined_mesh(const FoundationMeshType& foundation_mesh, const VertexMarker_& marker)
    {
      // Create subdivision levels
      const Index num_vertices = foundation_mesh.get_num_vertices();
      SubdivisionLevels sdls(num_vertices);
 
      for(Index i(0); i < num_vertices; ++i)
      {
        sdls[i] = std::uint64_t(marker(i));
      }
 
      // Create mesh
      AdaptiveMesh adaptive_mesh(foundation_mesh);
      adaptive_mesh.adapt(sdls, ImportBehaviour::All);
      return adaptive_mesh.to_conformal_mesh(Layer{adaptive_mesh.num_layers() - 1});
    }
 
    FoundationMeshType to_conformal_mesh(Layer layer) const
    {
      // Algorithm type
      using MeshWriter = Intern::ConformalMeshWriter<AdaptiveMesh, FoundationMeshType, ShapeType>;
 
      // Determine result mesh size
      std::array<Index, 4> entities = {};
      entities[0] = get_num_entities(layer, 0);
      entities[1] = get_num_entities(layer, 1);
      entities[2] = get_num_entities(layer, 2);
      entities[3] = get_num_entities(layer, 3);
      FoundationMeshType result(entities.data());
 
      // Write mesh
      MeshWriter::write_to_mesh(result, *this, layer);
      result.fill_neighbors();
 
      return result;
    }
 
    template<typename TargetMeshType_>
    MeshPart<TargetMeshType_> project_meshpart(Layer layer, const MeshPart<FoundationMeshType>& part) const
    {
      // Algorithm types
      using Projector = Intern::MeshPartProjector<AdaptiveMesh, TargetMeshType_, ShapeType>;
      using Visitor = Intern::MeshPartProjectorVisitor<MeshStorage>;
 
      // For projecting a mesh part, we need to collect all child keys on the
      // appropriate layers. This means both the keys of type zero entities
      // that are (grand)-children of any entity in the original mesh part, as
      // well as non-zero type elements on the target layer. The visitor will
      // do so in a PreOrder iteration through the refinement trees,
      // maintaining the relative order of entities.
      Visitor visitor(_storage, layer);
 
      // Collect keys of mesh part, performs a PreOrder iteration on all
      // refinement trees rooted at an entity of the original mesh part
      Projector::collect_keys(visitor, part, layer, _roots, _storage);
 
      // Create result mesh part
      std::array<Index, 4> num_elements =
        {visitor.vertices.size(), visitor.edges.size(), visitor.faces.size(), visitor.cells.size()};
      MeshPart<TargetMeshType_> result(num_elements.data(), false);
 
      // Convert collected keys to indices and write into result mesh part
      Projector::write_meshpart(visitor, result, _storage);
 
      return result;
    }
 
    template<int dim_>
    Index num_total_entities() const
    {
      return _storage.template num_total_entities<dim_>();
    }
 
    VertexType& vertex(Layer layer, Index vertex_index)
    {
      return _storage.get_vertex_by_index(layer, vertex_index).vertex;
    }
 
    const VertexType& vertex(Layer layer, Index vertex_index) const
    {
      return _storage.get_vertex_by_index(layer, vertex_index).vertex;
    }
 
    template<int dim_, int codim_>
    Index get_face_index(Layer layer, Index entity_idx, int face_idx) const
    {
      auto& entity = _storage.template get_by_index<dim_>(layer, entity_idx);
 
      return _storage.get_index(entity.topology.template key_by_dim<codim_>(face_idx));
    }
 
    Index get_num_entities(Layer layer, int dim) const
    {
      return _storage.num_entities(layer, dim);
    }
 
    template<int dim_>
    std::optional<Index> get_child(Layer layer, Index parent_idx, Index child) const
    {
      auto& element = _storage.template get_by_index<dim_>(layer, parent_idx);
 
      if(element.type.is_zero_refinement())
      {
        // If the parent has no children, it will be part of all further layers with the same Index.
        // For the purposes of transfer operators and similar, it is thus its own child
        return parent_idx;
      }
 
      if(child < TemplateSet::template num_children<dim_, dim_>(element.type))
      {
        // A child with the wanted index exists. Retrieve the key and convert to index.
        return _storage.get_index(element.children.template by_dim<dim_>()[child]);
      }
 
      // No child with the wanted index exists
      return std::nullopt;
    }
 
    template<int dim>
    Index get_num_children(Layer layer, Index elem)
    {
      auto type = _storage.template get_by_index<dim>(layer, elem).type;
      if(!type.is_zero_refinement())
      {
        return TemplateSet::template num_children<dim, dim>(type);
      }
 
      return 1;
    }
 
    bool has_vertex_changed(Layer layer, Index vertex_idx) const
    {
      auto& vertex = _storage.get_vertex_by_index(layer, vertex_idx);
      return vertex.layer <= layer && layer <= vertex.last_changed;
    }
 
    std::vector<std::pair<Index, Index>> bridge_pairs()
    {
      std::vector<std::pair<Index, Index>> pairs;
 
      for(auto& iter : _roots.vertices)
      {
        pairs.emplace_back(std::make_pair(iter.first, _storage.get_index(iter.second)));
      }
 
      return pairs;
    }
 
    FoundationMeshType& foundation_mesh()
    {
      return _foundation_mesh;
    }
 
    const FoundationMeshType& foundation_mesh() const
    {
      return _foundation_mesh;
    }
 
    template<int dim_>
    std::optional<Index> get_overlap_cell(Index foundation_index) const
    {
      const std::unordered_map<Index, ElementRef<dim_>>& root_map = _roots.template by_dim<dim_>();
 
      auto iter = root_map.find(foundation_index);
      if(iter != root_map.end())
      {
        return _storage.get_index(iter->second);
      }
 
      return std::nullopt;
    }
 
    std::vector<Index> overlap_mapping()
    {
      std::vector<Index> result(_storage.num_entities(0, 0));
 
      for(auto& entry : _roots.vertices)
      {
        result[_storage.get_index(entry.second)] = entry.first;
      }
 
      return result;
    }
 
    Index num_layers() const
    {
      return _storage.num_layers();
    }
 
    MeshPart<FoundationMeshType> exclusive_meshpart() const
    {
      // Get element roots
      const std::unordered_map<Index, ElementRef<shape_dim>>& root_map = _roots.template by_dim<ShapeType::dimension>();
 
      // Determine mesh part size
      const Index total_cells = _foundation_mesh.get_num_entities(shape_dim);
      const Index overlap_cells = root_map.size();
      const Index part_cells = total_cells - overlap_cells;
      std::array<Index, shape_dim + 1> part_size = {};
      part_size[shape_dim] = part_cells;
 
      // Create mesh part without topology
      MeshPart<FoundationMeshType> result(part_size.data(), false);
      auto& target_set = result.template get_target_set<shape_dim>();
 
      // Write foundation mesh indices into mesh part
      Index tset_index = 0;
      for(Index i(0); i < total_cells; ++i)
      {
        if(root_map.find(i) == root_map.end())
        {
          // Cell i is not a root of the adaptive mesh. Add it.
          target_set[tset_index] = i;
          tset_index++;
        }
      }
 
      return result;
    }
 
    MeshPart<FoundationMeshType> overlap_meshpart() const
    {
      // Get element roots
      const std::unordered_map<Index, ElementRef<shape_dim>>& root_map = _roots.template by_dim<ShapeType::dimension>();
 
      // Determine result target set size
      const Index overlap_cells = root_map.size();
      std::array<Index, shape_dim + 1> part_size = {};
      part_size[shape_dim] = overlap_cells;
 
      // Create mesh part without topology
      MeshPart<FoundationMeshType> result(part_size.data(), false);
      auto& target_set = result.template get_target_set<shape_dim>();
 
      // Write foundation mesh indices into mesh part
      Index tset_index = 0;
      for(auto& entry : root_map)
      {
        target_set[tset_index] = entry.first;
        tset_index++;
      }
 
      return result;
    }
 
    SubdivisionLevels transfer_sdls(SubdivisionLevels& lvls_in) const
    {
      const Index num_layers = _storage.num_layers();
      const Layer fine_layer = Layer{num_layers - 1};
      const Index num_base_vertices = _storage.num_entities(Layer{0}, 0);
      const Index num_fine_vertices = _storage.num_entities(fine_layer, 0);
 
      SubdivisionLevels result(num_base_vertices);
 
      for(Index v_fine_idx = 0; v_fine_idx < num_fine_vertices; v_fine_idx++)
      {
        // Find vertex in foundation mesh that is closest to the current fine vertex
        Real closest_distance = std::numeric_limits<Real>::max();
        Index closest_vertex = 0;
        for(Index v_base_idx = 0; v_base_idx < num_base_vertices; v_base_idx++)
        {
          const AdaptiveVertex& v_fine = _storage.get_vertex_by_index(fine_layer, v_fine_idx);
          const AdaptiveVertex& v_base = _storage.get_vertex_by_index(Layer{0}, v_base_idx);
          const Real distance = (v_fine.vertex - v_base.vertex).norm_euclid_sqr();
          if(distance < closest_distance)
          {
            closest_distance = distance;
            closest_vertex = v_base_idx;
          }
        }
 
        // Assign the subdivision levels of the current vertex to the closest vertex in the foundation mesh
        result[closest_vertex] = std::max(result[closest_vertex], lvls_in[v_fine_idx]);
      }
      return result;
    }
 
    void fill_neighbors()
    {
      _storage.fill_neighbors();
    }
 
    Index get_neighbor(Layer layer, Index element_idx, Index neighbor_idx) const
    {
      return _storage.get_neighbor(layer, element_idx, neighbor_idx);
    }
 
    void transform(const VertexType& origin, const VertexType& angles, const VertexType& offset)
    {
      _storage.transform(origin, angles, offset);
    }
 
  private:
    template<int dim_>
    Intern::ElementKey<dim_> _build(Index depth, const LevelTuple<dim_>& levels, const ElementTopology<dim_>& topology)
    {
      ASSERT(_dbg_check_topology_unique<dim_>(topology));
      ASSERT(_dbg_check_topology_orientations<dim_>(topology));
      ASSERT(_dbg_is_topology_consistent<dim_>(topology));
      ASSERT(_dbg_is_topology_layering_consistent<dim_>(Intern::Layer{depth}, topology));
 
      _stats.added_element<dim_>();
 
      auto type = typename TemplateSet::template RefinementTypeByDim<dim_>(levels);
 
      //std::cout << "Building element of type " << type << "\n";
      // Create tree root
      auto self = AdaptiveElement<dim_>(type, Layer{depth}, topology);
 
      if(type.is_zero_refinement())
      {
        // Zero type elements have no children. We are done here.
        auto key = _storage.insert(self);
        //std::cout << "Built element " << key << "\n";
        return key;
      }
 
      using TemplateShape = typename Shape::FaceTraits<ShapeType, dim_>::ShapeType;
 
      static constexpr int num_vertices = Shape::FaceTraits<TemplateShape, 0>::count;
 
      const RefinementTemplate<TemplateShape>& tmplt = TemplateSet::get_template(type);
 
      // Build child vertices
      _stats.added_element<0>(tmplt.template num_entities<0>());
 
      // Retrieve vertex coordinates of topology for following vertex interpolation
      std::array<VertexType, num_vertices> vertex_coordinates;
 
      for(int i(0); i < num_vertices; i++)
      {
        vertex_coordinates[std::size_t(i)] = _storage[topology.template key_by_dim<0>(i)].vertex;
      }
 
      // Produce child vertices by linear interpolation
      auto& child_vertices = self.children.template by_dim<0>();
      for(Index i = 0; i < tmplt.template num_entities<0>(); i++)
      {
        const VertexType vertex(Intern::interpolate(vertex_coordinates, tmplt.get_vertex_coefficients()[i]));
        child_vertices[i] = _storage.insert(AdaptiveVertex(vertex, Layer{depth + 1}));
      }
 
      // We are recording the last mesh layer on which any of the entities surrounding a vertex has changed.
      // This allows disregarding these vertices in a BPX solver later.
      // NOTE: This is quite specific. Can we move this out of the build routine somehow?
      if(tmplt.template num_entities<0>() > 0)
      {
        // Parent vertices created new vertices at this depth.
        // Increase referenced depth for BPX star-sum
        auto& topo_vertices = self.topology.template by_dim<0>();
        for(Index v = 0; v < topo_vertices.size(); v++)
        {
          auto& vertex = _storage[topo_vertices[v]];
          if(depth + 1 > vertex.last_changed.idx)
          {
            vertex.last_changed = Layer{depth + 1};
          }
        }
      }
 
      // Build child edges
      if constexpr(dim_ >= 1)
      {
        auto& children = self.children.template by_dim<1>();
        const auto& topo_templates = tmplt.template get_topologies<1>();
        for(Index i = 0; i < tmplt.template num_entities<1>(); i++)
        {
          ElementTopology<1> child_topo;
          _build_topology<dim_, 1>(topo_templates[i], child_topo, topology, self.children);
 
          children[i] = _build<1>(depth + 1, _levels<dim_, 1>(topo_templates[i], levels), child_topo);
        }
      }
 
      // Build faces
      if constexpr(dim_ >= 2)
      {
        auto& children = self.children.template by_dim<2>();
        const auto& topo_templates = tmplt.template get_topologies<2>();
        for(Index i = 0; i < tmplt.template num_entities<2>(); i++)
        {
          ElementTopology<2> child_topo;
          _build_topology<dim_, 2>(topo_templates[i], child_topo, topology, self.children);
 
          children[i] = _build<2>(depth + 1, _levels<dim_, 2>(topo_templates[i], levels), child_topo);
        }
      }
 
      // Build cells
      if constexpr(dim_ >= 3)
      {
        auto& children = self.children.template by_dim<3>();
        const auto& topo_templates = tmplt.template get_topologies<3>();
        for(Index i = 0; i < tmplt.template num_entities<3>(); i++)
        {
          ElementTopology<3> child_topo;
          _build_topology<dim_, 3>(topo_templates[i], child_topo, topology, self.children);
 
          children[i] = _build<3>(depth + 1, _levels<dim_, 3>(topo_templates[i], levels), child_topo);
        }
      }
 
      return _storage.insert(self);
    }
 
    template<int dim_>
    void _adapt_roots_of_dim(const Intern::MeshIndexSet& set, Intern::RefinementField<typename TemplateSet::VertexMarkerType>& sdls)
    {
      // Algorithm for collecting topologies from the foundation mesh. Entry
      // dimension is dim_ - 1, because topologies only contain entities of
      // smaller dimension than the entity the topology belongs to.
      using TopologyCollector = Intern::FoundationTopologyCollector<AdaptiveMesh, dim_, dim_ - 1>;
 
      // Vertices of foundation entities, for collecting subdivision levels
      auto& foundation_entities = _foundation_mesh.template get_index_set<dim_, 0>();
      auto& roots = _roots.template by_dim<dim_>();
 
      for(Index idx : set.template by_dim<dim_>())
      {
        // Determine refinement type of current entity
        auto markings = sdls.get_tuple(foundation_entities[idx]);
        auto type = typename TemplateSet::template RefinementTypeByDim<dim_>(markings);
 
        // Determine entities entry in mesh roots, if it exists
        auto iter = roots.find(idx);
        bool exists = iter != roots.end();
 
        if(exists && _storage.type(iter->second) == type)
        {
          // Element already exists and its type has not changed.
          // Adapt its children further.
          _adapt_children<dim_>(iter->second, markings);
        }
        else
        {
          //std::cout << "Rebuilding element with new type " << type << "\n";
          //std::cout << "Markings: ";
          for(Index i(0); i < markings.size; i++)
          {
            //std::cout << markings[i] << ", ";
          }
          //std::cout << "\n";
          // Element either does not exist or its type has changed.
          // (Re)build it.
 
          // We need the entity's topology to rebuild it.
          // NOTE: Even if the entity already exists in the mesh and is only
          // going to be replaced, we can not reuse its topology. If we entered
          // this branch the entity's refinement type has changed, which means
          // the refinement type of at least one sub-entity has also changed.
          // That sub-entity has itself been replaced at this point, and the
          // reference stored in this entity's topology is invalid.
 
          // OPTIMIZATION: While it is true that we can't reuse the previous
          // topology, if it exists, we are currently doing some needless work
          // in that case. The topology collection also determines the
          // orientation of any references in the topology. Those orientations
          // don't change, so it would be enough to only update the references,
          // if the topology has been collected before.
          ElementTopology<dim_> topology;
          TopologyCollector::collect(topology, idx, *this);
 
          if(exists)
          {
            // Entity already existed. Replace it.
            _erase(iter->second);
            iter->second = _build<dim_>(0, markings, topology);
          }
          else
          {
            // Entity is entirely new. Build it.
            roots[idx] = _build<dim_>(0, markings, topology);
          }
        }
      }
    }
 
    void _adapt_roots(const Intern::MeshIndexSet& set, Intern::RefinementField<typename TemplateSet::VertexMarkerType>& sdls)
    {
      if constexpr(shape_dim >= 0)
      {
        auto& vertex_set = _foundation_mesh.get_vertex_set();
        auto& vertex_roots = _roots.template by_dim<0>();
        for(const auto& vert : set.vertices)
        {
          if(vertex_roots.find(vert) == vertex_roots.end())
          {
            _stats.template added_element<0>();
            vertex_roots[vert] = _storage.insert(vertex_set[vert], Layer{0});
          }
          else
          {
            _stats.template kept_element<0>();
          }
        }
      }
      if constexpr(shape_dim >= 1)
      {
        _adapt_roots_of_dim<1>(set, sdls);
      }
      if constexpr(shape_dim >= 2)
      {
        _adapt_roots_of_dim<2>(set, sdls);
      }
      if constexpr(shape_dim >= 3)
      {
        _adapt_roots_of_dim<3>(set, sdls);
      }
    }
 
    template<int dim_>
    void _adapt_children(ElementRef<dim_> key, const LevelTuple<dim_>& levels)
    {
      // Calculate element type
      const auto type = _storage.type(key);
      const Index depth = _storage[key].layer.idx;
 
      //std::cout << "Adapting element of type " << type << "\n";
      // Adapting an entity's children requires a valid topology for the parent
      // element.  If this entity's type changed, then it should have been
      // rebuilt with a valid topology and calling _adapt_children is an error.
      XASSERTM(type == typename TemplateSet::template RefinementTypeByDim<dim_>(levels), "Called adapt_children on element with changed type.");
 
      // We are keeping the current entity as part of the mesh
      _stats.kept_element<dim_>();
 
      if(type.is_zero_refinement())
      {
        // Zero-type elements have no children that need further adaptation.
        return;
      }
 
      // Retrieve reference to current entity
      auto& self = _storage[key];
 
      using TemplateShape = typename Shape::FaceTraits<ShapeType, dim_>::ShapeType;
      const RefinementTemplate<TemplateShape>& tmplt = TemplateSet::get_template(type);
 
      // Retrieve vertex query from template set. Vertices do not need
      // adaptation, but we guarantee template sets that we will construct all
      // template queries in order.
      _stats.kept_element<0>(tmplt.template num_entities<0>());
 
      // Adapt child edges
      if constexpr(dim_ >= 1)
      {
        _adapt_sub_entities<dim_, 1>(depth, self, tmplt, levels);
 
        // Adapt child faces
        if constexpr(dim_ >= 2)
        {
          _adapt_sub_entities<dim_, 2>(depth, self, tmplt, levels);
 
          // Adapt child cells
          if constexpr(dim_ >= 3)
          {
            _adapt_sub_entities<dim_, 3>(depth, self, tmplt, levels);
          }
        }
      }
    }
 
  protected:
    template<int dim_>
    void _erase(Intern::ElementKey<dim_> key)
    {
      _stats.removed_elements(_storage.erase(key));
    }
 
    template<int parent_dim_, int child_dim_>
    void _adapt_sub_entities(
        Index depth,
        AdaptiveElement<parent_dim_>& self,
        const RefinementTemplate<typename Shape::FaceTraits<ShapeType, parent_dim_>::ShapeType>& tmplt,
        const LevelTuple<parent_dim_>& levels)
    {
      // children is an std::array of slotmap keys
      auto& children = self.children.template by_dim<child_dim_>();
      auto& topo_templates = tmplt.template get_topologies<child_dim_>();
 
      for(Index i = 0; i < tmplt.template num_entities<child_dim_>(); i++)
      {
        // Determine child SubdivisionLevels and type
        LevelTuple<child_dim_> child_levels = _levels<parent_dim_, child_dim_>(topo_templates[i], levels);
        auto child_type = typename TemplateSet::template RefinementTypeByDim<child_dim_>(child_levels);
 
        if(_storage.type(children[i]) == child_type)
        {
          // Type is unchanged. Adapt children recursively.
          _adapt_children<child_dim_>(children[i], child_levels);
        }
        else
        {
          // Type is changed. Erase child and rebuild it.
          _erase(children[i]);
 
          ElementTopology<child_dim_> child_topology;
          _build_topology<parent_dim_, child_dim_>(topo_templates[i], child_topology, self.topology, self.children);
 
          children[i] = _build<child_dim_>(depth + 1, child_levels, child_topology);
        }
      }
    }
 
    template<int parent_dim_, int entity_dim_>
    LevelTuple<entity_dim_> _levels(
        const TopologyTemplate<typename Shape::FaceTraits<ShapeType, entity_dim_>::ShapeType>& tmplt,
        const LevelTuple<parent_dim_>& lvls)
    {
      // OPTIMIZATION(mmuegge): Currenlty a new minimum is calculated for all
      // vertices. We can use that all sibling vertices share the same new
      // marking. As do all vertices on any of the boundary elements. Add a
      // MarkingSpreader class or something like that and cache the results.
      using ParentShape = typename Shape::FaceTraits<ShapeType, parent_dim_>::ShapeType;
 
      const auto& vertex_refs = tmplt.template get_references<0>();
 
      LevelTuple<entity_dim_> result;
 
      for(int i(0); i < result.size; ++i)
      {
        result[i] = TemplateSet::template spread_refinement_field<ParentShape>(vertex_refs[std::size_t(i)], lvls);
      }
 
      return result;
    }
 
    template<int parent_dim_, int topo_dim_, int dim_ = topo_dim_ - 1>
    void _build_topology(
        const TopologyTemplate<typename Shape::FaceTraits<ShapeType, topo_dim_>::ShapeType>& tmplt,
        ElementTopology<topo_dim_>& topo,
        const ElementTopology<parent_dim_>& parent_topo,
        const ElementChildren<parent_dim_>& siblings)
    {
      if constexpr(dim_ >= 0)
      {
        _build_topology<parent_dim_, topo_dim_, dim_ - 1>(tmplt, topo, parent_topo, siblings);
 
        using TopologyShape = typename Shape::FaceTraits<ShapeType, topo_dim_>::ShapeType;
        static constexpr int num_entities = Shape::FaceTraits<TopologyShape, dim_>::count;
 
        auto& target = topo.template by_dim<dim_>();
        for(int i(0); i < num_entities; i++)
        {
          EntityReference ref = tmplt.template get_reference<dim_>(i);
          target[std::size_t(i)] = _resolve_entity_reference<parent_dim_, dim_>(parent_topo, siblings, ref);
 
          if constexpr(dim_ >= 1)
          {
            if(ref.source == EntitySource::BoundaryEdge || ref.source == EntitySource::BoundaryFace)
            {
              // Calc orientation
              using ResultShape = typename Shape::FaceTraits<ShapeType, dim_>::ShapeType;
              using FaceMapping = Intern::FaceIndexMapping<TopologyShape, dim_, 0>;
              using CongruencySampler = Intern::CongruencySampler<ResultShape>;
 
              static constexpr int num_vertices = Shape::FaceTraits<ResultShape, 0>::count;
              std::array<ElementRef<0>, num_vertices> local_keys;
              std::array<ElementRef<0>, num_vertices> foreign_keys;
 
              const auto& entity = _storage[target[std::size_t(i)]];
              for(int vertex(0); vertex < num_vertices; vertex++)
              {
                int local_vertex = FaceMapping::map(i, vertex);
                local_keys[std::size_t(vertex)] = topo.template key_by_dim<0>(local_vertex).key;
                foreign_keys[std::size_t(vertex)] = entity.topology.template key_by_dim<0>(vertex).key;
              }
 
              target[std::size_t(i)].orientation = CongruencySampler::compare(local_keys.data(), foreign_keys.data());
            }
          }
        }
      }
    }
 
    /*
     * \brief Map a EntityReference to an OrientedElementRef
     *
     * \tparam entity_dim_ Dimension of the current entity
     * \tparam result_dim_ Dimension the reference refers to
     *
     * \param[in] topo Current topology
     * \param[in] siblings Current siblings
     * \param[in] ref EntityReference
     */
    template<int entity_dim_, int result_dim_>
    OrientedElementRef<result_dim_> _resolve_entity_reference(
        const ElementTopology<entity_dim_>& topo,
        const ElementChildren<entity_dim_>& siblings,
        const EntityReference& ref)
    {
      switch(ref.source)
      {
        case EntitySource::ParentTopology:
          return topo.template key_by_dim<result_dim_>(ref.index);
        case EntitySource::Sibling:
          return OrientedElementRef<result_dim_>(ref.orientation, siblings.template by_dim<result_dim_>()[std::size_t(ref.index)]);
        case EntitySource::BoundaryEdge:
        {
          if constexpr (entity_dim_ >= 2 && result_dim_ <= 1)
          {
            return _resolve_boundary_entity_reference<entity_dim_, result_dim_, 1>(topo, ref);
          }
          else
          {
            XABORTM("_resolve_entity_reference failed: EntitySource::BoundaryEdge: Either source or target has invalid dimension");
          }
        }
        case EntitySource::BoundaryFace:
        {
          if constexpr (entity_dim_ >= 3 && result_dim_ <= 2)
          {
            return _resolve_boundary_entity_reference<entity_dim_, result_dim_, 2>(topo, ref);
          }
          else
          {
            XABORTM("_resolve_entity_reference failed: EntitySource::BoundaryFace: Either source or target has invalid dimension");
          }
        }
        default:
          XABORTM("_resolve_entitiy_reference failed: Unkown EntitySource");
      }
    }
 
    template<int entity_dim_, int result_dim_, int boundary_dim_>
    OrientedElementRef<result_dim_> _resolve_boundary_entity_reference(
        const ElementTopology<entity_dim_>& topo,
        const EntityReference& ref)
    {
      ASSERT(ref.source == EntitySource::BoundaryEdge || ref.source == EntitySource::BoundaryFace);
 
      // Retrieve type of boundary element
      using BoundaryShape = typename Shape::FaceTraits<ShapeType, boundary_dim_>::ShapeType;
      OrientedElementRef<boundary_dim_> boundary_ref = topo.template key_by_dim<boundary_dim_>(ref.entity);
      typename TemplateSet::template RefinementTypeByDim<boundary_dim_> type = _storage.type(boundary_ref.key);
 
      // Use type and orientation to correct the index of the reference for orientation
      std::pair<Index, int> orientation_mapping =
        TemplateSet::template correct_for_orientation<BoundaryShape, result_dim_>(type, boundary_ref.orientation, ref.index);
 
      // Retrieve correct child based on orientation mapping
      AdaptiveElement<boundary_dim_>& boundary_element = _storage[boundary_ref];
      ElementRef<result_dim_> result_ref =
        boundary_element.children.template by_dim<result_dim_>()[orientation_mapping.first];
 
      // Return oriented reference
      return OrientedElementRef<result_dim_>(orientation_mapping.second, result_ref);
    }
 
    template<int dim_>
    void _delete_unused(std::unordered_map<Index, ElementRef<dim_>>& roots, const std::set<Index>& elements)
    {
      auto iter = roots.begin();
      while(iter != roots.end())
      {
        if(elements.find(iter->first) == elements.end())
        {
          _erase(iter->second);
          iter = roots.erase(iter);
        }
        else
        {
          iter++;
        }
      }
    }
 
    template<int dim_ = shape_dim>
    void _garbage_collect(Intern::MeshRoots<shape_dim>& roots, const Intern::MeshIndexSet& set)
    {
      if constexpr(dim_ >= 0)
      {
        _garbage_collect<dim_ -1>(roots, set);
        _delete_unused<dim_>(roots.template by_dim<dim_>(), set.by_dim<dim_>());
      }
    }
 
    template<int topo_dim_, int dim_ = topo_dim_ - 1>
    bool _dbg_check_topology_unique(const ElementTopology<topo_dim_>& topology)
    {
      if constexpr(dim_ >= 0)
      {
        using TopologyShape = typename Shape::FaceTraits<ShapeType, topo_dim_>::ShapeType;
 
        bool all_unique = true;
 
        static constexpr int num_entities = Shape::FaceTraits<TopologyShape, dim_>::count;
 
        const auto& entities = topology.template by_dim<dim_>();
        for(int i(0); i < num_entities; i++)
        {
          for(int j(0); j < num_entities; j++)
          {
            all_unique = all_unique && (i == j || (entities[std::size_t(i)] != entities[std::size_t(j)]));
          }
        }
 
        return all_unique && _dbg_check_topology_unique<topo_dim_, dim_ - 1>(topology);
      }
      else
      {
        return true;
      }
    }
 
    template<int topo_dim_, int dim_ = topo_dim_ - 1>
    bool _dbg_check_topology_orientations(const ElementTopology<topo_dim_>& topology)
    {
      if constexpr(dim_ == 0)
      {
        return true;
      }
      else
      {
        bool is_oriented = true;
 
        using TopologyShape = typename Shape::FaceTraits<ShapeType, topo_dim_>::ShapeType;
        using EntityShape = typename Shape::FaceTraits<ShapeType, dim_>::ShapeType;
 
        using FaceMapping = Intern::FaceIndexMapping<TopologyShape, dim_, 0>;
        using CongruencySampler = Intern::CongruencySampler<EntityShape>;
 
        static constexpr int num_entities = Shape::FaceTraits<TopologyShape, dim_>::count;
        static constexpr int num_vertices = Shape::FaceTraits<EntityShape, 0>::count;
 
        const auto& entities = topology.template by_dim<dim_>();
 
        for(int entity_idx(0); entity_idx < num_entities; entity_idx++)
        {
          const Intern::OrientedElement<dim_> o_ref = entities[std::size_t(entity_idx)];
          const auto& entity = _storage[o_ref];
 
          std::array<ElementRef<0>, num_vertices> local_keys;
          std::array<ElementRef<0>, num_vertices> foreign_keys;
 
          for(int vertex(0); vertex < num_vertices; vertex++)
          {
            const int local_vertex = int(FaceMapping::map(entity_idx, vertex));
            local_keys[std::size_t(vertex)] = topology.template key_by_dim<0>(local_vertex).key;
            foreign_keys[std::size_t(vertex)] = entity.topology.template key_by_dim<0>(vertex).key;
          }
 
          is_oriented = is_oriented && (CongruencySampler::compare(local_keys.data(), foreign_keys.data()) == o_ref.orientation);
        }
 
        return is_oriented && _dbg_check_topology_orientations<topo_dim_, dim_ - 1>(topology);
      }
    }
 
    template<int topo_dim_, int dim_ = topo_dim_ - 1>
    bool _dbg_is_topology_consistent(const ElementTopology<topo_dim_>& topology)
    {
      if constexpr(dim_ == 0)
      {
        return true;
      }
      else
      {
        bool is_consistent = true;
 
        using TopologyShape = typename Shape::FaceTraits<ShapeType, topo_dim_>::ShapeType;
        using EntityShape = typename Shape::FaceTraits<ShapeType, dim_>::ShapeType;
 
        using FaceMapping = Intern::FaceIndexMapping<TopologyShape, dim_, 0>;
        using CongruencyMapping = Intern::CongruencyMapping<EntityShape, 0>;
 
        static constexpr int num_entities = Shape::FaceTraits<TopologyShape, dim_>::count;
        static constexpr int num_vertices = Shape::FaceTraits<EntityShape, 0>::count;
 
        const auto& entities = topology.template by_dim<dim_>();
        for(int entity_idx(0); entity_idx < num_entities; entity_idx++)
        {
          const auto& entity = _storage[entities[std::size_t(entity_idx)]];
 
          for(int vertex(0); vertex < num_vertices; vertex++)
          {
            ElementRef<0> v_topo = topology.template by_dim<0>()[std::size_t(FaceMapping::map(entity_idx, vertex))].key;
            const int orientation = entities[std::size_t(entity_idx)].orientation;
            const Index mapped = Index(CongruencyMapping::map(orientation, vertex));
            ElementRef<0> v_entity = entity.topology.template by_dim<0>()[mapped].key;
            is_consistent = is_consistent && (v_topo == v_entity);
 
            if(v_topo != v_entity)
            {
              std::cout << "Inconsistent topology. Vertex " << vertex << " of entity " << entity_idx << " of dimension << " << dim_ << " mismatches!\n";
            }
          }
        }
 
        return is_consistent && _dbg_is_topology_consistent<topo_dim_, dim_ - 1>(topology);
      }
    }
 
    template<int topo_dim_, int dim_ = topo_dim_ - 1>
    bool _dbg_is_topology_layering_consistent(Intern::Layer layer, const ElementTopology<topo_dim_>& topology)
    {
      if constexpr(dim_ == 0)
      {
        bool is_consistent = true;
 
        // Vertices are layered consistently as long as we pick no vertices from a higher layer
        // NOTE(mmuegge): This can be merged with the dim_ != 0 case, if we
        // mark vertices as permanent by default. Which we should do anyway,
        // because they are
        for(const auto& entity : topology.template by_dim<0>())
        {
          is_consistent = is_consistent && entity.key.layer <= layer;
        }
 
        return is_consistent;
      }
      else
      {
        bool is_consistent = _dbg_is_topology_layering_consistent<topo_dim_, dim_ - 1>(layer, topology);
 
        for(const auto& entity : topology.template by_dim<dim_>())
        {
          is_consistent = is_consistent && (entity.key.layer == layer || (entity.key.layer < layer && entity.key.is_permanent));
        }
 
        return is_consistent;
      }
    }
  }; // class AdaptiveMesh
} // namespace FEAT::Geometry