diff --git a/src/Makefile.am b/src/Makefile.am
index 3a4555856fae0c9d5d9c2f1808d446f0d5d37f46..36a2c8cd68bd38999c935838cd4d0d8ec1823fa2 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -6,7 +6,7 @@ AM_CPPFLAGS = @AM_CPPFLAGS@ -I$(top_srcdir)/..
 
 srcincludedir = $(includedir)/dune/common
 srcinclude_HEADERS = averagedistanceassembler.hh localgeodesicfefunction.hh \
-                     quaternion.hh rodrefine.hh averageinterface.hh localstiffness.hh \
+                     quaternion.hh geodesicfefunctionadaptor.hh averageinterface.hh localstiffness.hh \
                      localgeodesicfestiffness.hh riemanniantrsolver.hh rodwriter.hh \
                      geodesicdifference.hh makestraightrod.hh rigidbodymotion.hh \
                      rotation.hh geodesicfeassembler.hh maxnormtrustregion.hh \
diff --git a/src/geodesicfefunctionadaptor.hh b/src/geodesicfefunctionadaptor.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b9340a6755c9682a7eab26441824aa5577a9e306
--- /dev/null
+++ b/src/geodesicfefunctionadaptor.hh
@@ -0,0 +1,85 @@
+#ifndef GEODESIC_FE_FUNCTION_ADAPTOR_HH
+#define GEODESIC_FE_FUNCTION_ADAPTOR_HH
+
+#include <vector>
+#include <map>
+
+#include "localgeodesicfefunction.hh"
+
+template <class GridType, class TargetSpace>
+void geodesicFEFunctionAdaptor(GridType& grid, std::vector<TargetSpace>& x)
+{
+    const int dim = GridType::dimension;
+
+    typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
+    typedef typename GridType::template Codim<dim>::LeafIterator VertexIterator;
+
+    // /////////////////////////////////////////////////////
+    //   Save leaf p1 data in a map
+    // /////////////////////////////////////////////////////
+
+    const typename GridType::Traits::LocalIdSet&   idSet    = grid.localIdSet();
+    const typename GridType::Traits::LeafIndexSet& indexSet = grid.leafIndexSet();
+
+    std::map<typename GridType::Traits::LocalIdSet::IdType, TargetSpace> dofMap;
+
+    VertexIterator vIt    = grid.template leafbegin<dim>();
+    VertexIterator vEndIt = grid.template leafend<dim>();
+
+    for (; vIt!=vEndIt; ++vIt)
+        dofMap.insert(std::make_pair(idSet.id(*vIt), x[indexSet.index(*vIt)]));
+
+
+
+    // /////////////////////////////////////////////////////
+    //   Globally refine the grid
+    // /////////////////////////////////////////////////////
+
+    grid.globalRefine(1);
+
+
+    // /////////////////////////////////////////////////////
+    //   Restore and interpolate the data
+    // /////////////////////////////////////////////////////
+
+    x.resize(grid.size(dim));
+
+    ElementIterator eIt    = grid.template leafbegin<0>();
+    ElementIterator eEndIt = grid.template leafend<0>();
+
+    for (; eIt!=eEndIt; ++eIt) {
+
+        // Set up a local gfe function on the father element
+        std::vector<TargetSpace> coefficients(eIt->father()->template count<dim>());
+
+        for (int i=0; i<eIt->father()->template count<dim>(); i++)
+            coefficients[i] = dofMap.find(idSet.subId(*eIt->father(),i,dim))->second;
+
+        LocalGeodesicFEFunction<dim,double,TargetSpace> fatherFunction(coefficients);
+
+        // The embedding of this element into the father geometry
+        const typename GridType::template Codim<0>::LocalGeometry& geometryInFather = eIt->geometryInFather();
+
+        for (int i=0; i<eIt->template count<dim>(); i++) {
+
+            if (dofMap.find(idSet.subId(*eIt,i,dim)) != dofMap.end()) {
+
+                // If the vertex exists on the coarser level we take the value from there.
+                // That should be faster and more accurate than interpolating
+                x[indexSet.subIndex(*eIt,i,dim)] = dofMap[idSet.subId(*eIt,i,dim)];
+
+            } else {
+
+                // Interpolate
+                x[indexSet.subIndex(*eIt,i,dim)] = fatherFunction.evaluate(geometryInFather.corner(i));
+
+            }
+
+        }
+
+    }
+
+
+}
+
+#endif
diff --git a/src/rodrefine.hh b/src/rodrefine.hh
deleted file mode 100644
index d847035cbfc1dfd4821311be14d78c401eff31c4..0000000000000000000000000000000000000000
--- a/src/rodrefine.hh
+++ /dev/null
@@ -1,78 +0,0 @@
-#ifndef ROD_REFINE_HH
-#define ROD_REFINE_HH
-
-#include <vector>
-#include <map>
-
-#include "rigidbodymotion.hh"
-
-template <class GridType>
-void globalRodRefine(GridType& grid, std::vector<RigidBodyMotion<3> >& x)
-{
-    // Should be a static assertion
-    assert(GridType::dimension == 1);
-
-    typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
-    typedef typename GridType::template Codim<1>::LeafIterator VertexIterator;
-
-    // /////////////////////////////////////////////////////
-    //   Save leaf p1 data in a map
-    // /////////////////////////////////////////////////////
-
-    const typename GridType::Traits::LocalIdSet&    idSet    = grid.localIdSet();
-    const typename GridType::Traits::LeafIndexSet& indexSet = grid.leafIndexSet();
-
-    std::map<typename GridType::Traits::LocalIdSet::IdType, RigidBodyMotion<3> > dofMap;
-
-    VertexIterator vIt    = grid.template leafbegin<1>();
-    VertexIterator vEndIt = grid.template leafend<1>();
-
-    for (; vIt!=vEndIt; ++vIt)
-        dofMap.insert(std::make_pair(idSet.id(*vIt), x[indexSet.index(*vIt)]));
-
-
-
-    // /////////////////////////////////////////////////////
-    //   Globally refine the grid
-    // /////////////////////////////////////////////////////
-
-    grid.globalRefine(1);
-
-
-    // /////////////////////////////////////////////////////
-    //   Restore and interpolate the data
-    // /////////////////////////////////////////////////////
-
-    x.resize(grid.size(1));
-
-    ElementIterator eIt    = grid.template leafbegin<0>();
-    ElementIterator eEndIt = grid.template leafend<0>();
-
-    for (; eIt!=eEndIt; ++eIt) {
-
-        for (int i=0; i<2; i++) {
-
-            if (dofMap.find(idSet.subId(*eIt,i,1)) != dofMap.end()) {
-
-                x[indexSet.subIndex(*eIt,i,1)] = dofMap[idSet.subId(*eIt,i,1)];
-
-            } else {
-                // Interpolate
-                RigidBodyMotion<3> p0 = dofMap[idSet.subId(*eIt->father(),0,1)];
-                RigidBodyMotion<3> p1 = dofMap[idSet.subId(*eIt->father(),1,1)];
-
-                x[indexSet.subIndex(*eIt,i,1)].r = (p0.r + p1.r);
-                x[indexSet.subIndex(*eIt,i,1)].r *= 0.5;
-                x[indexSet.subIndex(*eIt,i,1)].q
-                    = Rotation<3,double>::interpolate(p0.q, p1.q, 0.5);
-
-            }
-
-        }
-
-    }
-
-
-}
-
-#endif