adaptive_mesh_layer.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/geometry/common_factories.hpp"
#include "kernel/geometry/mesh_permutation.hpp"
#include "kernel/shape.hpp"
#include "kernel/util/string.hpp"
#include <iterator>
#include <kernel/adjacency/permutation.hpp>
#include <kernel/base_header.hpp>
#include <kernel/geometry/adaptive_mesh.hpp>
#include <kernel/util/tiny_algebra.hpp>
 
#include <memory>
#include <optional>
#include <utility>
 
namespace FEAT::Geometry
{
 
  template<typename AdaptiveMeshType_>
  class AdaptiveVertexSet
  {
  public:
    using AdaptiveMeshType = AdaptiveMeshType_;
 
    static constexpr int num_coords = AdaptiveMeshType::world_dim;
 
    using CoordType = typename AdaptiveMeshType::CoordType;
 
    using VertexType = typename AdaptiveMeshType::VertexType;
 
  private:
    std::shared_ptr<AdaptiveMeshType> _mesh;
 
    // Mesh layer we refer to
    Layer _layer;
 
    // Vertex permutation
    // Stored, because we can not apply it to the vertex set directly, as the
    // underlying memory is shared between all AdaptiveVertexSets.
    // Stored as an optional to avoid allocating a permutation if no permutation
    // is required by the user.
    std::optional<Adjacency::Permutation> _perm;
 
  public:
    AdaptiveVertexSet(std::shared_ptr<AdaptiveMeshType> mesh, Layer layer) : _mesh(std::move(mesh)), _layer(layer)
    {
    }
 
    AdaptiveVertexSet clone() const
    {
      return AdaptiveVertexSet(_mesh, _layer);
    }
 
    std::size_t bytes() const
    {
      // We are just a view into data stored elsewhere.
      return 0;
    }
 
    int get_num_coords() const
    {
      return num_coords;
    }
 
    Index get_num_vertices() const
    {
      return _mesh->get_num_entities(_layer, 0);
    }
 
    VertexType& operator[](Index idx)
    {
      if(_perm)
      {
        return _mesh->vertex(_layer, _perm.value().map(idx));
      }
      return _mesh->vertex(_layer, idx);
    }
 
    const VertexType& operator[](Index idx) const
    {
      if(_perm)
      {
        return _mesh->vertex(_layer, _perm.value().map(idx));
      }
      return _mesh->vertex(_layer, idx);
    }
 
    // NOTE: Unsupported on purpose. Transform underlying mesh if required.
    //void transform(const VertexType& origin, const VertexType& angles, const VertexType& offset)
 
    void permute(const Adjacency::Permutation& perm, bool invert = false)
    {
      if(!_perm)
      {
        _perm = std::optional{invert ? perm.inverse() : perm.clone()};
      }
      else
      {
        if(invert)
        {
          _perm.value().concat(perm.inverse());
        }
        else
        {
          _perm.value().concat(perm);
        }
      }
    }
 
    static String name()
    {
      return "AdaptiveVertexSet<...>";
    }
  };
 
  template<typename AdaptiveMeshType_, int cell_dim_, int face_dim_>
  class AdaptiveIndexTuple
  {
  public:
    using AdaptiveMeshType = AdaptiveMeshType_;
 
  private:
    std::shared_ptr<AdaptiveMeshType> _mesh;
    Layer _layer;
    Index _entity_idx;
    std::optional<Adjacency::Permutation> _inv_perm;
 
  public:
    // Wrapper Code
 
    static constexpr int num_indices = Shape::FaceTraits<Shape::Hypercube<cell_dim_>, face_dim_>::count;
 
    AdaptiveIndexTuple(std::shared_ptr<AdaptiveMeshType> mesh, Layer layer, Index entity_idx) :
      _mesh(std::move(mesh)),
      _layer(layer),
      _entity_idx(entity_idx),
      _inv_perm(std::nullopt)
    {
    }
 
    AdaptiveIndexTuple(const AdaptiveIndexTuple& other) = default;
    AdaptiveIndexTuple(AdaptiveIndexTuple&& other) = default;
 
    AdaptiveIndexTuple& operator=(const AdaptiveIndexTuple& other) = default;
    AdaptiveIndexTuple& operator=(AdaptiveIndexTuple&& other) = default;
 
    Index operator[](int idx) const
    {
      ASSERT(idx >= 0);
      ASSERT(idx < num_indices);
      Index result = _mesh->template get_face_index<cell_dim_, face_dim_>(_layer, _entity_idx, idx);
      if(_inv_perm)
      {
        return _inv_perm.value().map(result);
      }
      return result;
    }
 
    void permute_map(const Adjacency::Permutation& inv_perm)
    {
      if(_inv_perm)
      {
        _inv_perm.value().concat(inv_perm);
      }
      else
      {
        _inv_perm = inv_perm.clone();
      }
    }
  };
 
  template<typename AdaptiveMeshType_, int cell_dim_, int face_dim_>
  class AdaptiveIndexSet
  {
  public:
    using AdaptiveMeshType = AdaptiveMeshType_;
 
    static constexpr int num_indices = Shape::FaceTraits<Shape::Hypercube<cell_dim_>, face_dim_>::count;
 
    class AdaptiveIndexSetAdjactor
    {
      // Adaptive mesh reference
      std::shared_ptr<AdaptiveMeshType> _mesh;
      Layer _layer;
 
      // Iterator state
      Index _domain_node;
      Index _next;
 
    public:
      AdaptiveIndexSetAdjactor(std::shared_ptr<AdaptiveMeshType> mesh, Layer layer, Index domain_node, Index next) :
        _mesh(std::move(mesh)),
        _layer(layer),
        _domain_node(domain_node),
        _next(next)
      {
      }
 
