diff --git a/dune/solvers/transferoperators/genericmultigridtransfer.hh b/dune/solvers/transferoperators/genericmultigridtransfer.hh
index 8fe413ca273f7364de06c8451c3bc892a7c39ee6..d9ab41955ff34506d89780585980f21c8a7f0e7a 100644
--- a/dune/solvers/transferoperators/genericmultigridtransfer.hh
+++ b/dune/solvers/transferoperators/genericmultigridtransfer.hh
@@ -8,11 +8,7 @@
 
 #include <dune/grid/common/genericreferenceelements.hh>
 
-#include <dune/localfunctions/common/virtualinterface.hh>
-#include <dune/localfunctions/common/virtualwrappers.hh>
-#include <dune/localfunctions/p1.hh>
-#include <dune/localfunctions/q1.hh>
-#include <dune/localfunctions/prismp1.hh>
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
 
 #include "dune/ag-common/staticmatrixtools.hh"
 
@@ -30,37 +26,6 @@
 //template<class DiscFuncType, class MatrixBlock, class TransferOperatorType>
 class GenericMultigridTransfer {
 
-    template <class DT, class RT, int dim>
-    struct P1ElementFactory {
-
-    private:
-        // extract LocalBasisTraits from P1LocalFiniteElement
-        typedef typename Dune::P1LocalFiniteElement<DT, RT, dim>::Traits::LocalBasisType::Traits P1LocalBasisTraits;
-
-        // make sure the traits have order 0 since not all FE's do implement higher order partial derivatives
-        typedef typename Dune::FixedOrderLocalBasisTraits<P1LocalBasisTraits, 0>::Traits LocalBasisTraits;
-
-    public:
-        typedef Dune::LocalFiniteElementVirtualInterface<LocalBasisTraits> type;
-
-        const type* get(const Dune::GeometryType& gt) {
-            if (gt.isSimplex())
-                return &simplexBaseSet_;
-            else if (gt.isCube())
-                return &cubeBaseSet_;
-            else if (gt.isPrism())
-                // This cast is necessary because otherwise the code wouldn't compile for dim!=3
-                return (type*)&prismBaseSet_;
-            else
-                DUNE_THROW(Dune::NotImplemented, "transfer operators for " << gt);
-        }
-
-    private:
-        Dune::LocalFiniteElementVirtualImp<Dune::P1LocalFiniteElement<DT, RT, dim> > simplexBaseSet_;
-        Dune::LocalFiniteElementVirtualImp<Dune::Q1LocalFiniteElement<DT, RT, dim> > cubeBaseSet_;
-        Dune::LocalFiniteElementVirtualImp<Dune::PrismP1LocalFiniteElement<DT, RT> > prismBaseSet_;
-    };
-
 public:
 
     template<class MatrixType, class VectorType>
@@ -119,8 +84,9 @@ public:
             identity[i][i] = 1;
 
         // A factory for the shape functions
-        P1ElementFactory<ctype,field_type,dim> p1ElementFactory;
-        typedef typename P1ElementFactory<ctype,field_type,dim>::type FEType;
+        typedef typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
+        typedef typename P1FECache::FiniteElementType FEType;
+        P1FECache p1FECache;
 
         matrix.setSize(rows,cols);
         matrix.setBuildMode(TransferOperatorType::random);
@@ -142,9 +108,9 @@ public:
             typedef typename EntityType::HierarchicIterator HierarchicIterator;
 
             // Get local finite element
-            const FEType* coarseBaseSet = p1ElementFactory.get(cIt->type());
+            const FEType& coarseBaseSet = p1FECache.get(cIt->type());
 
-            const int numCoarseBaseFct = coarseBaseSet->localBasis().size();
+            const int numCoarseBaseFct = coarseBaseSet.localBasis().size();
 
             // preallocate vector for function evaluations
             std::vector<Dune::FieldVector<field_type,1> > values(numCoarseBaseFct);
@@ -163,14 +129,14 @@ public:
                 const typename EntityType::Geometry& fGeometryInFather = fIt->geometryInFather();
 
                 // Get local finite element
-                const FEType* fineBaseSet = p1ElementFactory.get(fIt->type());;
+                const FEType& fineBaseSet = p1FECache.get(fIt->type());;
 
-                const int numFineBaseFct = fineBaseSet->localBasis().size();
+                const int numFineBaseFct = fineBaseSet.localBasis().size();
 
 
                 for (int j=0; j<numFineBaseFct; j++)
                 {
-                    const Dune::LocalKey& jLocalKey = fineBaseSet->localCoefficients().localKey(j);
+                    const Dune::LocalKey& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
 
                     int globalFine = fineIndexSet.subIndex(*fIt, jLocalKey.subEntity(), jLocalKey.codim());
 
@@ -178,13 +144,13 @@ public:
                     Dune::FieldVector<ctype, dim> local = fGeometryInFather.global(fineBasePosition);
 
                     // Evaluate coarse grid base functions
-                    coarseBaseSet->localBasis().evaluateFunction(local, values);
+                    coarseBaseSet.localBasis().evaluateFunction(local, values);
 
                     for (int i=0; i<numCoarseBaseFct; i++)
                     {
                         if (values[i] > 0.001) 
                         {
-                            const Dune::LocalKey& iLocalKey = coarseBaseSet->localCoefficients().localKey(i);
+                            const Dune::LocalKey& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
                             int globalCoarse = coarseIndexSet.subIndex(*cIt, iLocalKey.subEntity(), iLocalKey.codim());
                             indices.add(globalFine, globalCoarse);
                         }
@@ -204,9 +170,9 @@ public:
         for (; cIt != cEndIt; ++cIt) {
 
             // Get local finite element
-            const FEType* coarseBaseSet = p1ElementFactory.get(cIt->type());
+            const FEType& coarseBaseSet = p1FECache.get(cIt->type());
 
-            const int numCoarseBaseFct = coarseBaseSet->localBasis().size();
+            const int numCoarseBaseFct = coarseBaseSet.localBasis().size();
 
             typedef typename GridType::template Codim<0>::Entity EntityType;
             typedef typename EntityType::HierarchicIterator HierarchicIterator;
@@ -228,13 +194,13 @@ public:
                 const typename EntityType::Geometry& fGeometryInFather = fIt->geometryInFather();
 
                 // Get local finite element
-                const FEType* fineBaseSet = p1ElementFactory.get(fIt->type());
+                const FEType& fineBaseSet = p1FECache.get(fIt->type());
 
-                const int numFineBaseFct = fineBaseSet->localBasis().size();
+                const int numFineBaseFct = fineBaseSet.localBasis().size();
 
                 for (int j=0; j<numFineBaseFct; j++)
                 {
-                    const Dune::LocalKey& jLocalKey = fineBaseSet->localCoefficients().localKey(j);
+                    const Dune::LocalKey& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
 
                     int globalFine = fineIndexSet.subIndex(*fIt, jLocalKey.subEntity(), jLocalKey.codim());
 
@@ -242,13 +208,13 @@ public:
                     Dune::FieldVector<ctype, dim> local = fGeometryInFather.global(fineBasePosition);
 
                     // Evaluate coarse grid base functions
-                    coarseBaseSet->localBasis().evaluateFunction(local, values);
+                    coarseBaseSet.localBasis().evaluateFunction(local, values);
 
                     for (int i=0; i<numCoarseBaseFct; i++)
                     {
                         if (values[i] > 0.001) 
                         {
-                            const Dune::LocalKey& iLocalKey = coarseBaseSet->localCoefficients().localKey(i);
+                            const Dune::LocalKey& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
                             int globalCoarse = coarseIndexSet.subIndex(*cIt, iLocalKey.subEntity(), iLocalKey.codim());
 
                             TransferMatrixBlock matValue = identity;