shapefunctioncache.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_SHAPEFUNCTIONCACHE_HH
#define DUNE_FUFEM_FORMS_SHAPEFUNCTIONCACHE_HH
 
#include <type_traits>
#include <utility>
#include <list>
#include <any>
#include <typeindex>
 
#include <dune/common/indices.hh>
#include <dune/common/tuplevector.hh>
#include <dune/common/typetree/nodeconcepts.hh>
#include <dune/common/typetree/treepath.hh>
 
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/fufem/quadraturerules/subentityquadraturerule.hh>
 
 
namespace Dune::Fufem::Forms::Impl {
 
  template<class C, class... T>
  static constexpr decltype(auto) accessByTreePath(C&& container, const Dune::TypeTree::TreePath<T...>& path)
  {
    if constexpr (sizeof...(T)==0)
      return container;
    else
      return accessByTreePath(container[path.front()], pop_front(path));
  }
 
  template<class ReferenceJacobian, class GeometryJacobianInverse>
  struct GlobalJacobianTraits
  {
    using GeometryField = typename GeometryJacobianInverse::field_type;
    using ReferenceJacobianField = typename ReferenceJacobian::field_type;
    using JacobianField = typename Dune::PromotionTraits<GeometryField, ReferenceJacobianField>::PromotedType;
    static constexpr int dimWorld = GeometryJacobianInverse::cols;
    using type = decltype(std::declval<ReferenceJacobian>()*std::declval<GeometryJacobianInverse>());
  };
 
  template<class GeometryJacobianInverse, class JF, int dimRange, int dim>
  struct GlobalJacobianTraits<Dune::FieldMatrix<JF, dimRange, dim>, GeometryJacobianInverse>
  {
    using GeometryField = typename GeometryJacobianInverse::field_type;
    using ReferenceJacobianField = JF;
    using JacobianField = typename Dune::PromotionTraits<GeometryField, ReferenceJacobianField>::PromotedType;
    static constexpr int dimWorld = GeometryJacobianInverse::cols;
    using type = Dune::FieldMatrix<JacobianField, dimRange, dimWorld>;
  };
 
 
 
} // end namespace Dune::Fufem::Forms::Impl
 
 
 
namespace Dune::Fufem::Forms {
 
 
 
  template<class Node, class CT=double>
  class ShapeFunctionCache;
 
 
 
  template<class Node, class CT>
  requires Dune::TypeTree::Concept::LeafTreeNode<Node>
  class ShapeFunctionCache<Node, CT>
  {
    using GeometryJacobianInverse = typename Node::Element::Geometry::JacobianInverse;
  public:
 
    using QuadratureRule = Dune::QuadratureRule<CT, Node::Element::dimension>;
 
    using FEValue = typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
    using FEJacobian = typename Node::FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
    using FEGlobalJacobian = typename Impl::GlobalJacobianTraits<FEJacobian, GeometryJacobianInverse>::type;
 
    using ValueCache = typename std::vector<std::vector<FEValue>>;
    using JacobianCache = typename std::vector<std::vector<FEJacobian>>;
    using GlobalJacobianCache = typename std::vector<std::vector<FEGlobalJacobian>>;
 
    ShapeFunctionCache(const QuadratureRule& rule, const Node& node)
    {
      setRule(rule);
      setTree(node);
    }
 
    ShapeFunctionCache(const QuadratureRule& rule) :
      node_(nullptr)
    {
      setRule(rule);
    }
 
    ShapeFunctionCache(const Node& node) :
      rule_(nullptr)
    {
      setTree(node);
    }
 
    ShapeFunctionCache(const ShapeFunctionCache& other) = default;
 
    ShapeFunctionCache() = default;
 
    void setRule(const QuadratureRule& rule)
    {
      rule_ = &rule;
      valueCache_.resize(rule_->size());
      jacobianCache_.resize(rule_->size());
      globalJacobianCache_.resize(rule_->size());
      valueCache_[0].clear();
      jacobianCache_[0].clear();
      globalJacobianCache_[0].clear();
    }
 
