discreteglobalbasisfunction.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FUNCTIONS_GRIDFUNCTIONS_DISCRETEGLOBALBASISFUNCTIONS_HH
#define DUNE_FUNCTIONS_GRIDFUNCTIONS_DISCRETEGLOBALBASISFUNCTIONS_HH
 
#include <memory>
#include <optional>
 
#include <dune/common/typetraits.hh>
 
#include <dune/typetree/treecontainer.hh>
 
#include <dune/functions/functionspacebases/hierarchicnodetorangemap.hh>
#include <dune/functions/functionspacebases/flatvectorview.hh>
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <dune/functions/gridfunctions/gridfunction.hh>
#include <dune/functions/backends/concepts.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
 
namespace Dune {
namespace Functions {
 
 
 
template<typename B, typename V,
  typename NTRE = HierarchicNodeToRangeMap,
  typename R = typename V::value_type>
class DiscreteGlobalBasisFunction
{
public:
  using Basis = B;
  using Vector = V;
 
  // In order to make the cache work for proxy-references
  // we have to use AutonomousValue<T> instead of std::decay_t<T>
  using Coefficient = Dune::AutonomousValue<decltype(std::declval<Vector>()[std::declval<typename Basis::MultiIndex>()])>;
 
  using GridView = typename Basis::GridView;
  using EntitySet = GridViewEntitySet<GridView, 0>;
  using Tree = typename Basis::LocalView::Tree;
  using NodeToRangeEntry = NTRE;
 
  using Domain = typename EntitySet::GlobalCoordinate;
  using Range = R;
 
  using LocalDomain = typename EntitySet::LocalCoordinate;
  using Element = typename EntitySet::Element;
 
  using Traits = Imp::GridFunctionTraits<Range(Domain), EntitySet, DefaultDerivativeTraits, 16>;
 
  class LocalFunction
  {
    using LocalView = typename Basis::LocalView;
    using size_type = typename Tree::size_type;
 
    template<class Node>
    using LocalBasisRange = typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
 
    template<class Node>
    using NodeData = typename std::vector<LocalBasisRange<Node>>;
 
    using PerNodeEvaluationBuffer = typename TypeTree::TreeContainer<NodeData,Tree>;
 
  public:
 
    using GlobalFunction = DiscreteGlobalBasisFunction;
    using Domain = LocalDomain;
    using Range = GlobalFunction::Range;
    using Element = GlobalFunction::Element;
 
    LocalFunction(const DiscreteGlobalBasisFunction& globalFunction)
      : globalFunction_(&globalFunction)
      , localView_(globalFunction.basis().localView())
      , evaluationBuffer_(localView_.tree())
      , bound_(false)
    {
      localDoFs_.reserve(localView_.maxSize());
    }
 
    LocalFunction(const LocalFunction& other)
      : globalFunction_(other.globalFunction_)
      , localView_(other.localView_)
      , evaluationBuffer_(localView_.tree())
      , bound_(other.bound_)
    {
      localDoFs_.reserve(localView_.maxSize());
      if (bound())
        localDoFs_ = other.localDoFs_;
    }
 
    LocalFunction& operator=(const LocalFunction& other)
    {
      globalFunction_ = other.globalFunction_;
      localView_ = other.localView_;
      bound_ = other.bound_;
      if (bound())
        localDoFs_ = other.localDoFs_;
      return *this;
    }
 
    void bind(const Element& element)
    {
      localView_.bind(element);
      bound_ = true;
      // Use cache of full local view size. For a subspace basis,
      // this may be larger than the number of local DOFs in the
      // tree. In this case only cache entries associated to local
      // DOFs in the subspace are filled. Cache entries associated
      // to local DOFs which not contained in the subspace will not
      // be touched.
      //
      // Alternatively one could use a cache that exactly fits
      // the size of the tree. However, this would require to
      // subtract an offset from localIndex(i) on each cache
      // access in operator().
      localDoFs_.resize(localView_.size());
      for (size_type i = 0; i < localView_.tree().size(); ++i)
      {
        // For a subspace basis the index-within-tree i
        // is not the same as the localIndex withn the
        // full local view.
        size_t localIndex = localView_.tree().localIndex(i);
        localDoFs_[localIndex] = globalFunction_->dofs()[localView_.index(localIndex)];
      }
    }
 
    void unbind()
    {
      localView_.unbind();
      bound_ = false;
    }
 
    bool bound() const {
      return bound_;
    }
 
