userfunctions.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_FORMS_USERFUNCTIONS_HH
#define DUNE_FUFEM_FORMS_USERFUNCTIONS_HH
 
#include <cstddef>
#include <type_traits>
#include <utility>
 
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>
 
#include <dune/istl/scaledidmatrix.hh>
 
#if DUNE_VERSION_GTE(DUNE_FUNCTIONS, 2, 11)
#include <dune/common/typetree/childaccess.hh>
#else
#include <dune/typetree/childextraction.hh>
#endif
 
#include <dune/fufem/forms/baseclass.hh>
#include <dune/fufem/forms/unaryoperators.hh>
#include <dune/fufem/forms/transformedoperator.hh>
#include <dune/fufem/forms/productoperator.hh>
#include <dune/fufem/forms/sumoperator.hh>
#include <dune/fufem/forms/boundunaryoperator.hh>
#include <dune/fufem/forms/outsideoperator.hh>
#include <dune/fufem/forms/localoperators.hh>
#include <dune/fufem/forms/geometryoperators.hh>
#include <dune/fufem/typetreefallback.hh>
 
 
 
 
namespace Dune::Fufem::Forms::Impl {
 
 
 
  // ************************************************************
  // Implementation of compose(outerOp, innerOp)
  // ************************************************************
 
  template<class OuterTransform, class InnerTransform, class InnerOp>
  requires isOperator_v<InnerOp>
  auto composeImpl(const OuterTransform& outerTransform, const TransformedOperator<InnerTransform, InnerOp>& op)
  {
    return TransformedOperator(LocalOperators::localCompose(outerTransform, op.transformation()), op.baseOperator());
  }
 
  template<class OuterOp, class InnerOp>
  requires isOperator_v<InnerOp>
  auto composeImpl(const OuterOp& outerOp, const InnerOp& op)
  {
    return TransformedOperator(outerOp, op);
  }
 
  template<class OuterOp, class... InnerOps>
  auto composeImpl(const OuterOp& outerOp, const Dune::Fufem::Forms::SumOperator<InnerOps...>& innerOps)
  {
    if constexpr(SumOperator<InnerOps...>::arity==0)
      return TransformedOperator(outerOp, innerOps);
    else
      return std::apply([&](const auto&...op_i) {
          return Dune::Fufem::Forms::SumOperator(composeImpl(outerOp, op_i) ...);
      }, innerOps.operators());
  }
 
 
 
  // ************************************************************
  // Implementation of sum(l,r).
  // ************************************************************
 
  // This has four special cases to collapse nested sums
  // and two for the case that one argument is a constant.
 
  template<class L, class R>
  requires (isOperator_v<L> and isOperator_v<R>)
  auto sum(const L& l, const R& r)
  {
    return SumOperator<L, R>(l, r);
  }
 
  template<class L, class... Ri>
  requires isOperator_v<L>
  auto sum(const L& l, const SumOperator<Ri...>& r)
  {
    return std::apply([&](const auto&... ri) {
        return SumOperator<L, Ri...>(l, ri...);
    }, r.operators());
  }
 
  template<class... Li, class R>
  requires isOperator_v<R>
  auto sum(const SumOperator<Li...>& l, const R& r)
  {
    return std::apply([&](const auto&... li) {
        return SumOperator<Li..., R>(li..., r);
    }, l.operators());
  }
 
  template<class... Li, class... Rj>
  auto sum(const SumOperator<Li...>& l, const SumOperator<Rj...>& r)
  {
    return std::apply([&](const auto&...li) {
      return std::apply([&](const auto&...rj) {
          return SumOperator<Li..., Rj...>(li..., rj...);
      }, r.operators());
    }, l.operators());
  }
 
  // Overload for for the case that R is a constant
  template<class L, class R>
  requires (isOperatorOrSumOperator_v<L> and not isOperatorOrSumOperator_v<R>)
  auto sum(const L& l, const R& r)
  {
    return Impl::composeImpl([r](auto&& l_value) {
      return Dune::Fufem::Forms::LocalOperators::SumOp()(l_value, r);
    }, l);
  }
 
  // Overload for for the case that L is a constant
  template<class L, class R>
  requires (not isOperatorOrSumOperator_v<L> and isOperatorOrSumOperator_v<R>)
  auto sum(const L& l, const R& r)
  {
    return Impl::composeImpl([l](auto&& r_value) {
      return Dune::Fufem::Forms::LocalOperators::SumOp()(l, r_value);
    }, r);
  }
 
 
 
