diff --git a/src/localgeodesicfefunction.hh b/src/localgeodesicfefunction.hh
index 6b63c8777b40d7a40e92ca5dd37d2f292f6aeeb9..f465e3fa8ac41790666b8026fb7924eacf613386 100644
--- a/src/localgeodesicfefunction.hh
+++ b/src/localgeodesicfefunction.hh
@@ -130,10 +130,10 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local)
 
     // the actual system matrix
     Dune::FieldVector<ctype,dim+1> w = barycentricCoordinates(local);
-
+    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+    
     Dune::FieldMatrix<ctype,dim+1,dim+1> dFdq(0);
-    for (int i=0; i<dim+1; i++)
-        dFdq.axpy(w[i], TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], q));
+    assembler.assembleHessian(q,dFdq);
 
     // ////////////////////////////////////
     //   solve the system