diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh
index 0f9cb017c0f9810bd221f3442412a53312ca74fe..d718cbed5c4ca67f1623322cf0c2a364346db887 100644
--- a/dune/gfe/localgeodesicfestiffness.hh
+++ b/dune/gfe/localgeodesicfestiffness.hh
@@ -220,6 +220,10 @@ public:
                                   const std::vector<TargetSpace>& solution,
                                   std::vector<Dune::FieldVector<double,blocksize> >& gradient) const;
 
+    virtual void assembleEmbeddedFDGradient(const Entity& element,
+                                  const std::vector<TargetSpace>& solution,
+                                  std::vector<typename TargetSpace::EmbeddedTangentVector>& gradient) const;
+
     // assembled data
     Dune::Matrix<Dune::FieldMatrix<double,blocksize,blocksize> > A_;
     
@@ -234,6 +238,14 @@ assembleGradient(const Entity& element,
     LocalGeodesicFEStiffnessImp<GridView,TargetSpace,globalIsometricCoordinates>::assembleGradient(element, localSolution, localGradient, this);
 }
 
+template <class GridView, class TargetSpace>
+void LocalGeodesicFEStiffness<GridView, TargetSpace>::
+assembleEmbeddedFDGradient(const Entity& element,
+                 const std::vector<TargetSpace>& localSolution,
+                 std::vector<typename TargetSpace::EmbeddedTangentVector>& localGradient) const
+{
+    LocalGeodesicFEStiffnessImp<GridView,TargetSpace,globalIsometricCoordinates>::assembleEmbeddedGradient(element, localSolution, localGradient, this);
+}
 
 template <class GridView, class TargetSpace>
 void LocalGeodesicFEStiffness<GridView,TargetSpace>::
@@ -423,14 +435,13 @@ public:
     virtual RT energy (const Entity& e,
                        const std::vector<TargetSpace>& localSolution) const = 0;
 
-#if 0
-    /** \brief Assemble the element gradient of the energy functional 
+    /** \brief Assemble the element gradient of the energy functional using a finite-difference approximation
 
-    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;
-#endif
+    This is mainly for debugging purposes.
+    */
+    virtual void assembleEmbeddedFDGradient(const Entity& element,
+                                            const std::vector<TargetSpace>& solution,
+                                            std::vector<typename TargetSpace::EmbeddedTangentVector>& gradient) const;
                                           
     /** \brief Assemble the element gradient of the energy functional 
 
@@ -507,6 +518,16 @@ assembleGradient(const Entity& element,
     LocalGeodesicFEStiffnessImp<GridView,TargetSpace,globalIsometricCoordinates>::assembleGradient(element, localSolution, localGradient,this);
 }
 
+template <class GridView, int dim>
+void LocalGeodesicFEStiffness<GridView, UnitVector<dim> >::
+assembleEmbeddedFDGradient(const Entity& element,
+                 const std::vector<TargetSpace>& localSolution,
+                 std::vector<typename TargetSpace::EmbeddedTangentVector>& localGradient) const
+{
+    LocalGeodesicFEStiffnessImp<GridView,TargetSpace,globalIsometricCoordinates>::assembleEmbeddedGradient(element, localSolution, localGradient, this);
+}
+
+
 // ///////////////////////////////////////////////////////////
 //   Compute gradient by finite-difference approximation
 // ///////////////////////////////////////////////////////////