diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 10f794906cb7fdd2d9a32959a8d3a090c9cdd392..af1fce745e7dafef0a3bfbb6b3db2b3a1d76627b 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -84,9 +84,6 @@ public:
     DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
                                       const TargetSpace& q) const;
 
-    /** \brief For debugging: Evaluate the derivative of the function using a finite-difference approximation*/
-    DerivativeType evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local) const;
-
     /** \brief Evaluate the derivative of the function value with respect to a coefficient */
     void evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
                                                  int coefficient,
@@ -301,34 +298,6 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local, const TargetSpace
     return result;
 }
 
-template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
-typename LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::DerivativeType
-LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
-evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local) const
-{
-    double eps = 1e-6;
-
-    Dune::FieldMatrix<ctype, embeddedDim, dim> result;
-
-    for (int i=0; i<dim; i++) {
-
-        Dune::FieldVector<ctype, dim> forward  = local;
-        Dune::FieldVector<ctype, dim> backward = local;
-
-        forward[i]  += eps;
-        backward[i] -= eps;
-
-        EmbeddedTangentVector fdDer = evaluate(forward).globalCoordinates() - evaluate(backward).globalCoordinates();
-        fdDer /= 2*eps;
-
-        for (int j=0; j<embeddedDim; j++)
-            result[j][i] = fdDer[j];
-
-    }
-
-    return result;
-}
-
 template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
 void LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
 evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
@@ -699,28 +668,6 @@ public:
         return result;
     }
 
-    /** \brief For debugging: Evaluate the derivative of the function using a finite-difference approximation*/
-    DerivativeType evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local) const
-    {
-        Dune::FieldMatrix<ctype, embeddedDim, dim> result(0);
-
-        // get translation part
-        std::vector<Dune::FieldMatrix<ctype,1,dim> > sfDer(translationCoefficients_.size());
-        localFiniteElement_.localBasis().evaluateJacobian(local, sfDer);
-
-        for (size_t i=0; i<translationCoefficients_.size(); i++)
-            for (int j=0; j<3; j++)
-                result[j].axpy(translationCoefficients_[i][j], sfDer[i][0]);
-
-        // get orientation part
-        Dune::FieldMatrix<ctype,4,dim> qResult = orientationFEFunction_->evaluateDerivativeFD(local);
-        for (int i=0; i<4; i++)
-            for (int j=0; j<dim; j++)
-                result[3+i][j] = qResult[i][j];
-
-        return result;
-    }
-
     /** \brief Evaluate the derivative of the function value with respect to a coefficient */
     void evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
                                                  int coefficient,
diff --git a/test/localgeodesicfefunctiontest.cc b/test/localgeodesicfefunctiontest.cc
index 41b2ed102d6a6949446225a47ee256958446907d..aea15a5580e3e36ef2a85cf7c732dc9a6fd4671b 100644
--- a/test/localgeodesicfefunctiontest.cc
+++ b/test/localgeodesicfefunctiontest.cc
@@ -34,6 +34,34 @@ double diameter(const std::vector<TargetSpace>& v)
     return d;
 }
 
+template <int dim, class ctype, class LocalFunction>
+auto
+evaluateDerivativeFD(const LocalFunction& f, const Dune::FieldVector<ctype, dim>& local)
+-> decltype(f.evaluateDerivative(local))
+{
+    double eps = 1e-8;
+    //static const int embeddedDim = LocalFunction::TargetSpace::embeddedDim;
+    decltype(f.evaluateDerivative(local)) result;
+
+    for (int i=0; i<dim; i++) {
+
+        Dune::FieldVector<ctype, dim> forward  = local;
+        Dune::FieldVector<ctype, dim> backward = local;
+
+        forward[i]  += eps;
+        backward[i] -= eps;
+
+        auto fdDer = f.evaluate(forward).globalCoordinates() - f.evaluate(backward).globalCoordinates();
+        fdDer /= 2*eps;
+
+        for (size_t j=0; j<result.N(); j++)
+            result[j][i] = fdDer[j];
+
+    }
+
+    return result;
+}
+
 
 template <int domainDim>
 void testDerivativeTangentiality(const RealTuple<double,1>& x,
@@ -155,7 +183,7 @@ void testDerivative(const LocalGeodesicFEFunction<domainDim,double,typename PQkL
         Dune::FieldMatrix<double, embeddedDim, domainDim> derivative = f.evaluateDerivative(quadPos);
 
         // evaluate fd approximation of derivative
-        Dune::FieldMatrix<double, embeddedDim, domainDim> fdDerivative = f.evaluateDerivativeFD(quadPos);
+        Dune::FieldMatrix<double, embeddedDim, domainDim> fdDerivative = evaluateDerivativeFD(f,quadPos);
 
         Dune::FieldMatrix<double, embeddedDim, domainDim> diff = derivative;
         diff -= fdDerivative;