Merge branch 'feature/dofvector_global_evaluation' into 'develop'

Feature/dofvector global evaluation

See merge request spraetor/dune-amdis!48

src/amdis/FileWriter.hpp

@@ -12,7 +12,7 @@
 #include <amdis/Initfile.hpp>
 #include <amdis/common/Size.hpp>
 #include <amdis/common/ValueCategory.hpp>
-#include <amdis/gridfunctions/DOFVectorView.hpp>
+#include <amdis/gridfunctions/DiscreteFunction.hpp>
 #include <amdis/io/FileWriterInterface.hpp>
 #include <amdis/utility/Filesystem.hpp>
 
@@ -43,14 +43,17 @@ namespace AMDiS
   constexpr std::size_t VTKFieldSize = Size<Range>;
 
 
-  template <class GlobalBasis, class RangeType, class TreePath>
+  template <class GB, class VT, class TP>
   class FileWriter
       : public FileWriterInterface
   {
-  private: // typedefs and static constants
+    using GlobalBasis = GB;
+    using ValueType = VT;
+    using TreePath = TP;
+
+  private:
     using GridView = typename GlobalBasis::GridView;
-    using Vector = DOFVectorConstView<GlobalBasis,RangeType,TreePath>;
-    using Range = typename Vector::Range;
+    using Range = typename DiscreteFunction<GB,VT,TP>::Range;
 
     /// Dimension of the mesh
     static constexpr int dim = GridView::dimension;
@@ -60,33 +63,32 @@ namespace AMDiS
 
   public:
     /// Constructor.
-    FileWriter(std::string const& baseName,
-               Vector const& dofvector)
-      : FileWriterInterface(baseName)
-      , dofvector_(dofvector)
+    FileWriter(std::string const& name, DiscreteFunction<GB,VT,TP> const& discreteFct)
+      : FileWriterInterface(name)
+      , discreteFct_(discreteFct)
       , animation_(false)
     {
-      Parameters::get(baseName + "->ParaView animation", animation_);
+      Parameters::get(name + "->ParaView animation", animation_);
 
-      int m = Parameters::get<int>(baseName + "->ParaView mode").value_or(0);
+      int m = Parameters::get<int>(name + "->ParaView mode").value_or(0);
       mode_ = (m == 0 ? Dune::VTK::ascii : Dune::VTK::appendedraw);
 
-      init(baseName, ValueCategory_t<Range>{});
+      init(name, ValueCategory_t<Range>{});
     }
 
     template <class ValueCat>
     void init(std::string const& baseName, ValueCat)
     {
-      int subSampling = Parameters::get<int>(baseName + "->subsampling").value_or(0);
-      if (subSampling > 0)
-        vtkWriter_ = std::make_shared<Dune::SubsamplingVTKWriter<GridView>>(gridView(), subSampling);
+      auto subSampling = Parameters::get<int>(baseName + "->subsampling");
+      if (subSampling)
+        vtkWriter_ = std::make_shared<Dune::SubsamplingVTKWriter<GridView>>(gridView(), subSampling.value());
       else
         vtkWriter_ = std::make_shared<Dune::VTKWriter<GridView>>(gridView());
 
       if (animation_)
         vtkSeqWriter_ = std::make_shared<Dune::VTKSequenceWriter<GridView>>(vtkWriter_, filename_, dir_, "");
 
-      vtkWriter_->addVertexData(dofvector_, Dune::VTK::FieldInfo(name_, VTKFieldType<Range>, VTKFieldSize<Range>));
+      vtkWriter_->addVertexData(discreteFct_, Dune::VTK::FieldInfo(name_, VTKFieldType<Range>, VTKFieldSize<Range>));
     }
 
     void init(std::string const&, tag::unknown) {}
@@ -105,11 +107,11 @@ namespace AMDiS
   protected:
     GridView const& gridView() const
     {
-      return dofvector_.basis().gridView();
+      return discreteFct_.basis().gridView();
     }
 
   private:
-    Vector dofvector_;
+    DiscreteFunction<GB,VT,TP> discreteFct_;
 
     std::shared_ptr<Dune::VTKWriter<GridView>> vtkWriter_;
     std::shared_ptr<Dune::VTKSequenceWriter<GridView>> vtkSeqWriter_;
@@ -122,12 +124,12 @@ namespace AMDiS
   };
 
 
-  template <class GlobalBasis, class Range, class TreePath>
-  std::shared_ptr<FileWriter<GlobalBasis,Range,TreePath>>
-  makeFileWriterPtr(std::string baseName,
-      DOFVectorConstView<GlobalBasis,Range,TreePath> const& dofvector)
+  template <class GlobalBasis, class ValueType, class TreePath>
+  std::shared_ptr<FileWriter<GlobalBasis,ValueType,TreePath>>
+  makeFileWriterPtr(std::string const& name,
+      DiscreteFunction<GlobalBasis,ValueType,TreePath> const& discreteFct)
   {
-    return std::make_shared<FileWriter<GlobalBasis,Range,TreePath>>(baseName, dofvector);
+    return std::make_shared<FileWriter<GlobalBasis,ValueType,TreePath>>(name, discreteFct);
   }
 
 } // end namespace AMDiS

