From bbbca3bd2bf886cab7bf7ce7f7a2cd75e21472ab Mon Sep 17 00:00:00 2001
From: "Praetorius, Simon" <simon.praetorius@tu-dresden.de>
Date: Fri, 26 Apr 2019 15:51:45 +0200
Subject: [PATCH] reimplement interpolate function with averaging
---
src/amdis/CMakeLists.txt | 1 +
src/amdis/common/FieldMatVec.hpp | 17 +
src/amdis/common/FieldMatVec.inc.hpp | 54 +++
src/amdis/functions/CMakeLists.txt | 6 +
.../functions/HierarchicNodeToRangeMap.hpp | 41 ++
src/amdis/functions/Interpolate.hpp | 436 ++++++++++++++++++
src/amdis/gridfunctions/DOFVectorView.hpp | 42 +-
7 files changed, 589 insertions(+), 8 deletions(-)
create mode 100644 src/amdis/functions/CMakeLists.txt
create mode 100644 src/amdis/functions/HierarchicNodeToRangeMap.hpp
create mode 100644 src/amdis/functions/Interpolate.hpp
diff --git a/src/amdis/CMakeLists.txt b/src/amdis/CMakeLists.txt
index c21899d2..dd4aa640 100644
--- a/src/amdis/CMakeLists.txt
+++ b/src/amdis/CMakeLists.txt
@@ -67,6 +67,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 dbc3f6c5..ceb547bf 100644
--- a/src/amdis/common/FieldMatVec.hpp
+++ b/src/amdis/common/FieldMatVec.hpp
@@ -398,15 +398,32 @@ namespace Dune
template <class T1, class T2, int M, int N, int L>
FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_ABt(FieldMatrix<T1, M, L> const& A, FieldMatrix<T2, N, L> const& B);
+
template <class T1, class T2, class T3, int M, int N, int L>
FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, L> const& A, FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,M,N>& C);
+ template <class T1, class T2, class T3, int N, int L>
+ FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, L> const& A, FieldMatrix<T2, N, L> const& B, FieldVector<T3,N>& C);
+
+ template <class T1, class T2, class T3, int N, int L>
+ FieldVector<T3,N>& multiplies_ABt(FieldVector<T1, L> const& A, FieldMatrix<T2, N, L> const& B, FieldVector<T3,N>& C);
+
+ template <class T1, class T2, class T3, int N, int L>
+ FieldMatrix<T3,1,N>& multiplies_ABt(FieldVector<T1, L> const& A, FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,1,N>& C);
+
+
template <class T1, class T2, class T3, int M, int N>
FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, N> const& A, DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,M,N>& C);
template <class T1, class T2, class T3, int N>
FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C);
+ template <class T1, class T2, class T3, int N>
+ FieldVector<T3,N>& multiplies_ABt(FieldVector<T1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C);
+
+ template <class T1, class T2, class T3, int N>
+ FieldMatrix<T3,1,N>& multiplies_ABt(FieldVector<T1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,1,N>& C);
+
// -----------------------------------------------------------------------------
template <class T, int N>
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<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, L> const& A, FieldMatrix
return C;
}
+template <class T1, class T2, class T3, int N, int L>
+FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, L> const& A, FieldMatrix<T2, N, L> const& B, FieldVector<T3,N>& 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 <class T1, class T2, class T3, int N, int L>
+FieldVector<T3,N>& multiplies_ABt(FieldVector<T1, L> const& A, FieldMatrix<T2, N, L> const& B, FieldVector<T3,N>& 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 <class T1, class T2, class T3, int N, int L>
+FieldMatrix<T3,1,N>& multiplies_ABt(FieldVector<T1, L> const& A, FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,1,N>& 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 <class T1, class T2, class T3, int M, int N>
FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, N> const& A, DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,M,N>& C)
{
@@ -537,6 +572,25 @@ FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, N> const& A, DiagonalMatri
return C;
}
+template <class T1, class T2, class T3, int N>
+FieldVector<T3,N>& multiplies_ABt(FieldVector<T1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C)
+{
+ for (int n = 0; n < N; ++n) {
+ C[n] = A[n] * B.diagonal(n);
+ }
+ return C;
+}
+
+template <class T1, class T2, class T3, int N>
+FieldMatrix<T3,1,N>& multiplies_ABt(FieldVector<T1, N> const& A, DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,1,N>& C)
+{
+ for (int n = 0; n < N; ++n) {
+ C[0][n] = A[n] * B.diagonal(n);
+ }
+ return C;
+}
+
+
template <class T, int N>
T const& at(FieldMatrix<T,N,1> 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 00000000..13244525
--- /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 00000000..25971440
--- /dev/null
+++ b/src/amdis/functions/HierarchicNodeToRangeMap.hpp
@@ -0,0 +1,41 @@
+#pragma once
+
+#include <utility>
+#include <type_traits>
+
+#include <dune/common/concept.hh>
+
+#include <dune/functions/functionspacebases/concepts.hh>
+#include <dune/functions/common/indexaccess.hh>
+
+#include <amdis/common/ConceptsBase.hpp>
+
+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 <class Node, class TreePath, class Range,
+ REQUIRES(Dune::models<Dune::Functions::Concept::HasIndexAccess, Range, Dune::index_constant<0>>())>
+ decltype(auto) operator()(const Node& node, const TreePath& treePath, Range&& y) const
+ {
+ return Dune::Functions::resolveStaticMultiIndex(y, treePath);
+ }
+
+ template <class Node, class TreePath, class Range,
+ REQUIRES(not Dune::models<Dune::Functions::Concept::HasIndexAccess, Range, Dune::index_constant<0>>())>
+ decltype(auto) operator()(const Node& node, const TreePath& treePath, Range&& y) const
+ {
+ return std::forward<Range>(y);
+ }
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/functions/Interpolate.hpp b/src/amdis/functions/Interpolate.hpp
new file mode 100644
index 00000000..5845b52f
--- /dev/null
+++ b/src/amdis/functions/Interpolate.hpp
@@ -0,0 +1,436 @@
+#pragma once
+
+#include <memory>
+#include <vector>
+
+#include <dune/common/tuplevector.hh>
+#include <dune/common/version.hh>
+#include <dune/common/std/optional.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+
+#include <amdis/common/Concepts.hpp>
+#include <amdis/common/FieldMatVec.hpp>
+#include <amdis/common/Logical.hpp>
+#include <amdis/functions/HierarchicNodeToRangeMap.hpp>
+#include <amdis/typetree/Traversal.hpp>
+
+namespace AMDiS {
+namespace Impl {
+
+ struct FakeCounter
+ {
+ struct Sink
+ {
+ template <class T> constexpr Sink& operator =(T const&) { return *this; }
+ template <class T> constexpr Sink& operator+=(T const&) { return *this; }
+ template <class T> constexpr Sink& operator-=(T const&) { return *this; }
+
+ constexpr operator int() const { return 1; }
+ };
+
+ template <class SizeInfo>
+ constexpr void resize(SizeInfo const& sizeInfo) { /* do nothing */ }
+ constexpr std::size_t size() const { return 0; }
+
+ template <class MI> constexpr Sink operator[](MI const&) { return Sink{}; }
+ template <class MI> 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 B, class Vec, class C, class BV, class LF, class NTRE, bool average>
+ 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 Node, class TreePath>
+ class LocalFunctionComponent
+ : public Dune::LocalFiniteElementFunctionBase<typename Node::FiniteElement>::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<LocalFunction, LocalDomain>,
+ "Function passed to LocalInterpolateVisitor does not model the Callable<LocalCoordinate> concept");
+ }
+
+ /// Apply the visitor to a node in the basis-tree (with corresponding treepath)
+ template <class Node, class TreePath>
+ 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<bool> 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<Node, TreePath> localFj(node, treePath, localF_, nodeToRangeEntry_);
+ thread_local std::vector<RangeField> 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 <class B, class Vec>
+ decltype(auto) toVectorBackend(B const& basis, Vec& vec)
+ {
+ return Dune::Hybrid::ifElse(Dune::models<Dune::Functions::Concept::VectorBackend<B>, Vec>(),
+ [&](auto id) -> decltype(auto) { return vec; },
+ [&](auto id) -> decltype(auto) { return Dune::Functions::istlVectorBackend(vec); });
+ }
+
+ template <class B, class Vec>
+ decltype(auto) toConstVectorBackend(B const& basis, Vec const& vec)
+ {
+ return Dune::Hybrid::ifElse(Dune::models<Dune::Functions::Concept::ConstVectorBackend<B>, 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 <class B, class TP, class Vec, class C, class BV, class GF, class NTRE, bool average>
+void interpolateTreeSubset(B const& basis, TP const& treePath, Vec& vec, C& count, BV const& bitVec,
+ GF const& gf, NTRE const& nodeToRangeEntry, bool_t<average>)
+{
+ 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<B, TYPEOF(vector), TYPEOF(counter), TYPEOF(bitVector), TYPEOF(lf), NTRE, average>;
+ for_each_leaf_node(subTree, Visitor{vector, counter, bitVector, lf, localView, nodeToRangeEntry});
+ }
+}
+
+
+template <class B, class... I, class Vec, class C, class BV, class GF, bool average,
+ REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename B::GridView>,B>()),
+ REQUIRES(Concepts::GridFunction<GF>)>
+void interpolateTreeSubset(B const& basis, Dune::TypeTree::HybridTreePath<I...> const& treePath,
+ Vec& vec, C& count, BV const& bitVec, GF const& gf, bool_t<average>)
+{
+ auto ntrm = AMDiS::HierarchicNodeToRangeMap();
+ AMDiS::interpolateTreeSubset(basis, treePath, vec, count, bitVec, gf, ntrm, bool_t<average>{});
+}
+
+
+template <class B, class... I, class Vec, class C, class GF, class NTRE, bool average,
+ REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename B::GridView>,B>()),
+ REQUIRES(Concepts::GridFunction<GF>)>
+void interpolateTree(B const& basis, Dune::TypeTree::HybridTreePath<I...> const& treePath,
+ Vec& vec, C& count, GF const& gf, NTRE const& nodeToRangeEntry, bool_t<average>)
+{
+ auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector();
+ AMDiS::interpolateTreeSubset(basis, treePath, vec, count, bitVec, gf, nodeToRangeEntry, bool_t<average>{});
+}
+
+
+template <class B, class... I, class Vec, class C, class GF, bool average,
+ REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename B::GridView>,B>()),
+ REQUIRES(Concepts::GridFunction<GF>)>
+void interpolateTree(B const& basis, Dune::TypeTree::HybridTreePath<I...> const& treePath,
+ Vec& vec, C& count, GF const& gf, bool_t<average>)
+{
+ auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector();
+ auto ntrm = AMDiS::HierarchicNodeToRangeMap();
+ AMDiS::interpolateTreeSubset(basis, treePath, vec, count, bitVec, gf, ntrm, bool_t<average>{});
+}
+
+
+/**
+ * \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 <class B, class... I, class Vec, class C, class GF, class BV, bool average,
+ REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename B::GridView>,B>()),
+ REQUIRES(Concepts::GridFunction<GF>)>
+void interpolateFiltered(B const& basis, Dune::TypeTree::HybridTreePath<I...> const& treePath,
+ Vec& vec, C& count, BV const& bitVec, GF const& gf, bool_t<average>)
+{
+ auto ntrm = AMDiS::HierarchicNodeToRangeMap();
+ AMDiS::interpolateTreeSubset(basis, treePath, vec, count, bitVec, gf, ntrm, bool_t<average>{});
+}
+
+
+/**
+ * \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 <class B, class Vec, class C, class BV, class GF, bool average,
+ REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename B::GridView>,B>()),
+ REQUIRES(not Dune::Functions::Imp::isHybridTreePath<Vec>()),
+ REQUIRES(Concepts::GridFunction<GF>)>
+void interpolateFiltered(B const& basis, Vec& vec, C& count, BV const& bitVec, GF const& gf, bool_t<average>)
+{
+ auto root = Dune::TypeTree::hybridTreePath();
+ auto ntrm = AMDiS::HierarchicNodeToRangeMap();
+ AMDiS::interpolateTreeSubset(basis, root, vec, count, bitVec, gf, ntrm, bool_t<average>{});
+}
+
+
+/**
+ * \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 <class B, class Vec, class C, class GF, bool average,
+ REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename B::GridView>,B>()),
+ REQUIRES(not Dune::Functions::Imp::isHybridTreePath<Vec>()),
+ REQUIRES(Concepts::GridFunction<GF>)>
+void interpolate(B const& basis, Vec& vec, C& count, GF const& gf, bool_t<average>)
+{
+ auto root = Dune::TypeTree::hybridTreePath();
+ auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector();
+ AMDiS::interpolateFiltered(basis, root, vec, count, bitVec, gf, bool_t<average>{});
+}
+
+template <class B, class... I, class Vec, class GF,
+ REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename B::GridView>,B>()),
+ REQUIRES(not Dune::Functions::Imp::isHybridTreePath<Vec>()),
+ REQUIRES(Concepts::GridFunction<GF>)>
+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 <class B, class... I, class Vec, class C, class GF, bool average,
+ REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename B::GridView>,B>()),
+ REQUIRES(Concepts::GridFunction<GF>)>
+void interpolate(B const& basis, Dune::TypeTree::HybridTreePath<I...> const& treePath,
+ Vec& vec, C& count, GF const& gf, bool_t<average>)
+{
+ auto bitVec = Dune::Functions::Imp::AllTrueBitSetVector();
+ AMDiS::interpolateFiltered(basis, treePath, vec, count, bitVec, gf, bool_t<average>{});
+}
+
+template <class B, class... I, class Vec, class GF,
+ REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename B::GridView>,B>()),
+ REQUIRES(Concepts::GridFunction<GF>)>
+void interpolate(B const& basis, Dune::TypeTree::HybridTreePath<I...> 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 f95677a9..edbd33f9 100644
--- a/src/amdis/gridfunctions/DOFVectorView.hpp
+++ b/src/amdis/gridfunctions/DOFVectorView.hpp
@@ -1,10 +1,18 @@
#pragma once
#include <amdis/GridFunctions.hpp>
+#include <amdis/functions/Interpolate.hpp>
#include <amdis/gridfunctions/DiscreteFunction.hpp>
namespace AMDiS
{
+ namespace tag
+ {
+ struct average {};
+ struct assign {};
+
+ } // end namespace tag
+
/// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
template <class GB, class VT, class TP>
class DOFVectorView
@@ -23,8 +31,9 @@ namespace AMDiS
{}
/// Constructor. Stores a pointer to the mutable `dofvector`.
- DOFVectorView(DOFVector<GB,VT>& dofVector, TP const& treePath)
- : Super(dofVector, treePath)
+ template <class PreTreePath>
+ DOFVectorView(DOFVector<GB,VT>& 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 <class Expr>
- void interpolate_noalias(Expr&& expr)
+ template <class Expr, class Tag = tag::average>
+ 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<Tag, tag::average>::value) {
+ thread_local std::vector<std::uint8_t> 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 <class Expr>
- void interpolate(Expr&& expr)
+ template <class Expr, class Tag = tag::average>
+ void interpolate(Expr&& expr, Tag strategy = {})
{
// create temporary copy of data
DOFVector<GB,VT> 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 <class GlobalBasis, class ValueType, class PreTreePath>
+ DOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector, PreTreePath const& preTreePath)
+ -> DOFVectorView<GlobalBasis, ValueType, TYPEOF(makeTreePath(preTreePath))>;
+
// Deduction guide for DOFVectorView class
template <class GlobalBasis, class ValueType>
DOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector)
--
GitLab