diff --git a/src/amdis/CMakeLists.txt b/src/amdis/CMakeLists.txt index 706ce91d9481db4232090b5ad6295ae0645adc1c..3f3e58613886b786a246520ea80bce5971735c24 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/common/FieldMatVec.hpp b/src/amdis/common/FieldMatVec.hpp index dbc3f6c56db6970aeea80d219a82fe16b6cef91a..ceb547bf0ab9f54c7caed3b0aa9dc5eaaf4a133d 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 340b24e68fdc6c6e2aff7cafbb41beed8b92dd10..93f0c219159bc36224dce8e2ddd2957f0764a9b8 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/CMakeLists.txt b/src/amdis/functions/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..13244525b3280e8c2abc2f98606792edf8d03e57 --- /dev/null +++ b/src/amdis/functions/CMakeLists.txt @@ -0,0 +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 0000000000000000000000000000000000000000..259714400aafa064e6f958f92c249f3cac1c7b75 --- /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 new file mode 100644 index 0000000000000000000000000000000000000000..5845b52fe8a3966ec40058cc3a06a8a004422d8d --- /dev/null +++ b/src/amdis/functions/Interpolate.hpp @@ -0,0 +1,436 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace AMDiS { +namespace Impl { + + struct FakeCounter + { + 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; } + }; + + + /// \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 + { + 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; + + /// Functor called in the LocalInterpolation + 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_, transformTreePath(treePath_), tmp_vec))[comp_]; + } + + void setComponent(std::size_t comp) + { + 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_; + LocalFunction const& localF_; + NodeToRangeEntry const& nodeToRangeEntry_; + + std::size_t comp_ = 0; + }; + + 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) + : 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"); + } + + /// Apply the visitor to a node in the basis-tree (with corresponding treepath) + 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); + + // 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) + { + 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; + + // 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]) + { + 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 = AMDiS::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 = AMDiS::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 = AMDiS::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 = AMDiS::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::FakeCounter(); + 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::FakeCounter(); + 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 f95677a960a72583a1807c9e38efae57feee6a27..edbd33f96db49a9c082837875955cbedb053e966 100644 --- a/src/amdis/gridfunctions/DOFVectorView.hpp +++ b/src/amdis/gridfunctions/DOFVectorView.hpp @@ -1,10 +1,18 @@ #pragma once #include +#include #include 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 @@ -23,8 +31,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 +47,27 @@ namespace AMDiS * v.interpolate_noalias([](auto const& x) { return x[0]; }); * ``` **/ - template - void interpolate_noalias(Expr&& expr) + 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()); - Dune::Functions::interpolate(basis, treePath, coefficients(), FWD(gridFct)); + if (std::is_same::value) { + 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 +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) + 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)); + tmpView.interpolate_noalias(FWD(expr), strategy); // move data from temporary vector into stored DOFVector coefficients().vector() = std::move(tmp.vector()); @@ -107,6 +128,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)