template_builder.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/base_header.hpp>
#include <kernel/geometry/intern/congruency_mapping.hpp>
#include <kernel/geometry/intern/congruency_sampler.hpp>
#include <kernel/geometry/intern/congruency_trafo.hpp>
#include <kernel/geometry/intern/face_index_mapping.hpp>
#include <kernel/geometry/raw_refinement_templates.hpp>
#include <kernel/geometry/refinement_types.hpp>
#include <kernel/shape.hpp>
#include <kernel/util/assertion.hpp>
#include <kernel/util/string.hpp>
#include <kernel/util/tiny_algebra.hpp>
 
#include <algorithm>
#include <array>
#include <cstdint>
#include <unordered_map>
#include <vector>
#include <iostream>
 
namespace FEAT::Geometry
{
  namespace Intern
  {
    template<typename Shape_>
    constexpr int orientations()
    {
      if constexpr(std::is_same_v<Shape_, Shape::Hypercube<1>>)
      {
        return 2;
      }
 
      if constexpr(std::is_same_v<Shape_, Shape::Hypercube<2>>)
      {
        return 8;
      }
 
      return 0;
    }
 
    template<int n_>
    bool is_same_vertex(const Tiny::Vector<Real, n_>& a, const Tiny::Vector<Real, n_>& b)
    {
      bool is_same = true;
 
      const Real tol = Math::pow(Math::eps<Real>(), Real(0.7));
      for(int i(0); i < n_; ++i)
      {
        if(Math::abs(a[i] - b[i]) > tol)
        {
          is_same = false;
        }
      }
 
      return is_same;
    }
 
    template<typename Shape_, int num_coords_ = Shape_::dimension>
    bool is_same_raw_entity(const RawEntity<Shape_, num_coords_>& a, const RawEntity<Shape_, num_coords_>& b)
    {
      return std::is_permutation(a.coords.begin(), a.coords.end(), b.coords.begin(), Intern::is_same_vertex<num_coords_>);
    }
 
    template<typename Shape_, int num_coords_ = Shape_::dimension>
    std::vector<Tiny::Vector<Real, num_coords_>> all_vertices(const std::vector<RawEntity<Shape_, num_coords_>>& entities)
    {
      std::vector<Tiny::Vector<Real, num_coords_>> vertices;
 
      for(const auto& entity : entities)
      {
        for(const auto& vertex : entity.coords)
        {
          vertices.push_back(vertex);
        }
      }
 
      return vertices;
    }
 
    template<int n_>
    bool is_internal_vertex(Tiny::Vector<Real, n_> vertex)
    {
      const Real tol = Math::pow(Math::eps<Real>(), Real(0.7));
 
      // A vertex is internal if all its coordinates are in (0, 1)
      bool is_internal = true;
      for(int i(0); i < n_; ++i)
      {
        if(vertex[i] < tol || (1 - tol) < vertex[i])
        {
          is_internal = false;
        }
      }
      return is_internal;
    }
 
    template<typename Shape_, int num_coords_ = Shape_::dimension>
    bool is_internal_entity(const RawEntity<Shape_, num_coords_>& entity)
    {
      Tiny::Vector<Real, num_coords_> centroid(0);
      for(const auto& vertex : entity.coords)
      {
        centroid += vertex;
      }
      return is_internal_vertex((1.0 / entity.coords.size()) * centroid);
    }
 
    template<typename T_, typename Comp>
    std::vector<T_> deduplicate(Comp comparison, const std::vector<T_>& list)
    {
      std::vector<T_> deduplicated;
      for(const T_& current : list)
      {
        auto pred = [&](T_& t) { return comparison(current, t); };
        bool is_new = !std::any_of(deduplicated.begin(), deduplicated.end(), pred);
        if(is_new)
        {
          deduplicated.push_back(current);
        }
      }
      return deduplicated;
    }
 
    template<typename T_, typename Predicate_>
    std::vector<T_> filter(Predicate_ predicate, const std::vector<T_>& list)
    {
      std::vector<T_> filtered;
      for(const T_& current : list)
      {
        if(predicate(current))
        {
          filtered.push_back(current);
        }
      }
      return filtered;
    }
 
    template<typename Shape_>
    std::vector<Tiny::Vector<Real, Shape_::dimension>> internal_vertices(const RawTemplate<Shape_>& tmplt)
    {
      std::vector<Tiny::Vector<Real, Shape_::dimension>> vertices;
 
      for(const auto& entity : tmplt.entities)
      {
        for(const auto& vertex : entity.coords)
        {
          vertices.push_back(vertex);
        }
      }
 
      vertices = filter([](auto vertex) { return is_internal_vertex(vertex); }, vertices);
      vertices = deduplicate([](auto a, auto b) { return Intern::is_same_vertex(a, b); }, vertices);
 
      return vertices;
    }
 
    template<typename Shape_, int face_dim_, int num_coords_ = Shape_::dimension>
    std::vector<RawEntity<typename Shape::FaceTraits<Shape_, face_dim_>::ShapeType, num_coords_>>
    internal_entities(const std::vector<RawEntity<Shape_, num_coords_>>& entities)
    {
      if constexpr(Shape_::dimension == face_dim_)
      {
        return entities;
      }
      else
      {
        static constexpr const int faces_per_entity = Shape::FaceTraits<Shape_, face_dim_>::count;
 
        std::vector<RawEntity<typename Shape::FaceTraits<Shape_, face_dim_>::ShapeType, num_coords_>> result;
 
        for(const auto& entity : entities)
        {
          for(int i(0); i < faces_per_entity; ++i)
          {
            result.push_back(entity.template face<face_dim_>(i));
          }
        }
 
        result = filter([](auto& entity) { return is_internal_entity(entity); }, result);
        result = deduplicate([](auto& a, auto& b) { return Intern::is_same_raw_entity(a, b); }, result);
 
        return result;
      }
    }
 
    template<typename Shape_>
    static Tiny::Vector<Real, Shape_::dimension>
    orient_vertex(const Tiny::Vector<Real, Shape_::dimension>& vertex, int orientation)
    {
      static constexpr int num_coords = Shape_::dimension;
      using Trafo = Intern::CongruencyTrafo<Shape_>;
 
      Tiny::Vector<Real, num_coords> offset(0.5);
 
      Tiny::Matrix<Real, num_coords, num_coords> transform(0.0);
      Tiny::Vector<Real, num_coords> unused(0.0);
      Trafo::compute(transform, unused, orientation);
 
      Tiny::Vector<Real, num_coords> result;
      result.set_mat_vec_mult(transform, (vertex - offset));
 
      return result + offset;
    }
 
