datahandleprovider.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_PDELAB_GRIDFUNCTIONSPACE_DATAHANDLEPROVIDER_HH
#define DUNE_PDELAB_GRIDFUNCTIONSPACE_DATAHANDLEPROVIDER_HH
 
#include <vector>
#include <stack>
 
#include <dune/common/typetraits.hh>
#include <dune/common/reservedvector.hh>
#include <dune/typetree/visitor.hh>
 
#include <dune/pdelab/ordering/utility.hh>
 
namespace Dune {
  namespace PDELab {
 
    namespace {
 
      template<typename EntityIndex>
      struct get_size_for_entity
        : public TypeTree::TreeVisitor
        , public TypeTree::DynamicTraversal
      {
 
        template<typename Ordering, typename TreePath>
        void leaf(const Ordering& ordering, TreePath tp)
        {
          _size += ordering.size(_entity_index);
        }
 
        get_size_for_entity(const EntityIndex& entity_index)
          : _size(0)
          , _entity_index(entity_index)
        {}
 
        std::size_t size() const
        {
          return _size;
        }
 
      private:
 
        std::size_t _size;
        const EntityIndex& _entity_index;
 
      };
 
 
      template<typename EntityIndex, typename OffsetIterator>
      struct get_leaf_offsets_for_entity
        : public TypeTree::TreeVisitor
        , public TypeTree::DynamicTraversal
      {
 
        template<typename Ordering, typename TreePath>
        void leaf(const Ordering& ordering, TreePath tp)
        {
          *(++_oit) = ordering.size(_entity_index);
        }
 
        get_leaf_offsets_for_entity(const EntityIndex& entity_index, OffsetIterator oit)
          : _oit(oit)
          , _entity_index(entity_index)
        {}
 
        OffsetIterator offsetIterator() const
        {
          return _oit;
        }
 
      private:
 
        OffsetIterator _oit;
        const EntityIndex& _entity_index;
 
      };
 
 
      template<typename DOFIndex, typename ContainerIndex, std::size_t tree_depth, bool map_dof_indices = false>
      struct indices_for_entity
        : public TypeTree::TreeVisitor
        , public TypeTree::DynamicTraversal
      {
 
        typedef std::size_t size_type;
        typedef typename DOFIndex::EntityIndex EntityIndex;
        typedef typename std::vector<ContainerIndex>::iterator CIIterator;
        typedef typename std::conditional<
          map_dof_indices,
          typename std::vector<DOFIndex>::iterator,
          DummyDOFIndexIterator
          >::type DIIterator;
 
 
        template<typename Ordering, typename Child, typename TreePath, typename ChildIndex>
        void beforeChild(const Ordering& ordering, const Child& child, TreePath tp, ChildIndex childIndex)
        {
          _stack.push(std::make_pair(_ci_it,_di_it));
        }
 
        template<typename Ordering, typename TreePath>
        void leaf(const Ordering& ordering, TreePath tp)
        {
          static_assert(tp.size()>=1, "Cannot call method 'indices_for_entity::leaf' for empty tree path");
          size_type size = ordering.extract_entity_indices(_entity_index,
                                                           tp.back(),
                                                           _ci_it,
                                                           _ci_end,
                                                           _di_it);
 
          _ci_end += size;
          _ci_it = _ci_end;
          _di_end += size;
          _di_it = _di_end;
        }
 
        template<typename Ordering, typename Child, typename TreePath, typename ChildIndex>
        void afterChild(const Ordering& ordering, const Child& child, TreePath tp, ChildIndex childIndex)
        {
          // pop
          ordering.extract_entity_indices(_entity_index,
                                          childIndex,
                                          _stack.top().first,
                                          _ci_end);
 
          if (Ordering::consume_tree_index)
            for (DIIterator it = _stack.top().second;
                 it != _di_end;
                 ++it)
              it->treeIndex().push_back(childIndex);
 
          _stack.pop();
        }
 
 
        indices_for_entity(const EntityIndex& entity_index,
                           CIIterator ci_begin,
                           DIIterator di_begin = DIIterator())
          : _entity_index(entity_index)
          , _ci_it(ci_begin)
          , _ci_end(ci_begin)
          , _di_it(di_begin)
          , _di_end(di_begin)
        {}
 
 
        // Exposed for multidomain support
        CIIterator ci_end() const
        {
          return _ci_end;
        }
 
        // Exposed for multidomain support
        DIIterator di_end() const
        {
          return _di_end;
        }
 
      private:
 
        const EntityIndex& _entity_index;
        CIIterator _ci_it;
        CIIterator _ci_end;
        DIIterator _di_it;
        DIIterator _di_end;
 
        std::stack<
          std::pair<
            CIIterator,
            DIIterator
            >,
          ReservedVector<
            std::pair<
              CIIterator,
              DIIterator
              >,
            tree_depth
            >
          > _stack;
      };
 
    } // anonymous namespace
 
 
    template<typename GFS>
    class DataHandleProvider
    {
 