      bool operator==(const AdaptiveIndexSetAdjactor& other) const
      {
        return _domain_node == other._domain_node && _next == other._next && _layer == other._layer &&
               _mesh == other._mesh;
      }
 
      bool operator!=(const AdaptiveIndexSetAdjactor& other) const
      {
        return !(*this == other);
      }
 
      AdaptiveIndexSetAdjactor& operator++()
      {
        _next++;
        return *this;
      }
 
      Index operator*() const
      {
        return _mesh->template get_face_index<cell_dim_, face_dim_>(_layer, _domain_node, _next);
      }
    };
 
    using ImageIterator = AdaptiveIndexSetAdjactor;
 
    using IndexTupleType = AdaptiveIndexTuple<AdaptiveMeshType, cell_dim_, face_dim_>;
 
  private:
    // Mesh reference
    std::shared_ptr<AdaptiveMeshType> _mesh;
    Layer _layer;
 
    std::optional<Adjacency::Permutation> _permutation;
    std::optional<Adjacency::Permutation> _inverse_face_permutation;
 
  public:
    AdaptiveIndexSet(std::shared_ptr<AdaptiveMeshType> mesh, Layer layer) : _mesh(std::move(mesh)), _layer(layer)
    {
    }
 
    AdaptiveIndexSet(const AdaptiveIndexSet& other) = default;
    AdaptiveIndexSet(AdaptiveIndexSet&& other) = default;
 
    AdaptiveIndexSet& operator=(const AdaptiveIndexSet& other) = default;
    AdaptiveIndexSet& operator=(AdaptiveIndexSet&& other) = default;
 
    AdaptiveIndexSet clone() const
    {
      return AdaptiveIndexSet(_mesh, _layer);
    }
 
    std::size_t bytes() const
    {
      return 0;
    }
 
    Index operator()(Index i, int j)
    {
      if(_permutation)
      {
        return _mesh->template get_face_index<cell_dim_, face_dim_>(
          _layer,
          _permutation.value().map(i),
          _inverse_face_permutation.value().map(j)
        );
      }
      return _mesh->template get_face_index<cell_dim_, face_dim_>(_layer, i, j);
    }
 
    Index operator()(Index i, int j) const
    {
      if(_permutation)
      {
        return _mesh->template get_face_index<cell_dim_, face_dim_>(
          _layer,
          _permutation.value().map(i),
          _inverse_face_permutation.value().map(j)
        );
      }
      return _mesh->template get_face_index<cell_dim_, face_dim_>(_layer, i, j);
    }
 
    IndexTupleType operator[](Index i)
    {
      ASSERT(i < get_num_entities());
      if(_permutation)
      {
        auto tuple = IndexTupleType(_mesh, _layer, _permutation.value().map(i));
        tuple.permute_map(_inverse_face_permutation.value());
        return tuple;
      }
      return IndexTupleType(_mesh, _layer, i);
    }
 
    IndexTupleType operator[](Index i) const
    {
      ASSERT(i < get_num_entities());
      if(_permutation)
      {
        auto tuple = IndexTupleType(_mesh, _layer, _permutation.value().map(i));
        tuple.permute_map(_inverse_face_permutation.value());
        return tuple;
      }
      return IndexTupleType(_mesh, _layer, i);
    }
 
    Index get_num_entities() const
    {
      return _mesh->get_num_entities(_layer, cell_dim_);
    }
 
    int get_num_indices() const
    {
      return num_indices;
    }
 
    Index get_index_bound() const
    {
      return _mesh->get_num_entities(_layer, face_dim_);
    }
 
    // Unsupported.
    // IndexSet API depends on pointer arithmetic that we can not support.
    // Rewrite IndexSet API first, if this is required somewhere
    //AdaptiveIndexTuple* get_indices()
    //const AdaptiveIndexTuple* get_indices() const
 
    void permute(const Adjacency::Permutation& perm, const Adjacency::Permutation& inv_perm_face)
    {
      // Permute index tuples
      if(_permutation)
      {
        _permutation.value().concat(perm);
      }
      else
      {
        _permutation = perm.clone();
      }
 
      // Permute indices by the inverse face permutation
      if(_inverse_face_permutation)
      {
        _inverse_face_permutation.value().concat(inv_perm_face);
      }
      else
      {
        _inverse_face_permutation = inv_perm_face.clone();
      }
    }
 
    static String name()
    {
      return "AdaptiveIndexSet<" + stringify(cell_dim_) + ", " + stringify(face_dim_) + ">";
    }
 
    // Adjactor Interface
 
    Index get_num_nodes_domain() const
    {
      return get_num_entities();
    }
 
    Index get_num_nodes_image() const
    {
      return _mesh->get_num_entities(_layer, face_dim_);
    }
 
    AdaptiveIndexSetAdjactor image_begin(Index domain_node) const
    {
      return AdaptiveIndexSetAdjactor(_mesh, _layer, domain_node, 0);
    }
 
    AdaptiveIndexSetAdjactor image_end(Index domain_node) const
    {
      return AdaptiveIndexSetAdjactor(_mesh, _layer, domain_node, num_indices);
    }
  };
 
