diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index e13acbd8d4934b9b0fd6010175be22eb0a519e71..eb4cab23ad962cd83c373964c46eadc0178b79cd 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -151,9 +151,9 @@ private:
     }
 
     /** \brief Compute derivate of F(w,q) (the derivative of the weighted distance fctl) wrt to w */
-    Dune::Matrix<Dune::FieldMatrix<RT,1,1> > computeDFdw(const TargetSpace& q) const
+    Dune::Matrix<RT> computeDFdw(const TargetSpace& q) const
     {
-        Dune::Matrix<Dune::FieldMatrix<RT,1,1> > dFdw(embeddedDim,localFiniteElement_.localBasis().size());
+        Dune::Matrix<RT> dFdw(embeddedDim,localFiniteElement_.localBasis().size());
         for (size_t i=0; i<localFiniteElement_.localBasis().size(); i++) {
             Dune::FieldVector<RT,embeddedDim> tmp = TargetSpace::derivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], q);
             for (int j=0; j<embeddedDim; j++)
@@ -250,20 +250,20 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local, const TargetSpace
     localFiniteElement_.localBasis().evaluateJacobian(local, B);
 
     // compute negative derivative of F(w,q) (the derivative of the weighted distance fctl) wrt to w
-    Dune::Matrix<Dune::FieldMatrix<RT,1,1> > dFdw = computeDFdw(q);
+    Dune::Matrix<RT> dFdw = computeDFdw(q);
     dFdw *= -1;
 
     // multiply the two previous matrices: the result is the right hand side
     // RHS = dFdw * B;
     // We need to write this multiplication by hand, because B is actually an (#coefficients) times 1 times dim matrix,
     // and we need to ignore the middle index.
-    Dune::Matrix<Dune::FieldMatrix<RT,1,1> > RHS(dFdw.N(), dim);
+    Dune::Matrix<RT> RHS(dFdw.N(), dim);
 
     for (size_t i=0; i<RHS.N(); i++)
       for (size_t j=0; j<RHS.M(); j++) {
         RHS[i][j] = 0;
         for (size_t k=0; k<dFdw.M(); k++)
-          RHS[i][j] += dFdw[i][k][0][0]*B[k][0][j];
+          RHS[i][j] += dFdw[i][k]*B[k][0][j];
       }
 
     // the actual system matrix
@@ -413,7 +413,7 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
     // the matrix that turns coordinates on the reference simplex into coordinates on the standard simplex
     std::vector<Dune::FieldMatrix<ctype,1,dim> > BNested(coefficients_.size());
     localFiniteElement_.localBasis().evaluateJacobian(local, BNested);
-    Dune::Matrix<Dune::FieldMatrix<RT,1,1> > B(coefficients_.size(), dim);
+    Dune::Matrix<RT> B(coefficients_.size(), dim);
     for (size_t i=0; i<coefficients_.size(); i++)
         for (size_t j=0; j<dim; j++)
             B[i][j] = BNested[i][0][j];
diff --git a/dune/gfe/tensorssd.hh b/dune/gfe/tensorssd.hh
index ac83b94ce4d57e7721b1e6bed113f62101043eca..4d1a256ab0502682a2c33bd0068bf02131a14230 100644
--- a/dune/gfe/tensorssd.hh
+++ b/dune/gfe/tensorssd.hh
@@ -71,7 +71,7 @@ public:
         return *this;
     }
 
-    friend TensorSSD<T,N1,N2> operator*(const TensorSSD<T,N1,N2>& a, const Dune::Matrix<Dune::FieldMatrix<T,1,1> >& b)
+    friend TensorSSD<T,N1,N2> operator*(const TensorSSD<T,N1,N2>& a, const Dune::Matrix<T>& b)
     {
         TensorSSD<T,N1,N2> result(b.M());
             
@@ -83,7 +83,7 @@ public:
                 for (size_t k=0; k<b.M(); k++) {
                     result.data_[i][j][k] = 0;
                     for (size_t l=0; l<N4; l++)
-                        result.data_[i][j][k] += a.data_[i][j][l]*b[l][k][0][0];
+                        result.data_[i][j][k] += a.data_[i][j][l]*b[l][k];
                 }
                     
         return result;