Added computation of global DOF IDs for pre-bases

src/amdis/common/TupleUtility.hpp

@@ -32,4 +32,28 @@ namespace AMDiS
 
   } // end namespace Impl
 
+  template <class Tuple, template <class> class Map>
+  struct MapTuple;
+
+  template <class Tuple, template <class> class Map>
+  using MapTuple_t = typename MapTuple<Tuple,Map>::type;
+
+  template <class... T, template <class> class Map>
+  struct MapTuple<std::tuple<T...>, Map>
+  {
+    using type = std::tuple<Map<T>...>;
+  };
+
+  template <class Indices, template <std::size_t> class Map>
+  struct IndexMapTuple;
+
+  template <class Indices, template <std::size_t> class Map>
+  using IndexMapTuple_t = typename IndexMapTuple<Indices,Map>::type;
+
+  template <std::size_t... I, template <std::size_t> class Map>
+  struct IndexMapTuple<std::index_sequence<I...>, Map>
+  {
+    using type = std::tuple<Map<I>...>;
+  };
+
 } // end namespace AMDiS

src/amdis/functions/CMakeLists.txt

 #install headers
 
 install(FILES
+    GlobalIdSet.hpp
     HierarchicNodeToRangeMap.hpp
     Interpolate.hpp
+    Nodes.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/functions)

src/amdis/functions/GlobalIdSet.hpp