  // ************************************************************
  // Implementation of product(p, left, right)
  // ************************************************************
 
  // The following are the generic functions for producing
  // product operators from a binary operator and two other
  // terms. At least one of those terms should be a operator.
 
  // Overload for two plain operators
  template<class Op, class L, class R>
  requires (isOperator_v<L> and isOperator_v<R>
    and Dune::IsCallable<Op(const typename L::Range, const typename R::Range)>::value)
  auto productImpl(const Op& op, const L& l, const R& r)
  {
    static_assert(L::arity+R::arity<=2, "Trying to construct multi-linear operator with arity>2");
 
    if constexpr((L::arity==1) and (L::arity==1))
    {
      if constexpr(L::argIndex==1)
        return ProductOperator(LocalOperators::TransposedBinaryOp(op), r, l);
      else
        return ProductOperator(op, l, r);
    }
    else
      return ProductOperator(op, l, r);
  }
 
  // Overload for one raw value and one operator
  // This will wrap the raw value into a operator.
  template<class Op, class L, class R>
  requires (not isOperatorOrSumOperator_v<L> and isOperator_v<R>
    and Dune::IsCallable<Op(const L, const typename R::Range)>::value)
  auto productImpl(const Op& op, const L& l, const R& r)
  {
    using namespace std::placeholders;
    return composeImpl(std::bind(op, l, _1), r);
  }
 
  // Overload for one operator and one raw value
  // This will wrap the raw value into a operator.
  template<class Op, class L, class R>
  requires (isOperator_v<L> and not isOperatorOrSumOperator_v<R>
    and Dune::IsCallable<Op(const typename L::Range, const R)>::value)
  auto productImpl(const Op& op, const L& l, const R& r)
  {
    using namespace std::placeholders;
    return composeImpl(std::bind(op, _1, r), l);
  }
 
  // Overload for one sum-operator and one plain term
  // This will multiply out the product for the sum-operator.
  template<class Op, class... Li, class R>
  requires requires (const Op op, const Li... li, const R r) { (productImpl(op, li, r),...); }
  auto productImpl(const Op& op, const SumOperator<Li...>& l, const R& r)
  {
    return std::apply([&](const auto&... li) {
        return (productImpl(op, li, r) + ...);
    }, l.operators());
  }
 
  // Overload for one plain term and one sum-operator
  // This will multiply out the product for the sum-operator.
  template<class Op, class L, class... Ri>
  requires requires (const Op op, const L l, const Ri... ri) { (productImpl(op, l, ri),...); }
  auto productImpl(const Op& op, const L& l, const SumOperator<Ri...>& r)
  {
    return std::apply([&](const auto&...ri) {
        return (productImpl(op, l, ri) + ...);
    }, r.operators());
  }
 
  // Overload for two sum-operators
  // This will multiply out the product for both sum-operators.
  template<class Op, class... Li, class... Rj>
  requires requires (const Op op, const Li... li, const SumOperator<Rj...> r) { (productImpl(op, li, r),...); }
  auto productImpl(const Op& op, const SumOperator<Li...>& l, const SumOperator<Rj...>& r)
  {
    return std::apply([&](const auto&...li) {
      return (productImpl(op, li, r) + ...);
    }, l.operators());
  }
 
 
 
  // ************************************************************
  // Implementation of dot(l, r).
  // ************************************************************
 
  // These specializations only exist for optimization of special cases.
 
  // This is the default implementation.
  template<class L, class R>
  requires ((isOperatorOrSumOperator_v<L> or isOperatorOrSumOperator_v<R>)
    and requires (LocalOperators::DotOp dot, const L l, const R r) { productImpl(dot, l, r); })
  auto dotImpl(const L& l, const R& r)
  {
    return productImpl(Dune::Fufem::Forms::LocalOperators::DotOp(), l, r);
  }
 
  template<class B, class TP, std::size_t argIndex_L, std::size_t argIndex_R>
  auto dotImpl(const FEFunctionOperator<B,TP,argIndex_L>& l, const FEFunctionOperator<B,TP,argIndex_R>& r);
 
  template<class B, class TP, std::size_t argIndex_L, std::size_t argIndex_R>
  auto dotImpl(const FEFunctionJacobianOperator<B,TP,argIndex_L>& l, const FEFunctionJacobianOperator<B,TP,argIndex_R>& r);
 
