unaryoperators.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_UNARYOPERATORS_HH
#define DUNE_FUFEM_FORMS_UNARYOPERATORS_HH
 
#include <cstddef>
#include <type_traits>
#include <vector>
 
#include <dune/common/typetraits.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/indices.hh>
#include <dune/common/tupleutility.hh>
#include <dune/common/typelist.hh>
#include <dune/common/version.hh>
 
#if DUNE_VERSION_GTE(DUNE_FUNCTIONS, 2, 11)
#include <dune/common/typetree/childaccess.hh>
#include <dune/common/typetree/treepath.hh>
#else
#include <dune/typetree/childextraction.hh>
#include <dune/typetree/treepath.hh>
#endif
 
#include <dune/functions/functionspacebases/subspacebasis.hh>
 
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
 
#include <dune/fufem/forms/baseclass.hh>
#include <dune/fufem/forms/tensors.hh>
#include <dune/fufem/forms/shapefunctioncache.hh>
 
#include <dune/fufem/typetreefallback.hh>
#include <dune/fufem/utilities/asnumber.hh>
 
 
 
namespace Dune::Fufem::Forms {
 
 
 
  template<class B, class TP, std::size_t argIndex>
  class FEOperatorBase : public UnaryOperator<argIndex>
  {
 
  protected:
    using LocalView = typename B::LocalView;
    using Tree = typename LocalView::Tree;
    using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TP>;
    using SubspaceBasis = typename Dune::Functions::SubspaceBasis<B, TP>;
 
    // Helper function to derive the leafTreePath associated
    // to a treePath. For treePaths referring to a leaf this is
    // the identity. For treePaths referring to the a power
    // node whose children are leaf, it's the first child.
    // This is used to implement vector valued ansatz functions
    // for power nodes.
    static auto leafTreePath(const TP& tp)
    {
      if constexpr(Dune::Fufem::Impl::Concept::LeafTreeNode<Node>)
        return tp;
      else if constexpr(Dune::Fufem::Impl::Concept::UniformInnerTreeNode<Node>
          and Dune::Fufem::Impl::Concept::LeafTreeNode<Dune::TypeTree::Child<Node, 0>>)
        return Dune::TypeTree::push_back(tp, Dune::Indices::_0);
    }
 
    using LeafTreePath = decltype(leafTreePath(std::declval<TP>()));
    using LeafNode = typename Dune::TypeTree::ChildForTreePath<Tree, LeafTreePath>;
    using Basis = B;
    using TreePath = TP;
 
  public:
 
    using Element = typename Basis::LocalView::Element;
    using Intersection = typename Basis::GridView::Intersection;
 
    FEOperatorBase(const Basis& basis, const TreePath& treePath, bool isAffine = true) :
      subspaceBasis_(basis, treePath),
      isAffine_(isAffine)
    {}
 
    class LocalOperator
    {
    protected:
 
      using TreeCache = ShapeFunctionCache<Tree>;
      using LeafNodeCache = ShapeFunctionCache<LeafNode>;
 
    public:
 
      using Element = typename FEOperatorBase::Element;
      using Intersection = typename FEOperatorBase::Intersection;
      using CacheManager = Dune::Fufem::Forms::CacheManager<typename Element::Geometry::ctype, Element::Geometry::mydimension>;
 
      LocalOperator(const SubspaceBasis& subspaceBasis, bool isAffine) :
        subspaceBasis_(subspaceBasis),
        leafTreePath_(leafTreePath(subspaceBasis_.prefixPath())),
        leafNode_(nullptr),
        localView_(nullptr),
        quadratureRuleKey_(),
        cacheIndex_(std::numeric_limits<std::size_t>::max()),
        isAffine_(isAffine)
      {}
 
      LocalOperator(const FEOperatorBase& op) :
        LocalOperator(op.subspaceBasis_, op.isAffine_)
      {}
 
      void unbind()
      {}
 
      auto quadratureRuleKey() const
      {
        return quadratureRuleKey_;
      }
 
      template<class... LV>
      void registerLocalViews(const LV&... lvs)
      {
        Impl::visitMatchingLocalView([&](const auto& localView) {
          localView_ = &localView;
          leafNode_ = & Dune::TypeTree::child(localView.tree(), leafTreePath_);
        }, subspaceBasis_.rootBasis(), lvs...);
      }
 
      template<class... LV>
      void registerOutsideLocalViews(const LV&... lvs)
      {}
 
      void registerCaches(CacheManager& cacheManager)
      {
        const auto& tree = localView_->tree();
        cacheIndex_ = cacheManager.registerCache(UniqueCacheId(tree), TreeCache(tree));
        if (not isAffine_)
          cacheManager.template getCache<TreeCache>(cacheIndex_)[subspaceBasis_.prefixPath()].setNonAffine();
      }
 
