diff --git a/CHANGELOG.md b/CHANGELOG.md
index e17db3871db4fd6dfde5bed5b34c134f5461cabb..663fa3503d2277a580c5b21362f6362585158746 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,7 +1,8 @@
 # Master
 
-- The class `LocalGeodesicFEFunction` now takes a `dune-functions` basis
-  as its first argument, instead of a `LocalFiniteElement`.
+- The classes `LocalGeodesicFEFunction` and `LocalProjectedFEFunction`
+  now take a `dune-functions` basis as their first arguments, instead
+  of a `LocalFiniteElement`.
 
 - The entire code is now in the `Dune::GFE` namespace.
 
diff --git a/dune/gfe/assemblers/simofoxenergy.hh b/dune/gfe/assemblers/simofoxenergy.hh
index a59fc9e445ab2fccb21b1bf41d40d1d14d01c616..424376beb1d0090b037dc20456f4fe5c6845b744 100644
--- a/dune/gfe/assemblers/simofoxenergy.hh
+++ b/dune/gfe/assemblers/simofoxenergy.hh
@@ -5,6 +5,7 @@
 #include <dune/common/parametertree.hh>
 #include <dune/common/transpose.hh>
 #include <dune/common/tuplevector.hh>
+#include <dune/functions/functionspacebases/subspacebasis.hh>
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/geometry/quadraturerules.hh>
 #include <dune/gfe/linearalgebra.hh>
@@ -48,7 +49,7 @@ namespace Dune::GFE
    * \tparam LocalFEFunction Decides which interpolation function is used for the interpolation of the midsurface
    * and director field.
    */
-  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type = double>
+  template <class Basis, template <typename, typename> typename LocalFEFunction, typename field_type = double>
   class SimoFoxEnergyLocalStiffness
     : public LocalEnergy<Basis, ProductManifold<RealTuple<field_type, 3>,
           UnitVector<field_type, 3> > >
@@ -65,22 +66,6 @@ namespace Dune::GFE
     static constexpr int gridDim  = GridView::dimension;
     static constexpr int dimworld = GridView::dimensionworld;
 
-    // The local finite element type used for midsurface position and displacement interpolation
-    using MidSurfaceElement =
-      typename std::result_of<decltype (&LocalFiniteElementFactory<Basis, 0>::get)(typename Basis::LocalView, decltype(Dune::Indices::_0))>::type;
-    // The local finite element type used for the director interpolation
-    using DirectorElement =
-      typename std::result_of<decltype (&LocalFiniteElementFactory<Basis, 1>::get)(typename Basis::LocalView, decltype(Dune::Indices::_1))>::type;
-
-    // The local finite element function type to evaluate the midsurface position
-    using LocalMidSurfaceFunctionType = LocalFEFunction<gridDim, DT, MidSurfaceElement, RealTuple<field_type, 3> >;
-    // The local finite element function type to evaluate the director
-    using LocalDirectorFunctionType = LocalFEFunction<gridDim, DT, DirectorElement, UnitVector<field_type, 3> >;
-
-    // Extra function type for reference quantities since they are unconditionally doubles, i.e. no ADOL-C types
-    using LocalMidSurfaceReferenceFunctionType = LocalFEFunction<gridDim, DT, MidSurfaceElement, RealTuple<double, 3> >;
-    using LocalDirectorReferenceFunctionType   = LocalFEFunction<gridDim, DT, DirectorElement, UnitVector<double, 3> >;
-
   public:
     /** \brief Constructor with a set of material parameters
      * \param parameters The material parameters
@@ -195,7 +180,7 @@ namespace Dune::GFE
    * b is similar to the second fundamental form of the surface, but the unit normal is replaced with the shell director
    * gamma_1 and gamma_2 are the transverse shear in the two parametric directions
    * Paper Equation 4.10 */