  template<
    typename AdaptiveMeshType_,
    typename Shape_,
    int face_dim_ = Shape_::dimension - 1>
  class AdaptiveIndexSetWrapper :
    public AdaptiveIndexSetWrapper<AdaptiveMeshType_, Shape_, face_dim_ - 1>
  {
    static_assert(face_dim_ < Shape_::dimension, "invalid face dimension");
    static_assert(face_dim_ > 0, "invalid face dimension");
 
  public:
    using AdaptiveMeshType = AdaptiveMeshType_;
 
    static constexpr int cell_dim = Shape_::dimension;
 
    using BaseClass =
      AdaptiveIndexSetWrapper<AdaptiveMeshType_, Shape_, face_dim_ - 1>;
 
    using IndexSetType =
      AdaptiveIndexSet<AdaptiveMeshType_, Shape_::dimension, face_dim_>;
 
 
  protected:
    template<int face_dim__>
    using AdaptiveIndexSetByFaceDim =
      AdaptiveIndexSet<AdaptiveMeshType, cell_dim, face_dim__>;
 
    std::shared_ptr<AdaptiveMeshType> _mesh;
    Layer _layer;
 
    IndexSetType _index_set;
 
  public:
    AdaptiveIndexSetWrapper(std::shared_ptr<AdaptiveMeshType> mesh, Layer layer) :
      BaseClass(mesh, layer),
      _mesh(mesh),
      _layer(layer),
      _index_set(mesh, layer)
    {
    }
 
    AdaptiveIndexSetWrapper(const AdaptiveIndexSetWrapper& other) = default;
    AdaptiveIndexSetWrapper(AdaptiveIndexSetWrapper&& other) = default;
 
    AdaptiveIndexSetWrapper& operator=(const AdaptiveIndexSetWrapper& other) = default;
    AdaptiveIndexSetWrapper& operator=(AdaptiveIndexSetWrapper&& other) = default;
 
    void clone(const AdaptiveIndexSetWrapper& other)
    {
      BaseClass::clone(other);
      _index_set = other._index_set;
    }
 
    AdaptiveIndexSetWrapper clone() const
    {
      return AdaptiveIndexSetWrapper(_mesh, _layer);
    }
 
    template<int face_dim__>
    AdaptiveIndexSetByFaceDim<face_dim__>& get_index_set()
    {
      static_assert(face_dim__ >= 0, "invalid face dimension");
      static_assert(face_dim__ < Shape_::dimension, "invalid face dimension");
 
      return AdaptiveIndexSetWrapper<AdaptiveMeshType, Shape_, face_dim__>::_index_set;
    }
 
    template<int face_dim__>
    const AdaptiveIndexSetByFaceDim<face_dim__>& get_index_set() const
    {
      static_assert(face_dim__ >= 0, "invalid face dimension");
      static_assert(face_dim__ < Shape_::dimension, "invalid face dimension");
 
      return AdaptiveIndexSetWrapper<AdaptiveMeshType, Shape_, face_dim__>::_index_set;
    }
 
    template<std::size_t np_>
    void permute(const Adjacency::Permutation& shape_perm,
      const std::array<Adjacency::Permutation, np_>& inv_perm)
    {
      BaseClass::permute(shape_perm, inv_perm);
      _index_set.permute(shape_perm, inv_perm.at(face_dim_));
    }
 
    static String name()
    {
      return "AdaptiveIndexSetWrapper<" + Shape_::name() + ", " + stringify(face_dim_) + ">";
    }
 
    std::size_t bytes()
    {
      return BaseClass::bytes() + _index_set.bytes();
    }
  };
 
  template<typename AdaptiveMeshType_, typename Shape_>
  class AdaptiveIndexSetWrapper<AdaptiveMeshType_, Shape_, 0>
  {
    static_assert(Shape_::dimension > 0, "invalid shape dimension");
 
  public:
    using AdaptiveMeshType = AdaptiveMeshType_;
 
    static constexpr int cell_dim = Shape_::dimension;
 
    using IndexSetType =
      AdaptiveIndexSet<AdaptiveMeshType_, Shape_::dimension, 0>;
 
 
  protected:
    template<int face_dim__>
    using AdaptiveIndexSetByFaceDim =
      AdaptiveIndexSet<AdaptiveMeshType, cell_dim, face_dim__>;
 
    std::shared_ptr<AdaptiveMeshType> _mesh;
    Layer _layer;
 
    IndexSetType _index_set;
 
  public:
    AdaptiveIndexSetWrapper(std::shared_ptr<AdaptiveMeshType> mesh, Layer layer) :
      _mesh(mesh),
      _layer(layer),
      _index_set(mesh, layer)
    {
    }
 
    AdaptiveIndexSetWrapper(const AdaptiveIndexSetWrapper& other) = default;
    AdaptiveIndexSetWrapper(AdaptiveIndexSetWrapper&& other) = default;
 
    AdaptiveIndexSetWrapper& operator=(const AdaptiveIndexSetWrapper& other) = default;
    AdaptiveIndexSetWrapper& operator=(AdaptiveIndexSetWrapper&& other) = default;
 
    void clone(const AdaptiveIndexSetWrapper& other)
    {
      _index_set = other._index_set;
    }
 
    AdaptiveIndexSetWrapper clone() const
    {
      return AdaptiveIndexSetWrapper(_mesh, _layer);
    }
 
    template<int face_dim__>
    AdaptiveIndexSetByFaceDim<face_dim__>& get_index_set()
    {
      static_assert(face_dim__ == 0, "invalid face dimension");
 
      return AdaptiveIndexSetWrapper<AdaptiveMeshType, Shape_, face_dim__>::_index_set;
    }
 
    template<int face_dim__>
    const AdaptiveIndexSetByFaceDim<face_dim__>& get_index_set() const
    {
      static_assert(face_dim__ == 0, "invalid face dimension");
 
      return AdaptiveIndexSetWrapper<AdaptiveMeshType, Shape_, face_dim__>::_index_set;
    }
 
