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/grid/utility/hierarchicsearch.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 {
 
 
namespace ImplDoc {
 
template<typename B, typename V, typename NTRE>
class DiscreteGlobalBasisFunctionBase
{
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 LocalDomain = typename EntitySet::LocalCoordinate;
  using Element = typename EntitySet::Element;
 
protected:
 
  // This collects all data that is shared by all related
  // global and local functions. This way we don't need to
  // keep track of it individually.
  struct Data
  {
    EntitySet entitySet;
    std::shared_ptr<const Basis> basis;
    std::shared_ptr<const Vector> coefficients;
    std::shared_ptr<const NodeToRangeEntry> nodeToRangeEntry;
  };
 
public:
  class LocalFunctionBase
  {
    using LocalView = typename Basis::LocalView;
    using size_type = typename Tree::size_type;
 
  public:
    using Domain = LocalDomain;
    using Element = typename EntitySet::Element;
 
  protected:
    LocalFunctionBase(const std::shared_ptr<const Data>& data)
      : data_(data)
      , localView_(data_->basis->localView())
    {
      localDoFs_.reserve(localView_.maxSize());
    }
 
    LocalFunctionBase(const LocalFunctionBase& other)
      : data_(other.data_)
      , localView_(other.localView_)
    {
      localDoFs_.reserve(localView_.maxSize());
      if (bound())
        localDoFs_ = other.localDoFs_;
    }
 
    LocalFunctionBase& operator=(const LocalFunctionBase& other)
    {
      data_ = other.data_;
      localView_ = other.localView_;
      if (bound())
        localDoFs_ = other.localDoFs_;
      return *this;
    }
 
  public:
    void bind(const Element& element)
    {
      localView_.bind(element);
      // 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 are 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());
      const auto& dofs = *data_->coefficients;
      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 within the
        // full local view.
        size_t localIndex = localView_.tree().localIndex(i);
        localDoFs_[localIndex] = dofs[localView_.index(localIndex)];
      }
    }
 
    void unbind()
    {
      localView_.unbind();
    }
 
    bool bound() const
    {
      return localView_.bound();
    }
 
    const Element& localContext() const
    {
      return localView_.element();
    }
 
  protected:
 
    template<class To, class From>
    void assignWith(To& to, const From& from) const
    {
        auto from_flat = flatVectorView(from);
        auto to_flat = flatVectorView(to);
        assert(from_flat.size() == to_flat.size());
        for (size_type i = 0; i < to_flat.size(); ++i)
          to_flat[i] = from_flat[i];
    }
 
    template<class Node, class TreePath, class Range>
    decltype(auto) nodeToRangeEntry(const Node& node, const TreePath& treePath, Range& y) const
    {
      return (*data_->nodeToRangeEntry)(node, treePath, y);
    }
 
    std::shared_ptr<const Data> data_;
    LocalView localView_;
    std::vector<Coefficient> localDoFs_;
  };
 
protected:
  DiscreteGlobalBasisFunctionBase(const std::shared_ptr<const Data>& data)
    : data_(data)
  {
    /* Nothing. */
  }
 
public:
 
  const Basis& basis() const
  {
    return *data_->basis;
  }
 
  const Vector& dofs() const
  {
    return *data_->coefficients;
  }
 
  const NodeToRangeEntry& nodeToRangeEntry() const
  {
    return *data_->nodeToRangeEntry;
  }
 
  const EntitySet& entitySet() const
  {
    return data_->entitySet;
  }
 
protected:
  std::shared_ptr<const Data> data_;
};
 
} // namespace ImplDoc
 
 
 
template<typename DGBF>
class DiscreteGlobalBasisFunctionDerivative;
 
template<typename B, typename V,
  typename NTRE = HierarchicNodeToRangeMap,
  typename R = typename V::value_type>