src/amdis/ProblemStat.hpp

@@ -32,6 +32,7 @@
 #include <amdis/common/Utility.hpp>
 
 #include <amdis/GridFunctions.hpp>
+#include <amdis/gridfunctions/DiscreteFunction.hpp>
 #include <amdis/gridfunctions/DOFVectorView.hpp>
 
 #include <amdis/io/FileWriterInterface.hpp>
@@ -198,7 +199,7 @@ namespace AMDiS
     auto getSolution(TreePath const& path = {}) const
     {
       auto&& tp = makeTreePath(path);
-      return makeDOFVectorView(*solution_, tp);
+      return makeDiscreteFunction(*solution_, tp);
     }
 
 

src/amdis/gridfunctions/CMakeLists.txt

@@ -5,8 +5,9 @@ install(FILES
     ConstantGridFunction.hpp
     CoordsGridFunction.hpp
     DerivativeGridFunction.hpp
+    DiscreteFunction.hpp
+    DiscreteFunction.inc.hpp
     DOFVectorView.hpp
-    DOFVectorView.inc.hpp
     FunctorGridFunction.hpp
     GridFunctionConcepts.hpp
     Integrate.hpp

src/amdis/gridfunctions/DOFVectorView.hpp

 #pragma once
 
-#include <vector>
-
-#include <dune/common/std/optional.hh>
-#include <dune/functions/common/defaultderivativetraits.hh>
-#include <dune/functions/common/treedata.hh>
-#include <dune/functions/functionspacebases/defaultnodetorangemap.hh>
-#include <dune/functions/functionspacebases/flatvectorview.hh>
-#include <dune/functions/gridfunctions/gridviewentityset.hh>
-#include <dune/typetree/childextraction.hh>
-
 #include <amdis/GridFunctions.hpp>
-#include <amdis/utility/FiniteElementType.hpp>
+#include <amdis/gridfunctions/DiscreteFunction.hpp>
 
 namespace AMDiS
 {
-  /**
-    * \addtogroup GridFunctions
-    * @{
-    **/
-
-  template <class GlobalBasisType, class RangeType, class TreePathType>
-  class DOFVectorConstView
+  /// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
+  template <class GB, class VT, class TP>
+  class DOFVectorView
+      : public DiscreteFunction<GB, VT, TP>
   {
-  public:
-    using GlobalBasis = GlobalBasisType;
-    using TreePath = TreePathType;
-
-    using Tree = typename GlobalBasis::LocalView::Tree;
-    using SubTree = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
-    using NodeToRangeEntry = Dune::Functions::DefaultNodeToRangeMap<SubTree>;
-
-    using GridView = typename GlobalBasis::GridView;
-    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
-
-    using Domain = typename EntitySet::GlobalCoordinate;
-    using Range = RangeType_t<SubTree>;
-    static_assert(std::is_arithmetic<RangeType>::value, "");
-    // Don't know how to determine Range with non-trivial RangeType
-
-    using RawSignature = typename Dune::Functions::SignatureTraits<Range(Domain)>::RawSignature;
-    using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>;
-    using DerivativeRange = typename DerivativeTraits::Range;
-
-    using LocalDomain = typename EntitySet::LocalCoordinate;
-    using Element = typename EntitySet::Element;
-    using Geometry = typename Element::Geometry;
-
-    enum { hasDerivative = false };
-
-  public: // a local view on the gradients
-
-    /// A LocalFunction representing the derivative of the DOFVector
-    class GradientLocalFunction
-    {
-    public:
-      using Domain = LocalDomain;
-      using Range = DerivativeRange;
-
-      enum { hasDerivative = false };
-
-    private:
-      using LocalView = typename GlobalBasis::LocalView;
-
-      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(DOFVectorConstView const& globalFunction)
-        : globalFunction_(&globalFunction)
-        , localView_(globalFunction_->basis().localView())
-        , subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_->treePath()))
-      {
-        referenceGradientContainer_.init(*subTree_);
-      }
-
-      void bind(Element const& element)
-      {
-        localView_.bind(element);
-        geometry_.emplace(element.geometry());
-        bound_ = true;
-      }
-
-      void unbind()
-      {
-        localView_.unbind();
-        geometry_.reset();
-        bound_ = false;
-      }
-
-      /// Evaluate Gradient at bound element in local coordinates
-      Range operator()(Domain const& x) const;
-
-      friend int order(GradientLocalFunction const& self)
-      {
-        assert( self.bound_ );
-        return std::max(0, polynomialDegree(*self.subTree_)-1);
-      }
-
-      /// Return the bound element
-      Element const& localContext() const
-      {
-        assert( bound_ );
-        return localView_.element();
-      }
-
-    private:
-      DOFVectorConstView const* globalFunction_;
-      LocalView localView_;
-
-      SubTree const* subTree_;
-      mutable ReferenceGradientContainer referenceGradientContainer_;
-
-      Dune::Std::optional<Geometry> geometry_;
-      bool bound_ = false;
-    };
-
-
-  public: // a local view on the values
-
-    /// A LocalFunction, i.e., an element local view on the DOFVector
-    class LocalFunction
-    {
-    public:
-      using Domain = typename DOFVectorConstView::LocalDomain;
-      using Range = typename DOFVectorConstView::Range;
-
-      enum { hasDerivative = true };
-
-    private:
-      using LocalView = typename GlobalBasis::LocalView;
-
-      template <class LeafNode>
-      using LocalBasisRange = RangeType_t<LeafNode>;
-        // = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
-
-      template <class Node>
-      using NodeData = typename std::vector<LocalBasisRange<Node>>;
-
-      using ShapeFunctionValueContainer = Dune::Functions::TreeData<SubTree, NodeData, true>;
-
-    public:
-      LocalFunction(DOFVectorConstView const& globalFunction)
-        : globalFunction_(&globalFunction)
-        , localView_(globalFunction_->basis().localView())
-        , subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_->treePath()))
-      {
-        shapeFunctionValueContainer_.init(*subTree_);
-      }
-
-      void bind(Element const& element)
-      {
-        localView_.bind(element);
-        bound_ = true;
-      }
-
-      void unbind()
-      {
-        localView_.unbind();
-        bound_ = false;
-      }
-
-      /// Evaluate LocalFunction at bound element in local coordinates
-      Range operator()(Domain const& x) const;
-
-      /// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction
-      friend GradientLocalFunction derivative(LocalFunction const& localFunction)
-      {
-        static_assert(isValidRange<DerivativeTraits>(),"Derivative of DOFVector not defined.");
-        return GradientLocalFunction{*localFunction.globalFunction_};
-      }
-
-      friend int order(LocalFunction const& self)
-      {
-        assert( self.bound_ );
-        return polynomialDegree(*self.subTree_);
-      }
-
-      /// Return the bound element
-      Element const& localContext() const
-      {
-        assert( bound_ );
-        return localView_.element();
-      }
-
-    private:
-      DOFVectorConstView const* globalFunction_;
-      LocalView localView_;
-      SubTree const* subTree_;
-
-      mutable ShapeFunctionValueContainer shapeFunctionValueContainer_;
-
-      bool bound_ = false;
-    };
+    using Self = DOFVectorView;
+    using Super = DiscreteFunction<GB, VT, TP>;
 