    template<std::size_t np_>
    void permute(const Adjacency::Permutation& shape_perm,
      const std::array<Adjacency::Permutation, np_>& inv_perm)
    {
      _index_set.permute(shape_perm, inv_perm.at(0));
    }
 
    static String name()
    {
      return "AdaptiveIndexSetWrapper<" + Shape_::name() + ",0>";
    }
 
    std::size_t bytes() const
    {
      return _index_set.bytes();
    }
  };
 
  template<typename AdaptiveMeshType_, typename Shape_>
  class AdaptiveIndexSetHolder :
    public AdaptiveIndexSetHolder<
      AdaptiveMeshType_,
      typename Shape::FaceTraits<Shape_, Shape_::dimension - 1>::ShapeType>
  {
    using AdaptiveMeshType = AdaptiveMeshType_;
 
    using BaseClass = AdaptiveIndexSetHolder<
      AdaptiveMeshType_,
      typename Shape::FaceTraits<Shape_, Shape_::dimension - 1>::ShapeType>;
 
    using IndexSetWrapperType = AdaptiveIndexSetWrapper<AdaptiveMeshType, Shape_>;
 
  protected:
    template<int shape_dim_>
    using AdaptiveIndexSetWrapperByShapeDim =
      AdaptiveIndexSetWrapper<
        AdaptiveMeshType,
        typename Shape::FaceTraits<Shape_, shape_dim_>::ShapeType>;
 
    std::shared_ptr<AdaptiveMeshType> _mesh;
    Layer _layer;
 
    IndexSetWrapperType _index_set_wrapper;
 
  public:
    AdaptiveIndexSetHolder(std::shared_ptr<AdaptiveMeshType> mesh, Layer layer) :
      BaseClass(mesh, layer),
      _mesh(mesh),
      _layer(layer),
      _index_set_wrapper(mesh, layer)
    {
    }
 
    AdaptiveIndexSetHolder(const AdaptiveIndexSetHolder& other) = default;
    AdaptiveIndexSetHolder(AdaptiveIndexSetHolder&& other) = default;
 
    AdaptiveIndexSetHolder& operator=(const AdaptiveIndexSetHolder& other) = default;
    AdaptiveIndexSetHolder& operator=(AdaptiveIndexSetHolder&& other) = default;
 
    void clone(const AdaptiveIndexSetHolder& other)
    {
      BaseClass::clone(other);
      _index_set_wrapper.clone(other._index_set_wrapper);
    }
 
    AdaptiveIndexSetHolder clone() const
    {
      return AdaptiveIndexSetHolder(_mesh, _layer);
    }
 
    template<int shape_dim_>
    AdaptiveIndexSetWrapperByShapeDim<shape_dim_>& get_index_set_wrapper()
    {
      static_assert(shape_dim_ > 0, "invalid shape dimension");
      static_assert(shape_dim_ <= Shape_::dimension, "invalid shape dimension");
      typedef typename Shape::FaceTraits<Shape_, shape_dim_>::ShapeType ShapeType;
      return AdaptiveIndexSetHolder<AdaptiveMeshType, ShapeType>::_index_set_wrapper;
    }
 
    template<int shape_dim_>
    const AdaptiveIndexSetWrapperByShapeDim<shape_dim_>& get_index_set_wrapper() const
    {
      static_assert(shape_dim_ > 0, "invalid shape dimension");
      static_assert(shape_dim_ <= Shape_::dimension, "invalid shape dimension");
      typedef typename Shape::FaceTraits<Shape_, shape_dim_>::ShapeType ShapeType;
      return AdaptiveIndexSetHolder<AdaptiveMeshType, ShapeType>::_index_set_wrapper;
    }
 
    template<int cell_dim_, int face_dim_>
    AdaptiveIndexSet<AdaptiveMeshType, cell_dim_, face_dim_>& get_index_set()
    {
      static_assert(cell_dim_ <= Shape_::dimension, "invalid cell dimension");
      static_assert(face_dim_ < cell_dim_, "invalid face/cell dimension");
      static_assert(face_dim_ >= 0, "invalid face dimension");
      return get_index_set_wrapper<cell_dim_>().template get_index_set<face_dim_>();
    }
 
    template<int cell_dim_, int face_dim_>
    const AdaptiveIndexSet<AdaptiveMeshType, cell_dim_, face_dim_>& get_index_set() const
    {
      static_assert(cell_dim_ <= Shape_::dimension, "invalid cell dimension");
      static_assert(face_dim_ < cell_dim_, "invalid face/cell dimension");
      static_assert(face_dim_ >= 0, "invalid face dimension");
      return get_index_set_wrapper<cell_dim_>().template get_index_set<face_dim_>();
    }
 
    template<std::size_t np_>
    void permute(const std::array<Adjacency::Permutation, np_>& perms,
      const std::array<Adjacency::Permutation, np_>& inv_perms)
    {
      BaseClass::permute(perms, inv_perms);
      _index_set_wrapper.permute(perms.at(Shape_::dimension), inv_perms);
    }
 
    static String name()
    {
      return "AdaptiveIndexSetHolder<" + Shape_::name() + ">";
    }
 
    std::size_t bytes() const
    {
      return BaseClass::bytes() + _index_set_wrapper.bytes();
    }
  };
 
  template<typename AdaptiveMeshType_>
  class AdaptiveIndexSetHolder<AdaptiveMeshType_, Shape::Vertex>
  {
  public:
    AdaptiveIndexSetHolder() = default;
 