    template<typename Shape_, int num_coords_ = Shape_::dimension>
    static RawEntity<Shape_, num_coords_> orient_raw_entity(const RawEntity<Shape_, num_coords_>& entity, int orientation)
    {
      RawEntity<Shape_, num_coords_> result;
      for(int i(0); i < entity.coords.size(); i++)
      {
        result.coords[i] = Intern::orient_vertex<Shape_>(entity.coords[i], orientation);
      }
    }
 
 
    template<typename Shape_, int source_dim_>
    struct VertexEmbedder
    {
    };
 
    template<>
    struct VertexEmbedder<Shape::Quadrilateral, 1>
    {
      using SourceVertex = Tiny::Vector<Real, 1>;
      using TargetVertex = Tiny::Vector<Real, 2>;
 
      static TargetVertex embed(int elem, const SourceVertex& vertex)
      {
        switch(elem)
        {
        case 0:
          return Tiny::Vector<Real, 2>{vertex[0], 0.0};
        case 1:
          return Tiny::Vector<Real, 2>{vertex[0], 1.0};
        case 2:
          return Tiny::Vector<Real, 2>{0.0, vertex[0]};
        case 3:
          return Tiny::Vector<Real, 2>{1.0, vertex[0]};
        default:
          XABORTM("Invalid face-edge index.");
        }
      }
    };
 
    template<>
    struct VertexEmbedder<Shape::Hexahedron, 1>
    {
      using SourceVertex = Tiny::Vector<Real, 1>;
      using TargetVertex = Tiny::Vector<Real, 3>;
 
      static TargetVertex embed(int elem, const SourceVertex& vertex)
      {
        switch(elem)
        {
        case 0:
          return Tiny::Vector<Real, 3>{vertex[0], 0.0, 0.0};
        case 1:
          return Tiny::Vector<Real, 3>{vertex[0], 1.0, 0.0};
        case 2:
          return Tiny::Vector<Real, 3>{vertex[0], 0.0, 1.0};
        case 3:
          return Tiny::Vector<Real, 3>{vertex[0], 1.0, 1.0};
        case 4:
          return Tiny::Vector<Real, 3>{0.0, vertex[0], 0.0};
        case 5:
          return Tiny::Vector<Real, 3>{1.0, vertex[0], 0.0};
        case 6:
          return Tiny::Vector<Real, 3>{0.0, vertex[0], 1.0};
        case 7:
          return Tiny::Vector<Real, 3>{1.0, vertex[0], 1.0};
        case 8:
          return Tiny::Vector<Real, 3>{0.0, 0.0, vertex[0]};
        case 9:
          return Tiny::Vector<Real, 3>{1.0, 0.0, vertex[0]};
        case 10:
          return Tiny::Vector<Real, 3>{0.0, 1.0, vertex[0]};
        case 11:
          return Tiny::Vector<Real, 3>{1.0, 1.0, vertex[0]};
        default:
          XABORTM("Invalid cell-edge index.");
        }
      }
    };
 
    template<>
    struct VertexEmbedder<Shape::Hexahedron, 2>
    {
      using SourceVertex = Tiny::Vector<Real, 2>;
      using TargetVertex = Tiny::Vector<Real, 3>;
 
      static TargetVertex embed(int elem, const SourceVertex& vertex)
      {
        switch(elem)
        {
        case 0:
          return Tiny::Vector<Real, 3>{vertex[0], vertex[1], 0.0};
        case 1:
          return Tiny::Vector<Real, 3>{vertex[0], vertex[1], 1.0};
        case 2:
          return Tiny::Vector<Real, 3>{vertex[0], 0.0, vertex[1]};
        case 3:
          return Tiny::Vector<Real, 3>{vertex[0], 1.0, vertex[1]};
        case 4:
          return Tiny::Vector<Real, 3>{0.0, vertex[0], vertex[1]};
        case 5:
          return Tiny::Vector<Real, 3>{1.0, vertex[0], vertex[1]};
        default:
          XABORTM("Invalid cell-edge index.");
        }
      }
    };
 
    template<typename Shape_>
    std::array<Tiny::Vector<Real, Shape_::dimension>, Shape::FaceTraits<Shape_, 0>::count> topology_vertices()
    {
      static constexpr int num_coords = Shape_::dimension;
      static constexpr int num_vertices = Shape::FaceTraits<Shape_, 0>::count;
 
      std::array<Tiny::Vector<Real, num_coords>, num_vertices> vertices;
 
      for(int i(0); i < num_vertices; i++)
      {
        // Build topology vertex from reference cell
        Tiny::Vector<Real, num_coords> vertex;
        for(int j(0); j < num_coords; ++j)
        {
          vertex[j] = 0.5 * (Shape::ReferenceCell<Shape_>::template vertex<Real>(i, j) + 1.0);
        }
        vertices[i] = vertex;
      }
 
      return vertices;
    }
 
  } // namespace Intern
 
  // Refinement template data types
 
  enum class EntitySource : std::uint8_t
  {
    ParentTopology, 
    Sibling,        
    BoundaryEdge,   
    BoundaryFace,   
  };
 
  inline std::ostream& operator<<(std::ostream& stream, EntitySource source)
  {
    static std::unordered_map<EntitySource, String> mapping{
      {EntitySource::ParentTopology, "EntitySource::ParentTopology"},
      {EntitySource::Sibling, "EntitySource::Sibling"},
      {EntitySource::BoundaryEdge, "EntitySource::BoundaryEdge"},
      {EntitySource::BoundaryFace, "EntitySource::BoundaryFace"},
    };
    stream << mapping[source];
    return stream;
  }
 
  struct EntityReference
  {
    EntitySource source;
    Index index;
    int orientation;
    int entity;
 
    EntityReference() = default;
    EntityReference(EntitySource s, Index i, int o, int e) : source(s), index(i), orientation(o), entity(e)
    {
    }
 
    // Equals operator for EntityReferences
    bool operator==(const EntityReference& rhs) const
    {
      return (source == rhs.source) && (index == rhs.index) && (orientation == rhs.orientation) &&
            (entity == rhs.entity);
    }
 