    void setTree(const Node& node)
    {
      node_ = &node;
    }
 
    void setNonAffine()
    {
      isNonAffine_ = true;
    }
 
    const QuadratureRule& rule() const
    {
      return *rule_;
    }
 
    void invalidate()
    {
      // We use the size at the first quadrature point to indicate
      // an invalidated cache. To this end we set it's size to zero.
      if (isNonAffine_)
      {
        valueCache_[0].clear();
        jacobianCache_[0].clear();
      }
      globalJacobianCache_[0].clear();
    }
 
    decltype(auto) operator[](const Dune::TypeTree::TreePath<>& treePath) const
    {
      return *this;
    }
 
    decltype(auto) operator[](const Dune::TypeTree::TreePath<>& treePath)
    {
      return *this;
    }
 
    const auto& getValues()
    {
      if (valueCache_[0].empty())
      {
        if (node_->size() == 0)
          return valueCache_;
        const auto& localBasis = node_->finiteElement().localBasis();
        for(std::size_t k=0; k<rule_->size(); ++k)
          localBasis.evaluateFunction((*rule_)[k].position(), valueCache_[k]);
      }
      return valueCache_;
    }
 
    const auto& getJacobians()
    {
      if (jacobianCache_[0].empty())
      {
        if (node_->size() == 0)
          return jacobianCache_;
        const auto& localBasis = node_->finiteElement().localBasis();
        for(std::size_t k=0; k<rule_->size(); ++k)
          localBasis.evaluateJacobian((*rule_)[k].position(), jacobianCache_[k]);
      }
      return jacobianCache_;
    }
 
    const auto& getGlobalJacobians()
    {
      if (globalJacobianCache_[0].empty())
      {
        if (node_->size() == 0)
          return globalJacobianCache_;
        const auto& geometry = node_->element().geometry();
        const auto& jacobians = getJacobians();
        for(std::size_t k=0; k<rule_->size(); ++k)
        {
          globalJacobianCache_[k].resize(jacobians[k].size());
          const auto& jacobianInverse = geometry.jacobianInverse((*rule_)[k].position());
          for(std::size_t i=0; i<jacobians[k].size(); ++i)
            globalJacobianCache_[k][i] = jacobians[k][i] * jacobianInverse;
        }
      }
      return globalJacobianCache_;
    }
 
  private:
    const QuadratureRule* rule_ = nullptr;
    const Node* node_ = nullptr;
    ValueCache valueCache_;
    JacobianCache jacobianCache_;
    GlobalJacobianCache globalJacobianCache_;
    bool isNonAffine_ = false;
  };
 
 
 
  template<class Node, class CT>
  requires Dune::TypeTree::Concept::UniformInnerTreeNode<Node>
  class ShapeFunctionCache<Node, CT>
  {
    using ChildNode = Dune::TypeTree::Child<Node,0>;
    using ChildCache = ShapeFunctionCache<ChildNode, CT>;
  public:
 
    using QuadratureRule = Dune::QuadratureRule<CT, Node::Element::dimension>;
 
    ShapeFunctionCache(const QuadratureRule& rule, const Node& node) :
      childCache_(rule, node.child(0))
    {}
 
    ShapeFunctionCache(const QuadratureRule& rule) :
      childCache_(rule)
    {}
 
    ShapeFunctionCache(const Node& node) :
      childCache_(node.child(0))
    {}
 
    ShapeFunctionCache(const ShapeFunctionCache& other) = default;
 
    ShapeFunctionCache() = default;
 
    void setRule(const QuadratureRule& rule)
    {
      childCache_.setRule(rule);
    }
 
    void setTree(const Node& node)
    {
      childCache_.setTree(node.child(0));
    }
 
    void setNonAffine()
    {
      childCache_.setNonAffine();
    }
 
    const QuadratureRule& rule() const
    {
      return childCache_.rule();
    }
 
    void invalidate()
    {
      childCache_.invalidate();
    }
 
    const ChildCache& operator[](std::size_t i) const
    {
      assert(i<Node::degree());
      return childCache_;
    }
 
