boundunaryoperator.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_BOUNDUNARYOPERATOR_HH
#define DUNE_FUFEM_FORMS_BOUNDUNARYOPERATOR_HH
 
#include <cstddef>
#include <type_traits>
#include <tuple>
#include <utility>
 
#include <dune/common/shared_ptr.hh>
#include <dune/common/typetree/childaccess.hh>
#include <dune/common/typetree/traversal.hh>
#include <dune/common/typetree/treecontainer.hh>
#include <dune/common/typetree/treepath.hh>
 
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <dune/functions/backends/concepts.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
 
#include <dune/fufem/forms/shapefunctioncache.hh>
#include <dune/fufem/forms/baseclass.hh>
#include <dune/fufem/forms/tensors.hh>
#include <dune/fufem/forms/localfunctionadaptor.hh>
#include <dune/fufem/forms/localoperators.hh>
 
 
 
namespace Dune::Fufem::Forms {
 
 
  template<class Op, class V>
  class BoundUnaryOperator : public MultilinearOperator<0>
  {
    using UnaryOperator = Op;
    using LocalUnaryOperator = decltype(localOperator(std::declval<UnaryOperator>()));
    using Tree = typename std::decay_t<decltype(std::get<0>(std::declval<UnaryOperator>().basis()))>::LocalView::Tree;
    using NodeFlags = typename TypeTree::UniformTreeContainer<char, Tree>;
 
    template<class OP, class Func>
    static void forEachOperator(OP&& op, Func&& f)
    {
      if constexpr (isOperator_v<UnaryOperator>)
        f(op);
      else if constexpr (isSumOperator_v<UnaryOperator>)
        Impl::forEachTupleEntry(op.operators(), f);
    }
 
    using MultOp = Dune::Fufem::Forms::LocalOperators::MultOp;
    using SumOp = Dune::Fufem::Forms::LocalOperators::SumOp;
    using AssignOp = Dune::Fufem::Forms::LocalOperators::AssignOp;
 
  public:
 
    // Types present in DiscreteGlobalBasisFunctionBase
    using Basis = std::decay_t<decltype(std::get<0>(std::declval<UnaryOperator>().basis()))>;
    using Vector = V;
    using Coefficient = Dune::AutonomousValue<decltype(std::declval<Vector>()[std::declval<typename Basis::MultiIndex>()])>;
    using GridView = typename Basis::GridView;
    using EntitySet = typename Dune::Functions::GridViewEntitySet<GridView, 0>;
    using Domain = typename EntitySet::GlobalCoordinate;
    using LocalDomain = typename EntitySet::LocalCoordinate;
 
    // Types required by Operator interface
    using Element = typename EntitySet::Element;
    using Range = std::decay_t<decltype(MultOp()(std::declval<Coefficient>(), std::declval<typename UnaryOperator::Range>()))>;
 
    BoundUnaryOperator(const UnaryOperator& unaryOperator, const Vector& coefficients) :
      cacheId_(*this),
      unaryOperator_(unaryOperator),
      coefficients_(Dune::stackobject_to_shared_ptr(coefficients)),
      entitySet_(std::get<0>(unaryOperator.basis()).gridView())
    {
      static_assert(isOperatorOrSumOperator_v<UnaryOperator>, "The type passed to BoundUnaryOperator is not a operator.");
      static_assert(UnaryOperator::arity==1, "The operator passed to BoundUnaryOperator is not unary.");
    }
 
    BoundUnaryOperator(const UnaryOperator& unaryOperator, Vector&& coefficients) :
      cacheId_(*this),
      unaryOperator_(unaryOperator),
      coefficients_(std::make_shared<const Vector>(std::move(coefficients))),
      entitySet_(std::get<0>(unaryOperator.basis()).gridView())
    {
      static_assert(isOperatorOrSumOperator_v<UnaryOperator>, "The type passed to BoundUnaryOperator is not a operator.");
      static_assert(UnaryOperator::arity==1, "The operator passed to BoundUnaryOperator is not unary.");
    }
 
