diff --git a/dune/gfe/geodesicfeassembler.hh b/dune/gfe/geodesicfeassembler.hh
index 9f55c197080d3036bf99709b21f2bc0737dbd784..bf1593c94c53608f9b850f05e103a78a9e627e07 100644
--- a/dune/gfe/geodesicfeassembler.hh
+++ b/dune/gfe/geodesicfeassembler.hh
@@ -22,7 +22,7 @@ class GeodesicFEAssembler {
     enum { gridDim = GridView::dimension };
     
     //! Dimension of a tangent space
-    enum { blocksize = TargetSpace::EmbeddedTangentVector::size };
+    enum { blocksize = TargetSpace::TangentVector::size };
     
     //!
     typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh
index 18e9884c33e50397845914617a4500ba36815c78..0558aba249c20f1c2ab31e5fbbe2628217c75c2f 100644
--- a/dune/gfe/localgeodesicfestiffness.hh
+++ b/dune/gfe/localgeodesicfestiffness.hh
@@ -316,7 +316,7 @@ class LocalGeodesicFEStiffness <GridView,UnitVector<dim> >
 public:
     
     //! Each block is x, y, theta in 2d, T (R^3 \times SO(3)) in 3d
-    enum { blocksize = TargetSpace::EmbeddedTangentVector::size };
+    enum { blocksize = TargetSpace::TangentVector::size };
 
     /** \brief Assemble the local stiffness matrix at the current position
 
@@ -330,10 +330,17 @@ public:
 
     /** \brief Assemble the element gradient of the energy functional 
 
+    The default implementation in this class uses a finite difference approximation */
+    virtual void assembleEmbeddedGradient(const Entity& element,
+                                          const std::vector<TargetSpace>& solution,
+                                          std::vector<typename TargetSpace::EmbeddedTangentVector>& gradient) const;
+
+    /** \brief Assemble the element gradient of the energy functional 
+
     The default implementation in this class uses a finite difference approximation */
     virtual void assembleGradient(const Entity& element,
                                   const std::vector<TargetSpace>& solution,
-                                  std::vector<Dune::FieldVector<double,blocksize> >& gradient) const;
+                                  std::vector<typename TargetSpace::TangentVector>& gradient) const;
 
     // assembled data
     Dune::Matrix<Dune::FieldMatrix<double,blocksize,blocksize> > A_;
@@ -342,11 +349,13 @@ public:
 
 template <class GridView, int dim>
 void LocalGeodesicFEStiffness<GridView, UnitVector<dim> >::
-assembleGradient(const Entity& element,
-                 const std::vector<TargetSpace>& localSolution,
-                 std::vector<Dune::FieldVector<double,blocksize> >& localGradient) const
+assembleEmbeddedGradient(const Entity& element,
+                         const std::vector<TargetSpace>& localSolution,
+                         std::vector<typename TargetSpace::EmbeddedTangentVector>& localGradient) const
 {
 
+    const int embeddedBlocksize = TargetSpace::EmbeddedTangentVector::size;
+
     // ///////////////////////////////////////////////////////////
     //   Compute gradient by finite-difference approximation
     // ///////////////////////////////////////////////////////////
@@ -360,7 +369,7 @@ assembleGradient(const Entity& element,
 
     for (size_t i=0; i<localSolution.size(); i++) {
         
-        for (int j=0; j<blocksize; j++) {
+        for (int j=0; j<embeddedBlocksize; j++) {
             
             // The return value does not have unit norm.  But assigning it to a UnitVector object
             // will normalize it.  This amounts to an extension of the energy functional 
@@ -382,6 +391,24 @@ assembleGradient(const Entity& element,
 
 }
 
+template <class GridView, int dim>
+void LocalGeodesicFEStiffness<GridView, UnitVector<dim> >::
+assembleGradient(const Entity& element,
+                 const std::vector<TargetSpace>& localSolution,
+                 std::vector<typename TargetSpace::TangentVector>& localGradient) const
+{
+    std::vector<typename TargetSpace::EmbeddedTangentVector> embeddedLocalGradient;
+
+    // first compute the gradient in embedded coordinates
+    assembleEmbeddedGradient(element, localSolution, embeddedLocalGradient);
+
+    // transform to coordinates on the tangent space
+    localGradient.resize(embeddedLocalGradient.size());
+
+    for (size_t i=0; i<localGradient.size(); i++)
+        localSolution[i].orthonormalFrame().mv(embeddedLocalGradient[i], localGradient[i]);
+
+}
 
 template <class GridType, int dim>
 void LocalGeodesicFEStiffness<GridType,UnitVector<dim> >::
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index ba09eaa6416513565a57981a06eaae0c9704b90a..fac86e406132c0fd04cb48e78a09e3f77118aa0c 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -18,9 +18,9 @@
 template <class GridType, class TargetSpace>
 class RiemannianTrustRegionSolver 
     : public IterativeSolver<std::vector<TargetSpace>,
-                             Dune::BitSetVector<TargetSpace::EmbeddedTangentVector::size> >
+                             Dune::BitSetVector<TargetSpace::TangentVector::size> >
 { 
-    const static int blocksize = TargetSpace::EmbeddedTangentVector::size;
+    const static int blocksize = TargetSpace::TangentVector::size;
 
     const static int gridDim = GridType::dimension;
 
diff --git a/harmonicmaps.cc b/harmonicmaps.cc
index 7d529c9d17fae47b4fc3654d357268595333d351..d276de77ae628c87b732e78f8c2e3629d97f76b5 100644
--- a/harmonicmaps.cc
+++ b/harmonicmaps.cc
@@ -45,7 +45,7 @@ typedef RealTuple<1> TargetSpace;
 #endif
 
 // Tangent vector of the image space
-const int blocksize = TargetSpace::EmbeddedTangentVector::size;
+const int blocksize = TargetSpace::TangentVector::size;
 
 using namespace Dune;
 
@@ -254,6 +254,8 @@ int main (int argc, char *argv[]) try
     writeRod(x, resultPath + "rod3d.result");
 #endif
 
+    exit(0);
+
     // //////////////////////////////////////////////////////////
     //   Recompute and compare against exact solution
     // //////////////////////////////////////////////////////////