+    using GlobalBasis = GB;
+    using TreePath = TP;
 
   public:
-    /// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
-    DOFVectorConstView(DOFVector<GlobalBasis,RangeType> const& dofVector, TreePath const& treePath)
-      : dofVector_(&dofVector)
-      , treePath_(treePath)
-      , entitySet_(dofVector.basis().gridView())
-      , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(dofVector.basis(), treePath))
+    /// Constructor forwards to the treePath constructor, with empty TreePath
+    DOFVectorView(DOFVector<GB,VT>& dofVector)
+      : DOFVectorView{dofVector, Dune::TypeTree::hybridTreePath()}
     {}
 
-    /// Evaluate the view on this DOFVector in global coordinates
-    Range operator()(Domain const& x) const
-    {
-      error_exit("Not implemented.");
-      return Range(0);
-    }
-
-    /// \brief Create a local function for this view on the DOFVector. \relates LocalFunction
-    friend LocalFunction localFunction(DOFVectorConstView const& self)
-    {
-      return LocalFunction{self};
-    }
-
-    EntitySet const& entitySet() const
-    {
-      return entitySet_;
-    }
-
-  public:
-    /// Return global basis
-    GlobalBasis const& basis() const
-    {
-      return dofVector_->basis();
-    }
-
-    /// Return treePath associated with this view
-    TreePath const& treePath() const
-    {
-      return treePath_;
-    }
-
-    /// Return const coefficient vector
-    DOFVector<GlobalBasis,RangeType> const& coefficients() const
-    {
-      return *dofVector_;
-    }
-
-  protected:
-    DOFVector<GlobalBasis,RangeType> const* dofVector_;
-    TreePath const treePath_;
-
-    EntitySet entitySet_;
-    NodeToRangeEntry nodeToRangeEntry_;
-  };
-
-
-  // A mutable version of DOFVectorView
-  template <class GlobalBasisType, class RangeType, class TreePathType>
-  class DOFVectorMutableView
-      : public DOFVectorConstView<GlobalBasisType, RangeType, TreePathType>
-  {
-    using Super = DOFVectorConstView<GlobalBasisType, RangeType, TreePathType>;
-
-    using GlobalBasis = GlobalBasisType;
-    using TreePath = TreePathType;
-
-  public:
     /// Constructor. Stores a pointer to the mutable `dofvector`.
-    DOFVectorMutableView(DOFVector<GlobalBasis,RangeType>& dofVector, TreePath const& treePath)
+    DOFVectorView(DOFVector<GB,VT>& dofVector, TP const& treePath)
       : Super(dofVector, treePath)
       , mutableDofVector_(&dofVector)
     {}
 
   public:
-    /// Interpolation of GridFunction to DOFVector
+    /// \brief Interpolation of GridFunction to DOFVector, assuming that there is no
+    /// reference to this DOFVector in the expression.
+    /**
+     * **Example:**
+     * ```
+     * auto v = makeDOFVectorView(prob.getSolutionVector(),0);
+     * v.interpolate([](auto const& x) { return x[0]; });
+     * ```
+     **/
     template <class Expr>
-    DOFVectorMutableView& interpolate(Expr&& expr)
+    void interpolate_noalias(Expr&& expr)
     {
-      auto const& basis = Super::basis();
-      auto const& treePath = Super::treePath();
+      auto const& basis = this->basis();
+      auto const& treePath = this->treePath();
 
       auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), basis.gridView());
 