    BoundUnaryOperator(const UnaryOperator& unaryOperator, std::shared_ptr<const Vector> coefficients) :
      cacheId_(*this),
      unaryOperator_(unaryOperator),
      coefficients_(std::move(coefficients)),
      entitySet_(std::get<0>(unaryOperator.basis()).gridView())
    {
      static_assert(isOperatorOrSumOperator_v<UnaryOperator>, "The type passed to BoundUnaryOperator is not a operator.");
      static_assert(UnaryOperator::arity==1, "The operator passed to BoundUnaryOperator is not unary.");
    }
 
 
 
    class LocalOperator
    {
    public:
 
      using Element = typename BoundUnaryOperator::Element;
      using Intersection = typename BoundUnaryOperator::LocalUnaryOperator::Intersection;
      using CacheManager = Dune::Fufem::Forms::CacheManager<typename Element::Geometry::ctype, Element::Geometry::mydimension>;
      using Range = typename BoundUnaryOperator::Range;
 
    private:
      static constexpr int dimension = EntitySet::Element::dimension;
 
      using LocalView = typename Basis::LocalView;
      using Tree = typename LocalView::Tree;
      using size_type = typename Basis::LocalView::Tree::size_type;
      using QuadratureRule = typename CacheManager::QuadratureRule;
      using QuadraturePoint = typename QuadratureRule::value_type;
      using ValueCache = SimpleCache<double, dimension, Range>;
 
      void loadLocalDOFs(const LocalView& localView, std::vector<Coefficient>& localDoFs) const
      {
        localDoFs.resize(localView.size());
        TypeTree::forEachLeafNode(localView.tree(), [&](auto&& node, auto&& treePath) {
          if (usedNodes_[treePath])
            for (auto i : Dune::range(node.size()))
            {
              auto localIndex = node.localIndex(i);
              localDoFs[localIndex] = (*coefficients_)[localView.index(localIndex)];
            }
        });
      }
 
    public:
 
      LocalOperator(LocalUnaryOperator&& localUnaryOperator, NodeFlags usedNodes, const Basis& basis, std::shared_ptr<const Vector> coefficients, UniqueCacheId cacheId) :
        basis_(basis),
        localUnaryOperator_(std::move(localUnaryOperator)),
        usedNodes_(usedNodes),
        coefficients_(std::move(coefficients)),
        internalLocalView_{basis.localView(), basis.localView()},
        cacheId_(cacheId)
      {
        localDoFs_[0].reserve(internalLocalView_[0].maxSize());
        localDoFs_[1].reserve(internalLocalView_[1].maxSize());
        localUnaryOperator_.registerLocalViews(internalLocalView_[0]);
        localUnaryOperator_.registerOutsideLocalViews(internalLocalView_[1]);
      }
 
      LocalOperator(const LocalOperator& other) :
        LocalOperator(LocalUnaryOperator(other.localUnaryOperator_), other.usedNodes_, other.basis_, other.coefficients_, other.cacheId_)
      {}
 
      LocalOperator(LocalOperator&& other) :
        LocalOperator(std::move(other.localUnaryOperator_), other.usedNodes_, other.basis_, other.coefficients_, other.cacheId_)
      {}
 
      void bind(const Element& element)
      {
        if (externalLocalView_[0] == nullptr)
        {
          internalLocalView_[0].bind(element);
          loadLocalDOFs(internalLocalView_[0], localDoFs_[0]);
        }
        else
          loadLocalDOFs(*(externalLocalView_[0]), localDoFs_[0]);
        localUnaryOperator_.bind(element);
      }
 
