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());

src/amdis/linearalgebra/istl/Communication.hpp

@@ -79,7 +79,6 @@ namespace AMDiS
   /// Dummy implementation for ISTL-specific communication when no MPI is found
   template <class Basis>
   class ISTLCommunication
-      : public DefaultCommunication<Basis>
   {
     using Self = ISTLCommunication;
 

src/amdis/linearalgebra/mtl/ITL_Solver.hpp

@@ -456,7 +456,7 @@ namespace AMDiS
       // default iterative solver
       Map::addCreator("default", gmres);
 
-      init_direct(std::is_same<typename Dune::FieldTraits<T>::real_type, double>{});
+      init_direct(std::is_same<typename Dune::FieldTraits<typename Matrix::value_type>::real_type, double>{});
     }
 
     static void init_direct(std::false_type) {}

src/amdis/linearalgebra/mtl/UmfpackRunner.hpp

@@ -48,7 +48,6 @@ namespace AMDiS
     /// Implementation of \ref RunnerInterface::init()
     void init(Matrix const& matrix, Comm&) override
     {
-      DUNE_UNUSED_PARAMETER(comm);
       try {
         solver_.reset(new SolverType(matrix, symmetricStrategy_, allocInit_));
       } catch (mtl::mat::umfpack::error const& e) {

src/amdis/utility/CMakeLists.txt

@@ -4,5 +4,6 @@ install(FILES
     LocalToGlobalAdapter.hpp
     MacroGridFactory.hpp
     QuadratureFactory.hpp
+    Twist.hpp
     UniqueBorderPartition.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/utility)

src/amdis/utility/Twist.hpp

+#pragma once
+
+#include <vector>
+
+#include <dune/geometry/referenceelements.hh>
+#include <dune/localfunctions/common/localkey.hh>
+
+#include <amdis/Output.hpp>
+
+namespace AMDiS
+{
+  /// \brief Permutate the dof indices w.r.t. a global entity orientation
+  template <class IdSet, int dim>
+  class Twist
+  {
+    using IdType = typename IdSet::IdType;
+
+    using RefElements = typename Dune::ReferenceElements<double,dim>;
+    using RefElement = typename RefElements::ReferenceElement;
+
+  public:
+    /// Constructor. Stores a reference to the passed idSet
+    Twist(IdSet const& idSet)
+      : idSet_(idSet)
+    {}
+
+    /// Bind the twist to a codim-0 element. Calculates the global ids of all the element vertices.
+    template <class Element>
+    void bind(Element const& element)
+    {
+      static_assert(dim == Element::mydimension, "");
+      refElem_ = &RefElements::general(element.type());
+
+      ids_.resize(refElem_->size(dim));
+      for (int i = 0; i < refElem_->size(dim); ++i)
+        ids_[i] = idSet_.subId(element, i, dim);
+    }
+
+    /// Get the permutated local dof index, living on a subEntity of the bound element.
+    /**
+     * If the global entity orientation differs from the local orientation, perform a
+     * permutation of the local dof indices on that entity. This is currently implemented
+     * for edge entities only.
+     *
+     * \param localKey The Dune::LocalKey of the local dof, identifying the entity and local index
+     * \param size     The number of local dofs on that entity.
+     */
+    unsigned int get(Dune::LocalKey const& localKey, unsigned int size) const
+    {
+      int subDim = dim - localKey.codim();
+      if (subDim == 1 && 1 < dim) // facet == edge
+        return id(localKey,0) < id(localKey,1) ? localKey.index() : size - localKey.index() - 1;
+      else if (subDim == 2 && 2 < dim) { // facet == face
+        test_exit(size == 1, "Bases with more then one DoF per subentity are not yet supported.");
+        return localKey.index();
+      }
+
+      return localKey.index();
+    }
+
+  private:
+    IdType const& id(Dune::LocalKey const& localKey, int ii) const
+    {
+      return ids_[refElem_->subEntity(localKey.subEntity(), localKey.codim(), ii, dim)];
+    }
+
+  private:
+    IdSet const& idSet_;
+    RefElement const* refElem_ = nullptr;
+    std::vector<IdType> ids_;
+  };
+
+} // end namespace AMDiS

test/ISTLCommTest.cpp

@@ -40,12 +40,15 @@ bool test(Basis& basis, std::string const& prefix)
 
   // Make communicator and exchange dofs
   auto comm = Comm::create(basis, prefix);
+
+#if HAVE_MPI
   comm->copyOwnerToAll(dofs, dofs);
 
   // Compare post-exchange data with reference
   for (std::size_t i = 0; i < dofs.size(); ++i)
     if (std::abs(dofs[i]-ref[i]) > AMDIS_TEST_TOL)
       passed = false;
+#endif
 
   return passed;
 }