diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 7e2ae2cb9dfd3b81db3c17abce0542957927984a..212b40b1022b740ad2aa9b13511a0f948e6f855f 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -5,7 +5,6 @@
 #include <dune/common/parametertree.hh>
 #include <dune/geometry/quadraturerules.hh>
 
-#include <dune/fufem/functions/virtualgridfunction.hh>
 #include <dune/fufem/boundarypatch.hh>
 
 #include <dune/gfe/localenergy.hh>
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index 48506937e535f63b93066fa090d40d7d444c42b2..09e4c0c22d009bccfe56bce3c7bccacdaf727aa0 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -7,7 +7,6 @@
 
 #include <dune/matrix-vector/crossproduct.hh>
 
-#include <dune/fufem/functions/virtualgridfunction.hh>
 #include <dune/fufem/boundarypatch.hh>
 
 #include <dune/gfe/localenergy.hh>
diff --git a/src/compute-disc-error.cc b/src/compute-disc-error.cc
index 6fa0ca0059ced22006ffd29436fdead6fab2a93f..b320954964b5da277f85105775714f4908743643 100644
--- a/src/compute-disc-error.cc
+++ b/src/compute-disc-error.cc
@@ -443,7 +443,7 @@ void measureAnalyticalEOC(const GridView gridView,
   //   Measure the discretization error
   /////////////////////////////////////////////////////////////////
 
-  // Read reference solution and its derivative into a PythonFunction
+  // Read reference solution and its derivative into a Python function
   typedef VirtualDifferentiableFunction<FieldVector<double, dimworld>, typename TargetSpace::CoordinateType> FBase;
 
   Python::Module module = Python::import(parameterSet.get<std::string>("referenceSolution"));
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 294d78a0beee1011538b0ef7109decd7cd7f1523..3ed30be55becc0c59b444c9ce426abdba0abe555 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -32,7 +32,6 @@
 
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/functiontools/boundarydofs.hh>
-#include <dune/fufem/functiontools/basisinterpolator.hh>
 #include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
 #include <dune/fufem/dunepython.hh>
 
@@ -230,11 +229,11 @@ int main (int argc, char *argv[]) try
     // Make Python function that computes which vertices are on the Dirichlet boundary,
     // based on the vertex positions.
     std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
-    PythonFunction<FieldVector<double,dimworld>, bool> pythonDirichletVertices(Python::evaluate(lambda));
+    auto pythonDirichletVertices = Python::make_function<bool>(Python::evaluate(lambda));
 
     // Same for the Neumann boundary
     lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
-    PythonFunction<FieldVector<double,dimworld>, bool> pythonNeumannVertices(Python::evaluate(lambda));
+    auto pythonNeumannVertices = Python::make_function<bool>(Python::evaluate(lambda));
 
     for (auto&& vertex : vertices(gridView))
     {
@@ -304,25 +303,15 @@ int main (int argc, char *argv[]) try
     x[_0].resize(deformationFEBasis.size());
 
     lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
-    PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
+    auto pythonInitialDeformation = Python::make_function<FieldVector<double,3> >(Python::evaluate(lambda));
 
     std::vector<FieldVector<double,3> > v;
-    ::Functions::interpolate(fufemDeformationFEBasis, v, pythonInitialDeformation);
+    Functions::interpolate(deformationFEBasis, v, pythonInitialDeformation);
 
     for (size_t i=0; i<x[_0].size(); i++)
       x[_0][i] = v[i];
 
     x[_1].resize(orientationFEBasis.size());
-#if 0
-    lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
-    PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
-
-    std::vector<FieldVector<double,3> > v;
-    Functions::interpolate(feBasis, v, pythonInitialDeformation);
-
-    for (size_t i=0; i<x.size(); i++)
-      xDisp[i] = v[i];
-#endif
 #else
     SolutionType x(feBasis.size());
 
@@ -364,10 +353,10 @@ int main (int argc, char *argv[]) try
 
     } else {
     lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
-    PythonFunction<FieldVector<double,dimworld>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
+      auto pythonInitialDeformation = Python::make_function<FieldVector<double,3> >(Python::evaluate(lambda));
 
     std::vector<FieldVector<double,3> > v;
-      ::Functions::interpolate(fufemFeBasis, v, pythonInitialDeformation);
+      Functions::interpolate(feBasis, v, pythonInitialDeformation);
 
     for (size_t i=0; i<x.size(); i++)
       x[i].r = v[i];
@@ -514,15 +503,15 @@ int main (int argc, char *argv[]) try
         Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
 
         // Extract object member functions as Dune functions
-        PythonFunction<FieldVector<double,dimworld>, FieldVector<double,3> >   deformationDirichletValues(dirichletValuesPythonObject.get("deformation"));
-        PythonFunction<FieldVector<double,dimworld>, FieldMatrix<double,3,3> > orientationDirichletValues(dirichletValuesPythonObject.get("orientation"));
+        auto deformationDirichletValues = Python::make_function<FieldVector<double,3> >(dirichletValuesPythonObject.get("deformation"));
+        auto orientationDirichletValues = Python::make_function<FieldMatrix<double,3,3> >(dirichletValuesPythonObject.get("orientation"));
 
         std::vector<FieldVector<double,3> > ddV;
         std::vector<FieldMatrix<double,3,3> > dOV;
 
 #ifdef MIXED_SPACE
-        ::Functions::interpolate(fufemDeformationFEBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
-        ::Functions::interpolate(fufemOrientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
+        Functions::interpolate(deformationFEBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
+        Functions::interpolate(orientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
 
         for (size_t j=0; j<x[_0].size(); j++)
           if (deformationDirichletNodes[j][0])
@@ -532,8 +521,8 @@ int main (int argc, char *argv[]) try
           if (orientationDirichletNodes[j][0])
             x[_1][j].set(dOV[j]);
 #else
-        ::Functions::interpolate(fufemFeBasis, ddV, deformationDirichletValues, dirichletDofs);
-        ::Functions::interpolate(fufemFeBasis, dOV, orientationDirichletValues, dirichletDofs);
+        Functions::interpolate(feBasis, ddV, deformationDirichletValues, dirichletDofs);
+        Functions::interpolate(feBasis, dOV, orientationDirichletValues, dirichletDofs);
 
         for (size_t j=0; j<x.size(); j++)
           if (dirichletNodes[j][0])
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 235be15b26ad2335940e2d56fc2dac0f875092f9..3921c142c19f218566896afbbd85f3f82f206df6 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -177,15 +177,15 @@ int main (int argc, char *argv[]) try
   // Make Python function that computes which vertices are on the Dirichlet boundary,
   // based on the vertex positions.
   std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, FieldVector<bool,dim>> pythonDirichletVertices(Python::evaluate(lambda));
+  auto pythonDirichletVertices = Python::make_function<FieldVector<bool,dim> >(Python::evaluate(lambda));
 
   // Same for the Neumann boundary
   lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
+  auto pythonNeumannVertices = Python::make_function<bool>(Python::evaluate(lambda));
 
   // Same for the Surface Shell Boundary
   lambda = std::string("lambda x: (") + parameterSet.get<std::string>("surfaceShellVerticesPredicate", "0") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, bool> pythonSurfaceShellVertices(Python::evaluate(lambda));
+  auto pythonSurfaceShellVertices = Python::make_function<bool>(Python::evaluate(lambda));
 
   while (numLevels > 0) {
     for (auto&& e : elements(grid->leafGridView())){
@@ -326,7 +326,7 @@ int main (int argc, char *argv[]) try
   //Initial deformation of the underlying substrate
   BlockVector<FieldVector<double,dim> > displacement(compositeBasis.size({0}));
   lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > pythonInitialDeformation(Python::evaluate(lambda));
+  auto pythonInitialDeformation = Python::make_function<FieldVector<double,dim> >(Python::evaluate(lambda));
   Dune::Functions::interpolate(deformationPowerBasis, displacement, pythonInitialDeformation);
 
   BlockVector<FieldVector<double,dim> > identity(compositeBasis.size({0}));
@@ -426,7 +426,7 @@ int main (int argc, char *argv[]) try
   } else {
     // Read grid deformation from deformation function
     auto gridDeformationLambda = std::string("lambda x: (") + parameterSet.get<std::string>("gridDeformation") + std::string(")");
-    PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > gridDeformation(Python::evaluate(gridDeformationLambda));
+    auto gridDeformation = Python::make_function<FieldVector<double,dim> >(Python::evaluate(gridDeformationLambda));
 
     //Iterate over boundary, each facet on the boundary has an element (boundaryElement.inside()) with a unique global id (idSet.subId);
     //we store the new geometry in the map with this id as reference
@@ -462,8 +462,8 @@ int main (int argc, char *argv[]) try
   Python::Reference surfaceShellClass = Python::import(materialParameters.get<std::string>("surfaceShellParameters"));
   Python::Callable surfaceShellCallable = surfaceShellClass.get("SurfaceShellParameters");
   Python::Reference pythonObject = surfaceShellCallable();
-  PythonFunction<Dune::FieldVector<double, dim>, double> fThickness(pythonObject.get("thickness"));
-  PythonFunction<Dune::FieldVector<double, dim>, Dune::FieldVector<double, 2>> fLame(pythonObject.get("lame"));
+  auto fThickness = Python::make_function<double>(pythonObject.get("thickness"));
+  auto fLame = Python::make_function<FieldVector<double, 2> >(pythonObject.get("lame"));
 
   for (int i = 0; i < numHomotopySteps; i++)
   {
@@ -546,8 +546,8 @@ int main (int argc, char *argv[]) try
     Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
 
     // Extract object member functions as Dune functions
-    PythonFunction<FieldVector<double,dim>, FieldVector<double,targetDim> >   deformationDirichletValues(dirichletValuesPythonObject.get("deformation"));
-    PythonFunction<FieldVector<double,dim>, FieldMatrix<double,targetDim,targetDim> > rotationalDirichletValues(dirichletValuesPythonObject.get("orientation"));
+    auto deformationDirichletValues = Python::make_function<FieldVector<double,targetDim> >   (dirichletValuesPythonObject.get("deformation"));
+    auto rotationalDirichletValues = Python::make_function<FieldMatrix<double,targetDim,targetDim> > (dirichletValuesPythonObject.get("orientation"));
 
     BlockVector<FieldVector<double,targetDim> > ddV;
     Dune::Functions::interpolate(deformationPowerBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
diff --git a/src/gradient-flow.cc b/src/gradient-flow.cc
index ac87cc85783c30796c7d2cceb34b054d895f5ed5..ea1be930de232f2fb8b5f7dd6ccaa777f195886f 100644
--- a/src/gradient-flow.cc
+++ b/src/gradient-flow.cc
@@ -23,10 +23,10 @@
 #include <dune/grid/io/file/vtk.hh>
 
 #include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 
 #include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/basisinterpolator.hh>
 #include <dune/fufem/functiontools/boundarydofs.hh>
 #include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
 #include <dune/fufem/discretizationerror.hh>
@@ -160,7 +160,7 @@ int main (int argc, char *argv[]) try
   // Make Python function that computes which vertices are on the Dirichlet boundary,
   // based on the vertex positions.
   std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
-  PythonFunction<FieldVector<double,dimworld>, bool> pythonDirichletVertices(Python::evaluate(lambda));
+  auto pythonDirichletVertices = Python::make_function<bool>(Python::evaluate(lambda));
 
   for (auto&& vertex : vertices(grid->leafGridView()))
   {
@@ -177,14 +177,12 @@ int main (int argc, char *argv[]) try
   //   Initial iterate
   ////////////////////////////
 
-  // Read initial iterate into a PythonFunction
-  typedef PythonFunction<FieldVector<double, dimworld>, TargetSpace::CoordinateType> FBase;
-
+  // Read initial iterate into a Python function
   Python::Module module = Python::import(parameterSet.get<std::string>("initialIterate"));
-  auto pythonInitialIterate = module.get("f").toC<std::shared_ptr<FBase>>();
+  auto pythonInitialIterate = Python::make_function<TargetSpace::CoordinateType>(module.get("f"));
 
   std::vector<TargetSpace::CoordinateType> v;
-  ::Functions::interpolate(fufemFeBasis, v, *pythonInitialIterate);
+  Functions::interpolate(feBasis, v, pythonInitialIterate);
 
   SolutionType x(feBasis.size());
 
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 68c123ebf6a6a324e6b0e665f84e5d420cd31c38..ee6abd62c7d61405222a42ba759918dda43b941f 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -220,7 +220,7 @@ int main (int argc, char *argv[])
     // Make Python function that computes which vertices are on the Dirichlet boundary,
     // based on the vertex positions.
     std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
-    PythonFunction<FieldVector<double,dimworld>, bool> pythonDirichletVertices(Python::evaluate(lambda));
+    auto pythonDirichletVertices = Python::make_function<bool>(Python::evaluate(lambda));
 
     for (auto&& vertex : vertices(gridView))
     {
@@ -243,7 +243,7 @@ int main (int argc, char *argv[])
     //   Initial iterate
     // //////////////////////////
 
-    // Read initial iterate into a PythonFunction
+    // Read initial iterate into a Python function
     Python::Module module = Python::import(parameterSet.get<std::string>("initialIterate"));
     auto pythonInitialIterate = Python::makeFunction<TargetSpace::CoordinateType(const FieldVector<double,dimworld>&)>(module.get("f"));