diff --git a/dune/gfe/harmonicenergystiffness.hh b/dune/gfe/harmonicenergystiffness.hh
index 3fb31da5cab1f4d50b8d939de1d8970b34b842ff..66cf878f0b008083ff881ffc33a3eb16b0b2e214 100644
--- a/dune/gfe/harmonicenergystiffness.hh
+++ b/dune/gfe/harmonicenergystiffness.hh
@@ -146,14 +146,15 @@ assembleEmbeddedGradient(const Entity& element,
         // loop over all the element's degrees of freedom and compute the gradient wrt it
         for (size_t i=0; i<localSolution.size(); i++) {
          
-            Tensor3<double, TargetSpace::EmbeddedTangentVector::size, gridDim,TargetSpace::EmbeddedTangentVector::size> derivativeDerivative;
+            Tensor3<double, TargetSpace::EmbeddedTangentVector::size,TargetSpace::EmbeddedTangentVector::size,gridDim> derivativeDerivative;
             localGeodesicFEFunction.evaluateDerivativeOfGradientWRTCoefficient(quadPos, i, derivativeDerivative);
         
             for (int j=0; j<derivative.rows; j++) {
                 
                 for (int k=0; k<derivative.cols; k++) {
                     
-                    localGradient[i].axpy(weight*derivative[j][k], derivativeDerivative[j][k]);
+                    for (int l=0; l<TargetSpace::EmbeddedTangentVector::size; l++)
+                        localGradient[i][l] += weight*derivative[j][k] * derivativeDerivative[l][j][k];
                     
                 }
                 
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index bb2098ef903f6078e7ac2cef83da3d6d1d7bd26f..2045b0def46220fc234d86e0a5e8f53226367444 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -88,7 +88,7 @@ public:
     /** \brief Evaluate the derivative of the gradient of the function with respect to a coefficient */
     void evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
                                                     int coefficient,
-                                                    Tensor3<double, TargetSpace::EmbeddedTangentVector::size,dim,TargetSpace::EmbeddedTangentVector::size>& result) const;
+                                                    Tensor3<double, TargetSpace::EmbeddedTangentVector::size,TargetSpace::EmbeddedTangentVector::size,dim>& result) const;
 
 private:
 
@@ -280,7 +280,7 @@ template <int dim, class ctype, class TargetSpace>
 void LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
 evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
                                            int coefficient,
-                                           Tensor3<double, TargetSpace::EmbeddedTangentVector::size,dim,TargetSpace::EmbeddedTangentVector::size>& result) const
+                                           Tensor3<double, TargetSpace::EmbeddedTangentVector::size,TargetSpace::EmbeddedTangentVector::size,dim>& result) const
 {
     const int embeddedDim = EmbeddedTangentVector::size;