diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 8f7c34528b63bd2aa66b62c272ec8c5dd58acc24..7e2ae2cb9dfd3b81db3c17abce0542957927984a 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -97,8 +97,8 @@ public:
      */
     CosseratEnergyLocalStiffness(const Dune::ParameterTree& parameters,
                                  const BoundaryPatch<GridView>* neumannBoundary,
-                                 const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > neumannFunction,
-                                 const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad)
+                                 const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction,
+                                 const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad)
     : neumannBoundary_(neumannBoundary),
       neumannFunction_(neumannFunction)
     {
@@ -271,10 +271,10 @@ public:
     const BoundaryPatch<GridView>* neumannBoundary_;
 
     /** \brief The function implementing the Neumann data */
-    const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > neumannFunction_;
+    const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction_;
 
     /** \brief The function implementing a volume load */
-    const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad_;
+    const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad_;
 };
 
 template <class Basis, int dim, class field_type>
@@ -367,12 +367,7 @@ energy(const typename Basis::LocalView& localView,
             continue;
 
         // Value of the volume load density at the current position
-        Dune::FieldVector<double,3> volumeLoadDensity;
-
-        if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_))
-            std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_)->evaluateLocal(element, quadPos, volumeLoadDensity);
-        else
-            volumeLoad_->evaluate(element.geometry().global(quad[pt].position()), volumeLoadDensity);
+        auto volumeLoadDensity = volumeLoad_(element.geometry().global(quad[pt].position()));
 
         // Only translational dofs are affected by the volume load
         for (size_t i=0; i<volumeLoadDensity.size(); i++)
@@ -406,12 +401,7 @@ energy(const typename Basis::LocalView& localView,
             RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
 
             // Value of the Neumann data at the current position
-            Dune::FieldVector<double,3> neumannValue;
-
-            if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_))
-                std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
-            else
-                neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
+            auto neumannValue = neumannFunction_(it.geometry().global(quad[pt].position()));
 
             // Only translational dofs are affected by the Neumann force
             for (size_t i=0; i<neumannValue.size(); i++)
@@ -550,12 +540,7 @@ energy(const typename Basis::LocalView& localView,
             RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
 
             // Value of the Neumann data at the current position
-            Dune::FieldVector<double,3> neumannValue;
-
-            if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_))
-                std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
-            else
-                neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
+            auto neumannValue = neumannFunction_(it.geometry().global(quad[pt].position()));
 
             // Only translational dofs are affected by the Neumann force
             for (size_t i=0; i<neumannValue.size(); i++)
diff --git a/dune/gfe/localenergy.hh b/dune/gfe/localenergy.hh
index ba6f8da46fcd2db1323418662ecf1b74dc2ff867..d80c0f0f5a01ea8c9645cac690a03ae3c0c9c4b9 100644
--- a/dune/gfe/localenergy.hh
+++ b/dune/gfe/localenergy.hh
@@ -8,19 +8,20 @@ namespace Dune {
 namespace GFE {
 
 /** \brief Base class for energies defined by integrating over one grid element */
-template<class Basis, class TargetSpace>
+template<class Basis, class... TargetSpaces>
 class LocalEnergy
 {
 public:
+  using RT = typename std::common_type<typename TargetSpaces::ctype...>::type;
 
   /** \brief Compute the energy
    *
    * \param localView Local view specifying the current element and the FE space there
    * \param coefficients The coefficients of a FE function on the current element
    */
-  virtual typename TargetSpace::ctype
+  virtual RT
   energy (const typename Basis::LocalView& localView,
-          const std::vector<TargetSpace>& localSolution) const = 0;
+          const std::vector<TargetSpaces>&... localSolution) const = 0;
 
   /** Empty virtual default destructor
    *
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index 102115dbc0a371ce6eb83a33e2e5dd739bd031d3..85e97321ea5719bec0602bcb7507fb9c45628922 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -42,8 +42,8 @@ public:
   NonplanarCosseratShellEnergy(const Dune::ParameterTree& parameters,
                                const std::vector<UnitVector<double,3> >& vertexNormals,
                                const BoundaryPatch<GridView>* neumannBoundary,
-                               const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > neumannFunction,
-                               const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad)
+                               const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction,
+                               const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad)
   : vertexNormals_(vertexNormals),
     neumannBoundary_(neumannBoundary),
     neumannFunction_(neumannFunction),
@@ -117,10 +117,10 @@ public:
   const BoundaryPatch<GridView>* neumannBoundary_;
 
   /** \brief The function implementing the Neumann data */
-  const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > neumannFunction_;
+  const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction_;
 
   /** \brief The function implementing a volume load */
-  const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad_;
+  const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad_;
 };
 
 template <class Basis, int dim, class field_type>
@@ -306,12 +306,7 @@ energy(const typename Basis::LocalView& localView,
       continue;
 
     // Value of the volume load density at the current position
-    Dune::FieldVector<double,3> volumeLoadDensity;
-
-    if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_))
-      std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_)->evaluateLocal(element, quadPos, volumeLoadDensity);
-    else
-      volumeLoad_->evaluate(geometry.global(quad[pt].position()), volumeLoadDensity);
+    Dune::FieldVector<double,3> volumeLoadDensity = volumeLoad_(geometry.global(quad[pt].position()));
 
     // Only translational dofs are affected by the volume load
     for (size_t i=0; i<volumeLoadDensity.size(); i++)
