Add generalized derivative and local-to-global adapter to handle the global derivatives uniformly

Summary

Cleanup of Gridfunctions, i.e. separate order, derivative and gridfunction concepts. Add generalized derivative with gradient, divergence and partial derivatives. Add a local-to-global adapter similar to the one found in dune-localfunctions.

Details

1. Cleanup of Gridfunctions: The GridFunctionConcepts is split into several parts, Order.hpp, Derivative.hpp and GridFunction.hpp for a cleaner structure. Thereby, the free functions order(), derivative() and localFunctions() now redirect to the corresponding member functions when those member functions are available.

2. Generalized derivatives: The function derivative() is extended by an argument characterizing the type of the derivative, i.e. tag::gradient, tag::divergence, or tag::partial where the last one accepts an argument for the component of the partial derivative. In order to determine the range type of the grid functions, a class DerivativeTraits is added, parametrized with the derivative type and the function signature. For the discrete function it is possible to calculate all three derivative types. This is tested in GradienTest finally.

3. Local-to-Global adapter: The calculation of global derivatives from local-basis and geometry Jacobians is abstracted and put into the new class LocalToGlobalBasisAdapter that performes cached evaluation of the local basis jacobians and geometry jacobians and the transforms those into global derivatives. This class can be used in the localoperators and is partially adapted already.

src/amdis/gridfunctions/DiscreteFunction.inc.hpp

@@ -200,6 +200,17 @@ public:
     return localView_.element();
   }
 
+  Geometry const& geometry() const
+  {
+    assert( bound_ );
+    return geometry_.value();
+  }
+
+  LocalView const& localView() const
+  {
+    return localView_;
+  }
+
 protected:
   DiscreteFunction globalFunction_;
   Type type_;
@@ -219,30 +230,27 @@ public:
   using Range = typename Super::Range;
   using Domain = typename Super::Domain;
 
-  using Super::DerivativeLocalFunctionBase;
-  using Super::globalFunction_;
-  using Super::geometry_;
-  using Super::subTree_;
-  using Super::localView_;
+  // adopt constructor from base class
+  using DerivativeLocalFunctionBase<tag::gradient>::DerivativeLocalFunctionBase;
 
   /// Evaluate Gradient at bound element in local coordinates
   Range operator()(Domain const& x) const
   {
-    assert( Super::bound_ );
+    assert( this->bound_ );
     Range dy(0);
 
-    auto&& coefficients = *globalFunction_.dofVector_;
-    auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
-    for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
+    auto&& coefficients = *this->globalFunction_.dofVector_;
+    auto&& nodeToRangeEntry = this->globalFunction_.nodeToRangeEntry_;
+    for_each_leaf_node(*this->subTree_, [&,this](auto const& node, auto const& tp)
     {
-      auto localBasis = makeLocalToGlobalBasisAdapter(node, geometry_.value());
+      auto localBasis = makeLocalToGlobalBasisAdapter(node, this->geometry());
       auto const& gradients = localBasis.gradientsAt(x);
 
       // Get range entry associated to this node
       auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
 
       for (std::size_t i = 0; i < localBasis.size(); ++i) {
-        auto&& multiIndex = localView_.index(node.localIndex(i));
+        auto&& multiIndex = this->localView().index(node.localIndex(i));
 
         // Get coefficient associated to i-th shape function
         auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);
@@ -276,11 +284,8 @@ public:
   using Range = typename Super::Range;
   using Domain = typename Super::Domain;
 
-  using Super::DerivativeLocalFunctionBase;
-  using Super::globalFunction_;
-  using Super::geometry_;
-  using Super::subTree_;
-  using Super::localView_;
+  // adopt constructor from base class
+  using DerivativeLocalFunctionBase<tag::divergence>::DerivativeLocalFunctionBase;
 
   /// Evaluate divergence at bound element in local coordinates
   Range operator()(Domain const& x) const
@@ -299,13 +304,13 @@ private:
   {
     static_assert(Size_v<Range> == 1, "Implemented for scalar coefficients only.");
 
-    assert( Super::bound_ );
+    assert( this->bound_ );
     Range dy(0);
 
-    auto&& coefficients = *globalFunction_.dofVector_;
-    auto&& node = *subTree_;
+    auto&& coefficients = *this->globalFunction_.dofVector_;
+    auto&& node = *this->subTree_;
 
-    auto localBasis = makeLocalToGlobalBasisAdapter(node.child(0), geometry_.value());
+    auto localBasis = makeLocalToGlobalBasisAdapter(node.child(0), this->geometry());
     auto const& gradients = localBasis.gradientsAt(x);
 
     auto re = Dune::Functions::flatVectorView(dy);
@@ -315,7 +320,7 @@ private:
 
       assert(int(grad.size()) == GridView::dimensionworld);
       for (std::size_t j = 0; j < GridView::dimensionworld; ++j) {
-        auto&& multiIndex = localView_.index(node.child(j).localIndex(i));
+        auto&& multiIndex = this->localView().index(node.child(j).localIndex(i));
         re[0] += coefficients[multiIndex] * grad[j];
       }
     }
@@ -335,32 +340,28 @@ public:
   using Range = typename Super::Range;
   using Domain = typename Super::Domain;
 
-  using Super::DerivativeLocalFunctionBase;
-  using Super::globalFunction_;
-  using Super::geometry_;
-  using Super::subTree_;
-  using Super::localView_;
+  using DerivativeLocalFunctionBase<tag::partial>::DerivativeLocalFunctionBase;
 
   /// Evaluate partial derivative at bound element in local coordinates
   Range operator()(Domain const& x) const
   {
-    assert( Super::bound_ );
+    assert( this->bound_ );
     Range dy(0);
 
-    std::size_t comp = Super::type_.comp;
+    std::size_t comp = this->type_.comp;
 
-    auto&& coefficients = *globalFunction_.dofVector_;
-    auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
-    for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
+    auto&& coefficients = *this->globalFunction_.dofVector_;
+    auto&& nodeToRangeEntry = this->globalFunction_.nodeToRangeEntry_;
+    for_each_leaf_node(*this->subTree_, [&,this](auto const& node, auto const& tp)
     {
-      auto localBasis = makeLocalToGlobalBasisAdapter(node, geometry_.value());
+      auto localBasis = makeLocalToGlobalBasisAdapter(node, this->geometry());
       auto const& partial = localBasis.partialsAt(x, comp);
 
       // Get range entry associated to this node
       auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
 
       for (std::size_t i = 0; i < localBasis.size(); ++i) {
-        auto&& multiIndex = localView_.index(node.localIndex(i));
+        auto&& multiIndex = this->localView().index(node.localIndex(i));
 
         // Get coefficient associated to i-th shape function
         auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);