localoperators.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_LOCALOPERATORS_HH
#define DUNE_FUFEM_FORMS_LOCALOPERATORS_HH
 
#include <array>
#include <cstddef>
#include <tuple>
#include <type_traits>
#include <utility>
 
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/indices.hh>
#include <dune/common/promotiontraits.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/tuplevector.hh>
#include <dune/common/typetraits.hh>
 
#include <dune/istl/scaledidmatrix.hh>
#include <dune/istl/multitypeblockvector.hh>
 
 
 
namespace Dune::Fufem::Forms::LocalOperators {
 
 
  struct AssignOp {
 
    template<class X, class Y>
    requires requires(X& x, const Y& y) { x = y; }
    void operator()(X& x, const Y& y) const
    {
      x = y;
    }
 
    template<class... Xi, class Y>
    requires (Dune::IsNumber<Y>::value)
    void operator()(std::tuple<Xi...>& x, const Y& y) const
    {
      return std::apply([&](auto&... xi) {
        ((xi = y),...);
      }, x);
    }
 
  };
 
  struct SumOp {
 
    template<class K1, class K2>
    requires requires(K1 x, const K2 y) { x += y; }
    static void addTo(K1& x, const K2& y)
    {
      x += y;
    }
 
    template<class K1, class K2, int n>
    static void addTo(Dune::FieldVector<K1,n>& x, const Dune::FieldVector<K2,n>& y)
    {
      x += y;
    }
 
    template<class K1, class K2, int n, int m>
    static void addTo(Dune::FieldMatrix<K1,n,m>& x, const Dune::FieldMatrix<K2,n,m>& y)
    {
      x += y;
    }
 
    template<class K1, class K2, int n>
    static void addTo(Dune::ScaledIdentityMatrix<K1,n>& x, const Dune::ScaledIdentityMatrix<K2,n>& y)
    {
      x.scalar() += y.scalar();
    }
 
    template<class K1, class K2, int n>
    static void addTo(Dune::FieldMatrix<K1,n,n>& x, const Dune::ScaledIdentityMatrix<K2,n>& y)
    {
      for(std::size_t i=0; i<n; ++i)
        x[i][i] += y.scalar();
    }
 
    template<class K1, class K2, int n>
    static void addTo(Dune::FieldMatrix<K1,n,1>& x, const Dune::FieldVector<K2,n>& y)
    {
      for(std::size_t i=0; i<n; ++i)
        x[i][0] += y[i];
    }
 
    template<class K1, class K2, int n>
    static void addTo(Dune::FieldMatrix<K1,1,n>& x, const Dune::FieldVector<K2,n>& y)
    {
      for(std::size_t i=0; i<n; ++i)
        x[0][i] += y[i];
    }
 
    template<class K1, class K2, int n>
    static void addTo(Dune::FieldVector<K1,n>& x, const Dune::FieldMatrix<K2,n,1>& y)
    {
      for(std::size_t i=0; i<n; ++i)
        x[i] += y[i][0];
    }
 
    template<class K1, class K2, int n>
    static void addTo(Dune::FieldVector<K1,n>& x, const Dune::FieldMatrix<K2,1,n>& y)
    {
      for(std::size_t i=0; i<n; ++i)
        x[i] += y[0][i];
    }
 
    template<class... Xi, class... Yi>
    requires(sizeof...(Xi) == sizeof...(Yi))
    static void addTo(std::tuple<Xi...>& x, const std::tuple<Yi...>& y)
    {
      constexpr auto indices = std::index_sequence_for<Xi...>{};
      Dune::unpackIntegerSequence([&](auto... i) {
        ((std::get<i>(x) += std::get<i>(y)),...);
      }, indices);
    }
 
    template<class K1, class K2>
    requires requires(const K1& x, const K2& y) { x + y; }
    auto operator()(const K1& x, const K2& y) const
    {
      return x + y;
    }
 
    template<class K1, class K2, int n>
    auto operator()(const Dune::FieldVector<K1,n>& x, const Dune::FieldVector<K2,n>& 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::FieldMatrix<K2,n,m>& y) const
    {
      return x + y;
    }
 
    template<class K1, class K2, int n>
    auto operator()(const Dune::ScaledIdentityMatrix<K1,n>& x, const Dune::ScaledIdentityMatrix<K2,n>& y) const
    {
      using K = typename Dune::PromotionTraits<K1,K2>::PromotedType;
      return Dune::ScaledIdentityMatrix<K, n>(x.scalar() + y.scalar());
    }
 
    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 Dune::PromotionTraits<K1,K2>::PromotedType;
      auto z = Dune::FieldMatrix<K,n,1>(x);
      for(std::size_t i=0; i<n; ++i)
        z[i][0] += y[i];
      return z;
    }
 
    template<class K1, class K2, int n>
    auto operator()(const Dune::FieldMatrix<K1,1,n>& x, const Dune::FieldVector<K2,n>& y) const
    {
      using K = typename Dune::PromotionTraits<K1,K2>::PromotedType;
      auto z = Dune::FieldMatrix<K,1,n>(x);
      for(std::size_t i=0; i<n; ++i)
        z[0][i] += y[i];
      return z;
    }
 