  template<class InnerOp>
  using TransposedOperator = TransformedOperator<Dune::Fufem::Forms::LocalOperators::TransposeOp, InnerOp>;
 
  template<class L, class R>
  auto dotImpl(const Impl::TransposedOperator<L>& l, const Impl::TransposedOperator<R>& r)
  {
    using namespace Dune::Indices;
    return Dune::Fufem::Forms::Impl::dotImpl(l.baseOperator(), r.baseOperator());
  }
 
  template<class B, class TP, std::size_t argIndex_L, std::size_t argIndex_R>
  auto dotImpl(const FEFunctionOperator<B,TP,argIndex_L>& l, const FEFunctionOperator<B,TP,argIndex_R>& r)
  {
    using Node = typename Dune::TypeTree::ChildForTreePath<typename B::LocalView::Tree, TP>;
    if constexpr (Dune::Fufem::Impl::Concept::UniformInnerTreeNode<Node>)
    {
      return unpackIntegerSequence(
        [&](auto... i) {
          return (Dune::Fufem::Forms::Impl::dotImpl(l.childOperator(i), r.childOperator(i)) + ...);
        },
        std::make_index_sequence<Node::degree()>{});
    }
    else
      return productImpl(Dune::Fufem::Forms::LocalOperators::DotOp(), l, r);
  }
 
  template<class B, class TP, std::size_t argIndex_L, std::size_t argIndex_R>
  auto dotImpl(const FEFunctionJacobianOperator<B,TP,argIndex_L>& l, const FEFunctionJacobianOperator<B,TP,argIndex_R>& r)
  {
    using Node = typename Dune::TypeTree::ChildForTreePath<typename B::LocalView::Tree, TP>;
    if constexpr (Dune::Fufem::Impl::Concept::UniformInnerTreeNode<Node>)
    {
      return unpackIntegerSequence(
        [&](auto... i) {
          return (Dune::Fufem::Forms::Impl::dotImpl(l.childOperator(i), r.childOperator(i)) + ...);
        },
        std::make_index_sequence<Node::degree()>{});
    }
    else
      return productImpl(Dune::Fufem::Forms::LocalOperators::DotOp(), l, r);
  }
 
 
 
} // namespace Dune::Fufem::Forms::Impl
 
 
 
namespace Dune::Fufem::Forms {
 
 
 
  // ************************************************************
  // Linear combination of operators of the same arity
  // ************************************************************
 
#ifdef DOXYGEN
  template<class L, class R>
  auto operator+ (const L& l, const R& r);
#else
  // Overload for two operators of the same arity
  template<class L, class R>
  requires (isOperatorOrSumOperator_v<L> or isOperatorOrSumOperator_v<R>)
  auto operator+ (const L& l, const R& r)
  {
    if constexpr (isOperatorOrSumOperator_v<L> and isOperatorOrSumOperator_v<R>)
      static_assert(L::arity==R::arity, "The operators passed to operator+ do not have the same arity");
    else if constexpr (isOperatorOrSumOperator_v<R>)
      static_assert(R::arity==0, "Plain values can only be added to 0-linear operators");
    else if constexpr (isOperatorOrSumOperator_v<L>)
      static_assert(L::arity==0, "Plain values can only be added to 0-linear operators");
    return Impl::sum(l, r);
  }
#endif
 
  template<class Op>
  requires isOperatorOrSumOperator_v<Op>
  auto operator- (const Op& op)
  {
    if constexpr (requires(typename Op::Range x) { -x; })
      return Impl::composeImpl(std::negate(), op);
    else if constexpr (requires(typename Op::Range x) { -1*x; })
      return Impl::composeImpl([](auto&& x) { return -1*x;}, op);
  }
 
#ifdef DOXYGEN
  template<class L, class R>
  auto operator- (const L& l, const R& r);
#else
  template<class L, class R>
  requires (isOperatorOrSumOperator_v<L> or isOperatorOrSumOperator_v<R>)
  auto operator- (const L& l, const R& r)
  {
    if constexpr (isOperatorOrSumOperator_v<L> and isOperatorOrSumOperator_v<R>)
      static_assert(L::arity==R::arity, "The operators passed to operator- do not have the same arity");
    else if constexpr (isOperatorOrSumOperator_v<R>)
      static_assert(R::arity==0, "Differences with plain values can only be used with 0-linear operators");
    else if constexpr (isOperatorOrSumOperator_v<L>)
      static_assert(L::arity==0, "Differences with plain values can only be used with 0-linear operators");
    return Impl::sum(l, -r);
  }
#endif
 
 
 