class DiscreteGlobalBasisFunction
  : public ImplDoc::DiscreteGlobalBasisFunctionBase<B, V, NTRE>
{
  using Base = ImplDoc::DiscreteGlobalBasisFunctionBase<B, V, NTRE>;
  using Data = typename Base::Data;
 
public:
  using Basis = typename Base::Basis;
  using Vector = typename Base::Vector;
 
  using Domain = typename Base::Domain;
  using Range = R;
 
  using Traits = Imp::GridFunctionTraits<Range(Domain), typename Base::EntitySet, DefaultDerivativeTraits, 16>;
 
private:
 
  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, typename Base::Tree>;
 
public:
  class LocalFunction
    : public Base::LocalFunctionBase
  {
    using LocalBase = typename Base::LocalFunctionBase;
    using size_type = typename Base::Tree::size_type;
    using LocalBase::nodeToRangeEntry;
 
  public:
 
    using GlobalFunction = DiscreteGlobalBasisFunction;
    using Domain = typename LocalBase::Domain;
    using Range = GlobalFunction::Range;
    using Element = typename LocalBase::Element;
 
    LocalFunction(const DiscreteGlobalBasisFunction& globalFunction)
      : LocalBase(globalFunction.data_)
      , evaluationBuffer_(this->localView_.tree())
    {
      /* Nothing. */
    }
 
    Range operator()(const Domain& x) const
    {
      Range y;
      istlVectorBackend(y) = 0;
 
      TypeTree::forEachLeafNode(this->localView_.tree(), [&](auto&& node, auto&& treePath) {
        const auto& fe = node.finiteElement();
        const auto& localBasis = fe.localBasis();
        auto& shapeFunctionValues = evaluationBuffer_[treePath];
 
        localBasis.evaluateFunction(x, shapeFunctionValues);
 
        // Compute linear combinations of basis function jacobian.
        // Non-scalar coefficients of dimension coeffDim are handled by
        // processing the coeffDim linear combinations independently
        // and storing them as entries of an array.
        using Value = LocalBasisRange< std::decay_t<decltype(node)> >;
        static constexpr auto coeffDim = decltype(flatVectorView(this->localDoFs_[node.localIndex(0)]).size())::value;
        auto values = std::array<Value, coeffDim>{};
        istlVectorBackend(values) = 0;
        for (size_type i = 0; i < localBasis.size(); ++i)
        {
          auto c = flatVectorView(this->localDoFs_[node.localIndex(i)]);
          for (std::size_t j = 0; j < coeffDim; ++j)
            values[j].axpy(c[j], shapeFunctionValues[i]);
        }
 
        // Assign computed values to node entry of range.
        // Types are matched using the lexicographic ordering provided by flatVectorView.
        LocalBase::assignWith(nodeToRangeEntry(node, treePath, y), values);
      });
 
      return y;
    }
 
    friend typename DiscreteGlobalBasisFunctionDerivative<DiscreteGlobalBasisFunction>::LocalFunction derivative(const LocalFunction& lf)
    {
      auto dlf = localFunction(DiscreteGlobalBasisFunctionDerivative<DiscreteGlobalBasisFunction>(lf.data_));
      if (lf.bound())
        dlf.bind(lf.localContext());
      return dlf;
    }
 
  private:
    mutable PerNodeEvaluationBuffer evaluationBuffer_;
  };
 
  template<class B_T, class V_T, class NTRE_T>
  DiscreteGlobalBasisFunction(B_T && basis, V_T && coefficients, NTRE_T&& nodeToRangeEntry)
    : Base(std::make_shared<Data>(Data{{basis.gridView()}, wrap_or_move(std::forward<B_T>(basis)), wrap_or_move(std::forward<V_T>(coefficients)), 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 typename Base::NodeToRangeEntry> nodeToRangeEntry)
    : Base(std::make_shared<Data>(Data{{basis->gridView()}, basis, coefficients, nodeToRangeEntry}))
  {}
 
  Range operator() (const Domain& x) const
  {
    HierarchicSearch search(this->data_->basis->gridView().grid(), this->data_->basis->gridView().indexSet());
 
    const auto e = search.findEntity(x);
    auto localThis = localFunction(*this);
    localThis.bind(e);
    return localThis(e.geometry().local(x));
  }
 
  friend DiscreteGlobalBasisFunctionDerivative<DiscreteGlobalBasisFunction> derivative(const DiscreteGlobalBasisFunction& f)
  {
    return DiscreteGlobalBasisFunctionDerivative<DiscreteGlobalBasisFunction>(f.data_);
  }
 
  friend LocalFunction localFunction(const DiscreteGlobalBasisFunction& t)
  {
    return LocalFunction(t);
  }
};
 
 
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());
}
 
 
template<typename DGBF>
class DiscreteGlobalBasisFunctionDerivative
  : public ImplDoc::DiscreteGlobalBasisFunctionBase<typename DGBF::Basis, typename DGBF::Vector, typename DGBF::NodeToRangeEntry>
{
  using Base = ImplDoc::DiscreteGlobalBasisFunctionBase<typename DGBF::Basis, typename DGBF::Vector, typename DGBF::NodeToRangeEntry>;
  using Data = typename Base::Data;
 
public:
  using DiscreteGlobalBasisFunction = DGBF;
 
  using Basis = typename Base::Basis;
  using Vector = typename Base::Vector;
 
  using Domain = typename Base::Domain;
  using Range = typename SignatureTraits<typename DiscreteGlobalBasisFunction::Traits::DerivativeInterface>::Range;
 
  using Traits = Imp::GridFunctionTraits<Range(Domain), typename Base::EntitySet, DefaultDerivativeTraits, 16>;
 
private:
 
  template<class Node>
  using LocalBasisRange = typename Node::FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
  template<class Node>
  using NodeData = typename std::vector< LocalBasisRange<Node> >;
  using PerNodeEvaluationBuffer = typename TypeTree::TreeContainer<NodeData, typename Base::Tree>;
 
public:
 
  class LocalFunction
    : public Base::LocalFunctionBase
  {
    using LocalBase = typename Base::LocalFunctionBase;
    using size_type = typename Base::Tree::size_type;
    using LocalBase::nodeToRangeEntry;
 
  public:
    using GlobalFunction = DiscreteGlobalBasisFunctionDerivative;
    using Domain = typename LocalBase::Domain;
    using Range = GlobalFunction::Range;
    using Element = typename LocalBase::Element;
 
    LocalFunction(const GlobalFunction& globalFunction)
      : LocalBase(globalFunction.data_)
      , evaluationBuffer_(this->localView_.tree())
    {
      /* Nothing. */
    }
 
    void bind(const Element& element)
    {
      LocalBase::bind(element);
      geometry_.emplace(element.geometry());
    }
 
    void unbind()
    {
      geometry_.reset();
      LocalBase::unbind();
    }
 
    Range operator()(const Domain& x) const
    {
      Range y;
      istlVectorBackend(y) = 0;
 
      const auto& jacobianInverse = geometry_->jacobianInverse(x);
 
      TypeTree::forEachLeafNode(this->localView_.tree(), [&](auto&& node, auto&& treePath) {
        const auto& fe = node.finiteElement();
        const auto& localBasis = fe.localBasis();
        auto& shapeFunctionJacobians = evaluationBuffer_[treePath];
 
        localBasis.evaluateJacobian(x, shapeFunctionJacobians);
 
        // Compute linear combinations of basis function jacobian.
        // Non-scalar coefficients of dimension coeffDim are handled by
        // processing the coeffDim linear combinations independently
        // and storing them as entries of an array.
        using RefJacobian = LocalBasisRange< std::decay_t<decltype(node)> >;
        static constexpr auto coeffDim = decltype(flatVectorView(this->localDoFs_[node.localIndex(0)]).size())::value;
        auto refJacobians = std::array<RefJacobian, coeffDim>{};
        istlVectorBackend(refJacobians) = 0;
        for (size_type i = 0; i < localBasis.size(); ++i)
        {
          auto c = flatVectorView(this->localDoFs_[node.localIndex(i)]);
          for (std::size_t j = 0; j < coeffDim; ++j)
            refJacobians[j].axpy(c[j], shapeFunctionJacobians[i]);
        }
 
        // Transform Jacobians form local to global coordinates.
        using Jacobian = decltype(refJacobians[0] * jacobianInverse);
        auto jacobians = std::array<Jacobian, coeffDim>{};
        std::transform(
          refJacobians.begin(), refJacobians.end(), jacobians.begin(),
          [&](const auto& refJacobian) { return refJacobian * jacobianInverse; });
 
        // Assign computed Jacobians to node entry of range.
        // Types are matched using the lexicographic ordering provided by flatVectorView.
        LocalBase::assignWith(nodeToRangeEntry(node, treePath, y), jacobians);
      });
 
      return y;
    }
 
    friend typename Traits::LocalFunctionTraits::DerivativeInterface derivative(const LocalFunction&)
    {
      DUNE_THROW(NotImplemented, "derivative of derivative is not implemented");
    }
 
  private:
    mutable PerNodeEvaluationBuffer evaluationBuffer_;
    std::optional<typename Element::Geometry> geometry_;
  };
 
  DiscreteGlobalBasisFunctionDerivative(const std::shared_ptr<const Data>& data)
    : Base(data)
  {
    /* Nothing. */
  }
 
  Range operator()(const Domain& x) const
  {
    HierarchicSearch search(this->data_->basis->gridView().grid(), this->data_->basis->gridView().indexSet());
 
    const auto e = search.findEntity(x);
    auto localThis = localFunction(*this);
    localThis.bind(e);
    return localThis(e.geometry().local(x));
  }
 
  friend typename Traits::DerivativeInterface derivative(const DiscreteGlobalBasisFunctionDerivative& f)
  {
    DUNE_THROW(NotImplemented, "derivative of derivative is not implemented");
  }
 
  friend LocalFunction localFunction(const DiscreteGlobalBasisFunctionDerivative& f)
  {
    return LocalFunction(f);
  }
};
 
 
} // namespace Functions
} // namespace Dune
 
#endif // DUNE_FUNCTIONS_GRIDFUNCTIONS_DISCRETEGLOBALBASISFUNCTIONS_HH