renamed DOFVectorConstView into DiscreteFunction and added some more interpolation methods

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>
 
@@ -49,8 +49,8 @@ namespace AMDiS
   {
   private: // typedefs and static constants
     using GridView = typename GlobalBasis::GridView;
-    using Vector = DOFVectorConstView<GlobalBasis,RangeType,TreePath>;
-    using Range = typename Vector::Range;
+    using DiscreteFunction = AMDiS::DiscreteFunction<GlobalBasis,RangeType,TreePath>;
+    using Range = typename DiscreteFunction::Range;
 
     /// Dimension of the mesh
     static constexpr int dim = GridView::dimension;
@@ -61,9 +61,9 @@ namespace AMDiS
   public:
     /// Constructor.
     FileWriter(std::string const& baseName,
-               Vector const& dofvector)
+               DiscreteFunction const& discreteFct)
       : FileWriterInterface(baseName)
-      , dofvector_(dofvector)
+      , discreteFct_(discreteFct)
       , animation_(false)
     {
       Parameters::get(baseName + "->ParaView animation", animation_);
@@ -86,7 +86,7 @@ namespace AMDiS
       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 +105,11 @@ namespace AMDiS
   protected:
     GridView const& gridView() const
     {
-      return dofvector_.basis().gridView();
+      return discreteFct_.basis().gridView();
     }
 
   private:
-    Vector dofvector_;
+    DiscreteFunction discreteFct_;
 
     std::shared_ptr<Dune::VTKWriter<GridView>> vtkWriter_;
     std::shared_ptr<Dune::VTKSequenceWriter<GridView>> vtkSeqWriter_;
@@ -125,9 +125,9 @@ 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)
+      DiscreteFunction<GlobalBasis,Range,TreePath> const& discreteFct)
   {
-    return std::make_shared<FileWriter<GlobalBasis,Range,TreePath>>(baseName, dofvector);
+    return std::make_shared<FileWriter<GlobalBasis,Range,TreePath>>(baseName, 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/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
+   **/
+  template <class GlobalBasisType, class RangeType, class TreePathType>
+  class DiscreteFunction
+  {
+  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;
+
+    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;
+    };
+
+
+  public: // a local view on the values
+
+    /// A LocalFunction, i.e., an element local view on the DOFVector
+    class LocalFunction
+    {
+    public:
+      using Domain = typename DiscreteFunction::LocalDomain;
+      using Range = typename DiscreteFunction::Range;
+
+      enum { hasDerivative = true };
+
+    private:
+      using LocalView = typename GlobalBasis::LocalView;
+
+    public:
+      LocalFunction(DiscreteFunction const& globalFunction)
+        : globalFunction_(&globalFunction)
+        , localView_(globalFunction_->basis().localView())
+        , subTree_(&child(localView_.tree(), globalFunction_->treePath()))
+      {}
+
+      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:
+      DiscreteFunction const* globalFunction_;
+      LocalView localView_;
+      SubTree const* subTree_;
+
+      bool bound_ = false;
+    };
+
+
+  public:
+    /// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
+    DiscreteFunction(DOFVector<GlobalBasis,RangeType> const& dofVector, TreePath const& treePath)
+      : dofVector_(&dofVector)
+      , treePath_(treePath)
+      , entitySet_(dofVector.basis().gridView())
+      , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(dofVector.basis(), treePath))
+    {}
+
+    /// 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};
+    }
+
+    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 Generator for a \ref DiscreteFunction
+  template <class GlobalBasis, class RangeType, class TreePath>
+  auto makeDiscreteFunction(DOFVector<GlobalBasis, RangeType> const& dofVector, TreePath const& treePath)
+  {
+    return DiscreteFunction<GlobalBasis, RangeType, TreePath>{dofVector, treePath};
+  }
+
+  /// A Generator for a \ref DiscreteFunction
+  template <class GlobalBasis, class RangeType>
+  auto makeDiscreteFunction(DOFVector<GlobalBasis, RangeType> const& dofVector)
+  {
+    auto treePath = Dune::TypeTree::hybridTreePath();
+    return DiscreteFunction<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
+  }
+
+} // end namespace AMDiS
+
+#include "DiscreteFunction.inc.hpp"

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

@@ -3,10 +3,12 @@
 #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>::
+typename DiscreteFunction<GB,RT,TP>::Range DiscreteFunction<GB,RT,TP>::
 LocalFunction::operator()(LocalDomain const& x) const
 {
   assert( bound_ );
@@ -50,7 +52,7 @@ LocalFunction::operator()(LocalDomain const& x) const
 
 
 template <class GB, class RT, class TP>
-typename DOFVectorConstView<GB,RT,TP>::DerivativeRange DOFVectorConstView<GB,RT,TP>::
+typename DiscreteFunction<GB,RT,TP>::DerivativeRange DiscreteFunction<GB,RT,TP>::
 GradientLocalFunction::operator()(LocalDomain const& x) const
 {
   assert( bound_ );
@@ -106,4 +108,22 @@ GradientLocalFunction::operator()(LocalDomain const& x) const
   return dy;
 }
 
+
+template <class GB, class RT, class TP>
+typename DiscreteFunction<GB,RT,TP>::Range DiscreteFunction<GB,RT,TP>::
+operator()(Domain const& x) const
+{
+  using Grid = typename GB::GridView::Grid;
+  using IS = typename GB::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

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
   CMD_ARGS "${CMAKE_SOURCE_DIR}/examples/init/ellipt.dat.2d")

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 <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();
+  }
+
+  AMDiS::finalize();
+  return 0;
+}