diff --git a/dune/gfe/averagedistanceassembler.hh b/dune/gfe/averagedistanceassembler.hh
index a940fe650bf04ad407374a8d96bea4d80225c7dd..e92ac0cef3df9f07a73ee24e7bc53e10962c98cb 100644
--- a/dune/gfe/averagedistanceassembler.hh
+++ b/dune/gfe/averagedistanceassembler.hh
@@ -3,7 +3,8 @@
 
 #include <vector>
 
-template <class TargetSpace>
+/** \tparam N Number of coefficients (i.e., simplex corners) */
+template <class TargetSpace, int N>
 class AverageDistanceAssembler
 {
     static const int size         = TargetSpace::TangentVector::size;
@@ -11,8 +12,8 @@ class AverageDistanceAssembler
 
 public:
 
-    AverageDistanceAssembler(const std::vector<TargetSpace>& coefficients,
-                             const std::vector<double>& weights)
+    AverageDistanceAssembler(const Dune::array<TargetSpace,N>& coefficients,
+                             const Dune::array<double,N>& weights)
         : coefficients_(coefficients),
           weights_(weights)
     {}
@@ -74,9 +75,9 @@ public:
         
     }
 
-    const std::vector<TargetSpace> coefficients_;
+    const Dune::array<TargetSpace,N> coefficients_;
 
