Commit bbbca3bd authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

reimplement interpolate function with averaging

parent e9859f1d
......@@ -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")
......
......@@ -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>
......
......@@ -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)
{
......
#install headers
install(FILES
HierarchicNodeToRangeMap.hpp
Interpolate.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/functions)
#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
#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>{});
}