    // Not-equals operator for EntityReferences
    bool operator!=(const EntityReference& rhs) const
    {
      return !(*this == rhs);
    }
  };
 
  inline std::ostream& operator<<(std::ostream& stream, const EntityReference& reference)
  {
    stream << "EntityReference {"
           << "source: " << reference.source << ", index: " << stringify(reference.index)
           << ", orientation: " << stringify(reference.orientation) << ", entity: " << stringify(reference.entity)
           << "}";
    return stream;
  }
 
  template<typename Shape_, int reference_dim_ = Shape_::dimension - 1>
  struct TopologyTemplate : TopologyTemplate<Shape_, reference_dim_ - 1>
  {
    static constexpr int num_entities = Shape::FaceTraits<Shape_, reference_dim_>::count;
 
    std::array<EntityReference, num_entities> references;
 
    template<int dim_>
    std::array<EntityReference, Shape::FaceTraits<Shape_, dim_>::count>& get_references()
    {
      return TopologyTemplate<Shape_, dim_>::references;
    }
 
    template<int dim_>
    const std::array<EntityReference, Shape::FaceTraits<Shape_, dim_>::count>& get_references() const
    {
      return TopologyTemplate<Shape_, dim_>::references;
    }
 
    template<int dim_>
    EntityReference get_reference(int idx) const
    {
      return TopologyTemplate<Shape_, dim_>::references[idx];
    }
  };
 
  template<typename Shape_>
  struct TopologyTemplate<Shape_, 0>
  {
    static constexpr int num_entities = Shape::FaceTraits<Shape_, 0>::count;
 
    std::array<EntityReference, num_entities> references;
 
    template<int dim_>
    std::array<EntityReference, num_entities>& get_references()
    {
      static_assert(dim_ == 0, "invalid reference dimension");
      return references;
    }
 
    template<int dim_>
    const std::array<EntityReference, num_entities>& get_references() const
    {
      static_assert(dim_ == 0, "invalid reference dimension");
      return references;
    }
 
    template<int dim_>
    EntityReference get_reference(Index idx) const
    {
      return TopologyTemplate<Shape_, dim_>::references[idx];
    }
  };
 
  template<typename Shape_, typename CoShape_ = Shape_>
  struct RefinementTemplate
    : RefinementTemplate<Shape_, typename Shape::FaceTraits<CoShape_, CoShape_::dimension - 1>::ShapeType>
  {
    template<int dim_>
    using TopologyTemplateByDim = TopologyTemplate<typename Shape::FaceTraits<Shape_, dim_>::ShapeType>;
 
    std::vector<TopologyTemplate<CoShape_>> topologies;
 
    template<int dim_>
    std::vector<TopologyTemplateByDim<dim_>>& get_topologies()
    {
      using QueryShape = typename Shape::FaceTraits<Shape_, dim_>::ShapeType;
      return RefinementTemplate<Shape_, QueryShape>::topologies;
    }
 
    template<int dim_>
    const std::vector<TopologyTemplateByDim<dim_>>& get_topologies() const
    {
      using QueryShape = typename Shape::FaceTraits<Shape_, dim_>::ShapeType;
      return RefinementTemplate<Shape_, QueryShape>::topologies;
    }
 
    template<int dim_>
    Index num_entities() const
    {
      if constexpr(dim_ > Shape_::dimension)
      {
        return 0;
      }
      else if constexpr(dim_ > 0 && dim_ <= Shape_::dimension)
      {
        using QueryShape = typename Shape::FaceTraits<Shape_, dim_>::ShapeType;
        return RefinementTemplate<Shape_, QueryShape>::topologies.size();
      }
      else
      {
        using QueryShape = typename Shape::FaceTraits<Shape_, dim_>::ShapeType;
        return RefinementTemplate<Shape_, QueryShape>::vertices.size();
      }
    }
  };
 
  template<typename Shape_>
  struct RefinementTemplate<Shape_, Shape::Vertex>
  {
    static constexpr int num_coefficients = Shape::FaceTraits<Shape_, 0>::count;
 
    std::vector<Tiny::Vector<Real, Shape_::dimension>> vertices;
    std::vector<std::array<Real, num_coefficients>> vertex_coefficients;
 
    std::vector<Tiny::Vector<Real, Shape_::dimension>>& get_vertices()
    {
      return vertices;
    }
 
    const std::vector<Tiny::Vector<Real, Shape_::dimension>>& get_vertices() const
    {
      return vertices;
    }
 
    std::vector<std::array<Real, num_coefficients>>& get_vertex_coefficients()
    {
      return vertex_coefficients;
    }
 
    const std::vector<std::array<Real, num_coefficients>>& get_vertex_coefficients() const
    {
      return vertex_coefficients;
    }
 
    template<int dim_>
    Index num_entities() const
    {
      static_assert(dim_ == 0);
 
      return vertices.size();
    }
  };
 
  class OrientationMapping
  {
    std::vector<std::pair<Index, int>> _mapping;
 
  public:
    std::pair<Index, int> operator()(Index idx) const
    {
      return _mapping.at(idx);
    }
 
    std::vector<std::pair<Index, int>>& get_mapping()
    {
      return _mapping;
    }
 
    const std::vector<std::pair<Index, int>>& get_mapping() const
    {
      return _mapping;
    }
  };
 
  template<typename RawData_, typename Shape_, int dim_ = Shape_::dimension>
  struct OrientationMappingWrapper : public OrientationMappingWrapper<RawData_, Shape_, dim_ - 1>
  {
    static constexpr int num_orientations = Intern::orientations<Shape_>();
 
    template<typename S>
    using RefinementTypeByShape = typename RawData_::template RefinementTypeByShape<S>;
 
    std::unordered_map<RefinementTypeByShape<Shape_>, std::array<OrientationMapping, num_orientations>> mappings = {};
 
    template<int n_>
    OrientationMapping& get_mapping(RefinementTypeByShape<Shape_> type, int orientation)
    {
      static_assert(n_ <= dim_);
 
      ASSERT(orientation < num_orientations);
 
      auto& mappings_n = OrientationMappingWrapper<RawData_, Shape_, n_>::mappings;
      if(mappings_n.find(type) == mappings_n.end())
      {
        mappings_n[type] = {};
      }
 
      return mappings_n.at(type)[orientation];
    }
 