  // ************************************************************
  // Bilinear combination of operators
  // ************************************************************
 
  template<class Op, class L, class R>
  auto product(const Op& op, const L& l, const R& r)
  requires ((isOperatorOrSumOperator_v<L> or isOperatorOrSumOperator_v<R>)
    and requires (const Op op, const L l, const R r) { Impl::productImpl(op, l, r); })
  {
    return Impl::productImpl(op, l, r);
  }
 
  template<class L, class R, class Dummy=void>
  requires ((isOperatorOrSumOperator_v<L> or isOperatorOrSumOperator_v<R>)
    and requires (LocalOperators::MultOp mult, const L l, const R r) { Impl::productImpl(mult, l, r); })
  auto operator* (const L& l, const R& r)
  {
    // This overload is recursively constrained to cases where `operator*` makes sense.
    // Otherwise Dune::dot(a,b) may be activated for operators because it relies on availability
    // of a*b. This is a problem, because it leads to ambiguity with Forms::dot(a,b).
    return Impl::productImpl(LocalOperators::MultOp{}, l, r);
  }
 
  template<class L, class R>
  requires ((isOperatorOrSumOperator_v<L> or isOperatorOrSumOperator_v<R>)
    and requires (const L l, const R r) { Impl::dotImpl(l, r); })
  auto dot(const L& l, const R& r)
  {
    return Impl::dotImpl(l, r);
  }
 
 
 
  // ************************************************************
  // Utilities for range transformations
  // ************************************************************
 
  template<class OuterOp, class InnerOp>
  requires (isOperatorOrSumOperator_v<InnerOp>
    and requires (const OuterOp outerOp, typename InnerOp::Range inner_range) { outerOp(inner_range); })
  auto compose(const OuterOp& outerOp, const InnerOp& innerOp)
  {
    return Impl::composeImpl(outerOp, innerOp);
  }
 
 
 
  template<class OuterOp>
  class RangeOperator
  {
  public:
    RangeOperator(const OuterOp& outerOp) : outerOp_(outerOp) {}
 
    template<class Op, class... Args>
    requires isOperatorOrSumOperator_v<Op>
    auto operator()(Op&& baseOperator, Args... args) const
    {
      using namespace std::placeholders;
      return compose(std::bind(outerOp_, _1, args...), std::forward<Op>(baseOperator));
    }
 
    template<class Op>
    requires isOperatorOrSumOperator_v<Op>
    auto operator()(Op&& baseOperator) const
    {
      return compose(outerOp_, std::forward<Op>(baseOperator));
    }
 
    template<class Arg0, class... Args>
    requires (not isOperatorOrSumOperator_v<Arg0>)
    auto operator()(Arg0&& arg0, Args&&... args) const
    {
      return outerOp_(std::forward<Arg0>(arg0), std::forward<Args>(args)...);
    }
 
  private:
    OuterOp outerOp_;
  };
 
 
 
  // ************************************************************
  // Predefined linear range transformations
  // ************************************************************
 
  template<class Op>
  auto transpose(const Op& op)
  {
    return RangeOperator(Dune::Fufem::Forms::LocalOperators::TransposeOp())(op);
  }
 
  template<class Op>
  auto symmetrize(const Op& op)
  {
    return RangeOperator(Dune::Fufem::Forms::LocalOperators::SymOp())(op);
  }
 
  template<class Op>
  auto trace(const Op& op)
  {
    return RangeOperator(Dune::Fufem::Forms::LocalOperators::TraceOp())(op);
  }
 
 
 
  // ************************************************************
  // Predefined non-linear range transformations
  // ************************************************************
 
  template<class Op>
  auto inv(const Op& op)
  {
    if constexpr(isOperatorOrSumOperator_v<Op>)
      static_assert(Op::arity==0, "inv(op) can only be applied to 0-linear operators");
    return RangeOperator(Dune::Fufem::Forms::LocalOperators::InvertOp())(op);
  }
 
