utility.hh

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=8 sw=2 sts=2:
 
#ifndef DUNE_PDELAB_ORDERING_UTILITY_HH
#define DUNE_PDELAB_ORDERING_UTILITY_HH
 
#include <vector>
#include <bitset>
 
#include <dune/pdelab/common/dofindex.hh>
#include <dune/pdelab/common/globaldofindex.hh>
#include <dune/pdelab/ordering/transformations.hh>
 
namespace Dune {
  namespace PDELab {
 
 
    struct MergeMode
    {
 
      enum type {
        lexicographic, 
        interleaved  
      };
 
    };
 
    enum class MultiIndexOrder {
      Inner2Outer,
      Outer2Inner,
    };
 
 
#ifndef DOXYGEN
 
    namespace ordering {
 
      // This is an implementation detail of the composite orderings, no need to confuse our users!
      struct update_direct_children
        : public TypeTree::DirectChildrenVisitor
        , public TypeTree::DynamicTraversal
      {
 
        template<typename GFS, typename Child, typename TreePath, typename ChildIndex>
        void afterChild(const GFS& gfs, Child& child, TreePath tp, ChildIndex childIndex) const
        {
          child.update();
        }
 
      };
 
    } // end namespace ordering
 
#endif // DOXYGEN
 
 
    struct DefaultDOFIndexAccessor
    {
 
      template<typename DOFIndex, typename SizeType, typename IndexType>
      static typename std::enable_if<
        std::is_integral<IndexType>::value
        >::type
      store(DOFIndex& dof_index, const GeometryType& gt, SizeType entity_index, IndexType tree_index)
      {
        dof_index.clear();
        dof_index.entityIndex()[0] = GlobalGeometryTypeIndex::index(gt);
        dof_index.entityIndex()[1] = entity_index;
        dof_index.treeIndex().push_back(tree_index);
      }
 
      template<typename DOFIndex, typename SizeType, typename IndexType>
      static typename std::enable_if<
        !std::is_integral<IndexType>::value
        >::type
      store(DOFIndex& dof_index, const GeometryType& gt, SizeType entity_index, IndexType tree_index)
      {
        dof_index.entityIndex()[0] = GlobalGeometryTypeIndex::index(gt);
        dof_index.entityIndex()[1] = entity_index;
        dof_index.treeIndex() = tree_index;
      }
 
      template<typename DOFIndex, typename SizeType, typename IndexType>
      static typename std::enable_if<
        std::is_integral<IndexType>::value
        >::type
      store(DOFIndex& dof_index, SizeType gt_index, SizeType entity_index, IndexType tree_index)
      {
        dof_index.clear();
        dof_index.entityIndex()[0] = gt_index;
        dof_index.entityIndex()[1] = entity_index;
        dof_index.treeIndex().push_back(tree_index);
      }
 
      template<typename DOFIndex, typename SizeType, typename IndexType>
      static typename std::enable_if<
        !std::is_integral<IndexType>::value
        >::type
      store(DOFIndex& dof_index, SizeType gt_index, SizeType entity_index, IndexType tree_index)
      {
        dof_index.entityIndex()[0] = gt_index;
        dof_index.entityIndex()[1] = entity_index;
        dof_index.treeIndex() = tree_index;
      }
 
 
      struct GeometryIndex
      {
 
        template<typename Index>
        static std::size_t geometryType(const Index& geometry_index)
        {
          return geometry_index[0];
        }
 
        template<typename Index>
        static std::size_t entityIndex(const Index& geometry_index)
        {
          return geometry_index[1];
        }
 
        template<typename Index, typename SizeType>
        static void store(Index& index, const GeometryType& gt, SizeType entity_index)
        {
          index[0] = GlobalGeometryTypeIndex::index(gt);
          index[1] = entity_index;
        }
 
        template<typename Index, typename SizeType>
        static void store(Index& index, SizeType geometry_index, SizeType entity_index)
        {
          index[0] = geometry_index;
          index[1] = entity_index;
        }
 
      };
 
      template<typename DOFIndex>
      static std::size_t geometryType(const DOFIndex& dof_index)
      {
        return GeometryIndex::geometryType(dof_index.entityIndex());
      }
 
      template<typename DOFIndex>
      static std::size_t entityIndex(const DOFIndex& dof_index)
      {
        return GeometryIndex::entityIndex(dof_index.entityIndex());
      }
 
    };
 
    struct SimpleDOFIndexAccessor
    {
 
      template<typename DOFIndex, typename SizeType>
      static void store(DOFIndex& dof_index, const GeometryType& gt, SizeType entity_index, SizeType tree_index)
      {
        dof_index = entity_index;
      }
 
    };
 
 
    template<typename DI, typename CI, MultiIndexOrder CIOrder = MultiIndexOrder::Inner2Outer>
    struct SimpleOrderingTraits
    {
 
      typedef DI DOFIndex;
 
      typedef CI ContainerIndex;
 
      typedef std::size_t SizeType;
 
      typedef DefaultDOFIndexAccessor DOFIndexAccessor;
 
      static constexpr MultiIndexOrder ContainerIndexOrder = CIOrder;
    };
 
 
    template<typename SizeType_, typename CI, MultiIndexOrder CIOrder>
    struct SimpleOrderingTraits<SimpleDOFIndex<SizeType_>,CI,CIOrder>
    {
 
      typedef SimpleDOFIndex<SizeType_> DOFIndex;
 
      typedef CI ContainerIndex;
 
      typedef SizeType_ SizeType;
 