    template<int n_>
    const OrientationMapping& get_mapping(RefinementTypeByShape<Shape_> type, int orientation) const
    {
      static_assert(n_ <= dim_);
 
      ASSERT(mappings.find(type) != mappings.end());
      ASSERT(orientation < num_orientations);
 
      auto& mappings_n = OrientationMappingWrapper<RawData_, Shape_, n_>::mappings;
      return mappings_n.at(type)[orientation];
    }
  };
 
  template<typename RawData_, typename Shape_>
  struct OrientationMappingWrapper<RawData_, Shape_, 0>
  {
    static constexpr int num_orientations = Intern::orientations<Shape_>();
 
    template<typename S>
    using RefinementTypeByShape = typename RawData_::template RefinementTypeByShape<S>;
 
    std::unordered_map<RefinementTypeByShape<Shape_>, std::array<OrientationMapping, num_orientations>> mappings = {};
 
    template<int n_>
    OrientationMapping& get_mapping(RefinementTypeByShape<Shape_> type, int orientation)
    {
      static_assert(n_ == 0);
 
      ASSERT(mappings.find(type) != mappings.end());
      ASSERT(orientation < num_orientations);
 
      return OrientationMappingWrapper<RawData_, Shape_, n_>::mappings.at(type)[orientation];
    }
 
    template<int n_>
    const OrientationMapping& get_mapping(RefinementTypeByShape<Shape_> type, int orientation) const
    {
      static_assert(n_ == 0);
 
      ASSERT(mappings.find(type) != mappings.end());
      ASSERT(orientation < num_orientations);
 
      return OrientationMappingWrapper<RawData_, Shape_, n_>::mappings.at(type)[orientation];
    }
  };
 
  template<typename RawData_, typename Shape_>
  struct OrientationMappingHolder
    : OrientationMappingHolder<RawData_, typename Shape::FaceTraits<Shape_, Shape_::dimension - 1>::ShapeType>
  {
    OrientationMappingWrapper<RawData_, Shape_> wrapper = {};
 
    template<int dim_>
    using RefinementTypeByDim = typename RawData_::template RefinementTypeByDim<dim_>;
 
    template<int dim_, int codim_>
    OrientationMapping&
    get_mapping(RefinementTypeByDim<dim_> type, int orientation)
    {
      using DimShape = typename Shape::FaceTraits<Shape_, dim_>::ShapeType;
      return OrientationMappingHolder<RawData_, DimShape>::wrapper.template get_mapping<codim_>(
        type,
        orientation);
    }
 
    template<int dim_, int codim_>
    const OrientationMapping&
    get_mapping(RefinementTypeByDim<dim_> type, int orientation) const
    {
      using DimShape = typename Shape::FaceTraits<Shape_, dim_>::ShapeType;
      return OrientationMappingHolder<RawData_, DimShape>::wrapper.template get_mapping<codim_>(
        type,
        orientation);
    }
  };
 
  template<typename RawData_>
  struct OrientationMappingHolder<RawData_, Shape::Vertex>
  {
  };
 
  template<typename Shape_, int num_coords_>
  struct EntitySearchEntry
  {
    RawEntity<Shape_, num_coords_> raw_entity;
 
    EntityReference reference;
 
    EntitySearchEntry(RawEntity<Shape_, num_coords_> ent, EntityReference ref) : raw_entity(ent), reference(ref)
    {
    }
  };
 
  template<typename Shape_, int num_coords_>
  std::ostream& operator<<(std::ostream& stream, const EntitySearchEntry<Shape_, num_coords_>& entry)
  {
    stream << "EntitySearchEntry<" << Shape_::name() << ", " << num_coords_ << "> { raw_entity: " << entry.raw_entity
           << ", reference: " << entry.reference << " }";
    return stream;
  }
 
