diff --git a/src/amdis/common/TupleUtility.hpp b/src/amdis/common/TupleUtility.hpp index 2f6aaead244295067097262ac55f60dcae9372ff..e7eca21df9dde8df0015259d41b9eab0627708f2 100644 --- a/src/amdis/common/TupleUtility.hpp +++ b/src/amdis/common/TupleUtility.hpp @@ -32,4 +32,28 @@ namespace AMDiS } // end namespace Impl + template class Map> + struct MapTuple; + + template class Map> + using MapTuple_t = typename MapTuple::type; + + template class Map> + struct MapTuple, Map> + { + using type = std::tuple...>; + }; + + template class Map> + struct IndexMapTuple; + + template class Map> + using IndexMapTuple_t = typename IndexMapTuple::type; + + template class Map> + struct IndexMapTuple, Map> + { + using type = std::tuple...>; + }; + } // end namespace AMDiS diff --git a/src/amdis/functions/CMakeLists.txt b/src/amdis/functions/CMakeLists.txt index 13244525b3280e8c2abc2f98606792edf8d03e57..3d89bca76501b246d622936a54864643890d54bb 100644 --- a/src/amdis/functions/CMakeLists.txt +++ b/src/amdis/functions/CMakeLists.txt @@ -1,6 +1,8 @@ #install headers install(FILES + GlobalIdSet.hpp HierarchicNodeToRangeMap.hpp Interpolate.hpp + Nodes.hpp DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/functions) diff --git a/src/amdis/functions/GlobalIdSet.hpp b/src/amdis/functions/GlobalIdSet.hpp new file mode 100644 index 0000000000000000000000000000000000000000..b2c089d0558f76a1b4bfa894164326c828bf58b9 --- /dev/null +++ b/src/amdis/functions/GlobalIdSet.hpp @@ -0,0 +1,508 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace Dune +{ + namespace Functions + { + // forward declarations... + template + class SubspaceBasis; + + template + class CompositePreBasis; + + template + class LagrangePreBasis; + + template + class PowerPreBasis; + + template + class LagrangeDGPreBasis; + + template + class TaylorHoodPreBasis; + } +} + + +namespace AMDiS +{ + // forward declaration + template + 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 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 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; + using PartitionType = Dune::PartitionType; + + using PreBasis = typename GlobalBasis::PreBasis; + using TreePath = typename GlobalBasis::PrefixPath; + using NodeIdSet = AMDiS::NodeIdSet; + + 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; + std::vector data_; + }; + + + template + class GlobalBasisIdSet> + : public GlobalBasisIdSet + { + public: + GlobalBasisIdSet(Dune::Functions::SubspaceBasis const& basis) + : GlobalBasisIdSet(basis.rootBasis()) + {} + }; + + + template + class NodeIdSet + { + public: + using PreBasis = PB; + using Node = Node_t; + 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 + 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{}, c, + [&](auto codim) { return node_->element().template subEntity(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 NodeIdSet, TP> + { + public: + using PreBasis = Dune::Functions::PowerPreBasis; + using Node = Node_t; + 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(), std::size_t(0))); + using SubNodeIdSet = NodeIdSet; + 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 + 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 NodeIdSet, TP> + { + public: + using PreBasis = Dune::Functions::CompositePreBasis; + using Node = Node_t; + 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; + + // The I'th SubPreBasis + template + using SubPreBasis = typename PreBasis::template SubPreBasis; + + template + using SubTreePath = decltype(Dune::TypeTree::push_back(std::declval(), Dune::index_constant{})); + + template + using SubNodeIdSet = NodeIdSet, SubTreePath>; + + // A tuple of NodeIdSets for the SubPreBases + using IdsTuple = IndexMapTuple_t, SubNodeIdSet>; + + public: + NodeIdSet(GridView const& gridView) + : idsTuple_(Tools::apply_indices([&](auto... i) { + return std::make_tuple(SubNodeIdSet(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(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 + 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 NodeIdSet, TP> + { + public: + using PreBasis = Dune::Functions::TaylorHoodPreBasis; + using Node = Node_t; + using GridView = typename PreBasis::GridView; + using size_type = std::size_t; + + private: + using PTP = decltype(Dune::TypeTree::push_back(std::declval(), Dune::index_constant<1>{})); + using VTP = decltype(Dune::TypeTree::push_back(std::declval(), Dune::index_constant<0>{})); + using V0TP = decltype(Dune::TypeTree::push_back(std::declval(), std::size_t(0))); + + // Note: PreBasis::PQ1PreBasis is private + using PQ1PreBasis = typename NodeIndexSet_t::PQ1NodeIndexSet::PreBasis; + using PQ2PreBasis = typename NodeIndexSet_t::PQ2NodeIndexSet::PreBasis; + + using PQ1NodeIdSet = NodeIdSet; // pressure + using PQ2NodeIdSet = NodeIdSet; // 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 + 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 NodeIdSet, TP> + { + public: + using PreBasis = Dune::Functions::LagrangeDGPreBasis; + using Node = Node_t; + 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 + 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 diff --git a/src/amdis/functions/Nodes.hpp b/src/amdis/functions/Nodes.hpp new file mode 100644 index 0000000000000000000000000000000000000000..5a247083d66b41fe028d9d0a2c3756c8c3e2ada1 --- /dev/null +++ b/src/amdis/functions/Nodes.hpp @@ -0,0 +1,37 @@ +#pragma once + +#include + +namespace AMDiS +{ + // dune version independent extraction of node type from preBasis + template + using Node_t = +#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) + typename PB::template Node; +#else + typename PB::Node; +#endif + + // dune version independent extraction of node indexset type from preBasis + template + using NodeIndexSet_t = +#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) + typename PB::template IndexSet; +#else + typename PB::IndexSet; +#endif + + // dune version independent creation of node from preBasis + template + 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 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 282214b3a68816de41ce0f583939d8f86a0bc0bf..9812673160d54ea53efffc0af0a9a881c013ef61 100644 --- a/test/CMakeLists.txt +++ b/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) diff --git a/test/GlobalIdSetTest.cpp b/test/GlobalIdSetTest.cpp new file mode 100644 index 0000000000000000000000000000000000000000..b332cf98fa44730c6d904d70dfb20f5d0096dbd9 --- /dev/null +++ b/test/GlobalIdSetTest.cpp @@ -0,0 +1,81 @@ +#include +#include +#include + +#include +#include + +#include "Tests.hpp" + +using namespace AMDiS; + +template +void checkBasisIds(const Basis& basis) +{ + using IdSet = GlobalBasisIdSet; + using IdType = typename IdSet::IdType; + + std::set cache{}; + + IdSet 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); + cache.insert(id); + + DUNE_UNUSED auto pt = idSet.partitionType(i); + } + + idSet.unbind(); + } + + AMDIS_TEST(cache.size() == basis.dimension()); +} + +template +void checkGrid(const Grid& grid) +{ + // create basis + using namespace Dune::Functions::BasisBuilder; + auto basis1 = makeBasis(grid.leafGridView(), + composite( + power<2>(lagrange<2>()), + lagrange<1>() + )); + + auto basis2 = makeBasis(grid.leafGridView(), + composite( + power<2>(lagrangeDG<2>()), + lagrangeDG<1>() + )); + + auto basis3 = makeBasis(grid.leafGridView(), taylorHood()); + + checkBasisIds(basis1); + checkBasisIds(basis2); + checkBasisIds(basis3); +} + +int main(int argc, char** argv) +{ + Environment env(argc, argv); + + // non-overlapping + Dune::YaspGrid<2> grid1({1.0, 1.0}, {8,8}, 0, 0); + Dune::YaspGrid<3> grid2({1.0, 1.0, 1.0}, {8,8,8}, 0, 0); + + // overlapping + Dune::YaspGrid<2> grid3({1.0, 1.0}, {8,8}, 0, 1); + Dune::YaspGrid<3> grid4({1.0, 1.0, 1.0}, {8,8,8}, 0, 1); + + checkGrid(grid1); + checkGrid(grid2); + checkGrid(grid3); + checkGrid(grid4); + + return report_errors(); +}