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 <cstddef>
#include <type_traits>
#include <utility>
#include <array>
 
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/promotiontraits.hh>
#include <dune/common/rangeutilities.hh>
 
#include <dune/istl/scaledidmatrix.hh>
 
 
 
// This namespace contains implementations of several local operators.
// Here 'local' means that these operators are applied poinwise to
// numbers, vectors, matrices, ... . The operators are implemented as non-templated
// imutable classes that provide a templated or overloaded operator().
// Some of the operations could be implemented as plain lambda functions,
// but using named callbacks makes the expression trees much more readable
// in error messages and allows to detect certain operators to provide
// shortcuts for special cases.
namespace Dune::Fufem::Forms::LocalOperators {
 
 
 
  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 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>
    auto operator()(const T1& x0, const T2& x1, const T3& x3, const Ts&... xs) const
    requires requires { (*this)( (*this)(x0, x1), x3, xs...); }
    {
      return (*this)( (*this)(x0, x1), x3, xs...);
    }
 
  };
 
 
 
  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;
    }
 
  };
 
  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 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[0][0]*y;
    }
 
    template<class K, class Y>
    auto operator()(const Dune::FieldVector<K,1>& x, const Y& y) const
    {
      return x[0]*y;
    }
 
  };
 
 
 
  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,
      class = std::void_t<decltype(1./std::declval<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