  template<typename Shape_, int num_coords_ = Shape_::dimension>
  class TemplateSearchSpace
    : public TemplateSearchSpace<typename Shape::FaceTraits<Shape_, Shape_::dimension - 1>::ShapeType, num_coords_>
  {
  public:
    using ShapeType = Shape_;
 
    using Entry = EntitySearchEntry<Shape_, num_coords_>;
    static const constexpr int num_coords = num_coords_;
 
    using BaseClass =
      TemplateSearchSpace<typename Shape::FaceTraits<Shape_, Shape_::dimension - 1>::ShapeType, num_coords_>;
 
    template<int dim_>
    using EntryByDim = EntitySearchEntry<typename Shape::FaceTraits<Shape_, dim_>::ShapeType, num_coords_>;
 
    template<int dim_>
    using RawTemplateByDimes = RawTemplate<typename Shape::FaceTraits<ShapeType, dim_>::ShapeType>;
 
    template<int dim_>
    std::vector<EntryByDim<dim_>>& entries()
    {
      static_assert(dim_ <= Shape_::dimension);
      static_assert(dim_ >= 0);
 
      using DimShape = typename Shape::FaceTraits<Shape_, dim_>::ShapeType;
      return TemplateSearchSpace<DimShape, num_coords_>::_entries;
    }
 
    template<int dim_>
    const std::vector<EntryByDim<dim_>>& entries() const
    {
      static_assert(dim_ <= Shape_::dimension);
      static_assert(dim_ >= 0);
 
      using DimShape = typename Shape::FaceTraits<Shape_, dim_>::ShapeType;
      return TemplateSearchSpace<DimShape, num_coords_>::_entries;
    }
 
    template<typename Shape__>
    EntityReference search(const RawEntity<Shape__, num_coords_>& entity) const
    {
      return TemplateSearchSpace<Shape__, num_coords>::_search(entity);
    }
 
    template<typename Shape__>
    void add_sibling(const RawEntity<Shape__, num_coords_>& entity, Index idx)
    {
      return TemplateSearchSpace<Shape__, num_coords>::_add_sibling(entity, idx);
    }
 
    template<typename EntityShape_, int dim_ = EntityShape_::dimension>
    void add_boundary_entity(int parent_index, const RawTemplate<EntityShape_>& tmplt)
    {
      static_assert(EntityShape_::dimension <= ShapeType::dimension);
 
      using DimShape = typename Shape::FaceTraits<EntityShape_, dim_>::ShapeType;
      using Embedder = Intern::VertexEmbedder<ShapeType, EntityShape_::dimension>;
 
      static constexpr int num_vertices = Shape::FaceTraits<DimShape, 0>::count;
 
      EntitySource source = EntityShape_::dimension == 1 ? EntitySource::BoundaryEdge : EntitySource::BoundaryFace;
 
      if constexpr(dim_ == 0)
      {
        auto& vertices = this->template entries<0>();
        Index vertex_index = 0;
        for(const auto& vertex : Intern::internal_vertices(tmplt))
        {
          RawEntity<Shape::Vertex, num_coords> entity(Embedder::embed(parent_index, vertex));
          vertices.emplace_back(entity, EntityReference{source, vertex_index++, 0, parent_index});
        }
      }
      else
      {
        add_boundary_entity<EntityShape_, dim_ - 1>(parent_index, tmplt);
 
        auto& entries = this->template entries<dim_>();
 
        // Add internal elements of dimension dim_ to search space
        Index idx = 0;
        for(const auto& entity : Intern::internal_entities<EntityShape_, dim_>(tmplt.entities))
        {
          RawEntity<DimShape, num_coords> embedded;
 
          for(int i(0); i < num_vertices; i++)
          {
            embedded.coords[i] = Embedder::embed(parent_index, entity.coords[i]);
          }
          entries.emplace_back(embedded, EntityReference{source, idx++, 0, parent_index});
        }
      }
    }
 
    template<int dim_ = ShapeType::dimension - 1>
    void add_topology()
    {
      static constexpr int num_vertices = Shape::FaceTraits<ShapeType, 0>::count;
 
      std::array<Tiny::Vector<Real, num_coords>, num_vertices> vertices = Intern::topology_vertices<ShapeType>();
      _add_topology_rec(vertices);
    }
 
 
  protected:
    template<int dim_ = ShapeType::dimension - 1>
    void _add_topology_rec(
      const std::array<Tiny::Vector<Real, num_coords>, Shape::FaceTraits<ShapeType, 0>::count>& vertices)
    {
      if constexpr(dim_ == 0)
      {
        static constexpr int num_vertices = Shape::FaceTraits<ShapeType, 0>::count;
 
        auto& entries = this->template entries<0>();
 
        for(int vertex(0); vertex < num_vertices; vertex++)
        {
          RawEntity<Shape::Vertex, num_coords> entity;
          entity.coords[0] = vertices[vertex];
          entries.emplace_back(entity, EntityReference{EntitySource::ParentTopology, static_cast<Index>(vertex), 0, 0});
        }
      }
      else
      {
        using FaceShape = typename Shape::FaceTraits<ShapeType, dim_>::ShapeType;
        using FaceMapping = Intern::FaceIndexMapping<ShapeType, dim_, 0>;
 
        static constexpr int num_faces = Shape::FaceTraits<ShapeType, dim_>::count;
        static constexpr int num_vertices = Shape::FaceTraits<FaceShape, 0>::count;
 
        auto& entries = this->template entries<dim_>();
        for(int face(0); face < num_faces; face++)
        {
          RawEntity<FaceShape, num_coords> entity;
 
          for(int vertex(0); vertex < num_vertices; vertex++)
          {
            entity.coords[vertex] = vertices[FaceMapping::map(face, vertex)];
          }
          entries.emplace_back(entity, EntityReference{EntitySource::ParentTopology, static_cast<Index>(face), 0, 0});
        }
 
        _add_topology_rec<dim_ - 1>(vertices);
      }
    }
 
    EntityReference _search(const RawEntity<Shape_, num_coords_>& entity) const
    {
      auto comp = [&](auto entry) { return Intern::is_same_raw_entity(entry.raw_entity, entity); };
      auto it = std::find_if(_entries.begin(), _entries.end(), comp);
 
      if(it == _entries.end())
      {
        std::cout << "Searching for entity " << entity << "\n";
        XABORTM("Entity not in search space");
      }
 
      // We found a valid reference. Next we need to determine the orientation of entity,
      // as compared to the found reference.
      // Because the CongurencySampler can not directly compare vertices, we produce index-lists here.
      using CongruencySampler = Intern::CongruencySampler<Shape_>;
 
      int orientation = CongruencySampler::compare_with(entity.coords.data(), it->raw_entity.coords.data(), Intern::is_same_vertex<num_coords_>);
      EntityReference ref = it->reference;
 
      ref.orientation = orientation;
 
      if constexpr(num_coords_ < 3)
      {
        // 3D or less templates should not produce boundary face references
        ASSERT(ref.source != EntitySource::BoundaryFace);
      }
 
      return ref;
    }
 
    void _add_sibling(const RawEntity<Shape_, num_coords_>& entity, Index idx)
    {
      _entries.push_back({entity, EntityReference{EntitySource::Sibling, idx, 0, 0}});
    }
 
    std::vector<Entry> _entries;
  };
 
  template<int num_coords_>
  class TemplateSearchSpace<Shape::Vertex, num_coords_>
  {
  public:
    using ShapeType = Shape::Vertex;
 
    using Entry = EntitySearchEntry<Shape::Vertex, num_coords_>;
 
    static const constexpr int num_coords = num_coords_;
 
    template<int dim_>
    std::vector<Entry>& entries()
    {
      static_assert(dim_ == 0);
      return _entries;
    }
 
    template<int dim_>
    const std::vector<Entry>& entries() const
    {
      static_assert(dim_ == 0);
 
      return _entries;
    }
 
    EntityReference search_vertex(const Tiny::Vector<Real, num_coords>& vertex) const
    {
      for(const auto& entry : _entries)
      {
        if(Intern::is_same_vertex(vertex, entry.raw_entity.coords[0]))
        {
          return entry.reference;
        }
      }
      std::cout << "Vertex: " << vertex << " not found\n";
      XABORTM("Vertex not in search space");
    }
 
    void add_sibling_vertex(const Tiny::Vector<Real, num_coords>& vertex, Index idx)
    {
      _entries.push_back(
        {RawEntity<Shape::Vertex, num_coords>(vertex), EntityReference{EntitySource::Sibling, idx, 0, 0}});
    }
 
