index_set.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/shape.hpp>
#include <kernel/adjacency/permutation.hpp>
 
// includes, system
#include <array>
#include <vector>
 
namespace FEAT
{
  namespace Geometry
  {
    template<int num_indices_>
    struct IndexTuple
    {
      static_assert(num_indices_ > 0, "invalid number of indices");
 
      static constexpr int num_indices = num_indices_;
 
      Index indices[num_indices];
 
      Index& operator[](int i)
      {
        ASSERT(i >= 0);
        ASSERT(i < num_indices);
        return indices[i];
      }
 
      const Index& operator[](int i) const
      {
        ASSERT(i >= 0);
        ASSERT(i < num_indices);
        return indices[i];
      }
 
      void permute_map(const Adjacency::Permutation& inv_perm)
      {
        for(int i(0); i < num_indices; ++i)
          indices[i] = inv_perm.map(indices[i]);
      }
    }; // struct IndexTuple
 
    template<int num_indices_>
    class IndexSet
    {
      static_assert(num_indices_ > 0, "invalid index count");
 
    public:
      static constexpr int num_indices = num_indices_;
 
      typedef IndexTuple<num_indices> IndexTupleType;
 
      typedef const Index* ImageIterator;
 
    protected:
      Index _index_bound;
 
      std::vector<IndexTupleType> _indices;
 
      explicit IndexSet(Index idx_bnd, const std::vector<IndexTupleType>& idx) :
        _index_bound(idx_bnd), _indices(idx)
      {
      }
 
    public:
      explicit IndexSet(
        Index num_entities = 0,
        Index index_bound = 0)
          :
        _index_bound(index_bound),
        _indices(std::size_t(num_entities))
      {
      }
 
      IndexSet(IndexSet&& other) :
        _index_bound(other._index_bound),
        _indices(std::forward<std::vector<IndexTupleType>>(other._indices))
      {
      }
 
      IndexSet& operator=(IndexSet&& other)
      {
        // avoid self-move
        if(this == &other)
          return *this;
 
        _index_bound = other._index_bound;
        _indices = std::forward<std::vector<IndexTupleType>>(other._indices);
 
        return *this;
      }
 
      virtual ~IndexSet()
      {
      }
 
      IndexSet clone() const
      {
        return IndexSet(this->_index_bound, this->_indices);
      }
 
      std::size_t bytes() const
      {
        return _indices.size() * sizeof(IndexTupleType);
      }
 
      int get_num_indices() const
      {
        return num_indices_;
      }
 
      Index get_num_entities() const
      {
        return Index(_indices.size());
      }
 
      Index get_index_bound() const
      {
        return _index_bound;
      }
 
      IndexTupleType* get_indices()
      {
        return _indices.data();
      }
 
      const IndexTupleType* get_indices() const
      {
        return _indices.data();
      }
 
      IndexTupleType& operator[](Index i)
      {
        ASSERT(i < get_num_entities());
        return _indices[i];
      }
 
      const IndexTupleType& operator[](Index i) const
      {
        ASSERT(i < get_num_entities());
        return _indices[i];
      }
 
      Index& operator()(Index i, int j)
      {
        ASSERT(i < get_num_entities());
        ASSERT(j < num_indices);
        return _indices[i][j];
      }
 
      // copy-paste documentation because \copydoc does not work with operators
      const Index& operator()(Index i, int j) const
      {
        ASSERT(i < get_num_entities());
        ASSERT(j < num_indices);
        return _indices[i][j];
      }
 
      void set_index_bound(Index bound)
      {
        _index_bound = bound;
      }
 
      void permute(const Adjacency::Permutation& perm, const Adjacency::Permutation& inv_perm_face)
      {
        if(_indices.empty())
          return;
 
        // first, permute the actual index tuples by the forward permutation
        if(!perm.empty())
        {
          XASSERT(perm.size() == get_num_entities());
          perm.apply(_indices.data());
        }
 
        // now, permute all indices by the inverse face permutation
        if(!inv_perm_face.empty())
        {
          XASSERT(inv_perm_face.size() == get_index_bound());
          const Index n = Index(_indices.size());
          for(Index i(0); i < n; ++i)
            _indices.at(i).permute_map(inv_perm_face);
        }
      }
 