      void registerOutsideCaches(CacheManager& cacheManager)
      {}
 
    protected:
 
      const LeafNode& leafNode() const
      {
        return *leafNode_;
      }
 
      auto& leafNodeCache(CacheManager& cacheManager)
      {
        return cacheManager.template getCache<TreeCache>(cacheIndex_)[leafTreePath_];
      }
 
      SubspaceBasis subspaceBasis_;
      const LeafTreePath leafTreePath_;
      const LeafNode* leafNode_;
      const LocalView* localView_;
      QuadratureRuleKey quadratureRuleKey_;
      std::size_t cacheIndex_;
      bool isAffine_;
    };
 
    auto basis() const
    {
      return std::tie(subspaceBasis_.rootBasis());
    }
 
    auto treePath() const
    {
      return std::tie(subspaceBasis_.prefixPath());
    }
 
    bool isAffine() const
    {
      return isAffine_;
    }
 
    template<std::size_t k>
    FEOperatorBase<B, TP, k> rebindArgIndex() const
    {
      return {std::get<0>(basis()), std::get<0>(treePath()), isAffine_};
    }
 
  protected:
    SubspaceBasis subspaceBasis_;
    bool isAffine_;
  };
 
 
 
  template<class B, class TP, std::size_t argIndex>
  class FEFunctionOperator : public FEOperatorBase<B, TP, argIndex>
  {
    using Base = FEOperatorBase<B, TP, argIndex>;
    using Node = typename Base::Node;
    using LeafLBTraits = typename Base::LeafNode::FiniteElement::Traits::LocalBasisType::Traits;
 
  public:
 
    using Element = typename Base::Element;
    using Range = std::conditional_t<
      Dune::Fufem::Impl::Concept::LeafTreeNode<Node>,
      typename LeafLBTraits::RangeType,
      Dune::FieldVector<typename LeafLBTraits::RangeFieldType, Node::degree()>>;
 
    using Base::Base;
 
    FEFunctionOperator(const Base& other)
      : Base(other)
    {}
 
    class LocalOperator : public Base::LocalOperator
    {
      using Base::LocalOperator::quadratureRuleKey_;
      using Base::LocalOperator::leafNode;
      using Base::LocalOperator::leafNodeCache;
      using Base::LocalOperator::LocalOperator;
 
    public:
 
      using Intersection = typename Base::LocalOperator::Intersection;
      using CacheManager = typename Base::LocalOperator::CacheManager;
      using Range = typename FEFunctionOperator::Range;
 
      void bind(const Element&)
      {
        if (leafNode().size()!=0)
          quadratureRuleKey_ = QuadratureRuleKey(leafNode().finiteElement());
        else
          quadratureRuleKey_ = QuadratureRuleKey(this->localView_->element().type(), 0);
      }
 
      void bind(const Intersection& intersection, const Element& element, const Element& otherElement)
      {
        bind(element);
      }
 
      template<class... OutsideCacheManager>
      void bindToCaches(CacheManager& cacheManager, OutsideCacheManager&... outsideCacheManager)
      {
        valueCache_ = &(leafNodeCache(cacheManager).getValues());
      }
 
      auto operator()(std::size_t index) const
      {
        using namespace Impl::Tensor;
        using namespace Dune::Indices;
        const auto& values = (*valueCache_)[index];
        std::size_t size = leafNode().size();
        std::size_t indexOffset = (size==0) ? 0 : leafNode().localIndex(0);
        if constexpr(Dune::Fufem::Impl::Concept::LeafTreeNode<Node>)
          return SparseTensor(_2,
            [&, size, indexOffset](auto&& f) {
              for(std::size_t i=0; i<size; ++i)
                f(values[i], _0, indexOffset+i);
            },
            size,
            Dune::MetaType<Range>{}
          );
        else if constexpr(Dune::Fufem::Impl::Concept::UniformInnerTreeNode<Node>
            and Dune::Fufem::Impl::Concept::LeafTreeNode<Dune::TypeTree::Child<Node, 0>>)
          return SparseTensor(_2,
            [&, size, indexOffset](auto&& f) {
              for(std::size_t j=0; j<Node::degree(); ++j)
                for(std::size_t i=0; i<size; ++i)
                {
                  auto result = Range(0);
                  result[j] = Dune::Fufem::asNumber(values[i]);
                  f(result, _0, indexOffset + j*size + i);
                }
            },
            Node::degree() * size,
            Dune::MetaType<Range&&>{}
          );
      }
 
    private:
      const typename Base::LocalOperator::LeafNodeCache::ValueCache* valueCache_ = nullptr;
    };
 
    friend LocalOperator localOperator(const FEFunctionOperator& op)
    {
      return LocalOperator(op);
    }
 