    Tiny::Vector<Real, num_coords_> vertex_by_reference(const EntityReference& reference) const
    {
      for(const auto& entry : _entries)
      {
        if(entry.reference == reference)
        {
          return entry.raw_entity.coords[0];
        }
      }
      XABORTM("EntifyReference not in search space!");
    }
 
    std::vector<Entry> _entries;
  };
 
 
  template<class RawData_>
  class TemplateBuilder
  {
    using MaxShape = typename RawData_::MaxShape;
 
    template<int dim_>
    using RawTemplateMapByDim = typename RawData_::template TemplateMapByDim<dim_>;
 
    using RawEdgeMap = typename RawData_::EdgeMap;
    using RawFaceMap = typename RawData_::FaceMap;
    using RawCellMap = typename RawData_::CellMap;
 
    template<int dim_>
    using RefinementTypeByDim = typename RawData_::template RefinementTypeByDim<dim_>;
 
    template<int dim_>
    using RefinementTemplateByDim = RefinementTemplate<typename Shape::FaceTraits<MaxShape, dim_>::ShapeType>;
 
    template<int dim_>
    using TemplateMapByDim = std::unordered_map<RefinementTypeByDim<dim_>, RefinementTemplate<Shape::Hypercube<dim_>>>;
 
    using EdgeMap = TemplateMapByDim<1>;
    using FaceMap = TemplateMapByDim<2>;
    using CellMap = TemplateMapByDim<3>;
 
    EdgeMap _edge_templates = {};
    FaceMap _face_templates = {};
    CellMap _cell_templates = {};
 
    OrientationMappingHolder<RawData_, Shape::Hexahedron> _orientation_mappings = {};
 
    std::array<std::array<std::uint64_t, 2>, 4> _edge_type_adjustments = {};
    std::array<std::array<std::uint64_t, 4>, 16> _face_type_adjustments = {};
    std::array<std::array<std::uint64_t, 8>, 256> _cell_type_adjustments = {};
 
    std::array<bool, 4> _type_adjustment_initialized = {false, false, false, false};
 
    const EdgeMap& edges()
    {
      if(_edge_templates.empty())
      {
        _edge_templates = _build<1>(RawData_::raw_edges());
      }
      return _edge_templates;
    }
 
    const FaceMap& faces()
    {
      if(_face_templates.empty())
      {
        _face_templates = _build<2>(RawData_::raw_faces());
      }
      return _face_templates;
    }
 
    const CellMap& cells()
    {
      if(_cell_templates.empty())
      {
        _cell_templates = _build<3>(RawData_::raw_cells());
      }
      return _cell_templates;
    }
 
  public:
    template<int dim_>
    const RefinementTemplateByDim<dim_>& get_template(RefinementTypeByDim<dim_> type)
    {
      if constexpr(dim_ == 1)
      {
        return edges().at(type);
      }
      else if constexpr(dim_ == 2)
      {
        return faces().at(type);
      }
      else if constexpr(dim_ == 3)
      {
        return cells().at(type);
      }
      XABORTM("Unsupported template dimension!");
    }
 
    template<int dim_>
    bool has_template(RefinementTypeByDim<dim_> type)
    {
      if constexpr(dim_ == 1)
      {
        return edges().find(type) != edges().end();
      }
      else if constexpr(dim_ == 2)
      {
        return faces().find(type) != faces().end();
      }
      else if constexpr(dim_ == 3)
      {
        return cells().find(type) != cells().end();
      }
      XABORTM("Unsupported template dimension!");
    }
 
    template<int dim_, int codim_>
    std::pair<Index, int> correct_for_orientation(RefinementTypeByDim<dim_> type, int orientation, Index idx)
    {
      return _orientation_mappings.template get_mapping<dim_, codim_>(type, orientation)(idx);
    }
 
    /*
     * \brief Accessor for type adjustments.
     *
     * \tparam dim_ Dimension to retrieve adjustments for
     *
     * \param[in] type Current markings of a mesh element
     *
     * \returns An array containing a delta of subdivision levels.
     *
     * A template set does not necessarily contain templates for all possible
     * markings of a mesh element. In that case the AdaptiveMesh corrects the
     * markings during the creation of the final refinement field. This method
     * supplies the necessary corrections.
     *
     * \note The adjustments will only ever add subdivision levels.
     */
    template<int dim_>
    const std::array<std::uint64_t, Shape::FaceTraits<Shape::Hypercube<dim_>, 0>::count>&
    type_adjustment(Index type)
    {
      if constexpr(dim_ == 1)
      {
        if(!_type_adjustment_initialized[1])
        {
          _create_type_adjustments<Shape::Hypercube<dim_>>(edges());
        }
        return _edge_type_adjustments.at(type);
      }
      else if constexpr(dim_ == 2)
      {
        if(!_type_adjustment_initialized[2])
        {
          _create_type_adjustments<Shape::Hypercube<dim_>>(faces());
        }
        return _face_type_adjustments.at(type);
      }
      else if constexpr(dim_ == 3)
      {
        if(!_type_adjustment_initialized[3])
        {
          _create_type_adjustments<Shape::Hypercube<dim_>>(cells());
        }
        return _cell_type_adjustments.at(type);
      }
      XABORTM("Unsupported template dimension!");
    }
 
    void stats()
    {
      std::array<Index, 4> maxima = {0, 0, 0, 0};
 
      for(const auto& entry : edges())
      {
        maxima[0] = Math::max(maxima[0], entry.second.template num_entities<0>());
        maxima[1] = Math::max(maxima[1], entry.second.template num_entities<1>());
      }
      std::cout << "Edges:\n";
      std::cout << "  Max Vertices: " << stringify(unsigned(maxima[0])) << "\n";
      std::cout << "  Max Edges: " << stringify(unsigned(maxima[1])) << "\n";
 
      maxima = {0, 0, 0, 0};
      for(const auto& entry : faces())
      {
        maxima[0] = Math::max(maxima[0], entry.second.template num_entities<0>());
        maxima[1] = Math::max(maxima[1], entry.second.template num_entities<1>());
        maxima[2] = Math::max(maxima[2], entry.second.template num_entities<2>());
      }
      std::cout << "Faces:\n";
      std::cout << "  Max Vertices: " << stringify(unsigned(maxima[0])) << "\n";
      std::cout << "  Max Edges: " << stringify(unsigned(maxima[1])) << "\n";
      std::cout << "  Max Faces: " << stringify(unsigned(maxima[2])) << "\n";
 
      maxima = {0, 0, 0, 0};
      for(const auto& entry : cells())
      {
        maxima[0] = Math::max(maxima[0], entry.second.template num_entities<0>());
        maxima[1] = Math::max(maxima[1], entry.second.template num_entities<1>());
        maxima[2] = Math::max(maxima[2], entry.second.template num_entities<2>());
        maxima[3] = Math::max(maxima[3], entry.second.template num_entities<3>());
      }
      std::cout << "Faces:\n";
      std::cout << "  Max Vertices: " << stringify(unsigned(maxima[0])) << "\n";
      std::cout << "  Max Edges: " << stringify(unsigned(maxima[1])) << "\n";
      std::cout << "  Max Faces: " << stringify(unsigned(maxima[2])) << "\n";
      std::cout << "  Max Cells: " << stringify(unsigned(maxima[3])) << "\n";
    }
 
