Replaced GlobalGeodesicFEFunction by a new one GlobalGFEFunction

There are two important changes:  First of all, the new class implements
the `dune-functions` interface rather than the deprecated one based
on inheritance from `VirtualGridViewFunction`. Secondly, the new class
does not hard-wire geodesic interpolation anymore.  Rather, you can
give it interpolation classes which then govern how interpolation is done.

This change was long overdue.  It was now needed to make
gradient-flow.cc compile with the post 2.9 version of dune-fufem again.

CHANGELOG.md

 # Master
 
+- Replaced the class `GlobalGeodesicFEFunction` by a new one called
+  `GlobalGFEFunction`.  There are two important changes:  First of all,
+  the new class implements the `dune-functions` interface rather than
+  the deprecated one based on inheritance from `VirtualGridViewFunction`.
+  Secondly, the new class does not hard-wire geodesic interpolation
+  anymore.  Rather, you can give it interpolation classes which then
+  govern how interpolation is done.
+
 - Added dune-gmsh4 as a dependency
 - Build cosserat-continuum for different combinations of LFE-orders and
   GFE-orders, the respective program is called

dune/gfe/globalgeodesicfefunction.hh

-#ifndef GLOBAL_GEODESIC_FINITE_ELEMENT_FUNCTION_HH
-#define GLOBAL_GEODESIC_FINITE_ELEMENT_FUNCTION_HH
-
-#include <vector>
-
-#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
-
-#include <dune/gfe/localgeodesicfefunction.hh>
-
-template <class dctype, int DimDomain, class rctype, int DimRange>
-struct DerivativeTypefier<Dune::FieldVector<dctype, DimDomain>, UnitVector<rctype,DimRange> >
-{
-  typedef Dune::FieldMatrix<rctype, DimRange, DimDomain> DerivativeType;
-};
-
-template <class dctype, int DimDomain, class rctype, int DimRange>
-struct DerivativeTypefier<Dune::FieldVector<dctype, DimDomain>, RealTuple<rctype,DimRange> >
-{
-  typedef Dune::FieldMatrix<rctype, DimRange, DimDomain> DerivativeType;
-};
-
-
-
-/** \brief Global geodesic finite element function.
- *
- *  \tparam B  - The global basis type.
- *  \tparam TargetSpace - The manifold that this functions takes its values in.
- */
-template<class B, class TargetSpace>
-class GlobalGeodesicFEFunction
-: public VirtualGridViewFunction<typename B::GridView, TargetSpace>
-{
-
-public:
-    typedef B Basis;
-
-    typedef typename Basis::LocalFiniteElement LocalFiniteElement;
-    typedef typename Basis::GridView GridView;
-    typedef typename GridView::template Codim<0>::Entity Element;
-    typedef typename GridView::Grid::ctype  ctype;
-
-    typedef LocalGeodesicFEFunction<GridView::dimension, ctype, LocalFiniteElement, TargetSpace> LocalGFEFunction;
-    typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
-
-    //! Dimension of the grid.
-    constexpr static int gridDim = GridView::dimension;
-
-    //! Dimension of the embedded tangent space
-    constexpr static int embeddedDim = EmbeddedTangentVector::dimension;
-
-
-    //! Create global function by a global basis and the corresponding coefficient vector
-    GlobalGeodesicFEFunction(const Basis& basis, const std::vector<TargetSpace>& coefficients) :
-        VirtualGridViewFunction<typename B::GridView, TargetSpace>(basis.getGridView()),
-        basis_(basis),
-        coefficients_(coefficients)
-    {}
-
-
-    /** \brief Evaluate the function at local coordinates. */
-    void evaluateLocal(const Element& element, const Dune::FieldVector<ctype,gridDim>& local, TargetSpace& out) const
-    {
-        int numOfBaseFct = basis_.getLocalFiniteElement(element).localBasis().size();
-
-        // Extract local coefficients
-        std::vector<TargetSpace> localCoeff(numOfBaseFct);
-
-        for (int i=0; i<numOfBaseFct; i++)
-            localCoeff[i] = coefficients_[basis_.index(element,i)];
-
-        // create local gfe function
-        LocalGFEFunction localGFE(basis_.getLocalFiniteElement(element),localCoeff);
-        out = localGFE.evaluate(local);
-    }
-
-    /** \brief Evaluate the derivative of the function at local coordinates. */
-    void evaluateDerivativeLocal(const Element& element, const Dune::FieldVector<ctype,gridDim>& local,
-                                 Dune::FieldMatrix<ctype, embeddedDim, gridDim>& out) const
-    {
-        int numOfBaseFct = basis_.getLocalFiniteElement(element).localBasis().size();
-
-        // Extract local coefficients
-        std::vector<TargetSpace> localCoeff(numOfBaseFct);
-
-        for (int i=0; i<numOfBaseFct; i++)
-            localCoeff[i] = coefficients_[basis_.index(element,i)];
-
-        // create local gfe function
-        LocalGFEFunction localGFE(basis_.getLocalFiniteElement(element),localCoeff);
-
-        // use it to evaluate the derivative
-        Dune::FieldMatrix<ctype, embeddedDim, gridDim> refJac = localGFE.evaluateDerivative(local);
-
-        out =0.0;
-        //transform the gradient
-        const auto jacInvTrans = element.geometry().jacobianInverseTransposed(local);
-        for (size_t k=0; k< refJac.N(); k++)
-            jacInvTrans.umv(refJac[k],out[k]);
-
-    }
-
-    /** \brief Export basis */
-    const Basis& basis() const
-    {
-        return basis_;
-    }
-
-    /** \brief Export coefficients. */
-    const std::vector<TargetSpace>& coefficients() const
-    {
-        return coefficients_;
-    }
-
-private:
-    //! The global basis
-    const Basis& basis_;
-    //! The coefficient vector
-    const std::vector<TargetSpace>& coefficients_;
-};
-#endif

