Commit 2d132a0e authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Allow std containers as coefficients in DiscreteFunction and allow DefaultGlobalBasis

add unified access for the implementation of DiscreteFunction

Improved test and range initialization of DiscreteFunction

extend the DiscreteFunctionTest to test for construction with std containers

Add range type to solution() and oldSolution() functions
parent f9a94a4e
......@@ -62,12 +62,16 @@ namespace AMDiS
}
/// Return a const view to a oldSolution component
template <class... Indices>
/**
* \tparam Range The range type return by evaluating the view in coordinates. If not specified,
* it is automatically selected using \ref RangeType_t template.
**/
template <class Range = void, class... Indices>
auto oldSolution(Indices... ii) const
{
test_exit_dbg(bool(oldSolution_),
"OldSolution need to be created. Call initialize with INIT_UH_OLD.");
return valueOf(*oldSolution_, ii...);
return valueOf<Range>(*oldSolution_, ii...);
}
/// Implementation of \ref ProblemTimeInterface::transferInitialSolution().
......
......@@ -393,19 +393,27 @@ namespace AMDiS
/// Return a mutable view to a solution component
template <class... Indices>
/**
* \tparam Range The range type return by evaluating the view in coordinates. If not specified,
* it is automatically selected using \ref RangeType_t template.
**/
template <class Range = void, class... Indices>
auto solution(Indices... ii)
{
assert(bool(solution_) && "You have to call initialize() before.");
return valueOf(*solution_, ii...);
return valueOf<Range>(*solution_, ii...);
}
/// Return a const view to a solution component
template <class... Indices>
/**
* \tparam Range The range type return by evaluating the view in coordinates. If not specified,
* it is automatically selected using \ref RangeType_t template.
**/
template <class Range = void, class... Indices>
auto solution(Indices... ii) const
{
assert(bool(solution_) && "You have to call initialize() before.");
return valueOf(*solution_, ii...);
return valueOf<Range>(*solution_, ii...);
}
......
......@@ -3,14 +3,16 @@
#include <optional>
#include <vector>
#include <dune/common/ftraits.hh>
#include <dune/common/typeutilities.hh>
#include <dune/functions/common/defaultderivativetraits.hh>
#include <dune/functions/functionspacebases/defaultnodetorangemap.hh>
#include <dune/functions/functionspacebases/flatvectorview.hh>
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <dune/typetree/childextraction.hh>
#include <amdis/common/Tags.hpp>
#include <amdis/functions/NodeCache.hpp>
#include <amdis/linearalgebra/Access.hpp>
#include <amdis/typetree/FiniteElementType.hpp>
#include <amdis/typetree/RangeType.hpp>
#include <amdis/typetree/TreePath.hpp>
......@@ -22,28 +24,28 @@ namespace AMDiS
/**
* \ingroup GridFunctions
*
* \tparam Coeff Const or mutable coefficient type of the DOFVector
* \tparam GB Thy type of the global basis
* \tparam Coeff Const or mutable coefficient vector
* \tparam GB The type of the global basis
* \tparam TreePath A realization of \ref Dune::TypeTree::HybridTreePath
* \tparam Range The range type for th evaluation of the discrete function
*
* **Requirements:**
* - GB models \ref Concepts::GlobalBasis
**/
template <class Coeff, class GB, class TreePath>
template <class Coeff, class GB, class TreePath, class Range = void>
class DiscreteFunction;
/// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
template <class Coeff, class GB, class TreePath>
template <class Coeff, class GB, class TreePath, class R>
class DiscreteFunction
: public DiscreteFunction<Coeff const,GB,TreePath>
: public DiscreteFunction<Coeff const,GB,TreePath,R>
{
using Self = DiscreteFunction<Coeff,GB,TreePath>;
using Super = DiscreteFunction<Coeff const,GB,TreePath>;
using Self = DiscreteFunction<Coeff,GB,TreePath,R>;
using Super = DiscreteFunction<Coeff const,GB,TreePath,R>;
using Coefficients = Coeff;
using GlobalBasis = GB;
using ValueType = typename Coeff::value_type;
public:
/// Constructor. Stores a pointer to the mutable `dofvector`.
......@@ -54,16 +56,17 @@ namespace AMDiS
{}
template <class... Path>
DiscreteFunction(Coefficients& dofVector, std::shared_ptr<GlobalBasis const> const& basis, Path... path)
: DiscreteFunction(dofVector, *basis, path...)
DiscreteFunction(Coefficients& coefficients, std::shared_ptr<GlobalBasis const> const& basis, Path... path)
: DiscreteFunction(coefficients, *basis, path...)
{}
/// Construct a DiscreteFunction directly from a DOFVector
template <class DV, class... Path,
Dune::disableCopyMove<DiscreteFunction,DV> = 0,
class Coeff_ = TYPEOF(std::declval<DV>().coefficients()),
class GB_ = TYPEOF(*std::declval<DV>().basis())>
DiscreteFunction(DV&& dofVector, Path... path)
: DiscreteFunction(dofVector.coefficients(), *dofVector.basis(), path...)
class GB_ = TYPEOF(std::declval<DV>().basis())>
explicit DiscreteFunction(DV&& dofVector, Path... path)
: DiscreteFunction(dofVector.coefficients(), dofVector.basis(), path...)
{}
public:
......@@ -125,11 +128,11 @@ namespace AMDiS
/// Return the const DOFVector
using Super::coefficients;
template <class... Indices>
template <class Range = void, class... Indices>
auto child(Indices... ii)
{
auto tp = cat(this->treePath_, makeTreePath(ii...));
return DiscreteFunction<Coeff, GB, TYPEOF(makeTreePath(tp))>{*mutableCoeff_, this->basis(), tp};
return DiscreteFunction<Coeff, GB, TYPEOF(makeTreePath(tp)), Range>{*mutableCoeff_, this->basis(), tp};
}
using Super::child;
......@@ -139,22 +142,20 @@ namespace AMDiS
};
/// A Const DiscreteFunction
template <class Coeff, class GB, class TreePath>
class DiscreteFunction<Coeff const,GB,TreePath>
template <class Coeff, class GB, class TreePath, class R>
class DiscreteFunction<Coeff const,GB,TreePath,R>
{
private:
using Coefficients = std::remove_const_t<Coeff>;
using GlobalBasis = GB;
using LocalView = typename GlobalBasis::LocalView;
using GridView = typename GlobalBasis::GridView;
using ValueType = typename Coefficients::value_type;
using Coefficient = TYPEOF(std::declval<Traits::Access>()(std::declval<Coeff const>(), std::declval<typename GB::MultiIndex>()));
using Tree = typename GlobalBasis::LocalView::Tree;
using SubTree = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
using TreeCache = NodeCache_t<Tree>;
using SubTreeCache = typename Dune::TypeTree::ChildForTreePath<TreeCache, TreePath>;
public:
/// Set of entities the DiscreteFunction is defined on
using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
......@@ -162,14 +163,14 @@ namespace AMDiS
/// Global coordinates of the EntitySet
using Domain = typename EntitySet::GlobalCoordinate;
/// Range type of this DiscreteFunction
using Range = RangeType_t<SubTree,ValueType>;
/// Range type of this DiscreteFunction. If R=void deduce the Range automatically, using RangeType_t
using Range = std::conditional_t<std::is_same_v<R,void>, RangeType_t<SubTree,Coefficient>, R>;
/// \brief This GridFunction has no derivative function, it can be created
/// by \ref DiscreteGridFunction.
enum { hasDerivative = false };
public:
private:
/// A LocalFunction representing the localfunction and derivative of the DOFVector on a bound element
template <class Type>
class LocalFunction;
......@@ -191,9 +192,10 @@ namespace AMDiS
/// Construct a DiscreteFunction directly from a DOFVector
template <class DV, class... Path,
Dune::disableCopyMove<DiscreteFunction,DV> = 0,
class Coeff_ = TYPEOF(std::declval<DV>().coefficients()),
class GB_ = TYPEOF(std::declval<DV>().basis())>
DiscreteFunction(DV const& dofVector, Path... path)
explicit DiscreteFunction(DV const& dofVector, Path... path)
: DiscreteFunction(dofVector.coefficients(), dofVector.basis(), path...)
{}
......@@ -230,11 +232,11 @@ namespace AMDiS
return *coefficients_;
}
template <class... Indices>
template <class Range = void, class... Indices>
auto child(Indices... ii) const
{
auto tp = cat(this->treePath_, makeTreePath(ii...));
return DiscreteFunction<Coeff const, GB, TYPEOF(makeTreePath(tp))>{*coefficients_, *basis_, tp};
return DiscreteFunction<Coeff const, GB, TYPEOF(makeTreePath(tp)), Range>{*coefficients_, *basis_, tp};
}
protected:
......@@ -247,45 +249,47 @@ namespace AMDiS
// deduction guides
template <class Coeff, class Basis, class... Path,
class TP = TYPEOF(makeTreePath(std::declval<Path>()...)),
template <class C, class Basis, class... Indices,
class GB = Underlying_t<Basis>,
class TP = TYPEOF(makeTreePath(std::declval<Indices>()...)),
REQUIRES(Concepts::GlobalBasis<GB>)>
DiscreteFunction(Coeff&, Basis const&, Path...)
-> DiscreteFunction<Coeff,GB,TP>;
DiscreteFunction(C&, Basis const&, Indices...)
-> DiscreteFunction<C,GB,TP>;
template <class DV, class... Path,
class Coeff = decltype(std::declval<DV>().coefficients()),
class GB = decltype(std::declval<DV>().basis()),
class TP = TYPEOF(makeTreePath(std::declval<Path>()...))>
DiscreteFunction(DV&, Path...)
-> DiscreteFunction<std::remove_reference_t<Coeff>,Underlying_t<GB>,TP>;
template <class DV, class... Indices,
class C = decltype(std::declval<DV>().coefficients()),
class GB = decltype(std::declval<DV>().basis()),
class TP = TYPEOF(makeTreePath(std::declval<Indices>()...))>
DiscreteFunction(DV&, Indices...)
-> DiscreteFunction<std::remove_reference_t<C>,Underlying_t<GB>,TP>;
// grid functions representing the DOFVector
// -----------------------------------------
/// A Generator for the childs of a mutable \ref DiscreteFunction
template <class Coeff, class GB, class... Path, class... Indices>
auto valueOf(DiscreteFunction<Coeff,GB,Path...>& df, Indices... ii)
template <class Range = void, class C, class GB, class TP, class R, class... Indices>
auto valueOf(DiscreteFunction<C,GB,TP,R>& df, Indices... ii)
{
return df.child(ii...);
return df.template child<Range>(ii...);
}
/// A Generator for the childs of a const \ref DiscreteFunction
template <class Coeff, class GB, class... Path, class... Indices>
auto valueOf(DiscreteFunction<Coeff,GB,Path...> const& df, Indices... ii)
template <class Range = void, class C, class GB, class TP, class R, class... Indices>
auto valueOf(DiscreteFunction<C,GB,TP,R> const& df, Indices... ii)
{
return df.child(ii...);
return df.template child<Range>(ii...);
}
/// A Generator to transform a DOFVector into a \ref DiscreteFunction
template <class DV, class... Indices,
class = decltype(std::declval<DV>().coefficients()),
class = decltype(std::declval<DV>().basis())>
template <class Range = void, class DV, class... Indices,
class C = decltype(std::declval<DV>().coefficients()),
class GB = decltype(std::declval<DV>().basis()),
class TP = TYPEOF(makeTreePath(std::declval<Indices>()...))>
auto valueOf(DV& dofVec, Indices... ii)
{
return DiscreteFunction{dofVec.coefficients(), dofVec.basis(), makeTreePath(ii...)};
using DF = DiscreteFunction<std::remove_reference_t<C>,Underlying_t<GB>,TP,Range>;
return DF{dofVec.coefficients(), dofVec.basis(), makeTreePath(ii...)};
}
} // end namespace AMDiS
......
......@@ -11,8 +11,8 @@
namespace AMDiS {
// Evaluate DiscreteFunction in global coordinates
template <class Coeff, class GB, class TP>
typename DiscreteFunction<Coeff const,GB,TP>::Range DiscreteFunction<Coeff const,GB,TP>::
template <class C, class GB, class TP, class R>
typename DiscreteFunction<C const,GB,TP,R>::Range DiscreteFunction<C const,GB,TP,R>::
operator()(Domain const& x) const
{
using Grid = typename GlobalBasis::GridView::Grid;
......@@ -30,9 +30,9 @@ typename DiscreteFunction<Coeff const,GB,TP>::Range DiscreteFunction<Coeff const
// Interpolation of GridFunction to DOFVector
template <class Coeff, class GB, class TP>
template <class C, class GB, class TP, class R>
template <class Expr, class Tag>
void DiscreteFunction<Coeff,GB,TP>::
void DiscreteFunction<C,GB,TP,R>::
interpolate_noalias(Expr&& expr, Tag strategy)
{
auto const& basis = this->basis();
......@@ -41,7 +41,7 @@ void DiscreteFunction<Coeff,GB,TP>::
auto&& gf = makeGridFunction(FWD(expr), basis.gridView());
if (std::is_same_v<Tag, tag::average>) {
VectorType_t<short,Coeff> counter(basis);
VectorType_t<short,Coefficients> counter(basis);
AMDiS::interpolate(basis, coefficients(), gf, treePath, counter);
coefficients().forEach([&counter](auto dof, auto& coeff)
......@@ -55,13 +55,13 @@ void DiscreteFunction<Coeff,GB,TP>::
// Interpolation of GridFunction to DOFVector
template <class Coeff, class GB, class TP>
template <class C, class GB, class TP, class R>
template <class Expr, class Tag>
void DiscreteFunction<Coeff,GB,TP>::
void DiscreteFunction<C,GB,TP,R>::
interpolate(Expr&& expr, Tag strategy)
{
// create temporary copy of data
Coeff tmp(coefficients());
Coefficients tmp(coefficients());
Self tmpView{tmp, this->basis(), this->treePath()};
tmpView.interpolate_noalias(FWD(expr), strategy);
......
......@@ -4,6 +4,8 @@
#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/common/FieldMatVec.hpp>
#include <amdis/common/Math.hpp>
#include <amdis/common/RecursiveForEach.hpp>
#include <amdis/functions/HierarchicNodeToRangeMap.hpp>
#include <amdis/functions/NodeIndices.hpp>
#include <amdis/functions/Order.hpp>
#include <amdis/typetree/FiniteElementType.hpp>
......@@ -14,7 +16,7 @@
namespace AMDiS {
namespace Impl {
// specialization of Coeff has gather method
// specialization for containers with gather method
template <class Coeff, class LocalView, class LocalCoeff,
class = decltype(std::declval<Coeff>().gather(std::declval<LocalView>(), std::declval<LocalCoeff&>()))>
void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoeff, Dune::PriorityTag<2>)
......@@ -22,8 +24,9 @@ void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoe
coeff.gather(localView, localCoeff);
}
// fallback implementation
template <class Coeff, class LocalView, class LocalCoeff>
// implementation for containers with multi-index access
template <class Coeff, class LocalView, class LocalCoeff,
class = decltype(std::declval<Coeff>()[std::declval<typename LocalView::MultiIndex>()])>
void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoeff, Dune::PriorityTag<1>)
{
localCoeff.resize(localView.size());
......@@ -32,12 +35,24 @@ void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoe
*it++ = coeff[idx];
}
// fallback implementation
template <class Coeff, class LocalView, class LocalCoeff,
class = decltype(std::declval<Coeff>()[std::integral_constant<std::size_t,0>{}])>
void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoeff, Dune::PriorityTag<0>)
{
auto backend = Dune::Functions::istlVectorBackend(coeff);
localCoeff.resize(localView.size());
auto it = localCoeff.begin();
for (auto const& idx : nodeIndices(localView))
*it++ = backend[idx];
}
} // end namespace Impl
template <class Coeff, class GB, class TP>
template <class C, class GB, class TP, class R>
template <class Type>
class DiscreteFunction<Coeff const,GB,TP>::LocalFunction
class DiscreteFunction<C const,GB,TP,R>::LocalFunction
{
using DomainType = typename DiscreteFunction::Domain;
using RangeType = typename DiscreteFunction::Range;
......@@ -55,17 +70,16 @@ private:
using LocalView = typename GlobalBasis::LocalView;
using Element = typename EntitySet::Element;
using Geometry = typename Element::Geometry;
using NodeToRangeMap = Dune::Functions::DefaultNodeToRangeMap<SubTree>;
using NodeToRangeMap = HierarchicNodeToRangeMap;
public:
/// Constructor. Stores a copy of the DiscreteFunction.
template <class LV>
LocalFunction(std::shared_ptr<LV> localView, TP const& treePath, Coeff const& coefficients, Type type)
LocalFunction(std::shared_ptr<LV> localView, TP const& treePath, C const& coefficients, Type type)
: localView_(std::move(localView))
, treePath_(treePath)
, coefficients_(coefficients)
, type_(type)
, nodeToRangeMap_(subTree())
{}
/// \brief Bind the LocalView to the element
......@@ -103,7 +117,7 @@ public:
return evaluateDivergence(local);
}
return Range(0);
return Range{};
}
/// \brief Create a LocalFunction representing the derivative.
......@@ -119,7 +133,7 @@ public:
-> decltype(AMDiS::order(std::declval<SubTree>()))
{
assert(bound_);
return AMDiS::order(subTreeCache());
return AMDiS::order(subTree());
}
/// \brief Return the bound element
......@@ -143,42 +157,62 @@ private:
template <class LocalCoordinate>
auto evaluatePartialDerivative (const LocalCoordinate& local) const;
// get the geometry of the bound element, that is constructed in the bind() method
Geometry const& geometry() const
{
assert(bound_);
return *geometry_;
}
// return the sub-tree of the basis-tree associated to the tree-path
SubTree const& subTree() const
{
return Dune::TypeTree::child(localView_->tree(), treePath_);
}
SubTreeCache const& subTreeCache() const
// return the sub-tree of the basis-tree-cache associated to the tree-path
// NOTE: this is only return in case the local-view does provide a treeCache() function
template <class LV,
class = decltype(std::declval<LV>().treeCache())>
auto const& subTreeCache(Dune::PriorityTag<1>) const
{
return Dune::TypeTree::child(localView_->treeCache(), treePath_);
}
// construct a new tree-cache that stores a reference to the sub-tree
template <class>
auto subTreeCache(Dune::PriorityTag<0>) const
{
return makeNodeCache(subTree());
}
decltype(auto) subTreeCache() const
{
return subTreeCache<LocalView>(Dune::PriorityTag<2>());
}
private:
std::shared_ptr<LocalView> localView_;
TP treePath_;
Coeff const& coefficients_;
Coefficients const& coefficients_;
Type type_;
NodeToRangeMap nodeToRangeMap_;
std::optional<Geometry> geometry_;
std::vector<ValueType> localCoefficients_;
std::vector<Coefficient> localCoefficients_;
bool bound_ = false;
};
template <class Coeff, class GB, class TP>
template <class C, class GB, class TP, class R>
template <class Type>
template <class LocalCoordinate>
auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
::evaluateFunction(LocalCoordinate const& local) const
{
Range y(0);
Range y;
Recursive::forEach(y, [](auto& v) { v = 0; });
Traversal::forEachLeafNode(this->subTreeCache(), [&](auto const& node, auto const& tp)
{
auto const& shapeFunctionValues = node.localBasisValuesAt(local);
......@@ -209,13 +243,15 @@ auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
}
template <class Coeff, class GB, class TP>
template <class C, class GB, class TP, class R>
template <class Type>
template <class LocalCoordinate>
auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
::evaluateJacobian(LocalCoordinate const& local) const
{
Range dy(0);
Range dy;
Recursive::forEach(dy, [](auto& v) { v = 0; });
Traversal::forEachLeafNode(this->subTreeCache(), [&](auto const& node, auto const& tp)
{
LocalToGlobalBasisAdapter localBasis(node, this->geometry());
......@@ -246,15 +282,16 @@ auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
}
template <class Coeff, class GB, class TP>
template <class C, class GB, class TP, class R>
template <class Type>
template <class LocalCoordinate>
auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
::evaluateDivergence(LocalCoordinate const& local) const
{
static_assert(static_size_v<Range> == 1, "Divergence implemented for scalar coefficients only.");
Range dy(0);
Range dy;
Recursive::forEach(dy, [](auto& v) { v = 0; });
auto&& node = this->subTreeCache();
......@@ -275,15 +312,17 @@ auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
}
template <class Coeff, class GB, class TP>
template <class C, class GB, class TP, class R>
template <class Type>
template <class LocalCoordinate>
auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
::evaluatePartialDerivative(LocalCoordinate const& local) const
{
std::size_t comp = this->type_.comp;
Range dy(0);
Range dy;
Recursive::forEach(dy, [](auto& v) { v = 0; });
Traversal::forEachLeafNode(this->subTreeCache(), [&](auto const& node, auto const& tp)
{
LocalToGlobalBasisAdapter localBasis(node, this->geometry());
......
#pragma once
#include <dune/common/typeutilities.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
namespace AMDiS {
namespace Traits {
/// Unified accessor class
struct Access
{
private:
// if C implements the vector facade interface
template <class C, class MI>
auto operator()(C const& c, MI const& mi, Dune::PriorityTag<2>) const
-> decltype(c.at(mi))
{
return c.at(mi);
}
// if C supports multi-index access
template <class C, class MI>
auto operator()(C const& c, MI const& mi, Dune::PriorityTag<1>) const
-> decltype(c[mi])
{
return c[mi];
}
// fallback implementation for simple vectors
template <class C, class MI>
auto operator()(C const& c, MI const& mi, Dune::PriorityTag<0>) const
-> decltype(Dune::Functions::istlVectorBackend(c)[mi])
{
return Dune::Functions::istlVectorBackend(c)[mi];
}
public:
template <class C, class MI>
auto operator()(C const& c, MI const& mi) const