derivative of gridFunction and functor applied to gridFunctions added

dune/amdis/DOFVectorView.hpp

 #pragma once
 
+#include <vector>
+
+#include <dune/common/std/optional.hh>
+#include <dune/functions/common/defaultderivativetraits.hh>
+#include <dune/functions/functionspacebases/defaultnodetorangemap.hh>
+#include <dune/functions/functionspacebases/flatvectorbackend.hh>
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+#include <dune/typetree/childextraction.hh>
+
+#include <dune/amdis/terms/GridViewExpression.hpp>
+#include <dune/amdis/terms/TermGenerator.hpp>
+
 namespace AMDiS
 {
+  template <class GlobalBasisType, class TreePathType, bool isConst = true>
+  class DOFVectorView;
+
   template <class GlobalBasisType, class TreePathType>
-  class DOFVectorView
+  class DOFVectorView<GlobalBasisType, TreePathType, true>
   {
   public:
     using GlobalBasis = GlobalBasisType;
@@ -11,25 +26,103 @@ namespace AMDiS
     using Vector = DOFVector<GlobalBasis>;
 
     using Tree = typename GlobalBasis::LocalView::Tree;
-    using SubTree = typename TypeTree::ChildForTreePath<Tree, TreePath>;
-    using NodeToRangeEntry = Dune::Functions::DefaultNodeToRangeMap<SubTree>
+    using SubTree = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
+    using NodeToRangeEntry = Dune::Functions::DefaultNodeToRangeMap<SubTree>;
 
     using GridView = typename GlobalBasis::GridView;
-    using EntitySet = GridViewEntitySet<GridView, 0>;
+    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
 
     using Domain = typename EntitySet::GlobalCoordinate;
     using Range = typename Vector::value_type;
 
+    using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<Range(Domain)>;
+    using DerivativeRange = typename DerivativeTraits::Range;
+
     using LocalDomain = typename EntitySet::LocalCoordinate;
     using Element = typename EntitySet::Element;
+    using Geometry = typename Element::Geometry;
+
+    template <class Block>
+    using Flat = Dune::Functions::FlatVectorBackend<Block>;
+
+  public: // a local view on the gradients
+
+    class GradientLocalFunction
+    {
+    public:
+      using Domain = LocalDomain;
+      using Range = DerivativeRange;
+
+    private:
+      using LocalBasisView = typename GlobalBasis::LocalView;
+      using LocalIndexSet = typename GlobalBasis::LocalIndexSet;
+
+      template <class LeafNode>
+      using LocalBasisJacobian = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
+
+      template <class Node>
+      using NodeData = typename std::vector<LocalBasisJacobian<Node>>;
+
+      using ReferenceGradientContainer = Dune::Functions::TreeData<SubTree, NodeData, true>;
+
+    public:
+      GradientLocalFunction(DOFVectorView const& globalFunction)
+        : globalFunction_(&globalFunction)
+        , localBasisView_(globalFunction_->basis().localView())
+        , localIndexSet_(globalFunction_->basis().localIndexSet())
+        , subTree_(&Dune::TypeTree::child(localBasisView_.tree(), globalFunction_->treePath()))
+      {
+        referenceGradientContainer_.init(*subTree_);
+      }
+
+      void bind(Element const& element)
+      {
+        localBasisView_.bind(element);
+        localIndexSet_.bind(localBasisView_);
+        geometry_.emplace(element.geometry());
+        bound_ = true;
+      }
+
+      void unbind()
+      {
+        localIndexSet_.unbind();
+        localBasisView_.unbind();
+        geometry_.reset();
+        bound_ = false;
+      }
+
+      /// Evaluate Gradient at bound element in local coordinates
+      Range operator()(Domain const& x) const;
+
+      /// Return the bound element
+      Element const& localContext() const
+      {
+        assert( bound_ );
+        return localBasisView_.element();
+      }
+
+    private:
+      DOFVectorView const* globalFunction_;
+      LocalBasisView localBasisView_;
+      LocalIndexSet localIndexSet_;
+
+      SubTree const* subTree_;
+      mutable ReferenceGradientContainer referenceGradientContainer_;
+
+      Dune::Std::optional<Geometry> geometry_;
+      bool bound_ = false;
+    };
+
+
+  public: // a local view on the values
 
     class LocalFunction
     {
     public:
       using Domain = typename DOFVectorView::LocalDomain;
       using Range = typename DOFVectorView::Range;
-      using Element = typename DOFVectorView::Element;
 
+    private:
       using LocalBasisView = typename GlobalBasis::LocalView;
       using LocalIndexSet = typename GlobalBasis::LocalIndexSet;
 
@@ -39,9 +132,6 @@ namespace AMDiS
       template <class Node>
       using NodeData = typename std::vector<LocalBasisRange<Node>>;
 
-      template <class Block>
-      using Flat = Dune::Functions::FlatVectorBackend<Block>;
-
       using ShapeFunctionValueContainer = Dune::Functions::TreeData<SubTree, NodeData, true>;
 
     public:
@@ -49,9 +139,9 @@ namespace AMDiS
         : globalFunction_(&globalFunction)
         , localBasisView_(globalFunction_->basis().localView())
         , localIndexSet_(globalFunction_->basis().localIndexSet())
-        , subTree_(TypeTree::child(localBasisView_.tree(), globalFunction_->treePath())
+        , subTree_(&Dune::TypeTree::child(localBasisView_.tree(), globalFunction_->treePath()))
       {
-        shapeFunctionValueContainer_.init(subTree_);
+        shapeFunctionValueContainer_.init(*subTree_);
       }
 
       void bind(Element const& element)
@@ -68,19 +158,19 @@ namespace AMDiS
         bound_ = false;
       }
 
-      /**
-      * \brief Evaluate LocalFunction at bound element.
-      *
-      * The result of this method is undefined if you did
-      * not call bind() beforehand or changed the coefficient
-      * vector after the last call to bind(). In the latter case
-      * you have to call bind() again in order to make operator()
-      * usable.
-      */
+      /// Evaluate LocalFunction at bound element in local coordinates
       Range operator()(Domain const& x) const;
 
+      /// Create a LocalFunction representing the gradient
+      friend GradientLocalFunction derivative(LocalFunction const& localFunction)
+      {
+        return GradientLocalFunction{*localFunction.globalFunction_};
+      }
+
+      /// Return the bound element
       Element const& localContext() const
       {
+        assert( bound_ );
         return localBasisView_.element();
       }
 
@@ -88,32 +178,41 @@ namespace AMDiS
       DOFVectorView const* globalFunction_;
       LocalBasisView localBasisView_;
       LocalIndexSet localIndexSet_;
+      SubTree const* subTree_;
 
       mutable ShapeFunctionValueContainer shapeFunctionValueContainer_;
-      SubTree const& subTree_;
 
       bool bound_ = false;
     };
 
+
   public:
 
-    DOFVectorView(DOFVector<GlobalBasis>& dofVector, TreePath const& treePath)
-      : basis_(basis)
+    /// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
+    DOFVectorView(DOFVector<GlobalBasis> const& dofVector, TreePath const& treePath)
+      : dofVector_(&dofVector)
       , treePath_(treePath)
-      , coefficients_(coefficients)
-      , entitySet_(basis.gridView())
-      , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(basis, treePath))
+      , entitySet_(dofVector.getFeSpace().gridView())
+      , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(dofVector.getFeSpace(), treePath))
     {}
 
