diff --git a/dune/gfe/averagedistanceassembler.hh b/dune/gfe/averagedistanceassembler.hh
index 33319bc7ab4bd39b76d02aa2815e2b2d465954ae..bd675e9e1580ca7a896f484c62a9d96262020a8a 100644
--- a/dune/gfe/averagedistanceassembler.hh
+++ b/dune/gfe/averagedistanceassembler.hh
@@ -3,8 +3,8 @@
 
 #include <vector>
 
-/** \tparam N Number of coefficients (i.e., simplex corners) */
-template <class TargetSpace, int N>
+/** \tparam TargetSpace The manifold that we are mapping to */
+template <class TargetSpace>
 class AverageDistanceAssembler
 {
     static const int size         = TargetSpace::TangentVector::dimension;
@@ -12,12 +12,8 @@ class AverageDistanceAssembler
 
 public:
 
-    AverageDistanceAssembler(const Dune::array<TargetSpace,N>& coefficients,
-                             const Dune::array<double,N>& weights) DUNE_DEPRECATED
-        : coefficients_(coefficients.begin(), coefficients.end()),
-          weights_(weights.begin(), weights.end())
-    {}
-
+    /** \brief Constructor with given coefficients \f$ v_i \f$ and weights \f$ w_i \f$
+     */
     AverageDistanceAssembler(const std::vector<TargetSpace>& coefficients,
                              const std::vector<double>& weights)
         : coefficients_(coefficients),
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 43ee577b77a302046b63322f0c7e4beaa16025e8..5f4f729996b10db302598437fe5fc7cad09fbe39 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -156,10 +156,10 @@ evaluate(const Dune::FieldVector<ctype, dim>& local) const
     for (size_t i=0; i<w.size(); i++)
         w[i] = wNested[i][0];
 
-    /** \todo The 'dim+1' parameter is a dummy here */
-    AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
+    // The energy functional whose mimimizer is the value of the geodesic interpolation
+    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
 
-    TargetSpaceRiemannianTRSolver<TargetSpace,dim+1> solver;
+    TargetSpaceRiemannianTRSolver<TargetSpace> solver;
 
     solver.setup(&assembler,
                  coefficients_[0],   // initial iterate
@@ -223,7 +223,7 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
     for (size_t i=0; i<w.size(); i++)
         w[i] = wNested[i][0];
     
-    AverageDistanceAssembler<TargetSpace,1> assembler(coefficients_, w);
+    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
     
     Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
     assembler.assembleEmbeddedHessian(q,dFdq);
@@ -314,8 +314,14 @@ evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& loc
     TargetSpace q = evaluate(local);
 
     // dFdq
-    Dune::array<ctype,dim+1> w = barycentricCoordinates(local);
-    AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
+    std::vector<Dune::FieldVector<ctype,1> > wNested;
+    localFiniteElement_.localBasis().evaluateFunction(local,wNested);
+
+    std::vector<ctype> w(wNested.size());
+    for (size_t i=0; i<w.size(); i++)
+        w[i] = wNested[i][0];
+
+    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
     
     Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
     assembler.assembleEmbeddedHessian(q,dFdq);
@@ -429,7 +435,7 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
     for (size_t i=0; i<w.size(); i++)
         w[i] = wNested[i][0];
 
-    AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
+    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
     
     Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
     assembler.assembleEmbeddedHessian(q,dFdq);
diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index aa2a8ade6b03c4082c9e5fd176f95390e6b70c28..d59f290067b0be7bba9456d1a340e3da1d4983c4 100644
--- a/dune/gfe/targetspacertrsolver.cc
+++ b/dune/gfe/targetspacertrsolver.cc
@@ -6,9 +6,9 @@
 
 #include "maxnormtrustregion.hh"
 
-template <class TargetSpace, int N>
-void TargetSpaceRiemannianTRSolver<TargetSpace,N>::
-setup(const AverageDistanceAssembler<TargetSpace,N>* assembler,
+template <class TargetSpace>
+void TargetSpaceRiemannianTRSolver<TargetSpace>::
+setup(const AverageDistanceAssembler<TargetSpace>* assembler,
       const TargetSpace& x,
       double tolerance,
       int maxTrustRegionSteps,
@@ -45,8 +45,8 @@ setup(const AverageDistanceAssembler<TargetSpace,N>* assembler,
 }
 
 
-template <class TargetSpace, int N>
-void TargetSpaceRiemannianTRSolver<TargetSpace,N>::solve()
+template <class TargetSpace>
+void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
 {
     MaxNormTrustRegion<blocksize> trustRegion(1,   // we have only one block
                                               initialTrustRegionRadius_);
diff --git a/dune/gfe/targetspacertrsolver.hh b/dune/gfe/targetspacertrsolver.hh
index b721ae7454080bbd0c9145576170b69eef3c6d47..16e776a8bf0561f416b0e2f869b88c394f1c443f 100644
--- a/dune/gfe/targetspacertrsolver.hh
+++ b/dune/gfe/targetspacertrsolver.hh
@@ -10,9 +10,8 @@
 
 /** \brief Riemannian trust-region solver for geodesic finite-element problems 
  \tparam TargetSpace The manifold that our functions take values in
- \tparam N Number of simplex vertices
  */
-template <class TargetSpace, int N>
+template <class TargetSpace>
 class TargetSpaceRiemannianTRSolver 
     : public NumProc
 { 
@@ -34,7 +33,7 @@ public:
     {}
 
     /** \brief Set up the solver using a monotone multigrid method as the inner solver */
-    void setup(const AverageDistanceAssembler<TargetSpace,N>* assembler,
+    void setup(const AverageDistanceAssembler<TargetSpace>* assembler,
                const TargetSpace& x,
                double tolerance,
                int maxTrustRegionSteps,
@@ -71,7 +70,7 @@ protected:
     double innerTolerance_;
 
     /** \brief The assembler for the average-distance functional */
-    const AverageDistanceAssembler<TargetSpace,N>* assembler_;
+    const AverageDistanceAssembler<TargetSpace>* assembler_;
 
     /** \brief The solver for the quadratic inner problems */
     std::auto_ptr< ::LoopSolver<CorrectionType> > innerSolver_;