    Range operator()(const Domain& x) const
    {
      auto y = Range(0);
 
      TypeTree::forEachLeafNode(localView_.tree(), [&](auto&& node, auto&& treePath) {
        const auto& nodeToRangeEntry = globalFunction_->nodeToRangeEntry();
        const auto& fe = node.finiteElement();
        const auto& localBasis = fe.localBasis();
        auto& shapeFunctionValues = evaluationBuffer_[treePath];
 
        localBasis.evaluateFunction(x, shapeFunctionValues);
 
        // Get range entry associated to this node
        auto re = flatVectorView(nodeToRangeEntry(node, treePath, y));
 
        for (size_type i = 0; i < localBasis.size(); ++i)
        {
          // Get coefficient associated to i-th shape function
          auto c = flatVectorView(localDoFs_[node.localIndex(i)]);
 
          // Get value of i-th shape function
          auto v = flatVectorView(shapeFunctionValues[i]);
 
          // Notice that the range entry re, the coefficient c, and the shape functions
          // value v may all be scalar, vector, matrix, or general container valued.
          // The matching of their entries is done via the multistage procedure described
          // in the class documentation of DiscreteGlobalBasisFunction.
          auto&& dimC = c.size();
          auto dimV = v.size();
          assert(dimC*dimV == re.size());
          for(size_type j=0; j<dimC; ++j)
          {
            auto&& c_j = c[j];
            for(size_type k=0; k<dimV; ++k)
              re[j*dimV + k] += c_j*v[k];
          }
        }
      });
 
      return y;
    }
 
    const Element& localContext() const
    {
      return localView_.element();
    }
 
    friend typename Traits::LocalFunctionTraits::DerivativeInterface derivative(const LocalFunction& t)
    {
      DUNE_THROW(NotImplemented,"not implemented");
    }
 
  private:
 
    const DiscreteGlobalBasisFunction* globalFunction_;
    LocalView localView_;
    mutable PerNodeEvaluationBuffer evaluationBuffer_;
    bool bound_ = false;
    std::vector<Coefficient> localDoFs_;
  };
 
  template<class B_T, class V_T, class NTRE_T>
  DiscreteGlobalBasisFunction(B_T && basis, V_T && coefficients, NTRE_T&& nodeToRangeEntry) :
    entitySet_(basis.gridView()),
    basis_(wrap_or_move(std::forward<B_T>(basis))),
    coefficients_(wrap_or_move(std::forward<V_T>(coefficients))),
    nodeToRangeEntry_(wrap_or_move(std::forward<NTRE_T>(nodeToRangeEntry)))
  {}
 
  DiscreteGlobalBasisFunction(std::shared_ptr<const Basis> basis, std::shared_ptr<const V> coefficients, std::shared_ptr<const NodeToRangeEntry> nodeToRangeEntry) :
    entitySet_(basis->gridView()),
    basis_(basis),
    coefficients_(coefficients),
    nodeToRangeEntry_(nodeToRangeEntry)
  {}
 
  const Basis& basis() const
  {
    return *basis_;
  }
 
  const V& dofs() const
  {
    return *coefficients_;
  }
 
  const NodeToRangeEntry& nodeToRangeEntry() const
  {
    return *nodeToRangeEntry_;
  }
 
  // TODO: Implement this using hierarchic search
  Range operator() (const Domain& x) const
  {
    DUNE_THROW(NotImplemented,"not implemented");
  }
 
  friend typename Traits::DerivativeInterface derivative(const DiscreteGlobalBasisFunction& t)
  {
    DUNE_THROW(NotImplemented,"not implemented");
  }
 
  friend LocalFunction localFunction(const DiscreteGlobalBasisFunction& t)
  {
    return LocalFunction(t);
  }
 
  const EntitySet& entitySet() const
  {
    return entitySet_;
  }
 
private:
 
  EntitySet entitySet_;
  std::shared_ptr<const Basis> basis_;
  std::shared_ptr<const V> coefficients_;
  std::shared_ptr<const NodeToRangeEntry> nodeToRangeEntry_;
};
 
 
 
template<typename... TT>
void localFunction(DiscreteGlobalBasisFunction<TT...>&& t) = delete;
 
 
 
template<typename R, typename B, typename V>
auto makeDiscreteGlobalBasisFunction(B&& basis, V&& vector)
{
  using Basis = std::decay_t<B>;
  using NTREM = HierarchicNodeToRangeMap;
 
  // Small helper functions to wrap vectors using istlVectorBackend
  // if they do not already satisfy the VectorBackend interface.
  auto toConstVectorBackend = [&](auto&& v) -> decltype(auto) {
    if constexpr (models<Concept::ConstVectorBackend<Basis>, decltype(v)>()) {
      return std::forward<decltype(v)>(v);
    } else {
      return istlVectorBackend(v);
    }
  };
 
  using Vector = std::decay_t<decltype(toConstVectorBackend(std::forward<V>(vector)))>;
  return DiscreteGlobalBasisFunction<Basis, Vector, NTREM, R>(
      std::forward<B>(basis),
      toConstVectorBackend(std::forward<V>(vector)),
      HierarchicNodeToRangeMap());
}
 
 
 
} // namespace Functions
} // namespace Dune
 
#endif // DUNE_FUNCTIONS_GRIDFUNCTIONS_DISCRETEGLOBALBASISFUNCTIONS_HH