+    /// Evaluate the view on this DOFVector in global coordinates
     Range operator()(Domain const& x) const
     {
       error_exit("Not implemented.");
       return Range(0);
     }
 
+    /// Create a local function for this view on the DOFVector
     friend LocalFunction localFunction(DOFVectorView const& self)
     {
-      return LocalExpr{self};
+      return LocalFunction{self};
+    }
+
+    /// Create a local function implementing the gradient of the DOFVector
+    friend GradientLocalFunction derivative(DOFVectorView const& self)
+    {
+      return GradientLocalFunction{self};
     }
 
     EntitySet const& entitySet() const
@@ -124,38 +223,134 @@ namespace AMDiS
 
   public:
 
+    /// Return global basis
+    GlobalBasis const& basis() const
+    {
+      return dofVector_->getFeSpace();
+    }
+
+    /// Return treePath associated with this view
+    TreePath const& treePath() const
+    {
+      return treePath_;
+    }
+
+    /// Return const coefficient vector
+    DOFVector<GlobalBasis> const& coefficients() const
+    {
+      return *dofVector_;
+    }
+
+
+  protected:
+
+    DOFVector<GlobalBasis> const* dofVector_;
+    TreePath const treePath_;
+
+    EntitySet entitySet_;
+    NodeToRangeEntry nodeToRangeEntry_;
+  };
+
+
+  // A mutable version of DOFVectorView
+  template <class GlobalBasisType, class TreePathType>
+  class DOFVectorView<GlobalBasisType, TreePathType, false>
+      : public DOFVectorView<GlobalBasisType, TreePathType, true>
+  {
+    using DOFVectorConstView = DOFVectorView<GlobalBasisType, TreePathType, true>;
+
+    using GlobalBasis = GlobalBasisType;
+    using TreePath = TreePathType;
+
+  public:
+
+    /// Constructor. Stores a pointer to the mutable `dofvector`.
+    DOFVectorView(DOFVector<GlobalBasis>& dofVector, TreePath const& treePath)
+      : DOFVectorConstView(dofVector, treePath)
+      , mutableDofVector_(&dofVector)
+    {}
+
+
+  public:
+
+    /// Interpolation of expr to DOFVector
     template <class PreExpr>
-    DOFVectorView& operator<<(PreExpr const& preExpr)
+    DOFVectorView& operator<<(PreExpr&& preExpr)
     {
+      // create expression from passed parameter
       typename ToTerm<PreExpr>::type expr = toTerm(std::forward<PreExpr>(preExpr));
 
-      auto gridViewExpr = makeGridViewExpression(expr, basis_.gridView());
-      DOFVector<GlobalBasis> tmp(basis(), "tmp");
-      Dune::Functions::interpolate(basis_, treePath_, tmp, gridViewExpr);
+      auto const& basis = DOFVectorConstView::basis();
+      auto const& treePath = DOFVectorConstView::treePath();
+
+      auto gridViewExpr = makeGridViewExpression(expr, basis.gridView());
+      DOFVector<GlobalBasis> tmp(basis, "tmp");
+      Dune::Functions::interpolate(basis, treePath, tmp, gridViewExpr);
 
-      coefficients_.getVector() = std::move(vec.getVector());
+      // move data from temporary vector into stored DOFVector
+      mutableDofVector_->getVector() = std::move(tmp.getVector());
       return *this;
     }
 
+    /// Interpolation of expr to DOFVector
+    template <class GridFct>
+    DOFVectorView& interpol(GridFct const& gridFct)
+    {
+      auto const& basis = DOFVectorConstView::basis();
+      auto const& treePath = DOFVectorConstView::treePath();
 
-  public:
+      DOFVector<GlobalBasis> tmp(basis, "tmp");
+      Dune::Functions::interpolate(basis, treePath, tmp, gridFct);
 
-    GlobalBasis const& basis() const { return dofVector_.getFeSpace(); }
-    TreePath const& treePath() const { return treePath_; }
+      // move data from temporary vector into stored DOFVector
+      mutableDofVector_->getVector() = std::move(tmp.getVector());
+      return *this;
+    }
 
-    DOFVector<GlobalBasis>& coefficients() { return dofVector_; }
-    DOFVector<GlobalBasis> const& coefficients() const { return dofVector_; }
+    /// Return the mutable DOFVector
+    DOFVector<GlobalBasis>& coefficients() { return *mutableDofVector_; }
 
+    /// Return the const DOFVector
+    using DOFVectorConstView::coefficients;
 
-  private:
 
-    DOFVector<GlobalBasis>& dofVector_;
-    TreePath const treePath_;
+  protected:
 
-    EntitySet entitySet_;
-    NodeToRangeEntry nodeToRangeEntry_;
+    DOFVector<GlobalBasis>* mutableDofVector_;
   };
 