    AdaptiveIndexSetHolder(std::shared_ptr<AdaptiveMeshType_> /*mesh*/, Layer /*layer*/){
    }
 
    AdaptiveIndexSetHolder(const AdaptiveIndexSetHolder& other) = default;
    AdaptiveIndexSetHolder(AdaptiveIndexSetHolder&& other) = default;
 
    AdaptiveIndexSetHolder& operator=(const AdaptiveIndexSetHolder& other) = default;
    AdaptiveIndexSetHolder& operator=(AdaptiveIndexSetHolder&& other) = default;
 
    void clone(const AdaptiveIndexSetHolder& /*unused*/)
    {
    }
 
    AdaptiveIndexSetHolder clone() const
    {
      return AdaptiveIndexSetHolder();
    }
 
    template<std::size_t np_>
    void permute(const std::array<Adjacency::Permutation, np_>&,
      const std::array<Adjacency::Permutation, np_>&)
    {
    }
 
    static String name()
    {
      return "AdaptiveIndexSetHolder<Vertex>";
    }
 
    std::size_t bytes() const
    {
      return std::size_t(0);
    }
  };
 
  template<typename AdaptiveMeshType_>
  class AdaptiveNeighbors
  {
  protected:
    using AdaptiveMeshType = AdaptiveMeshType_;
 
    std::shared_ptr<AdaptiveMeshType> _mesh;
    Layer _layer;
 
  public:
 
    static constexpr int num_indices =
      Shape::FaceTraits<typename AdaptiveMeshType::ShapeType, AdaptiveMeshType_::shape_dim - 1>::count;
 
    class AdaptiveNeighborsAdjactor
    {
      // Adaptive mesh reference
      std::shared_ptr<AdaptiveMeshType> _mesh;
      Layer _layer;
 
      // Iterator state
      Index _domain_node;
      Index _next;
 
      std::optional<Adjacency::Permutation> _inverse_face_permutation = std::nullopt;
 
    public:
      AdaptiveNeighborsAdjactor(std::shared_ptr<AdaptiveMeshType> mesh, Layer layer, Index domain_node, Index next) :
        _mesh(std::move(mesh)),
        _layer(layer),
        _domain_node(domain_node),
        _next(next)
      {
      }
 
      AdaptiveNeighborsAdjactor(
        std::shared_ptr<AdaptiveMeshType> mesh,
        Layer layer,
        Index domain_node,
        Index next,
        Adjacency::Permutation inverse_face_permutation) :
        _mesh(std::move(mesh)),
        _layer(layer),
        _domain_node(domain_node),
        _next(next),
        _inverse_face_permutation(std::move(inverse_face_permutation))
      {
      }
 
      bool operator==(const AdaptiveNeighborsAdjactor& other) const
      {
        return _domain_node == other._domain_node && _next == other._next && _layer == other._layer &&
               _mesh == other._mesh;
      }
 
      bool operator!=(const AdaptiveNeighborsAdjactor& other) const
      {
        return !(*this == other);
      }
 
      AdaptiveNeighborsAdjactor& operator++()
      {
        _next++;
        return *this;
      }
 
      Index operator*() const
      {
        if(_inverse_face_permutation)
        {
          return _mesh->get_neighbor(_layer, _domain_node, _inverse_face_permutation.value().map(_next));
        }
        return _mesh->get_neighbor(_layer, _domain_node, _next);
      }
    };
 
    using ImageIterator = AdaptiveNeighborsAdjactor;
 
  private:
    // Mesh reference
    std::optional<Adjacency::Permutation> _permutation;
    std::optional<Adjacency::Permutation> _inverse_face_permutation;
 
  public:
    AdaptiveNeighbors(std::shared_ptr<AdaptiveMeshType> mesh, Layer layer) : _mesh(std::move(mesh)), _layer(layer)
    {
    }
 
    AdaptiveNeighbors(const AdaptiveNeighbors& other) = default;
    AdaptiveNeighbors(AdaptiveNeighbors&& other) = default;
 
    AdaptiveNeighbors& operator=(const AdaptiveNeighbors& other) = default;
    AdaptiveNeighbors& operator=(AdaptiveNeighbors&& other) = default;
 
    AdaptiveNeighbors clone() const
    {
      return AdaptiveNeighbors(_mesh, _layer);
    }
 
    std::size_t bytes() const
    {
      return 0;
    }
 
    Index operator()(Index i, int j)
    {
      if(_permutation)
      {
        return _mesh->get_neighbor(
          _layer,
          _permutation.value().map(i),
          Index(_inverse_face_permutation.value().map(j))
        );
      }
      return _mesh->get_neighbor(_layer, i, Index(j));
    }
 
    Index operator()(Index i, int j) const
    {
      if(_permutation)
      {
        return _mesh->get_neighbor(
          _layer,
          _permutation.value().map(i),
          Index(_inverse_face_permutation.value().map(j))
        );
      }
      return _mesh->get_neighbor(_layer, i, Index(j));
    }
 
    IndexTuple<num_indices> operator[](Index i) const
    {
      IndexTuple<num_indices> result;
      for(int j(0); j < num_indices; j++)
      {
        result[j] = (*this)(i, j);
      }
      return result;
    }
 
    Index get_num_entities() const
    {
      return _mesh->get_num_entities(_layer, AdaptiveMeshType::ShapeType::dimension);
    }
 
    int get_num_indices() const
    {
      return num_indices;
    }
 
    Index get_index_bound() const
    {
      return _mesh->get_num_entities(_layer, AdaptiveMeshType::ShapeType::dimension);
    }
 