-      DOFVector<GlobalBasis,RangeType> tmp(*mutableDofVector_);
-      Dune::Functions::interpolate(basis, treePath, tmp, std::forward<decltype(gridFct)>(gridFct));
+      Dune::Functions::interpolate(basis, treePath, coefficients(),
+        std::forward<decltype(gridFct)>(gridFct));
+    }
+
+    /// \brief Interpolation of GridFunction to DOFVector
+    /**
+     * **Example:**
+     * ```
+     * auto v = makeDOFVectorView(prob.getSolutionVector(),0);
+     * v.interpolate(v + [](auto const& x) { return x[0]; });
+     * ```
+     * Allows to have a reference to the DOFVector in the expression, e.g. as
+     * \ref DiscreteFunction or \ref gradientAtQP() of a DiscreteFunction.
+     **/
+    template <class Expr>
+    void interpolate(Expr&& expr)
+    {
+      // create temporary copy of data
+      DOFVector<GB,VT> tmp(coefficients());
+      Self tmpView{tmp, this->treePath()};
+      tmpView.interpolate_noalias(std::forward<Expr>(expr));
 
       // move data from temporary vector into stored DOFVector
-      mutableDofVector_->vector() = std::move(tmp.vector());
+      coefficients().vector() = std::move(tmp.vector());
+    }
+
+    /// \brief Interpolation of GridFunction to DOFVector, alias to \ref interpolate()
+    template <class Expr>
+    DOFVectorView& operator<<(Expr&& expr)
+    {
+      interpolate(expr);
       return *this;
     }
 
+    /// \brief interpolate `(*this) + expr` to DOFVector
     template <class Expr>
-    DOFVectorMutableView& operator<<(Expr&& expr)
+    DOFVectorView& operator+=(Expr&& expr)
     {
-      return interpolate(expr);
+      interpolate((*this) + expr);
+      return *this;
     }
 
+    /// \brief interpolate `(*this) - expr` to DOFVector
+    template <class Expr>
+    DOFVectorView& operator-=(Expr&& expr)
+    {
+      interpolate((*this) - expr);
+      return *this;
+    }
 
     /// Return the mutable DOFVector
-    DOFVector<GlobalBasis,RangeType>& coefficients() { return *mutableDofVector_; }
+    DOFVector<GB,VT>& coefficients() { return *mutableDofVector_; }
 
     /// Return the const DOFVector
     using Super::coefficients;
 
   protected:
-    DOFVector<GlobalBasis,RangeType>* mutableDofVector_;
+    DOFVector<GB,VT>* mutableDofVector_;
   };
 
-  /** @} **/
-
-
-#ifndef DOXYGEN
-  // A Generator for a const \ref DOFVectorView.
-  template <class GlobalBasis, class RangeType, class TreePath>
-  auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType> const& dofVector, TreePath const& treePath)
-  {
-    return DOFVectorConstView<GlobalBasis, RangeType, TreePath>{dofVector, treePath};
-  }
-
-  // A Generator for a mutable \ref DOFVectorView.
-  template <class GlobalBasis, class RangeType, class TreePath>
-  auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType>& dofVector, TreePath const& treePath)
-  {
-    return DOFVectorMutableView<GlobalBasis, RangeType, TreePath>{dofVector, treePath};
-  }
 
+#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
+  // Deduction guide for DOFVectorView class
+  template <class GlobalBasis, class ValueType>
+  DOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector)
+    -> DOFVectorView<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>;
+#endif
 