dune/gfe/globalgfefunction.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_GFE_GLOBALGFEFUNCTION_HH
+#define DUNE_GFE_GLOBALGFEFUNCTION_HH
+
+#include <memory>
+#include <optional>
+#include <vector>
+
+#include <dune/common/typetraits.hh>
+
+#include <dune/grid/utility/hierarchicsearch.hh>
+
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+#include <dune/functions/gridfunctions/gridfunction.hh>
+#include <dune/functions/backends/concepts.hh>
+
+namespace Dune::GFE {
+
+namespace Impl {
+
+/** \brief Common base class for GlobalGFEFunction and its derivative
+ *
+ * \tparam B Scalar(!) function-space basis
+ * \tparam V Container of coefficients
+ * \tparam LocalInterpolationRule How to interpolate manifold-valued data
+ */
+template<typename B, typename V, typename LocalInterpolationRule>
+class GlobalGFEFunctionBase
+{
+public:
+  using Basis = B;
+  using Vector = V;
+
+  // In order to make the cache work for proxy-references
+  // we have to use AutonomousValue<T> instead of std::decay_t<T>
+  using Coefficient = Dune::AutonomousValue<decltype(std::declval<Vector>()[std::declval<typename Basis::MultiIndex>()])>;
+
+  using GridView = typename Basis::GridView;
+  using EntitySet = Functions::GridViewEntitySet<GridView, 0>;
+  using Tree = typename Basis::LocalView::Tree;
+
+  using Domain = typename EntitySet::GlobalCoordinate;
+
+  using LocalDomain = typename EntitySet::LocalCoordinate;
+  using Element = typename EntitySet::Element;
+
+protected:
+
+  // This collects all data that is shared by all related
+  // global and local functions. This way we don't need to
+  // keep track of it individually.
+  struct Data
+  {
+    EntitySet entitySet;
+    std::shared_ptr<const Basis> basis;
+    std::shared_ptr<const Vector> coefficients;
+  };
+
+public:
+  class LocalFunctionBase
+  {
+    using LocalView = typename Basis::LocalView;
+    using size_type = typename Tree::size_type;
+
+  public:
+    using Domain = LocalDomain;
+    using Element = typename EntitySet::Element;
+
+  protected:
+    LocalFunctionBase(const std::shared_ptr<const Data>& data)
+      : data_(data)
+      , localView_(data_->basis->localView())
+    {
+      localDoFs_.reserve(localView_.maxSize());
+    }
+
+    /**
+     * \brief Copy-construct the local-function.
+     *
+     * This copy-constructor copies the cached local DOFs only
+     * if the `other` local-function is bound to an element.
+     **/
+    LocalFunctionBase(const LocalFunctionBase& other)
+      : data_(other.data_)
+      , localView_(other.localView_)
+    {
+      localDoFs_.reserve(localView_.maxSize());
+      if (bound())
+        localDoFs_ = other.localDoFs_;
+    }
+
+    /**
+     * \brief Copy-assignment of the local-function.
+     *
+     * Assign all members from `other` to `this`, except the
+     * local DOFs. Those are copied only if the `other`
+     * local-function is bound to an element.
+     **/
+    LocalFunctionBase& operator=(const LocalFunctionBase& other)
+    {
+      data_ = other.data_;
+      localView_ = other.localView_;
+      if (bound())
+        localDoFs_ = other.localDoFs_;
+      return *this;
+    }
+
+  public:
+    /**
+     * \brief Bind LocalFunction to grid element.
+     *
+     * You must call this method before `operator()`
+     * and after changes to the coefficient vector.
+     */
+    void bind(const Element& element)
+    {
+      localView_.bind(element);
+
+      localDoFs_.resize(localView_.size());
+      const auto& dofs = *data_->coefficients;
+      for (size_type i = 0; i < localView_.tree().size(); ++i)
+      {
+        // For a subspace basis the index-within-tree i
+        // is not the same as the localIndex within the
+        // full local view.
+        size_t localIndex = localView_.tree().localIndex(i);
+        localDoFs_[localIndex] = dofs[localView_.index(localIndex)];
+      }
+
+      // create local GFE function
+      // TODO Store this object by value
+      localInterpolationRule_ = std::make_unique<LocalInterpolationRule>(this->localView_.tree().finiteElement(),localDoFs_);
+    }
+
+    //! Unbind the local-function.
+    void unbind()
+    {
+      localView_.unbind();
+    }
+
+    //! Check if LocalFunction is already bound to an element.
+    bool bound() const
+    {
+      return localView_.bound();
+    }
+
+    //! Return the element the local-function is bound to.
+    const Element& localContext() const
+    {
+      return localView_.element();
+    }
+
+  protected:
+
+    std::shared_ptr<const Data> data_;
+    LocalView localView_;
+    std::vector<Coefficient> localDoFs_;
+    std::unique_ptr<LocalInterpolationRule> localInterpolationRule_;
+  };
+
+protected:
+  GlobalGFEFunctionBase(const std::shared_ptr<const Data>& data)
+    : data_(data)
+  {
+    /* Nothing. */
+  }
+
+public:
+
+  //! Return a const reference to the stored basis.
+  const Basis& basis() const
+  {
+    return *data_->basis;
+  }
+
+  //! Return the coefficients of this discrete function by reference.
+  const Vector& dofs() const
+  {
+    return *data_->coefficients;
+  }
+
+  //! Get associated set of entities the local-function can be bound to.
+  const EntitySet& entitySet() const
+  {
+    return data_->entitySet;
+  }
+
+protected:
+  std::shared_ptr<const Data> data_;
+};
+
+} // namespace Impl
+
+
+
+template<typename GGF>
+class GlobalGFEFunctionDerivative;
+
+/**
+ * \brief A global geometric finite element function
+ *
+ * \tparam B Type of global basis
+ * \tparam LIR Local interpolation rule for manifold-valued data
+ * \tparam TargetSpace Range type of this function
+ */
+template<typename B, typename LIR, typename TargetSpace>
+class GlobalGFEFunction
+  : public Impl::GlobalGFEFunctionBase<B, std::vector<TargetSpace>, LIR>
+{
+  using Base = Impl::GlobalGFEFunctionBase<B, std::vector<TargetSpace>, LIR>;
+  using Data = typename Base::Data;
+
+public:
+  using Basis = typename Base::Basis;
+  using Vector = typename Base::Vector;
+  using LocalInterpolationRule = LIR;
+
+  using Domain = typename Base::Domain;
+  using Range = typename TargetSpace::CoordinateType;
+
+  using Traits = Functions::Imp::GridFunctionTraits<Range(Domain), typename Base::EntitySet, Functions::DefaultDerivativeTraits, 16>;
+
+  class LocalFunction
+    : public Base::LocalFunctionBase
+  {
+    using LocalBase = typename Base::LocalFunctionBase;
+    using size_type = typename Base::Tree::size_type;
+
+  public:
+
+    using GlobalFunction = GlobalGFEFunction;
+    using Domain = typename LocalBase::Domain;
+    using Range = GlobalFunction::Range;
+    using Element = typename LocalBase::Element;
+
+    //! Create a local-function from the associated grid-function
+    LocalFunction(const GlobalGFEFunction& globalFunction)
+      : LocalBase(globalFunction.data_)
+    {
+      /* Nothing. */
+    }
+
+    /**
+     * \brief Evaluate this local-function in coordinates `x` in the 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.
+     */
+    Range operator()(const Domain& x) const
+    {
+      return this->localInterpolationRule_->evaluate(x).globalCoordinates();
+    }
+
+    //! Local function of the derivative
+    friend typename GlobalGFEFunctionDerivative<GlobalGFEFunction>::LocalFunction derivative(const LocalFunction& lf)
+    {
+      auto dlf = localFunction(GlobalGFEFunctionDerivative<GlobalGFEFunction>(lf.data_));
+      if (lf.bound())
+        dlf.bind(lf.localContext());
+      return dlf;
+    }
+  };
+
+  //! Create a grid-function, by wrapping the arguments in `std::shared_ptr`.
+  template<class B_T, class V_T>
+  GlobalGFEFunction(B_T && basis, V_T && coefficients)
+    : Base(std::make_shared<Data>(Data{{basis.gridView()}, wrap_or_move(std::forward<B_T>(basis)), wrap_or_move(std::forward<V_T>(coefficients))}))
+  {}
+
+  //! Create a grid-function, by moving the arguments in `std::shared_ptr`.
+  GlobalGFEFunction(std::shared_ptr<const Basis> basis, std::shared_ptr<const Vector> coefficients)
+    : Base(std::make_shared<Data>(Data{{basis->gridView()}, basis, coefficients}))
+  {}
+
+  /** \brief Evaluate at a point given in world coordinates
+   *
+   * \warning This has to find the element that the evaluation point is in.
+   *   It is therefore very slow.
+   */
+  Range operator() (const Domain& x) const
+  {
+    HierarchicSearch search(this->data_->basis->gridView().grid(), this->data_->basis->gridView().indexSet());
+
+    const auto e = search.findEntity(x);
+    auto localThis = localFunction(*this);
+    localThis.bind(e);
+    return localThis(e.geometry().local(x));
+  }
+
+  //! Derivative of the `GlobalGFEFunction`
+  friend GlobalGFEFunctionDerivative<GlobalGFEFunction> derivative(const GlobalGFEFunction& f)
+  {
+    return GlobalGFEFunctionDerivative<GlobalGFEFunction>(f.data_);
+  }
+
+  /**
+   * \brief Construct local function from a GlobalGFEFunction.
+   *
+   * The obtained local function satisfies the concept
+   * `Dune::Functions::Concept::LocalFunction`. It must be bound
+   * to an entity from the entity set of the GlobalGFEFunction
+   * before it can be used.
+   */
+  friend LocalFunction localFunction(const GlobalGFEFunction& t)
+  {
+    return LocalFunction(t);
+  }
+};
+
+
+/**
+ * \brief Derivative of a `GlobalGFEFunction`
+ *
+ * Function returning the derivative of the given `GlobalGFEFunction`
+ * with respect to global coordinates.
+ *
+ * \tparam GGF instance of the `GlobalGFEFunction` this is a derivative of
+ */
+template<typename GGF>
+class GlobalGFEFunctionDerivative
+  : public Impl::GlobalGFEFunctionBase<typename GGF::Basis, typename GGF::Vector, typename GGF::LocalInterpolationRule>
+{
+  using Base = Impl::GlobalGFEFunctionBase<typename GGF::Basis, typename GGF::Vector, typename GGF::LocalInterpolationRule>;
+  using Data = typename Base::Data;
+
+public:
+  using GlobalGFEFunction = GGF;
+
+  using Basis = typename Base::Basis;
+  using Vector = typename Base::Vector;
+
+  using Domain = typename Base::Domain;
+  using Range = typename Functions::SignatureTraits<typename GlobalGFEFunction::Traits::DerivativeInterface>::Range;
+
+  using Traits = Functions::Imp::GridFunctionTraits<Range(Domain), typename Base::EntitySet, Functions::DefaultDerivativeTraits, 16>;
+
+  /**
+   * \brief local function evaluating the derivative in reference coordinates
+   *
+   * Note that the function returns the derivative with respect to global
+   * coordinates even when the point is given in reference coordinates on
+   * an element.
+   */
+  class LocalFunction
+    : public Base::LocalFunctionBase
+  {
+    using LocalBase = typename Base::LocalFunctionBase;
+    using size_type = typename Base::Tree::size_type;
+
+  public:
+    using GlobalFunction = GlobalGFEFunctionDerivative;
+    using Domain = typename LocalBase::Domain;
+    using Range = GlobalFunction::Range;
+    using Element = typename LocalBase::Element;
+
+    //! Create a local function from the associated grid function
+    LocalFunction(const GlobalFunction& globalFunction)
+      : LocalBase(globalFunction.data_)
+    {
+      /* Nothing. */
+    }
+
+    /**
+     * \brief Bind LocalFunction to grid element.
+     *
+     * You must call this method before `operator()`
+     * and after changes to the coefficient vector.
+     */
+    void bind(const Element& element)
+    {
+      LocalBase::bind(element);
+      geometry_.emplace(element.geometry());
+    }
+
+    //! Unbind the local-function.
+    void unbind()
+    {
+      geometry_.reset();
+      LocalBase::unbind();
+    }
+
+    /**
+     * \brief Evaluate this local-function in coordinates `x` in the 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.
+     *
+     * Note that the function returns the derivative with respect to global
+     * coordinates even though the evaluation point is given in reference coordinates
+     * on the current element.
+     */
+    Range operator()(const Domain& x) const
+    {
+        // Jacobian with respect to local coordinates
+        auto refJac = this->localInterpolationRule_->evaluateDerivative(x);
+
+        // Transform to world coordinates
+        return refJac * geometry_->jacobianInverse(x);
+    }
+
+    //! Not implemented
+    friend typename Traits::LocalFunctionTraits::DerivativeInterface derivative(const LocalFunction&)
+    {
+      DUNE_THROW(NotImplemented, "derivative of derivative is not implemented");
+    }
+
+  private:
+    std::optional<typename Element::Geometry> geometry_;
+  };
+
+  /**
+   * \brief create object from `GlobalGFEFunction` data
+   *
+   * Please call `derivative(globalGFEFunction)` to create an instance
+   * of this class.
+   */
+  GlobalGFEFunctionDerivative(const std::shared_ptr<const Data>& data)
+    : Base(data)
+  {
+    /* Nothing. */
+  }
+
+  /** \brief Evaluate the discrete grid-function derivative in global coordinates
+   *
+   * \warning This has to find the element that the evaluation point is in.
+   *   It is therefore very slow.
+   */
+  Range operator()(const Domain& x) const
+  {
+    HierarchicSearch search(this->data_->basis->gridView().grid(), this->data_->basis->gridView().indexSet());
+
+    const auto e = search.findEntity(x);
+    auto localThis = localFunction(*this);
+    localThis.bind(e);
+    return localThis(e.geometry().local(x));
+  }
+
+  friend typename Traits::DerivativeInterface derivative(const GlobalGFEFunctionDerivative& f)
+  {
+    DUNE_THROW(NotImplemented, "derivative of derivative is not implemented");
+  }
+
+  //! Construct local function from a `GlobalGFEFunctionDerivative`
+  friend LocalFunction localFunction(const GlobalGFEFunctionDerivative& f)
+  {
+    return LocalFunction(f);
+  }
+};
+
+} // namespace Dune::GFE
+
+#endif // DUNE_GFE_GLOBALGFEFUNCTION_HH