      static String name()
      {
        return "IndexSet<" + stringify(num_indices_) + ">";
      }
 
      /* *************************************************************************************** */
      /*          A D J A C T O R     I N T E R F A C E     I M P L E M E N T A T I O N          */
      /* *************************************************************************************** */
      Index get_num_nodes_domain() const
      {
        return get_num_entities();
      }
 
      Index get_num_nodes_image() const
      {
        return _index_bound;
      }
 
      ImageIterator image_begin(Index domain_node) const
      {
        ASSERT(domain_node < get_num_entities());
        return reinterpret_cast<const Index*>(&_indices.data()[domain_node]);
      }
 
      ImageIterator image_end(Index domain_node) const
      {
        ASSERT(domain_node < get_num_entities());
        return reinterpret_cast<const Index*>(&_indices.data()[domain_node+1]);
      }
    }; // class IndexSet<...>
 
    /* ***************************************************************************************** */
    /* ***************************************************************************************** */
    /* ***************************************************************************************** */
 
    template<
      typename Shape_,
      int face_dim_ = Shape_::dimension - 1>
    class IndexSetWrapper :
      public IndexSetWrapper<Shape_, face_dim_ - 1>
    {
      static_assert(face_dim_ < Shape_::dimension, "invalid face dimension");
      static_assert(face_dim_ > 0, "invalid face dimension");
 
    public:
      typedef IndexSetWrapper<Shape_, face_dim_ - 1> BaseClass;
      typedef IndexSet<Shape::FaceTraits<Shape_, face_dim_>::count> IndexSetType;
 
    protected:
      IndexSetType _index_set;
 
    public:
      IndexSetWrapper() :
        BaseClass(),
        _index_set()
      {
      }
 
      explicit IndexSetWrapper(const Index num_entities[]) :
        BaseClass(num_entities),
        _index_set(
          num_entities[Shape_::dimension],
          num_entities[face_dim_])
      {
      }
 
      template<std::size_t m_>
      explicit IndexSetWrapper(const std::array<Index, m_>& num_entities) :
        BaseClass(num_entities),
        _index_set(
          num_entities[Shape_::dimension],
          num_entities[face_dim_])
      {
        static_assert(int(m_) > Shape_::dimension, "invalid array size");
      }
 
      IndexSetWrapper(IndexSetWrapper&& other) :
        BaseClass(std::forward<BaseClass>(other)),
        _index_set(std::forward<IndexSetType>(other._index_set))
      {
      }
 
      IndexSetWrapper& operator=(IndexSetWrapper&& other)
      {
        if(this == &other)
          return *this;
 
        BaseClass::operator=(std::forward<BaseClass>(other));
        _index_set = std::forward<IndexSetType>(other._index_set);
 
        return *this;
      }
 
      virtual ~IndexSetWrapper()
      {
      }
 
      void clone(const IndexSetWrapper& other)
      {
        BaseClass::clone(other);
        this->_index_set = other._index_set.clone();
      }
 
      IndexSetWrapper clone() const
      {
        IndexSetWrapper isw;
        isw.clone(*this);
        return isw;
      }
 
      template<int face_dim__>
      IndexSet<Shape::FaceTraits<Shape_, face_dim__>::count>& get_index_set()
      {
        static_assert(face_dim__ >= 0, "invalid face dimension");
        static_assert(face_dim__ < Shape_::dimension, "invalid face dimension");
        return IndexSetWrapper<Shape_, face_dim__>::_index_set;
      }
 
      template<int face_dim__>
      const IndexSet<Shape::FaceTraits<Shape_, face_dim__>::count>& get_index_set() const
      {
        static_assert(face_dim__ >= 0, "invalid face dimension");
        static_assert(face_dim__ < Shape_::dimension, "invalid face dimension");
        return IndexSetWrapper<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 "IndexSetWrapper<" + Shape_::name() + "," + stringify(face_dim_) + ">";
      }
 
      std::size_t bytes() const
      {
        return BaseClass::bytes() + _index_set.bytes();
      }
    };
 