      void bind(const Intersection& intersection, const Element& element, const Element& otherElement)
      {
        if (externalLocalView_[0] == nullptr)
        {
          internalLocalView_[0].bind(element);
          loadLocalDOFs(internalLocalView_[0], localDoFs_[0]);
        }
        else
          loadLocalDOFs(*(externalLocalView_[0]), localDoFs_[0]);
        if (externalLocalView_[1] == nullptr)
        {
          internalLocalView_[1].bind(otherElement);
          loadLocalDOFs(internalLocalView_[1], localDoFs_[1]);
        }
        else
          loadLocalDOFs(*(externalLocalView_[1]), localDoFs_[1]);
        localUnaryOperator_.bind(intersection, element, otherElement);
      }
 
      void unbind()
      {
        if (externalLocalView_[0] == nullptr)
          internalLocalView_[0].unbind();
        if (externalLocalView_[1] == nullptr)
          internalLocalView_[1].unbind();
      }
 
      // Additional interface of local nullary operator
      auto quadratureRuleKey() const
      {
        return localUnaryOperator_.quadratureRuleKey();
      }
 
      auto operator()(std::size_t index) const
      {
        auto& valueCache = cacheManager_->template getCache<ValueCache>(cacheIndex_);
        auto& values = valueCache.getValues();
        if (valueCache.isEmpty())
        {
          const auto& rule = valueCache.rule();
          for (auto k : Dune::range(rule.size()))
          {
            auto& y = values[k];
            AssignOp()(y, 0);
            forEachOperator(localUnaryOperator_, [&](auto& op) {
              sparseForEach(op(k), [&](const auto& value, auto stencilIndex, auto localIndex) {
                SumOp::addTo(y, MultOp()(localDoFs_[stencilIndex][localIndex], value));
              });
            });
          }
          valueCache.setEmpty(false);
        }
        return Impl::Tensor::RankZeroTensor(values[index]);
      }
 
      template<class... LV>
      void registerLocalViews(const LV&... lvs)
      {
        localUnaryOperator_.registerLocalViews(lvs...);
        Impl::visitMatchingLocalView([&](const auto& localView) {
          externalLocalView_[0] = &localView;
          localDoFs_[0].reserve(externalLocalView_[0]->maxSize());
        }, basis_, lvs...);
      }
 
      template<class... LV>
      void registerOutsideLocalViews(const LV&... lvs)
      {
        localUnaryOperator_.registerOutsideLocalViews(lvs...);
        Impl::visitMatchingLocalView([&](const auto& localView) {
          externalLocalView_[1] = &localView;
          localDoFs_[1].reserve(externalLocalView_[1]->maxSize());
        }, basis_, lvs...);
      }
 
      void registerCaches(CacheManager& cacheManager)
      {
        localUnaryOperator_.registerCaches(cacheManager);
        cacheIndex_ = cacheManager.registerCache(cacheId_, ValueCache(true));
      }
 
      void registerOutsideCaches(CacheManager& cacheManager)
      {
        localUnaryOperator_.registerOutsideCaches(cacheManager);
      }
 
      template<class... OutsideCacheManager>
      void bindToCaches(CacheManager& cacheManager, OutsideCacheManager&... outsideCacheManager)
      {
        localUnaryOperator_.bindToCaches(cacheManager, outsideCacheManager...);
        cacheManager_ = &cacheManager;
      }
 
    private:
      const Basis& basis_;
      mutable LocalUnaryOperator localUnaryOperator_;
      NodeFlags usedNodes_;
      std::shared_ptr<const Vector> coefficients_;
      mutable std::array<std::vector<Coefficient>,2> localDoFs_;
      std::array<LocalView, 2> internalLocalView_;
      std::array<const LocalView*,2> externalLocalView_ = {nullptr, nullptr};
      mutable CacheManager* cacheManager_ = nullptr;
      UniqueCacheId cacheId_;
      std::size_t cacheIndex_ = 0;
    };
 
    auto basis() const
    {
      return std::tuple<>();
    }
 
    auto treePath() const
    {
      return std::tuple<>();
    }
 