  private:
    template<int n_>
    Index _find_vertex(const std::vector<Tiny::Vector<Real, n_>>& vertices, const Tiny::Vector<Real, n_>& vertex)
    {
      auto pred = [&](auto v) { return Intern::is_same_vertex(vertex, v); };
      return std::distance(vertices.begin(), std::find_if(vertices.begin(), vertices.end(), pred));
    }
 
    template<int dim_>
    auto& type_adjustments()
    {
      if constexpr(dim_ == 1)
      {
        return _edge_type_adjustments;
      }
      else if constexpr(dim_ == 2)
      {
        return _face_type_adjustments;
      }
      else if constexpr(dim_ == 3)
      {
        return _cell_type_adjustments;
      }
      XABORTM("Unsupported template dimension!");
    }
 
    template<typename TemplateShape_, typename TopologyShape_, int dim_ = TopologyShape_::dimension - 1>
    void _build_topology(
      const RawEntity<TopologyShape_, TemplateShape_::dimension>& raw_entity,
      const TemplateSearchSpace<TemplateShape_>& search_space,
      TopologyTemplate<TopologyShape_>& topo)
    {
      static constexpr int num_entities = Shape::FaceTraits<TopologyShape_, dim_>::count;
      if constexpr(dim_ == 0)
      {
        auto& vertex_refs = topo.template get_references<0>();
        for(int i(0); i < num_entities; i++)
        {
          vertex_refs[i] = search_space.search_vertex(raw_entity.coords[i]);
        }
      }
      else
      {
        auto& refs = topo.template get_references<dim_>();
        for(int i(0); i < num_entities; i++)
        {
          refs[i] = search_space.search(raw_entity.template face<dim_>(i));
        }
 
        _build_topology<TemplateShape_, TopologyShape_, dim_ - 1>(raw_entity, search_space, topo);
      }
    }
 
    template<typename TemplateShape_, int dim_ = TemplateShape_::dimension>
    void _build_template(
      const RawTemplate<TemplateShape_>& raw_template,
      RefinementTemplate<TemplateShape_>& tmplt,
      TemplateSearchSpace<TemplateShape_>& search_space)
    {
      if constexpr(dim_ == 0)
      {
        for(const auto& vertex : Intern::internal_vertices(raw_template))
        {
          // Add vertex and its coefficients to template
          tmplt.get_vertices().push_back(vertex);
          tmplt.get_vertex_coefficients().push_back(Intern::vertex_coefficients<TemplateShape_>(vertex));
 
          // Add entry to search space for building topologies later
          search_space.add_sibling_vertex(vertex, tmplt.get_vertices().size() - 1);
        }
      }
      else
      {
        _build_template<TemplateShape_, dim_ - 1>(raw_template, tmplt, search_space);
 
        using TopologyShape = typename Shape::FaceTraits<TemplateShape_, dim_>::ShapeType;
 
        auto& topologies = tmplt.template get_topologies<dim_>();
        for(const auto& entity : Intern::internal_entities<TemplateShape_, dim_>(raw_template.entities))
        {
          TopologyTemplate<TopologyShape> topology{};
          _build_topology<TemplateShape_, TopologyShape>(entity, search_space, topology);
          topologies.push_back(topology);
 
          search_space.add_sibling(entity, topologies.size() - 1);
        }
      }
    }
 
    template<typename TemplateShape_, int dim_ = TemplateShape_::dimension>
    void _fill_orientation_mappings(
      const RefinementTypeByDim<TemplateShape_::dimension>& oriented_type,
      const int orientation,
      const RefinementTemplate<TemplateShape_>& local_template,
      const TemplateSearchSpace<TemplateShape_>& local_search_space,
      const RefinementTemplate<TemplateShape_>& oriented_template,
      const TemplateSearchSpace<TemplateShape_>& oriented_search_space)
    {
 
      using VertexType = Tiny::Vector<Real, TemplateShape_::dimension>;
 
      if constexpr(dim_ == 0)
      {
        OrientationMapping& vertex_mapping =
          _orientation_mappings.template get_mapping<TemplateShape_::dimension, 0>(oriented_type, orientation);
 
        // For each vertex of the local template
        // * map it into the oriented frame of reference
        // * search for the corresponding vertex of the oriented template
        // * save the found index in the orientation mapping
        for(const VertexType& local_vertex : local_template.get_vertices())
        {
          VertexType oriented_vertex = Intern::orient_vertex<TemplateShape_>(local_vertex, orientation);
          Index mapped = _find_vertex(oriented_template.get_vertices(), oriented_vertex);
          vertex_mapping.get_mapping().emplace_back(mapped, 0);
        }
      }
      else
      {
        _fill_orientation_mappings<TemplateShape_, dim_ - 1>(
          oriented_type,
          orientation,
          local_template,
          local_search_space,
          oriented_template,
          oriented_search_space);
 
        using CurrentShape = typename Shape::FaceTraits<TemplateShape_, dim_>::ShapeType;
        static constexpr int num_vertices = Shape::FaceTraits<CurrentShape, 0>::count;
 
        OrientationMapping& mapping =
          _orientation_mappings.template get_mapping<TemplateShape_::dimension, dim_>(oriented_type, orientation);
 
        for(const auto& local_topo : local_template.template get_topologies<dim_>())
        {
          // Collect local vertices from the search space
          std::array<VertexType, num_vertices> local = {};
 
          const auto& local_refs = local_topo.template get_references<0>();
          for(int i(0); i < num_vertices; i++)
          {
            local[i] =
              Intern::orient_vertex<TemplateShape_>(local_search_space.vertex_by_reference(local_refs[i]), orientation);
          }
 
          Index entity_idx = 0;
          for(const auto& oriented_topo : oriented_template.template get_topologies<dim_>())
          {
            // Collect oriented vertices from the oriented search space
            std::array<VertexType, num_vertices> oriented = {};
 
            const auto& oriented_refs = oriented_topo.template get_references<0>();
            for(int i(0); i < num_vertices; i++)
            {
              oriented[i] = oriented_search_space.vertex_by_reference(oriented_refs[i]);
            }
 
            // Compare local and oriented entities
            RawEntity<CurrentShape, TemplateShape_::dimension> ent_local;
            ent_local.coords = local;
 
            RawEntity<CurrentShape, TemplateShape_::dimension> ent_oriented;
            ent_oriented.coords = oriented;
 
            if(Intern::is_same_raw_entity(ent_local, ent_oriented))
            {
              // The entities are identical. This means local_topo and
              // oriented_topo are the same mesh element and we can save this
              // mapping in the orientation mapping. Additionally we can
              // pre-compute their relative orientation at this point, because
              // we know how both of these entities will be constructed in the
              // final mesh and how their parents are oriented relative to each
              // other.
              using CongruencySampler = Intern::CongruencySampler<CurrentShape>;
              int o = CongruencySampler::compare_with(
                ent_local.coords.data(),
                ent_oriented.coords.data(),
                Intern::is_same_vertex<TemplateShape_::dimension>);
 
              mapping.get_mapping().emplace_back(entity_idx, o);
              break;
            }
            entity_idx++;
          }
        }
      }
    }
 