    template<typename Shape_>
    class IndexSetWrapper<Shape_, 0>
    {
      static_assert(Shape_::dimension > 0, "invalid shape dimension");
 
    public:
      typedef IndexSet<Shape::FaceTraits<Shape_, 0>::count> IndexSetType;
 
    protected:
      IndexSetType _index_set;
 
    public:
      IndexSetWrapper() :
        _index_set()
      {
      }
 
      explicit IndexSetWrapper(const Index num_entities[]) :
        _index_set(
          num_entities[Shape_::dimension],
          num_entities[0])
      {
      }
 
      template<std::size_t m_>
      explicit IndexSetWrapper(const std::array<Index, m_>& num_entities) :
        _index_set(
          num_entities[Shape_::dimension],
          num_entities[0])
      {
        static_assert(int(m_) > Shape_::dimension, "invalid array size");
      }
 
      IndexSetWrapper(IndexSetWrapper&& other) :
        _index_set(std::forward<IndexSetType>(other._index_set))
      {
      }
 
      IndexSetWrapper& operator=(IndexSetWrapper&& other)
      {
        if(this == &other)
          return *this;
 
        _index_set = std::forward<IndexSetType>(other._index_set);
 
        return *this;
      }
 
      virtual ~IndexSetWrapper()
      {
      }
 
      void clone(const IndexSetWrapper& other)
      {
        this->_index_set = other._index_set.clone();
      }
 
      IndexSetWrapper clone() const
      {
        IndexSetWrapper isw;
        isw.clone(*this);
        return isw;
      }
 
      template<int face_dim__>
      IndexSet<Shape::FaceTraits<Shape_, face_dim__>::count>& get_index_set()
      {
        static_assert(face_dim__ == 0, "invalid face dimension");
        return IndexSetWrapper<Shape_, face_dim__>::_index_set;
      }
 
      template<int face_dim__>
      const IndexSet<Shape::FaceTraits<Shape_, face_dim__>::count>& get_index_set() const
      {
        static_assert(face_dim__ == 0, "invalid face dimension");
        return IndexSetWrapper<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 "IndexSetWrapper<" + Shape_::name() + ",0>";
      }
 
      std::size_t bytes() const
      {
        return _index_set.bytes();
      }
    };
 
    /* ***************************************************************************************** */
    /* ***************************************************************************************** */
    /* ***************************************************************************************** */
 
    template<typename Shape_>
    class IndexSetHolder :
      public IndexSetHolder<typename Shape::FaceTraits<Shape_, Shape_::dimension - 1>::ShapeType>
    {
    public:
      typedef IndexSetHolder<typename Shape::FaceTraits<Shape_, Shape_::dimension - 1>::ShapeType> BaseClass;
      typedef IndexSetWrapper<Shape_> IndexSetWrapperType;
 
      template<
        int cell_dim_,
        int face_dim_>
      struct IndexSet
      {
        typedef Geometry::IndexSet
        <
          Shape::FaceTraits
          <
            typename Shape::FaceTraits<Shape_, cell_dim_> ::ShapeType,
            face_dim_
          >
          ::count
        > Type;
      }; // struct IndexSet<...>
 
    protected:
      IndexSetWrapperType _index_set_wrapper;
 
    public:
      IndexSetHolder() :
        BaseClass(),
        _index_set_wrapper()
      {
      }
 
      explicit IndexSetHolder(const Index num_entities[]) :
        BaseClass(num_entities),
        _index_set_wrapper(num_entities)
      {
      }
 
      template<std::size_t m_>
      explicit IndexSetHolder(const std::array<Index, m_>& num_entities) :
        BaseClass(num_entities),
        _index_set_wrapper(num_entities)
      {
        static_assert(int(m_) > Shape_::dimension, "invalid array size");
      }
 
      IndexSetHolder(IndexSetHolder&& other) :
        BaseClass(std::forward<BaseClass>(other)),
        _index_set_wrapper(std::forward<IndexSetWrapperType>(other._index_set_wrapper))
      {
      }
 
      IndexSetHolder& operator=(IndexSetHolder&& other)
      {
        if(this == &other)
          return *this;
 
        BaseClass::operator=(std::forward<BaseClass>(other));
        _index_set_wrapper = std::forward<IndexSetWrapperType>(other._index_set_wrapper);
 
        return *this;
      }
 
      virtual ~IndexSetHolder()
      {
      }
 
      void clone(const IndexSetHolder& other)
      {
        BaseClass::clone(other);
        this->_index_set_wrapper.clone(other._index_set_wrapper);
      }
 
      IndexSetHolder clone() const
      {
        IndexSetHolder ish;
        ish.clone(*this);
        return ish;
      }
 
      template<int shape_dim_>
      IndexSetWrapper<typename Shape::FaceTraits<Shape_, shape_dim_>::ShapeType>& 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 IndexSetHolder<ShapeType>::_index_set_wrapper;
      }
 
