make AMDiS GridFunction differentiable w.r.t. dune-functions definitions. Make...

make AMDiS GridFunction differentiable w.r.t. dune-functions definitions. Make DiscreteLocalFunction copyable

amdis/BiLinearForm.hpp

@@ -84,8 +84,8 @@ namespace AMDiS
       : BiLinearForm(Dune::wrap_or_move(FWD(rowBasis)))
     {}
 
-    std::shared_ptr<RowBasis const> const& rowBasis() const { return rowBasis_; }
-    std::shared_ptr<ColBasis const> const& colBasis() const { return colBasis_; }
+    RowBasis const& rowBasis() const { return *rowBasis_; }
+    ColBasis const& colBasis() const { return *colBasis_; }
 
     /// \brief Associate a local operator with this BiLinearForm
     /**

amdis/BiLinearForm.inc.hpp

@@ -21,13 +21,13 @@ addOperator(ContextTag contextTag, Expr const& expr,
       "col must be a valid treepath, or an integer/index-constant");
 
   auto i = makeTreePath(row);
-  auto node_i = child(this->rowBasis()->localView().treeCache(), i);
+  auto node_i = child(this->rowBasis().localView().treeCache(), i);
   auto j = makeTreePath(col);
-  auto node_j = child(this->colBasis()->localView().treeCache(), j);
+  auto node_j = child(this->colBasis().localView().treeCache(), j);
 
   using LocalContext = typename ContextTag::type;
   using Tr = DefaultAssemblerTraits<LocalContext, ElementMatrix>;
-  auto op = makeLocalOperator<LocalContext>(expr, this->rowBasis()->gridView());
+  auto op = makeLocalOperator<LocalContext>(expr, this->rowBasis().gridView());
   auto localAssembler = makeUniquePtr(makeAssembler<Tr>(std::move(op), node_i, node_j));
 
   operators_[i][j].push(contextTag, std::move(localAssembler));
@@ -45,7 +45,7 @@ assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView)
   elementMatrix_.resize(rowLocalView.size(), colLocalView.size());
   elementMatrix_ = 0;
 
-  auto const& gv = this->rowBasis()->gridView();
+  auto const& gv = this->rowBasis().gridView();
   auto const& element = rowLocalView.element();
   auto geometry = element.geometry();
 
@@ -68,11 +68,11 @@ template <class RB, class CB, class T, class Traits>
 void BiLinearForm<RB,CB,T,Traits>::
 assemble()
 {
-  auto rowLocalView = this->rowBasis()->localView();
-  auto colLocalView = this->colBasis()->localView();
+  auto rowLocalView = this->rowBasis().localView();
+  auto colLocalView = this->colBasis().localView();
 
   this->init();
-  for (auto const& element : elements(this->rowBasis()->gridView(), typename Traits::PartitionSet{})) {
+  for (auto const& element : elements(this->rowBasis().gridView(), typename Traits::PartitionSet{})) {
     rowLocalView.bind(element);
     if (this->rowBasis() == this->colBasis())
       this->assemble(rowLocalView, rowLocalView);

amdis/DOFVector.hpp

@@ -82,9 +82,9 @@ namespace AMDiS
     {}
 
     /// Return the global basis
-    std::shared_ptr<GlobalBasis const> const& basis() const
+    GlobalBasis const& basis() const
     {
-      return basis_;
+      return *basis_;
     }
 
     Coefficients const& coefficients() const

amdis/DOFVector.inc.hpp

@@ -16,12 +16,12 @@ backup(std::string const& filename)
 {
   std::ofstream out(filename, std::ios::binary);
 
-  std::int64_t numElements = this->basis()->gridView().size(0);
+  std::int64_t numElements = this->basis().gridView().size(0);
   out.write((char*)&numElements, sizeof(std::int64_t));
 
-  auto localView = this->basis()->localView();
+  auto localView = this->basis().localView();
   std::vector<value_type> data;
-  for (auto const& element : elements(this->basis()->gridView()))
+  for (auto const& element : elements(this->basis().gridView()))
   {
     localView.bind(element);
     this->gather(localView, data);
@@ -43,13 +43,13 @@ restore(std::string const& filename)
 
   std::int64_t numElements = 0;
   in.read((char*)&numElements, sizeof(std::int64_t));
-  assert(numElements == this->basis()->gridView().size(0));
+  assert(numElements == this->basis().gridView().size(0));
 
   // assume the order of element traversal is not changed
-  auto localView = this->basis()->localView();
+  auto localView = this->basis().localView();
   std::vector<value_type> data;
-  this->init(sizeInfo(*this->basis()), true);
-  for (auto const& element : elements(this->basis()->gridView()))
+  this->init(sizeInfo(this->basis()), true);
+  for (auto const& element : elements(this->basis().gridView()))
   {
     std::uint64_t len = 0;
     in.read((char*)&len, sizeof(std::uint64_t));

amdis/LinearForm.hpp

@@ -55,7 +55,7 @@ namespace AMDiS
       : LinearForm(Dune::wrap_or_move(FWD(basis)))
     {}
 
-    std::shared_ptr<GlobalBasis const> const& basis() const { return basis_; }
+    GlobalBasis const& basis() const { return *basis_; }
 
     /// \brief Associate a local operator with this LinearForm
     /**

amdis/LinearForm.inc.hpp

@@ -20,11 +20,11 @@ addOperator(ContextTag contextTag, Expr const& expr, TreePath path)
       "path must be a valid treepath, or an integer/index-constant");
 
   auto i = makeTreePath(path);
-  auto node = child(this->basis()->localView().treeCache(), i);
+  auto node = child(this->basis().localView().treeCache(), i);
 
   using LocalContext = typename ContextTag::type;
   using Tr = DefaultAssemblerTraits<LocalContext, ElementVector>;
-  auto op = makeLocalOperator<LocalContext>(expr, this->basis()->gridView());
+  auto op = makeLocalOperator<LocalContext>(expr, this->basis().gridView());
   auto localAssembler = makeUniquePtr(makeAssembler<Tr>(std::move(op), node));
 
   operators_[i].push(contextTag, std::move(localAssembler));
@@ -39,7 +39,7 @@ assemble(LocalView const& localView)
   elementVector_.resize(localView.size());
   elementVector_ = 0;
 
-  auto const& gv = this->basis()->gridView();
+  auto const& gv = this->basis().gridView();
   auto const& element = localView.element();
   auto geometry = element.geometry();
 
@@ -60,10 +60,10 @@ template <class GB, class T, class Traits>
 void LinearForm<GB,T,Traits>::
 assemble()
 {
-  auto localView = this->basis()->localView();
+  auto localView = this->basis().localView();
 
   this->init(sizeInfo(*this->basis()), true);
-  for (auto const& element : elements(this->basis()->gridView(), typename Traits::PartitionSet{})) {
+  for (auto const& element : elements(this->basis().gridView(), typename Traits::PartitionSet{})) {
     localView.bind(element);
     this->assemble(localView);
     localView.unbind();

amdis/common/DerivativeTraits.hpp

@@ -15,6 +15,7 @@ namespace AMDiS
 {
   namespace tag
   {
+    struct value {};
     struct jacobian {};
     struct gradient {};
     struct divergence {};
@@ -27,6 +28,12 @@ namespace AMDiS
   template <class Sig, class Type>
   struct DerivativeTraits;
 
+  template <class R, class D>
+  struct DerivativeTraits<R(D), tag::value>
+  {
+    using Range = R;
+  };
+
   template <class R, class D>
   struct DerivativeTraits<R(D), tag::jacobian>
       : public Dune::Functions::DefaultDerivativeTraits<R(D)>

amdis/functions/GlobalIdSet.hpp

@@ -124,11 +124,11 @@ namespace AMDiS
 
   public:
     GlobalBasisIdSet(GlobalBasis const& globalBasis)
-      : tree_(makeNode(globalBasis.preBasis(), TreePath{}))
+      : tree_(globalBasis.localView().tree())
       , nodeIdSet_(globalBasis.gridView())
       , twist_(globalBasis.gridView().grid().globalIdSet())
     {
-      Dune::Functions::initializeTree(tree_);
+      initializeTree(tree_);
     }
 
     /// \brief Bind the IdSet to a grid element.
@@ -139,7 +139,7 @@ namespace AMDiS
      **/
     void bind(Element const& element)
     {
-      Dune::Functions::bindTree(tree_, element);
+      bindTree(tree_, element);
       nodeIdSet_.bind(tree_);
       twist_.bind(element);
       data_.resize(size());

amdis/functions/LocalView.hpp

@@ -7,6 +7,7 @@
 #include <dune/common/concept.hh>
 #include <dune/functions/functionspacebases/concepts.hh>
 
+#include <amdis/common/TypeTraits.hpp>
 #include <amdis/functions/NodeCache.hpp>
 #include <amdis/functions/Nodes.hpp>
 
@@ -46,8 +47,8 @@ namespace AMDiS
     using MultiIndex = typename NodeIndexSet::MultiIndex;
 
   private:
-    template <class NIS>
-    using hasIndices = decltype(std::declval<NIS>().indices(std::declval<std::vector<typename NIS::MultiIndex>>().begin()));
+    template <class NIS, class Iter>
+    using hasIndices = decltype(std::declval<NIS>().indices(std::declval<Iter>()));
 
   public:
     /// \brief Construct local view for a given global finite element basis
@@ -61,6 +62,30 @@ namespace AMDiS
       initializeTree(tree_);
     }
 