    public:
 
      typedef std::size_t size_type;
 
      //------------------------------
      // generic data handle interface
      //------------------------------
 
      bool dataHandleContains (int codim) const
      {
        return gfs().ordering().contains(codim);
      }
 
      bool dataHandleFixedSize (int codim) const
      {
        return gfs().ordering().fixedSize(codim);
      }
 
 
      constexpr bool sendLeafSizes() const
      {
        return false;
      }
 
      template<typename Entity>
      size_type dataHandleSize (const Entity& e) const
      {
        typedef typename GFS::Ordering Ordering;
 
        typedef typename Ordering::Traits::DOFIndex::EntityIndex EntityIndex;
        EntityIndex ei;
 
        Ordering::Traits::DOFIndexAccessor::GeometryIndex::store(
          ei,
          e.type(),
          gfs().gridView().indexSet().index(e)
        );
 
        get_size_for_entity<EntityIndex> get_size(ei);
        TypeTree::applyToTree(gfs().ordering(),get_size);
 
        return get_size.size();
      }
 
      template<typename V, typename EntityIndex>
      void setup_dof_indices(V& v, size_type n, const EntityIndex& ei, std::integral_constant<bool,true>) const
      {
        v.resize(n);
        for (typename V::iterator it = v.begin(),
               endit = v.end();
             it != endit;
             ++it)
          {
            it->treeIndex().clear();
            it->entityIndex() = ei;
          }
      }
 
      template<typename V, typename EntityIndex>
      void setup_dof_indices(V& v, size_type n, const EntityIndex& ei, std::integral_constant<bool,false>) const
      {}
 
      template<typename V>
      typename V::iterator dof_indices_begin(V& v, std::integral_constant<bool,true>) const
      {
        return v.begin();
      }
 
      template<typename V>
      DummyDOFIndexIterator dof_indices_begin(V& v, std::integral_constant<bool,false>) const
      {
        return DummyDOFIndexIterator();
      }
 
      template<typename Entity, typename ContainerIndex, typename DOFIndex, typename OffsetIterator, bool map_dof_indices>
      void dataHandleIndices (const Entity& e,
                              std::vector<ContainerIndex>& container_indices,
                              std::vector<DOFIndex>& dof_indices,
                              OffsetIterator oit,
                              std::integral_constant<bool,map_dof_indices> map_dof_indices_value
                              ) const
      {
        typedef typename GFS::Ordering Ordering;
 
        static_assert((std::is_same<ContainerIndex,typename Ordering::Traits::ContainerIndex>::value),
                      "dataHandleContainerIndices() called with invalid ContainerIndex type.");
 
        typedef typename Ordering::Traits::DOFIndex::EntityIndex EntityIndex;
        EntityIndex ei;
 
        Ordering::Traits::DOFIndexAccessor::GeometryIndex::store(
          ei,
          e.type(),
          gfs().entitySet().indexSet().index(e)
        );
 
        get_leaf_offsets_for_entity<EntityIndex,OffsetIterator> get_offsets(ei,oit);
        TypeTree::applyToTree(gfs().ordering(),get_offsets);
        OffsetIterator end_oit = oit + (TypeTree::TreeInfo<Ordering>::leafCount + 1);
 
        // convert sizes to offsets - last entry contains total size
        std::partial_sum(oit,end_oit,oit);
        size_type size = *(oit + TypeTree::TreeInfo<Ordering>::leafCount);
 
        container_indices.resize(size);
        // Clear index state
        for (typename std::vector<ContainerIndex>::iterator it = container_indices.begin(),
               endit = container_indices.end();
             it != endit;
             ++it)
          it->clear();
 
        setup_dof_indices(dof_indices,size,ei,map_dof_indices_value);
 
#ifndef NDEBUG
        // this checks that gfs has common entity set on leaf nodes (throw if that's not the case)
        // in case of multidomain gfs, indices_for_entity will fail
        // because DOFIndex will have wrong entity indices
        auto es_visitor = impl::common_entity_set<typename GFS::Traits::EntitySet>{};
        TypeTree::applyToTree(gfs(), es_visitor);
#endif
 
        indices_for_entity<
          DOFIndex,
          ContainerIndex,
          TypeTree::TreeInfo<Ordering>::depth,
          map_dof_indices
          > extract_indices(ei,container_indices.begin(),dof_indices_begin(dof_indices,map_dof_indices_value));
        TypeTree::applyToTree(gfs().ordering(),extract_indices);
 
      }
 
    protected:
 
      const GFS& gfs() const
      {
        return static_cast<const GFS&>(*this);
      }
 
    };
 
  } // namespace PDELab
} // namespace Dune
 
#endif // DUNE_PDELAB_GRIDFUNCTIONSPACE_DATAHANDLEPROVIDER_HH