      template<int shape_dim_>
      const IndexSetWrapper<typename Shape::FaceTraits<Shape_, shape_dim_>::ShapeType>& 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 IndexSetHolder<ShapeType>::_index_set_wrapper;
      }
 
      template<
        int cell_dim_,
        int face_dim_>
      Geometry::IndexSet<
        Shape::FaceTraits<
          typename Shape::FaceTraits<
            Shape_,
            cell_dim_>
          ::ShapeType,
          face_dim_>
        ::count>&
      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 Geometry::IndexSet<
        Shape::FaceTraits<
          typename Shape::FaceTraits<
            Shape_,
            cell_dim_>
          ::ShapeType,
          face_dim_>
        ::count>&
      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 "IndexSetHolder<" + Shape_::name() + ">";
      }
 
      std::size_t bytes() const
      {
        return BaseClass::bytes() + _index_set_wrapper.bytes();
      }
    };
 
    template<>
    class IndexSetHolder<Shape::Vertex>
    {
    public:
      IndexSetHolder()
      {
      }
 
      explicit IndexSetHolder(const Index* /*num_entities*/)
      {
      }
 
      template<std::size_t m_>
      explicit IndexSetHolder(const std::array<Index, m_>&)
      {
      }
 
      IndexSetHolder(IndexSetHolder&&)
      {
      }
 
      IndexSetHolder& operator=(IndexSetHolder&&)
      {
        return *this;
      }
 
      virtual ~IndexSetHolder()
      {
      }
 
      void clone(const IndexSetHolder&)
      {
      }
 
      IndexSetHolder clone() const
      {
        return IndexSetHolder();
      }
 
      template<std::size_t np_>
      void permute(const std::array<Adjacency::Permutation, np_>&,
        const std::array<Adjacency::Permutation, np_>&)
      {
      }
 
      static String name()
      {
        return "IndexSetHolder<Vertex>";
      }
 
      std::size_t bytes() const
      {
        return std::size_t(0);
      }
    };
 
    template<int shape_dim_, int co_dim_>
    struct IndexSetHolderDimensionUpdater
    {
      template<typename IndexSetHolderType_>
      static void update(IndexSetHolderType_& DOXY(ish))
      {
        //dummy
      }
    };
 
    template<int shape_dim_>
    struct IndexSetHolderDimensionUpdater<shape_dim_,1>
    {
      template<typename IndexSetHolderType_>
      static void update(IndexSetHolderType_& ish)
      {
        auto& vert_at_subshape_index_set(ish.template get_index_set<shape_dim_-1,0>());
 
        Index num_shapes(ish.template get_index_set<shape_dim_,0>().get_num_entities());
        // Get the type of the IndexSet to dimensions <Shape_::dimension, face_dim_>
        typename std::remove_reference<decltype( (ish.template get_index_set<shape_dim_, shape_dim_-1>()) )>::type
          // The output IndexSet has num_shapes entities and the maximum index is the number of subshapes
          subshape_at_shape_index_set(num_shapes, vert_at_subshape_index_set.get_num_entities());
 
        // Replace the corresponding IndexSet in the IndexSetHolder by the just generated IndexSet of the right
        // size
        ish.template get_index_set<shape_dim_, shape_dim_-1>() = std::move(subshape_at_shape_index_set);
      }
    };
 
    template<int shape_dim_>
    struct IndexSetHolderDimensionUpdater<shape_dim_,2>
    {
      template<typename IndexSetHolderType_>
      static void update(IndexSetHolderType_& ish)
      {
        // Update edge\@cell IndexSet dimensions
        // vert\@edge IndexSet
        auto& vert_at_subshape_index_set(ish.template get_index_set<shape_dim_-2,0>());
        // Number of cells
        Index num_shapes(ish.template get_index_set<shape_dim_,0>().get_num_entities());
        // Get the type of the IndexSet to dimensions <Shape_::dimension, Shape::dimension-2>
        typename std::remove_reference<decltype( (ish.template get_index_set<shape_dim_, shape_dim_-2>()) )>::type
          // The output IndexSet has num_shapes entities and the maximum index is the number of subshapes (=edges)
          subshape_at_shape_index_set(num_shapes, vert_at_subshape_index_set.get_num_entities());
        // Replace the corresponding IndexSet in the IndexSetHolder by the just generated IndexSet of the right
        // size
        ish.template get_index_set<shape_dim_, shape_dim_-2>() = std::move(subshape_at_shape_index_set);
 
        // Now do it again for the edge\@face IndexSet
        // Number of faces
        Index num_intermediate_shapes(ish.template get_index_set<shape_dim_-1,0>().get_num_entities());
        // Get the type of the IndexSet to dimensions <Shape_::dimension-1, Shape::dimension-2>
        typename std::remove_reference<decltype( (ish.template get_index_set<shape_dim_-1, shape_dim_-2>()) )>::type
          // The output IndexSet has num_intermediate_shapes entities and the maximum index is the number of subshapes
          subshape_at_intermediate_index_set(
            num_intermediate_shapes, vert_at_subshape_index_set.get_num_entities());
        // Replace the corresponding IndexSet in the IndexSetHolder by the just generated IndexSet of the right
        // size
        ish.template get_index_set<shape_dim_-1, shape_dim_-2>() = std::move(subshape_at_intermediate_index_set);
      }
    };
 
    template<int dim_>
    struct NumEntitiesExtractor
    {
      template<typename IndexSetHolderType_>
      static void set_num_entities(const IndexSetHolderType_& ish, Index* num_entities)
      {
        // The number of entities of dimension dim_ is the length of the vertex\@shape[dim_] IndexSet
        num_entities[dim_] = ish.template get_index_set<dim_,0>().get_num_entities();
        // Recurse down
        NumEntitiesExtractor<dim_-1>::set_num_entities(ish, num_entities);
      }
 
      template<typename IndexSetHolderType_>
      static void set_num_entities_with_numverts(const IndexSetHolderType_& ish, Index* num_entities)
      {
        // The number of entities of dimension dim_ is the length of the vertex\@shape[dim_] IndexSet
        num_entities[dim_] = ish.template get_index_set<dim_,0>().get_num_entities();
        // Recurse down
        NumEntitiesExtractor<dim_-1>::set_num_entities_with_numverts(ish, num_entities);
      }
    };
 
    template<>
    struct NumEntitiesExtractor<1>
    {
      template<typename IndexSetHolderType_>
      static void set_num_entities(const IndexSetHolderType_& ish, Index* num_entities)
      {
        num_entities[1] = ish.template get_index_set<1,0>().get_num_entities();
      }
 
      template<typename IndexSetHolderType_>
      static void set_num_entities_with_numverts(const IndexSetHolderType_& ish, Index* num_entities)
      {
        num_entities[0] = ish.template get_index_set<1,0>().get_index_bound();
        num_entities[1] = ish.template get_index_set<1,0>().get_num_entities();
      }
    };
  } // namespace Geometry
} // namespace FEAT