-  // A Generator for a const \ref DOFVectorView.
-  template <class GlobalBasis, class RangeType>
-  auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType> const& dofVector)
+  /// A Generator for a mutable \ref DOFVectorView
+  template <class GlobalBasis, class ValueType, class TreePath>
+  auto makeDOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector, TreePath const& treePath)
   {
-    auto treePath = Dune::TypeTree::hybridTreePath();
-    return DOFVectorConstView<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
+    return DOFVectorView<GlobalBasis, ValueType, TreePath>{dofVector, treePath};
   }
 
-  // A Generator for a mutable \ref DOFVectorView.
-  template <class GlobalBasis, class RangeType>
-  auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType>& dofVector)
+  /// A Generator for a mutable \ref DOFVectorView
+  template <class GlobalBasis, class ValueType>
+  auto makeDOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector)
   {
     auto treePath = Dune::TypeTree::hybridTreePath();
-    return DOFVectorMutableView<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
+    return DOFVectorView<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>{dofVector, treePath};
   }
-#endif
 
 } // end namespace AMDiS
-
-#include "DOFVectorView.inc.hpp"

src/amdis/gridfunctions/DerivativeGridFunction.hpp

@@ -3,6 +3,7 @@
 #include <type_traits>
 
 #include <dune/functions/common/defaultderivativetraits.hh>
+#include <dune/grid/utility/hierarchicsearch.hh>
 
 #include <amdis/gridfunctions/GridFunctionConcepts.hpp>
 
@@ -52,11 +53,22 @@ namespace AMDiS
       static_assert(isValidRange<DerivativeTraits>(), "Derivative of GridFunction not defined");
     }
 
-    /// NOTE: no global derivative available
+    /// Evaluate derivative in global coordinates. NOTE: expensive
     Range operator()(Domain const& x) const
     {
-      error_exit("Not implemented");
-      return Range(0);
+      auto gv = entitySet().gridView();
+
+      using GridView = decltype(gv);
+      using Grid = typename GridView::Grid;
+      using IS = typename GridView::IndexSet;
+
+      Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};
+
+      auto element = hsearch.findEntity(x);
+      auto geometry = element.geometry();
+      auto localFct = derivative(localFunction(gridFct_));
+      localFct.bind(element);
+      return localFct(geometry.local(x));
     }
 
     /// Return the derivative-localFunction of the GridFunction.

src/amdis/gridfunctions/DiscreteFunction.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/flatvectorview.hh>
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+#include <dune/typetree/childextraction.hh>
+
+#include <amdis/GridFunctions.hpp>
+#include <amdis/utility/FiniteElementType.hpp>
+
+namespace AMDiS
+{
+  /// \brief A view on a subspace of a \ref DOFVector
+  /**
+   * \ingroup GridFunctions
+   *
+   * \tparam GB  GlobalBasis type that models \ref Dune::Functions::Concept::GlobalBasis
+   * \tparam VT  Coefficient type of the DOFVector
+   * \tparam TP  TreePath type a realization of \ref Dune::TypeTree::HybridTreePath
+   **/
+  template <class GB, class VT, class TP>
+  class DiscreteFunction
+  {
+    static_assert(std::is_arithmetic<VT>::value, "");
+
+  private:
+    using GlobalBasis = GB;
+    using TreePath = TP;
+
+    using Tree = typename GlobalBasis::LocalView::Tree;
+    using SubTree = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
+    using NodeToRangeEntry = Dune::Functions::DefaultNodeToRangeMap<SubTree>;
+
+    using GridView = typename GlobalBasis::GridView;
+
+  public:
+    /// Set of entities the DiscreteFunction is defined on
+    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
+
+    /// Global coordinates of the EntitySet
+    using Domain = typename EntitySet::GlobalCoordinate;
+
+    /// Range type of this DiscreteFunction
+    using Range = RangeType_t<SubTree,VT>;
+
+    /// \brief This GridFunction has no derivative function, it can be created
+    /// by \ref DiscreteGridFunction.
+    enum { hasDerivative = false };
+
+  public:
+    /// A LocalFunction representing the derivative of the DOFVector on a bound element
+    class GradientLocalFunction;
+
+    /// A LocalFunction representign the value the DOFVector on a bound element
+    class LocalFunction;
+
+  public:
+    /// Constructor forwards to the treePath constructor, with empty TreePath
+    DiscreteFunction(DOFVector<GB,VT> const& dofVector)
+      : DiscreteFunction{dofVector, Dune::TypeTree::hybridTreePath()}
+    {}
+
+    /// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
+    DiscreteFunction(DOFVector<GB,VT> const& dofVector, TP const& treePath)
+      : dofVector_(&dofVector)
+      , treePath_(treePath)
+      , entitySet_(dofVector.basis().gridView())
+      , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(dofVector.basis(), treePath))
+    {}
+
+    /// \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
+    friend LocalFunction localFunction(DiscreteFunction const& self)
+    {
+      return LocalFunction{self};
+    }
+
+    /// \brief Return a \ref Dune::Functions::GridViewEntitySet
+    EntitySet const& entitySet() const
+    {
+      return entitySet_;
+    }
+
+  public:
+    /// \brief Return global basis bound to the DOFVector
+    GlobalBasis const& basis() const
+    {
+      return dofVector_->basis();
+    }
+
+    /// \brief Return treePath associated with this discrete function
+    TreePath const& treePath() const
+    {
+      return treePath_;
+    }
+
+    /// \brief Return const coefficient vector
+    DOFVector<GB,VT> const& coefficients() const
+    {
+      return *dofVector_;
+    }
+
+  protected:
+    DOFVector<GB,VT> const* dofVector_;
+    TreePath treePath_;
+    EntitySet entitySet_;
+    NodeToRangeEntry nodeToRangeEntry_;
+  };
+
+
+#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
+  // Deduction guide for DiscreteFunction class
+  template <class GlobalBasis, class ValueType>
+  DiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector)
+    -> DiscreteFunction<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>;
+#endif
+
+  /// A Generator for a \ref DiscreteFunction
+  template <class GlobalBasis, class ValueType, class TreePath>
+  auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector, TreePath const& treePath)
+  {
+    return DiscreteFunction<GlobalBasis, ValueType, TreePath>{dofVector, treePath};
+  }
+
+  /// A Generator for a \ref DiscreteFunction
+  template <class GlobalBasis, class ValueType>
+  auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector)
+  {
+    auto treePath = Dune::TypeTree::hybridTreePath();
+    return DiscreteFunction<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>{dofVector, treePath};
+  }
+
+} // end namespace AMDiS
+
+#include "DiscreteFunction.inc.hpp"