+
+  // A Generator for a const \ref DOFVectorView.
+  template <class GlobalBasis, class TreePath>
+  auto makeDOFVectorView(DOFVector<GlobalBasis> const& dofVector, TreePath const& treePath)
+  {
+    return DOFVectorView<GlobalBasis, TreePath, true>{dofVector, treePath};
+  }
+
+  // A Generator for a mutable \ref DOFVectorView.
+  template <class GlobalBasis, class TreePath>
+  auto makeDOFVectorView(DOFVector<GlobalBasis>& dofVector, TreePath const& treePath)
+  {
+    return DOFVectorView<GlobalBasis, TreePath, false>{dofVector, treePath};
+  }
+
+
+  // A Generator for a const \ref DOFVectorView.
+  template <class GlobalBasis>
+  auto makeDOFVectorView(DOFVector<GlobalBasis> const& dofVector)
+  {
+    auto treePath = Dune::TypeTree::hybridTreePath();
+    return DOFVectorView<GlobalBasis, decltype(treePath), true>{dofVector, treePath};
+  }
+
+  // A Generator for a mutable \ref DOFVectorView.
+  template <class GlobalBasis>
+  auto makeDOFVectorView(DOFVector<GlobalBasis>& dofVector)
+  {
+    auto treePath = Dune::TypeTree::hybridTreePath();
+    return DOFVectorView<GlobalBasis, decltype(treePath), false>{dofVector, treePath};
+  }
+
 } // end namespace AMDiS
 
 #include "DOFVectorView.inc.hpp"

