From 449811fe4e94002a8eba719e014cebef5e61d35b Mon Sep 17 00:00:00 2001
From: "Praetorius, Simon" <simon.praetorius@tu-dresden.de>
Date: Wed, 28 Aug 2019 12:32:46 +0200
Subject: [PATCH] added twist utilities to flip the order of edge DOFs
---
src/amdis/functions/GlobalIdSet.hpp | 67 +++++++++++------
.../linearalgebra/istl/Communication.hpp | 1 -
src/amdis/linearalgebra/mtl/ITL_Solver.hpp | 2 +-
src/amdis/linearalgebra/mtl/UmfpackRunner.hpp | 1 -
src/amdis/utility/CMakeLists.txt | 1 +
src/amdis/utility/Twist.hpp | 73 +++++++++++++++++++
test/ISTLCommTest.cpp | 3 +
7 files changed, 123 insertions(+), 25 deletions(-)
create mode 100644 src/amdis/utility/Twist.hpp
diff --git a/src/amdis/functions/GlobalIdSet.hpp b/src/amdis/functions/GlobalIdSet.hpp
index ee0674c2..29f56a69 100644
--- a/src/amdis/functions/GlobalIdSet.hpp
+++ b/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());
diff --git a/src/amdis/linearalgebra/istl/Communication.hpp b/src/amdis/linearalgebra/istl/Communication.hpp
index eb63684f..fc56b2b6 100644
--- a/src/amdis/linearalgebra/istl/Communication.hpp
+++ b/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;
diff --git a/src/amdis/linearalgebra/mtl/ITL_Solver.hpp b/src/amdis/linearalgebra/mtl/ITL_Solver.hpp
index fd800d0c..bc2360a6 100644
--- a/src/amdis/linearalgebra/mtl/ITL_Solver.hpp
+++ b/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) {}
diff --git a/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp b/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
index 546d25d0..368b3f41 100644
--- a/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
+++ b/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) {
diff --git a/src/amdis/utility/CMakeLists.txt b/src/amdis/utility/CMakeLists.txt
index 33add258..ec1450c4 100644
--- a/src/amdis/utility/CMakeLists.txt
+++ b/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)
diff --git a/src/amdis/utility/Twist.hpp b/src/amdis/utility/Twist.hpp
new file mode 100644
index 00000000..d6ef68b1
--- /dev/null
+++ b/src/amdis/utility/Twist.hpp
@@ -0,0 +1,73 @@
+#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
diff --git a/test/ISTLCommTest.cpp b/test/ISTLCommTest.cpp
index e3604680..16fcd63f 100644
--- a/test/ISTLCommTest.cpp
+++ b/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;
}
--
GitLab