src/amdis/gridfunctions/DOFVectorView.inc.hpp → src/amdis/gridfunctions/DiscreteFunction.inc.hpp

@@ -3,11 +3,144 @@
 #include <amdis/common/FieldMatVec.hpp>
 #include <amdis/LocalBasisCache.hpp>
 
+#include <dune/grid/utility/hierarchicsearch.hh>
+
 namespace AMDiS {
 
-template <class GB, class RT, class TP>
-typename DOFVectorConstView<GB,RT,TP>::Range DOFVectorConstView<GB,RT,TP>::
-LocalFunction::operator()(LocalDomain const& x) const
+template <class GB, class VT, class TP>
+class DiscreteFunction<GB,VT,TP>::LocalFunction
+{
+public:
+  using Domain = typename EntitySet::LocalCoordinate;
+  using Range = typename DiscreteFunction::Range;
+
+  enum { hasDerivative = true };
+
+private:
+  using LocalView = typename GlobalBasis::LocalView;
+  using Element = typename EntitySet::Element;
+  using Geometry = typename Element::Geometry;
+
+public:
+  LocalFunction(DiscreteFunction const& globalFunction)
+    : globalFunction_(&globalFunction)
+    , localView_(globalFunction_->basis().localView())
+    , subTree_(&child(localView_.tree(), globalFunction_->treePath()))
+  {}
+
+  /// \brief Bind the LocalView to the element
+  void bind(Element const& element)
+  {
+    localView_.bind(element);
+    bound_ = true;
+  }
+
+  /// \brief Unbind the LocalView from the element
+  void unbind()
+  {
+    localView_.unbind();
+    bound_ = false;
+  }
+
+  /// \brief Evaluate LocalFunction at bound element in local coordinates
+  Range operator()(Domain const& x) const;
+
+  /// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction
+  friend GradientLocalFunction derivative(LocalFunction const& localFunction)
+  {
+    return GradientLocalFunction{*localFunction.globalFunction_};
+  }
+
+  /// \brief The \ref polynomialDegree() of the LocalFunctions
+  friend int order(LocalFunction const& self)
+  {
+    assert( self.bound_ );
+    return polynomialDegree(*self.subTree_);
+  }
+
+  /// \brief Return the bound element
+  Element const& localContext() const
+  {
+    assert( bound_ );
+    return localView_.element();
+  }
+
+private:
+  DiscreteFunction const* globalFunction_;
+  LocalView localView_;
+  SubTree const* subTree_;
+  bool bound_ = false;
+};
+
+
+template <class GB, class VT, class TP>
+class DiscreteFunction<GB,VT,TP>::GradientLocalFunction
+{
+  using R = typename DiscreteFunction::Range;
+  using D = typename DiscreteFunction::Domain;
+  using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
+  using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>;
+
+public:
+  using Domain = typename EntitySet::LocalCoordinate;
+  using Range = typename DerivativeTraits::Range;
+
+  enum { hasDerivative = false };
+
+private:
+  using LocalView = typename GlobalBasis::LocalView;
+  using Element = typename EntitySet::Element;
+  using Geometry = typename Element::Geometry;
+
+public:
+  GradientLocalFunction(DiscreteFunction const& globalFunction)
+    : globalFunction_(&globalFunction)
+    , localView_(globalFunction_->basis().localView())
+    , subTree_(&child(localView_.tree(), globalFunction_->treePath()))
+  {}
+
+  void bind(Element const& element)
+  {
+    localView_.bind(element);
+    geometry_.emplace(element.geometry());
+    bound_ = true;
+  }
+
+  void unbind()
+  {
+    localView_.unbind();
+    geometry_.reset();
+    bound_ = false;
+  }
+
+  /// Evaluate Gradient at bound element in local coordinates
+  Range operator()(Domain const& x) const;
+
+  friend int order(GradientLocalFunction const& self)
+  {
+    assert( self.bound_ );
+    return std::max(0, polynomialDegree(*self.subTree_)-1);
+  }
+
+  /// Return the bound element
+  Element const& localContext() const
+  {
+    assert( bound_ );
+    return localView_.element();
+  }
+
+private:
+  DiscreteFunction const* globalFunction_;
+  LocalView localView_;
+  SubTree const* subTree_;
+  Dune::Std::optional<Geometry> geometry_;
+  bool bound_ = false;
+};
+
+
+template <class GB, class VT, class TP>
+typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
+LocalFunction::operator()(Domain const& x) const
 {
   assert( bound_ );
   auto y = Range(0);
@@ -49,9 +182,9 @@ LocalFunction::operator()(LocalDomain const& x) const
 }
 
 
-template <class GB, class RT, class TP>
-typename DOFVectorConstView<GB,RT,TP>::DerivativeRange DOFVectorConstView<GB,RT,TP>::
-GradientLocalFunction::operator()(LocalDomain const& x) const
+template <class GB, class VT, class TP>
+typename DiscreteFunction<GB,VT,TP>::GradientLocalFunction::Range DiscreteFunction<GB,VT,TP>::
+GradientLocalFunction::operator()(Domain const& x) const
 {
   assert( bound_ );
   Range dy;
@@ -106,4 +239,22 @@ GradientLocalFunction::operator()(LocalDomain const& x) const
   return dy;
 }
 
+
+template <class GB, class VT, class TP>
+typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
+operator()(Domain const& x) const
+{
+  using Grid = typename GlobalBasis::GridView::Grid;
+  using IS = typename GlobalBasis::GridView::IndexSet;
+
+  auto const& gv = this->basis().gridView();
+  Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};
+
+  auto element = hsearch.findEntity(x);
+  auto geometry = element.geometry();
+  auto localFct = localFunction(*this);
+  localFct.bind(element);
+  return localFct(geometry.local(x));
+}
+
 } // end namespace AMDiS

