added twist utilities to flip the order of edge DOFs

src/amdis/functions/GlobalIdSet.hpp

@@ -11,11 +11,12 @@
 #include <dune/grid/common/gridenums.hh>
 #include <dune/localfunctions/common/localkey.hh>
 
+#include <amdis/Output.hpp>
 #include <amdis/common/Apply.hpp>
 #include <amdis/common/ForEach.hpp>
 #include <amdis/common/TupleUtility.hpp>
 #include <amdis/functions/Nodes.hpp>
-#include <amdis/Output.hpp>
+#include <amdis/utility/Twist.hpp>
 
 namespace Dune
 {
@@ -88,9 +89,20 @@ namespace AMDiS
     {
       using Super = std::pair<EntityIdType, std::size_t>;
 
-      IdType(int i = 0) : Super() {};
+      IdType(std::size_t i = 0) : Super(EntityIdType{}, i) {};
       using Super::Super;
 
+      IdType& operator++()
+      {
+        ++this->second;
+        return *this;
+      }
+
+      IdType operator+(std::size_t shift) const
+      {
+        return IdType{this->first, this->second + shift};
+      }
+
       friend std::ostream& operator<<(std::ostream& os, IdType const& id)
       {
         os << "(" << id.first << "," << id.second << ")";
@@ -101,11 +113,13 @@ namespace AMDiS
     using PreBasis = typename GlobalBasis::PreBasis;
     using TreePath = typename GlobalBasis::PrefixPath;
     using NodeIdSet = AMDiS::NodeIdSet<PreBasis, TreePath>;
+    using Twist = AMDiS::Twist<typename Grid::GlobalIdSet, GridView::dimension>;
 
   public:
     GlobalBasisIdSet(GlobalBasis const& globalBasis)
       : tree_(makeNode(globalBasis.preBasis(), TreePath{}))
       , nodeIdSet_(globalBasis.gridView())
+      , twist_(globalBasis.gridView().grid().globalIdSet())
     {
       Dune::Functions::initializeTree(tree_);
     }
@@ -120,8 +134,9 @@ namespace AMDiS
     {
       Dune::Functions::bindTree(tree_, element);
       nodeIdSet_.bind(tree_);
+      twist_.bind(element);
       data_.resize(size());
-      nodeIdSet_.fillIn(data_.begin());
+      nodeIdSet_.fillIn(twist_, data_.begin());
     }
 
     /// \brief unbind from the element
@@ -162,6 +177,7 @@ namespace AMDiS
   protected:
     Tree tree_;
     NodeIdSet nodeIdSet_;
+    Twist twist_;
     using Data = std::pair<IdType, PartitionType>;
     std::vector<Data> data_;
   };
@@ -195,6 +211,7 @@ namespace AMDiS
   public:
     NodeIdSet(GridView const& gridView)
       : gridView_(gridView)
+      , sizes_{}
     {}
 
     /// \brief Bind the idset to a tree node
@@ -202,6 +219,15 @@ namespace AMDiS
     {
       node_ = &node;
       size_ = node_->finiteElement().size();
+
+      std::fill(std::begin(sizes_), std::end(sizes_), 0u);
+      for (size_type i = 0; i < size_ ; ++i) {
+        Dune::LocalKey const& localKey = node_->finiteElement().localCoefficients().localKey(i);
+        sizes_[localKey.codim()]++;
+      }
+      auto refElem = Dune::referenceElement<double,GridView::dimension>(node_->element().type());
+      for (size_type c = 0; c <= GridView::dimension ; ++c)
+        sizes_[c] /= refElem.size(c);
     }
 
     /// \brief Unbind the idset
@@ -218,8 +244,8 @@ namespace AMDiS
 
     /// \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
+    template <class Twist, class It>
+    It fillIn(Twist const& twist, It it, size_type shift = 0) const
     {
       const auto& gridIdSet = gridView_.grid().globalIdSet();
 
@@ -227,11 +253,7 @@ namespace AMDiS
         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->first = {gridIdSet.subId(node_->element(), s, c), shift + twist.get(localKey,sizes_[c])};
 
         it->second = Dune::Hybrid::switchCases(std::make_index_sequence<dim+1>{}, c,
           [&](auto codim) { return node_->element().template subEntity<codim>(s).partitionType(); },
@@ -248,6 +270,7 @@ namespace AMDiS
     GridView gridView_;
     const Node* node_ = nullptr;
     size_type size_ = 0;
+    std::array<unsigned int,GridView::dimension+1> sizes_;
   };
 
 
@@ -293,13 +316,13 @@ namespace AMDiS
 
     /// \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
+    template <class Twist, class It>
+    It fillIn(Twist const& twist, 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);
+        it = subIds_.fillIn(twist, it, shift);
         shift += subTreeSize;
       }
       return it;
@@ -370,13 +393,13 @@ namespace AMDiS
 
     /// \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
+    template <class Twist, class It>
+    It fillIn(Twist const& twist, It it, size_type shift = 0) const
     {
       Tools::for_each(idsTuple_, [&](auto const& ids)
       {
         size_type subTreeSize = ids.size();
-        it = ids.fillIn(it, shift);
+        it = ids.fillIn(twist, it, shift);
         shift += subTreeSize;
       });
       return it;
@@ -443,15 +466,15 @@ namespace AMDiS
 
     /// \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
+    template <class Twist, class It>
+    It fillIn(Twist const& twist, It it, size_type shift = 0) const
     {
       for (int child = 0; child < dow; ++child) {
         size_type subTreeSize = pq2NodeIdSet_.size();
-        it = pq2NodeIdSet_.fillIn(it, shift);
+        it = pq2NodeIdSet_.fillIn(twist, it, shift);
         shift += subTreeSize;
       }
-      it = pq1NodeIdSet_.fillIn(it, shift);
+      it = pq1NodeIdSet_.fillIn(twist, it, shift);
       return it;
     }
 
@@ -496,8 +519,8 @@ namespace AMDiS
 
     /// \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
+    template <class Twist, class It>
+    It fillIn(Twist const& /*twist*/, It it, size_type shift = 0) const
     {
       const auto& gridIdSet = gridView_.grid().globalIdSet();
       auto elementId = gridIdSet.id(node_->element());