    template<int dim_>
    TemplateMapByDim<dim_> _build(const RawTemplateMapByDim<dim_>& raw_templates)
    {
      using TemplateShape = typename Shape::FaceTraits<Shape::Hypercube<dim_>, dim_>::ShapeType;
 
      TemplateMapByDim<dim_> result;
 
      // We require the full template search spaces to create orientation mappings with.
      std::unordered_map<RefinementTypeByDim<dim_>, TemplateSearchSpace<TemplateShape>> search_spaces;
 
      for(const auto& entry : raw_templates)
      {
        RefinementTemplate<Shape::Hypercube<dim_>> tmplt;
        TemplateSearchSpace<TemplateShape> search_space = _setup_search_space<Shape::Hypercube<dim_>>(entry.first);
 
        _build_template(entry.second, tmplt, search_space);
 
        result.insert({entry.first, tmplt});
        search_spaces.insert({entry.first, search_space});
      }
 
      // We have now created all templates for this dimension.
      // Next, we create the corresponding orientation mappings.
 
      // NOTE(mmuegge): _fill_orientation_mappings is fully dimension-generic,
      // apart from CongruencySamplers for higher dimensions. We thus disable
      // orientation mapping creation for all dimensions above 2
      if constexpr(dim_ <= 2)
      {
        for(const auto& entry : raw_templates)
        {
          RefinementTypeByDim<dim_> local_type = entry.first;
          auto& search_space = search_spaces[local_type];
 
          for(int orientation = 0; orientation < Intern::orientations<Shape::Hypercube<dim_>>(); orientation++)
          {
            const RefinementTypeByDim<dim_> oriented_type = local_type.orient(orientation);
            auto& oriented_search_space = search_spaces[oriented_type];
 
            const auto& local_tmplt = result.at(local_type);
            const auto& oriented_tmplt = result.at(oriented_type);
 
            _fill_orientation_mappings<TemplateShape>(
              oriented_type,
              orientation,
              local_tmplt,
              search_space,
              oriented_tmplt,
              oriented_search_space);
          }
        }
      }
      return result;
    }
 
    template<typename Shape_>
    TemplateSearchSpace<Shape_> _setup_search_space(RefinementTypeByDim<Shape_::dimension> type) const
    {
      // Setup result.
      TemplateSearchSpace<Shape_> result = {};
 
      // For faces and cells we also have to add the boundary edges and internal edges
      if constexpr(Shape_::dimension > 1)
      {
        static const constexpr int num_edges = Shape::FaceTraits<Shape_, 1>::count;
        for(int i(0); i < num_edges; ++i)
        {
          auto edge_type = type.template face_type<1>(i);
          result.add_boundary_entity(i, RawData_::raw_edges().at(edge_type));
        }
      }
 
      // For cells, we have to add boundary faces and internal faces
      if constexpr(Shape_::dimension > 2)
      {
        static const constexpr int num_faces = Shape::FaceTraits<Shape_, 2>::count;
        for(int i(0); i < num_faces; ++i)
        {
          auto face_type = type.template face_type<2>(i);
          result.add_boundary_entity(i, RawData_::raw_faces().at(face_type));
        }
      }
 
      result.add_topology();
 
      return result;
    }
 
    template<typename Shape_, typename TemplateMap_>
    void _create_type_adjustments(const TemplateMap_& available_templates)
    {
      auto& adjustments = type_adjustments<Shape_::dimension>();
 
      static const constexpr int num_verts = Shape::FaceTraits<Shape_, 0>::count;
      for(Index type(0); type < (1 << num_verts); ++type)
      {
        auto type_marking = VertexMarking<Shape_>(type);
 
        std::size_t closest_count = num_verts;
        auto closest_marking = VertexMarking<Shape_>::all_marked();
 
        for(auto& entry : available_templates)
        {
          VertexMarking<Shape_> entry_marking(entry.first.to_vertex_marking());
 
          // Is the current type a subtype of this type?
          if(!entry_marking.covers(type_marking))
          {
            continue;
          }
          std::size_t distance = type_marking.distance(entry_marking);
 
          if(distance < closest_count)
          {
            closest_count = distance;
            closest_marking = entry_marking;
          }
        }
 
        // Found the type that requires the least changes
        // Next, determine exact required changes
 
        auto& adjustment = adjustments.at(type);
        for(int i(0); i < num_verts; ++i)
        {
          if(closest_marking.is_vertex_marked(i) && !type_marking.is_vertex_marked(i))
          {
            adjustment[i] = 1;
          }
        }
      }
    }
  };
} // namespace FEAT::Geometry