    ChildCache& operator[](std::size_t i)
    {
      assert(i<Node::degree());
      return childCache_;
    }
 
    template<class... T>
    decltype(auto) operator[](const Dune::TypeTree::TreePath<T...>& treePath) const
    {
      return Impl::accessByTreePath(*this, treePath);
    }
 
    template<class... T>
    decltype(auto) operator[](const Dune::TypeTree::TreePath<T...>& treePath)
    {
      return Impl::accessByTreePath(*this, treePath);
    }
 
    constexpr std::size_t size() const noexcept
    {
      return Node::degree();
    }
 
  private:
    ChildCache childCache_;
  };
 
 
 
  namespace Impl {
 
    template<class Node, class CT, class Indices>
    struct TupleVectorOfChildShapeFunctionCaches
    {};
 
    template<class Node, class CT, std::size_t... indices>
    struct TupleVectorOfChildShapeFunctionCaches<Node, CT, std::index_sequence<indices...>>
    {
      using type = Dune::TupleVector<ShapeFunctionCache<Dune::TypeTree::Child<Node, indices>, CT>...>;
    };
 
    template<class Node>
    using ChildIndices = std::make_index_sequence<decltype(Node::degree())::value>;
 
    template<class Node, class CT>
    using TupleVectorOfChildShapeFunctionCaches_t = typename TupleVectorOfChildShapeFunctionCaches<Node, CT, ChildIndices<Node>>::type;
 
  } // end namespace Imp
 
 
 
  template<class Node, class CT>
  requires (Dune::TypeTree::Concept::StaticDegreeInnerTreeNode<Node>
            and not Dune::TypeTree::Concept::UniformInnerTreeNode<Node>)
  class ShapeFunctionCache<Node, CT>
    : public Impl::TupleVectorOfChildShapeFunctionCaches_t<Node, CT>
  {
    using Base = Impl::TupleVectorOfChildShapeFunctionCaches_t<Node, CT>;
  public:
 
    using QuadratureRule = Dune::QuadratureRule<CT, Node::Element::dimension>;
 
    ShapeFunctionCache(const QuadratureRule& rule, const Node& node)
    {
      setRule(rule);
      setTree(node);
    }
 
    ShapeFunctionCache(const QuadratureRule& rule)
    {
      setRule(rule);
    }
 
    ShapeFunctionCache(const Node& node)
    {
      setTree(node);
    }
 
 
    ShapeFunctionCache(const ShapeFunctionCache& other) = default;
 
    ShapeFunctionCache() = default;
 
    void setRule(const QuadratureRule& rule)
    {
      Hybrid::forEach(*this, [&](auto& childCache) {
        childCache.setRule(rule);
      });
    }
 
    void setTree(const Node& node)
    {
      Hybrid::forEach(Dune::range(Node::degree()), [&](const auto& i){
        (*this)[i].setTree(node.child(i));
      });
    }
 
    void setNonAffine()
    {
      Hybrid::forEach(Dune::range(Node::degree()), [&](const auto& i){
        (*this)[i].setNonAffine();
      });
    }
 
    const QuadratureRule& rule() const
    {
      return (*this)[Dune::Indices::_0].rule();
    }
 
    using Base::operator[];
 
    template<class... T>
    decltype(auto) operator[](const Dune::TypeTree::TreePath<T...>& treePath) const
    {
      return Impl::accessByTreePath(*this, treePath);
    }
 
    template<class... T>
    decltype(auto) operator[](const Dune::TypeTree::TreePath<T...>& treePath)
    {
      return Impl::accessByTreePath(*this, treePath);
    }
 
    void invalidate() {
      Hybrid::forEach(*this, [&](auto& childCache) {
        childCache.invalidate();
      });
    }
 
  };
 
 
 
  template<class C>
  class MultipleQuadratureCache
  {
    using Cache = C;
    using CT = typename Cache::QuadratureRule::CoordType;
    constexpr static int dimension = Cache::QuadratureRule::d;
 
