From 3a6b1fea825a1c41dfd3d87fa4efef24a3d43184 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@mi.fu-berlin.de>
Date: Sun, 19 Mar 2023 18:31:58 +0100
Subject: [PATCH] Cleanup and improve
* Split header
* Detangle dependencies
* Use better names in user interface
* `FormOperator` -> `RangeOperator`
* `transformedForm(op, form)` -> `transformRange(form, op)`
* Remove `makeSumForm()`, using plain `SumForm(...)` does the trick due to CTAD
* Add operator `D` for computing partial derivatives as `D(u, i)`
* Further cleaup of form interface. Less magic, more explicit.
* Add support for passing more arguemnts to `RangeOperator`
---
doc/CMakeLists.txt | 1 +
dune/fufem-forms/CMakeLists.txt | 11 +-
dune/fufem-forms/baseclass.hh | 105 +++
dune/fufem-forms/forms.hh | 1347 +--------------------------
dune/fufem-forms/fufem-forms.hh | 6 -
dune/fufem-forms/nullaryforms.hh | 74 ++
dune/fufem-forms/operators.hh | 213 +++++
dune/fufem-forms/productform.hh | 248 +++++
dune/fufem-forms/sumform.hh | 156 ++++
dune/fufem-forms/transformedform.hh | 158 ++++
dune/fufem-forms/unaryforms.hh | 400 ++++++++
src/fracture.cc | 2 +-
12 files changed, 1398 insertions(+), 1323 deletions(-)
create mode 100644 dune/fufem-forms/baseclass.hh
delete mode 100644 dune/fufem-forms/fufem-forms.hh
create mode 100644 dune/fufem-forms/nullaryforms.hh
create mode 100644 dune/fufem-forms/operators.hh
create mode 100644 dune/fufem-forms/productform.hh
create mode 100644 dune/fufem-forms/sumform.hh
create mode 100644 dune/fufem-forms/transformedform.hh
create mode 100644 dune/fufem-forms/unaryforms.hh
diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt
index be52cfc..4828ea3 100644
--- a/doc/CMakeLists.txt
+++ b/doc/CMakeLists.txt
@@ -1 +1,2 @@
add_subdirectory("doxygen")
+add_subdirectory("manual")
diff --git a/dune/fufem-forms/CMakeLists.txt b/dune/fufem-forms/CMakeLists.txt
index bc3da27..343ec23 100644
--- a/dune/fufem-forms/CMakeLists.txt
+++ b/dune/fufem-forms/CMakeLists.txt
@@ -1,2 +1,11 @@
#install headers
-install(FILES fufem-forms.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/fufem-forms)
+install(FILES
+ baseclass.hh
+ forms.hh
+ nullaryforms.hh
+ operators.hh
+ productform.hh
+ sumform.hh
+ transformedform.hh
+ unaryforms.hh
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/fufem-forms)
diff --git a/dune/fufem-forms/baseclass.hh b/dune/fufem-forms/baseclass.hh
new file mode 100644
index 0000000..c51df77
--- /dev/null
+++ b/dune/fufem-forms/baseclass.hh
@@ -0,0 +1,105 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FUFEM_FORMS_BASECLASS_HH
+#define DUNE_FUFEM_FORMS_BASECLASS_HH
+
+#include <cstddef>
+#include <type_traits>
+#include <tuple>
+#include <utility>
+
+#include <dune/common/typetraits.hh>
+#include <dune/common/indices.hh>
+
+#include <dune/typetree/childextraction.hh>
+#include <dune/typetree/treepath.hh>
+
+
+
+namespace Dune::Fufem::Forms::Impl {
+
+ template<class Tuple, class F>
+ void forEachTupleEntry(Tuple&& tuple, F&& f){
+ std::apply([&](auto&&... entry){
+ (f(std::forward<decltype(entry)>(entry)), ...);
+ }, tuple);
+ }
+
+ // Helper function to derive the leafTreePath associated
+ // to a treePath. For treePaths refering to a leaf this is
+ // the identity. For treePaths referring to the a power
+ // node whose children are leaf, it's the first child.
+ // This is used to implement vector valued ansatz functions
+ // for power nodes.
+ template<class Tree, class TreePath>
+ static auto leafTreePath(const TreePath& tp)
+ {
+ using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
+ if constexpr(Node::isLeaf)
+ return tp;
+ else if constexpr(Node::isPower and Node::ChildType::isLeaf)
+ return Dune::TypeTree::push_back(tp, Dune::Indices::_0);
+ }
+
+ template<class Tree, class TreePath>
+ using LeafTreePath = decltype(leafTreePath<Tree>(std::declval<TreePath>()));
+
+ template<class F>
+ using LazyFormRange = decltype(std::declval<F>()(std::declval<typename F::Coordinate>(), 0));
+
+ // Helper template to determine the range after bind, evaluation, and index access.
+ template<class F>
+ using FormRange = decltype(std::declval<LazyFormRange<F>>()(std::tuple<int,int>{}));
+
+} // namespace Impl
+
+
+
+namespace Dune::Fufem::Forms {
+
+
+
+ // Base class of form implementations. The term _form_ follows the
+ // terminalogy in UFL, altough it does stricly speaking not refer
+ // to a multilinear form here. Instead all forms are multilinear
+ // maps from finite element spaces into the space of functions
+ // from the domain into some finite dimensional vector spaces.
+ // Poinwise products (scalar-vector, matrix-vector, dot,...) of
+ // multilinear maps induces multilinear maps of higher order.
+ // The product of a k-linear map and an l-linear is a (k+l) linear
+ // map into a suitable function space.
+ //
+ // Example:
+ // Denote by {A->B} the set of all functions f:B->A. Now let
+ // D a domain in R^d, V \subset {D->R} and W \subset {D->R}
+ // two FE-spaces, and f \in {D->R^(d,d)} a fixed matrix valued
+ // function.
+ //
+ // Then f is a 0-linear map into {D->R^(d,d)}.
+ // The gradient operator \nabla:V -> {D->R^d} is a 1-linear map into {D->R^d}.
+ // The gradient operator \nabla:W -> {D->R^d} is a 1-linear map into {D->R^d}.
+ // The pointwise matrix-vector product f\nabla:V -> {D->R^d} is a 1-linear map into {D->R^d}.
+ // The pointwise dot product dot(f\nabla(.), \nabla(.)):V\times W -> {D->R}
+ // is a 2-linear map into {D->R}.
+ //
+ // Integrating a k-linear map F:V_1 \times ... \times V_k -> {D->R} over D
+ // results in a k-linear form \int_D F(...)dx : V_1 \times ... \times V_k -> D
+ class Form {};
+
+
+
+ // Check if a class is a form
+
+ template<class F>
+ using isForm = std::is_base_of<Form, std::decay_t<F>>;
+
+ template<class F>
+ inline constexpr bool isForm_v = isForm<F>::value;
+
+
+
+} // namespace Dune::Fufem::Forms
+
+
+
+#endif // DUNE_FUFEM_FORMS_BASECLASS_HH
diff --git a/dune/fufem-forms/forms.hh b/dune/fufem-forms/forms.hh
index 5cb108b..d00c496 100644
--- a/dune/fufem-forms/forms.hh
+++ b/dune/fufem-forms/forms.hh
@@ -3,1317 +3,36 @@
#ifndef DUNE_FUFEM_ASSEMBLERS_FORMS_HH
#define DUNE_FUFEM_ASSEMBLERS_FORMS_HH
+#include <cstddef>
#include <type_traits>
#include <tuple>
#include <utility>
-#include <dune/common/typetraits.hh>
-
-#include <dune/istl/matrix.hh>
-#include <dune/istl/matrixindexset.hh>
-
-#include <dune/functions/functionspacebases/subspacebasis.hh>
+#include <dune/typetree/childextraction.hh>
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem-forms/shapefunctioncache.hh>
+#include <dune/fufem-forms/baseclass.hh>
+#include <dune/fufem-forms/nullaryforms.hh>
+#include <dune/fufem-forms/transformedform.hh>
+#include <dune/fufem-forms/unaryforms.hh>
+#include <dune/fufem-forms/sumform.hh>
+#include <dune/fufem-forms/productform.hh>
+#include <dune/fufem-forms/operators.hh>
-namespace Dune::Fufem::Forms {
-
-
-namespace Impl {
-
- 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 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 {
- template<class K, int n, int m>
- auto operator()(const Dune::FieldMatrix<K,n,m>& x) const
- {
- Dune::FieldMatrix<K,m,n> y;
- for(std::size_t i=0; i<x.N(); ++i)
- for(std::size_t j=0; j<x.M(); ++j)
- y[j][i] = x[i][j];
- return y;
- }
- };
-
- struct DotOp {
-
- template<class K1, class K2, int n>
- auto operator()(const Dune::FieldVector<K1,n>& x, const Dune::FieldVector<K2,n>& y) const
- {
- using K = typename PromotionTraits<K1,K2>::PromotedType;
- auto result = K(0);
- for(std::size_t i=0; i<n; ++i)
- result += x[i]*y[i];
- return result;
- }
-
- template<class K1, class K2, int n, int m>
- auto operator()(const Dune::FieldMatrix<K1,n,m>& x, const Dune::FieldMatrix<K2,n,m>& y) const
- {
- using K = typename PromotionTraits<K1,K2>::PromotedType;
- auto result = K(0);
- for(std::size_t i=0; i<n; ++i)
- for(std::size_t j=0; j<m; ++j)
- result += x[i][j]*y[i][j];
- return result;
- }
-
- template<class K1, class K2, int n>
- auto operator()(const Dune::FieldMatrix<K1,n,1>& x, const Dune::FieldVector<K2,n>& y) const
- {
- using K = typename PromotionTraits<K1,K2>::PromotedType;
- auto result = K(0);
- for(std::size_t i=0; i<n; ++i)
- result += x[i][0]*y[i];
- return result;
- }
-
- template<class K1, class K2, int n>
- auto operator()(const Dune::FieldVector<K1,n>& x, const Dune::FieldMatrix<K2,n,1>& y) const
- {
- using K = typename PromotionTraits<K1,K2>::PromotedType;
- auto result = K(0);
- for(std::size_t i=0; i<n; ++i)
- result += x[i]*y[i][0];
- return result;
- }
-
- template<class K1, class K2, int n>
- auto operator()(const Dune::ScaledIdentityMatrix<K1,n>& x, const Dune::FieldMatrix<K2,n,n>& y) const
- {
- using K = typename PromotionTraits<K1,K2>::PromotedType;
- auto result = K(0);
- for(std::size_t i=0; i<x.N(); ++i)
- result += y[i][i];
- return x.scalar()*result;
- }
-
- template<class K1, class K2, int n>
- auto operator()(const Dune::FieldMatrix<K1,n,n>& x, const Dune::ScaledIdentityMatrix<K2,n>& y) const
- {
- using K = typename PromotionTraits<K1,K2>::PromotedType;
- auto result = K(0);
- for(std::size_t i=0; i<x.N(); ++i)
- result += x[i][i];
- return result*y.scalar();
- }
-
- template<class K1, class K2, int n>
- auto operator()(const Dune::ScaledIdentityMatrix<K1,n>& x, const Dune::ScaledIdentityMatrix<K2,n>& y) const
- {
- return n*x.scalar()*y.scalar();
- }
-
- };
-
- struct MultOp {
-
- template<class X, class Y,
- class = std::void_t<decltype(std::declval<X>()*std::declval<Y>())>>
- auto operator()(const X& x, const Y& y) const
- {
- return x*y;
- }
-
- template<class K1, class K2, int n, int m>
- auto operator()(const Dune::FieldMatrix<K1,n,m>& x, const Dune::FieldVector<K2,m>& y) const
- {
- using K = typename PromotionTraits<K1,K2>::PromotedType;
- Dune::FieldVector<K,n> result;
- x.mv(y, result);
- return result;
- }
-
- template<class K, class Y>
- auto operator()(const Dune::FieldMatrix<K,1,1>& x, const Y& y) const
- {
- return x[0][0]*y;
- }
-
- template<class K, class Y>
- auto operator()(const Dune::FieldVector<K,1>& x, const Y& y) const
- {
- return x[0]*y;
- }
-
- };
-
-
- // Helper function to derive the leafTreePath associated
- // to a treePath. For treePaths refering to a leaf this is
- // the identity. For treePaths referring to the a power
- // node whose children are leaf, it's the first child.
- // This is used to implement vector valued ansatz functions
- // for power nodes.
- template<class Tree, class TreePath>
- static auto leafTreePath(const TreePath& tp)
- {
- using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
- if constexpr(Node::isLeaf)
- return tp;
- else if constexpr(Node::isPower and Node::ChildType::isLeaf)
- return Dune::TypeTree::push_back(tp, Dune::Indices::_0);
- }
-
- template<class Tree, class TreePath>
- using LeafTreePath = decltype(leafTreePath<Tree>(std::declval<TreePath>()));
-
- template<class F>
- using LazyFormRange = decltype(std::declval<F>()(std::declval<typename F::Coordinate>(), 0));
-
- // Helper template to determine the range after bind, evaluation, and index access.
- template<class F>
- using FormRange = decltype(std::declval<LazyFormRange<F>>()(std::tuple<int,int>{}));
-
-} // namespace Impl
-
-
-
- // Base class of form implementations. The term _form_ follows the
- // terminalogy in UFL, altough it does stricly speaking not refer
- // to a multilinear form here. Instead all forms are multilinear
- // maps from finite element spaces into the space of functions
- // from the domain into some finite dimensional vector spaces.
- // Poinwise products (scalar-vector, matrix-vector, dot,...) of
- // multilinear maps induces multilinear maps of higher order.
- // The product of a k-linear map and an l-linear is a (k+l) linear
- // map into a suitable function space.
- //
- // Example:
- // Denote by {A->B} the set of all functions f:B->A. Now let
- // D a domain in R^d, V \subset {D->R} and W \subset {D->R}
- // two FE-spaces, and f \in {D->R^(d,d)} a fixed matrix valued
- // function.
- //
- // Then f is a 0-linear map into {D->R^(d,d)}.
- // The gradient operator \nabla:V -> {D->R^d} is a 1-linear map into {D->R^d}.
- // The gradient operator \nabla:W -> {D->R^d} is a 1-linear map into {D->R^d}.
- // The pointwise matrix-vector product f\nabla:V -> {D->R^d} is a 1-linear map into {D->R^d}.
- // The pointwise dot product dot(f\nabla(.), \nabla(.)):V\times W -> {D->R}
- // is a 2-linear map into {D->R}.
- //
- // Integrating a k-linear map F:V_1 \times ... \times V_k -> {D->R} over D
- // results in a k-linear form \int_D F(...)dx : V_1 \times ... \times V_k -> D
- class Form {
- public:
-
- 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 {}
- };
-
-
-
- // Check if a class is a form
-
- template<class F>
- using isForm = std::is_base_of<Form, std::decay_t<F>>;
-
- template<class F>
- 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
-
- template<class F,
- std::enable_if_t<isForm_v<F>, int> = 0>
- auto transpose(const F& f);
-
- template<class B, class TP, std::size_t argIndex>
- class FEFunctionJacobianForm;
-
- template<class B, class TP, std::size_t argIndex>
- class FEFunctionDivergenceForm;
-
- template<class B, class TP, std::size_t argIndex>
- class FEFunctionForm;
-
-
-
-
-
- // Zero-linear map induced by fixed coefficient function
- template<class E, class C>
- class ConstCoefficient : public Form
- {
- public:
- static constexpr std::size_t arity = 0;
-
- using Element = E;
- using Coordinate = typename Element::Geometry::LocalCoordinate;
-
- ConstCoefficient(const C& c) :
- c_(c)
- {}
-
- auto quadratureRuleKey() const
- {
- return QuadratureRuleKey(Element::dimension, 0);
- }
-
- auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
- {
- return [=](const auto& args) {
- return c_;
- };
- }
-
- private:
- C c_;
- };
-
-
-
- // Zero-linear map induced by fixed coefficient function
- template<class F>
- class Coefficient : public Form
- {
- using LocalFunction = std::decay_t<decltype(localFunction(std::declval<F&>()))>;
-
- public:
- static constexpr std::size_t arity = 0;
-
- using Element = typename F::EntitySet::Element;
- using Coordinate = typename Element::Geometry::LocalCoordinate;
-
- Coefficient(const F& f) :
- f_(f),
- localF_(localFunction(f))
- {}
-
- Coefficient(const F& f, std::size_t order) :
- f_(f),
- localF_(localFunction(f))
- {}
-
- auto quadratureRuleKey() const
- {
- return QuadratureRuleKey(Element::dimension, 0);
- }
-
- void bindToElement(const Element& element)
- {
- localF_.bind(element);
- }
-
- auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
- {
- auto value = localF_(x);
- return [value](const auto& args) {
- return value;
- };
- }
-
- private:
- F f_;
- LocalFunction localF_;
- };
-
-
-
- // Linear map on FE-space that maps v onto itself
- template<class B, class TP, std::size_t argIndex>
- class FEFunctionForm : public UnaryForm<FEFunctionForm<B,TP,argIndex>, argIndex>
- {
- public:
- using Basis = B;
- using TreePath = TP;
- using LocalView = typename Basis::LocalView;
- using Tree = typename LocalView::Tree;
- using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
-
- private:
- using LeafNode = typename Dune::TypeTree::ChildForTreePath<Tree, Impl::LeafTreePath<Tree,TreePath>>;
-
- using Buffer = std::vector<typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::RangeType>;
- using ValueCache = typename ShapeFunctionCache<Tree>::template NodalValueCache<LeafNode>;
-
- public:
-
- using Element = typename Basis::LocalView::Element;
- using Coordinate = typename Element::Geometry::LocalCoordinate;
-
- FEFunctionForm(const Basis& basis, const TreePath& treePath) :
- basis_(basis),
- treePath_(treePath),
- leafTreePath_(Impl::leafTreePath<Tree>(treePath)),
- leafNode_(nullptr)
- {}
-
- auto quadratureRuleKey() const
- {
- return QuadratureRuleKey(leafNode_->finiteElement());
- }
-
- void bindToLocalView(const LocalView& localView)
- {
- leafNode_ = & Dune::TypeTree::child(localView.tree(), leafTreePath_);
- }
-
- template<class Cache>
- void bindToCache(Cache& cache)
- {
- valueCache_ = &(cache.getValues(leafTreePath_));
- }
-
- const auto& treePath() const
- {
- return treePath_;
- }
-
- const auto& basis() const
- {
- return basis_;
- }
-
- auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
- {
- const auto& values = (*valueCache_)[quadPointIndex];
- if constexpr(Node::isLeaf)
- return [&](const auto& args) {
- return values[std::get<argIndex>(args)];
- };
- else if constexpr(Node::isPower and Node::ChildType::isLeaf)
- return [&](const auto& args) {
- using Field = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
- constexpr auto components = Node::degree();
- auto componentIndex = std::get<argIndex>(args) / leafNode_->size();
- auto shapeFunctionIndex = std::get<argIndex>(args) % leafNode_->size();
- auto result = Dune::FieldVector<Field,components>(0);
- result[componentIndex] = values[shapeFunctionIndex];
- return result;
- };
- }
-
- friend auto jacobian(const FEFunctionForm& f)
- {
- return FEFunctionJacobianForm<Basis, TreePath, argIndex>(f.basis_, f.treePath_);
- }
-
- friend auto gradient(const FEFunctionForm& f)
- {
- return Fufem::Forms::transpose(jacobian(f));
- }
-
- friend auto grad(const FEFunctionForm& f)
- {
- return Fufem::Forms::transpose(jacobian(f));
- }
-
- friend auto divergence(const FEFunctionForm& f)
- {
- return FEFunctionDivergenceForm<Basis, TreePath, argIndex>(f.basis_, f.treePath_);
- }
-
- friend auto div(const FEFunctionForm& f)
- {
- return divergence(f);
- }
-
- template<std::size_t i>
- auto childForm(Dune::index_constant<i> childIndex) const
- {
- auto childTP = Dune::TypeTree::push_back(treePath_, childIndex);
- return FEFunctionForm<Basis, decltype(childTP), argIndex>(basis_, childTP);
- }
-
- private:
- const Basis& basis_;
- const TreePath treePath_;
- Impl::LeafTreePath<Tree,TreePath> leafTreePath_;
- const LeafNode* leafNode_;
- mutable Buffer evaluationBuffer_;
- const ValueCache* valueCache_;
- };
-
-
-
- // Linear map on FE-space that maps v onto its Jacobian
- template<class B, class TP, std::size_t argIndex>
- class FEFunctionJacobianForm : public UnaryForm<FEFunctionJacobianForm<B,TP,argIndex>, argIndex>
- {
- public:
- using Basis = B;
- using TreePath = TP;
- using LocalView = typename Basis::LocalView;
- using Tree = typename LocalView::Tree;
- using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
-
- private:
- using LeafNode = typename Dune::TypeTree::ChildForTreePath<Tree, Impl::LeafTreePath<Tree,TreePath>>;
-
- using Buffer = std::vector<typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::JacobianType>;
- using GlobalJacobianCache = typename ShapeFunctionCache<Tree>::template NodalGlobalJacobianCache<LeafNode>;
-
- public:
-
- using Element = typename Basis::LocalView::Element;
- using Coordinate = typename Element::Geometry::LocalCoordinate;
-
- FEFunctionJacobianForm(const Basis& basis, const TreePath& treePath) :
- basis_(basis),
- treePath_(treePath),
- leafTreePath_(Impl::leafTreePath<Tree>(treePath)),
- leafNode_(nullptr)
- {}
-
- auto quadratureRuleKey() const
- {
- return QuadratureRuleKey(leafNode_->finiteElement()).derivative();
- }
-
- void bindToLocalView(const LocalView& localView)
- {
- leafNode_ = & Dune::TypeTree::child(localView.tree(), leafTreePath_);
- }
-
- template<class Cache>
- void bindToCache(Cache& cache)
- {
- globalJacobianCache_ = &(cache.getGlobalJacobians(leafTreePath_));
- }
-
- const auto& treePath() const
- {
- return treePath_;
- }
-
- const auto& basis() const
- {
- return basis_;
- }
-
- auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
- {
- const auto& globalJacobians = (*globalJacobianCache_)[quadPointIndex];
- if constexpr(Node::isLeaf)
- return [&](const auto& args) {
- return globalJacobians[std::get<argIndex>(args)];
- };
- else if constexpr(Node::isPower and Node::ChildType::isLeaf)
- return [&](const auto& args) {
- using F = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
- 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);
- result[componentIndex] = globalJacobians[shapeFunctionIndex][0];
- return result;
- };
- }
-
- template<std::size_t i>
- auto childForm(Dune::index_constant<i> childIndex) const
- {
- auto childTP = Dune::TypeTree::push_back(treePath_, childIndex);
- return FEFunctionJacobianForm<Basis, decltype(childTP), argIndex>(basis_, childTP);
- }
-
- private:
- const Basis& basis_;
- const TreePath treePath_;
- Impl::LeafTreePath<Tree,TreePath> leafTreePath_;
- const LeafNode* leafNode_;
- mutable Buffer evaluationBuffer_;
- const GlobalJacobianCache* globalJacobianCache_;
- };
-
-
-
- // Linear map on FE-space that maps v onto its divergence
- template<class B, class TP, std::size_t argIndex>
- class FEFunctionDivergenceForm : public UnaryForm<FEFunctionDivergenceForm<B,TP,argIndex>, argIndex>
- {
- public:
- using Basis = B;
- using TreePath = TP;
- using LocalView = typename Basis::LocalView;
- using Tree = typename LocalView::Tree;
- using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
-
- private:
- using LeafNode = typename Dune::TypeTree::ChildForTreePath<Tree, Impl::LeafTreePath<Tree,TreePath>>;
-
- using Buffer = std::vector<typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::JacobianType>;
- using GlobalJacobianCache = typename ShapeFunctionCache<Tree>::template NodalGlobalJacobianCache<LeafNode>;
-
- public:
- using Element = typename Basis::LocalView::Element;
- using Coordinate = typename Element::Geometry::LocalCoordinate;
-
- FEFunctionDivergenceForm(const Basis& basis, const TreePath& treePath) :
- basis_(basis),
- treePath_(treePath),
- leafTreePath_(Impl::leafTreePath<Tree>(treePath)),
- leafNode_(nullptr)
- {}
-
- auto quadratureRuleKey() const
- {
- return QuadratureRuleKey(leafNode_->finiteElement()).derivative();
- }
-
- void bindToLocalView(const LocalView& localView)
- {
- leafNode_ = & Dune::TypeTree::child(localView.tree(), leafTreePath_);
- }
-
- template<class Cache>
- void bindToCache(Cache& cache)
- {
- globalJacobianCache_ = &(cache.getGlobalJacobians(leafTreePath_));
- }
-
- const auto& treePath() const
- {
- return treePath_;
- }
-
- const auto& basis() const
- {
- return basis_;
- }
-
- auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
- {
- const auto& globalJacobians = (*globalJacobianCache_)[quadPointIndex];
- if constexpr(Node::isLeaf)
- return [&](const auto& args) {
- using Field = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
- constexpr auto dimension = LeafNode::FiniteElement::Traits::LocalBasisType::Traits::dimDomain;
- auto div = Field{0};
- for(std::size_t i=0; i<dimension; ++i)
- div += globalJacobians[std::get<argIndex>(args)][i][i];
- return div;
- };
- else if constexpr(Node::isPower and Node::ChildType::isLeaf)
- return [&](const auto& args) {
- // Here we should compute the trace of the global Jacobian J
- // that we would need to build first. But since only the
- // componentIndex'th row of J is nonzero, we can simply
- // return the diagonal entry of this row which coincides
- // with the componentIndex'th entry of the gradient of
- // the respective scalar basis function.
- auto componentIndex = std::get<argIndex>(args) / leafNode_->size();
- auto shapeFunctionIndex = std::get<argIndex>(args) % leafNode_->size();
- return globalJacobians[shapeFunctionIndex][0][componentIndex];
- };
- }
-
- private:
- const Basis& basis_;
- const TreePath treePath_;
- Impl::LeafTreePath<Tree,TreePath> leafTreePath_;
- const LeafNode* leafNode_;
- mutable Buffer evaluationBuffer_;
- const GlobalJacobianCache* globalJacobianCache_;
- };
-
-
-
- // Generic product of multilinear forms. This combines
- // a k_1-linear map F_1 with a k_2-linear map F_2
- // into one (k_1+k_2)-linear map using the pointwise
- // bilinear map P on the range spaces.
- template<class P, class Form0, class Form1>
- class ProductForm : public Form
- {
- using Product = P;
-
- public:
-
- static constexpr std::size_t arity = Form0::arity + Form1::arity;
-
- using Element = typename Form0::Element;
- using Coordinate = typename Element::Geometry::LocalCoordinate;
-
- ProductForm(const Product& product, const Form0& form0, const Form1& form1) :
- product_(product),
- form0_(form0),
- form1_(form1)
- {}
-
- 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
- {
- if constexpr (std::is_same_v<decltype(form1_.testBasis()), void>)
- return form0_.testBasis();
- else
- return form1_.testBasis();
- }
-
- decltype(auto) trialBasis() const
- {
- if constexpr (std::is_same_v<decltype(form1_.trialBasis()), void>)
- return form0_.trialBasis();
- else
- return form1_.trialBasis();
- }
-
- decltype(auto) testTreePath() const
- {
- if constexpr (std::is_same_v<decltype(form1_.testTreePath()), void>)
- return form0_.testTreePath();
- else
- return form1_.testTreePath();
- }
-
- decltype(auto) trialTreePath() const
- {
- 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
- {
- auto result0 = form0_(x, quadPointIndex);
- auto result1 = form1_(x, quadPointIndex);
- return [=](const auto& args) {
- return product_(result0(args), result1(args));
- };
- }
-
- private:
- Product product_;
- Form0 form0_;
- Form1 form1_;
- };
-
-
-
- template<class T, class Form0>
- class TransformedForm : public Form
- {
- using Transformation = T;
-
- public:
- static constexpr std::size_t arity = Form0::arity;
-
- using Element = typename Form0::Element;
- using Coordinate = typename Element::Geometry::LocalCoordinate;
-
- TransformedForm(const Transformation& transformation, const Form0& form0) :
- transformation_(transformation),
- form0_(form0)
- {}
-
- 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
- {
- return form0_.testBasis();
- }
-
- decltype(auto) trialBasis() const
- {
- return form0_.trialBasis();
- }
-
- decltype(auto) testTreePath() const
- {
- return form0_.testTreePath();
- }
-
- decltype(auto) trialTreePath() const
- {
- return form0_.trialTreePath();
- }
-
- auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
- {
- auto result0 = form0_(x, quadPointIndex);
- return [=](const auto& args) {
- return transformation_(result0(args));
- };
- }
-
- auto subForm() const
- {
- return form0_;
- }
-
- private:
- Transformation transformation_;
- Form0 form0_;
- };
-
-
-
- template<class F>
- class TransposedForm
- : public TransformedForm<Dune::Fufem::Forms::Impl::TransposeOp, F>
- {
- using Base = TransformedForm<Dune::Fufem::Forms::Impl::TransposeOp, F>;
- public:
- using Base::arity;
-
- TransposedForm(const F& form) :
- Base(Dune::Fufem::Forms::Impl::TransposeOp(), form)
- {}
- };
-
-
-
- // Sum of multilinear forms of the same arity k gives another multilinear
- // of arity k.
- template<class Form0, class... Forms>
- class SumForm : public Form
- {
- public:
- static constexpr std::size_t arity = Form0::arity;
-
- using Element = typename Form0::Element;
- using Coordinate = typename Element::Geometry::LocalCoordinate;
-
- SumForm(const Form0& form0, const Forms&... forms) :
- forms_(form0, forms...)
- {}
-
- 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.bindToTrialLocalView(localView);
- });
- }
-
- template<class Cache>
- void bindToTestCache(Cache& cache)
- {
- Impl::forEachTupleEntry(forms_, [&](auto& form) {
- form.bindToTestCache(cache);
- });
- }
-
- template<class Cache>
- void bindToTrialCache(Cache& cache)
- {
- Impl::forEachTupleEntry(forms_, [&](auto& form) {
- 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_;
- }
-
- const auto& forms() const
- {
- return forms_;
- }
-
- private:
- std::tuple<Form0, Forms...> forms_;
- };
-
-
-
- template<class Form>
- auto makeSumForm(Form form)
- {
- return SumForm<Form>(form);
- }
-
- template<class... Forms>
- auto makeSumForm(SumForm<Forms...> sumForm)
- {
- return sumForm;
- }
-
-
-
- template<class Basis>
- auto testFunction(const Basis& basis)
- {
- using RootBasis = std::decay_t<decltype(basis.rootBasis())>;
- return FEFunctionForm<RootBasis, typename Basis::PrefixPath, 0>(basis.rootBasis(), basis.prefixPath());
- }
-
- template<class Basis, class... Args, std::enable_if_t<(sizeof...(Args)>0), int> = 0>
- auto testFunction(const Basis& basis, Args&&... args)
- {
- return testFunction(Functions::subspaceBasis(basis, args...));
- }
-
-
-
-
- template<class Basis>
- auto trialFunction(const Basis& basis)
- {
- using RootBasis = std::decay_t<decltype(basis.rootBasis())>;
- return FEFunctionForm<RootBasis, typename Basis::PrefixPath, 1>(basis.rootBasis(), basis.prefixPath());
- }
-
- template<class Basis, class... Args, std::enable_if_t<(sizeof...(Args)>0), int> = 0>
- auto trialFunction(const Basis& basis, Args&&... args)
- {
- return trialFunction(Functions::subspaceBasis(basis, args...));
- }
-
-
-
- template<class L, class R,
- std::enable_if_t<isForm_v<L> and isForm_v<R>, int> = 0>
- auto sum(const L& l, const R& r)
- {
- return SumForm<L, R>(l, r);
- }
-
- template<class L, class... Ri,
- std::enable_if_t<isForm_v<L>, int> = 0>
- auto sum(const L& l, const SumForm<Ri...>& r)
- {
- return std::apply([&](const auto&... ri) {
- return SumForm<L, Ri...>(l, ri...);
- }, r.forms());
- }
-
- template<class... Li, class R,
- std::enable_if_t<isForm_v<R>, int> = 0>
- auto sum(const SumForm<Li...>& l, const R& r)
- {
- return std::apply([&](const auto&... li) {
- return SumForm<Li..., R>(li..., r);
- }, l.forms());
- }
-
- template<class... Li, class... Rj>
- auto sum(const SumForm<Li...>& l, const SumForm<Rj...>& r)
- {
- return std::apply([&](const auto&...li) {
- return std::apply([&](const auto&...rj) {
- return SumForm<Li..., Rj...>(li..., rj...);
- }, r.forms());
- }, l.forms());
- }
-
-
-
- template<class Op, class F>
- auto transformedForm(const Op& op, const F& f)
- {
- return TransformedForm(op, f);
- }
-
- template<class Op, class... Fi>
- auto transformedForm(const Op& op, const SumForm<Fi...>& f)
- {
- return std::apply([&](const auto&...fi) {
- return (transformedForm(op, fi) + ...);
- }, f.forms());
- }
-
-
-
- template<class Op>
- class FormOperator
- {
- public:
- FormOperator(const Op& op) : op_(op) {}
-
- template<class Arg, std::enable_if_t<isForm_v<Arg>, int> = 0>
- auto operator()(Arg&& arg) const
- {
- return transformedForm(op_, std::forward<Arg>(arg));
- }
-
- template<class Arg, std::enable_if_t<not isForm_v<Arg>, int> = 0>
- auto operator()(Arg&& arg) const
- {
- return op_(std::forward<Arg>(arg));
- }
-
- private:
- Op op_;
- };
-
-
-
- // The following are the generic functions for producing
- // product forms from a binary operator and two other
- // terms. At least one of those terms should be a form.
-
- // 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,
- 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);
- }
-
- // Overload for one raw value and one form
- // 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,
- 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);
- }
-
- // Overload for one form and one raw value
- // 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,
- 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));
- }
-
- // Overload for one sum-form and one plain term
- // This will multiply out the product for the sum-form.
- template<class Op, class... Li, class R,
- class = std::void_t<decltype((productForm(std::declval<Op>(), std::declval<Li>(), std::declval<R>()),...))>>
- auto productForm(const Op& op, const SumForm<Li...>& l, const R& r)
- {
- return std::apply([&](const auto&... li) {
- return (productForm(op, li, r) + ...);
- }, l.forms());
- }
-
- // Overload for one plain term and one sum-form
- // This will multiply out the product for the sum-form.
- template<class Op, class L, class... Ri,
- class = std::void_t<decltype((productForm(std::declval<Op>(), std::declval<L>(), std::declval<Ri>()),...))>>
- auto productForm(const Op& op, const L& l, const SumForm<Ri...>& r)
- {
- return std::apply([&](const auto&...ri) {
- return (productForm(op, l, ri) + ...);
- }, r.forms());
- }
-
- // Overload for two sum-forms
- // This will multiply out the product for both sum-forms.
- template<class Op, class... Li, class... Rj,
- class = std::void_t<decltype((productForm(std::declval<Op>(), std::declval<Li>(), std::declval<SumForm<Rj...>>()),...))>>
- auto productForm(const Op& op, const SumForm<Li...>& l, const SumForm<Rj...>& r)
- {
- return std::apply([&](const auto&...li) {
- return (productForm(op, li, r) + ...);
- }, l.forms());
- }
-
-
-
- // This is the default implementation. The following specializations
- // are optimized versions for certain special cases.
- template<class L, class R,
- std::enable_if_t<isForm_v<L> or isForm_v<R>, int> = 0,
- class = std::void_t<decltype(productForm(std::declval<Impl::DotOp>(), std::declval<L>(), std::declval<R>()))>>
- auto dot(const L& l, const R& r)
- {
- return productForm(Dune::Fufem::Forms::Impl::DotOp(), l, r);
- }
-
- template<class B, class TP, std::size_t argIndex_L, std::size_t argIndex_R>
- auto dot(const FEFunctionForm<B,TP,argIndex_L>& l, const FEFunctionForm<B,TP,argIndex_R>& r);
-
- template<class B, class TP, std::size_t argIndex_L, std::size_t argIndex_R>
- auto dot(const FEFunctionJacobianForm<B,TP,argIndex_L>& l, const FEFunctionJacobianForm<B,TP,argIndex_R>& r);
-
- template<class L, class R>
- auto dot(const TransposedForm<L>& l, const TransposedForm<R>& r)
- {
- using namespace Dune::Indices;
- return Dune::Fufem::Forms::dot(l.subForm(), r.subForm());
- }
-
- template<class B, class TP, std::size_t argIndex_L, std::size_t argIndex_R>
- auto dot(const FEFunctionForm<B,TP,argIndex_L>& l, const FEFunctionForm<B,TP,argIndex_R>& r)
- {
- using Node = typename FEFunctionForm<B,TP,argIndex_L>::Node;
- if constexpr (Node::isPower)
- {
- return unpackIntegerSequence(
- [&](auto... i) {
- return (Dune::Fufem::Forms::dot(l.childForm(i), r.childForm(i)) + ...);
- },
- std::make_index_sequence<Node::degree()>{});
- }
- else
- return productForm(Dune::Fufem::Forms::Impl::DotOp(), l, r);
- }
-
- template<class B, class TP, std::size_t argIndex_L, std::size_t argIndex_R>
- auto dot(const FEFunctionJacobianForm<B,TP,argIndex_L>& l, const FEFunctionJacobianForm<B,TP,argIndex_R>& r)
- {
- using Node = typename FEFunctionJacobianForm<B,TP,argIndex_L>::Node;
- if constexpr (Node::isPower)
- {
- return unpackIntegerSequence(
- [&](auto... i) {
- return (Dune::Fufem::Forms::dot(l.childForm(i), r.childForm(i)) + ...);
- },
- std::make_index_sequence<Node::degree()>{});
- }
- else
- return productForm(Dune::Fufem::Forms::Impl::DotOp(), l, r);
- }
-
-
-
- // This overload is recursively constrained to cases where `operator*` makes sense.
- // Otherwise Dune::dot(a,b) may be activated for forms because it relies on availability
- // of a*b. This is a problem, because it leads to ambiguity with Forms::dot(a,b).
- template<class L, class R,
- std::enable_if_t<isForm_v<L> or isForm_v<R>, int> = 0,
- class = std::void_t<decltype(productForm(std::declval<Impl::MultOp>(), std::declval<L>(), std::declval<R>()))>>
- auto operator* (const L& l, const R& r)
- {
- return productForm(Impl::MultOp{}, l, r);
- }
-
-
-
- template<class F,
- std::enable_if_t<isForm_v<F>, int> = 0>
- auto operator- (const F& f)
- {
- auto op = [](const auto& x) {
- return -x;
- };
- return transformedForm(op, f);
- }
-
- template<class L, class R,
- std::enable_if_t<isForm_v<L> and isForm_v<R>, int> = 0>
- auto operator+ (const L& l, const R& r)
- {
- static_assert(L::arity==R::arity, "The forms passed to operator+ do not have the same arity");
- return sum(l, r);
- }
-
- template<class L, class R,
- std::enable_if_t<isForm_v<L> and isForm_v<R>, int> = 0>
- auto operator- (const L& l, const R& r)
- {
- static_assert(L::arity==R::arity, "The forms passed to operator- do not have the same arity");
- return sum(l, -r);
- }
-
-
-
- template<class F,
- std::enable_if_t<isForm_v<F>, int>>
- auto transpose(const F& f)
- {
- return TransposedForm(f);
- }
+namespace Dune::Fufem::Forms {
template<class Form, bool isSemiAffine=true>
class IntegratedLinearForm
{
- using TestRootBasis = std::decay_t<decltype(std::declval<Form>().testBasis())>;
+ using TestRootBasis = std::decay_t<decltype(std::declval<Form>().template basis<0>())>;
using TestTree = typename TestRootBasis::LocalView::Tree;
@@ -1322,7 +41,7 @@ namespace Impl {
IntegratedLinearForm(const Form& sumForm) :
sumForm_(sumForm),
- testRootBasis_(sumForm_.testBasis())
+ testRootBasis_(sumForm_.template basis<0>())
{}
template<class LocalVector, class TestLocalView>
@@ -1333,13 +52,13 @@ namespace Impl {
const auto& geometry = element.geometry();
sumForm_.bindToElement(testLocalView.element());
- sumForm_.bindToTestLocalView(testLocalView);
+ sumForm_.bindToLocalView(testLocalView);
auto& cacheForRule = cache_[sumForm_.quadratureRuleKey()];
cacheForRule.bind(testLocalView.tree());
cacheForRule.invalidate(isSemiAffine);
- sumForm_.bindToTestCache(cacheForRule);
+ sumForm_.bindToCache(cacheForRule);
const auto& quad = cacheForRule.rule();
@@ -1352,7 +71,7 @@ namespace Impl {
Impl::forEachTupleEntry(sumForm_.forms(), [&](auto& form) {
auto evaluatedForm = form(quadPoint.position(), k);
- const auto& testPrefix = form.testTreePath();
+ const auto& testPrefix = form.template treePath<0>();
const auto& testNode = Dune::TypeTree::child(testLocalView.tree(), testPrefix);
for (std::size_t i=0; i<testNode.size(); ++i)
{
@@ -1379,8 +98,8 @@ namespace Impl {
template<class Form, bool isSemiAffine=true>
class IntegratedBilinearForm
{
- using TestRootBasis = std::decay_t<decltype(std::declval<Form>().testBasis())>;
- using AnsatzRootBasis = std::decay_t<decltype(std::declval<Form>().trialBasis())>;
+ using TestRootBasis = std::decay_t<decltype(std::declval<Form>().template basis<0>())>;
+ using AnsatzRootBasis = std::decay_t<decltype(std::declval<Form>().template basis<1>())>;
using TestTree = typename TestRootBasis::LocalView::Tree;
@@ -1389,8 +108,8 @@ namespace Impl {
IntegratedBilinearForm(const Form& sumForm) :
sumForm_(sumForm),
- testRootBasis_(sumForm_.testBasis()),
- ansatzRootBasis_(sumForm_.trialBasis())
+ testRootBasis_(sumForm_.template basis<0>()),
+ ansatzRootBasis_(sumForm_.template basis<1>())
{}
template<class LocalMatrix, class TestLocalView, class AnsatzLocalView>
@@ -1402,15 +121,13 @@ namespace Impl {
const auto& geometry = element.geometry();
sumForm_.bindToElement(testLocalView.element());
- sumForm_.bindToTestLocalView(testLocalView);
- sumForm_.bindToTrialLocalView(ansatzLocalView);
+ sumForm_.bindToLocalView(testLocalView, ansatzLocalView);
auto& cacheForRule = cache_[sumForm_.quadratureRuleKey()];
cacheForRule.bind(testLocalView.tree());
cacheForRule.invalidate(isSemiAffine);
- sumForm_.bindToTestCache(cacheForRule);
- sumForm_.bindToTrialCache(cacheForRule);
+ sumForm_.bindToCache(cacheForRule, cacheForRule);
const auto& quad = cacheForRule.rule();
@@ -1423,9 +140,9 @@ namespace Impl {
Impl::forEachTupleEntry(sumForm_.forms(), [&](auto& form) {
auto evaluatedForm = form(quadPoint.position(), k);
- const auto& testPrefix = form.testTreePath();
+ const auto& testPrefix = form.template treePath<0>();
const auto& testNode = Dune::TypeTree::child(testLocalView.tree(), testPrefix);
- const auto& ansatzPrefix = form.trialTreePath();
+ const auto& ansatzPrefix = form.template treePath<1>();
const auto& ansatzTree = Dune::TypeTree::child(ansatzLocalView.tree(), ansatzPrefix);
for (std::size_t i=0; i<testNode.size(); ++i)
{
@@ -1457,7 +174,7 @@ namespace Impl {
template<class Form, class Patch, bool isSemiAffine=true>
class IntegratedBoundaryLinearForm
{
- using TestRootBasis = std::decay_t<decltype(std::declval<Form>().testBasis())>;
+ using TestRootBasis = std::decay_t<decltype(std::declval<Form>().template basis<0>())>;
using TestTree = typename TestRootBasis::LocalView::Tree;
@@ -1468,7 +185,7 @@ namespace Impl {
IntegratedBoundaryLinearForm(const Form& sumForm, const Patch& patch) :
sumForm_(sumForm),
- testRootBasis_(sumForm_.testBasis()),
+ testRootBasis_(sumForm_.template basis<0>()),
patch_(patch)
{}
@@ -1491,13 +208,13 @@ namespace Impl {
const auto& geometry = intersection.geometry();
sumForm_.bindToElement(testLocalView.element());
- sumForm_.bindToTestLocalView(testLocalView);
+ sumForm_.bindToLocalView(testLocalView);
auto& cacheForRule = cache_[std::pair(sumForm_.quadratureRuleKey(), facet)];
cacheForRule.bind(testLocalView.tree());
cacheForRule.invalidate(isSemiAffine);
- sumForm_.bindToTestCache(cacheForRule);
+ sumForm_.bindToCache(cacheForRule);
const auto& quad = cacheForRule.rule();
@@ -1511,7 +228,7 @@ namespace Impl {
Impl::forEachTupleEntry(sumForm_.forms(), [&](auto& form) {
auto evaluatedForm = form(quadPoint.position(), k);
- const auto& testPrefix = form.testTreePath();
+ const auto& testPrefix = form.template treePath<0>();
const auto& testNode = Dune::TypeTree::child(testLocalView.tree(), testPrefix);
for (std::size_t i=0; i<testNode.size(); ++i)
{
@@ -1590,7 +307,7 @@ namespace Impl {
std::enable_if_t<isForm_v<Form>, int> = 0>
auto integrate(Form form)
{
- auto sumForm = makeSumForm(form);
+ auto sumForm = SumForm(form);
if constexpr(Form::arity==1)
return IntegratedLinearForm(sumForm);
else if constexpr(Form::arity==2)
@@ -1601,7 +318,7 @@ namespace Impl {
std::enable_if_t<isForm_v<Form>, int> = 0>
auto integrate(Form form, const BoundaryPatch<GridView>& patch)
{
- auto sumForm = makeSumForm(form);
+ auto sumForm = SumForm(form);
if constexpr(Form::arity==1)
return IntegratedBoundaryLinearForm(sumForm, patch);
}
@@ -1614,7 +331,7 @@ namespace Impl {
std::enable_if_t<isForm_v<Form>, int> = 0>
auto integrate(Form form, NonAffineFamiliyTag)
{
- auto sumForm = makeSumForm(form);
+ auto sumForm = SumForm(form);
using SumForm = decltype(sumForm);
if constexpr(Form::arity==1)
return IntegratedLinearForm<SumForm, false>(sumForm);
@@ -1626,7 +343,7 @@ namespace Impl {
std::enable_if_t<isForm_v<Form>, int> = 0>
auto integrate(Form form, const BoundaryPatch<GridView>& patch, NonAffineFamiliyTag)
{
- auto sumForm = makeSumForm(form);
+ auto sumForm = SumForm(form);
using SumForm = decltype(sumForm);
using Patch = BoundaryPatch<GridView>;
if constexpr(Form::arity==1)
diff --git a/dune/fufem-forms/fufem-forms.hh b/dune/fufem-forms/fufem-forms.hh
deleted file mode 100644
index 6262192..0000000
--- a/dune/fufem-forms/fufem-forms.hh
+++ /dev/null
@@ -1,6 +0,0 @@
-#ifndef DUNE_FUFEM_FORMS_HH
-#define DUNE_FUFEM_FORMS_HH
-
-// add your classes here
-
-#endif // DUNE_FUFEM_FORMS_HH
diff --git a/dune/fufem-forms/nullaryforms.hh b/dune/fufem-forms/nullaryforms.hh
new file mode 100644
index 0000000..d2ce26a
--- /dev/null
+++ b/dune/fufem-forms/nullaryforms.hh
@@ -0,0 +1,74 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FUFEM_FORMS_NULLARYFORMS_HH
+#define DUNE_FUFEM_FORMS_NULLARYFORMS_HH
+
+#include <type_traits>
+#include <tuple>
+#include <utility>
+
+#include <dune/common/typetraits.hh>
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/matrixindexset.hh>
+
+#include <dune/fufem-forms/baseclass.hh>
+
+
+
+namespace Dune::Fufem::Forms {
+
+
+
+ // Zero-linear map induced by fixed coefficient function
+ template<class F>
+ class Coefficient : public Form
+ {
+ using LocalFunction = std::decay_t<decltype(localFunction(std::declval<F&>()))>;
+
+ public:
+ static constexpr std::size_t arity = 0;
+
+ using Element = typename F::EntitySet::Element;
+ using Coordinate = typename Element::Geometry::LocalCoordinate;
+
+ Coefficient(const F& f) :
+ f_(f),
+ localF_(localFunction(f))
+ {}
+
+ Coefficient(const F& f, std::size_t order) :
+ f_(f),
+ localF_(localFunction(f))
+ {}
+
+ auto quadratureRuleKey() const
+ {
+ return QuadratureRuleKey(Element::dimension, 0);
+ }
+
+ void bindToElement(const Element& element)
+ {
+ localF_.bind(element);
+ }
+
+ auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
+ {
+ auto value = localF_(x);
+ return [value](const auto& args) {
+ return value;
+ };
+ }
+
+ private:
+ F f_;
+ LocalFunction localF_;
+ };
+
+
+
+} // namespace Dune::Fufem::Forms
+
+
+
+#endif // DUNE_FUFEM_FORMS_NULLARYFORMS_HH
diff --git a/dune/fufem-forms/operators.hh b/dune/fufem-forms/operators.hh
new file mode 100644
index 0000000..92be0c5
--- /dev/null
+++ b/dune/fufem-forms/operators.hh
@@ -0,0 +1,213 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FUFEM_FORMS_OPERATORS_HH
+#define DUNE_FUFEM_FORMS_OPERATORS_HH
+
+#include <cstddef>
+#include <type_traits>
+#include <utility>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/istl/scaledidmatrix.hh>
+
+#include <dune/fufem-forms/baseclass.hh>
+#include <dune/fufem-forms/transformedform.hh>
+#include <dune/fufem-forms/unaryforms.hh>
+#include <dune/fufem-forms/productform.hh>
+
+
+
+namespace Dune::Fufem::Forms::Impl {
+
+
+
+ struct DotOp {
+
+ template<class K1, class K2, int n>
+ auto operator()(const Dune::FieldVector<K1,n>& x, const Dune::FieldVector<K2,n>& y) const
+ {
+ using K = typename PromotionTraits<K1,K2>::PromotedType;
+ auto result = K(0);
+ for(std::size_t i=0; i<n; ++i)
+ result += x[i]*y[i];
+ return result;
+ }
+
+ template<class K1, class K2, int n, int m>
+ auto operator()(const Dune::FieldMatrix<K1,n,m>& x, const Dune::FieldMatrix<K2,n,m>& y) const
+ {
+ using K = typename PromotionTraits<K1,K2>::PromotedType;
+ auto result = K(0);
+ for(std::size_t i=0; i<n; ++i)
+ for(std::size_t j=0; j<m; ++j)
+ result += x[i][j]*y[i][j];
+ return result;
+ }
+
+ template<class K1, class K2, int n>
+ auto operator()(const Dune::FieldMatrix<K1,n,1>& x, const Dune::FieldVector<K2,n>& y) const
+ {
+ using K = typename PromotionTraits<K1,K2>::PromotedType;
+ auto result = K(0);
+ for(std::size_t i=0; i<n; ++i)
+ result += x[i][0]*y[i];
+ return result;
+ }
+
+ template<class K1, class K2, int n>
+ auto operator()(const Dune::FieldVector<K1,n>& x, const Dune::FieldMatrix<K2,n,1>& y) const
+ {
+ using K = typename PromotionTraits<K1,K2>::PromotedType;
+ auto result = K(0);
+ for(std::size_t i=0; i<n; ++i)
+ result += x[i]*y[i][0];
+ return result;
+ }
+
+ template<class K1, class K2, int n>
+ auto operator()(const Dune::ScaledIdentityMatrix<K1,n>& x, const Dune::FieldMatrix<K2,n,n>& y) const
+ {
+ using K = typename PromotionTraits<K1,K2>::PromotedType;
+ auto result = K(0);
+ for(std::size_t i=0; i<x.N(); ++i)
+ result += y[i][i];
+ return x.scalar()*result;
+ }
+
+ template<class K1, class K2, int n>
+ auto operator()(const Dune::FieldMatrix<K1,n,n>& x, const Dune::ScaledIdentityMatrix<K2,n>& y) const
+ {
+ using K = typename PromotionTraits<K1,K2>::PromotedType;
+ auto result = K(0);
+ for(std::size_t i=0; i<x.N(); ++i)
+ result += x[i][i];
+ return result*y.scalar();
+ }
+
+ template<class K1, class K2, int n>
+ auto operator()(const Dune::ScaledIdentityMatrix<K1,n>& x, const Dune::ScaledIdentityMatrix<K2,n>& y) const
+ {
+ return n*x.scalar()*y.scalar();
+ }
+
+ };
+
+ struct MultOp {
+
+ template<class X, class Y,
+ class = std::void_t<decltype(std::declval<X>()*std::declval<Y>())>>
+ auto operator()(const X& x, const Y& y) const
+ {
+ return x*y;
+ }
+
+ template<class K1, class K2, int n, int m>
+ auto operator()(const Dune::FieldMatrix<K1,n,m>& x, const Dune::FieldVector<K2,m>& y) const
+ {
+ using K = typename PromotionTraits<K1,K2>::PromotedType;
+ Dune::FieldVector<K,n> result;
+ x.mv(y, result);
+ return result;
+ }
+
+ template<class K, class Y>
+ auto operator()(const Dune::FieldMatrix<K,1,1>& x, const Y& y) const
+ {
+ return x[0][0]*y;
+ }
+
+ template<class K, class Y>
+ auto operator()(const Dune::FieldVector<K,1>& x, const Y& y) const
+ {
+ return x[0]*y;
+ }
+
+ };
+
+
+
+} // namespace Impl
+
+
+
+namespace Dune::Fufem::Forms {
+
+
+
+ // This is the default implementation. The following specializations
+ // are optimized versions for certain special cases.
+ template<class L, class R,
+ std::enable_if_t<isForm_v<L> or isForm_v<R>, int> = 0,
+ class = std::void_t<decltype(productForm(std::declval<Impl::DotOp>(), std::declval<L>(), std::declval<R>()))>>
+ auto dot(const L& l, const R& r)
+ {
+ return productForm(Dune::Fufem::Forms::Impl::DotOp(), l, r);
+ }
+
+ template<class B, class TP, std::size_t argIndex_L, std::size_t argIndex_R>
+ auto dot(const FEFunctionForm<B,TP,argIndex_L>& l, const FEFunctionForm<B,TP,argIndex_R>& r);
+
+ template<class B, class TP, std::size_t argIndex_L, std::size_t argIndex_R>
+ auto dot(const FEFunctionJacobianForm<B,TP,argIndex_L>& l, const FEFunctionJacobianForm<B,TP,argIndex_R>& r);
+
+ template<class L, class R>
+ auto dot(const TransposedForm<L>& l, const TransposedForm<R>& r)
+ {
+ std::cout << "eliminating transpose" << std::endl;
+ using namespace Dune::Indices;
+ return Dune::Fufem::Forms::dot(l.baseForm(), r.baseForm());
+ }
+
+ template<class B, class TP, std::size_t argIndex_L, std::size_t argIndex_R>
+ auto dot(const FEFunctionForm<B,TP,argIndex_L>& l, const FEFunctionForm<B,TP,argIndex_R>& r)
+ {
+ using Node = typename FEFunctionForm<B,TP,argIndex_L>::Node;
+ if constexpr (Node::isPower)
+ {
+ return unpackIntegerSequence(
+ [&](auto... i) {
+ return (Dune::Fufem::Forms::dot(l.childForm(i), r.childForm(i)) + ...);
+ },
+ std::make_index_sequence<Node::degree()>{});
+ }
+ else
+ return productForm(Dune::Fufem::Forms::Impl::DotOp(), l, r);
+ }
+
+ template<class B, class TP, std::size_t argIndex_L, std::size_t argIndex_R>
+ auto dot(const FEFunctionJacobianForm<B,TP,argIndex_L>& l, const FEFunctionJacobianForm<B,TP,argIndex_R>& r)
+ {
+ using Node = typename FEFunctionJacobianForm<B,TP,argIndex_L>::Node;
+ if constexpr (Node::isPower)
+ {
+ return unpackIntegerSequence(
+ [&](auto... i) {
+ return (Dune::Fufem::Forms::dot(l.childForm(i), r.childForm(i)) + ...);
+ },
+ std::make_index_sequence<Node::degree()>{});
+ }
+ else
+ return productForm(Dune::Fufem::Forms::Impl::DotOp(), l, r);
+ }
+
+
+
+ // This overload is recursively constrained to cases where `operator*` makes sense.
+ // Otherwise Dune::dot(a,b) may be activated for forms because it relies on availability
+ // of a*b. This is a problem, because it leads to ambiguity with Forms::dot(a,b).
+ template<class L, class R,
+ std::enable_if_t<isForm_v<L> or isForm_v<R>, int> = 0,
+ class = std::void_t<decltype(productForm(std::declval<Impl::MultOp>(), std::declval<L>(), std::declval<R>()))>>
+ auto operator* (const L& l, const R& r)
+ {
+ return productForm(Impl::MultOp{}, l, r);
+ }
+
+
+
+} // namespace Dune::Fufem::Forms
+
+
+#endif // DUNE_FUFEM_FORMS_OPERATORS_HH
diff --git a/dune/fufem-forms/productform.hh b/dune/fufem-forms/productform.hh
new file mode 100644
index 0000000..14ecebc
--- /dev/null
+++ b/dune/fufem-forms/productform.hh
@@ -0,0 +1,248 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FUFEM_FORMS_PRODUCTFORM_HH
+#define DUNE_FUFEM_FORMS_PRODUCTFORM_HH
+
+#include <cstddef>
+#include <type_traits>
+#include <tuple>
+#include <utility>
+
+#include <dune/fufem-forms/baseclass.hh>
+#include <dune/fufem-forms/transformedform.hh>
+#include <dune/fufem-forms/sumform.hh>
+
+
+
+namespace Dune::Fufem::Forms::Impl {
+
+
+ template<class Op>
+ class TransposedBinaryOp
+ {
+ public:
+ TransposedBinaryOp(const Op& op) : op_(op) {}
+
+ template<class X, class Y>
+ decltype(auto) operator()(X&& x, Y&& y) const
+ {
+ return op_(std::forward<X>(x), std::forward<Y>(y));
+ }
+
+ private:
+ Op op_;
+ };
+
+} // namespace Dune::Fufem::Forms::Impl
+
+namespace Dune::Fufem::Forms {
+
+
+ template<class C, class BaseForm, class NullaryForm>
+ class ProductWithNullaryForm : public BaseForm
+ {
+ using Contraction = C;
+
+ public:
+
+ static constexpr std::size_t arity = NullaryForm::arity + BaseForm::arity;
+
+ using Element = typename BaseForm::Element;
+ using Coordinate = typename BaseForm::Coordinate;
+
+ ProductWithNullaryForm(const Contraction& product, const BaseForm& form0, const NullaryForm& form1) :
+ BaseForm(form0),
+ product_(product),
+ form1_(form1)
+ {}
+
+ auto quadratureRuleKey() const
+ {
+ return BaseForm::quadratureRuleKey().product(form1_.quadratureRuleKey());
+ }
+
+ void bindToElement(const Element& element)
+ {
+ form1_.bindToElement(element);
+ if constexpr (BaseForm::arity==0)
+ BaseForm::bindToElement(element);
+ }
+
+ auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
+ {
+ auto result0 = BaseForm::operator()(x, quadPointIndex);
+ auto result1 = form1_(x, quadPointIndex);
+ return [=](const auto& args) {
+ return product_(result0(args), result1(args));
+ };
+ }
+
+ private:
+ Contraction product_;
+ NullaryForm form1_;
+ };
+
+
+
+ // Generic product of multilinear forms. This combines
+ // a k_1-linear map F_1 with a k_2-linear map F_2
+ // into one (k_1+k_2)-linear map using the pointwise
+ // bilinear map P on the range spaces.
+ template<class P, class Form0, class Form1>
+ class UnaryUnaryProductForm : public Dune::Fufem::Forms::Form
+ {
+ using Contraction = P;
+
+ public:
+
+ static constexpr std::size_t arity = Form0::arity + Form1::arity;
+
+ using Element = typename Form0::Element;
+ using Coordinate = typename Element::Geometry::LocalCoordinate;
+
+ UnaryUnaryProductForm(const Contraction& product, const Form0& form0, const Form1& form1) :
+ product_(product),
+ form0_(form0),
+ form1_(form1)
+ {}
+
+ auto quadratureRuleKey() const
+ {
+ return form0_.quadratureRuleKey().product(form1_.quadratureRuleKey());
+ }
+
+ void bindToElement(const Element& element)
+ {
+ form0_.bindToElement(element);
+ form1_.bindToElement(element);
+ }
+
+ template<class LocalView0, class LocalView1>
+ void bindToLocalView(const LocalView0& localView0, const LocalView1& localView1)
+ {
+ form0_.bindToLocalView(localView0);
+ form1_.bindToLocalView(localView1);
+ }
+
+ template<class Cache0, class Cache1>
+ void bindToCache(Cache0& cache0, Cache1& cache1)
+ {
+ form0_.bindToCache(cache0);
+ form1_.bindToCache(cache1);
+ }
+
+ template<int k>
+ decltype(auto) basis() const
+ {
+ if constexpr(k==0)
+ return form0_.template basis<0>();
+ if constexpr(k==1)
+ return form1_.template basis<0>();
+ }
+
+ template<int k>
+ decltype(auto) treePath() const
+ {
+ if constexpr(k==0)
+ return form0_.template treePath<0>();
+ if constexpr(k==1)
+ return form1_.template treePath<0>();
+ }
+
+ auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
+ {
+ auto result0 = form0_(x, quadPointIndex);
+ auto result1 = form1_(x, quadPointIndex);
+ return [=](const auto& args) {
+ return product_(result0(args), result1(args));
+ };
+ }
+
+ private:
+ Contraction product_;
+ Form0 form0_;
+ Form1 form1_;
+ };
+
+
+
+ // The following are the generic functions for producing
+ // product forms from a binary operator and two other
+ // terms. At least one of those terms should be a form.
+
+ // 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,
+ 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)
+ {
+ static_assert(L::arity+R::arity<=2, "Trying to construct multi-linear operator with rank>2");
+ if constexpr(R::arity==0)
+ return ProductWithNullaryForm(op, l, r);
+ else if constexpr(L::arity==0)
+ return ProductWithNullaryForm(Impl::TransposedBinaryOp(op), r, l);
+ else if constexpr((L::argIndex==0) and (R::argIndex==1))
+ return UnaryUnaryProductForm(op, l, r);
+ else if constexpr((L::argIndex==1) and (R::argIndex==0))
+ return UnaryUnaryProductForm(Impl::TransposedBinaryOp(op), r, l);
+ }
+
+ // Overload for one raw value and one form
+ // 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,
+ 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 transformRange(r, [op, l](auto&& r_value) { return op(l, r_value); });
+ }
+
+ // Overload for one form and one raw value
+ // 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,
+ 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 transformRange(l, [op, r](auto&& l_value) { return op(l_value, r); });
+ }
+
+ // Overload for one sum-form and one plain term
+ // This will multiply out the product for the sum-form.
+ template<class Op, class... Li, class R,
+ class = std::void_t<decltype((productForm(std::declval<Op>(), std::declval<Li>(), std::declval<R>()),...))>>
+ auto productForm(const Op& op, const SumForm<Li...>& l, const R& r)
+ {
+ return std::apply([&](const auto&... li) {
+ return (productForm(op, li, r) + ...);
+ }, l.forms());
+ }
+
+ // Overload for one plain term and one sum-form
+ // This will multiply out the product for the sum-form.
+ template<class Op, class L, class... Ri,
+ class = std::void_t<decltype((productForm(std::declval<Op>(), std::declval<L>(), std::declval<Ri>()),...))>>
+ auto productForm(const Op& op, const L& l, const SumForm<Ri...>& r)
+ {
+ return std::apply([&](const auto&...ri) {
+ return (productForm(op, l, ri) + ...);
+ }, r.forms());
+ }
+
+ // Overload for two sum-forms
+ // This will multiply out the product for both sum-forms.
+ template<class Op, class... Li, class... Rj,
+ class = std::void_t<decltype((productForm(std::declval<Op>(), std::declval<Li>(), std::declval<SumForm<Rj...>>()),...))>>
+ auto productForm(const Op& op, const SumForm<Li...>& l, const SumForm<Rj...>& r)
+ {
+ return std::apply([&](const auto&...li) {
+ return (productForm(op, li, r) + ...);
+ }, l.forms());
+ }
+
+
+
+} // namespace Dune::Fufem::Forms
+
+
+#endif // DUNE_FUFEM_FORMS_PRODUCTFORM_HH
diff --git a/dune/fufem-forms/sumform.hh b/dune/fufem-forms/sumform.hh
new file mode 100644
index 0000000..089e564
--- /dev/null
+++ b/dune/fufem-forms/sumform.hh
@@ -0,0 +1,156 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FUFEM_FORMS_SUMFORM_HH
+#define DUNE_FUFEM_FORMS_SUMFORM_HH
+
+#include <cstddef>
+#include <type_traits>
+#include <tuple>
+#include <utility>
+
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
+
+#include <dune/fufem-forms/baseclass.hh>
+
+
+
+namespace Dune::Fufem::Forms {
+
+
+
+ // Sum of multilinear forms of the same arity k gives another multilinear
+ // of arity k.
+ template<class Form0, class... Forms>
+ class SumForm : public Dune::Fufem::Forms::Form
+ {
+ public:
+ static constexpr std::size_t arity = Form0::arity;
+
+ using Element = typename Form0::Element;
+ using Coordinate = typename Element::Geometry::LocalCoordinate;
+
+ SumForm(Form0 form0, Forms... forms) :
+ forms_(form0, forms...)
+ {}
+
+ 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... LocalViews>
+ void bindToLocalView(const LocalViews&... localViews)
+ {
+ Impl::forEachTupleEntry(forms_, [&](auto& form) {
+ form.bindToLocalView(localViews...);
+ });
+ }
+
+ template<class... Caches>
+ void bindToCache(Caches&... caches)
+ {
+ Impl::forEachTupleEntry(forms_, [&](auto& form) {
+ form.bindToCache(caches...);
+ });
+ }
+
+ template<int k>
+ decltype(auto) basis() const
+ {
+ return std::get<0>(forms_).template basis<k>();
+ }
+
+ auto& forms()
+ {
+ return forms_;
+ }
+
+ const auto& forms() const
+ {
+ return forms_;
+ }
+
+ private:
+ std::tuple<Form0, Forms...> forms_;
+ };
+
+
+
+// template<class... Forms>
+// SumForm(Forms...) -> SumForm<Forms...>;
+
+// template<class... Forms>
+// SumForm(SumForm<Forms...>) -> SumForm<Forms...>;
+
+
+
+ template<class L, class R,
+ std::enable_if_t<isForm_v<L> and isForm_v<R>, int> = 0>
+ auto sum(const L& l, const R& r)
+ {
+ return SumForm<L, R>(l, r);
+ }
+
+ template<class L, class... Ri,
+ std::enable_if_t<isForm_v<L>, int> = 0>
+ auto sum(const L& l, const SumForm<Ri...>& r)
+ {
+ return std::apply([&](const auto&... ri) {
+ return SumForm<L, Ri...>(l, ri...);
+ }, r.forms());
+ }
+
+ template<class... Li, class R,
+ std::enable_if_t<isForm_v<R>, int> = 0>
+ auto sum(const SumForm<Li...>& l, const R& r)
+ {
+ return std::apply([&](const auto&... li) {
+ return SumForm<Li..., R>(li..., r);
+ }, l.forms());
+ }
+
+ template<class... Li, class... Rj>
+ auto sum(const SumForm<Li...>& l, const SumForm<Rj...>& r)
+ {
+ return std::apply([&](const auto&...li) {
+ return std::apply([&](const auto&...rj) {
+ return SumForm<Li..., Rj...>(li..., rj...);
+ }, r.forms());
+ }, l.forms());
+ }
+
+
+
+ template<class L, class R,
+ std::enable_if_t<isForm_v<L> and isForm_v<R>, int> = 0>
+ auto operator+ (const L& l, const R& r)
+ {
+ static_assert(L::arity==R::arity, "The forms passed to operator+ do not have the same arity");
+ return sum(l, r);
+ }
+
+ template<class L, class R,
+ std::enable_if_t<isForm_v<L> and isForm_v<R>, int> = 0>
+ auto operator- (const L& l, const R& r)
+ {
+ static_assert(L::arity==R::arity, "The forms passed to operator- do not have the same arity");
+ return sum(l, -r);
+ }
+
+
+
+} // namespace Dune::Fufem::Forms
+
+
+#endif // DUNE_FUFEM_FORMS_SUMFORM_HH
diff --git a/dune/fufem-forms/transformedform.hh b/dune/fufem-forms/transformedform.hh
new file mode 100644
index 0000000..e240492
--- /dev/null
+++ b/dune/fufem-forms/transformedform.hh
@@ -0,0 +1,158 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FUFEM_FORMS_TRANSFORMEDFORM_HH
+#define DUNE_FUFEM_FORMS_TRANSFORMEDFORM_HH
+
+#include <cstddef>
+#include <type_traits>
+#include <tuple>
+#include <utility>
+
+
+
+namespace Dune::Fufem::Forms {
+
+
+
+ template<class BaseForm, class Op>
+ class TransformedForm : public BaseForm
+ {
+ using Transformation = Op;
+
+ public:
+
+ using Element = typename BaseForm::Element;
+ using Coordinate = typename BaseForm::Coordinate;
+
+ TransformedForm(const BaseForm& form) :
+ BaseForm(form),
+ transformation_()
+ {}
+
+ TransformedForm(const BaseForm& form, const Transformation& transformation) :
+ BaseForm(form),
+ transformation_(transformation)
+ {}
+
+ auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
+ {
+ auto y = BaseForm::operator()(x, quadPointIndex);
+ return [=](const auto& args) {
+ return transformation_(y(args));
+ };
+ }
+
+ auto baseForm() const
+ {
+ return (const BaseForm&)(*this);
+ }
+
+ private:
+ Transformation transformation_;
+ };
+
+
+
+ template<class F, class Op>
+ auto transformRange(const F& f, const Op& op)
+ {
+ return TransformedForm(f, op);
+ }
+
+
+ template<class Form0, class... Forms>
+ class SumForm;
+
+ template<class Op, class... Fi>
+ auto transformRange(const Dune::Fufem::Forms::SumForm<Fi...>& f, const Op& op)
+ {
+ return std::apply([&](const auto&...fi) {
+ return Dune::Fufem::Forms::SumForm(transformRange(fi, op) ...);
+ }, f.forms());
+ }
+
+
+
+ template<class Op>
+ class RangeOperator
+ {
+ public:
+ RangeOperator(const Op& op) : op_(op) {}
+
+ template<class F, class... Args,
+ std::enable_if_t<isForm_v<F>, int> = 0>
+ auto operator()(F&& form, Args... args) const
+ {
+ auto captureOp = [op=op_, argTuple = std::tuple(args...)](auto&& arg) {
+ return std::apply([&](auto&&... args) {
+ return op(std::forward<decltype(arg)>(arg), args...);
+ }, argTuple);
+ };
+ return transformRange(std::forward<F>(form), captureOp);
+ }
+
+ template<class Arg0, class... Args,
+ std::enable_if_t<not isForm_v<Arg0>, int> = 0>
+ auto operator()(Arg0&& arg0, Args&&... args) const
+ {
+ return op_(std::forward<Arg0>(arg0), std::forward<Args...>(args...));
+ }
+
+ private:
+ Op op_;
+ };
+
+
+
+ namespace Impl {
+
+ struct TransposeOp {
+ template<class K, int n, int m>
+ auto operator()(const Dune::FieldMatrix<K,n,m>& x) const
+ {
+ Dune::FieldMatrix<K,m,n> y;
+ for(std::size_t i=0; i<x.N(); ++i)
+ for(std::size_t j=0; j<x.M(); ++j)
+ y[j][i] = x[i][j];
+ return y;
+ }
+ };
+
+ } // namespace Impl
+
+
+
+ template<class F>
+ using TransposedForm = TransformedForm<F, Dune::Fufem::Forms::Impl::TransposeOp>;
+
+ template<class F,
+ std::enable_if_t<isForm_v<F>, int> = 0>
+ auto transpose(const F& f)
+ {
+ return TransposedForm<F>(f);
+ }
+
+
+
+ template<class F,
+ std::enable_if_t<isForm_v<F>, int> = 0>
+ auto operator- (const F& f)
+ {
+ return transformRange(f, [](const auto& x) { return -x; });
+ }
+
+
+
+ template<class F,
+ std::enable_if_t<isForm_v<F>, int> = 0>
+ auto D(const F& f, std::size_t k)
+ {
+ return transformRange(grad(f), [k](auto&& g) { return g[k]; });
+ }
+
+
+
+} // namespace Dune::Fufem::Forms
+
+
+#endif // DUNE_FUFEM_FORMS_TRANSFORMEDFORM_HH
diff --git a/dune/fufem-forms/unaryforms.hh b/dune/fufem-forms/unaryforms.hh
new file mode 100644
index 0000000..ce2b992
--- /dev/null
+++ b/dune/fufem-forms/unaryforms.hh
@@ -0,0 +1,400 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FUFEM_FORMS_UNARYFORMS_HH
+#define DUNE_FUFEM_FORMS_UNARYFORMS_HH
+
+#include <cstddef>
+#include <type_traits>
+#include <tuple>
+#include <vector>
+
+#include <dune/common/typetraits.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/indices.hh>
+
+#include <dune/typetree/childextraction.hh>
+#include <dune/typetree/treepath.hh>
+
+#include <dune/functions/functionspacebases/subspacebasis.hh>
+
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
+
+#include <dune/fufem-forms/baseclass.hh>
+#include <dune/fufem-forms/transformedform.hh>
+
+
+namespace Dune::Fufem::Forms {
+
+
+
+ template<std::size_t k>
+ struct UnaryForm : public Form
+ {
+ static constexpr std::size_t argIndex = k;
+ static constexpr std::size_t arity = 1;
+
+ template<class Element>
+ void bindToElement(const Element&) {}
+ };
+
+
+
+ // Forward declarations
+
+ template<class B, class TP, std::size_t argIndex>
+ class FEFunctionJacobianForm;
+
+ template<class B, class TP, std::size_t argIndex>
+ class FEFunctionDivergenceForm;
+
+ template<class B, class TP, std::size_t argIndex>
+ class FEFunctionForm;
+
+
+
+ // Linear map on FE-space that maps v onto itself
+ template<class B, class TP, std::size_t argIndex>
+ class FEFunctionForm : public UnaryForm<argIndex>
+ {
+ public:
+ using Basis = B;
+ using TreePath = TP;
+ using LocalView = typename Basis::LocalView;
+ using Tree = typename LocalView::Tree;
+ using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
+
+ private:
+ using LeafNode = typename Dune::TypeTree::ChildForTreePath<Tree, Impl::LeafTreePath<Tree,TreePath>>;
+
+ using Buffer = std::vector<typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::RangeType>;
+ using ValueCache = typename ShapeFunctionCache<Tree>::template NodalValueCache<LeafNode>;
+
+ public:
+
+ using Element = typename Basis::LocalView::Element;
+ using Coordinate = typename Element::Geometry::LocalCoordinate;
+
+ FEFunctionForm(const Basis& basis, const TreePath& treePath) :
+ basis_(basis),
+ treePath_(treePath),
+ leafTreePath_(Impl::leafTreePath<Tree>(treePath)),
+ leafNode_(nullptr)
+ {}
+
+ auto quadratureRuleKey() const
+ {
+ return QuadratureRuleKey(leafNode_->finiteElement());
+ }
+
+ void bindToLocalView(const LocalView& localView)
+ {
+ leafNode_ = & Dune::TypeTree::child(localView.tree(), leafTreePath_);
+ }
+
+ template<class Cache>
+ void bindToCache(Cache& cache)
+ {
+ valueCache_ = &(cache.getValues(leafTreePath_));
+ }
+
+ template<int k>
+ const auto& treePath() const
+ {
+ return treePath_;
+ }
+
+ template<int k>
+ const auto& basis() const
+ {
+ return basis_;
+ }
+
+ auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
+ {
+ const auto& values = (*valueCache_)[quadPointIndex];
+ if constexpr(Node::isLeaf)
+ return [&](const auto& args) {
+ return values[std::get<argIndex>(args)];
+ };
+ else if constexpr(Node::isPower and Node::ChildType::isLeaf)
+ return [&](const auto& args) {
+ using Field = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
+ constexpr auto components = Node::degree();
+ auto componentIndex = std::get<argIndex>(args) / leafNode_->size();
+ auto shapeFunctionIndex = std::get<argIndex>(args) % leafNode_->size();
+ auto result = Dune::FieldVector<Field,components>(0);
+ result[componentIndex] = values[shapeFunctionIndex];
+ return result;
+ };
+ }
+
+ friend auto jacobian(const FEFunctionForm& f)
+ {
+ return FEFunctionJacobianForm<Basis, TreePath, argIndex>(f.basis_, f.treePath_);
+ }
+
+ friend auto gradient(const FEFunctionForm& f)
+ {
+ return Fufem::Forms::transpose(jacobian(f));
+ }
+
+ friend auto grad(const FEFunctionForm& f)
+ {
+ return Fufem::Forms::transpose(jacobian(f));
+ }
+
+ friend auto divergence(const FEFunctionForm& f)
+ {
+ return FEFunctionDivergenceForm<Basis, TreePath, argIndex>(f.basis_, f.treePath_);
+ }
+
+ friend auto div(const FEFunctionForm& f)
+ {
+ return divergence(f);
+ }
+
+ template<std::size_t i>
+ auto childForm(Dune::index_constant<i> childIndex) const
+ {
+ auto childTP = Dune::TypeTree::push_back(treePath_, childIndex);
+ return FEFunctionForm<Basis, decltype(childTP), argIndex>(basis_, childTP);
+ }
+
+ private:
+ const Basis& basis_;
+ const TreePath treePath_;
+ Impl::LeafTreePath<Tree,TreePath> leafTreePath_;
+ const LeafNode* leafNode_;
+ mutable Buffer evaluationBuffer_;
+ const ValueCache* valueCache_;
+ };
+
+
+
+ // Linear map on FE-space that maps v onto its Jacobian
+ template<class B, class TP, std::size_t argIndex>
+ class FEFunctionJacobianForm : public UnaryForm<argIndex>
+ {
+ public:
+ using Basis = B;
+ using TreePath = TP;
+ using LocalView = typename Basis::LocalView;
+ using Tree = typename LocalView::Tree;
+ using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
+
+ private:
+ using LeafNode = typename Dune::TypeTree::ChildForTreePath<Tree, Impl::LeafTreePath<Tree,TreePath>>;
+
+ using Buffer = std::vector<typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::JacobianType>;
+ using GlobalJacobianCache = typename ShapeFunctionCache<Tree>::template NodalGlobalJacobianCache<LeafNode>;
+
+ public:
+
+ using Element = typename Basis::LocalView::Element;
+ using Coordinate = typename Element::Geometry::LocalCoordinate;
+
+ FEFunctionJacobianForm(const Basis& basis, const TreePath& treePath) :
+ basis_(basis),
+ treePath_(treePath),
+ leafTreePath_(Impl::leafTreePath<Tree>(treePath)),
+ leafNode_(nullptr)
+ {}
+
+ auto quadratureRuleKey() const
+ {
+ return QuadratureRuleKey(leafNode_->finiteElement()).derivative();
+ }
+
+ void bindToLocalView(const LocalView& localView)
+ {
+ leafNode_ = & Dune::TypeTree::child(localView.tree(), leafTreePath_);
+ }
+
+ template<class Cache>
+ void bindToCache(Cache& cache)
+ {
+ globalJacobianCache_ = &(cache.getGlobalJacobians(leafTreePath_));
+ }
+
+ template<int k>
+ const auto& treePath() const
+ {
+ return treePath_;
+ }
+
+ template<int k>
+ const auto& basis() const
+ {
+ return basis_;
+ }
+
+ auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
+ {
+ const auto& globalJacobians = (*globalJacobianCache_)[quadPointIndex];
+ if constexpr(Node::isLeaf)
+ return [&](const auto& args) {
+ return globalJacobians[std::get<argIndex>(args)];
+ };
+ else if constexpr(Node::isPower and Node::ChildType::isLeaf)
+ return [&](const auto& args) {
+ using F = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
+ 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);
+ result[componentIndex] = globalJacobians[shapeFunctionIndex][0];
+ return result;
+ };
+ }
+
+ template<std::size_t i>
+ auto childForm(Dune::index_constant<i> childIndex) const
+ {
+ auto childTP = Dune::TypeTree::push_back(treePath_, childIndex);
+ return FEFunctionJacobianForm<Basis, decltype(childTP), argIndex>(basis_, childTP);
+ }
+
+ private:
+ const Basis& basis_;
+ const TreePath treePath_;
+ Impl::LeafTreePath<Tree,TreePath> leafTreePath_;
+ const LeafNode* leafNode_;
+ mutable Buffer evaluationBuffer_;
+ const GlobalJacobianCache* globalJacobianCache_;
+ };
+
+
+
+ // Linear map on FE-space that maps v onto its divergence
+ template<class B, class TP, std::size_t argIndex>
+ class FEFunctionDivergenceForm : public UnaryForm<argIndex>
+ {
+ public:
+ using Basis = B;
+ using TreePath = TP;
+ using LocalView = typename Basis::LocalView;
+ using Tree = typename LocalView::Tree;
+ using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
+
+ private:
+ using LeafNode = typename Dune::TypeTree::ChildForTreePath<Tree, Impl::LeafTreePath<Tree,TreePath>>;
+
+ using Buffer = std::vector<typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::JacobianType>;
+ using GlobalJacobianCache = typename ShapeFunctionCache<Tree>::template NodalGlobalJacobianCache<LeafNode>;
+
+ public:
+
+ using Element = typename Basis::LocalView::Element;
+ using Coordinate = typename Element::Geometry::LocalCoordinate;
+
+ FEFunctionDivergenceForm(const Basis& basis, const TreePath& treePath) :
+ basis_(basis),
+ treePath_(treePath),
+ leafTreePath_(Impl::leafTreePath<Tree>(treePath)),
+ leafNode_(nullptr)
+ {}
+
+ auto quadratureRuleKey() const
+ {
+ return QuadratureRuleKey(leafNode_->finiteElement()).derivative();
+ }
+
+ void bindToLocalView(const LocalView& localView)
+ {
+ leafNode_ = & Dune::TypeTree::child(localView.tree(), leafTreePath_);
+ }
+
+ template<class Cache>
+ void bindToCache(Cache& cache)
+ {
+ globalJacobianCache_ = &(cache.getGlobalJacobians(leafTreePath_));
+ }
+
+ template<int k>
+ const auto& treePath() const
+ {
+ return treePath_;
+ }
+
+ template<int k>
+ const auto& basis() const
+ {
+ return basis_;
+ }
+
+ auto operator()(const Coordinate& x, std::size_t quadPointIndex) const
+ {
+ const auto& globalJacobians = (*globalJacobianCache_)[quadPointIndex];
+ if constexpr(Node::isLeaf)
+ return [&](const auto& args) {
+ using Field = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
+ constexpr auto dimension = LeafNode::FiniteElement::Traits::LocalBasisType::Traits::dimDomain;
+ auto div = Field{0};
+ for(std::size_t i=0; i<dimension; ++i)
+ div += globalJacobians[std::get<argIndex>(args)][i][i];
+ return div;
+ };
+ else if constexpr(Node::isPower and Node::ChildType::isLeaf)
+ return [&](const auto& args) {
+ // Here we should compute the trace of the global Jacobian J
+ // that we would need to build first. But since only the
+ // componentIndex'th row of J is nonzero, we can simply
+ // return the diagonal entry of this row which coincides
+ // with the componentIndex'th entry of the gradient of
+ // the respective scalar basis function.
+ auto componentIndex = std::get<argIndex>(args) / leafNode_->size();
+ auto shapeFunctionIndex = std::get<argIndex>(args) % leafNode_->size();
+ return globalJacobians[shapeFunctionIndex][0][componentIndex];
+ };
+ }
+
+ private:
+ const Basis& basis_;
+ const TreePath treePath_;
+ Impl::LeafTreePath<Tree,TreePath> leafTreePath_;
+ const LeafNode* leafNode_;
+ mutable Buffer evaluationBuffer_;
+ const GlobalJacobianCache* globalJacobianCache_;
+ };
+
+
+
+
+ template<class Basis>
+ auto testFunction(const Basis& basis)
+ {
+ using RootBasis = std::decay_t<decltype(basis.rootBasis())>;
+ return FEFunctionForm<RootBasis, typename Basis::PrefixPath, 0>(basis.rootBasis(), basis.prefixPath());
+ }
+
+ template<class Basis, class... Args, std::enable_if_t<(sizeof...(Args)>0), int> = 0>
+ auto testFunction(const Basis& basis, Args&&... args)
+ {
+ return testFunction(Dune::Functions::subspaceBasis(basis, args...));
+ }
+
+
+
+
+ template<class Basis>
+ auto trialFunction(const Basis& basis)
+ {
+ using RootBasis = std::decay_t<decltype(basis.rootBasis())>;
+ return FEFunctionForm<RootBasis, typename Basis::PrefixPath, 1>(basis.rootBasis(), basis.prefixPath());
+ }
+
+ template<class Basis, class... Args, std::enable_if_t<(sizeof...(Args)>0), int> = 0>
+ auto trialFunction(const Basis& basis, Args&&... args)
+ {
+ return trialFunction(Dune::Functions::subspaceBasis(basis, args...));
+ }
+
+
+
+} // namespace Dune::Fufem::Forms
+
+
+
+#endif // DUNE_FUFEM_FORMS_UNARYFORMS_HH
diff --git a/src/fracture.cc b/src/fracture.cc
index ebcf7f7..2f9417e 100644
--- a/src/fracture.cc
+++ b/src/fracture.cc
@@ -167,7 +167,7 @@ int main (int argc, char *argv[]) try
auto displacementBasis = Dune::Functions::subspaceBasis(basis, _0);
auto damageBasis = Dune::Functions::subspaceBasis(basis, _1);
- auto trace = FormOperator([](const auto& M) {
+ auto trace = RangeOperator([](const auto& M) {
double tr = 0;
for (int ii=0; ii<dim; ii++)
tr += M[ii][ii];
--
GitLab