dune/gfe/l2distancesquaredenergy.hh

@@ -5,6 +5,7 @@
 
 #include <dune/geometry/quadraturerules.hh>
 
+#include <dune/gfe/globalgfefunction.hh>
 #include <dune/gfe/localenergy.hh>
 #include <dune/gfe/localgeodesicfefunction.hh>
 
@@ -20,10 +21,12 @@ class L2DistanceSquaredEnergy
   // some other sizes
   constexpr static int gridDim = GridView::dimension;
 
+  using LocalInterpolationRule = LocalGeodesicFEFunction<gridDim, DT, typename Basis::LocalView::Tree::FiniteElement, typename TargetSpace::template rebind<double>::other>;
+
 public:
 
   // This is the function that we are computing the L2-distance to
-  std::shared_ptr<VirtualGridViewFunction<GridView,typename TargetSpace::template rebind<double>::other > > origin_;
+  std::shared_ptr<Dune::GFE::GlobalGFEFunction<Basis, LocalInterpolationRule, typename TargetSpace::template rebind<double>::other > > origin_;
 
   /** \brief Assemble the energy for a single element */
   RT energy (const typename Basis::LocalView& localView,
@@ -36,11 +39,13 @@ public:
     typedef LocalGeodesicFEFunction<gridDim, double, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
     LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
 
+    const auto element = localView.element();
+    auto localOrigin = localFunction(*origin_);
+    localOrigin.bind(element);
+
     // Just guessing an appropriate quadrature order
     auto quadOrder = localFiniteElement.localBasis().order() * 2 * gridDim;
 
-    const auto element = localView.element();
-
     const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
 
     for (size_t pt=0; pt<quad.size(); pt++)
@@ -54,21 +59,14 @@ public:
 
       // The function value
       auto value = localGeodesicFEFunction.evaluate(quadPos);
-      typename TargetSpace::template rebind<double>::other originValue;
-      origin_->evaluateLocal(element,quadPos, originValue);
+      auto originValue = localOrigin(quadPos);
 
       // The derivative of the 'origin' function
       // First: as function defined on the reference element
-      typename VirtualGridViewFunction<GridView,typename TargetSpace::template rebind<double>::other>::DerivativeType originReferenceDerivative;
-      origin_->evaluateDerivativeLocal(element,quadPos,originReferenceDerivative);
+      auto originReferenceDerivative = derivative(localOrigin)(quadPos);
 
       // The derivative of the function defined on the actual element
-      typename VirtualGridViewFunction<GridView,typename TargetSpace::template rebind<double>::other >::DerivativeType originDerivative(0);
-
-      auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
-
-      for (size_t comp=0; comp<originReferenceDerivative.N(); comp++)
-        jacobianInverseTransposed.umv(originReferenceDerivative[comp], originDerivative[comp]);
+      auto originDerivative = originReferenceDerivative * element.geometry().jacobianInverse(quadPos);
 
       double weightFactor = originDerivative.frobenius_norm();
       // un-comment the following line to switch off the weight factor

src/gradient-flow.cc

@@ -25,11 +25,10 @@
 #include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
 #include <dune/functions/functionspacebases/interpolate.hh>
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
 
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/functiontools/boundarydofs.hh>
-#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
-#include <dune/fufem/discretizationerror.hh>
 #include <dune/fufem/dunepython.hh>
 
 #include <dune/solvers/solvers/iterativesolver.hh>
@@ -39,7 +38,7 @@
 #include <dune/gfe/localgeodesicfeadolcstiffness.hh>
 #include <dune/gfe/geodesicfeassembler.hh>
 #include <dune/gfe/riemanniantrsolver.hh>
-#include <dune/gfe/globalgeodesicfefunction.hh>
+#include <dune/gfe/globalgfefunction.hh>
 #include <dune/gfe/embeddedglobalgfefunction.hh>
 #include <dune/gfe/harmonicenergy.hh>
 #include <dune/gfe/l2distancesquaredenergy.hh>
@@ -137,22 +136,19 @@ int main (int argc, char *argv[]) try
   }
 
   grid->globalRefine(numLevels-1);