    template<std::size_t i>
    auto childOperator(Dune::index_constant<i> childIndex) const
    {
      using Basis = typename Base::Basis;
      auto childTP = Dune::TypeTree::push_back(Base::subspaceBasis_.prefixPath(), childIndex);
      return FEFunctionOperator<Basis, decltype(childTP), argIndex>(Base::subspaceBasis_.rootBasis(), childTP, Base::isAffine_);
    }
 
    template<std::size_t k>
    FEFunctionOperator<B, TP, k> rebindArgIndex() const
    {
      return {Base::template rebindArgIndex<k>()};
    }
 
  };
 
 
 
 
  template<class B, class TP, std::size_t argIndex>
  class FEFunctionJacobianOperator : public FEOperatorBase<B, TP, argIndex>
  {
    using Base = FEOperatorBase<B, TP, argIndex>;
    using Node = typename Base::Node;
    using LeafLBTraits = typename Base::LeafNode::FiniteElement::Traits::LocalBasisType::Traits;
 
    using ReferenceJacobian = std::conditional_t<
      Dune::Fufem::Impl::Concept::LeafTreeNode<Node>,
      typename LeafLBTraits::JacobianType,
      Dune::FieldMatrix<typename LeafLBTraits::RangeFieldType, Node::degree(), LeafLBTraits::dimDomain>>;
 
    using GeometryJacobianInverse = typename Node::Element::Geometry::JacobianInverse;
 
  public:
 
    using Element = typename Base::Element;
    using Range = typename Impl::GlobalJacobianTraits<ReferenceJacobian, GeometryJacobianInverse>::type;
 
    using Base::Base;
 
    FEFunctionJacobianOperator(const Base& other)
      : Base(other)
    {}
 
    class LocalOperator : public Base::LocalOperator
    {
      using Base::LocalOperator::quadratureRuleKey_;
      using Base::LocalOperator::leafNode;
      using Base::LocalOperator::leafNodeCache;
      using Base::LocalOperator::LocalOperator;
 
    public:
 
      using Intersection = typename Base::LocalOperator::Intersection;
      using CacheManager = typename Base::LocalOperator::CacheManager;
      using Range = typename FEFunctionJacobianOperator::Range;
 
      void bind(const Element&)
      {
        if (leafNode().size()!=0)
          quadratureRuleKey_ = QuadratureRuleKey(leafNode().finiteElement()).derivative();
        else
          quadratureRuleKey_ = QuadratureRuleKey(this->localView_->element().type(), 0);
      }
 
      void bind(const Intersection& intersection, const Element& element, const Element& otherElement)
      {
        bind(element);
      }
 
      template<class... OutsideCacheManager>
      void bindToCaches(CacheManager& cacheManager, OutsideCacheManager&... outsideCacheManager)
      {
        globalJacobianCache_ = &(leafNodeCache(cacheManager).getGlobalJacobians());
      }
 
      auto operator()(std::size_t index) const
      {
        using namespace Impl::Tensor;
        using namespace Dune::Indices;
        const auto& globalJacobians = (*globalJacobianCache_)[index];
        std::size_t size = leafNode().size();
        std::size_t indexOffset = (size==0) ? 0 : leafNode().localIndex(0);
        if constexpr(Dune::Fufem::Impl::Concept::LeafTreeNode<Node>)
          return SparseTensor(_2,
            [&, size, indexOffset](auto&& f) {
              for(std::size_t i=0; i<size; ++i)
                f(globalJacobians[i], _0, indexOffset+i);
            },
            size,
            Dune::MetaType<Range>{}
          );
        else if constexpr(Dune::Fufem::Impl::Concept::UniformInnerTreeNode<Node>
            and Dune::Fufem::Impl::Concept::LeafTreeNode<Dune::TypeTree::Child<Node, 0>>)
          return SparseTensor(_2,
            [&, size, indexOffset](auto&& f) {
              for(std::size_t j=0; j<Node::degree(); ++j)
                for(std::size_t i=0; i<size; ++i)
                {
                  auto result = Range(0);
                  result[j] = globalJacobians[i][0];
                  f(result, _0, indexOffset + j*size + i);
                }
            },
            Node::degree() * size,
            Dune::MetaType<Range&&>{}
          );
      }
 
    private:
      const typename Base::LocalOperator::LeafNodeCache::GlobalJacobianCache* globalJacobianCache_ = nullptr;
    };
 
    friend LocalOperator localOperator(const FEFunctionJacobianOperator& op)
    {
      return LocalOperator(op);
    }
 
    template<std::size_t i>
    auto childOperator(Dune::index_constant<i> childIndex) const
    {
      using Basis = typename Base::Basis;
      auto childTP = Dune::TypeTree::push_back(Base::subspaceBasis_.prefixPath(), childIndex);
      return FEFunctionJacobianOperator<Basis, decltype(childTP), argIndex>(Base::subspaceBasis_.rootBasis(), childTP, Base::isAffine_);
    }
 