+#pragma once
+
+#include <tuple>
+#include <utility>
+#include <vector>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/version.hh>
+#include <dune/functions/functionspacebases/nodes.hh>
+#include <dune/grid/common/gridenums.hh>
+#include <dune/localfunctions/common/localkey.hh>
+
+#include <amdis/common/Apply.hpp>
+#include <amdis/common/ForEach.hpp>
+#include <amdis/common/TupleUtility.hpp>
+#include <amdis/functions/Nodes.hpp>
+#include <amdis/Output.hpp>
+
+namespace Dune
+{
+  namespace Functions
+  {
+    // forward declarations...
+    template <class RB, class TP>
+    class SubspaceBasis;
+
+    template <class MI, class IMS, class... SPB>
+    class CompositePreBasis;
+
+    template <class GV, int k, class MI>
+    class LagrangePreBasis;
+
+    template <class MI, class IMS, class SPB, std::size_t C>
+    class PowerPreBasis;
+
+    template <class GV, int k, class MI>
+    class LagrangeDGPreBasis;
+
+    template<typename GV, class MI, bool HI>
+    class TaylorHoodPreBasis;
+  }
+}
+
+
+namespace AMDiS
+{
+  // forward declaration
+  template <class PreBasis, class TP>
+  class NodeIdSet;
+
+  /// \brief PRovide global ids for all DOFs in a global basis
+  /**
+   * A GlobalBasisIdSet provides and interface to retrieve unique IDs
+   * (over all partitions) for all DOFs in the basis. It is implemented
+   * in terms of the \ref NodeIdSet for each basis node that locally
+   * on each element calculates a global unique id.
+   *
+   * **Examples:**
+   * ```
+    GlobalBasisIdSet<Basis> idSet(basis);
+    for (const auto& e : elements(basis.gridView())) {
+      idSet.bind(e);
+      for (std::size_t i = 0; i < idSet.size(); ++i) {
+        auto id = idSet.id(i); // global unique ID of i'th local DOF
+        auto pt = idSet.partitionType(i); // partition type of i'th local DOF
+      }
+      idSet.unbind();
+    }
+   * ```
+   *
+   **/
+  template <class GB>
+  class GlobalBasisIdSet
+  {
+  public:
+    using GlobalBasis = GB;
+    using Tree = typename GB::LocalView::Tree;
+    using size_type = std::size_t;
+
+    using GridView = typename GlobalBasis::GridView;
+    using Grid = typename GridView::Grid;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using EntityIdType = typename Grid::GlobalIdSet::IdType; 
+    using IdType = std::pair<EntityIdType, std::size_t>;
+    using PartitionType = Dune::PartitionType;
+
+    using PreBasis = typename GlobalBasis::PreBasis;
+    using TreePath = typename GlobalBasis::PrefixPath;
+    using NodeIdSet = AMDiS::NodeIdSet<PreBasis, TreePath>;
+
+  public:
+    GlobalBasisIdSet(GlobalBasis const& globalBasis)
+      : tree_(makeNode(globalBasis.preBasis(), TreePath{}))
+      , nodeIdSet_(globalBasis.gridView())
+    {
+      Dune::Functions::initializeTree(tree_);
+    }
+
+    /// \brief Bind the IdSet to a grid element.
+    /**
+     * Binding the IdSet to an element collects the ids from all DOFs
+     * on that element and stores it in a local vector. Thus, the global
+     * DOF ids can be extracted using \ref id() afterwards cheaply.
+     **/
+    void bind(Element const& element)
+    {
+      Dune::Functions::bindTree(tree_, element);
+      nodeIdSet_.bind(tree_);
+      data_.resize(size());
+      nodeIdSet_.fillIn(data_.begin());
+    }
+
+    /// \brief unbind from the element
+    void unbind()
+    {
+      nodeIdSet_.unbind();
+    }
+
+    /// \brief The number of DOFs on the current element in the basis tree
+    size_type size() const
+    {
+      return tree_.size();
+    }
+
+    /// \brief Return the global unique ID of the `i`th DOF on the
+    /// currently bound element in the basis tree.
+    IdType id(size_type i) const
+    {
+      assert( i < data_.size() );
+      return data_[i].first;
+    }
+
+    /// \brief Return the global unique ID of the entity associated to the
+    /// `i`th DOF on the currently bound element in the basis tree.
+    EntityIdType entityId(size_type i) const
+    {
+      return id(i).first;
+    }
+
+    /// \brief Return the partition type of the `i`th DOF on the
+    /// currently bound element in the basis tree.
+    PartitionType partitionType(size_type i) const
+    {
+      assert( i < data_.size() );
+      return data_[i].second;
+    }
+
+  protected:
+    Tree tree_;
+    NodeIdSet nodeIdSet_;
+    using Data = std::pair<IdType, PartitionType>;
+    std::vector<Data> data_;
+  };
+
+
+  template <class RB, class TP>
+  class GlobalBasisIdSet<Dune::Functions::SubspaceBasis<RB,TP>>
+      : public GlobalBasisIdSet<RB>
+  {
+  public:
+    GlobalBasisIdSet(Dune::Functions::SubspaceBasis<RB,TP> const& basis)
+      : GlobalBasisIdSet<RB>(basis.rootBasis())
+    {}
+  };
+
+
+  template <class PB, class TP>
+  class NodeIdSet
+  {
+  public:
+    using PreBasis = PB;
+    using Node = Node_t<PreBasis,TP>;
+    using GridView = typename PreBasis::GridView;
+    using size_type = std::size_t;
+
+    static_assert(Node::isLeaf, "Generic NodeIdSet implemented for LeafNodes only. Provide a spcialization for your node!");
+
+  private:
+    static constexpr int dim = GridView::template Codim<0>::Entity::mydimension;
+
+  public:
+    NodeIdSet(GridView const& gridView)
+      : gridView_(gridView)
+    {}
+
+    /// \brief Bind the idset to a tree node
+    void bind(const Node& node)
+    {
+      node_ = &node;
+      size_ = node_->finiteElement().size();
+    }
+
+    /// \brief Unbind the idset
+    void unbind()
+    {
+      node_ = nullptr;
+    }
+
+    /// \brief Number of DOFs in this node
+    size_type size() const
+    {
+      return size_;
+    }
+
+    /// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis
+    // [[expects: node_ != nullptr]]
+    template <class It>
+    It fillIn(It it, size_type shift = 0) const
+    {
+      const auto& gridIdSet = gridView_.grid().globalIdSet();
+
+      for (size_type i = 0; i < size_ ; ++i, ++it) {
+        Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
+        unsigned int s = localKey.subEntity();
+        unsigned int c = localKey.codim();
+        unsigned int idx = localKey.index();
+        if (!(c == GridView::dimension || c == 0 || idx == 0))
+          DUNE_THROW(Dune::NotImplemented, "Bases with more then one DoF per subentity are not supported.");
+
+        it->first = {gridIdSet.subId(node_->element(), s, c), shift + idx};
+
+        it->second = Dune::Hybrid::switchCases(std::make_index_sequence<dim+1>{}, c,
+          [&](auto codim) { return node_->element().template subEntity<codim>(s).partitionType(); },
+          [&]() {
+            error_exit("Invalid codimension c = {}", c);
+            return Dune::PartitionType{};
+          });
+      }
+
+      return it;
+    }
+
+  protected:
+    GridView gridView_;
+    const Node* node_ = nullptr;
+    size_type size_ = 0;
+  };
+
+
+  template <class MI, class IMS, class SPB, std::size_t C, class TP>
+  class NodeIdSet<Dune::Functions::PowerPreBasis<MI,IMS,SPB,C>, TP>
+  {
+  public:
+    using PreBasis = Dune::Functions::PowerPreBasis<MI,IMS,SPB,C>;
+    using Node = Node_t<PreBasis,TP>;
+    using GridView = typename PreBasis::GridView;
+    using size_type = std::size_t;
+
+  protected:
+    using SubPreBasis = typename PreBasis::SubPreBasis;
+    using SubTreePath = decltype(Dune::TypeTree::push_back(std::declval<TP>(), std::size_t(0)));
+    using SubNodeIdSet = NodeIdSet<SubPreBasis, SubTreePath>;
+    static const std::size_t children = C;
+
+  public:
+    NodeIdSet(GridView const& gridView)
+      : subIds_(gridView)
+    {}
+
+    /// \brief Bind the idset to a tree node
+    void bind(const Node& node)
+    {
+      node_ = &node;
+      subIds_.bind(node.child(0));
+    }
+
+    /// \brief Unbind the idset
+    void unbind()
+    {
+      node_ = nullptr;
+      subIds_.unbind();
+    }
+
+    /// \brief Number of DOFs in this node
+    size_type size() const
+    {
+      return node_->size();
+    }
+
+    /// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis
+    // [[expects: node_ != nullptr]]
+    template <class It>
+    It fillIn(It it, size_type shift = 0) const
+    {
+      for (std::size_t child = 0; child < children; ++child)
+      {
+        size_type subTreeSize = subIds_.size();
+        it = subIds_.fillIn(it, shift);
+        shift += subTreeSize;
+      }
+      return it;
+    }
+
+  protected:
+    SubNodeIdSet subIds_;
+    const Node* node_ = nullptr;
+  };
+
+
+  template <class MI, class IMS, class... SPB, class TP>
+  class NodeIdSet<Dune::Functions::CompositePreBasis<MI,IMS,SPB...>, TP>
+  {
+  public:
+    using PreBasis = Dune::Functions::CompositePreBasis<MI,IMS,SPB...>;
+    using Node = Node_t<PreBasis,TP>;
+    using GridView = typename PreBasis::GridView;
+    using size_type = std::size_t;
+
+  protected:
+    static const std::size_t children = sizeof...(SPB);
+    using ChildIndices = std::make_index_sequence<children>;
+
+    // The I'th SubPreBasis
+    template <std::size_t I>
+    using SubPreBasis = typename PreBasis::template SubPreBasis<I>;
+
+    template <std::size_t I>
+    using SubTreePath = decltype(Dune::TypeTree::push_back(std::declval<TP>(), Dune::index_constant<I>{}));
+
+    template <std::size_t I>
+    using SubNodeIdSet = NodeIdSet<SubPreBasis<I>, SubTreePath<I>>;
+
+    // A tuple of NodeIdSets for the SubPreBases
+    using IdsTuple = IndexMapTuple_t<std::make_index_sequence<children>, SubNodeIdSet>;
+
+  public:
+    NodeIdSet(GridView const& gridView)
+      : idsTuple_(Tools::apply_indices<children>([&](auto... i) {
+          return std::make_tuple(SubNodeIdSet<VALUE(i)>(gridView)...);
+          }))
+    {}
+
+    /// \brief Bind the idset to a tree node
+    void bind(const Node& node)
+    {
+      node_ = &node;
+      Tools::for_range<0,children>([&](auto i) {
+        std::get<VALUE(i)>(idsTuple_).bind(node.child(i));
+      });
+    }
+
+    /// \brief Unbind the idset
+    void unbind()
+    {
+      node_ = nullptr;
+      Tools::for_each(idsTuple_, [](auto& ids) {
+        ids.unbind();
+      });
+    }
+
+    /// \brief Number of DOFs in this node
+    size_type size() const
+    {
+      return node_->size();
+    }
+
+    /// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis
+    // [[expects: node_ != nullptr]]
+    template <class It>
+    It fillIn(It it, size_type shift = 0) const
+    {
+      Tools::for_each(idsTuple_, [&](auto const& ids)
+      {
+        size_type subTreeSize = ids.size();
+        it = ids.fillIn(it, shift);
+        shift += subTreeSize;
+      });
+      return it;
+    }
+
+  private:
+    IdsTuple idsTuple_;
+    const Node* node_ = nullptr;
+  };
+
+
+  template <class GV, class MI, bool HI, class TP>
+  class NodeIdSet<Dune::Functions::TaylorHoodPreBasis<GV,MI,HI>, TP>
+  {
+  public:
+    using PreBasis = Dune::Functions::TaylorHoodPreBasis<GV,MI,HI>;
+    using Node = Node_t<PreBasis,TP>;
+    using GridView = typename PreBasis::GridView;
+    using size_type = std::size_t;
+
+  private:
+    using PTP = decltype(Dune::TypeTree::push_back(std::declval<TP>(), Dune::index_constant<1>{}));
+    using VTP = decltype(Dune::TypeTree::push_back(std::declval<TP>(), Dune::index_constant<0>{}));
+    using V0TP = decltype(Dune::TypeTree::push_back(std::declval<VTP>(), std::size_t(0)));
+
+    // Note: PreBasis::PQ1PreBasis is private
+    using PQ1PreBasis = typename NodeIndexSet_t<PreBasis,TP>::PQ1NodeIndexSet::PreBasis;
+    using PQ2PreBasis = typename NodeIndexSet_t<PreBasis,TP>::PQ2NodeIndexSet::PreBasis;
+
+    using PQ1NodeIdSet = NodeIdSet<PQ1PreBasis, PTP>; // pressure
+    using PQ2NodeIdSet = NodeIdSet<PQ2PreBasis, V0TP>; // velocity
+
+  private:
+    static constexpr int dow = GridView::dimensionworld;
+
+  public:
+    NodeIdSet(GridView const& gridView)
+      : pq1NodeIdSet_(gridView)
+      , pq2NodeIdSet_(gridView)
+    {}
+
+    /// \brief Bind the idset to a tree node
+    void bind(const Node& node)
+    {
+      node_ = &node;
+      using namespace Dune::Indices;
+      pq1NodeIdSet_.bind(node.child(_1));
+      pq2NodeIdSet_.bind(node.child(_0, 0));
+    }
+
+    /// \brief Unbind the idset
+    void unbind()
+    {
+      pq1NodeIdSet_.unbind();
+      pq2NodeIdSet_.unbind();
+      node_ = nullptr;
+    }
+
+    /// \brief Number of DOFs in this node
+    size_type size() const
+    {
+      return node_->size();
+    }
+
+    /// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis
+    // [[expects: node_ != nullptr]]
+    template <class It>
+    It fillIn(It it, size_type shift = 0) const
+    {
+      for (int child = 0; child < dow; ++child) {
+        size_type subTreeSize = pq2NodeIdSet_.size();
+        it = pq2NodeIdSet_.fillIn(it, shift);
+        shift += subTreeSize;
+      }
+      it = pq1NodeIdSet_.fillIn(it, shift);
+      return it;
+    }
+
+  protected:
+    PQ1NodeIdSet pq1NodeIdSet_;
+    PQ2NodeIdSet pq2NodeIdSet_;
+    const Node* node_ = nullptr;
+  };
+
+
+  template <class GV, int k, class MI, class TP>
+  class NodeIdSet<Dune::Functions::LagrangeDGPreBasis<GV, k, MI>, TP>
+  {
+  public:
+    using PreBasis = Dune::Functions::LagrangeDGPreBasis<GV, k, MI>;
+    using Node = Node_t<PreBasis, TP>;
+    using GridView = typename PreBasis::GridView;
+    using size_type = std::size_t;
+
+    NodeIdSet(GridView const& gridView)
+      : gridView_(gridView)
+    {}
+
+    /// \brief Bind the idset to a tree node
+    void bind(const Node& node)
+    {
+      node_ = &node;
+      size_ = node_->finiteElement().size();
+    }
+
+    /// \brief Unbind the idset
+    void unbind()
+    {
+      node_ = nullptr;
+    }
+
+    /// \brief Number of DOFs in this node
+    size_type size() const
+    {
+      return size_;
+    }
+
+    /// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis
+    // [[expects: node_ != nullptr]]
+    template <class It>
+    It fillIn(It it, size_type shift = 0) const
+    {
+      const auto& gridIdSet = gridView_.grid().globalIdSet();
+      auto elementId = gridIdSet.id(node_->element());
+
+      for (size_type i = 0; i < size_; ++i, ++it)
+      {
+        it->first = {elementId, shift + i};
+        it->second = node_->element().partitionType();
+      }
+
+      return it;
+    }
+
+  protected:
+    GridView gridView_;
+    const Node* node_ = nullptr;
+    size_type size_ = 0;
+  };
+
+
+} // end namespace AMDiS