      typedef SimpleDOFIndexAccessor DOFIndexAccessor;
 
      static constexpr MultiIndexOrder ContainerIndexOrder = CIOrder;
    };
 
    template <typename DI, typename CI,
              MultiIndexOrder CIOrder = MultiIndexOrder::Inner2Outer>
    struct OrderingTraits : public SimpleOrderingTraits<DI, CI, CIOrder> {
 
      // The maximum dimension supported (length of bitsets)
      // 32 dimensions should probably be fine for now... ;-)
      static const std::size_t max_dim = 32;
 
      typedef std::bitset<max_dim> CodimFlag;
 
      typedef typename DI::TreeIndex TreeIndex;
 
      typedef typename DI::View DOFIndexView;
      typedef typename DI::View::TreeIndex TreeIndexView;
 
      typedef typename DI::size_type SizeType;
      typedef typename DI::size_type size_type;
    };
 
    template <typename ES, typename DI, typename CI,
              MultiIndexOrder CIOrder = MultiIndexOrder::Outer2Inner>
    struct LocalOrderingTraits : public OrderingTraits<DI, CI, CIOrder> {
 
      using EntitySet = ES;
      using GridView = typename ES::GridView;
    };
 
    template<typename ES, typename DI, typename CI>
    struct GridViewOrderingTraits
      : public LocalOrderingTraits<ES,DI,CI>
    {
 
      typedef typename DI::EntityIndex EntityIndex;
      typedef typename DI::View::EntityIndex EntityIndexView;
 
    };
 
 
    template<typename DI, typename CI>
    class VirtualOrderingBase
    {
    public:
 
      typedef OrderingTraits<DI,CI> Traits;
 
      VirtualOrderingBase() {}
      virtual ~VirtualOrderingBase() = default;
 
      virtual void map_index_dynamic(typename Traits::DOFIndexView di, typename Traits::ContainerIndex& ci) const = 0;
    };
 
 
    template<typename child_type>
    struct extract_child_bases
      : public TypeTree::DirectChildrenVisitor
      , public TypeTree::DynamicTraversal
    {
 
      template<typename Node, typename Child, typename TreePath, typename ChildIndex>
      void afterChild(const Node& node, Child& child, TreePath tp, ChildIndex child_index)
      {
        extract_child(node,child,child_index);
      }
 
      template<typename Node, typename Child, typename ChildIndex>
      typename std::enable_if<Node::has_dynamic_ordering_children>::type
      extract_child(const Node& node, Child& child, ChildIndex child_index)
      {
        _children[child_index] = &child;
      }
 
      template<typename Node, typename Child, typename ChildIndex>
      typename std::enable_if<!Node::has_dynamic_ordering_children>::type
      extract_child(const Node& node, Child& child, ChildIndex child_index)
      {
      }
 
      extract_child_bases(std::vector<child_type*>& children)
        : _children(children)
      {}
 
    private:
      std::vector<child_type*>& _children;
 
    };
 
 
 
    struct DummyDOFIndexIterator
    {
 
      typedef std::size_t size_type;
 
      DummyDOFIndexIterator& operator++()
      {
        return *this;
      }
 
      DummyDOFIndexIterator& operator+=(size_type i)
      {
        return *this;
      }
 
      DummyDOFIndexIterator& operator*()
      {
        return *this;
      }
 
      DummyDOFIndexIterator* operator->()
      {
        return this;
      }
 
      DummyDOFIndexIterator& treeIndex()
      {
        return *this;
      }
 
      bool operator==(const DummyDOFIndexIterator& r) const
      {
        return true;
      }
 
      bool operator!=(const DummyDOFIndexIterator& r) const
      {
        return !operator==(r);
      }
 
      void push_back(size_type i)
      {}
 
    };
 
    template <class Size, class ContainerIndex_, MultiIndexOrder OriginOrder>
    struct SizeProviderAdapter {
 
      template <class Ordering>
      SizeProviderAdapter(const std::shared_ptr<const Ordering> &ordering)
          : _size_provider([=](const ContainerIndex_ &partial_multiindex) {
              return ordering->size(partial_multiindex);
            }) {
        static_assert(Ordering::Traits::ContainerIndexOrder == OriginOrder);
      }
 
      using ContainerIndex = ContainerIndex_;
 
      using SizePrefix = ContainerIndex;
 
      using SizeType = Size;
 
      static constexpr MultiIndexOrder ContainerIndexOrder = MultiIndexOrder::Outer2Inner;
 
      SizeType size(const SizePrefix &prefix) const {
        if constexpr (OriginOrder == MultiIndexOrder::Inner2Outer) {
          // reversing Outer2Inner prefix into a Inner2Outer suffix
          ContainerIndex suffix;
          suffix.resize(prefix.size());
          std::reverse_copy(prefix.begin(), prefix.end(),suffix.begin());
          // forward size request to ordering with new Inner2Outer suffix
          return _size_provider(suffix);
        } else {
          // prefix is already Outer2Inner, forward to size provider
          return _size_provider(prefix);
        }
      }
 
    private:
 
      const std::function<SizeType(const ContainerIndex &)> _size_provider;
    };
 
    template<class Ordering>
    SizeProviderAdapter(const std::shared_ptr<const Ordering>& ordering)
        -> SizeProviderAdapter<typename Ordering::Traits::SizeType,
                              typename Ordering::Traits::ContainerIndex,
                              Ordering::Traits::ContainerIndexOrder>;
 
 
  } // namespace PDELab
} // namespace Dune
 
#endif // DUNE_PDELAB_ORDERING_UTILITY_HH