+    /// Copy constructor.
+    // NOTE: needs to be implemented manually, because of the reference dependency
+    // between members tree_ and treeCache_
+    LocalView (LocalView const& other)
+      : globalBasis_(other.globalBasis_)
+      , tree_(other.tree_)
+      , treeCache_(makeNodeCache(tree_))
+      , nodeIndexSet_(other.nodeIndexSet_)
+      , element_(other.element_)
+      , indices_(other.indices_)
+    {}
+
+    // Move constructor
+    // NOTE: needs to be implemented manually, because of the reference dependency
+    // between members tree_ and treeCache_
+    LocalView (LocalView&& other)
+      : globalBasis_(std::move(other.globalBasis_))
+      , tree_(std::move(other.tree_))
+      , treeCache_(makeNodeCache(tree_))
+      , nodeIndexSet_(std::move(other.nodeIndexSet_))
+      , element_(std::move(other.element_))
+      , indices_(std::move(other.indices_))
+    {}
+
     /// \brief Bind the view to a grid element
     /**
      * Having to bind the view to an element before being able to actually access any of
@@ -74,7 +99,7 @@ namespace AMDiS
       nodeIndexSet_.bind(tree_);
       indices_.resize(size());
 
-      if constexpr (Dune::Std::is_detected<hasIndices,NodeIndexSet>{})
+      if constexpr (Dune::Std::is_detected_v<hasIndices, NodeIndexSet, TYPEOF(indices_.begin())>)
         nodeIndexSet_.indices(indices_.begin());
       else
         for (size_type i = 0; i < size(); ++i)
@@ -151,10 +176,11 @@ namespace AMDiS
 
   protected:
     GlobalBasis const* globalBasis_;
-    std::optional<Element> element_;
     Tree tree_;
     TreeCache treeCache_;
     NodeIndexSet nodeIndexSet_;
+
+    std::optional<Element> element_;
     std::vector<MultiIndex> indices_;
   };
 

amdis/gridfunctions/Derivative.hpp

@@ -17,6 +17,16 @@ namespace AMDiS
     return lf.makeDerivative(type);
   }
 
+  /// Implementation of local-function derivative interface of dune-functions.
+  /// Implements the jacobian derivative type.
+  template <class LocalFunction>
+  auto derivative(LocalFunction const& lf)
+    -> decltype(lf.makeDerivative(tag::gradient{}))
+  {
+    return lf.makeDerivative(tag::gradient{});
+  }
+
+
   namespace Concepts
   {
     /** \addtogroup Concepts

amdis/gridfunctions/DiscreteFunction.hpp

@@ -48,9 +48,9 @@ namespace AMDiS
   public:
     /// Constructor. Stores a pointer to the mutable `dofvector`.
     template <class... Path>
-    DiscreteFunction(Coefficients& dofVector, GlobalBasis const& basis, Path... path)
-      : Super(dofVector, basis, path...)
-      , mutableCoeff_(&dofVector)
+    DiscreteFunction(Coefficients& coefficients, GlobalBasis const& basis, Path... path)
+      : Super(coefficients, basis, path...)
+      , mutableCoeff_(&coefficients)
     {}
 
     template <class... Path>
@@ -145,6 +145,7 @@ namespace AMDiS
   private:
     using Coefficients = std::remove_const_t<Coeff>;
     using GlobalBasis = GB;
+    using LocalView = typename GlobalBasis::LocalView;
     using GridView = typename GlobalBasis::GridView;
     using ValueType = typename Coefficients::value_type;
 
@@ -154,8 +155,6 @@ namespace AMDiS
     using TreeCache = NodeCache_t<Tree>;
     using SubTreeCache = typename Dune::TypeTree::ChildForTreePath<TreeCache, TreePath>;
 
-    using NodeToRangeEntry = Dune::Functions::DefaultNodeToRangeMap<SubTree>;
-
   public:
     /// Set of entities the DiscreteFunction is defined on
     using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
@@ -171,15 +170,8 @@ namespace AMDiS
     enum { hasDerivative = false };
 
   public:
-    /// A LocalFunction representing the derivative of the DOFVector on a bound element
+    /// A LocalFunction representing the localfunction and derivative of the DOFVector on a bound element
     template <class Type>
-    class DerivativeLocalFunctionBase;
-
-    class GradientLocalFunction;
-    class PartialLocalFunction;
-    class DivergenceLocalFunction;
-
-    /// A LocalFunction representing the value the DOFVector on a bound element
     class LocalFunction;
 
   public:
@@ -190,7 +182,6 @@ namespace AMDiS
       , basis_(&basis)
       , treePath_(makeTreePath(path...))
       , entitySet_(basis_->gridView())
-      , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(*basis_, treePath_))
     {}
 
     template <class... Path>
@@ -206,14 +197,13 @@ namespace AMDiS
       : DiscreteFunction(dofVector.coefficients(), dofVector.basis(), path...)
     {}
 
-
     /// \brief Evaluate DiscreteFunction in global coordinates. NOTE: expensive
     Range operator()(Domain const& x) const;
 
     /// \brief Create a local function for this view on the DOFVector. \relates LocalFunction
-    LocalFunction makeLocalFunction() const
+    LocalFunction<tag::value> makeLocalFunction() const
     {
-      return LocalFunction{*this};
+      return {std::make_shared<LocalView>(basis().localView()), treePath(), coefficients(), tag::value{}};
     }
 
     /// \brief Return a \ref Dune::Functions::GridViewEntitySet
@@ -252,9 +242,9 @@ namespace AMDiS
     GlobalBasis const* basis_;
     TreePath treePath_;
     EntitySet entitySet_;
-    NodeToRangeEntry nodeToRangeEntry_;
   };
 
+
   // deduction guides
 
   template <class Coeff, class Basis, class... Path,

amdis/gridfunctions/DiscreteLocalFunction.inc.hpp

@@ -36,215 +36,136 @@ void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoe
 
 
 template <class Coeff, class GB, class TP>
+  template <class Type>
 class DiscreteFunction<Coeff const,GB,TP>::LocalFunction
 {
+  using DomainType = typename DiscreteFunction::Domain;
+  using RangeType = typename DiscreteFunction::Range;
+
+  template <class T>
+  using DerivativeRange = typename DerivativeTraits<RangeType(DomainType), T>::Range;
+
 public:
   using Domain = typename EntitySet::LocalCoordinate;
-  using Range = typename DiscreteFunction::Range;
+  using Range = DerivativeRange<Type>;
 
-  enum { hasDerivative = true };
+  enum { hasDerivative = std::is_same<Type, tag::value>::value };
 
 private:
   using LocalView = typename GlobalBasis::LocalView;
   using Element = typename EntitySet::Element;
   using Geometry = typename Element::Geometry;
+  using NodeToRangeMap = Dune::Functions::DefaultNodeToRangeMap<SubTree>;
 
 public:
   /// Constructor. Stores a copy of the DiscreteFunction.
-  LocalFunction(DiscreteFunction const& globalFunction)
-    : globalFunction_(globalFunction)
-    , localView_(globalFunction_.basis().localView())
-    , treeCache_(makeNodeCache(localView_.tree()))
-    , subTreeCache_(&Dune::TypeTree::child(treeCache_, globalFunction_.treePath()))
-  {}
-
-  /// Copy constructor.
-  LocalFunction(LocalFunction const& other)
-    : globalFunction_(other.globalFunction_)
-    , localView_(globalFunction_.basis().localView())
-    , treeCache_(makeNodeCache(localView_.tree()))
-    , subTreeCache_(&Dune::TypeTree::child(treeCache_, globalFunction_.treePath()))
+  template <class LV>
+  LocalFunction(std::shared_ptr<LV> localView, TP const& treePath, Coeff const& coefficients, Type type)
+    : localView_(std::move(localView))
+    , treePath_(treePath)
+    , coefficients_(coefficients)
+    , type_(type)
+    , nodeToRangeMap_(subTree())
   {}
 
   /// \brief Bind the LocalView to the element
   void bind(Element const& element)
   {
-    localView_.bind(element);
-    Impl::gather(globalFunction_.coefficients(), localView_, localCoefficients_, Dune::PriorityTag<4>{});
+    localView_->bind(element);
+    geometry_.emplace(element.geometry());
+    Impl::gather(coefficients_, *localView_, localCoefficients_, Dune::PriorityTag<4>{});
     bound_ = true;
   }
 
   /// \brief Unbind the LocalView from the element
   void unbind()
   {
-    localView_.unbind();
+    localView_->unbind();
+    geometry_.reset();
     bound_ = false;
   }
 
   /// \brief Evaluate LocalFunction at bound element in local coordinates
-  Range operator()(Domain const& x) const
-  {
-    assert( bound_ );
-    Range y(0);
-
-    auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
-    for_each_leaf_node(*subTreeCache_, [&,this](auto const& node, auto const& tp)
-    {
-      auto const& shapeFunctionValues = node.localBasisValuesAt(x);
-      std::size_t size = node.finiteElement().size();
-
-      // Get range entry associated to this node
-      auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y));
-
-      for (std::size_t i = 0; i < size; ++i) {
-        // Get coefficient associated to i-th shape function
-        auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);
-
-        // Get value of i-th shape function
-        auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);
-
-        std::size_t dimC = c.size();
-        std::size_t dimV = v.size();
-        assert(dimC*dimV == std::size_t(re.size()));
-        for(std::size_t j = 0; j < dimC; ++j) {
-          auto&& c_j = c[j];
-          for(std::size_t k = 0; k < dimV; ++k)
-            re[j*dimV + k] += c_j*v[k];
-        }
-      }
-    });
-
-    return y;
-  }
-
-  /// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction
-  GradientLocalFunction makeDerivative(tag::gradient type) const
+  Range operator() (const Domain& local) const
   {
-    return GradientLocalFunction{globalFunction_, type};
-  }
+    assert(bound_);
+    if constexpr (std::is_same_v<Type, tag::value>)
+      return evaluateFunction(local);
+    else if constexpr (std::is_same_v<Type, tag::gradient>)
+      return evaluateJacobian(local);
+    else if constexpr (std::is_same_v<Type, tag::partial>)
+      return evaluatePartialDerivative(local);
+    else if constexpr (std::is_same_v<Type, tag::divergence>) {
+      static constexpr bool isVectorNode
+        = SubTree::isPower && SubTree::degree() == GridView::dimensionworld;
+      static_assert(isVectorNode);
+      if (isVectorNode)
+        return evaluateDivergence(local);
+    }
 
-  DivergenceLocalFunction makeDerivative(tag::divergence type) const
-  {
-    return DivergenceLocalFunction{globalFunction_, type};
+    return Range(0);
   }
 
-  PartialLocalFunction makeDerivative(tag::partial type) const
+  /// \brief Create a LocalFunction representing the derivative.
+  template <class DerType>
+  LocalFunction<DerType> makeDerivative(DerType type) const
   {
-    return PartialLocalFunction{globalFunction_, type};
+    static_assert(std::is_same_v<Type, tag::value>);
+    return {localView_, treePath_, coefficients_, type};
   }
 
   /// \brief The polynomial degree of the LocalFunctions basis-functions
   auto order() const
     -> decltype(AMDiS::order(std::declval<SubTree>()))
   {
-    assert( bound_ );
-    return AMDiS::order(*subTreeCache_);
+    assert(bound_);
+    return AMDiS::order(subTreeCache());
   }
 
   /// \brief Return the bound element
   Element const& localContext() const
   {
-    assert( bound_ );
-    return localView_.element();
+    assert(bound_);
+    return localView_->element();
   }
 
 private:
-  DiscreteFunction globalFunction_;
-  LocalView localView_;
-  TreeCache treeCache_;
-  SubTreeCache const* subTreeCache_;
-
-  std::vector<ValueType> localCoefficients_;
-  bool bound_ = false;
-};
-
-
-template <class Coeff, class GB, class TP>
-  template <class Type>
-class DiscreteFunction<Coeff const,GB,TP>::DerivativeLocalFunctionBase
-{
-  using R = typename DiscreteFunction::Range;
-  using D = typename DiscreteFunction::Domain;
-  using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
-  using Traits = DerivativeTraits<RawSignature, Type>;
-
-public:
-  using Domain = typename EntitySet::LocalCoordinate;
-  using Range = typename Traits::Range;
-
-  enum { hasDerivative = false };
-
-private:
-  using LocalView = typename GlobalBasis::LocalView;
-  using Element = typename EntitySet::Element;
-  using Geometry = typename Element::Geometry;
-
-public:
-  /// Constructor. Stores a copy of the DiscreteFunction.
-  DerivativeLocalFunctionBase(DiscreteFunction const& globalFunction, Type const& type)
-    : globalFunction_(globalFunction)
-    , type_(type)
-    , localView_(globalFunction_.basis().localView())
-    , treeCache_(makeNodeCache(localView_.tree()))
-    , subTreeCache_(&Dune::TypeTree::child(treeCache_, globalFunction_.treePath()))
-  {}
 
-  /// Copy constructor.
-  DerivativeLocalFunctionBase(DerivativeLocalFunctionBase const& other)
-    : globalFunction_(other.globalFunction_)
-    , type_(other.type_)
-    , localView_(globalFunction_.basis().localView())
-    , treeCache_(makeNodeCache(localView_.tree()))
-    , subTreeCache_(&Dune::TypeTree::child(treeCache_, globalFunction_.treePath()))
-  {}
+  template <class LocalCoordinate>
+  auto evaluateFunction (const LocalCoordinate& local) const;
 
-  void bind(Element const& element)
-  {
-    localView_.bind(element);
-    geometry_.emplace(element.geometry());
+  template <class LocalCoordinate>
+  auto evaluateJacobian (const LocalCoordinate& local) const;
 
-    Impl::gather(globalFunction_.coefficients(), localView_, localCoefficients_, Dune::PriorityTag<4>{});
-    bound_ = true;
-  }
-
-  void unbind()
-  {
-    localView_.unbind();
-    geometry_.reset();
-    bound_ = false;
-  }
+  template <class LocalCoordinate>
+  auto evaluateDivergence (const LocalCoordinate& local) const;
 
-  auto order() const
-    -> decltype(AMDiS::order(std::declval<SubTree>()))
-  {