diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 989e5e50b7b0047dec566d3a7950df250c37fc7d..748852951f674febcc017a4ec230f9e6b9628c7b 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -189,7 +189,7 @@ public:
      */
     CosseratEnergyLocalStiffness(const Dune::ParameterTree& parameters,
                                  const BoundaryPatch<GridView>* neumannBoundary,
-                                 const 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> > > neumannFunction,
                                  const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad)
     : neumannBoundary_(neumannBoundary),
       neumannFunction_(neumannFunction)
@@ -363,7 +363,7 @@ public:
     const BoundaryPatch<GridView>* neumannBoundary_;
 
     /** \brief The function implementing the Neumann data */
-    const 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> > > neumannFunction_;
 
     /** \brief The function implementing a volume load */
     const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad_;
@@ -501,8 +501,8 @@ energy(const typename Basis::LocalView& localView,
             // Value of the Neumann data at the current position
             Dune::FieldVector<double,3> neumannValue;
 
-            if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_))
-                dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_)->evaluateLocal(element, quadPos, 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);
 
@@ -646,8 +646,8 @@ energy(const typename Basis::LocalView& localView,
             // Value of the Neumann data at the current position
             Dune::FieldVector<double,3> neumannValue;
 
-            if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_))
-                dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_)->evaluateLocal(element, quadPos, 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);
 
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 75eea0156923c331b665229540474f24317d4647..c086218ef05372efd42aaab6e000e74e55486714 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -293,7 +293,7 @@ int main (int argc, char *argv[]) try
     CosseratEnergyLocalStiffness<FEBasis,
                                  3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
                                                                               &neumannBoundary,
-                                                                              neumannFunction.get(),
+                                                                              neumannFunction,
                                                                               volumeLoad);
     LocalGeodesicFEADOLCStiffness<FEBasis,
                                   TargetSpace> localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
diff --git a/src/mixed-cosserat-continuum.cc b/src/mixed-cosserat-continuum.cc
index 420f254f6c67313a3127a067b421f0a82f83a9e1..1204d07c0c655baadfeb95fe8def3b8444d8af3b 100644
--- a/src/mixed-cosserat-continuum.cc
+++ b/src/mixed-cosserat-continuum.cc
@@ -297,7 +297,8 @@ int main (int argc, char *argv[]) try
     CosseratEnergyLocalStiffness<decltype(compositeBasis),
                         3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
                                                                      &neumannBoundary,
-                                                                     neumannFunction.get());
+                                                                     neumannFunction,
+                                                                     nullptr);
 
     MixedLocalGFEADOLCStiffness<decltype(compositeBasis),
                                 RealTuple<double,3>,