src/amdis/utility/RangeType.hpp

@@ -11,31 +11,43 @@ namespace AMDiS
 {
   namespace Impl
   {
-    template <class Node, bool isLeaf, bool isPower, bool isComposite>
+    template <class Node, class R, bool isLeaf, bool isPower, bool isComposite>
     struct RangeTypeImpl
     {
       static_assert( isLeaf || isPower || isComposite, "Unknown node-type for range definition" );
     };
   }
 
-  template <class Node>
+  /// \brief Range type of a node in the basis tree, composed of the leaf basis
+  /// range types.
+  /**
+   * Generate the range type by recursively combining leaf range types to a
+   * hybrid node range type. Range types for PowerNodes are thereby constructed
+   * as Dune::FieldVector of child range types. CompositeNodes produce a
+   * MultiTypeVector of the difference child range types.
+   *
+   * \tparam Node  Type of a basis-tree node
+   * \tparam R     Coefficient type [double]
+   **/
+  template <class Node, class R = double>
   using RangeType_t =
-    typename Impl::RangeTypeImpl<Node, Node::isLeaf, Node::isPower, Node::isComposite>::type;
+    typename Impl::RangeTypeImpl<Node, R, Node::isLeaf, Node::isPower, Node::isComposite>::type;
 
 
   namespace Impl
   {
     // Leaf node
-    template <class Node>
-    struct RangeTypeImpl<Node, true, false, false>
+    template <class Node, class R>
+    struct RangeTypeImpl<Node, R, true, false, false>
     {
       using LocalBasis = typename Node::FiniteElement::Traits::LocalBasisType;
-      using type = typename LocalBasis::Traits::RangeType;
+      using T = typename LocalBasis::Traits::RangeType;
+      using type = std::decay_t<decltype(std::declval<R>() * std::declval<T>())>;
     };
 
     // Power node
-    template <class Node>
-    struct RangeTypeImpl<Node, false, true, false>
+    template <class Node, class R>
+    struct RangeTypeImpl<Node, R, false, true, false>
     {
       using ChildNode = typename Node::template Child<0>::type;
 
@@ -49,12 +61,12 @@ namespace AMDiS
         using type = Dune::FieldVector<T, int(Node::CHILDREN)>;
       };
 
-      using type = typename FlatType<ChildNode::isLeaf, RangeType_t<ChildNode>>::type;
+      using type = typename FlatType<ChildNode::isLeaf, RangeType_t<ChildNode,R>>::type;
     };
 
     // Composite node
-    template <class Node>
-    struct RangeTypeImpl<Node, false, false, true>
+    template <class Node, class R>
+    struct RangeTypeImpl<Node, R, false, false, true>
     {
       template <class Idx> struct RangeTypeGenerator;
       template <std::size_t... I>
@@ -62,7 +74,7 @@ namespace AMDiS
       {
         template <std::size_t J>
         using ChildNode = typename Node::template Child<J>::type;
-        using type = MultiTypeVector<RangeType_t<ChildNode<I>>...>;
+        using type = MultiTypeVector<RangeType_t<ChildNode<I>,R>...>;
       };
 
       using type = typename RangeTypeGenerator<MakeSeq_t<Node::CHILDREN>>::type;

test/CMakeLists.txt

@@ -11,6 +11,10 @@ dune_add_test(SOURCES ConceptsTest.cpp
 dune_add_test(SOURCES DOFVectorTest.cpp
   LINK_LIBRARIES amdis)
 
+dune_add_test(SOURCES DiscreteFunctionTest.cpp
+  LINK_LIBRARIES amdis
+  CMD_ARGS "${CMAKE_SOURCE_DIR}/examples/init/ellipt.dat.2d")
+
 dune_add_test(SOURCES ExpressionsTest.cpp
   LINK_LIBRARIES amdis)
 

test/DiscreteFunctionTest.cpp

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#include <iostream>
+#include <type_traits>
+
+#include <amdis/AMDiS.hpp>
+#include <amdis/ProblemStat.hpp>
+#include <amdis/Operators.hpp>
+#include <amdis/common/Literals.hpp>
+#include <amdis/common/FieldMatVec.hpp>
+#include <amdis/gridfunctions/Integrate.hpp>
+
+#include "Tests.hpp"
+
+using namespace AMDiS;
+
+using ElliptParam   = YaspGridBasis<2, 2>;
+using ElliptProblem = ProblemStat<ElliptParam>;
+
+template <class GB, class T>
+bool comp(DOFVector<GB,T> const& U, DOFVector<GB,T> const& V)
+{
+  if (U.size() != V.size())
+    return false;
+
+  using Int = typename DOFVector<GB,T>::size_type;
+  for (Int i = 0; i < U.size(); ++i) {
+    if (U[i] != V[i])
+      return false;
+  }
+  return true;
+}
+
+int main(int argc, char** argv)
+{
+  AMDiS::init(argc, argv);
+
+  using namespace Dune::Indices;
+
+  ElliptProblem prob("ellipt");
+  prob.initialize(INIT_ALL);
+
+  // create 3 copies of the solution vector
+  auto U0 = prob.getSolutionVector();
+  auto U1 = prob.getSolutionVector();
+  auto U2 = prob.getSolutionVector();
+
+  auto u0 = makeDOFVectorView(U0);
+  auto u1 = makeDOFVectorView(U1);
+  auto u2 = makeDOFVectorView(U2);
+
+  auto expr = [](auto const& x) { return x[0] + x[1]; };
+  u0.interpolate_noalias(expr);
+  u1.interpolate(expr);
+  u2 << expr;
+
+  AMDIS_TEST( comp(U0, U1) );
+  AMDIS_TEST( comp(U0, U2) );
+
+  auto expr2 = [](auto const& x) { return 2*x[0] - 3*x[1]; };
+
+  u0.interpolate_noalias(u2 + expr2);
+  u1.interpolate(u1 + expr2);
+  u2 += expr2;
+
+  AMDIS_TEST( comp(U0, U1) );
+  AMDIS_TEST( comp(U0, U2) );
+
+  auto expr3 = [](auto const& x) { return -0.5*x[0] - 2*x[1]; };
+
+  u0.interpolate_noalias(u2 - expr3);
+  u1.interpolate(u1 - expr3);
+  u2 -= expr3;
+
+  AMDIS_TEST( comp(U0, U1) );
+  AMDIS_TEST( comp(U0, U2) );
+
+  auto du0 = derivative(u0);
+  auto localFct = localFunction(u0);
+  auto dlocalFct = derivative(localFct);
+  for (auto const& e : elements(prob.gridView()))
+  {
+    localFct.bind(e);
+    auto geo = e.geometry();
+    auto local = referenceElement(geo).position(0,0);
+    auto y0 = localFct(local);
+    auto y1 = u0(geo.global(local));
+
+    AMDIS_TEST(std::abs(y0 - y1) < 1.e-10);
+    localFct.unbind();
+
+    dlocalFct.bind(e);
+    auto g0 = dlocalFct(local);
+    auto g1 = du0(geo.global(local));
+
+    AMDIS_TEST(two_norm(g0 - g1) < 1.e-10);
+    dlocalFct.unbind();
+  }
+
+  auto V0 = makeDOFVector<float>(prob.globalBasis());
+  auto v0 = makeDOFVectorView(V0);
+  v0 << expr;
+
+  AMDiS::finalize();
+  return 0;
+}