From 502e1367e704a6bdafa900d2bb3be01979281f75 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Wed, 21 Aug 2019 13:55:30 +0200 Subject: [PATCH 1/2] added twist utilities to flip the order of edge DOFs --- src/amdis/functions/GlobalIdSet.hpp | 67 +++++++++++++++++--------- src/amdis/utility/CMakeLists.txt | 1 + src/amdis/utility/Twist.hpp | 73 +++++++++++++++++++++++++++++ 3 files changed, 119 insertions(+), 22 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 #include +#include #include #include #include #include -#include +#include namespace Dune { @@ -88,9 +89,20 @@ namespace AMDiS { using Super = std::pair; - 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; + using Twist = AMDiS::Twist; 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; std::vector 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(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 - It fillIn(It it, size_type shift = 0) const + template + 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{}, c, [&](auto codim) { return node_->element().template subEntity(s).partitionType(); }, @@ -248,6 +270,7 @@ namespace AMDiS GridView gridView_; const Node* node_ = nullptr; size_type size_ = 0; + std::array 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 - It fillIn(It it, size_type shift = 0) const + template + 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 - It fillIn(It it, size_type shift = 0) const + template + 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 - It fillIn(It it, size_type shift = 0) const + template + 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 - It fillIn(It it, size_type shift = 0) const + template + 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/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..4488b3da --- /dev/null +++ b/src/amdis/utility/Twist.hpp @@ -0,0 +1,73 @@ +#pragma once + +#include + +#include +#include + +#include + +namespace AMDiS +{ + /// \brief Permutate the dof indices w.r.t. a global entity orientation + template + class Twist + { + using IdType = typename IdSet::IdType; + + using RefElements = typename Dune::ReferenceElements; + 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 + 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 ids_; + }; + +} // end namespace AMDiS -- GitLab From 850473a6eef198f601893a882fdcdf64cdddd7d4 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Wed, 28 Aug 2019 11:06:59 +0200 Subject: [PATCH 2/2] Include corrected and some connected compile errors fixed --- src/amdis/linearalgebra/istl/Communication.hpp | 1 - src/amdis/linearalgebra/mtl/ITL_Solver.hpp | 2 +- src/amdis/linearalgebra/mtl/UmfpackRunner.hpp | 1 - src/amdis/utility/Twist.hpp | 2 +- test/ISTLCommTest.cpp | 3 +++ 5 files changed, 5 insertions(+), 4 deletions(-) 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 ISTLCommunication - : public DefaultCommunication { 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::real_type, double>{}); + init_direct(std::is_same::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 0f5729a6..16a88f8c 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/Twist.hpp b/src/amdis/utility/Twist.hpp index 4488b3da..d6ef68b1 100644 --- a/src/amdis/utility/Twist.hpp +++ b/src/amdis/utility/Twist.hpp @@ -5,7 +5,7 @@ #include #include -#include +#include 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