-  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
+  template <class Basis, template <typename, typename> typename LocalFEFunction, typename field_type>
   Dune::FieldVector<field_type, 8> SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::calculateGreenLagrangianStrains(
     const KinematicVariables &kin)
   {
@@ -222,7 +207,7 @@ namespace Dune::GFE
    *
    * Paper Equations 6.4 - 6.8
    */
-  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
+  template <class Basis, template <typename, typename> typename LocalFEFunction, typename field_type>
   template <typename Element, typename LocalDirectorFunction, typename LocalMidSurfaceFunction, typename LocalDirectorReferenceFunction,
       typename LocalMidSurfaceReferenceFunction, typename IntegrationPointPosition>
   auto SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::kinematicVariablesFactory(
@@ -245,7 +230,7 @@ namespace Dune::GFE
     return kin;
   }
 
-  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
+  template <class Basis, template <typename, typename> typename LocalFEFunction, typename field_type>
   auto SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::getReferenceLocalConfigurations(
     const typename Basis::LocalView &localView) const {
     using namespace Dune::Indices;
@@ -283,7 +268,7 @@ namespace Dune::GFE
    *  E are the components of the Green-Lagrangian strains [membrane strains, bending, transverse shear]
    *  see for details Paper Equation 4.11,4.10 and 10.1
    */
-  template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
+  template <class Basis, template <typename, typename> typename LocalFEFunction, typename field_type>
   typename SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::RT
   SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::energy(
     const typename Basis::LocalView &localView,
@@ -305,6 +290,19 @@ namespace Dune::GFE
     for (size_t i = 0; i < localMidSurfaceConfiguration.size(); ++i)
       displacements.emplace_back(localMidSurfaceConfiguration[i].globalCoordinates() - localRefMidSurfaceConfiguration[i].globalCoordinates());
 
+    // The local finite element type used for midsurface position and displacement interpolation
+    const auto midSurfaceBasis = Functions::subspaceBasis(localView.globalBasis(),_0,0);
+    const auto directorBasis = Functions::subspaceBasis(localView.globalBasis(),_1,0);
+
+    // The local finite element function type to evaluate the midsurface position
+    using LocalMidSurfaceFunctionType = LocalFEFunction<decltype(midSurfaceBasis), RealTuple<field_type, 3> >;
+    // The local finite element function type to evaluate the director
+    using LocalDirectorFunctionType = LocalFEFunction<decltype(directorBasis), UnitVector<field_type, 3> >;
+
+    // Extra function type for reference quantities since they are unconditionally doubles, i.e. no ADOL-C types
+    using LocalMidSurfaceReferenceFunctionType = LocalFEFunction<decltype(midSurfaceBasis), RealTuple<double, 3> >;
+    using LocalDirectorReferenceFunctionType   = LocalFEFunction<decltype(directorBasis), UnitVector<double, 3> >;
+
     const LocalMidSurfaceFunctionType localMidSurfaceFunction(midSurfaceElement, localMidSurfaceConfiguration);
     const LocalMidSurfaceFunctionType localMidSurfaceDisplacementFunction(midSurfaceElement, displacements);
     const LocalMidSurfaceReferenceFunctionType localMidSurfaceReferenceFunction(midSurfaceElement, localRefMidSurfaceConfiguration);
diff --git a/dune/gfe/bendingisometryhelper.hh b/dune/gfe/bendingisometryhelper.hh
index 8db7bf173becd76960361798c0c037329eb91de5..6fcc20ab15d2f94cd732e40c9a697d2143b2cd31 100644
--- a/dune/gfe/bendingisometryhelper.hh
+++ b/dune/gfe/bendingisometryhelper.hh
@@ -200,7 +200,6 @@ namespace Dune::GFE::Impl
     auto localView = discreteKirchhoffBasis.localView();
     auto localViewCoefficients = coefficientBasis.localView();
 
-    using GridView = typename DiscreteKirchhoffBasis::GridView;
     auto gridView = discreteKirchhoffBasis.gridView();
 
     /** Create an EmbeddedGlobalGFEFunction that serves as a GridViewFunction later on.
@@ -211,7 +210,7 @@ namespace Dune::GFE::Impl
         where the first three entries correspond to the deformation and the
         last four entries correspond to a (quaternion) rotation.
      */
-    typedef Dune::GFE::LocalProjectedFEFunction<GridView::dimension, double, typename CoefficientBasis::LocalView::Tree::FiniteElement, Dune::GFE::ProductManifold<RealTuple<double,3>, Rotation<double,3> > > LocalInterpolationRule;
+    typedef GFE::LocalProjectedFEFunction<CoefficientBasis, GFE::ProductManifold<RealTuple<double,3>, Rotation<double,3> > > LocalInterpolationRule;
     Dune::GFE::EmbeddedGlobalGFEFunction<CoefficientBasis, LocalInterpolationRule, Dune::GFE::ProductManifold<RealTuple<double,3>, Rotation<double,3> > > embeddedGlobalFunction(coefficientBasis, isometryCoefficients);
 
     auto productSpaceGridViewFunction = Dune::Functions::makeAnalyticGridViewFunction(embeddedGlobalFunction, gridView);
diff --git a/dune/gfe/functions/discretekirchhoffbendingisometry.hh b/dune/gfe/functions/discretekirchhoffbendingisometry.hh
index 98e022be979cb93b344bd35d94810cf605846fdf..25fcceb71e90935edf17157d329782d3f3638ef0 100644
--- a/dune/gfe/functions/discretekirchhoffbendingisometry.hh
+++ b/dune/gfe/functions/discretekirchhoffbendingisometry.hh
@@ -44,7 +44,7 @@ namespace Dune::GFE
 
     static constexpr int components_ = 3;
 
-    typedef Dune::GFE::LocalProjectedFEFunction<GridView::dimension, ctype, typename CoefficientBasis::LocalView::Tree::FiniteElement, Dune::GFE::ProductManifold<RealTuple<RT,components_>, Rotation<RT,components_> > > LocalInterpolationRule;
+    typedef GFE::LocalProjectedFEFunction<CoefficientBasis, GFE::ProductManifold<RealTuple<RT,components_>, Rotation<RT,components_> > > LocalInterpolationRule;
 
     /** \brief The type used for derivatives */
     typedef Dune::FieldMatrix<RT, components_, gridDim> DerivativeType;
diff --git a/dune/gfe/functions/interpolationderivatives.hh b/dune/gfe/functions/interpolationderivatives.hh
index 646e9c7a6ccbc608829a2f7e8ae5e0f7a429f6eb..4a0cd5d5defa87d97be17fb9ac1e3cabc9c3ad05 100644
--- a/dune/gfe/functions/interpolationderivatives.hh
+++ b/dune/gfe/functions/interpolationderivatives.hh
@@ -576,11 +576,16 @@ namespace Dune::GFE
    * standard FE interpolation, and the derivatives with respect to the coefficients can be
    * computed much simpler and faster than for the general case.
    */
-  template <int gridDim, typename field_type, typename LocalFiniteElement,int dim>
-  class InterpolationDerivatives<LocalProjectedFEFunction<gridDim, field_type, LocalFiniteElement, RealTuple<field_type,dim> > >
+  template <typename Basis, typename field_type, int dim>
+  class InterpolationDerivatives<LocalProjectedFEFunction<Basis, RealTuple<field_type,dim> > >
   {
     // TODO: The implementation here would be identical to the geodesic FE case
-    using LocalInterpolationRule = LocalProjectedFEFunction<gridDim, field_type, LocalFiniteElement, RealTuple<field_type,dim> >;
+    using LocalInterpolationRule = LocalProjectedFEFunction<Basis, RealTuple<field_type,dim> >;
+
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto gridDim = LocalCoordinate::size();
+
     using TargetSpace = typename LocalInterpolationRule::TargetSpace;
 
     constexpr static auto blocksize = TargetSpace::TangentVector::dimension;
@@ -626,10 +631,9 @@ namespace Dune::GFE
      *  \param[out] derivative      The derivative of the interpolation function
      *                              with respect to the evaluation point
      */
-    template <typename Element>
     void bind(short tapeNumber,
               const Element& element,
-              const typename Element::Geometry::LocalCoordinate& localPos,
+              const LocalCoordinate& localPos,
               typename TargetSpace::CoordinateType& valueGlobalCoordinates,
               typename LocalInterpolationRule::DerivativeType& derivative)
     {
@@ -692,7 +696,7 @@ namespace Dune::GFE
       // First derivatives of the function gradient wrt to the FE coefficients
       for (size_t i=0; i<nDofs; ++i)
         for (int j=0; j<blocksize; ++j)
-          for (int k=0; k<gridDim; ++k)
+          for (std::size_t k=0; k<gridDim; ++k)
             firstDerivative[blocksize + j*gridDim + k][i*blocksize+j] = shapeFunctionGradients_[i][0][k];
 
       // For RealTuple, firstDerivative and embeddedFirstDerivative coincide
@@ -712,10 +716,19 @@ namespace Dune::GFE
    * interpolation.  No matter what the target space is, the interpolation is always
    * Euclidean in the surrounding space.
    */
-  template <int gridDim, typename ctype, typename LocalFiniteElement, typename TargetSpace>
-  class InterpolationDerivatives<LocalProjectedFEFunction<gridDim, ctype, LocalFiniteElement, TargetSpace, false> >
+  template <typename Basis, typename TargetSpace>
+  class InterpolationDerivatives<LocalProjectedFEFunction<Basis, TargetSpace, false> >
   {
-    using LocalInterpolationRule = LocalProjectedFEFunction<gridDim, ctype, LocalFiniteElement, TargetSpace, false>;
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto gridDim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using FERangeType = typename LocalBasis::Traits::RangeType;
+    using FEJacobianType = typename LocalBasis::Traits::JacobianType;
+
+    using LocalInterpolationRule = LocalProjectedFEFunction<Basis, TargetSpace, false>;
     using CoordinateType = typename TargetSpace::CoordinateType;
 
     constexpr static auto blocksize = TargetSpace::TangentVector::dimension;
@@ -735,11 +748,10 @@ namespace Dune::GFE
     const bool doDerivative_;
 
     // Values of all scalar shape functions at the point we are bound to
-    std::vector<FieldVector<double,1> > shapeFunctionValues_;
+    std::vector<FERangeType> shapeFunctionValues_;
 
     // Gradients of all scalar shape functions at the point we are bound to
-    // TODO: The second dimension must be WorldDim
-    std::vector<FieldMatrix<double,1,gridDim> > shapeFunctionGradients_;
+    std::vector<FEJacobianType> shapeFunctionGradients_;
 
     // TODO: Don't hardcode FieldMatrix
     std::vector<FieldMatrix<double,blocksize,embeddedBlocksize> > orthonormalFrames_;
@@ -771,10 +783,9 @@ namespace Dune::GFE
      *  \param[out] derivative      The derivative of the interpolation function
      *                              with respect to the evaluation point
      */
-    template <typename Element>
     void bind(short tapeNumber,
               const Element& element,
-              const typename Element::Geometry::LocalCoordinate& localPos,
+              const LocalCoordinate& localPos,
               typename TargetSpace::CoordinateType& valueGlobalCoordinates,
               typename LocalInterpolationRule::DerivativeType& derivative)
     {
@@ -841,7 +852,7 @@ namespace Dune::GFE
       // First derivatives of the function gradient wrt to the FE coefficients
       for (size_t i=0; i<nDofs; ++i)
         for (int j=0; j<embeddedBlocksize; ++j)
-          for (int k=0; k<gridDim; ++k)
+          for (std::size_t k=0; k<gridDim; ++k)
             partialDerivative[embeddedBlocksize + j*gridDim + k][i*embeddedBlocksize+j] = shapeFunctionGradients_[i][0][k];
 
       // Euclidean derivative: Derivatives in the direction of the projected canonical vectors
@@ -1041,16 +1052,25 @@ namespace Dune::GFE
    * by hand with reasonable effort.  The corresponding implementation is faster than the one
    * based on ADOL-C.
    */
-  template <int gridDim, typename field_type, typename LocalFiniteElement,int dim>
-  class InterpolationDerivatives<LocalProjectedFEFunction<gridDim, field_type, LocalFiniteElement, UnitVector<field_type,dim> > >
+  template <typename Basis, typename field_type,int dim>
+  class InterpolationDerivatives<LocalProjectedFEFunction<Basis, UnitVector<field_type,dim> > >
   {
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto gridDim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using FERangeType = typename LocalBasis::Traits::RangeType;
+    using FEJacobianType = typename LocalBasis::Traits::JacobianType;
+
     using TargetSpace = UnitVector<field_type,dim>;
     using CoordinateType = typename TargetSpace::CoordinateType;
 
     constexpr static auto blocksize = TargetSpace::TangentVector::dimension;
     constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
 
-    using LocalInterpolationRule = LocalProjectedFEFunction<gridDim, field_type, LocalFiniteElement, TargetSpace>;
+    using LocalInterpolationRule = LocalProjectedFEFunction<Basis, TargetSpace>;
 
     /** \brief A vector that can be contracted with a derivative of the projection onto the unit sphere
      *
@@ -1234,11 +1254,10 @@ namespace Dune::GFE
     const bool doDerivative_;
 
     // Values of all scalar shape functions at the point we are bound to
-    std::vector<FieldVector<double,1> > shapeFunctionValues_;
+    std::vector<FERangeType> shapeFunctionValues_;
 
     // Gradients of all scalar shape functions at the point we are bound to
-    // TODO: The second dimension must be WorldDim
-    std::vector<FieldMatrix<double,1,gridDim> > shapeFunctionGradients_;
+    std::vector<FEJacobianType> shapeFunctionGradients_;
 
     // TODO: Is this needed at all?
     typename LocalInterpolationRule::DerivativeType euclideanDerivative_;
@@ -1291,10 +1310,9 @@ namespace Dune::GFE
      *  \param[out] derivative      The derivative of the interpolation function
      *                              with respect to the evaluation point
      */
-    template <typename Element>
     void bind(short tapeNumber,
               const Element& element,
-              const typename Element::Geometry::LocalCoordinate& localPos,
+              const LocalCoordinate& localPos,
               typename TargetSpace::CoordinateType& valueGlobalCoordinates,
               typename LocalInterpolationRule::DerivativeType& derivative)
     {
@@ -1486,7 +1504,7 @@ namespace Dune::GFE
 
               if (doDerivative_)
               {
-                for (int beta=0; beta<gridDim; beta++)
+                for (std::size_t beta=0; beta<gridDim; beta++)
                 {
                   VectorToContractWith derivativeWeightsBeta(derivativeWeightsTransposed[beta]);
                   derivativeWeightsBeta.bind(projection_.x_);
diff --git a/dune/gfe/functions/localprojectedfefunction.hh b/dune/gfe/functions/localprojectedfefunction.hh
index 0dbe01fe5f0600ac39ea4add0e8ddd93d9a73788..f17e3229923a8f126e14eb6feb4c1c7362cf9c23 100644
--- a/dune/gfe/functions/localprojectedfefunction.hh
+++ b/dune/gfe/functions/localprojectedfefunction.hh
@@ -18,19 +18,27 @@ namespace Dune::GFE
 
   /** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold
    *
-   * \tparam dim Dimension of the reference element
-   * \tparam ctype Type used for coordinates on the reference element
-   * \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
-   * \tparam TargetSpace The manifold that the function takes its values in
+   * \tparam Basis The basis used for the interpolation in the embedding space.
+   * \tparam TS The manifold that the function takes its values in
    * \tparam conforming Geometrical conformity of the functions (omits projections when set to false)
+   *
+   * \note Even though the embedding space is typically more than one-dimensional,
+   * we allow this basis to be scalar-valued.
    */
-  template <int dim, class ctype, class LocalFiniteElement, class TS, bool conforming=true>
+  template <class Basis, class TS, bool conforming=true>
   class LocalProjectedFEFunction
   {
   public:
     using TargetSpace=TS;
   private:
-    using LocalCoordinate = FieldVector<ctype,dim>;
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto dim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using ScalarFERange = typename LocalBasis::Traits::RangeType;
+    using ScalarFEJacobian = typename LocalBasis::Traits::JacobianType;
 
     typedef typename TargetSpace::ctype RT;
 
@@ -58,7 +66,7 @@ namespace Dune::GFE
     template<class U>
     struct rebind
     {
-      using other = LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,U,conforming>;
+      using other = LocalProjectedFEFunction<Basis,U,conforming>;
     };
 
     /** \brief The number of Lagrange points */
@@ -129,18 +137,19 @@ namespace Dune::GFE
 
   };
 
-  template <int dim, class ctype, class LocalFiniteElement, class TargetSpace, bool conforming>
-  auto LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
+  template <class Basis, class TargetSpace, bool conforming>
+  auto LocalProjectedFEFunction<Basis,TargetSpace,conforming>::
   evaluate(const LocalCoordinate& local) const
   {
-    // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
-    std::vector<Dune::FieldVector<ctype,1> > w;
+    // Compute value of FE function in embedding space
+    std::vector<ScalarFERange> w;
     localFiniteElement_.localBasis().evaluateFunction(local,w);
 
     typename TargetSpace::CoordinateType c(0);
     for (size_t i=0; i<coefficients_.size(); i++)
       c.axpy(w[i][0], coefficients_[i].globalCoordinates());
 
+    // Possibly project onto target space
     if constexpr (conforming)
       return TargetSpace::projectOnto(c);
     else
@@ -148,9 +157,9 @@ namespace Dune::GFE
 
   }
 
-  template <int dim, class ctype, class LocalFiniteElement, class TargetSpace,bool conforming>
-  typename LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::DerivativeType
-  LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
+  template <class Basis, class TargetSpace,bool conforming>
+  typename LocalProjectedFEFunction<Basis,TargetSpace,conforming>::DerivativeType
+  LocalProjectedFEFunction<Basis,TargetSpace,conforming>::
   evaluateDerivative(const LocalCoordinate& local) const
   {
     if constexpr(conforming)
@@ -162,7 +171,7 @@ namespace Dune::GFE
       return evaluateDerivative(local, q);
     }
     else {
-      std::vector<Dune::FieldMatrix<ctype, 1, dim> > wDer;
+      std::vector<ScalarFEJacobian> wDer;
       localFiniteElement_.localBasis().evaluateJacobian(local, wDer);
 
       Dune::FieldMatrix<RT, embeddedDim, dim> derivative(0);
@@ -175,16 +184,16 @@ namespace Dune::GFE
     }
   }
 
-  template <int dim, class ctype, class LocalFiniteElement, class TargetSpace,bool conforming>
-  typename LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::DerivativeType
-  LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
+  template <class Basis, class TargetSpace,bool conforming>
+  typename LocalProjectedFEFunction<Basis,TargetSpace,conforming>::DerivativeType
+  LocalProjectedFEFunction<Basis,TargetSpace,conforming>::
   evaluateDerivative(const LocalCoordinate& local, const TargetSpace& q) const
   {
-    // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
-    std::vector<Dune::FieldVector<ctype,1> > w;
+    // Compute value and derivative of FE approximation
+    std::vector<ScalarFERange> w;
     localFiniteElement_.localBasis().evaluateFunction(local,w);
 
-    std::vector<Dune::FieldMatrix<ctype,1,dim> > wDer;
+    std::vector<ScalarFEJacobian> wDer;
     localFiniteElement_.localBasis().evaluateJacobian(local,wDer);
 
     typename TargetSpace::CoordinateType embeddedInterpolation(0);
@@ -202,9 +211,9 @@ namespace Dune::GFE
     return derivativeOfProjection*derivative;
   }
 
-  template <int dim, class ctype, class LocalFiniteElement, class TargetSpace,bool conforming>
+  template <class Basis, class TargetSpace,bool conforming>
   auto
-  LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
+  LocalProjectedFEFunction<Basis,TargetSpace,conforming>::
   evaluateValueAndDerivative(const LocalCoordinate& local) const
   {
     // Construct the type of the result -- it depends on whether the interpolation
@@ -217,7 +226,7 @@ namespace Dune::GFE
     //  Compute the value of the interpolation function
     ///////////////////////////////////////////////////////////
 
-    std::vector<Dune::FieldVector<ctype,1> > w;
+    std::vector<ScalarFERange> w;
     localFiniteElement_.localBasis().evaluateFunction(local,w);
 
     typename TargetSpace::CoordinateType embeddedInterpolation(0);
@@ -234,7 +243,7 @@ namespace Dune::GFE
     ///////////////////////////////////////////////////////////
 
     // Compute the interpolation in the surrounding space
-    std::vector<Dune::FieldMatrix<ctype,1,dim> > wDer;
+    std::vector<ScalarFEJacobian> wDer;
     localFiniteElement_.localBasis().evaluateJacobian(local,wDer);
 
     result.second = 0.0;
@@ -253,17 +262,27 @@ namespace Dune::GFE
 
   /** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold -- specialization for SO(3)
    *
-   * \tparam dim Dimension of the reference element
-   * \tparam ctype Type used for coordinates on the reference element
-   * \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
+   * \tparam Basis The basis used for the interpolation in the embedding space.
+   * \tparam field_type The number type used for function values
+   * \tparam conforming Geometrical conformity of the functions (omits projections when set to false)
+   *
+   * \note Even though the embedding space is 9-dimensional, we allow this basis
+   * to be scalar-valued.
    */
-  template <int dim, class ctype, class LocalFiniteElement, class field_type, bool conforming>
-  class LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3>,conforming>
+  template <class Basis, class field_type, bool conforming>
+  class LocalProjectedFEFunction<Basis,Rotation<field_type,3>,conforming>
   {
   public:
     typedef Rotation<field_type,3> TargetSpace;
   private:
-    using LocalCoordinate = FieldVector<ctype,dim>;
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto dim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using ScalarFERange = typename LocalBasis::Traits::RangeType;
+    using ScalarFEJacobian = typename LocalBasis::Traits::JacobianType;
 
     typedef typename TargetSpace::ctype RT;
 
@@ -355,7 +374,7 @@ namespace Dune::GFE
     template<class U>
     struct rebind
     {
-      using other = LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,U>;
+      using other = LocalProjectedFEFunction<Basis,U>;
     };
 
     /** \brief The number of Lagrange points */
@@ -391,7 +410,7 @@ namespace Dune::GFE
       Rotation<field_type,3> result;
 
       // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
-      std::vector<Dune::FieldVector<ctype,1> > w;
+      std::vector<ScalarFERange> w;
       localFiniteElement_.localBasis().evaluateFunction(local,w);
 
       // Interpolate in R^{3x3}
@@ -427,10 +446,10 @@ namespace Dune::GFE
                                       const TargetSpace& q) const
     {
       // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
-      std::vector<Dune::FieldVector<ctype,1> > w;
+      std::vector<ScalarFERange> w;
       localFiniteElement_.localBasis().evaluateFunction(local,w);
 
-      std::vector<Dune::FieldMatrix<ctype,1,dim> > wDer;
+      std::vector<ScalarFEJacobian> wDer;
       localFiniteElement_.localBasis().evaluateJacobian(local,wDer);
 
       // Compute matrix representations for all coefficients (we only have them in quaternion representation)
@@ -492,7 +511,7 @@ namespace Dune::GFE
       std::pair<Value,DerivativeType> result;
 
       // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
-      std::vector<Dune::FieldVector<ctype,1> > w;
+      std::vector<ScalarFERange> w;
       localFiniteElement_.localBasis().evaluateFunction(local,w);
 
       // Interpolate in R^{3x3}
@@ -533,17 +552,26 @@ namespace Dune::GFE
 
   /** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold -- specialization for R^3 x SO(3)
    *
-   * \tparam dim Dimension of the reference element
-   * \tparam ctype Type used for coordinates on the reference element
-   * \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
+   * \tparam Basis The basis used for the interpolation in the embedding space.
+   * \tparam field_type The number type used for function values
+   *
+   * \note Even though the embedding space is 12-dimensional, we allow this basis
+   * to be scalar-valued.
    */
-  template <int dim, class ctype, class LocalFiniteElement, class field_type>
-  class LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
+  template <class Basis, class field_type>
+  class LocalProjectedFEFunction<Basis,ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
   {
   public:
     using TargetSpace = ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> >;
   private:
-    using LocalCoordinate = FieldVector<ctype,dim>;
+    using Element = typename Basis::GridView::template Codim<0>::Entity;
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+    static constexpr auto dim = LocalCoordinate::size();
+
+    using LocalFiniteElement = typename Basis::LocalView::Tree::FiniteElement;
+    using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+    using ScalarFERange = typename LocalBasis::Traits::RangeType;
+    using ScalarFEJacobian = typename LocalBasis::Traits::JacobianType;
 
     typedef typename TargetSpace::ctype RT;
 
@@ -575,14 +603,14 @@ namespace Dune::GFE
       for (size_t i=0; i<coefficients.size(); i++)
         orientationCoefficients[i] = coefficients[i][_1];
 
-      orientationFunction_ = std::make_unique<LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3> > > (localFiniteElement,orientationCoefficients);
+      orientationFunction_ = std::make_unique<LocalProjectedFEFunction<Basis,Rotation<field_type,3> > > (localFiniteElement,orientationCoefficients);
     }
 
     /** \brief Rebind the FEFunction to another TargetSpace */
     template<class U>
     struct rebind
     {
-      using other = LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,U>;
+      using other = LocalProjectedFEFunction<Basis,U>;
     };
 
     /** \brief The number of Lagrange points */
@@ -605,7 +633,7 @@ namespace Dune::GFE
       TargetSpace result;
 
       // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
-      std::vector<Dune::FieldVector<ctype,1> > w;
+      std::vector<ScalarFERange> w;
       localFiniteElement_.localBasis().evaluateFunction(local,w);
 
       result[_0] = FieldVector<field_type,3>(0.0);
@@ -639,7 +667,7 @@ namespace Dune::GFE
       DerivativeType result(0);
 
       // get translation part
-      std::vector<Dune::FieldMatrix<ctype,1,dim> > sfDer(translationCoefficients_.size());
+      std::vector<ScalarFEJacobian> sfDer(translationCoefficients_.size());
       localFiniteElement_.localBasis().evaluateJacobian(local, sfDer);
 
       for (size_t i=0; i<translationCoefficients_.size(); i++)
@@ -650,7 +678,7 @@ namespace Dune::GFE
       Dune::FieldMatrix<field_type,4,dim> qResult = orientationFunction_->evaluateDerivative(local,q[_1]);
 
       for (int i=0; i<4; i++)
-        for (int j=0; j<dim; j++)
+        for (std::size_t j=0; j<dim; j++)
           result[3+i][j] = qResult[i][j];
 
       return result;
@@ -675,7 +703,7 @@ namespace Dune::GFE
     // we need access to the coefficients for the various factor spaces separately.
     std::vector<Dune::FieldVector<field_type,3> > translationCoefficients_;
 
-    std::unique_ptr<LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type, 3> > > orientationFunction_;
+    std::unique_ptr<LocalProjectedFEFunction<Basis,Rotation<field_type, 3> > > orientationFunction_;
 
 
   };
diff --git a/src/compute-disc-error.cc b/src/compute-disc-error.cc
index f740428c49d17b4745bd92c09d28de892f6c1db4..839aa2fa50c4c136faa68da9ce4e2dbb399518b7 100644
--- a/src/compute-disc-error.cc
+++ b/src/compute-disc-error.cc
@@ -165,7 +165,7 @@ void measureDiscreteEOC(const GridView gridView,
   //typedef LocalGeodesicFEFunction<FEBasis, TargetSpace> LocalInterpolationRule;
   //if (parameterSet["interpolationMethod"] != "geodesic")
   //  DUNE_THROW(Exception, "Inconsistent choice of interpolation method");
-  typedef GFE::LocalProjectedFEFunction<GridView::dimension, double, typename FEBasis::LocalView::Tree::FiniteElement, TargetSpace> LocalInterpolationRule;
+  typedef GFE::LocalProjectedFEFunction<FEBasis, TargetSpace> LocalInterpolationRule;
   if (parameterSet["interpolationMethod"] != "projected")
     DUNE_THROW(Exception, "Inconsistent choice of interpolation method");
   std::cout << "Using local interpolation: " << className<LocalInterpolationRule>() << std::endl;
@@ -479,7 +479,7 @@ void measureAnalyticalEOC(const GridView gridView,
 
   // ONly used if parameterSet["interpolationMethod"] == "projected"
   auto numericalSolutionProjected = GFE::EmbeddedGlobalGFEFunction<FEBasis,
-      GFE::LocalProjectedFEFunction<dim, double, typename FEBasis::LocalView::Tree::FiniteElement, TargetSpace>,
+      GFE::LocalProjectedFEFunction<FEBasis, TargetSpace>,
       TargetSpace> (feBasis, x);
   auto localNumericalSolutionProjected = localFunction(numericalSolutionProjected);
 #else
@@ -498,7 +498,7 @@ void measureAnalyticalEOC(const GridView gridView,
 
   if (parameterSet["interpolationMethod"] == "projected")
     numericalSolution = std::make_unique<GFE::EmbeddedGlobalGFEFunction<FEBasis,
-        GFE::LocalProjectedFEFunction<dim, double, typename FEBasis::LocalView::Tree::FiniteElement, TargetSpace>,
+        GFE::LocalProjectedFEFunction<FEBasis, TargetSpace>,
         TargetSpace> > (feBasis, x);
 #endif
 
diff --git a/src/cosserat-rod.cc b/src/cosserat-rod.cc
index c4cd386d3cbdb8958b2f28eb701dc16f9695cd82..3df2a0ccc2c2bf1e9e539d2347fe983086dc56ff 100644
--- a/src/cosserat-rod.cc
+++ b/src/cosserat-rod.cc
@@ -221,7 +221,7 @@ int main (int argc, char *argv[]) try
 
   using ATargetSpace = TargetSpace::rebind<adouble>::other;
   using GeodesicInterpolationRule  = GFE::LocalGeodesicFEFunction<ScalarBasis, ATargetSpace>;
-  using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<1, double, ScalarBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+  using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<ScalarBasis, ATargetSpace>;
 
   // Assembler using ADOL-C
   std::shared_ptr<GFE::LocalEnergy<ScalarBasis,ATargetSpace> > localRodEnergy;
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 763e871b18fb3da4de5cff1b0f05b8be82512d0b..667b184b833d0e3d7a8748a67f9afe0abcdfd22b 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -289,7 +289,7 @@ int main (int argc, char *argv[])
 
   // Next: The local energy, i.e., the integral of the density over one element
   using GeodesicInterpolationRule  = GFE::LocalGeodesicFEFunction<FEBasis, ATargetSpace>;
-  using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+  using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<FEBasis, ATargetSpace>;
 
   std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy;
 
diff --git a/src/simofoxshell.cc b/src/simofoxshell.cc
index 3c47af5c798a6074b3bd55c02a9ba927be9a305e..f437dab615561da64944b1d6e55aa1db5af0cf64 100644
--- a/src/simofoxshell.cc
+++ b/src/simofoxshell.cc
@@ -42,9 +42,9 @@
 
 #include <dune/vtk/vtkreader.hh>
 
-template <int dim, class ctype, class LocalFiniteElement, class TS>
-using LocalFEFunction = Dune::GFE::LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TS>;
-//using LocalFEFunction = LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TS>;
+template <class Basis, class TS>
+using LocalFEFunction = Dune::GFE::LocalProjectedFEFunction<Basis,TS>;
+//using LocalFEFunction = LocalGeodesicFEFunction<Basis,TS>;
 
 // Order of the approximation space for the midsurface position
 const int midsurfaceOrder = 1;
diff --git a/test/cosseratrodtest.cc b/test/cosseratrodtest.cc
index 5b367c3644dc2f1044cd0fbaec0276d69323ab4a..51b42cd5ccc97273e438b8ae77c6fad7c6087cb5 100644
--- a/test/cosseratrodtest.cc
+++ b/test/cosseratrodtest.cc
@@ -133,7 +133,7 @@ int main (int argc, char *argv[]) try
 
   using ATargetSpace = TargetSpace::rebind<adouble>::other;
   using GeodesicInterpolationRule  = GFE::LocalGeodesicFEFunction<ScalarBasis, ATargetSpace>;
-  using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<1, double, ScalarBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+  using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<ScalarBasis, ATargetSpace>;
 
   // Assembler using ADOL-C
   std::shared_ptr<GFE::LocalEnergy<ScalarBasis,ATargetSpace> > localRodEnergy;
diff --git a/test/harmonicmaptest.cc b/test/harmonicmaptest.cc
index 71c157f85d784e104dc0de9472913f6f76a79fa0..10b8c715e734a138a8e3f377bce39dc6bdf46ab8 100644
--- a/test/harmonicmaptest.cc
+++ b/test/harmonicmaptest.cc
@@ -153,9 +153,9 @@ int main (int argc, char *argv[])
   using InterpolationRule = GFE::LocalGeodesicFEFunction<FEBasis, TargetSpace>;
 #else
 #if CONFORMING
-  using InterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace,true>;
+  using InterpolationRule = GFE::LocalProjectedFEFunction<FEBasis, TargetSpace,true>;
 #else
-  using InterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace,false>;
+  using InterpolationRule = GFE::LocalProjectedFEFunction<FEBasis, TargetSpace,false>;
 #endif
 #endif
 
diff --git a/test/interpolationderivativestest.cc b/test/interpolationderivativestest.cc
index 309da6264d292e266c711597cd94dac466e267bd..3860337a9116fdac0390663ea3183e3d419149f9 100644
--- a/test/interpolationderivativestest.cc
+++ b/test/interpolationderivativestest.cc
@@ -408,10 +408,7 @@ TestSuite checkDerivatives()
   using GeodesicInterpolationRule = GFE::LocalGeodesicFEFunction<decltype(scalarBasis),
       TargetSpace>;
 
-  using ProjectionBasedInterpolationRule = GFE::LocalProjectedFEFunction<domainDim,
-      typename Grid::ctype,
-      decltype(scalarBasis.localView().tree().finiteElement()),
-      TargetSpace>;
+  using ProjectionBasedInterpolationRule = GFE::LocalProjectedFEFunction<decltype(scalarBasis),TargetSpace>;
 
   // Select the one to test
   using LocalInterpolationRule = std::conditional_t<interpolationType==InterpolationType::Geodesic,
diff --git a/test/localintegralstiffnesstest.cc b/test/localintegralstiffnesstest.cc
index 0b3f8bd200f4169bcf3f3ecd2a3995eac326afde..f44ef299ff98f896ff67aceddd5def6b4f7e3c89 100644
--- a/test/localintegralstiffnesstest.cc
+++ b/test/localintegralstiffnesstest.cc
@@ -129,14 +129,10 @@ int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
     else
       std::cout << "Using nonconforming interpolation" << std::endl;
 
-    using LocalInterpolationRule = GFE::LocalProjectedFEFunction<dim,
-        typename GridView::ctype,
-        decltype(feBasis.localView().tree().finiteElement()),
+    using LocalInterpolationRule = GFE::LocalProjectedFEFunction<FEBasis,
         TargetSpace, interpolationType!=Nonconforming>;
 
-    using ALocalInterpolationRule = GFE::LocalProjectedFEFunction<dim,
-        typename GridView::ctype,
-        decltype(feBasis.localView().tree().finiteElement()),
+    using ALocalInterpolationRule = GFE::LocalProjectedFEFunction<FEBasis,
         ATargetSpace, interpolationType!=Nonconforming>;
 
     // Assemble using the old assembler
@@ -341,13 +337,13 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
   {
     std::cout << "Using projection-based interpolation" << std::endl;
 
-    using LocalDeformationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), GFE::RealTuple<double,dim> >;
-    using LocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), GFE::Rotation<double,dim> >;
+    using LocalDeformationInterpolationRule = GFE::LocalProjectedFEFunction<decltype(deformationFEBasis), GFE::RealTuple<double,dim> >;
+    using LocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<decltype(orientationFEBasis), GFE::Rotation<double,dim> >;
 
     using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
 
-    using ALocalDeformationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), GFE::RealTuple<adouble,dim> >;
-    using ALocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), GFE::Rotation<adouble,dim> >;
+    using ALocalDeformationInterpolationRule = GFE::LocalProjectedFEFunction<decltype(deformationFEBasis), GFE::RealTuple<adouble,dim> >;
+    using ALocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<decltype(orientationFEBasis), GFE::Rotation<adouble,dim> >;
 
     using ALocalInterpolationRule = std::tuple<ALocalDeformationInterpolationRule,ALocalOrientationInterpolationRule>;
 
diff --git a/test/localprojectedfefunctiontest.cc b/test/localprojectedfefunctiontest.cc
index 9ad47b0b0ce1375ddf7bf5beaf58f3cb2d81f330..8b5473bf3bcc0e8727c4a9b1227b1f4bb39646ee 100644
--- a/test/localprojectedfefunctiontest.cc
+++ b/test/localprojectedfefunctiontest.cc
@@ -10,7 +10,10 @@
 #include <dune/geometry/type.hh>
 #include <dune/geometry/referenceelements.hh>
 
-#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
+#include <dune/grid/onedgrid.hh>
+#include <dune/grid/uggrid.hh>
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
 
 #include <dune/gfe/functions/localprojectedfefunction.hh>
 #include <dune/gfe/spaces/productmanifold.hh>
@@ -122,12 +125,21 @@ void testPermutationInvariance(const std::vector<TargetSpace>& corners)
   if (domainDim!=2)
     return;
 
-  LagrangeLocalFiniteElementCache<double,double,domainDim,1> feCache;
-  typedef typename LagrangeLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement;
+  // Make a test grid
+  using Grid = UGGrid<domainDim>;
+  GridFactory<Grid> gridFactory;
+
+  gridFactory.insertVertex({0.0,0.0});
+  gridFactory.insertVertex({1.0,0.0});
+  gridFactory.insertVertex({0.0,1.0});
 
-  GeometryType simplex = GeometryTypes::simplex(domainDim);
+  gridFactory.insertElement(GeometryTypes::simplex(domainDim), {0,1,2});
 
-  //
+  const auto grid = gridFactory.createGrid();
+  const auto gridView = grid->leafGridView();
+  using GridView = decltype(gridView);
+
+  // Create the test configurations
   std::vector<TargetSpace> cornersRotated1(domainDim+1);
   std::vector<TargetSpace> cornersRotated2(domainDim+1);
 
@@ -135,15 +147,22 @@ void testPermutationInvariance(const std::vector<TargetSpace>& corners)
   cornersRotated1[1] = cornersRotated2[0] = corners[2];
   cornersRotated1[2] = cornersRotated2[1] = corners[0];
 
-  GFE::LocalProjectedFEFunction<2,double,LocalFiniteElement,TargetSpace> f0(feCache.get(simplex), corners);
-  GFE::LocalProjectedFEFunction<2,double,LocalFiniteElement,TargetSpace> f1(feCache.get(simplex), cornersRotated1);
-  GFE::LocalProjectedFEFunction<2,double,LocalFiniteElement,TargetSpace> f2(feCache.get(simplex), cornersRotated2);
+  using InterpolationBasis = Functions::LagrangeBasis<GridView,1>;
+  InterpolationBasis interpolationBasis(gridView);
+
+  auto localBasisView = interpolationBasis.localView();
+  localBasisView.bind(*gridView.template begin<0>());
+  const auto& localFiniteElement = localBasisView.tree().finiteElement();
+
+  GFE::LocalProjectedFEFunction<InterpolationBasis,TargetSpace> f0(localFiniteElement, corners);
+  GFE::LocalProjectedFEFunction<InterpolationBasis,TargetSpace> f1(localFiniteElement, cornersRotated1);
+  GFE::LocalProjectedFEFunction<InterpolationBasis,TargetSpace> f2(localFiniteElement, cornersRotated2);
 
   // A quadrature rule as a set of test points
   int quadOrder = 3;
 
-  const Dune::QuadratureRule<double, domainDim>& quad
-    = Dune::QuadratureRules<double, domainDim>::rule(simplex, quadOrder);
+  const auto& quad
+    = QuadratureRules<double, domainDim>::rule(GeometryTypes::simplex(domainDim), quadOrder);
 
   for (size_t pt=0; pt<quad.size(); pt++) {
 
@@ -171,9 +190,10 @@ void testPermutationInvariance(const std::vector<TargetSpace>& corners)
 
 }
 
-template <int domainDim, class TargetSpace, bool conforming=true>
-void testDerivative(const GFE::LocalProjectedFEFunction<domainDim,double,typename LagrangeLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType, TargetSpace, conforming>& f)
+template <class Basis, class TargetSpace, bool conforming=true>
+void testDerivative(const GFE::LocalProjectedFEFunction<Basis, TargetSpace, conforming>& f)
 {
+  static constexpr auto domainDim = Basis::GridView::template Codim<0>::Entity::mydimension;
   static const int embeddedDim = TargetSpace::EmbeddedTangentVector::dimension;
 
   // A quadrature rule as a set of test points
@@ -213,6 +233,26 @@ void test(const GeometryType& element)
 {
   std::cout << " --- Testing " << className<TargetSpace>() << ", domain dimension: " << element.dim() << " ---" << std::endl;
 
+  // Make a test grid
+  using Grid = std::conditional_t<domainDim==1,OneDGrid,UGGrid<domainDim> >;
+  GridFactory<Grid> gridFactory;
+
+  const auto referenceElement = Dune::referenceElement<double,domainDim>(element);
+
+  std::vector<unsigned int> refElementCorners(referenceElement.size(domainDim));
+  for (int i=0; i<referenceElement.size(domainDim); i++)
+  {
+    gridFactory.insertVertex(referenceElement.position(i,domainDim));
+    refElementCorners[i] = i;
+  }
+
+  gridFactory.insertElement(element, refElementCorners);
+
+  const auto grid = gridFactory.createGrid();
+  const auto gridView = grid->leafGridView();
+  using GridView = decltype(gridView);
+
+  // Make test point sets
   std::vector<TargetSpace> testPoints;
   GFE::ValueFactory<TargetSpace>::get(testPoints);
 
@@ -234,22 +274,28 @@ void test(const GeometryType& element)
       continue;
 
     // Make local gfe function to be tested
-    LagrangeLocalFiniteElementCache<double,double,domainDim,1> feCache;
-    typedef typename LagrangeLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement;
+    using InterpolationBasis = Functions::LagrangeBasis<GridView,1>;
+    InterpolationBasis interpolationBasis(gridView);
+
+    auto localBasisView = interpolationBasis.localView();
+    localBasisView.bind(*gridView.template begin<0>());
+    const auto& localFiniteElement = localBasisView.tree().finiteElement();
 
-    GFE::LocalProjectedFEFunction<domainDim,double,LocalFiniteElement,TargetSpace> f(feCache.get(element),corners);
-    GFE::LocalProjectedFEFunction<domainDim, double, LocalFiniteElement, TargetSpace,false> f_nonconforming(feCache.get(element), corners);
+    GFE::LocalProjectedFEFunction<InterpolationBasis,TargetSpace> f(localFiniteElement,corners);
+    GFE::LocalProjectedFEFunction<InterpolationBasis, TargetSpace,false> f_nonconforming(localFiniteElement,corners);
 
     //testPermutationInvariance(corners);
-    testDerivative<domainDim>(f);
-    testDerivative<domainDim>(f_nonconforming);
+    testDerivative(f);
+    testDerivative(f_nonconforming);
   }
 
 }
 
 
-int main()
+int main(int argc, char *argv[])
 {
+  MPIHelper::instance(argc, argv);
+
   // choke on NaN -- don't enable this by default, as there are
   // a few harmless NaN in the loopsolver
   //feenableexcept(FE_INVALID);