  template<class L, class R>
  requires ((isOperatorOrSumOperator_v<L> or isOperatorOrSumOperator_v<R>)
    and requires (LocalOperators::MultOp mult, const L l, const R r) { product(mult, l, inv(r)); })
  auto operator/ (const L& l, const R& r)
  {
    // Since inv() does not exist if R is a non-nullary operator,
    // we don't need to guard this further here.
    return product(LocalOperators::MultOp{}, l, inv(r));
  }
 
 
 
  // ************************************************************
  // Convenience functions for creating operators
  // ************************************************************
 
  class NonAffineFamily {};
 
  template<class Basis>
  auto testFunction(const Basis& basis)
  {
    using RootBasis = std::decay_t<decltype(basis.rootBasis())>;
    return FEFunctionOperator<RootBasis, typename Basis::PrefixPath, 0>(basis.rootBasis(), basis.prefixPath());
  }
 
  template<class Basis>
  auto testFunction(const Basis& basis, NonAffineFamily tag)
  {
    using RootBasis = std::decay_t<decltype(basis.rootBasis())>;
    return FEFunctionOperator<RootBasis, typename Basis::PrefixPath, 0>(basis.rootBasis(), basis.prefixPath(), false);
  }
 
  template<class Basis, class... Args>
  requires (sizeof...(Args)>0)
  auto testFunction(const Basis& basis, Args&&... args)
  {
    using LastArg = typename Dune::Functions::LastType<std::decay_t<Args>...>::type;
    if constexpr (std::is_same_v<LastArg, NonAffineFamily>)
    {
      return Dune::applyPartial([&](auto&&... treePathArgs){
          return testFunction(Dune::Functions::subspaceBasis(basis, treePathArgs...), NonAffineFamily{});
        },
        std::forward_as_tuple(std::forward<Args>(args)...),
        std::make_index_sequence<sizeof...(Args)-1>{}
      );
    }
    else
      return testFunction(Dune::Functions::subspaceBasis(basis, args...));
  }
 
  template<class Basis>
  auto trialFunction(const Basis& basis)
  {
    using RootBasis = std::decay_t<decltype(basis.rootBasis())>;
    return FEFunctionOperator<RootBasis, typename Basis::PrefixPath, 1>(basis.rootBasis(), basis.prefixPath());
  }
 
  template<class Basis>
  auto trialFunction(const Basis& basis, NonAffineFamily tag)
  {
    using RootBasis = std::decay_t<decltype(basis.rootBasis())>;
    return FEFunctionOperator<RootBasis, typename Basis::PrefixPath, 1>(basis.rootBasis(), basis.prefixPath(), false);
  }
 
  template<class Basis, class... Args>
  requires (sizeof...(Args)>0)
  auto trialFunction(const Basis& basis, Args&&... args)
  {
    using LastArg = typename Dune::Functions::LastType<std::decay_t<Args>...>::type;
    if constexpr (std::is_same_v<LastArg, NonAffineFamily>)
    {
      return Dune::applyPartial([&](auto&&... treePathArgs){
          return trialFunction(Dune::Functions::subspaceBasis(basis, treePathArgs...), NonAffineFamily{});
        },
        std::forward_as_tuple(std::forward<Args>(args)...),
        std::make_index_sequence<sizeof...(Args)-1>{}
      );
    }
    else
      return trialFunction(Dune::Functions::subspaceBasis(basis, args...));
  }
 
  template<class Op>
  requires (isOperatorOrSumOperator_v<Op> and (Op::arity<=1))
  auto outside(const Op& op)
  {
    if constexpr (isOperator_v<Op>)
      return OutsideOperator(op);
    if constexpr (isSumOperator_v<Op>)
      return std::apply([&](const auto&...op_i) {
          return Dune::Fufem::Forms::SumOperator(OutsideOperator(op_i) ...);
      }, op.operators());
  }
 
  template<class Op>
  requires (isOperatorOrSumOperator_v<Op> and (Op::arity<=1))
  auto avg(const Op& op)
  {
    return 0.5*(op+outside(op));
  }
 
  template<class Op>
  requires (isOperatorOrSumOperator_v<Op> and (Op::arity<=1))
  auto jump(const Op& op)
  {
    return op-outside(op);
  }
 
  template<class GridView>
  auto meshSize(const GridView& gridView)
  {
    return MeshSizeOperator(gridView);
  }
 
  template<class GridView>
  auto meshSize()
  {
    return MeshSizeOperator<GridView>();
  }
 
  template<class GridView>
  auto faceNormal(const GridView& gridView)
  {
    return FaceNormalOperator(gridView);
  }
 