@@ -344,12 +339,7 @@ energy(const typename Basis::LocalView& localView,
       RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
 
       // Value of the Neumann data at the current position
-      Dune::FieldVector<double,3> neumannValue;
-
-      if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_))
-        std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
-      else
-        neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
+      Dune::FieldVector<double,3> neumannValue = neumannFunction_(it.geometry().global(quad[pt].position()));
 
       // Only translational dofs are affected by the Neumann force
       for (size_t i=0; i<neumannValue.size(); i++)
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 122c7f44021af786cdc273a8a1d5bc1cb690f6c2..a64eefe9aac68578bd3b64fb6de8474c524b0706 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -81,45 +81,6 @@ const int blocksize = TargetSpace::TangentVector::dimension;
 
 using namespace Dune;
 
-/** \brief A constant vector-valued function, for simple Neumann boundary values */
-struct NeumannFunction
-    : public Dune::VirtualFunction<FieldVector<double,dimworld>, FieldVector<double,3> >
-{
-    NeumannFunction(const FieldVector<double,3> values,
-                    double homotopyParameter)
-    : values_(values),
-      homotopyParameter_(homotopyParameter)
-    {}
-
-    void evaluate(const FieldVector<double, dimworld>& x, FieldVector<double,3>& out) const {
-        out = 0;
-        out.axpy(homotopyParameter_, values_);
-    }
-
-    FieldVector<double,3> values_;
-    double homotopyParameter_;
-};
-
-/** \brief A constant vector-valued function, for simple volume loads */
-struct VolumeLoad
-    : public Dune::VirtualFunction<FieldVector<double,dimworld>, FieldVector<double,3> >
-{
-    VolumeLoad(const FieldVector<double,3> values,
-               double homotopyParameter)
-    : values_(values),
-      homotopyParameter_(homotopyParameter)
-    {}
-
-    void evaluate(const FieldVector<double, dimworld>& x, FieldVector<double,3>& out) const {
-        out = 0;
-        out.axpy(homotopyParameter_, values_);
-    }
-
-    FieldVector<double,3> values_;
-    double homotopyParameter_;
-};
-
-
 int main (int argc, char *argv[]) try
 {
     // initialize MPI, finalize is done automatically on exit
@@ -438,15 +399,25 @@ int main (int argc, char *argv[]) try
     // ////////////////////////////////////////////////////////////
 
     const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
-    std::shared_ptr<NeumannFunction> neumannFunction;
+    FieldVector<double,3> neumannValues {0,0,0};
     if (parameterSet.hasKey("neumannValues"))
-        neumannFunction = std::make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
-                                                       homotopyParameter);
+        neumannValues = parameterSet.get<FieldVector<double,3> >("neumannValues");
 