    // Unsupported.
    // IndexSet API depends on pointer arithmetic that we can not support.
    // Rewrite IndexSet API first, if this is required somewhere
    //AdaptiveIndexTuple* get_indices()
    //const AdaptiveIndexTuple* get_indices() const
 
    void permute(const Adjacency::Permutation& perm, const Adjacency::Permutation& inv_perm_face)
    {
      // Permute index tuples
      if(_permutation)
      {
        _permutation.value().concat(perm);
      }
      else
      {
        _permutation = perm.clone();
      }
 
      // Permute indices by the inverse face permutation
      if(_inverse_face_permutation)
      {
        _inverse_face_permutation.value().concat(inv_perm_face);
      }
      else
      {
        _inverse_face_permutation = inv_perm_face.clone();
      }
    }
 
    static String name()
    {
      return "AdaptiveNeighbors";
    }
 
    // Adjactor Interface
 
    Index get_num_nodes_domain() const
    {
      return get_num_entities();
    }
 
    Index get_num_nodes_image() const
    {
      return _mesh->get_num_entities(_layer, AdaptiveMeshType::ShapeType::dimension);
    }
 
    AdaptiveNeighborsAdjactor image_begin(Index domain_node) const
    {
      if(_permutation)
      {
        return AdaptiveIndexSetAdjactor(_mesh, _layer, _permutation.value().map(domain_node), 0, _inverse_face_permutation.value());
      }
      return AdaptiveIndexSetAdjactor(_mesh, _layer, domain_node, 0);
    }
 
    AdaptiveNeighborsAdjactor image_end(Index domain_node) const
    {
      if(_permutation)
      {
        return AdaptiveIndexSetAdjactor(_mesh, _layer, _permutation.value().map(domain_node), num_indices, _inverse_face_permutation.value());
      }
      return AdaptiveIndexSetAdjactor(_mesh, _layer, domain_node, num_indices);
    }
  };
 
  template<typename AdaptiveMeshType_>
  class AdaptiveMeshLayer
  {
  public:
    using AdaptiveMeshType = AdaptiveMeshType_;
 
    using AdaptiveNeighborsType = AdaptiveNeighbors<AdaptiveMeshType>;
 
    // Conformal Mesh Interface
 
    using ShapeType = typename AdaptiveMeshType::ShapeType;
 
    using VertexSetType = AdaptiveVertexSet<AdaptiveMeshType>;
 
    using VertexType = typename AdaptiveMeshType::VertexType;
 
    using CoordType = typename AdaptiveMeshType::CoordType;
 
    using IndexSetHolderType = AdaptiveIndexSetHolder<AdaptiveMeshType, ShapeType>;
 
    static constexpr int shape_dim = AdaptiveMeshType::shape_dim;
 
    static constexpr int world_dim = AdaptiveMeshType::world_dim;
 
    static constexpr bool is_structured = false;
 
    using MeshPermutationType = MeshPermutation<ShapeType>;
 
    template<int cell_dim_, int face_dim_>
    struct IndexSet
    {
      using Type = AdaptiveIndexSet<AdaptiveMeshType, cell_dim_, face_dim_>;
    };
 
  private:
    std::shared_ptr<AdaptiveMeshType> _mesh;
    Layer _layer;
 
    VertexSetType _vertex_set;
    IndexSetHolderType _index_set_holder;
 
    MeshPermutation<ShapeType> _permutation;
 
  public:
    AdaptiveMeshLayer(std::shared_ptr<AdaptiveMeshType> mesh, Layer layer) :
      _mesh(mesh),
      _layer(layer),
      _vertex_set(VertexSetType(_mesh, layer)),
      _index_set_holder(_mesh, layer),
      _permutation()
    {
    }
 
    AdaptiveMeshLayer(const AdaptiveMeshLayer& other) = default;
    AdaptiveMeshLayer(AdaptiveMeshLayer&& other) = default;
 
    AdaptiveMeshLayer& operator=(const AdaptiveMeshLayer& other) = default;
    AdaptiveMeshLayer& operator=(AdaptiveMeshLayer&& other) = default;
 
    void clone(const AdaptiveMeshLayer& other)
    {
      _mesh = other._mesh;
      _layer = other._layer;
      _vertex_set = other._vertex_set;
      _index_set_holder = other._index_set_holder;
      _permutation = other._permutation.clone();
    }
 
    AdaptiveMeshLayer clone() const
    {
      AdaptiveMeshLayer result(_mesh, _layer);
      result._vertex_set = _vertex_set;
      result._index_set_holder = _index_set_holder;
      result._permutation = _permutation;
      return result;
    }
 
    std::size_t bytes() const
    {
      return _vertex_set.bytes() + _index_set_holder.bytes() + _permutation.bytes();
    }
 
    Index get_num_entities(int dim) const
    {
      return _mesh->get_num_entities(_layer, dim);
    }
 
    Index get_num_vertices() const
    {
      return get_num_entities(0);
    }
 
    Index get_num_elements() const
    {
      return get_num_entities(ShapeType::dimension);
    }
 
    bool is_permuted() const
    {
      return !_permutation.empty();
    }
 
    const MeshPermutation<ShapeType>& get_mesh_permutation() const
    {
      return _permutation;
    }
 
    void create_permutation(PermutationStrategy strategy)
    {
      // make sure that we don't already have a permutation
      XASSERTM(this->_permutation.empty(), "mesh is already permuted!");
 
      // create the permutation
      this->_permutation.create(strategy, this->_index_set_holder, this->_vertex_set);
 
      // permute vertex set
      this->_vertex_set.permute(this->_permutation.get_perm(0));
 
      // permute index sets
      this->_index_set_holder.permute(this->_permutation.get_perms(), this->_permutation.get_inv_perms());
    }
 