    friend LocalOperator localOperator(const BoundUnaryOperator& boundUnaryOperator)
    {
      const auto& rootBasis = std::get<0>(boundUnaryOperator.unaryOperator_.basis());
      // Mark all leaf nodes whose DOFs need to be loaded
      auto localView = rootBasis.localView();
      auto usedNodes = TypeTree::makeTreeContainer(localView.tree(), [](auto&&) { return char(false); });
      forEachOperator(boundUnaryOperator.unaryOperator_, [&](auto& op) {
        const auto& treePath = std::get<0>(op.treePath());
        const auto& node = Dune::TypeTree::child(localView.tree(), treePath);
        TypeTree::forEachLeafNode(node, [&](auto&& leafNode, auto&& leafTreePath) {
          usedNodes[TypeTree::join(treePath,leafTreePath)] = true;
        });
      });
      return LocalOperator(localOperator(boundUnaryOperator.unaryOperator_), usedNodes, rootBasis, boundUnaryOperator.coefficients_, boundUnaryOperator.cacheId_);
    }
 
    friend LocalFunctionAdaptor<LocalOperator> localFunction(const BoundUnaryOperator& boundUnaryOperator)
    requires (BoundUnaryOperator::arity==0)
    {
      return LocalFunctionAdaptor<LocalOperator>(localOperator(boundUnaryOperator));
    }
 
    Range operator()(const Domain& x) const
    {
      DUNE_THROW(Dune::NotImplemented, "Evaluation of BoundUnaryOperator in global coordinates is not implemented");
    }
 
    const EntitySet& entitySet() const
    {
      return entitySet_;
    }
 
    const auto& unaryOperator() const
    {
      return unaryOperator_;
    }
 
    const auto& coefficients() const
    {
      return coefficients_;
    }
 
    friend auto jacobian(const BoundUnaryOperator& f)
    {
      static_assert(isOperator_v<UnaryOperator>, "jacobian(BoundUnaryOperator<SumOperator<...>,V>(...)) is not implemented.");
      return Dune::Fufem::Forms::BoundUnaryOperator(jacobian(f.unaryOperator_), f.coefficients_);
    }
 
    friend auto gradient(const BoundUnaryOperator& f)
    {
      static_assert(isOperator_v<UnaryOperator>, "gradient(BoundUnaryOperator<SumOperator<...>,V>(...)) is not implemented.");
      return Dune::Fufem::Forms::BoundUnaryOperator(gradient(f.unaryOperator_), f.coefficients_);
    }
 
    friend auto grad(const BoundUnaryOperator& f)
    {
      static_assert(isOperator_v<UnaryOperator>, "grad(BoundUnaryOperator<SumOperator<...>,V>(...)) is not implemented.");
      return Dune::Fufem::Forms::BoundUnaryOperator(grad(f.unaryOperator_), f.coefficients_);
    }
 
    friend auto divergence(const BoundUnaryOperator& f)
    {
      static_assert(isOperator_v<UnaryOperator>, "divergence(BoundUnaryOperator<SumOperator<...>,V>(...)) is not implemented.");
      return Dune::Fufem::Forms::BoundUnaryOperator(divergence(f.unaryOperator_), f.coefficients_);
    }
 
    friend auto div(const BoundUnaryOperator& f)
    {
      static_assert(isOperator_v<UnaryOperator>, "div(BoundUnaryOperator<SumOperator<...>,V>(...)) is not implemented.");
      return Dune::Fufem::Forms::BoundUnaryOperator(div(f.unaryOperator_), f.coefficients_);
    }
 
  private:
    UniqueCacheId cacheId_;
    UnaryOperator unaryOperator_;
    std::shared_ptr<const Vector> coefficients_;
    EntitySet entitySet_;
  };
 
 
 
} // namespace Dune::Fufem::Forms
 
 
#endif // DUNE_FUFEM_FORMS_BOUNDUNARYOPERATOR_HH