-    std::shared_ptr<VolumeLoad> volumeLoad;
+    auto neumannFunction = [&]( FieldVector<double,dim> ) {
+      auto nV = neumannValues;
+      nV *= homotopyParameter;
+      return nV;
+    };
+
+    FieldVector<double,3> volumeLoadValues {0,0,0};
     if (parameterSet.hasKey("volumeLoad"))
-        volumeLoad = std::make_shared<VolumeLoad>(parameterSet.get<FieldVector<double,3> >("volumeLoad"),
-                                                                                          homotopyParameter);
+        volumeLoadValues = parameterSet.get<FieldVector<double,3> >("volumeLoad");
+
+    auto volumeLoad = [&]( FieldVector<double,dim>) {
+      auto vL = volumeLoadValues;
+      vL *= homotopyParameter;
+      return vL;
+    };
 
     if (mpiHelper.rank() == 0) {
         std::cout << "Material parameters:" << std::endl;
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 0f263551e83035ae31499a861edae80575a059d7..f0843d81f4e46dca10975d882c9b15ff7e45390a 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -93,27 +93,6 @@ struct Dune::MathematicalConstants<adouble>
 typedef adouble ValueType;
 
 using namespace Dune;
-#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 8)
-/** \brief A constant vector-valued function, for simple Neumann boundary values */
-struct NeumannFunction
-    : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,dim> >
-{
-  NeumannFunction(const FieldVector<double,dim> values,
-                  double homotopyParameter)
-  : values_(values),
-    homotopyParameter_(homotopyParameter)
-  {}
-
-  void evaluate(const FieldVector<double, dim>& x, FieldVector<double,dim>& out) const
-  {
-    out = 0;
-    out.axpy(-homotopyParameter_, values_);
-  }
-
-  FieldVector<double,dim> values_;
-  double homotopyParameter_;
-};
-#endif
 
 int main (int argc, char *argv[]) try
 {
@@ -472,9 +451,7 @@ int main (int argc, char *argv[]) try
       // A constant vector-valued function, for simple Neumann boundary values
     std::shared_ptr<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<double,dim>)>> neumannFunctionPtr;
     neumannFunctionPtr = std::make_shared<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<double,dim>)>>([&](FieldVector<double,dim> ) {
-      auto nv = neumannValues;
-      nv *= (-homotopyParameter);
-      return nv;
+      return neumannValues * (-homotopyParameter);
     });
 #endif
 
diff --git a/test/geodesicfeassemblerwrappertest.cc b/test/geodesicfeassemblerwrappertest.cc
index 4a44e4892e30d65562d11940702528800d60e285..3e78813537d2fc681af3723c304491920ccfc1f0 100644
--- a/test/geodesicfeassemblerwrappertest.cc
+++ b/test/geodesicfeassemblerwrappertest.cc
@@ -65,18 +65,6 @@ using MatrixType = MultiTypeBlockMatrix<MatrixRow0,MatrixRow1>;
 //Types for the Non-mixed space
 using RBM = RigidBodyMotion<double, dim>;
 const static int blocksize = RBM::TangentVector::dimension;
-
-struct NeumannFunction
-    : public Dune::VirtualFunction<FieldVector<double,gridDim>, FieldVector<double,dim> >
-{
-    NeumannFunction(){}
-
-    void evaluate(const FieldVector<double, gridDim>& x, FieldVector<double,dim>& out) const {
-        out = 0;
-        out.axpy(1.0, values_);
-    }
-    FieldVector<double,3> values_ = {3e4,2e4,1e4};
-};  
 using CorrectionTypeWrapped = BlockVector<FieldVector<double, blocksize> >;
 using MatrixTypeWrapped = BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >;
 
@@ -143,7 +131,11 @@ int main (int argc, char *argv[])
   parameters["kappa"] = "1";
 
 
-  auto neumannFunction = std::make_shared<NeumannFunction>();
+  FieldVector<double,dim> values_ = {3e4,2e4,1e4};
+  auto neumannFunction = [&](FieldVector<double, gridDim>){
+    return values_;
+  };
+
   CosseratEnergyLocalStiffness<decltype(compositeBasis), dim,adouble> cosseratEnergy(parameters,
                                                                      &neumannBoundary,
                                                                      neumannFunction,