    void set_permutation(MeshPermutationType&& mesh_perm)
    {
      // make sure that we don't already have a permutation
      XASSERTM(this->_permutation.empty(), "mesh is already permuted!");
 
      // check the dimensions
      const std::array<Index, 4> mesh_size {
        get_num_entities(0),
        get_num_entities(1),
        get_num_entities(2),
        get_num_entities(3)
      };
      XASSERTM(mesh_perm.validate_sizes(mesh_size.data()) == 0, "mesh permutation has invalid size!");
 
      // save the permutation
      this->_permutation = std::forward<MeshPermutationType>(mesh_perm);
 
      // permute vertex set
      this->_vertex_set.permute(this->_permutation.get_perm(0));
 
      // permute index sets
      this->_index_set_holder.permute(this->_permutation.get_perms(), this->_permutation.get_inv_perms());
    }
 
    bool validate_element_coloring() const
    {
      // no coloring?
      const std::vector<Index>& coloring = this->get_mesh_permutation().get_element_coloring();
      if(coloring.empty())
        return true;
 
      // get vertices-at-element index set
      const auto& verts_at_elem = this->template get_index_set<shape_dim, 0>();
 
      // render transpose
      Adjacency::Graph elems_at_vert(Adjacency::RenderType::transpose, verts_at_elem);
 
      // loop over all color blocks
      for(std::size_t icol(0); icol+1u < coloring.size(); ++icol)
      {
        // get the bounds of our current color block
        const Index iel_beg = coloring[icol];
        const Index iel_end = coloring[icol+1u];
 
        // loop over all elements in the current color block
        for(Index iel(iel_beg); iel < iel_end; ++iel)
        {
          // loop over all vertices adjacent to this element
          for(int ivt(0); ivt < verts_at_elem.num_indices; ++ivt)
          {
            // loop over all elements adjacent to this vertex
            const Index ivtx = verts_at_elem(iel, ivt);
            for(auto it = elems_at_vert.image_begin(ivtx); it != elems_at_vert.image_end(ivtx); ++it)
            {
              // two adjacent element must not be in the same color block
              if((iel_beg <= *it) && (*it < iel_end) && (*it != iel))
                return false; // invalid coloring
            }
          }
        }
      } // next color block
 
      // ok, coloring is valid
      return true;
    }
 
    bool validate_element_layering() const
    {
      // no layering?
      const std::vector<Index>& layering = this->get_mesh_permutation().get_element_layering();
      if(layering.empty())
        return true;
 
      // get vertices-at-element index set
      const auto& verts_at_elem = this->template get_index_set<shape_dim, 0>();
 
      // render transpose
      Adjacency::Graph elems_at_vert(Adjacency::RenderType::transpose, verts_at_elem);
 
      // loop over all layers
      for(std::size_t ilay(0); ilay+1u < layering.size(); ++ilay)
      {
        // get the bounds of our current layer
        const Index iel_beg = layering[ilay];
        const Index iel_end = layering[ilay+1u];
 
        // get the lower bound for valid neighbors of our current layer = beginning of previous layer
        const Index iel_lower = layering[Math::max(ilay, std::size_t(1)) - 1u];
 
        // get the upper bound for valid neighbors of our current layer = end of next layer
        const Index iel_upper = layering[Math::min(ilay+2u, layering.size()-1u)];
 
        // loop over all elements in the current layer
        for(Index iel(iel_beg); iel < iel_end; ++iel)
        {
          // loop over all vertices adjacent to this element
          for(int ivt(0); ivt < verts_at_elem.num_indices; ++ivt)
          {
            // loop over all elements adjacent to this vertex
            const Index ivtx = verts_at_elem(iel, ivt);
            for(auto it = elems_at_vert.image_begin(ivtx); it != elems_at_vert.image_end(ivtx); ++it)
            {
              // adjacent element outside of adjacent layers?
              if(!((iel_lower <= *it) && (*it < iel_upper)))
                return false; // invalid layer
            }
          }
        }
      } // next color block
 
      // ok, layering is valid
      return true;
    }
 
    void fill_neighbours()
    {
      _mesh->fill_neighbors();
    }
 
    AdaptiveNeighborsType get_neighbors()
    {
      return AdaptiveNeighbors(_mesh, _layer);
    }
 
    AdaptiveNeighborsType get_neighbors() const
    {
      return AdaptiveNeighbors(_mesh, _layer);
    }
 
    VertexSetType& get_vertex_set()
    {
      return _vertex_set;
    }
 
    const VertexSetType& get_vertex_set() const
    {
      return _vertex_set;
    }
 
    template<int cell_dim_, int face_dim_>
    auto& get_index_set()
    {
      return _index_set_holder.template get_index_set<cell_dim_, face_dim_>();
    }
 
    template<int cell_dim_, int face_dim_>
    const auto& get_index_set() const
    {
      return _index_set_holder.template get_index_set<cell_dim_, face_dim_>();
    }
 
    IndexSetHolderType& get_index_set_holder()
    {
      return _index_set_holder;
    }
 
    const IndexSetHolderType& get_index_set_holder() const
    {
      return _index_set_holder;
    }
 
    IndexSetHolderType& get_topology()
    {
      return _index_set_holder;
    }
 
    const IndexSetHolderType& get_topology() const
    {
      return _index_set_holder;
    }
 
    // Wrapper specific code
 
    std::optional<Index> get_child(Index elem, Index child) const
    {
      return _mesh->template get_child<ShapeType::dimension>(_layer, elem, child);
    }
 