    using QuadratureRule = Dune::QuadratureRule<CT, dimension>;
    using FacetQuadratureRule = Dune::Fufem::SubEntityQuadratureRule<CT , dimension, 1>;
 
  public:
 
    using ElementKey = QuadratureRuleKey;
    using FacetKey = std::pair<QuadratureRuleKey, std::size_t>;
 
  private:
 
    template<class Key>
    const QuadratureRule& getRule(const Key& key)
    {
      if constexpr (std::is_same_v<Key, ElementKey>)
        return QuadratureRuleCache<CT, dimension>::rule(key);
      else
      {
        auto facet = key.second;
        auto&& type = key.first.geometryType();
        auto facetGeometryType = Dune::referenceElement<CT,dimension>(type).type(facet, 1);
        auto facetQuadratureRuleKey = key.first;
        facetQuadratureRuleKey.setGeometryType(facetGeometryType);
        const auto& rule = QuadratureRuleCache<CT, dimension-1>::rule(facetQuadratureRuleKey);
        // We need to store the face quadrature rule.
        // By using an std::list, we can store pointers
        // to the stored rules, because they never get invalidated.
        facetQuadratureRules_.push_back(FacetQuadratureRule(type, facet, rule));
        return facetQuadratureRules_.back();
      }
    }
 
    template<class Key>
    Cache makeCache(const Key& key)
    {
      auto cache = Cache(prototype_);
      cache.setRule(getRule(key));
      return cache;
    }
 
  public:
 
    Cache& operator[](const QuadratureRuleKey& key) {
      auto it = elementCacheMap_.find(key);
      if (it == elementCacheMap_.end())
        return (elementCacheMap_.insert(std::make_pair(key, makeCache(key)))).first->second;
      else
        return it->second;
    }
 
    Cache& operator[](const FacetKey& facetKey) {
      auto it = facetCacheMap_.find(facetKey);
      if (it == facetCacheMap_.end())
        return (facetCacheMap_.insert(std::make_pair(facetKey, makeCache(facetKey)))).first->second;
      else
        return it->second;
    }
 
    Cache& prototype()
    {
      return prototype_;
    }
 
    void clear()
    {
      prototype_.clear();
      elementCacheMap_.clear();
      facetCacheMap_.clear();
      facetQuadratureRules_.clear();
    }
 
  private:
    std::map<ElementKey, Cache> elementCacheMap_;
    std::map<FacetKey, Cache> facetCacheMap_;
    Cache prototype_;
    std::list<FacetQuadratureRule> facetQuadratureRules_;
  };
 
 
 
  class UniqueCacheId : public std::pair<const void*, std::type_index>
  {
    using Base = std::pair<const void*, std::type_index>;
  public:
 
    template<class T>
    UniqueCacheId(const T& t) :
      Base(&t, std::type_index(typeid(T)))
    {}
  };
 
  template<class CT, int dimension>
  class CacheManager
  {
  public:
    using QuadratureRule = Dune::QuadratureRule<CT, dimension>;
    using size_type = std::size_t;
 
  private:
 
    // Simple type erasure for stored caches.
    class TypeErasedCache
    {
      struct {
        std::any impl_;
        void*(*toRawPtr_)(std::any&);
        void(*invalidate_)(void*);
        void(*setRule_)(void*, const QuadratureRule&);
      } members_;
      void* rawPtr_;
 
      void* toRawPtr()
      {
        return members_.toRawPtr_(members_.impl_);
      }
 
    public:
 
      template<class Impl>
      TypeErasedCache(Impl&& impl)
        : members_{
            std::move(impl),
            [](std::any& impl) -> void* { return &std::any_cast<Impl&>(impl); },
            [](void* impl) { static_cast<Impl*>(impl)->invalidate(); },
            [](void* impl, const QuadratureRule& rule) { static_cast<Impl*>(impl)->setRule(rule); },
          }
        , rawPtr_(toRawPtr())
      {}
 
      TypeErasedCache(TypeErasedCache&& other)
        : members_(std::move(other.members_))
        , rawPtr_(toRawPtr())
      {}
 
