diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh
index 186b454c5d7459b74c2cd7d89734817bda2c2313..ca0205e0ab4d49532a6d00f9cd8f8bf4415929b6 100644
--- a/dune/gfe/localgeodesicfestiffness.hh
+++ b/dune/gfe/localgeodesicfestiffness.hh
@@ -27,22 +27,12 @@ public:
     //! Dimension of the embedding space
     enum { embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension };
 
-    /** \brief Assemble the local stiffness matrix at the current position
-
-    This default implementation used finite-difference approximations to compute the second derivatives
-
-    The formula for the Riemannian Hessian has been taken from Absil, Mahony, Sepulchre:
-    'Optimization algorithms on matrix manifolds', page 107.  There it says that
-    \f[
-        \langle Hess f(x)[\xi], \eta \rangle
-            = \frac 12 \frac{d^2}{dt^2} \Big(f(\exp_x(t(\xi + \eta))) - f(\exp_x(t\xi)) - f(\exp_x(t\eta))\Big)\Big|_{t=0}.
-    \f]
-    We compute that using a finite difference approximation.
+    /** \brief Assemble the local gradient and stiffness matrix at the current position
 
     */
     virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
                                  const std::vector<TargetSpace>& localSolution,
-                                 std::vector<typename TargetSpace::TangentVector>& localGradient);
+                                 std::vector<typename TargetSpace::TangentVector>& localGradient) = 0;
 
     /** \brief Compute the energy at the current configuration */
     virtual RT energy (const typename Basis::LocalView& localView,
@@ -53,35 +43,12 @@ public:
     The default implementation in this class uses a finite difference approximation */
     virtual void assembleGradient(const typename Basis::LocalView& localView,
                                   const std::vector<TargetSpace>& solution,
-                                  std::vector<typename TargetSpace::TangentVector>& gradient) const;
+                                  std::vector<typename TargetSpace::TangentVector>& gradient) const = 0;
 
     // assembled data
     Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> > A_;
 
 };
 
-
-template <class Basis, class TargetSpace>
-void LocalGeodesicFEStiffness<Basis, TargetSpace>::
-assembleGradient(const typename Basis::LocalView& localView,
-                 const std::vector<TargetSpace>& localSolution,
-                 std::vector<typename TargetSpace::TangentVector>& localGradient) const
-{
-  DUNE_THROW(Dune::NotImplemented, "!");
-}
-
-
-// ///////////////////////////////////////////////////////////
-//   Compute gradient by finite-difference approximation
-// ///////////////////////////////////////////////////////////
-template <class Basis, class TargetSpace>
-void LocalGeodesicFEStiffness<Basis, TargetSpace>::
-assembleGradientAndHessian(const typename Basis::LocalView& localView,
-                const std::vector<TargetSpace>& localSolution,
-                std::vector<typename TargetSpace::TangentVector>& localGradient)
-{
-  DUNE_THROW(Dune::NotImplemented, "!");
-}
-
 #endif