-    const std::vector<double> weights_;
+    const Dune::array<double,N> weights_;
 
 };
 
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index b8f53ac64fb74ca6784fe8f1bbeae3086a758d2c..4dcc4d5f98d09f819da42e05b0e62b6efc5d5e1a 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -95,7 +95,7 @@ public:
 
     /** \brief Assemble the energy for a single element */
     RT energy (const Entity& e,
-               const std::vector<TargetSpace>& localSolution) const;
+               const Dune::array<TargetSpace, gridDim+1>& localSolution) const;
                
     RT quadraticMembraneEnergy(const Dune::FieldMatrix<double,3,3>& U) const
     {
@@ -162,7 +162,7 @@ public:
 template <class GridView, int dim>
 typename CosseratEnergyLocalStiffness<GridView, dim>::RT CosseratEnergyLocalStiffness<GridView, dim>::
 energy(const Entity& element,
-       const std::vector<RigidBodyMotion<dim> >& localSolution) const
+       const Dune::array<RigidBodyMotion<dim>, gridDim+1>& localSolution) const
 {
     RT energy = 0;
 
diff --git a/dune/gfe/geodesicfeassembler.hh b/dune/gfe/geodesicfeassembler.hh
index 2a13170ad9641e9c25ee4f8fe624b2637fa0d38d..288d07e0853f80be9ab5a269eee9eada461e8e32 100644
--- a/dune/gfe/geodesicfeassembler.hh
+++ b/dune/gfe/geodesicfeassembler.hh
@@ -120,7 +120,7 @@ assembleMatrix(const std::vector<TargetSpace>& sol,
         const int numOfBaseFct = it->template count<gridDim>();  
         
         // Extract local solution
-        std::vector<TargetSpace> localSolution(numOfBaseFct);
+        Dune::array<TargetSpace,gridDim+1> localSolution;
         
         for (int i=0; i<numOfBaseFct; i++)
             localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
@@ -168,7 +168,7 @@ assembleGradient(const std::vector<TargetSpace>& sol,
         const int nDofs = it->template count<gridDim>();
 
         // Extract local solution
-        std::vector<TargetSpace> localSolution(nDofs);
+        Dune::array<TargetSpace,gridDim+1> localSolution;
         
         for (int i=0; i<nDofs; i++)
             localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
@@ -198,7 +198,7 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
     if (sol.size()!=indexSet.size(gridDim))
         DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
 
-    std::vector<TargetSpace> localSolution;
+    Dune::array<TargetSpace,gridDim+1> localSolution;
 
     ElementIterator it    = gridView_.template begin<0>();
     ElementIterator endIt = gridView_.template end<0>();
@@ -206,8 +206,6 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
     // Loop over all elements
     for (; it!=endIt; ++it) {
 
-        localSolution.resize(it->template count<gridDim>());
-
         for (int i=0; i<it->template count<gridDim>(); i++)
             localSolution[i]               = sol[indexSet.subIndex(*it,i,gridDim)];
 
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index bf7f0f670e0f5990eab25fde3949f2947a2e4c89..2a6760e2390a8a6059416f16d31c3a44daf95767 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -31,11 +31,9 @@ class LocalGeodesicFEFunction
 public:
 
     /** \brief Constructor */
-    LocalGeodesicFEFunction(const std::vector<TargetSpace>& coefficients)
+    LocalGeodesicFEFunction(const Dune::array<TargetSpace,dim+1>& coefficients)
         : coefficients_(coefficients)
-    {
-        assert(coefficients_.size() == dim+1);
-    }
+    {}
 
     /** \brief Evaluate the function */
     TargetSpace evaluate(const Dune::FieldVector<ctype, dim>& local) const;
@@ -80,8 +78,8 @@ private:
         return B;
     }
         
-    static std::vector<ctype> barycentricCoordinates(const Dune::FieldVector<ctype,dim>& local) {
-        std::vector<ctype> result(dim+1);
+    static Dune::array<ctype,dim+1> barycentricCoordinates(const Dune::FieldVector<ctype,dim>& local) {
+        Dune::array<ctype,dim+1> result;
         result[0] = 1;
         for (int i=0; i<dim; i++) {
             result[0]  -= local[i];
@@ -144,7 +142,7 @@ private:
     }
     
     /** \brief The coefficient vector */
-    std::vector<TargetSpace> coefficients_;
+    Dune::array<TargetSpace,dim+1> coefficients_;
 
 };
 
@@ -153,11 +151,11 @@ TargetSpace LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
 evaluate(const Dune::FieldVector<ctype, dim>& local) const
 {
     // First compute the coordinates on the standard simplex (in R^{n+1})
-    std::vector<ctype> w = barycentricCoordinates(local);
+    Dune::array<ctype,dim+1> w = barycentricCoordinates(local);
 
-    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+    AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
 
-    TargetSpaceRiemannianTRSolver<TargetSpace> solver;
+    TargetSpaceRiemannianTRSolver<TargetSpace,dim+1> solver;
 
     solver.setup(&assembler,
                  coefficients_[0],   // initial iterate
@@ -209,8 +207,8 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
     Dune::FieldMatrix<ctype,embeddedDim,dim> RHS = dFdw * B;
 
     // the actual system matrix
-    std::vector<ctype> w = barycentricCoordinates(local);
-    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+    Dune::array<ctype,dim+1> w = barycentricCoordinates(local);
+    AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
     
     Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
     assembler.assembleEmbeddedHessian(q,dFdq);
@@ -302,7 +300,7 @@ evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& loc
 
     // dFdq
     std::vector<ctype> w = barycentricCoordinates(local);
-    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+    AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
     
     Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
     assembler.assembleEmbeddedHessian(q,dFdq);
@@ -405,7 +403,7 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
     
     // the actual system matrix
     std::vector<ctype> w = barycentricCoordinates(local);
-    AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+    AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
     
     Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
     assembler.assembleEmbeddedHessian(q,dFdq);
diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh
index 566630e61463582d1555bc50b4fe29e7006a77a7..27ec09ed12d96e4ff302119577c375197aca4dfe 100644
--- a/dune/gfe/localgeodesicfestiffness.hh
+++ b/dune/gfe/localgeodesicfestiffness.hh
@@ -40,17 +40,17 @@ public:
 
     */
     virtual void assembleHessian(const Entity& e,
-                  const std::vector<TargetSpace>& localSolution);
+                  const Dune::array<TargetSpace,gridDim+1>& localSolution);
 
     /** \brief Compute the energy at the current configuration */
     virtual RT energy (const Entity& e,
-                       const std::vector<TargetSpace>& localSolution) const = 0;
+                       const Dune::array<TargetSpace,gridDim+1>& localSolution) const = 0;
 
     /** \brief Assemble the element gradient of the energy functional 
 
     The default implementation in this class uses a finite difference approximation */
     virtual void assembleGradient(const Entity& element,
-                                  const std::vector<TargetSpace>& solution,
+                                  const Dune::array<TargetSpace,gridDim+1>& solution,
                                   std::vector<typename TargetSpace::TangentVector>& gradient) const;
     
     // assembled data
@@ -62,7 +62,7 @@ public:
 template <class GridView, class TargetSpace>
 void LocalGeodesicFEStiffness<GridView, TargetSpace>::
 assembleGradient(const Entity& element,
-                 const std::vector<TargetSpace>& localSolution,
+                 const Dune::array<TargetSpace,gridDim+1>& localSolution,
                  std::vector<typename TargetSpace::TangentVector>& localGradient) const
 {
 
@@ -74,8 +74,8 @@ assembleGradient(const Entity& element,
 
     localGradient.resize(localSolution.size());
 
-    std::vector<TargetSpace> forwardSolution = localSolution;
-    std::vector<TargetSpace> backwardSolution = localSolution;
+    Dune::array<TargetSpace,gridDim+1> forwardSolution = localSolution;
+    Dune::array<TargetSpace,gridDim+1> backwardSolution = localSolution;
 
     for (size_t i=0; i<localSolution.size(); i++) {
         
@@ -111,7 +111,7 @@ assembleGradient(const Entity& element,
 template <class GridType, class TargetSpace>
 void LocalGeodesicFEStiffness<GridType, TargetSpace>::
 assembleHessian(const Entity& element,
-         const std::vector<TargetSpace>& localSolution)
+         const Dune::array<TargetSpace,gridDim+1>& localSolution)
 {
     // 1 degree of freedom per element vertex
     size_t nDofs = element.template count<gridDim>();
@@ -144,8 +144,8 @@ assembleHessian(const Entity& element,
             Dune::FieldVector<double,embeddedBlocksize> minusEpsXi = epsXi;
             minusEpsXi  *= -1;
             
-            std::vector<TargetSpace> forwardSolution  = localSolution;
-            std::vector<TargetSpace> backwardSolution = localSolution;
+            Dune::array<TargetSpace,gridDim+1> forwardSolution  = localSolution;
+            Dune::array<TargetSpace,gridDim+1> backwardSolution = localSolution;
             
             forwardSolution[i]  = TargetSpace::exp(localSolution[i],epsXi);
             backwardSolution[i] = TargetSpace::exp(localSolution[i],minusEpsXi);
@@ -169,8 +169,8 @@ assembleHessian(const Entity& element,
                     Dune::FieldVector<double,embeddedBlocksize> minusEpsXi  = epsXi;   minusEpsXi  *= -1;
                     Dune::FieldVector<double,embeddedBlocksize> minusEpsEta = epsEta;  minusEpsEta *= -1;
                         
-                    std::vector<TargetSpace> forwardSolutionXiEta  = localSolution;
-                    std::vector<TargetSpace> backwardSolutionXiEta  = localSolution;
+                    Dune::array<TargetSpace,gridDim+1> forwardSolutionXiEta  = localSolution;
+                    Dune::array<TargetSpace,gridDim+1> backwardSolutionXiEta  = localSolution;
             
                     if (i==j)
                         forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi+epsEta);
diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index e05bcbf546fa7abb5d83c57dff61f901a381423c..4cc9caa8eb3313f912bab2799282e98a013980ca 100644
--- a/dune/gfe/targetspacertrsolver.cc
+++ b/dune/gfe/targetspacertrsolver.cc
@@ -6,9 +6,9 @@
 
 #include "maxnormtrustregion.hh"
 
-template <class TargetSpace>
-void TargetSpaceRiemannianTRSolver<TargetSpace>::
-setup(const AverageDistanceAssembler<TargetSpace>* assembler,
+template <class TargetSpace, int N>
+void TargetSpaceRiemannianTRSolver<TargetSpace,N>::
+setup(const AverageDistanceAssembler<TargetSpace,N>* assembler,
       const TargetSpace& x,
       double tolerance,
       int maxTrustRegionSteps,
@@ -45,8 +45,8 @@ setup(const AverageDistanceAssembler<TargetSpace>* assembler,
 }
 
 
-template <class TargetSpace>
-void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
+template <class TargetSpace, int N>
+void TargetSpaceRiemannianTRSolver<TargetSpace,N>::solve()
 {
     MaxNormTrustRegion<blocksize> trustRegion(1,   // we have only one block
                                               initialTrustRegionRadius_);
diff --git a/dune/gfe/targetspacertrsolver.hh b/dune/gfe/targetspacertrsolver.hh
index 3799700194f42df16163b7548aae4d9750630be2..55c3e48369c7f143b0a35558f1ea9a411010633d 100644
--- a/dune/gfe/targetspacertrsolver.hh
+++ b/dune/gfe/targetspacertrsolver.hh
@@ -8,8 +8,11 @@
 #include <dune/solvers/iterationsteps/trustregiongsstep.hh>
 #include <dune/solvers/norms/energynorm.hh>
 
-/** \brief Riemannian trust-region solver for geodesic finite-element problems */
-template <class TargetSpace>
+/** \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>
 class TargetSpaceRiemannianTRSolver 
     : public NumProc
 //     : public IterativeSolver<std::vector<TargetSpace>,
@@ -33,7 +36,7 @@ public:
     {}
 
     /** \brief Set up the solver using a monotone multigrid method as the inner solver */
-    void setup(const AverageDistanceAssembler<TargetSpace>* assembler,
+    void setup(const AverageDistanceAssembler<TargetSpace,N>* assembler,
                const TargetSpace& x,
                double tolerance,
                int maxTrustRegionSteps,
@@ -70,7 +73,7 @@ protected:
     double innerTolerance_;
 
     /** \brief The assembler for the average-distance functional */
-    const AverageDistanceAssembler<TargetSpace>* assembler_;
+    const AverageDistanceAssembler<TargetSpace,N>* assembler_;
 
     /** \brief The solver for the quadratic inner problems */
     std::auto_ptr< ::LoopSolver<CorrectionType> > innerSolver_;