From 28549e8784c670f48dbb80c23df6f23036ee7701 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Mon, 25 Mar 2019 22:45:08 +0100 Subject: [PATCH 1/4] reimplement interpolate function with averaging --- src/amdis/common/FieldMatVec.hpp | 17 + src/amdis/common/FieldMatVec.inc.hpp | 54 +++ src/amdis/functions/Interpolate.hpp | 401 ++++++++++++++++++++++ src/amdis/gridfunctions/DOFVectorView.hpp | 35 +- 4 files changed, 499 insertions(+), 8 deletions(-) create mode 100644 src/amdis/functions/Interpolate.hpp diff --git a/src/amdis/common/FieldMatVec.hpp b/src/amdis/common/FieldMatVec.hpp index dbc3f6c5..ceb547bf 100644 --- a/src/amdis/common/FieldMatVec.hpp +++ b/src/amdis/common/FieldMatVec.hpp @@ -398,15 +398,32 @@ namespace Dune template FieldMatrix,M,N> multiplies_ABt(FieldMatrix const& A, FieldMatrix const& B); + template FieldMatrix& multiplies_ABt(FieldMatrix const& A, FieldMatrix const& B, FieldMatrix& C); + template + FieldVector& multiplies_ABt(FieldMatrix const& A, FieldMatrix const& B, FieldVector& C); + + template + FieldVector& multiplies_ABt(FieldVector const& A, FieldMatrix const& B, FieldVector& C); + + template + FieldMatrix& multiplies_ABt(FieldVector const& A, FieldMatrix const& B, FieldMatrix& C); + + template FieldMatrix& multiplies_ABt(FieldMatrix const& A, DiagonalMatrix const& B, FieldMatrix& C); template FieldVector& multiplies_ABt(FieldMatrix const& A, DiagonalMatrix const& B, FieldVector& C); + template + FieldVector& multiplies_ABt(FieldVector const& A, DiagonalMatrix const& B, FieldVector& C); + + template + FieldMatrix& multiplies_ABt(FieldVector const& A, DiagonalMatrix const& B, FieldMatrix& C); + // ----------------------------------------------------------------------------- template diff --git a/src/amdis/common/FieldMatVec.inc.hpp b/src/amdis/common/FieldMatVec.inc.hpp index 340b24e6..93f0c219 100644 --- a/src/amdis/common/FieldMatVec.inc.hpp +++ b/src/amdis/common/FieldMatVec.inc.hpp @@ -517,6 +517,41 @@ FieldMatrix& multiplies_ABt(FieldMatrix const& A, FieldMatrix return C; } +template +FieldVector& multiplies_ABt(FieldMatrix const& A, FieldMatrix const& B, FieldVector& C) +{ + for (int n = 0; n < N; ++n) { + C[n] = 0; + for (int l = 0; l < L; ++l) + C[n] += A[0][l] * B[n][l]; + } + return C; +} + +template +FieldVector& multiplies_ABt(FieldVector const& A, FieldMatrix const& B, FieldVector& C) +{ + for (int n = 0; n < N; ++n) { + C[n] = 0; + for (int l = 0; l < L; ++l) + C[n] += A[l] * B[n][l]; + } + return C; +} + +template +FieldMatrix& multiplies_ABt(FieldVector const& A, FieldMatrix const& B, FieldMatrix& C) +{ + for (int n = 0; n < N; ++n) { + C[0][n] = 0; + for (int l = 0; l < L; ++l) + C[0][n] += A[l] * B[n][l]; + } + return C; +} + + + template FieldMatrix& multiplies_ABt(FieldMatrix const& A, DiagonalMatrix const& B, FieldMatrix& C) { @@ -537,6 +572,25 @@ FieldVector& multiplies_ABt(FieldMatrix const& A, DiagonalMatri return C; } +template +FieldVector& multiplies_ABt(FieldVector const& A, DiagonalMatrix const& B, FieldVector& C) +{ + for (int n = 0; n < N; ++n) { + C[n] = A[n] * B.diagonal(n); + } + return C; +} + +template +FieldMatrix& multiplies_ABt(FieldVector const& A, DiagonalMatrix const& B, FieldMatrix& C) +{ + for (int n = 0; n < N; ++n) { + C[0][n] = A[n] * B.diagonal(n); + } + return C; +} + + template T const& at(FieldMatrix const& vec, std::size_t i) { diff --git a/src/amdis/functions/Interpolate.hpp b/src/amdis/functions/Interpolate.hpp new file mode 100644 index 00000000..de02d40a --- /dev/null +++ b/src/amdis/functions/Interpolate.hpp @@ -0,0 +1,401 @@ +#pragma once + +#include +#include + +#include +#include + +#include +#include +#include +#include + +namespace AMDiS { +namespace Impl { + + struct DummyCounter + { + struct Sink + { + template constexpr Sink& operator =(T const&) { return *this; } + template constexpr Sink& operator+=(T const&) { return *this; } + template constexpr Sink& operator-=(T const&) { return *this; } + + constexpr operator int() const { return 1; } + }; + + template + constexpr void resize(SizeInfo const& sizeInfo) { /* do nothing */ } + constexpr std::size_t size() const { return 0; } + + template constexpr Sink operator[](MI const&) { return Sink{}; } + template constexpr int operator[](MI const&) const { return 1; } + }; + + + template + class LocalInterpolateVisitor + { + public: + using Basis = B; + using LocalView = typename Basis::LocalView; + using MultiIndex = typename LocalView::MultiIndex; + + using NodeToRangeEntry = NTRE; + using VectorBackend = Vec; + using CounterBackend = C; + using BitVectorBackend = BV; + using LocalFunction = LF; + + using GridView = typename Basis::GridView; + using Element = typename GridView::template Codim<0>::Entity; + using LocalDomain = typename Element::Geometry::LocalCoordinate; + + template + class LocalFunctionComponent + : public Dune::LocalFiniteElementFunctionBase::type + { + using FiniteElement = typename Node::FiniteElement; + using Range = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType; + + public: + LocalFunctionComponent(Node const& node, TreePath const& treePath, LF const& localF, NTRE const& nodeToRangeEntry) + : node_(node) + , treePath_(treePath) + , localF_(localF) + , nodeToRangeEntry_(nodeToRangeEntry) + {} + + void evaluate(const LocalDomain& x, Range& y) const + { + const auto& tmp = localF_(x); + const auto& tmp_vec = Dune::MatVec::as_vector(tmp); + y = Dune::Functions::flatVectorView(nodeToRangeEntry_(node_, treePath_, tmp_vec))[comp_]; + } + + void setComponent(std::size_t comp) + { + comp_ = comp; + } + + private: + Node const& node_; + TreePath const& treePath_; + LocalFunction const& localF_; + NodeToRangeEntry const& nodeToRangeEntry_; + + std::size_t comp_ = 0; + }; + + public: + LocalInterpolateVisitor(Vec& vector, C& counter, BV const& bitVector, + LF const& localF, LocalView const& localView, + NodeToRangeEntry const& nodeToRangeEntry) + : vector_(vector) + , counter_(counter) + , bitVector_(bitVector) + , localF_(localF) + , localView_(localView) + , nodeToRangeEntry_(nodeToRangeEntry) + { + static_assert(Concepts::Callable, + "Function passed to LocalInterpolateVisitor does not model the Callable concept"); + } + + template + void operator()(Node const& node, TreePath const& treePath) + { + using FiniteElement = typename Node::FiniteElement; + using RangeField = typename FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType; + + auto&& fe = node.finiteElement(); + std::size_t feSize = fe.localBasis().size(); + + thread_local std::vector visit; + visit.resize(feSize, false); + + bool visit_any = false; + for (std::size_t i = 0; i < feSize; ++i) + { + auto multiIndex = localView_.index(node.localIndex(i)); + if (bitVector_[multiIndex]) { + visit[i] = true; + visit_any = true; + if (average) + counter_[multiIndex] += 1; + } else { + visit[i] = false; + } + } + + if (!visit_any) + return; + + LocalFunctionComponent localFj(node, treePath, localF_, nodeToRangeEntry_); + thread_local std::vector interpolationCoefficients; + + std::size_t blockSize = Dune::Functions::flatVectorView(vector_[localView_.index(0)]).size(); + for (std::size_t j = 0; j < blockSize; ++j) + { + localFj.setComponent(j); + fe.localInterpolation().interpolate(localFj, interpolationCoefficients); + assert(interpolationCoefficients.size() == feSize); + for (std::size_t i = 0; i < feSize; ++i) + { + if (visit[i]) + { + auto multiIndex = localView_.index(node.localIndex(i)); + auto vectorBlock = Dune::Functions::flatVectorView(vector_[multiIndex]); + if (average && counter_[multiIndex] > 1) + vectorBlock[j] += interpolationCoefficients[i]; + else + vectorBlock[j] = interpolationCoefficients[i]; + } + } + } + } + + protected: + VectorBackend& vector_; + CounterBackend& counter_; + BitVectorBackend const& bitVector_; + + LocalFunction const& localF_; + LocalView const& localView_; + NodeToRangeEntry const& nodeToRangeEntry_; + }; + + // Small helper functions to wrap vectors using istlVectorBackend + // if they do not already satisfy the VectorBackend interface. + template + decltype(auto) toVectorBackend(B const& basis, Vec& vec) + { + return Dune::Hybrid::ifElse(Dune::models, Vec>(), + [&](auto id) -> decltype(auto) { return vec; }, + [&](auto id) -> decltype(auto) { return Dune::Functions::istlVectorBackend(vec); }); + } + + template + decltype(auto) toConstVectorBackend(B const& basis, Vec const& vec) + { + return Dune::Hybrid::ifElse(Dune::models, Vec>(), + [&](auto id) -> decltype(auto) { return vec; }, + [&](auto id) -> decltype(auto) { return Dune::Functions::istlVectorBackend(vec); }); + } + +} // namespace Impl + + +/** + * \brief Interpolate given function in discrete function space + * + * Interpolation is done wrt the leaf node of the ansatz tree + * corresponding to the given tree path. + * + * Notice that this will only work if the range type of f and + * the block type of coeff are compatible and supported by + * flatVectorView. + * + * \param basis Global function space basis of discrete function space + * \param treePath Tree path specifying the part of the ansatz tree to use + * \param coeff Coefficient vector to represent the interpolation + * \param f Function to interpolate + * \param nodeToRangeEntry Polymorphic functor mapping local ansatz nodes to range-indices of given function + * \param bitVector A vector with flags marking all DOFs that should be interpolated + */ +template +void interpolateTreeSubset(B const& basis, TP const& treePath, Vec& vec, C& count, BV const& bitVec, + GF const& gf, NTRE const& nodeToRangeEntry, bool_t) +{ + auto&& vector = Impl::toVectorBackend(basis,vec); + auto&& counter = Impl::toVectorBackend(basis,count); + auto&& bitVector = Impl::toConstVectorBackend(basis,bitVec); + vector.resize(sizeInfo(basis)); + counter.resize(sizeInfo(basis)); + + // Obtain a local view of f + auto lf = localFunction(gf); + auto localView = basis.localView(); + + for (const auto& e : elements(basis.gridView())) + { + localView.bind(e); + lf.bind(e); + + auto&& subTree = Dune::TypeTree::child(localView.tree(),treePath); + + using Visitor + = Impl::LocalInterpolateVisitor; + for_each_leaf_node(subTree, Visitor{vector, counter, bitVector, lf, localView, nodeToRangeEntry}); + } +} + + +template ,B>()), + REQUIRES(Concepts::GridFunction)> +void interpolateTreeSubset(B const& basis, Dune::TypeTree::HybridTreePath const& treePath, + Vec& vec, C& count, BV const& bitVec, GF const& gf, bool_t) +{ + auto ntrm = Dune::Functions::HierarchicNodeToRangeMap(); + AMDiS::interpolateTreeSubset(basis, treePath, vec, count, bitVec, gf, ntrm, bool_t{}); +} + + +template ,B>()), + REQUIRES(Concepts::GridFunction)> +void interpolateTree(B const& basis, Dune::TypeTree::HybridTreePath const& treePath, + Vec& vec, C& count, GF const& gf, NTRE const& nodeToRangeEntry, bool_t) +{ + auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector(); + AMDiS::interpolateTreeSubset(basis, treePath, vec, count, bitVec, gf, nodeToRangeEntry, bool_t{}); +} + + +template ,B>()), + REQUIRES(Concepts::GridFunction)> +void interpolateTree(B const& basis, Dune::TypeTree::HybridTreePath const& treePath, + Vec& vec, C& count, GF const& gf, bool_t) +{ + auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector(); + auto ntrm = Dune::Functions::HierarchicNodeToRangeMap(); + AMDiS::interpolateTreeSubset(basis, treePath, vec, count, bitVec, gf, ntrm, bool_t{}); +} + + +/** + * \brief Interpolate given function in discrete function space + * + * Interpolation is done wrt the leaf node of the ansatz tree + * corresponding to the given tree path. + * + * Notice that this will only work if the range type of f and + * the block type of coeff are compatible and supported by + * flatVectorView. + * + * \param basis Global function space basis of discrete function space + * \param treePath Tree path specifying the part of the ansatz tree to use + * \param vec Coefficient vector to represent the interpolation + * \param count Vector that counts for the number of value assignments + * \param bitVec A vector with flags marking all DOFs that should be interpolated + * \param gf GridFunction to interpolate + */ +template ,B>()), + REQUIRES(Concepts::GridFunction)> +void interpolateFiltered(B const& basis, Dune::TypeTree::HybridTreePath const& treePath, + Vec& vec, C& count, BV const& bitVec, GF const& gf, bool_t) +{ + auto ntrm = Dune::Functions::HierarchicNodeToRangeMap(); + AMDiS::interpolateTreeSubset(basis, treePath, vec, count, bitVec, gf, ntrm, bool_t{}); +} + + +/** + * \brief Interpolate given function in discrete function space + * + * Interpolation is done wrt the leaf node of the ansatz tree + * corresponding to the given tree path. Only vector coefficients marked as 'true' in the + * bitVector argument are interpolated. Use this, e.g., to interpolate Dirichlet boundary values. + * + * Notice that this will only work if the range type of f and + * the block type of coeff are compatible and supported by + * flatVectorView. + * + * \param basis Global function space basis of discrete function space + * \param vec Coefficient vector to represent the interpolation + * \param count Vector that counts for the number of value assignments + * \param bitVec A vector with flags marking all DOFs that should be interpolated + * \param gf GridFunction to interpolate + */ +template ,B>()), + REQUIRES(not Dune::Functions::Imp::isHybridTreePath()), + REQUIRES(Concepts::GridFunction)> +void interpolateFiltered(B const& basis, Vec& vec, C& count, BV const& bitVec, GF const& gf, bool_t) +{ + auto root = Dune::TypeTree::hybridTreePath(); + auto ntrm = Dune::Functions::HierarchicNodeToRangeMap(); + AMDiS::interpolateTreeSubset(basis, root, vec, count, bitVec, gf, ntrm, bool_t{}); +} + + +/** + * \brief Interpolate given function in discrete function space + * + * Notice that this will only work if the range type of f and + * the block type of coeff are compatible and supported by + * flatVectorView. + * + * This function will only work, if the local ansatz tree of + * the basis is trivial, i.e., a single leaf node. + * + * \param basis Global function space basis of discrete function space + * \param vec Coefficient vector to represent the interpolation + * \param count Vector that counts for the number of value assignments + * \param gf GridFunction to interpolate + */ +template ,B>()), + REQUIRES(not Dune::Functions::Imp::isHybridTreePath()), + REQUIRES(Concepts::GridFunction)> +void interpolate(B const& basis, Vec& vec, C& count, GF const& gf, bool_t) +{ + auto root = Dune::TypeTree::hybridTreePath(); + auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector(); + AMDiS::interpolateFiltered(basis, root, vec, count, bitVec, gf, bool_t{}); +} + +template ,B>()), + REQUIRES(not Dune::Functions::Imp::isHybridTreePath()), + REQUIRES(Concepts::GridFunction)> +void interpolate(B const& basis, Vec& vec, GF const& gf) +{ + auto root = Dune::TypeTree::hybridTreePath(); + auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector(); + auto count = Impl::DummyCounter(); + AMDiS::interpolateFiltered(basis, root, vec, count, bitVec, gf, std::false_type{}); +} + + +/** + * \brief Interpolate given function in discrete function space + * + * Interpolation is done wrt the leaf node of the ansatz tree + * corresponding to the given tree path. + * + * Notice that this will only work if the range type of f and + * the block type of corresponding coeff entries are compatible + * and supported by flatVectorView. + * + * \param basis Global function space basis of discrete function space + * \param treePath Tree path specifying the part of the ansatz tree to use + * \param vec Coefficient vector to represent the interpolation + * \param count Vector that counts for the number of value assignments + * \param gf GridFunction to interpolate + */ +template ,B>()), + REQUIRES(Concepts::GridFunction)> +void interpolate(B const& basis, Dune::TypeTree::HybridTreePath const& treePath, + Vec& vec, C& count, GF const& gf, bool_t) +{ + auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector(); + AMDiS::interpolateFiltered(basis, treePath, vec, count, bitVec, gf, bool_t{}); +} + +template ,B>()), + REQUIRES(Concepts::GridFunction)> +void interpolate(B const& basis, Dune::TypeTree::HybridTreePath const& treePath, Vec& vec, GF const& gf) +{ + auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector(); + auto count = Impl::DummyCounter(); + AMDiS::interpolateFiltered(basis, treePath, vec, count, bitVec, gf, std::false_type{}); +} + +} // end namespace AMDiS diff --git a/src/amdis/gridfunctions/DOFVectorView.hpp b/src/amdis/gridfunctions/DOFVectorView.hpp index f95677a9..722c9098 100644 --- a/src/amdis/gridfunctions/DOFVectorView.hpp +++ b/src/amdis/gridfunctions/DOFVectorView.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include #include namespace AMDiS @@ -23,8 +24,9 @@ namespace AMDiS {} /// Constructor. Stores a pointer to the mutable `dofvector`. - DOFVectorView(DOFVector& dofVector, TP const& treePath) - : Super(dofVector, treePath) + template + DOFVectorView(DOFVector& dofVector, PreTreePath const& preTreePath) + : Super(dofVector, makeTreePath(preTreePath)) , mutableDofVector_(&dofVector) {} @@ -38,15 +40,27 @@ namespace AMDiS * v.interpolate_noalias([](auto const& x) { return x[0]; }); * ``` **/ - template - void interpolate_noalias(Expr&& expr) + template + void interpolate_noalias(Expr&& expr, bool_t = {}) { auto const& basis = this->basis(); auto const& treePath = this->treePath(); auto&& gridFct = makeGridFunction(FWD(expr), basis.gridView()); - Dune::Functions::interpolate(basis, treePath, coefficients(), FWD(gridFct)); + if (average) { + thread_local std::vector counter; + counter.clear(); + AMDiS::interpolate(basis, treePath, coefficients(), counter, FWD(gridFct), std::true_type{}); + + auto& coeff = coefficients().vector(); + for (std::size_t i = 0; i < counter.size(); ++i) { + if (counter[i] > 0) + coeff[i] /= double(counter[i]); + } + } else { + AMDiS::interpolate(basis, treePath, coefficients(), FWD(gridFct)); + } } /// \brief Interpolation of GridFunction to DOFVector @@ -59,13 +73,13 @@ namespace AMDiS * Allows to have a reference to the DOFVector in the expression, e.g. as * \ref DiscreteFunction or \ref gradientAtQP() of a DiscreteFunction. **/ - template - void interpolate(Expr&& expr) + template + void interpolate(Expr&& expr, bool_t = {}) { // create temporary copy of data DOFVector tmp(coefficients()); Self tmpView{tmp, this->treePath()}; - tmpView.interpolate_noalias(FWD(expr)); + tmpView.interpolate_noalias(FWD(expr), bool_t{}); // move data from temporary vector into stored DOFVector coefficients().vector() = std::move(tmp.vector()); @@ -107,6 +121,11 @@ namespace AMDiS #if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION + // Deduction guide for DOFVectorView class + template + DOFVectorView(DOFVector& dofVector, PreTreePath const& preTreePath) + -> DOFVectorView; + // Deduction guide for DOFVectorView class template DOFVectorView(DOFVector& dofVector) -- GitLab From 22049de1a88f0bc72a1ac3997a94ff5100ed7d26 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Tue, 26 Mar 2019 10:31:43 +0100 Subject: [PATCH 2/4] add amdis/functions to CMakeLists --- src/amdis/CMakeLists.txt | 1 + src/amdis/functions/CMakeLists.txt | 5 +++++ 2 files changed, 6 insertions(+) create mode 100644 src/amdis/functions/CMakeLists.txt diff --git a/src/amdis/CMakeLists.txt b/src/amdis/CMakeLists.txt index 706ce91d..3f3e5861 100644 --- a/src/amdis/CMakeLists.txt +++ b/src/amdis/CMakeLists.txt @@ -65,6 +65,7 @@ install(FILES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis) add_subdirectory("common") +add_subdirectory("functions") add_subdirectory("gridfunctions") add_subdirectory("linearalgebra") add_subdirectory("localoperators") diff --git a/src/amdis/functions/CMakeLists.txt b/src/amdis/functions/CMakeLists.txt new file mode 100644 index 00000000..98ce0646 --- /dev/null +++ b/src/amdis/functions/CMakeLists.txt @@ -0,0 +1,5 @@ +#install headers + +install(FILES + Interpolate.hpp +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/functions) -- GitLab From 658ccbc0c366f31f45111bb6f3d231768e647072 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Wed, 27 Mar 2019 15:29:30 +0100 Subject: [PATCH 3/4] interpolate using own HierarchicNodeToRangeMap and a wrapper for TreePath types --- src/amdis/functions/CMakeLists.txt | 1 + .../functions/HierarchicNodeToRangeMap.hpp | 41 +++++++++++++++++ src/amdis/functions/Interpolate.hpp | 45 ++++++++++++++++--- 3 files changed, 82 insertions(+), 5 deletions(-) create mode 100644 src/amdis/functions/HierarchicNodeToRangeMap.hpp diff --git a/src/amdis/functions/CMakeLists.txt b/src/amdis/functions/CMakeLists.txt index 98ce0646..13244525 100644 --- a/src/amdis/functions/CMakeLists.txt +++ b/src/amdis/functions/CMakeLists.txt @@ -1,5 +1,6 @@ #install headers install(FILES + HierarchicNodeToRangeMap.hpp Interpolate.hpp DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/functions) diff --git a/src/amdis/functions/HierarchicNodeToRangeMap.hpp b/src/amdis/functions/HierarchicNodeToRangeMap.hpp new file mode 100644 index 00000000..25971440 --- /dev/null +++ b/src/amdis/functions/HierarchicNodeToRangeMap.hpp @@ -0,0 +1,41 @@ +#pragma once + +#include +#include + +#include + +#include +#include + +#include + +namespace AMDiS +{ + /** + * \brief A simple node to range map using the nested tree indices + * + * This map directly usses the tree path entries of the given + * node to access the nested container. + * + * If the container does not provide any operator[] access, + * it is simply forwarded for all nodes. + */ + struct HierarchicNodeToRangeMap + { + template >())> + decltype(auto) operator()(const Node& node, const TreePath& treePath, Range&& y) const + { + return Dune::Functions::resolveStaticMultiIndex(y, treePath); + } + + template >())> + decltype(auto) operator()(const Node& node, const TreePath& treePath, Range&& y) const + { + return std::forward(y); + } + }; + +} // end namespace AMDiS diff --git a/src/amdis/functions/Interpolate.hpp b/src/amdis/functions/Interpolate.hpp index de02d40a..09d2ed06 100644 --- a/src/amdis/functions/Interpolate.hpp +++ b/src/amdis/functions/Interpolate.hpp @@ -3,12 +3,15 @@ #include #include +#include +#include #include #include #include #include #include +#include #include namespace AMDiS { @@ -34,6 +37,16 @@ namespace Impl { }; + /// \brief Visitor evalued on the leaf nodes of basis-tree + /** + * \tparam B GlobalBasis + * \tparam Vec The Coefficient vector + * \tparam C A counter vector(-like) datastructure + * \tparam BV BitVector indicating which DOFs to visit + * \tparam LF LocalFunction to evaluate in the local interpolation + * \tparam NTRE A node-to-range-map, by default \ref HierarchicNodeToRangeMap + * \tparam average Indicates whether to do value averaging on shared DOFs (true), or simple assignment. + **/ template class LocalInterpolateVisitor { @@ -52,6 +65,7 @@ namespace Impl { using Element = typename GridView::template Codim<0>::Entity; using LocalDomain = typename Element::Geometry::LocalCoordinate; + /// Functor called in the LocalInterpolation template class LocalFunctionComponent : public Dune::LocalFiniteElementFunctionBase::type @@ -71,7 +85,7 @@ namespace Impl { { const auto& tmp = localF_(x); const auto& tmp_vec = Dune::MatVec::as_vector(tmp); - y = Dune::Functions::flatVectorView(nodeToRangeEntry_(node_, treePath_, tmp_vec))[comp_]; + y = Dune::Functions::flatVectorView(nodeToRangeEntry_(node_, transformTreePath(treePath_), tmp_vec))[comp_]; } void setComponent(std::size_t comp) @@ -79,6 +93,20 @@ namespace Impl { comp_ = comp; } +#if DUNE_VERSION_GT(DUNE_FUNCTIONS,2,6) + static TreePath const& transformTreePath(TreePath const& treePath) + { + return treePath; + } +#else + // NOTE: due to a bug in dune-functions <= 2.6, a hybrid-treepath can not be passed to a HierarchicNodeToRangeMap, + // i.e. the HybridTreePath missed a size() function + static auto transformTreePath(TreePath const& treePath) + { + return Tools::apply([](auto... i) { return Dune::makeTupleVector(i...); }, treePath._data); + } +#endif + private: Node const& node_; TreePath const& treePath_; @@ -89,6 +117,7 @@ namespace Impl { }; public: + /// Constructor. Stores references to all passed objects. LocalInterpolateVisitor(Vec& vector, C& counter, BV const& bitVector, LF const& localF, LocalView const& localView, NodeToRangeEntry const& nodeToRangeEntry) @@ -103,6 +132,7 @@ namespace Impl { "Function passed to LocalInterpolateVisitor does not model the Callable concept"); } + /// Apply the visitor to a node in the basis-tree (with corresponding treepath) template void operator()(Node const& node, TreePath const& treePath) { @@ -115,6 +145,8 @@ namespace Impl { thread_local std::vector visit; visit.resize(feSize, false); + // Create a bitfield which DOFs to interpolate, using the global bitVector + // Here also the counter might be incremented bool visit_any = false; for (std::size_t i = 0; i < feSize; ++i) { @@ -135,12 +167,15 @@ namespace Impl { LocalFunctionComponent localFj(node, treePath, localF_, nodeToRangeEntry_); thread_local std::vector interpolationCoefficients; + // Traverse the range-components of the coefficient vector std::size_t blockSize = Dune::Functions::flatVectorView(vector_[localView_.index(0)]).size(); for (std::size_t j = 0; j < blockSize; ++j) { localFj.setComponent(j); fe.localInterpolation().interpolate(localFj, interpolationCoefficients); assert(interpolationCoefficients.size() == feSize); + + // Traverse all local DOFs (only if marked for visit with the bitVector) for (std::size_t i = 0; i < feSize; ++i) { if (visit[i]) @@ -238,7 +273,7 @@ template const& treePath, Vec& vec, C& count, BV const& bitVec, GF const& gf, bool_t) { - auto ntrm = Dune::Functions::HierarchicNodeToRangeMap(); + auto ntrm = AMDiS::HierarchicNodeToRangeMap(); AMDiS::interpolateTreeSubset(basis, treePath, vec, count, bitVec, gf, ntrm, bool_t{}); } @@ -261,7 +296,7 @@ void interpolateTree(B const& basis, Dune::TypeTree::HybridTreePath const& Vec& vec, C& count, GF const& gf, bool_t) { auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector(); - auto ntrm = Dune::Functions::HierarchicNodeToRangeMap(); + auto ntrm = AMDiS::HierarchicNodeToRangeMap(); AMDiS::interpolateTreeSubset(basis, treePath, vec, count, bitVec, gf, ntrm, bool_t{}); } @@ -289,7 +324,7 @@ template const& treePath, Vec& vec, C& count, BV const& bitVec, GF const& gf, bool_t) { - auto ntrm = Dune::Functions::HierarchicNodeToRangeMap(); + auto ntrm = AMDiS::HierarchicNodeToRangeMap(); AMDiS::interpolateTreeSubset(basis, treePath, vec, count, bitVec, gf, ntrm, bool_t{}); } @@ -318,7 +353,7 @@ template ) { auto root = Dune::TypeTree::hybridTreePath(); - auto ntrm = Dune::Functions::HierarchicNodeToRangeMap(); + auto ntrm = AMDiS::HierarchicNodeToRangeMap(); AMDiS::interpolateTreeSubset(basis, root, vec, count, bitVec, gf, ntrm, bool_t{}); } -- GitLab From 2946bc64599f0f928abfc0ee6b21bacf2cc3d107 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Fri, 26 Apr 2019 10:10:20 +0200 Subject: [PATCH 4/4] use strategy-tag instead of bool parameter for interpolate --- src/amdis/functions/Interpolate.hpp | 6 +++--- src/amdis/gridfunctions/DOFVectorView.hpp | 19 +++++++++++++------ 2 files changed, 16 insertions(+), 9 deletions(-) diff --git a/src/amdis/functions/Interpolate.hpp b/src/amdis/functions/Interpolate.hpp index 09d2ed06..5845b52f 100644 --- a/src/amdis/functions/Interpolate.hpp +++ b/src/amdis/functions/Interpolate.hpp @@ -17,7 +17,7 @@ namespace AMDiS { namespace Impl { - struct DummyCounter + struct FakeCounter { struct Sink { @@ -392,7 +392,7 @@ void interpolate(B const& basis, Vec& vec, GF const& gf) { auto root = Dune::TypeTree::hybridTreePath(); auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector(); - auto count = Impl::DummyCounter(); + auto count = Impl::FakeCounter(); AMDiS::interpolateFiltered(basis, root, vec, count, bitVec, gf, std::false_type{}); } @@ -429,7 +429,7 @@ template const& treePath, Vec& vec, GF const& gf) { auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector(); - auto count = Impl::DummyCounter(); + auto count = Impl::FakeCounter(); AMDiS::interpolateFiltered(basis, treePath, vec, count, bitVec, gf, std::false_type{}); } diff --git a/src/amdis/gridfunctions/DOFVectorView.hpp b/src/amdis/gridfunctions/DOFVectorView.hpp index 722c9098..edbd33f9 100644 --- a/src/amdis/gridfunctions/DOFVectorView.hpp +++ b/src/amdis/gridfunctions/DOFVectorView.hpp @@ -6,6 +6,13 @@ namespace AMDiS { + namespace tag + { + struct average {}; + struct assign {}; + + } // end namespace tag + /// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction template class DOFVectorView @@ -40,15 +47,15 @@ namespace AMDiS * v.interpolate_noalias([](auto const& x) { return x[0]; }); * ``` **/ - template - void interpolate_noalias(Expr&& expr, bool_t = {}) + template + void interpolate_noalias(Expr&& expr, Tag strategy = {}) { auto const& basis = this->basis(); auto const& treePath = this->treePath(); auto&& gridFct = makeGridFunction(FWD(expr), basis.gridView()); - if (average) { + if (std::is_same::value) { thread_local std::vector counter; counter.clear(); AMDiS::interpolate(basis, treePath, coefficients(), counter, FWD(gridFct), std::true_type{}); @@ -73,13 +80,13 @@ namespace AMDiS * Allows to have a reference to the DOFVector in the expression, e.g. as * \ref DiscreteFunction or \ref gradientAtQP() of a DiscreteFunction. **/ - template - void interpolate(Expr&& expr, bool_t = {}) + template + void interpolate(Expr&& expr, Tag strategy = {}) { // create temporary copy of data DOFVector tmp(coefficients()); Self tmpView{tmp, this->treePath()}; - tmpView.interpolate_noalias(FWD(expr), bool_t{}); + tmpView.interpolate_noalias(FWD(expr), strategy); // move data from temporary vector into stored DOFVector coefficients().vector() = std::move(tmp.vector()); -- GitLab