src/amdis/functions/Nodes.hpp

+#pragma once
+
+#include <dune/common/version.hh>
+
+namespace AMDiS
+{
+  // dune version independent extraction of node type from preBasis
+  template <class PB, class TP>
+  using Node_t =
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+    typename PB::template Node<TP>;
+#else
+    typename PB::Node;
+#endif
+
+  // dune version independent extraction of node indexset type from preBasis
+  template <class PB, class TP>
+  using NodeIndexSet_t =
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+    typename PB::template IndexSet<TP>;
+#else
+    typename PB::IndexSet;
+#endif
+
+  // dune version independent creation of node from preBasis
+  template <class PB, class TP>
+  auto makeNode(PB const& preBasis, TP const& treePath)
+  {
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+    return preBasis.node(treePath);
+#else
+    DUNE_UNUSED_PARAMETER(treePath);
+    return preBasis.makeNode();
+#endif
+  }
+
+} // end namespace AMDiS

test/CMakeLists.txt

@@ -46,6 +46,12 @@ dune_add_test(SOURCES FiniteElementTypeTest.cpp
 dune_add_test(SOURCES FilesystemTest.cpp
   LINK_LIBRARIES amdis)
 
+dune_add_test(SOURCES GlobalIdSetTest.cpp
+  LINK_LIBRARIES amdis
+  MPI_RANKS 2
+  TIMEOUT 300
+  CMAKE_GUARD MPI_FOUND)
+
 dune_add_test(SOURCES GradientTest.cpp
   LINK_LIBRARIES amdis)
 

test/GlobalIdSetTest.cpp

+#include <amdis/AMDiS.hpp>
+#include <amdis/ProblemStatTraits.hpp>
+#include <amdis/functions/GlobalIdSet.hpp>
+
+#include <dune/functions/functionspacebases/lagrangedgbasis.hh>
+#include <dune/functions/functionspacebases/taylorhoodbasis.hh>
+
+#include "Tests.hpp"
+
+using namespace AMDiS;
+
+template <class Basis>
+void checkBasisIds(const Basis& basis)