      TypeErasedCache(const TypeErasedCache& other)
        : members_(other.members_)
        , rawPtr_(toRawPtr())
      {}
 
      TypeErasedCache(TypeErasedCache& other)
        : members_(other.members_)
        , rawPtr_(toRawPtr())
      {}
 
      void invalidate()
      {
        members_.invalidate_(rawPtr_);
      }
 
      void setRule(const QuadratureRule& rule)
      {
        members_.setRule_(rawPtr_, rule);
      }
 
      template<class Impl>
      Impl& get()
      {
        return *static_cast<Impl*>(rawPtr_);
      }
    };
 
    const QuadratureRule* rule_ = nullptr;
    std::map<UniqueCacheId, size_type> cacheIndex_;
    std::vector<TypeErasedCache> caches_;
 
  public:
 
    CacheManager() = default;
    CacheManager(CacheManager&& other) = default;
    CacheManager(const CacheManager& other) = default;
 
    CacheManager& operator= (const CacheManager&) = delete;
    CacheManager& operator= (CacheManager&&) = delete;
 
    void clear()
    {
      rule_ = nullptr;
      cacheIndex_.clear();
      caches_.clear();
    }
 
    const QuadratureRule& rule() const
    {
      return *rule_;
    }
 
    void setRule(const QuadratureRule& rule)
    {
      rule_ = &rule;
      for(auto index : Dune::range(caches_.size()))
        caches_[index].setRule(rule);
    }
 
    void invalidate()
    {
      for(auto index : Dune::range(caches_.size()))
        caches_[index].invalidate();
    }
 
    // For diagnostics: Print number of stored caches
    void report() const
    {
      std::cout << "Number of managed caches " << caches_.size() << std::endl;
    }
 
    template<class Cache>
    size_type registerCache(UniqueCacheId uniqueCacheId, Cache&& cache)
    {
      auto it = cacheIndex_.find(uniqueCacheId);
      if (it == cacheIndex_.end())
      {
        caches_.push_back(std::forward<Cache>(cache));
        cacheIndex_[uniqueCacheId] = caches_.size()-1;
        return caches_.size()-1;
      }
      else
        return it->second;
    }
 
    template<class Cache>
    auto& getCache(size_type index)
    {
      return caches_[index].template get<Cache>();
    }
 
  };
 
 
 
  template<class CT, int dimension>
  using MultipleQuadratureCacheManager = MultipleQuadratureCache<CacheManager<CT, dimension>>;
 
 
 
  template<class CT, int dimension, class V>
  class SimpleCache
  {
  public:
 
    using QuadratureRule = Dune::QuadratureRule<CT, dimension>;
 
    using Value = V;
    using ValueCache = typename std::vector<Value>;
 
    SimpleCache(bool isNonAffine) :
      isNonAffine_(isNonAffine)
    {}
 
    SimpleCache(const SimpleCache& other) = default;
 
    void setRule(const QuadratureRule& rule)
    {
      rule_ = &rule;
      valueCache_.resize(rule_->size());
      isEmpty_ = true;
    }
 
    const QuadratureRule& rule() const
    {
      return *rule_;
    }
 
    void setNonAffine()
    {
      isNonAffine_ = true;
    }
 
    void invalidate()
    {
      if (isNonAffine_)
        isEmpty_ = true;
    }
 
    bool isEmpty()
    {
      return isEmpty_;
    }
 
    void setEmpty(bool isEmpty)
    {
      isEmpty_ = isEmpty;
    }
 
    auto& getValues()
    {
      return valueCache_;
    }
 
    const auto& getValues() const
    {
      return valueCache_;
    }
 
  private:
    const QuadratureRule* rule_ = nullptr;
    ValueCache valueCache_;
    bool isEmpty_ = true;
    bool isNonAffine_ = false;
  };
 
 
 
} // namespace Dune::Fufem::Forms
 
 
#endif // DUNE_FUFEM_FORMS_SHAPEFUNCTIONCACHE_HH