  template<class GridView>
  auto faceNormal()
  {
    return FaceNormalOperator<GridView>();
  }
 
  template<class Op, class V>
  requires isOperatorOrSumOperator_v<Op>
  auto bindToCoefficients(const Op& op, V&& vector)
  {
    using Basis = std::decay_t<decltype(std::get<0>(op.basis()))>;
    // Small helper functions to wrap vectors using istlVectorBackend
    // if they do not already satisfy the VectorBackend interface.
    auto toConstVectorBackend = [&](auto&& v) -> decltype(auto) {
      if constexpr (models<Dune::Functions::Concept::ConstVectorBackend<Basis>, decltype(v)>()) {
        return std::forward<decltype(v)>(v);
      } else {
        return Dune::Functions::istlVectorBackend(v);
      }
    };
    return BoundUnaryOperator(op, toConstVectorBackend(std::forward<V>(vector)));
  }
 
 
 
  template<class Op, class V>
  requires (isOperatorOrSumOperator_v<Op> and not isOperatorOrSumOperator_v<std::decay_t<V>>)
  auto operator|(const Op& op, V&& v)
  {
    return bindToCoefficients(op, std::forward<V>(v));
  }
 
 
 
  // ************************************************************
  // Differential operators
  // ************************************************************
 
  // Overloads for elementary operators
 
  template<class B, class TP, std::size_t argIndex>
  auto jacobian(const FEFunctionOperator<B, TP, argIndex>& op)
  {
    return FEFunctionJacobianOperator<B, TP, argIndex>(std::get<0>(op.basis()), std::get<0>(op.treePath()), op.isAffine());
  }
 
  template<class B, class TP, std::size_t argIndex>
  auto divergence(const FEFunctionOperator<B, TP, argIndex>& op)
  {
    return FEFunctionDivergenceOperator<B, TP, argIndex>(std::get<0>(op.basis()), std::get<0>(op.treePath()), op.isAffine());
  }
 
  template<class B, class TP, std::size_t argIndex>
  auto curl(const FEFunctionOperator<B, TP, argIndex>& op)
  {
    constexpr auto dim = FEFunctionOperator<B, TP, argIndex>::Element::Geometry::mydimension;
    static_assert(dim==2 or dim==3);
    return RangeOperator(Dune::Fufem::Forms::LocalOperators::CurlFromJacobianOp())(jacobian(op));
  }
 
  // Overloads for sums of operators
 
  template<class... Ops>
  requires requires(const Ops&... ops) { SumOperator(jacobian(ops)...); }
  auto jacobian(const SumOperator<Ops...>& op)
  {
    return std::apply([&](const auto&...op) {
        return SumOperator(jacobian(op) ...);
    }, op.operators());
  }
 
  template<class... Ops>
  requires requires(const Ops&... ops) { SumOperator(jacobian(ops)...); }
  auto divergence(const SumOperator<Ops...>& op)
  {
    return std::apply([&](const auto&...op) {
        return SumOperator(divergence(op) ...);
    }, op.operators());
  }
 
  template<class... Ops>
  requires requires(const Ops&... ops) { SumOperator(curl(ops)...); }
  auto curl(const SumOperator<Ops...>& op)
  {
    return std::apply([&](const auto&...op) {
        return SumOperator(curl(op) ...);
    }, op.operators());
  }
 
  // Other operators that are essetially shot cuts to the others
 
  template<class Op>
  auto gradient(const Op& op)
  requires requires {Fufem::Forms::transpose(jacobian(op));}
  {
    return Fufem::Forms::transpose(jacobian(op));
  }
 
  template<class Op>
  auto grad(const Op& op)
  requires requires {Fufem::Forms::transpose(jacobian(op));}
  {
    return Fufem::Forms::transpose(jacobian(op));
  }
 
  template<class Op>
  auto div(const Op& op)
  requires requires {divergence(op);}
  {
    return divergence(op);
  }
 
  template<class Op>
  auto rot(const Op& op)
  requires requires {curl(op);}
  {
    return curl(op);
  }
 
  template<class Op>
  auto D(const Op& op, std::size_t k)
  requires requires {grad(op);}
  {
    return compose([k](auto&& g) { return g[k]; }, grad(op));
  }
 
 
 
} // namespace Dune::Fufem::Forms
 
 
#endif // DUNE_FUFEM_FORMS_USERFUNCTIONS_HH