    Index get_num_children(Index elem) const
    {
      return _mesh->template get_num_children<ShapeType::dimension>(_layer, elem);
    }
 
    static String name()
    {
      return "AdaptiveMeshLayer<...>";
    }
 
    bool has_vertex_changed(Index vertex_idx) const
    {
      return _mesh->has_vertex_changed(_layer, vertex_idx);
    }
 
    typename AdaptiveMeshType::FoundationMeshType& foundation_mesh()
    {
      return _mesh->foundation_mesh();
    }
 
    const typename AdaptiveMeshType::FoundationMeshType& foundation_mesh() const
    {
      return _mesh->foundation_mesh();
    }
 
    AdaptiveMeshType& adaptive_mesh()
    {
      return *_mesh;
    }
 
    const AdaptiveMeshType& adaptive_mesh() const
    {
      return *_mesh;
    }
 
    std::shared_ptr<AdaptiveMeshType> adaptive_mesh_ptr() const
    {
      return _mesh;
    }
  };
 
  template<typename AdaptiveMeshLayerType>
  class AdaptiveChildMapping
  {
  public:
    class ChildCellIterator
    {
    public:
      ChildCellIterator(const AdaptiveMeshLayerType* coarse, Index cell, Index child) :
        _coarse(coarse),
        _cell(cell),
        _child(child)
      {
      }
 
      ChildCellIterator() = default;
 
      ChildCellIterator(const ChildCellIterator& other) = default;
      ChildCellIterator(ChildCellIterator&& other) = default;
 
      ChildCellIterator& operator=(const ChildCellIterator& other) = default;
      ChildCellIterator& operator=(ChildCellIterator&& other) = default;
 
      ~ChildCellIterator() = default;
 
      Index operator*() const
      {
        return _coarse->get_child(_cell, _child).value_or(0);
      }
 
      ChildCellIterator& operator++()
      {
        _child++;
        return *this;
      }
 
      bool operator!=(const ChildCellIterator& other)
      {
        return _cell != other._cell || _child != other._child;
      }
 
    private:
      static constexpr int dim = AdaptiveMeshLayerType::ShapeType::dimension;
 
      const AdaptiveMeshLayerType* _coarse;
      Index _cell;
      Index _child;
    };
 
    using ImageIterator = ChildCellIterator;
 
  private:
    // TODO: Should these be std::shared_ptrs?
    const AdaptiveMeshLayerType& _coarse;
    const AdaptiveMeshLayerType& _fine;
 
  public:
    AdaptiveChildMapping(const AdaptiveMeshLayerType& coarse, const AdaptiveMeshLayerType& fine) : _coarse(coarse), _fine(fine)
    {
    }
 
    Index get_num_nodes_domain() const
    {
      return _coarse.get_num_elements();
    }
 
    Index get_num_nodes_image() const
    {
      return _fine.get_num_elements();
    }
 
    ChildCellIterator image_begin(Index domain_node) const
    {
      return ChildCellIterator(&_coarse, domain_node, 0);
    }
 
    ChildCellIterator image_end(Index domain_node) const
    {
      Index num_children = _coarse.get_num_children(domain_node);
      return ChildCellIterator(&_coarse, domain_node, num_children);
    }
  };
 
  template<typename AdaptiveMeshType_, int dim_ = AdaptiveMeshType_::ShapeType::dimension>
  class FoundationAdaptiveAdjactor
  {
    const AdaptiveMeshType_& _amesh;
  public:
    class OnceIterator
    {
      Index _idx;
      bool _done;
 
    public:
 
      using iterator_category = std::input_iterator_tag;
      using difference_type = std::size_t;
      using value_type = Index;
      using pointer = Index*;
      using reference = Index&;
 
      OnceIterator() : _idx(0), _done(true)
      {
      }
 
      explicit OnceIterator(Index i) : _idx(i), _done(false)
      {
      }
 
      OnceIterator(const OnceIterator& other) = default;
      OnceIterator(OnceIterator&& other) = default;
 
      OnceIterator& operator=(const OnceIterator& other) = default;
      OnceIterator& operator=(OnceIterator&& other) = default;
 
      ~OnceIterator() = default;
 
 
      Index operator*()
      {
        return _idx;
      }
 
      Index operator->()
      {
        return _idx;
      }
 
      OnceIterator& operator++()
      {
        _done = true;
        return *this;
      }
 
      OnceIterator& operator++(int)
      {
        OnceIterator tmp = *this;
        ++(*this);
        return tmp;
      }
 
      friend bool operator==(const OnceIterator& a, const OnceIterator& b)
      {
        return (a._done && b._done);
      }
 
      friend bool operator!=(const OnceIterator& a, const OnceIterator& b)
      {
        return !(a == b);
      }
    };
 
    using ImageIterator = OnceIterator;
 
    explicit FoundationAdaptiveAdjactor(const AdaptiveMeshType_& amesh) : _amesh(amesh)
    {
    }
 
    Index get_num_nodes_domain() const
    {
      return _amesh.foundation_mesh().get_num_entities(dim_);
    }
 
    Index get_num_nodes_image() const
    {
      return _amesh.get_num_entities(Geometry::Layer{0}, dim_);
    }
 
    ImageIterator image_begin(Index domain_idx) const
    {
      if(auto mapping = _amesh.template get_overlap_cell<dim_>(domain_idx))
      {
        return OnceIterator(*mapping);
      }
      else
      {
        return OnceIterator{};
      }
    }
 
    ImageIterator image_end(Index /*domain_idx*/) const
    {
      return OnceIterator{};
    }
  };
} // namespace FEAT::Geometry