Skip to content
Snippets Groups Projects
Commit 05e93b66 authored by graeser's avatar graeser
Browse files

Simplify internal of form mechanism significantly 

While the new version is slighly longer, it uses far less magic
and is much easier to understand.

(E.g. no empty `std::tuple<>` return values, no methods that peel
out the argument they are respondible for from a passed `tuple`, ...).
parent 74685ddc
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
dune/fufem-forms/forms.hh
+ 331
− 234
View file @ 05e93b66
......@@ -7,6 +7,8 @@
#include <tuple>
#include <utility>
#include <dune/common/typetraits.hh>
#include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>
......@@ -24,14 +26,14 @@ namespace Dune::Fufem::Forms {
namespace Impl {
template<class... T, class F, std::size_t... i>
void forEachTupleEntryHelper(std::tuple<T...>& tuple, F&& f, std::index_sequence<i...>){
template<class Tuple, class F, std::size_t... i>
void forEachTupleEntryHelper(Tuple&& tuple, F&& f, std::index_sequence<i...>){
(f(std::get<i>(tuple)),...);
}
template<class... T, class F>
void forEachTupleEntry(std::tuple<T...>& tuple, F&& f){
forEachTupleEntryHelper(tuple, f, std::make_index_sequence<sizeof...(T)>{});
template<class Tuple, class F>
void forEachTupleEntry(Tuple&& tuple, F&& f){
forEachTupleEntryHelper(tuple, f, std::make_index_sequence<std::tuple_size_v<std::decay_t<Tuple>>>{});
}
struct TransposeOp {
......@@ -176,9 +178,6 @@ namespace Impl {
template<class F>
using FormRange = decltype(std::declval<LazyFormRange<F>>()(std::tuple<int,int>{}));
template<class P, class T1, class T2>
using ProductConcept = std::void_t<decltype(std::declval<P>()(std::declval<T1>(),std::declval<T2>()))>;
} // namespace Impl
......@@ -210,12 +209,29 @@ namespace Impl {
// results in a k-linear form \int_D F(...)dx : V_1 \times ... \times V_k -> D
class Form {
public:
const auto& quadratureRuleKey() const
{
return quadKey_;
}
protected:
QuadratureRuleKey quadKey_;
template<class Element>
void bindToElement(const Element&) {}
template<class LocalView>
void bindToTestLocalView(const LocalView&) {}
template<class LocalView>
void bindToTrialLocalView(const LocalView&) {}
template<class Cache>
void bindToTestCache(const Cache&) {}
template<class Cache>
void bindToTrialCache(const Cache&) {}
void testBasis() const {}
void trialBasis() const {}
void testTreePath() const {}
void trialTreePath() const {}
};
......@@ -229,6 +245,69 @@ namespace Impl {
inline constexpr bool isForm_v = isForm<F>::value;
template<class Impl, std::size_t argIndex>
struct UnaryForm;
template<class Impl>
struct UnaryForm<Impl, 0> : public Form
{
static constexpr std::size_t arity = 1;
template<class LocalView>
void bindToTestLocalView(const LocalView& localView)
{
static_cast<Impl*>(this)->bindToLocalView(localView);
}
template<class Cache>
void bindToTestCache(Cache& cache)
{
static_cast<Impl*>(this)->bindToCache(cache);
}
const auto& testBasis() const
{
return static_cast<const Impl*>(this)->basis();
}
const auto& testTreePath() const
{
return static_cast<const Impl*>(this)->treePath();
}
};
template<class Impl>
struct UnaryForm<Impl, 1> : public Form
{
static constexpr std::size_t arity = 1;
template<class LocalView>
void bindToTrialLocalView(const LocalView& localView)
{
static_cast<Impl*>(this)->bindToLocalView(localView);
}
template<class Cache>
void bindToTrialCache(Cache& cache)
{
static_cast<Impl*>(this)->bindToCache(cache);
}
const auto& trialBasis() const
{
return static_cast<const Impl*>(this)->basis();
}
const auto& trialTreePath() const
{
return static_cast<const Impl*>(this)->treePath();
}
};
// Forward declarations
......@@ -254,7 +333,6 @@ namespace Impl {
class ConstCoefficient : public Form
{
public:
static constexpr std::size_t arity = 0;
using Element = E;
......@@ -262,28 +340,11 @@ namespace Impl {
ConstCoefficient(const C& c) :
c_(c)
{
Form::quadKey_ = QuadratureRuleKey(Element::dimension, 0);
}
template<class LocalViews>
void bind(const Element& element, const LocalViews& localViews)
{}
template<class LocalViews, class Cache>
void bindCache(const LocalViews& localViews, Cache& cache)
{}
template<std::size_t i>
auto treePathForArgument() const
auto quadratureRuleKey() const
{
return std::make_tuple();
}
template<std::size_t i>
auto basisForArgument() const
{
return std::make_tuple();
return QuadratureRuleKey(Element::dimension, 0);
}
auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
......@@ -306,7 +367,6 @@ namespace Impl {
using LocalFunction = std::decay_t<decltype(localFunction(std::declval<F&>()))>;
public:
static constexpr std::size_t arity = 0;
using Element = typename F::EntitySet::Element;
......@@ -315,37 +375,21 @@ namespace Impl {
Coefficient(const F& f) :
f_(f),
localF_(localFunction(f))
{
Form::quadKey_ = QuadratureRuleKey(Element::dimension, 0);
}
{}
Coefficient(const F& f, std::size_t order) :
f_(f),
localF_(localFunction(f))
{
Form::quadKey_ = QuadratureRuleKey(Element::dimension, order);
}
template<class LocalViews>
void bind(const Element& element, const LocalViews& localViews)
{
localF_.bind(element);
}
template<class LocalViews, class Cache>
void bindCache(const LocalViews& localViews, Cache& cache)
{}
template<std::size_t i>
auto treePathForArgument() const
auto quadratureRuleKey() const
{
return std::make_tuple();
return QuadratureRuleKey(Element::dimension, 0);
}
template<std::size_t i>
auto basisForArgument() const
void bindToElement(const Element& element)
{
return std::make_tuple();
localF_.bind(element);
}
auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
......@@ -365,7 +409,7 @@ namespace Impl {
// Linear map on FE-space that maps v onto itself
template<class B, class TP, std::size_t argIndex>
class FEFunctionForm : public Form
class FEFunctionForm : public UnaryForm<FEFunctionForm<B,TP,argIndex>, argIndex>
{
public:
using Basis = B;
......@@ -385,8 +429,6 @@ namespace Impl {
using Element = typename Basis::LocalView::Element;
using Coordinate = typename Element::Geometry::LocalCoordinate;
static constexpr std::size_t arity = 1;
FEFunctionForm(const Basis& basis, const TreePath& treePath) :
basis_(basis),
treePath_(treePath),
......@@ -394,35 +436,30 @@ namespace Impl {
leafNode_(nullptr)
{}
template<class LocalViews>
void bind(const Element& element, const LocalViews& localViews)
auto quadratureRuleKey() const
{
leafNode_ = & Dune::TypeTree::child(std::get<argIndex>(localViews).tree(), leafTreePath_);
Form::quadKey_ = QuadratureRuleKey(leafNode_->finiteElement());
return QuadratureRuleKey(leafNode_->finiteElement());
}
template<class LocalViews, class Cache>
void bindCache(const LocalViews& localViews, Cache& cache)
void bindToLocalView(const LocalView& localView)
{
valueCache_ = &(cache.getValues(std::get<argIndex>(localViews).tree(), leafTreePath_));
leafNode_ = & Dune::TypeTree::child(localView.tree(), leafTreePath_);
}
template<std::size_t i>
auto treePathForArgument() const
template<class Cache>
void bindToCache(Cache& cache)
{
if constexpr (i==argIndex)
return std::tie(treePath_);
else
return std::make_tuple();
valueCache_ = &(cache.getValues(leafTreePath_));
}
template<std::size_t i>
auto basisForArgument() const
const auto& treePath() const
{
if constexpr (i==argIndex)
return std::tie(basis_);
else
return std::make_tuple();
return treePath_;
}
const auto& basis() const
{
return basis_;
}
auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
......@@ -489,7 +526,7 @@ namespace Impl {
// Linear map on FE-space that maps v onto its Jacobian
template<class B, class TP, std::size_t argIndex>
class FEFunctionJacobianForm : public Form
class FEFunctionJacobianForm : public UnaryForm<FEFunctionJacobianForm<B,TP,argIndex>, argIndex>
{
public:
using Basis = B;
......@@ -509,8 +546,6 @@ namespace Impl {
using Element = typename Basis::LocalView::Element;
using Coordinate = typename Element::Geometry::LocalCoordinate;
static constexpr std::size_t arity = 1;
FEFunctionJacobianForm(const Basis& basis, const TreePath& treePath) :
basis_(basis),
treePath_(treePath),
......@@ -518,35 +553,30 @@ namespace Impl {
leafNode_(nullptr)
{}
template<class LocalViews>
void bind(const Element& element, const LocalViews& localViews)
auto quadratureRuleKey() const
{
leafNode_ = & Dune::TypeTree::child(std::get<argIndex>(localViews).tree(), leafTreePath_);
Form::quadKey_ = QuadratureRuleKey(leafNode_->finiteElement()).derivative();
return QuadratureRuleKey(leafNode_->finiteElement()).derivative();
}
template<class LocalViews, class Cache>
void bindCache(const LocalViews& localViews, Cache& cache)
void bindToLocalView(const LocalView& localView)
{
globalJacobianCache_ = &(cache.getGlobalJacobians(std::get<argIndex>(localViews).tree(), leafTreePath_));
leafNode_ = & Dune::TypeTree::child(localView.tree(), leafTreePath_);
}
template<std::size_t i>
auto treePathForArgument() const
template<class Cache>
void bindToCache(Cache& cache)
{
if constexpr (i==argIndex)
return std::tie(treePath_);
else
return std::make_tuple();
globalJacobianCache_ = &(cache.getGlobalJacobians(leafTreePath_));
}
template<std::size_t i>
auto basisForArgument() const
const auto& treePath() const
{
if constexpr (i==argIndex)
return std::tie(basis_);
else
return std::make_tuple();
return treePath_;
}
const auto& basis() const
{
return basis_;
}
auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
......@@ -559,8 +589,8 @@ namespace Impl {
else if constexpr(Node::isPower and Node::ChildType::isLeaf)
return [&](const auto& args) {
using F = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
constexpr auto components = Node::degree();
constexpr auto dimension = LeafNode::FiniteElement::Traits::LocalBasisType::Traits::dimDomain;
static constexpr auto components = Node::degree();
static constexpr auto dimension = LeafNode::FiniteElement::Traits::LocalBasisType::Traits::dimDomain;
auto componentIndex = std::get<argIndex>(args) / leafNode_->size();
auto shapeFunctionIndex = std::get<argIndex>(args) % leafNode_->size();
auto result = Dune::FieldMatrix<F,components,dimension>(0);
......@@ -589,7 +619,7 @@ namespace Impl {
// Linear map on FE-space that maps v onto its divergence
template<class B, class TP, std::size_t argIndex>
class FEFunctionDivergenceForm : public Form
class FEFunctionDivergenceForm : public UnaryForm<FEFunctionDivergenceForm<B,TP,argIndex>, argIndex>
{
public:
using Basis = B;
......@@ -609,8 +639,6 @@ namespace Impl {
using Element = typename Basis::LocalView::Element;
using Coordinate = typename Element::Geometry::LocalCoordinate;
static constexpr std::size_t arity = 1;
FEFunctionDivergenceForm(const Basis& basis, const TreePath& treePath) :
basis_(basis),
treePath_(treePath),
......@@ -618,35 +646,30 @@ namespace Impl {
leafNode_(nullptr)
{}
template<class LocalViews>
void bind(const Element& element, const LocalViews& localViews)
auto quadratureRuleKey() const
{
leafNode_ = & Dune::TypeTree::child(std::get<argIndex>(localViews).tree(), leafTreePath_);
Form::quadKey_ = QuadratureRuleKey(leafNode_->finiteElement()).derivative();
return QuadratureRuleKey(leafNode_->finiteElement()).derivative();
}
template<class LocalViews, class Cache>
void bindCache(const LocalViews& localViews, Cache& cache)
void bindToLocalView(const LocalView& localView)
{
globalJacobianCache_ = &(cache.getGlobalJacobians(std::get<argIndex>(localViews).tree(), leafTreePath_));
leafNode_ = & Dune::TypeTree::child(localView.tree(), leafTreePath_);
}
template<std::size_t i>
auto treePathForArgument() const
template<class Cache>
void bindToCache(Cache& cache)
{
if constexpr (i==argIndex)
return std::tie(treePath_);
else
return std::make_tuple();
globalJacobianCache_ = &(cache.getGlobalJacobians(leafTreePath_));
}
template<std::size_t i>
auto basisForArgument() const
const auto& treePath() const
{
if constexpr (i==argIndex)
return std::tie(basis_);
else
return std::make_tuple();
return treePath_;
}
const auto& basis() const
{
return basis_;
}
auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
......@@ -708,31 +731,75 @@ namespace Impl {
form1_(form1)
{}
template<class LocalViews>
void bind(const Element& element, const LocalViews& localViews)
auto quadratureRuleKey() const
{
return form0_.quadratureRuleKey().product(form1_.quadratureRuleKey());
}
void bindToElement(const Element& element)
{
form0_.bindToElement(element);
form1_.bindToElement(element);
}
template<class LocalView>
void bindToTestLocalView(const LocalView& localView)
{
form0_.bindToTestLocalView(localView);
form1_.bindToTestLocalView(localView);
}
template<class LocalView>
void bindToTrialLocalView(const LocalView& localView)
{
form0_.bindToTrialLocalView(localView);
form1_.bindToTrialLocalView(localView);
}
template<class Cache>
void bindToTestCache(Cache& cache)
{
form0_.bindToTestCache(cache);
form1_.bindToTestCache(cache);
}
template<class Cache>
void bindToTrialCache(Cache& cache)
{
form0_.bindToTrialCache(cache);
form1_.bindToTrialCache(cache);
}
decltype(auto) testBasis() const
{
form0_.bind(element, localViews);
form1_.bind(element, localViews);
Form::quadKey_ = form0_.quadratureRuleKey().product(form1_.quadratureRuleKey());
if constexpr (std::is_same_v<decltype(form1_.testBasis()), void>)
return form0_.testBasis();
else
return form1_.testBasis();
}
template<class LocalViews, class Cache>
void bindCache(const LocalViews& localViews, Cache& cache)
decltype(auto) trialBasis() const
{
form0_.bindCache(localViews, cache);
form1_.bindCache(localViews, cache);
if constexpr (std::is_same_v<decltype(form1_.trialBasis()), void>)
return form0_.trialBasis();
else
return form1_.trialBasis();
}
template<std::size_t k>
auto treePathForArgument() const
decltype(auto) testTreePath() const
{
return std::tuple_cat(form0_.template treePathForArgument<k>(), form1_.template treePathForArgument<k>());
if constexpr (std::is_same_v<decltype(form1_.testTreePath()), void>)
return form0_.testTreePath();
else
return form1_.testTreePath();
}
template<std::size_t k>
auto basisForArgument() const
decltype(auto) trialTreePath() const
{
return std::tuple_cat(form0_.template basisForArgument<k>(), form1_.template basisForArgument<k>());
if constexpr (std::is_same_v<decltype(form1_.trialTreePath()), void>)
return form0_.trialTreePath();
else
return form1_.trialTreePath();
}
auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
......@@ -758,7 +825,6 @@ namespace Impl {
using Transformation = T;
public:
static constexpr std::size_t arity = Form0::arity;
using Element = typename Form0::Element;
......@@ -769,29 +835,58 @@ namespace Impl {
form0_(form0)
{}
template<class LocalViews>
void bind(const Element& element, const LocalViews& localViews)
auto quadratureRuleKey() const
{
return form0_.quadratureRuleKey();
}
void bindToElement(const Element& element)
{
form0_.bindToElement(element);
}
template<class LocalView>
void bindToTestLocalView(const LocalView& localView)
{
form0_.bindToTestLocalView(localView);
}
template<class LocalView>
void bindToTrialLocalView(const LocalView& localView)
{
form0_.bindToTrialLocalView(localView);
}
template<class Cache>
void bindToTestCache(Cache& cache)
{
form0_.bindToTestCache(cache);
}
template<class Cache>
void bindToTrialCache(Cache& cache)
{
form0_.bindToTrialCache(cache);
}
decltype(auto) testBasis() const
{
form0_.bind(element, localViews);
Form::quadKey_ = form0_.quadratureRuleKey();
return form0_.testBasis();
}
template<class LocalViews, class Cache>
void bindCache(const LocalViews& localViews, Cache& cache)
decltype(auto) trialBasis() const
{
form0_.bindCache(localViews, cache);
return form0_.trialBasis();
}
template<std::size_t k>
auto treePathForArgument() const
decltype(auto) testTreePath() const
{
return form0_.template treePathForArgument<k>();
return form0_.testTreePath();
}
template<std::size_t k>
auto basisForArgument() const
decltype(auto) trialTreePath() const
{
return form0_.template basisForArgument<k>();
return form0_.trialTreePath();
}
auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
......@@ -820,6 +915,8 @@ namespace Impl {
{
using Base = TransformedForm<Dune::Fufem::Forms::Impl::TransposeOp, F>;
public:
using Base::arity;
TransposedForm(const F& form) :
Base(Dune::Fufem::Forms::Impl::TransposeOp(), form)
{}
......@@ -833,7 +930,6 @@ namespace Impl {
class SumForm : public Form
{
public:
static constexpr std::size_t arity = Form0::arity;
using Element = typename Form0::Element;
......@@ -843,26 +939,64 @@ namespace Impl {
forms_(form0, forms...)
{}
template<class LocalViews>
void bind(const Element& element, const LocalViews& localViews)
auto quadratureRuleKey() const
{
auto quadKey = QuadratureRuleKey();
Impl::forEachTupleEntry(forms_, [&](const auto& form) {
quadKey = form.quadratureRuleKey().sum(quadKey);
});
return quadKey;
}
void bindToElement(const Element& element)
{
Impl::forEachTupleEntry(forms_, [&](auto& form) {
form.bindToElement(element);
});
}
template<class LocalView>
void bindToTestLocalView(const LocalView& localView)
{
Impl::forEachTupleEntry(forms_, [&](auto& form) {
form.bindToTestLocalView(localView);
});
}
template<class LocalView>
void bindToTrialLocalView(const LocalView& localView)
{
Impl::forEachTupleEntry(forms_, [&](auto& form) {
form.bind(element, localViews);
form.bindToTrialLocalView(localView);
});
Form::quadKey_ = std::get<0>(forms_).quadratureRuleKey();
}
template<class Cache>
void bindToTestCache(Cache& cache)
{
Impl::forEachTupleEntry(forms_, [&](auto& form) {
Form::quadKey_ = form.quadratureRuleKey().sum(Form::quadKey_);
form.bindToTestCache(cache);
});
}
template<class LocalViews, class Cache>
void bindCache(const LocalViews& localViews, Cache& cache)
template<class Cache>
void bindToTrialCache(Cache& cache)
{
Impl::forEachTupleEntry(forms_, [&](auto& form) {
form.bindCache(localViews, cache);
form.bindToTrialCache(cache);
});
}
decltype(auto) testBasis() const
{
return std::get<0>(forms_).testBasis();
}
decltype(auto) trialBasis() const
{
return std::get<0>(forms_).trialBasis();
}
auto& forms()
{
return forms_;
......@@ -1008,7 +1142,7 @@ namespace Impl {
// Overload for two plain forms
template<class Op, class L, class R,
std::enable_if_t<isForm_v<L> and isForm_v<R>, int> = 0,
class = Impl::ProductConcept<Op, Impl::FormRange<L>, Impl::FormRange<R>>>
std::enable_if_t<Dune::IsCallable<Op(Impl::FormRange<L>, Impl::FormRange<R>)>::value, int> = 0>
auto productForm(const Op& op, const L& l, const R& r)
{
return ProductForm(op, l, r);
......@@ -1018,7 +1152,7 @@ namespace Impl {
// This will wrap the raw value into a form.
template<class Op, class L, class R,
std::enable_if_t<(not isForm_v<L>) and isForm_v<R>, int> = 0,
class = Impl::ProductConcept<Op, L, Impl::FormRange<R>>>
std::enable_if_t<Dune::IsCallable<Op(L, Impl::FormRange<R>)>::value, int> = 0>
auto productForm(const Op& op, const L& l, const R& r)
{
return ProductForm(op, ConstCoefficient<typename R::Element,L>(l), r);
......@@ -1028,7 +1162,7 @@ namespace Impl {
// This will wrap the raw value into a form.
template<class Op, class L, class R,
std::enable_if_t<isForm_v<L> and (not isForm_v<R>), int> = 0,
class = Impl::ProductConcept<Op, Impl::FormRange<L>, R>>
std::enable_if_t<Dune::IsCallable<Op(Impl::FormRange<L>, R)>::value, int> = 0>
auto productForm(const Op& op, const L& l, const R& r)
{
return ProductForm(op, l, ConstCoefficient<typename L::Element,R>(r));
......@@ -1176,50 +1310,10 @@ namespace Impl {
template<std::size_t argIndex, class Form>
auto&& treePathForArgument(const Form& form)
{
// For each form treePathForArgument<k>() is either empty
// or a treePath. This is statically known. The simplest way
// to implement a static 'optional' is to collect
// the results for sub-forms is a tuple. The final form
// is required to have a unique treePath for each argument
// in the range from 0 to arity-1. Hence the tuple must be
// nonempty and we can simply pick its first entry.
return std::get<0>(form.template treePathForArgument<argIndex>());
}
template<std::size_t argIndex, class Form>
auto&& basisForArgument(const Form& form)
{
// For each form basisForArgument<k>() is either empty
// or a basis. This is statically known. The simplest way
// to implement a static 'optional' is to collect
// the results for sub-forms is a tuple. The final form
// is required to have a unique basis for each argument
// in the range from 0 to arity-1. Hence the tuple must be
// nonempty and we can simply pick its first entry.
return std::get<0>(form.template basisForArgument<argIndex>());
}
template<std::size_t argIndex, class... Forms>
auto&& basisForArgument(const SumForm<Forms...>& sumForm)
{
// For each form basisForArgument<k>() is either empty
// or a basis. This is statically known. The simplest way
// to implement a static 'optional' is to collect
// the results for sub-forms is a tuple. The final form
// is required to have a unique basis for each argument
// in the range from 0 to arity-1. Hence the tuple must be
// nonempty and we can simply pick its first entry.
return std::get<0>(std::get<0>(sumForm.forms()).template basisForArgument<argIndex>());
}
template<class Form, bool isSemiAffine=true>
class IntegratedLinearForm
{
using TestRootBasis = std::decay_t<decltype(basisForArgument<0>(std::declval<Form>()))>;
using TestRootBasis = std::decay_t<decltype(std::declval<Form>().testBasis())>;
using TestTree = typename TestRootBasis::LocalView::Tree;
......@@ -1228,7 +1322,7 @@ namespace Impl {
IntegratedLinearForm(const Form& sumForm) :
sumForm_(sumForm),
testRootBasis_(basisForArgument<0>(sumForm_))
testRootBasis_(sumForm_.testBasis())
{}
template<class LocalVector, class TestLocalView>
......@@ -1238,13 +1332,14 @@ namespace Impl {
const auto& geometry = element.geometry();
sumForm_.bind(element, std::tie(testLocalView));
sumForm_.bindToElement(testLocalView.element());
sumForm_.bindToTestLocalView(testLocalView);
auto& cacheForRule = cache_[sumForm_.quadratureRuleKey()];
cacheForRule.invalidateElementWiseCaches(testLocalView.tree());
if (not isSemiAffine)
cacheForRule.invalidateAllCaches(testLocalView.tree());
sumForm_.bindCache(std::tie(testLocalView), cacheForRule);
cacheForRule.bind(testLocalView.tree());
cacheForRule.invalidate(isSemiAffine);
sumForm_.bindToTestCache(cacheForRule);
const auto& quad = cacheForRule.rule();
......@@ -1257,7 +1352,7 @@ namespace Impl {
Impl::forEachTupleEntry(sumForm_.forms(), [&](auto& form) {
auto evaluatedForm = form(quadPoint.position(), k);
const auto& testPrefix = treePathForArgument<0>(form);
const auto& testPrefix = form.testTreePath();
const auto& testNode = Dune::TypeTree::child(testLocalView.tree(), testPrefix);
for (std::size_t i=0; i<testNode.size(); ++i)
{
......@@ -1284,8 +1379,8 @@ namespace Impl {
template<class Form, bool isSemiAffine=true>
class IntegratedBilinearForm
{
using TestRootBasis = std::decay_t<decltype(basisForArgument<0>(std::declval<Form>()))>;
using AnsatzRootBasis = std::decay_t<decltype(basisForArgument<1>(std::declval<Form>()))>;
using TestRootBasis = std::decay_t<decltype(std::declval<Form>().testBasis())>;
using AnsatzRootBasis = std::decay_t<decltype(std::declval<Form>().trialBasis())>;
using TestTree = typename TestRootBasis::LocalView::Tree;
......@@ -1294,8 +1389,8 @@ namespace Impl {
IntegratedBilinearForm(const Form& sumForm) :
sumForm_(sumForm),
testRootBasis_(basisForArgument<0>(sumForm_)),
ansatzRootBasis_(basisForArgument<1>(sumForm_))
testRootBasis_(sumForm_.testBasis()),
ansatzRootBasis_(sumForm_.trialBasis())
{}
template<class LocalMatrix, class TestLocalView, class AnsatzLocalView>
......@@ -1306,13 +1401,16 @@ namespace Impl {
const auto& geometry = element.geometry();
sumForm_.bind(element, std::tie(testLocalView, ansatzLocalView));
sumForm_.bindToElement(testLocalView.element());
sumForm_.bindToTestLocalView(testLocalView);
sumForm_.bindToTrialLocalView(ansatzLocalView);
auto& cacheForRule = cache_[sumForm_.quadratureRuleKey()];
cacheForRule.invalidateElementWiseCaches(testLocalView.tree());
if (not isSemiAffine)
cacheForRule.invalidateAllCaches(testLocalView.tree());
sumForm_.bindCache(std::tie(testLocalView,ansatzLocalView), cacheForRule);
cacheForRule.bind(testLocalView.tree());
cacheForRule.invalidate(isSemiAffine);
sumForm_.bindToTestCache(cacheForRule);
sumForm_.bindToTrialCache(cacheForRule);
const auto& quad = cacheForRule.rule();
......@@ -1325,9 +1423,9 @@ namespace Impl {
Impl::forEachTupleEntry(sumForm_.forms(), [&](auto& form) {
auto evaluatedForm = form(quadPoint.position(), k);
const auto& testPrefix = treePathForArgument<0>(form);
const auto& testPrefix = form.testTreePath();
const auto& testNode = Dune::TypeTree::child(testLocalView.tree(), testPrefix);
const auto& ansatzPrefix = treePathForArgument<1>(form);
const auto& ansatzPrefix = form.trialTreePath();
const auto& ansatzTree = Dune::TypeTree::child(ansatzLocalView.tree(), ansatzPrefix);
for (std::size_t i=0; i<testNode.size(); ++i)
{
......@@ -1359,7 +1457,7 @@ namespace Impl {
template<class Form, class Patch, bool isSemiAffine=true>
class IntegratedBoundaryLinearForm
{
using TestRootBasis = std::decay_t<decltype(basisForArgument<0>(std::declval<Form>()))>;
using TestRootBasis = std::decay_t<decltype(std::declval<Form>().testBasis())>;
using TestTree = typename TestRootBasis::LocalView::Tree;
......@@ -1370,7 +1468,7 @@ namespace Impl {
IntegratedBoundaryLinearForm(const Form& sumForm, const Patch& patch) :
sumForm_(sumForm),
testRootBasis_(basisForArgument<0>(sumForm_)),
testRootBasis_(sumForm_.testBasis()),
patch_(patch)
{}
......@@ -1387,19 +1485,19 @@ namespace Impl {
template<class LocalVector, class TestLocalView>
void operator()(const Intersection& intersection, LocalVector& localVector, const TestLocalView& testSubspaceLocalView)
{
const auto& element = intersection.inside();
auto facet = intersection.indexInInside();
const auto& testLocalView = testSubspaceLocalView.rootLocalView();
const auto& geometry = intersection.geometry();
sumForm_.bind(element, std::tie(testLocalView));
sumForm_.bindToElement(testLocalView.element());
sumForm_.bindToTestLocalView(testLocalView);
auto& cacheForRule = cache_[std::pair(sumForm_.quadratureRuleKey(), facet)];
cacheForRule.invalidateElementWiseCaches(testLocalView.tree());
if (not isSemiAffine)
cacheForRule.invalidateAllCaches(testLocalView.tree());
sumForm_.bindCache(std::tie(testLocalView), cacheForRule);
cacheForRule.bind(testLocalView.tree());
cacheForRule.invalidate(isSemiAffine);
sumForm_.bindToTestCache(cacheForRule);
const auto& quad = cacheForRule.rule();
......@@ -1413,7 +1511,7 @@ namespace Impl {
Impl::forEachTupleEntry(sumForm_.forms(), [&](auto& form) {
auto evaluatedForm = form(quadPoint.position(), k);
const auto& testPrefix = treePathForArgument<0>(form);
const auto& testPrefix = form.testTreePath();
const auto& testNode = Dune::TypeTree::child(testLocalView.tree(), testPrefix);
for (std::size_t i=0; i<testNode.size(); ++i)
{
......@@ -1578,4 +1676,3 @@ namespace Impl {
#endif // DUNE_FUFEM_ASSEMBLERS_FORMS_HH
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment