productoperator.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_PRODUCTOPERATOR_HH
#define DUNE_FUFEM_FORMS_PRODUCTOPERATOR_HH
 
#include <cstddef>
#include <type_traits>
#include <tuple>
#include <utility>
#include <functional>
 
#include <dune/common/typelist.hh>
 
#include <dune/fufem/forms/baseclass.hh>
#include <dune/fufem/forms/tensors.hh>
#include <dune/fufem/forms/localfunctionadaptor.hh>
 
 
 
namespace Dune::Fufem::Forms::Impl {
 
  template<class Operator0, class Operator1>
  struct ProductOperatorBaseTraits
  {
    static auto baseType() {
      if constexpr ((Operator0::arity==1) and (Operator1::arity==0))
        return Dune::Fufem::Forms::UnaryOperator<Operator0::argIndex>();
      else if constexpr ((Operator0::arity==0) and (Operator1::arity==1))
        return Dune::Fufem::Forms::UnaryOperator<Operator1::argIndex>();
      else
        return Dune::Fufem::Forms::MultilinearOperator<Operator0::arity+Operator1::arity>();
    }
    using type = std::decay_t<decltype(baseType())>;
  };
 
 
} // namespace Dune::Fufem::Forms::Impl
 
 
 
namespace Dune::Fufem::Forms {
 
 
 
  template<class P, class Operator0, class Operator1>
  class ProductOperator : public Impl::ProductOperatorBaseTraits<Operator0, Operator1>::type
  {
    using Contraction = P;
    using LocalOperator0 = decltype(localOperator(std::declval<Operator0>()));
    using LocalOperator1 = decltype(localOperator(std::declval<Operator1>()));
 
  public:
 
    using Element = typename Operator0::Element;
    using Range = decltype(std::declval<Contraction>()(std::declval<typename Operator0::Range>(), std::declval<typename Operator1::Range>()));
 
    ProductOperator(const Contraction& product, const Operator0& operator0, const Operator1& operator1) :
      product_(product),
      operator0_(operator0),
      operator1_(operator1)
    {
      static_assert(isOperator_v<Operator0>, "ProductOperator: First factor must be an operator.");
      static_assert(isOperator_v<Operator1>, "ProductOperator: Second factor must be an operator.");
      static_assert(Operator0::arity+Operator1::arity <= 2, "ProductOperator: Arity sum of factors must not exceed 2");
      if constexpr((Operator0::arity==1) and (Operator1::arity==1))
        static_assert((Operator0::argIndex==0) and (Operator1::argIndex==1), "ProductOperator: Factors must be ordered according to argument index.");
    }
 
    class LocalOperator
    {
      using Op1Cache = std::vector<std::vector<typename Operator1::Range>>;
      using Op1Tensor = decltype(std::declval<LocalOperator0>()(0));
      static constexpr bool op1IsLazy = Impl::Tensor::IsLazyTensor<Op1Tensor>::value;
    public:
 
      using Element = typename LocalOperator0::Element;
      using Intersection = typename LocalOperator0::Intersection;
      using CacheManager = typename LocalOperator0::CacheManager;
      using Range = typename ProductOperator::Range;
 
      LocalOperator(const Contraction& product, LocalOperator0&& localOperator0, LocalOperator1&& localOperator1):
        product_(product),
        localOperator0_(std::move(localOperator0)),
        localOperator1_(std::move(localOperator1))
      {}
 
      auto quadratureRuleKey() const
      {
        return quadratureRuleKey_;
      }
 
      void bind(const Element& element)
      {
        localOperator0_.bind(element);
        localOperator1_.bind(element);
        quadratureRuleKey_ = localOperator0_.quadratureRuleKey().product(localOperator1_.quadratureRuleKey());
      }
 
      void bind(const Intersection& intersection, const Element& element, const Element& otherElement)
      {
        localOperator0_.bind(intersection, element, otherElement);
        localOperator1_.bind(intersection, element, otherElement);
        quadratureRuleKey_ = localOperator0_.quadratureRuleKey().product(localOperator1_.quadratureRuleKey());
      }
 
