diff --git a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
index 66b8f85eff794ac87c2a533f64834953357ac4ca..ba0830b58373be5faac8dfe7263db44103a9da26 100644
--- a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
+++ b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
@@ -559,7 +559,7 @@ linearizedRodNeumannToDirichletMap(const std::string& rodName,
     const RodLeafBoundaryPatch& dirichletBoundary = this->complex_.rod(rodName).dirichletBoundary_;
 
     typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,6,6> > MatrixType;
-    GeodesicFEAssembler<typename RodGridType::LeafGridView, RigidBodyMotion<3> > assembler(dirichletBoundary.gridView(),
+    GeodesicFEAssembler<P1NodalBasis<typename RodGridType::LeafGridView>, RigidBodyMotion<3> > assembler(dirichletBoundary.gridView(),
                                                                     rod(rodName).localStiffness_);
     MatrixType stiffnessMatrix;
     assembler.assembleMatrix(currentIterate, 
diff --git a/dune/gfe/rodassembler.cc b/dune/gfe/rodassembler.cc
index dc356b2e0f90d4a5d58809bb9c6dc7c2c38cc871..6ab87a730b25587e38821837790817e25b5faf3f 100644
--- a/dune/gfe/rodassembler.cc
+++ b/dune/gfe/rodassembler.cc
@@ -18,7 +18,7 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
 {
     using namespace Dune;
 
-    const typename GridView::Traits::IndexSet& indexSet = this->gridView_.indexSet();
+    const typename GridView::Traits::IndexSet& indexSet = this->basis_.getGridView().indexSet();
 
     if (sol.size()!=indexSet.size(gridDim))
         DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
@@ -26,8 +26,8 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
     grad.resize(sol.size());
     grad = 0;
 
-    ElementIterator it    = this->gridView_.template begin<0>();
-    ElementIterator endIt = this->gridView_.template end<0>();
+    ElementIterator it    = this->basis_.getGridView().template begin<0>();
+    ElementIterator endIt = this->basis_.getGridView().template end<0>();
 
     // Loop over all elements
     for (; it!=endIt; ++it) {
@@ -36,7 +36,7 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
         static const int nDofs = 2;
 
         // Extract local solution
-        Dune::array<RigidBodyMotion<3>, nDofs> localSolution;
+        std::vector<RigidBodyMotion<3> > localSolution(nDofs);
         
         for (int i=0; i<nDofs; i++)
             localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
@@ -44,7 +44,10 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
         // Assemble local gradient
         std::vector<FieldVector<double,blocksize> > localGradient(nDofs);
 
-        this->localStiffness_->assembleGradient(*it, localSolution, localGradient);
+        this->localStiffness_->assembleGradient(*it, 
+                                                this->basis_.getLocalFiniteElement(*it),
+                                                localSolution, 
+                                                localGradient);
 
         // Add to global gradient
         for (int i=0; i<nDofs; i++)
@@ -150,7 +153,7 @@ getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
     //    if (gridView_ != &boundary.gridView())
     //        DUNE_THROW(Dune::Exception, "The boundary patch has to match the grid view of the assembler!");
 
-    const typename GridView::Traits::IndexSet& indexSet = this->gridView_.indexSet();
+    const typename GridView::Traits::IndexSet& indexSet = this->basis_.getGridView().indexSet();
 
     if (sol.size()!=indexSet.size(gridDim))
         DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
diff --git a/dune/gfe/rodassembler.hh b/dune/gfe/rodassembler.hh
index cd49bd4276745fda70a3a08261405f9951d96f4d..77fac7d8fc2a7101ea964dabdc2d9191b7f74c9f 100644
--- a/dune/gfe/rodassembler.hh
+++ b/dune/gfe/rodassembler.hh
@@ -24,7 +24,7 @@ class RodAssembler
 /** \brief The FEM operator for an extensible, shearable rod in 3d
  */
 template <class GridView>
-class RodAssembler<GridView,3> : public GeodesicFEAssembler<GridView, RigidBodyMotion<3> >
+class RodAssembler<GridView,3> : public GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<3> >
 {
         
     //typedef typename GridType::template Codim<0>::Entity EntityType;
@@ -46,7 +46,7 @@ public:
         //! ???
     RodAssembler(const GridView &gridView,
                  RodLocalStiffness<GridView,double>* localStiffness) 
-        : GeodesicFEAssembler<GridView, RigidBodyMotion<3> >(gridView,localStiffness)
+        : GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<3> >(gridView,localStiffness)
         { 
             std::vector<RigidBodyMotion<3> > referenceConfiguration(gridView.size(gridDim));
 
diff --git a/dune/gfe/rodlocalstiffness.hh b/dune/gfe/rodlocalstiffness.hh
index ac64fff625e039e0a82435a15827257d9ececc98..89296c027c2027b4fca3458c80d30fe6caac8f21 100644
--- a/dune/gfe/rodlocalstiffness.hh
+++ b/dune/gfe/rodlocalstiffness.hh
@@ -11,7 +11,7 @@
 
 template<class GridView, class RT>
 class RodLocalStiffness 
-    : public LocalGeodesicFEStiffness<GridView,RigidBodyMotion<3> >
+    : public LocalGeodesicFEStiffness<GridView, typename P1NodalBasis<GridView>::LocalFiniteElement, RigidBodyMotion<3> >
 {
     typedef RigidBodyMotion<3> TargetSpace;
 
@@ -95,9 +95,19 @@ public:
         referenceConfiguration_ = referenceConfiguration;
     }
     
+    /** \brief Local element energy for a P1 element */
     virtual RT energy (const Entity& e,
                        const Dune::array<RigidBodyMotion<3>, dim+1>& localSolution) const;
 
+    virtual RT energy (const Entity& e,
+                       const typename P1NodalBasis<GridView>::LocalFiniteElement& localFiniteElement,
+                       const std::vector<RigidBodyMotion<3> >& localSolution) const
+    {
+        assert(localSolution.size()==2);
+        Dune::array<RigidBodyMotion<3>, 2> localSolutionArray = {localSolution[0], localSolution[1]};
+        return energy(e,localSolutionArray);
+    }
+
     /** \brief Assemble the element gradient of the energy functional */
     void assembleGradient(const Entity& element,
                           const std::vector<RigidBodyMotion<3> >& solution,