    template<class K1, class K2, int n>
    auto operator()(const Dune::FieldVector<K1,n>& x, const Dune::FieldMatrix<K2,n,1>& y) const
    {
      return (*this)(y,x);
    }
 
    template<class K1, class K2, int n>
    auto operator()(const Dune::FieldVector<K1,n>& x, const Dune::FieldMatrix<K2,1,n>& y) const
    {
      return (*this)(y,x);
    }
 
    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 Dune::PromotionTraits<K1,K2>::PromotedType;
      auto z = Dune::FieldMatrix<K,n,n>(x);
      addTo(z, y);
      return z;
    }
 
    template<class K1, class K2, int n>
    auto operator()(const Dune::ScaledIdentityMatrix<K1,n>& x, const Dune::FieldMatrix<K2,n,n>& y) const
    {
      return (*this)(y,x);
    }
 
    template<class T1>
    auto operator()(const T1& x0) const
    {
      return x0;
    }
 
    template<class T1, class T2, class T3, class... Ts>
    requires requires(SumOp sumOp, const T1 x1, const T2 x2, const T3 x3, const Ts... xs) { sumOp( sumOp(x1, x2), x3, xs...); }
    auto operator()(const T1& x1, const T2& x2, const T3& x3, const Ts&... xs) const
    {
      return (*this)( (*this)(x1, x2), x3, xs...);
    }
 
    template<class... L, class... R>
    requires(sizeof...(L) == sizeof...(R))
    auto operator()(const Dune::MultiTypeBlockVector<L...>& x, const Dune::MultiTypeBlockVector<R...>& y) const
    {
      constexpr auto indices = std::index_sequence_for<L...>{};
      using Z = Dune::MultiTypeBlockVector<std::decay_t<decltype(std::declval<L>() + std::declval<R>())>...>;
      return Dune::unpackIntegerSequence([&](auto... i) {
        return Z((*this)(x[i],y[i]) ...);
      }, indices);
    }
 
    template<class... Xi, class... Yi>
    requires(sizeof...(Xi) == sizeof...(Yi))
    auto operator()(const Dune::TupleVector<Xi...>& x, const Dune::TupleVector<Yi...>& y) const
    {
      constexpr auto indices = std::index_sequence_for<Xi...>{};
      using Z = Dune::TupleVector<std::decay_t<decltype(std::declval<Xi>() + std::declval<Yi>())>...>;
      return Dune::unpackIntegerSequence([&](auto... i) {
        return Z((*this)(x[i],y[i]) ...);
      }, indices);
    }
 
  };
 
 
 
  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 Dune::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 Dune::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, int m>
    auto operator()(const Dune::FieldMatrix<K1,n,m>& x, const Dune::FieldVector<K2,n>& y) const
    {
      using K = typename Dune::PromotionTraits<K1,K2>::PromotedType;
      auto result = Dune::FieldVector<K,m>(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::FieldVector<K1,n>& x, const Dune::FieldMatrix<K2,n,m>& y) const
    {
      using K = typename Dune::PromotionTraits<K1,K2>::PromotedType;
      auto result = Dune::FieldVector<K,m>(0);
      for(std::size_t i=0; i<n; ++i)
        result += y[i]*x[i];
      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 Dune::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 Dune::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();
    }
 
    template<class L, class R, long unsigned int n>
    auto operator()(const std::array<L, n>& x, const std::array<R, n>& y) const
    {
      auto result = DotOp()(x[0], y[0]);
      for(std::size_t i=1; i<n; ++i)
        result += DotOp()(x[i], y[i]);
      return result;
    }
 
    template<class... L, class... R>
    requires(sizeof...(L) == sizeof...(R))
    auto operator()(const std::tuple<L...>& x, const std::tuple<R...>& y) const
    {
      constexpr auto indices = std::index_sequence_for<L...>{};
      return Dune::unpackIntegerSequence([&](auto... i) {
        auto dot = DotOp();
        return (dot(std::get<i>(x), std::get<i>(y)) + ...);
      }, indices);
    }
 
  };
 
  struct MultOp {
 
    template<class X, class Y>
    requires requires(const X x, const Y y) { x*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 Dune::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*y;
    }
 
    template<class K, class Y>
    auto operator()(const Dune::FieldVector<K,1>& x, const Y& y) const
    {
      return x*y;
    }
 
    template<class X, class... Yi>
    requires(Dune::IsNumber<X>::value)
    auto operator()(const X& x, const Dune::TupleVector<Yi...>& y) const
    {
      constexpr auto indices = std::index_sequence_for<Yi...>{};
      using Z = Dune::TupleVector<std::decay_t<decltype(x * std::declval<Yi>())>...>;
      return Dune::unpackIntegerSequence([&](auto... i) {
        return Z(x*y[i] ...);
      }, indices);
    }
 
    template<class... Xi, class Y>
    requires(Dune::IsNumber<Y>::value)
    auto operator()(const Dune::TupleVector<Xi...>& x, const Y& y) const
    {
      constexpr auto indices = std::index_sequence_for<Xi...>{};
      using Z = Dune::TupleVector<std::decay_t<decltype(std::declval<Xi>() * y)>...>;
      return Dune::unpackIntegerSequence([&](auto... i) {
        return Z(x[i]*y ...);
      }, indices);
    }
 
    template<class X, class... Yi>
    requires(Dune::IsNumber<X>::value)
    auto operator()(const X& x, const Dune::MultiTypeBlockVector<Yi...>& y) const
    {
      constexpr auto indices = std::index_sequence_for<Yi...>{};
      using Z = Dune::MultiTypeBlockVector<std::decay_t<decltype(x * std::declval<Yi>())>...>;
      return Dune::unpackIntegerSequence([&](auto... i) {
        return Z(x*y[i] ...);
      }, indices);
    }
 
    template<class... Xi, class Y>
    requires(Dune::IsNumber<Y>::value)
    auto operator()(const Dune::MultiTypeBlockVector<Xi...>& x, const Y& y) const
    {
      constexpr auto indices = std::index_sequence_for<Xi...>{};
      using Z = Dune::MultiTypeBlockVector<std::decay_t<decltype(std::declval<Xi>() * y)>...>;
      return Dune::unpackIntegerSequence([&](auto... i) {
        return Z(x[i]*y ...);
      }, indices);
    }
 
  };
 
 
 
  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 CurlFromJacobianOp {
    template<class K, int n>
    auto operator()(const Dune::FieldMatrix<K,n,n>& J) const
    {
      if constexpr (n==2)
        return Dune::FieldVector<K, 1>(J[1][0] - J[0][1]);
      if constexpr (n==3)
        return Dune::FieldVector<K, 3>({J[2][1] - J[1][2], J[0][2] - J[2][0], J[1][0] - J[0][1]});
    }
  };
 
 
 
  struct SymOp {
    template<class Matrix>
    auto operator()(const Matrix& M) const
    {
      return 0.5*(M + M.transposed());
    }
  };
 
 
 
  struct TraceOp {
    template<class Matrix>
    auto operator()(const Matrix& M) const
    {
      using T = std::decay_t<decltype(M[0][0])>;
      T tr = 0;
      for (auto k : Dune::range(M.N()))
        tr += M[k][k];
      return tr;
    }
  };
 
 
 
  struct InvertOp {
    template<class T>
    requires requires(const T t) { 1./t; }
    auto operator()(const T& t) const
    {
      return 1./t;
    }
  };
 
 
 
  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<Y>(y), std::forward<X>(x));
    }
 
  private:
    Op op_;
  };
 
 
 
  template<class... Ops>
  class ComposedOp;
 
  template<class OuterOp, class InnerOp>
  class ComposedOp<OuterOp, InnerOp>
  {
  public:
    ComposedOp(const OuterOp& outerOp, const InnerOp& innerOp) :
      outerOp_(outerOp),
      innerOp_(innerOp)
    {}
 
    template<class... Args>
    decltype(auto) operator()(Args&&... args) const
    {
      return outerOp_(innerOp_(std::forward<Args>(args)...));
    }
 
  private:
    OuterOp outerOp_;
    InnerOp innerOp_;
  };
 
  template<class OuterOp, class InnerOp0, class InnerOp1, class... InnerOps>
  class ComposedOp<OuterOp, InnerOp0, InnerOp1, InnerOps...>
  {
    using ComposedInnerOps = ComposedOp<InnerOp0, InnerOp1, InnerOps...>;
  public:
    ComposedOp(const OuterOp& outerOp, const InnerOp0& innerOp0, const InnerOp1& innerOp1, const InnerOps&... innerOps) :
      outerOp_(outerOp),
      innerOps_(innerOp0, innerOp1, innerOps...)
    {}
 
    ComposedOp(const OuterOp& outerOp, const ComposedInnerOps& innerOps) :
      outerOp_(outerOp),
      innerOps_(innerOps)
    {}
    template<class... Args>
    decltype(auto) operator()(Args&&... args) const
    {
      return outerOp_(innerOps_(std::forward<Args>(args)...));
    }
 
  private:
    OuterOp outerOp_;
    ComposedInnerOps innerOps_;
  };
 
 
  template<class OuterOp, class InnerOp>
  auto localCompose(const OuterOp& outerOp, const InnerOp& innerOp)
  {
    return ComposedOp<OuterOp, InnerOp>(outerOp, innerOp);
  }
 
  template<class OuterOp, class... InnerOps>
  auto localCompose(const OuterOp& outerOp, const ComposedOp<InnerOps...>& innerOp)
  {
    return ComposedOp<OuterOp, InnerOps...>(outerOp, innerOp);
  }
 
  auto localCompose(const TraceOp& outerOp, const SymOp& innerOp)
  {
    return outerOp;
  }
 
  auto localCompose(const SymOp& outerOp, const TransposeOp& innerOp)
  {
    return outerOp;
  }
 
  auto localCompose(const TransposeOp& outerOp, const SymOp& innerOp)
  {
    return innerOp;
  }
 
 
 
 
} // namespace Dune::Fufem::Forms::LocalOperators
 
 
 
#endif // DUNE_FUFEM_FORMS_LOCALOPERATORS_HH