      void unbind()
      {
        localOperator0_.unbind();
        localOperator1_.unbind();
      }
 
      template<class... LV>
      void registerLocalViews(const LV&... lvs)
      {
        localOperator0_.registerLocalViews(lvs...);
        localOperator1_.registerLocalViews(lvs...);
      }
 
      template<class... LV>
      void registerOutsideLocalViews(const LV&... lvs)
      {
        localOperator0_.registerOutsideLocalViews(lvs...);
        localOperator1_.registerOutsideLocalViews(lvs...);
      }
 
      void registerCaches(CacheManager& cacheManager)
      {
        localOperator0_.registerCaches(cacheManager);
        localOperator1_.registerCaches(cacheManager);
      }
 
      void registerOutsideCaches(CacheManager& cacheManager)
      {
        localOperator0_.registerOutsideCaches(cacheManager);
        localOperator1_.registerOutsideCaches(cacheManager);
      }
 
      template<class... OutsideCacheManager>
      void bindToCaches(CacheManager& cacheManager, OutsideCacheManager&... outsideCacheManager)
      {
        localOperator0_.bindToCaches(cacheManager, outsideCacheManager...);
        localOperator1_.bindToCaches(cacheManager, outsideCacheManager...);
        // For a binary product cache evaluation of the second factor
        if constexpr ((Operator0::arity==1) and (Operator1::arity==1) and (op1IsLazy))
        {
          const auto& quadRule = cacheManager.rule();
          if (op1Cache_.size()<quadRule.size())
            op1Cache_.resize(quadRule.size());
          for (auto k : Dune::range(quadRule.size()))
          {
            auto y = localOperator1_(k);
            op1Cache_[k].resize(y.nnz());
            std::size_t l = 0;
            sparseForEach(y, [&](const auto& y_j, auto... j) {
              op1Cache_[k][l++] = y_j;
            });
          }
        }
      }
 
      auto operator()(std::size_t index) const
      {
        using namespace Impl::Tensor;
        using namespace Dune::Indices;
        auto x = localOperator0_(index);
        auto y = localOperator1_(index);
        if constexpr ((Operator0::arity==1) and (Operator1::arity==1) and (op1IsLazy))
        {
          auto y_cached = SparseTensor(_2, [&, y, index](auto&& f) {
              std::size_t l = 0;
              sparseForEach(y, [&](const auto&, auto... j) {
                f(op1Cache_[index][l++], j...);
              });
            },
            op1Cache_.size(),
            Dune::MetaType<typename LocalOperator1::Range>{}
          );
          return interleavedOuterProduct(std::cref(product_), x, y_cached);
        }
        else
          return interleavedOuterProduct(std::cref(product_), x, y);
      }
 
    private:
      Contraction product_;
      LocalOperator0 localOperator0_;
      LocalOperator1 localOperator1_;
      QuadratureRuleKey quadratureRuleKey_;
      Op1Cache op1Cache_;
    };
 
    friend LocalOperator localOperator(const ProductOperator& productOperator)
    {
      return LocalOperator(productOperator.contraction(), localOperator(productOperator.operator0()), localOperator(productOperator.operator1()));
    }
 
    template<bool dummy=true, std::enable_if_t<dummy and (ProductOperator::arity==0), int> = 0>
    friend LocalFunctionAdaptor<LocalOperator> localFunction(const ProductOperator& productOperator)
    {
      return LocalFunctionAdaptor<LocalOperator>(localOperator(productOperator));
    }
 
    const Operator0& operator0() const
    {
      return operator0_;
    }
 
    const Operator1& operator1() const
    {
      return operator1_;
    }
 
    const Contraction& contraction() const
    {
      return product_;
    }
 
    auto basis() const
    {
      return std::tuple_cat(operator0_.basis(), operator1_.basis());
    }
 
    auto treePath() const
    {
      return std::tuple_cat(operator0_.treePath(), operator1_.treePath());
    }
 
  private:
    Contraction product_;
    Operator0 operator0_;
    Operator1 operator1_;
  };
 
 
 
} // namespace Dune::Fufem::Forms
 
#endif // DUNE_FUFEM_FORMS_PRODUCTOPERATOR_HH