+  auto gridView = grid->leafGridView();
 
   //////////////////////////////////////////////////////////////////////////////////
   //  Construct the scalar function space basis corresponding to the GFE space
   //////////////////////////////////////////////////////////////////////////////////
 
   typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, order> FEBasis;
-  //typedef Dune::Functions::Periodic1DPQ1NodalBasis<typename GridType::LeafGridView> FEBasis;
   FEBasis feBasis(grid->leafGridView());
 
   ///////////////////////////////////////////
   //   Read Dirichlet values
   ///////////////////////////////////////////
 
-  typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
-  FufemFEBasis fufemFeBasis(feBasis);
-
   BitSetVector<1> dirichletVertices(grid->leafGridView().size(dim), false);
 
   const auto& indexSet = grid->leafGridView().indexSet();
@@ -171,7 +167,7 @@ int main (int argc, char *argv[]) try
   BoundaryPatch<GridType::LeafGridView> dirichletBoundary(grid->leafGridView(), dirichletVertices);
 
   BitSetVector<blocksize> dirichletNodes(feBasis.size(), false);
-  constructBoundaryDofs(dirichletBoundary,fufemFeBasis,dirichletNodes);
+  constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
 
   ////////////////////////////
   //   Initial iterate
@@ -182,7 +178,16 @@ int main (int argc, char *argv[]) try
   auto pythonInitialIterate = Python::make_function<TargetSpace::CoordinateType>(module.get("f"));
 
   std::vector<TargetSpace::CoordinateType> v;
-  Functions::interpolate(feBasis, v, pythonInitialIterate);
+  using namespace Functions::BasisFactory;
+
+  auto powerBasis = makeBasis(
+    gridView,
+    power<TargetSpace::CoordinateType::dimension>(
+      lagrange<order>(),
+      blockedInterleaved()
+  ));
+
+  Functions::interpolate(powerBasis, v, pythonInitialIterate);
 
   SolutionType x(feBasis.size());
 
@@ -259,7 +264,7 @@ int main (int argc, char *argv[]) try
 
   for (int i=0; i<numTimeSteps; i++)
   {
-    auto previousTimeStepFct = std::make_shared<GlobalGeodesicFEFunction<FufemFEBasis,TargetSpace> >(fufemFeBasis,previousTimeStep);
+    auto previousTimeStepFct = std::make_shared<GFE::GlobalGFEFunction<FEBasis,GeodesicInterpolationRule::rebind<TargetSpace>::other,TargetSpace> >(feBasis,previousTimeStep);
     l2DistanceSquaredEnergy->origin_ = previousTimeStepFct;
 
     solver.setInitialIterate(x);