dune/amdis/DOFVectorView.inc.hpp

 #pragma once
 
-namespace AMDiS
+namespace AMDiS {
+
+template <class GlobalBasis, class TreePath>
+typename DOFVectorView<GlobalBasis, TreePath, true>::Range DOFVectorView<GlobalBasis, TreePath, true>::
+LocalFunction::operator()(LocalDomain const& x) const
+{
+  assert( bound_ );
+  auto y = Range(0);
+
+  auto&& coefficients = *globalFunction_->dofVector_;
+  auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
+  forEachLeafNode(*subTree_, [&,this](auto const& node, auto)
+  {
+    using Node = std::decay_t<decltype(node)>;
+    using LocalBasisRange = typename LocalFunction::template LocalBasisRange<Node>;
+    using MultiIndex = typename LocalIndexSet::MultiIndex;
+    using CoefficientBlock = typename std::decay<decltype(std::declval<Vector>()[std::declval<MultiIndex>()])>::type;
+    using RangeBlock = std::decay_t<decltype(std::declval<NodeToRangeEntry>()(node, y))>;
+
+    auto&& fe = node.finiteElement();
+    auto&& localBasis = fe.localBasis();
+
+    auto&& shapeFunctionValues = shapeFunctionValueContainer_[node];
+    localBasis.evaluateFunction(x, shapeFunctionValues);
+
+    // Get range entry associated to this node
+    auto&& re = nodeToRangeEntry(node, y);
+
+    for (std::size_t i = 0; i < localBasis.size(); ++i) {
+      auto&& multiIndex = localIndexSet_.index(node.localIndex(i));
+
+      // Get coefficient associated to i-th shape function
+      auto&& c = coefficients[multiIndex];
+
+      // Get value of i-th shape function
+      auto&& v = 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.
+      std::size_t dimC = Flat<CoefficientBlock>::size(c);
+      std::size_t dimV = Flat<LocalBasisRange>::size(v);
+      assert(dimC*dimV == Flat<RangeBlock>::size(re));
+      for(std::size_t j = 0; j < dimC; ++j) {
+        auto&& c_j = Flat<CoefficientBlock>::getEntry(c, j);
+        for(std::size_t k = 0; k < dimV; ++k) {
+          auto&& v_k = Flat<LocalBasisRange>::getEntry(v, k);
+          Flat<RangeBlock>::getEntry(re, j*dimV + k) += c_j*v_k;
+        }
+      }
+    }
+  });
+
+  return y;
+}
+
+template <class GlobalBasis, class TreePath>
+typename DOFVectorView<GlobalBasis, TreePath, true>::DerivativeRange DOFVectorView<GlobalBasis, TreePath, true>::
+GradientLocalFunction::operator()(LocalDomain const& x) const
 {
+  assert( bound_ );
+  Range dy;
+  for (std::size_t j = 0; j < dy.size(); ++j)
+    dy[j] = 0;
+
+  auto&& coefficients = *globalFunction_->dofVector_;
+  auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
+  forEachLeafNode(*subTree_, [&,this](auto const& node, auto)
+  {
+    // TODO: may DOFVectorView::Range to FieldVector type if necessary
+    using LocalDerivativeTraits = Dune::Functions::DefaultDerivativeTraits<Dune::FieldVector<double,1>(Domain)>;
+    using GradientBlock = typename LocalDerivativeTraits::Range;
+    using MultiIndex = typename LocalIndexSet::MultiIndex;
+    using CoefficientBlock = typename std::decay<decltype(std::declval<Vector>()[std::declval<MultiIndex>()])>::type;
+    using RangeBlock = std::decay_t<decltype(std::declval<NodeToRangeEntry>()(node, dy))>;
+
+    auto&& fe = node.finiteElement();
+    auto&& localBasis = fe.localBasis();
+
+    // The transposed inverse Jacobian of the map from the reference element to the element
+    auto&& jacobian = geometry_.value().jacobianInverseTransposed(x);
+
+    auto&& referenceGradients = referenceGradientContainer_[node];
+    localBasis.evaluateJacobian(x, referenceGradients);
+
+    // Compute the shape function gradients on the real element
+    std::vector<GradientBlock> gradients(referenceGradients.size());
+    for (std::size_t i = 0; i < gradients.size(); ++i)
+      multiplies_ABt(referenceGradients[i], jacobian, gradients[i]); // D[phi] * J^(-1) -> grad
+
+    // Get range entry associated to this node
+    auto&& re = nodeToRangeEntry(node, dy);
+
+    for (std::size_t i = 0; i < localBasis.size(); ++i) {
+      auto&& multiIndex = localIndexSet_.index(node.localIndex(i));
+
+      // Get coefficient associated to i-th shape function
+      auto&& c = coefficients[multiIndex];
+
+      // Get value of i-th transformed reference gradient
+      auto&& grad = gradients[i];
 
-      {
-        auto y = Range(0);
-
-        auto const& coefficients = globalFunction_.coefficients();
-        auto const& nodeToRangeEntry = globalFunction_.nodeToRangeEntry();
-        forEachLeafNode(subTree_, [&,this](auto const& node, auto)
-        {
-          using Node = std::decay_t<decltype(node)>;
-          using LocalBasisRange = typename LocalFunction::template LocalBasisRange<Node>;
-          using MultiIndex = typename LocalIndexSet::MultiIndex;
-          using CoefficientBlock = Range;
-          using RangeBlock = std::decay_t<decltype(nodeToRangeEntry(node, y))>;
-
-          auto&& fe = node.finiteElement();
-          auto&& localBasis = fe.localBasis();
-
-          auto&& shapeFunctionValues = shapeFunctionValueContainer_[node];
-          localBasis.evaluateFunction(x, shapeFunctionValues);
-
-          // Get range entry associated to this node
-          auto&& re = nodeToRangeEntry(node, y);
-
-          for (std::size_t i = 0; i < localBasis.size(); ++i) {
-            auto&& multiIndex = localIndexSet_.index(node.localIndex(i));
-
-            // Get coefficient associated to i-th shape function
-            auto&& c = coefficients[multiIndex];
-
-            // Get value of i-th shape function
-            auto&& v = 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.
-            auto dimC = Flat<CoefficientBlock>::size(c);
-            auto dimV = Flat<LocalBasisRange>::size(v);
-            assert(dimC*dimV == Flat<RangeBlock>::size(re));
-            for(std::size_t j = 0; j < dimC; ++j) {
-              auto&& c_j = Flat<CoefficientBlock>::getEntry(c, j);
-              for(std::size_t k = 0; k < dimV; ++k) {
-                auto&& v_k = Flat<LocalBasisRange>::getEntry(v, k);
-                Flat<RangeBlock>::getEntry(re, j*dimV + k) += c_j*v_k;
-              }
-            }
-          }
-        });
-
-        return y;
+      // Notice that the range entry re, the coefficient c, and the transformed
+      // reference gradient grad may all be scalar, vector, matrix, or general container valued.
+      std::size_t dimC = Flat<CoefficientBlock>::size(c);
+      std::size_t dimV = Flat<GradientBlock>::size(grad);
+      assert(dimC*dimV == Flat<RangeBlock>::size(re));
+      for(std::size_t j = 0; j < dimC; ++j) {
+        auto&& c_j = Flat<CoefficientBlock>::getEntry(c, j);
+        for(std::size_t k = 0; k < dimV; ++k) {
+          auto&& v_k = Flat<GradientBlock>::getEntry(grad, k);
+          Flat<RangeBlock>::getEntry(re, j*dimV + k) += c_j*v_k;
+        }
       }
+    }
+  });
 
+  return dy;
+}
 
 } // end namespace AMDiS

dune/amdis/common/FieldMatVec.hpp

@@ -107,7 +107,7 @@ namespace AMDiS
   template <class T, int N>
   auto unary_dot(Dune::FieldVector<T, N> const& x)
   {
-    auto op = [](auto const& a, auto const& b) { return a + sqr(std::abs(b)); };
+    auto op = [](auto const& a, auto const& b) { return a + Math::sqr(std::abs(b)); };
     return Impl::accumulate(x, op);
   }
 
@@ -187,7 +187,7 @@ namespace AMDiS
   {
     T result = 0;
     for (int i = 0; i < N; ++i)
-      result += sqr(lhs[i] - rhs[i]);
+      result += Math::sqr(lhs[i] - rhs[i]);
     return std::sqrt(result);
   }
 
@@ -274,4 +274,67 @@ namespace AMDiS
     return sqrt(dot(DF[0], DF[0]));
   }
 
+  template <class T, int M, int N>
+  Dune::FieldMatrix<T,N,M> trans(Dune::FieldMatrix<T, M, N> const& A)
+  {
+    Dune::FieldMatrix<T,N,M> At;
+    for (int i = 0; i < M; ++i)
+      for (int j = 0; j < N; ++j)
+        At[j][i] = A[i][j];
+
+    return At;
+  }
+
+
+  template <class T, int M, int N, int L>
+  Dune::FieldMatrix<T,M,N> multiplies(Dune::FieldMatrix<T, M, L> const& A,  Dune::FieldMatrix<T, L, N> const& B)
+  {
+    return A.rightmultiplyany(B);
+  }
+
+
+
+  template <class T, int M, int N, int L>
+  Dune::FieldMatrix<T,M,N> multiplies_AtB(Dune::FieldMatrix<T, L, M> const& A,  Dune::FieldMatrix<T, N, L> const& B)
+  {
+    Dune::FieldMatrix<T,M,N> C;
+
+    for (int m = 0; m < M; ++m) {
+      for (int n = 0; n < N; ++n) {
+        C[m][n] = 0;
+        for (int l = 0; l < L; ++l)
+          C[m][n] += A[l][m] * B[n][l];
+      }
+    }