    template<std::size_t k>
    FEFunctionJacobianOperator<B, TP, k> rebindArgIndex() const
    {
      return {Base::template rebindArgIndex<k>()};
    }
 
  };
 
 
 
  template<class B, class TP, std::size_t argIndex>
  class FEFunctionDivergenceOperator : public FEOperatorBase<B, TP, argIndex>
  {
    using Base = FEOperatorBase<B, TP, argIndex>;
    using Node = typename Base::Node;
    using LeafLBTraits = typename Base::LeafNode::FiniteElement::Traits::LocalBasisType::Traits;
 
  public:
 
    using Element = typename Base::Element;
    using Range = std::conditional_t<
      Dune::Fufem::Impl::Concept::LeafTreeNode<Node>,
      typename LeafLBTraits::RangeFieldType,
      typename LeafLBTraits::JacobianType::block_type>;
 
    using Base::Base;
 
    FEFunctionDivergenceOperator(const Base& other)
      : Base(other)
    {}
 
    class LocalOperator : public Base::LocalOperator
    {
      using Base::LocalOperator::quadratureRuleKey_;
      using Base::LocalOperator::leafNode;
      using Base::LocalOperator::leafNodeCache;
      using Base::LocalOperator::LocalOperator;
 
    public:
 
      using Intersection = typename Base::LocalOperator::Intersection;
      using CacheManager = typename Base::LocalOperator::CacheManager;
      using Range = typename FEFunctionDivergenceOperator::Range;
 
      void bind(const Element&)
      {
        if (leafNode().size()!=0)
          quadratureRuleKey_ = QuadratureRuleKey(leafNode().finiteElement()).derivative();
        else
          quadratureRuleKey_ = QuadratureRuleKey(this->localView_->element().type(), 0);
      }
 
      void bind(const Intersection& intersection, const Element& element, const Element& otherElement)
      {
        bind(element);
      }
 
      template<class... OutsideCacheManager>
      void bindToCaches(CacheManager& cacheManager, OutsideCacheManager&... outsideCacheManager)
      {
        globalJacobianCache_ = &(leafNodeCache(cacheManager).getGlobalJacobians());
      }
 
      auto operator()(std::size_t index) const
      {
        using namespace Impl::Tensor;
        using namespace Dune::Indices;
        using LocalBasisTraits = typename Base::LeafNode::FiniteElement::Traits::LocalBasisType::Traits;
        using Field = typename LocalBasisTraits::RangeFieldType;
        const auto& globalJacobians = (*globalJacobianCache_)[index];
        std::size_t size = leafNode().size();
        std::size_t indexOffset = (size==0) ? 0 : leafNode().localIndex(0);
        if constexpr(Dune::Fufem::Impl::Concept::LeafTreeNode<Node>)
          return SparseTensor(_2,
            [&, size, indexOffset](auto&& f) {
              for(std::size_t i=0; i<size; ++i)
              {
                auto div = Field{0};
                for(std::size_t j=0; j<globalJacobians[i].N(); ++j)
                  div += globalJacobians[i][j][j];
                f(div, _0, indexOffset+i);
              }
            },
            size,
            Dune::MetaType<Range&&>{}
          );
        else if constexpr(Dune::Fufem::Impl::Concept::UniformInnerTreeNode<Node>
            and Dune::Fufem::Impl::Concept::LeafTreeNode<Dune::TypeTree::Child<Node, 0>>)
          return SparseTensor(_2,
            [&, size, indexOffset](auto&& f) {
              // Here we should compute the trace of the global Jacobian J
              // that we would need to build first. But since only the
              // j-th row of J is nonzero, we can simply
              // return the diagonal entry of this row which coincides
              // with the j-th entry of the gradient of
              // the respective scalar basis function.
              for(std::size_t j=0; j<Node::degree(); ++j)
                for(std::size_t i=0; i<size; ++i)
                  f(globalJacobians[i][0][j], _0, indexOffset + j*size + i);
            },
            Node::degree() * size,
            Dune::MetaType<Range>{}
          );
      }
 
    private:
      const typename Base::LocalOperator::LeafNodeCache::GlobalJacobianCache* globalJacobianCache_ = nullptr;
    };
 
    friend LocalOperator localOperator(const FEFunctionDivergenceOperator& op)
    {
      return LocalOperator(op);
    }
 
    template<std::size_t k>
    FEFunctionDivergenceOperator<B, TP, k> rebindArgIndex() const
    {
      return {Base::template rebindArgIndex<k>()};
    }
 
  };
 
 
 
} // namespace Dune::Fufem::Forms
 
 
 
#endif // DUNE_FUFEM_FORMS_UNARYOPERATORS_HH