diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 094b58ff7bfa2bdef151c40f0a25a0844fc21cee..cc2d2ad6a8a46932cd334f0faedd36f31572a4da 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -433,7 +433,7 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
     dvDwF = 0;
     for (int i=0; i<embeddedDim; i++)
         for (int j=0; j<embeddedDim; j++)
-            dvDwF.data_[i][j][coefficient] = mixedDerivative[i][j];
+            dvDwF(i, j, coefficient) = mixedDerivative[i][j];
     
     
     // dFDq is not invertible, if the target space is embedded into a higher-dimensional
@@ -462,16 +462,16 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
         Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> hesse = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[k], q);
         for (int i=0; i<embeddedDim; i++)
             for (int j=0; j<embeddedDim; j++)
-                dqdwF.data_[i][j][k] = hesse[i][j];
+                dqdwF(i, j, k) = hesse[i][j];
     }
 
     TensorSSD<double, embeddedDim,embeddedDim> dqdwF_times_dvq(coefficients_.size());
     for (int i=0; i<embeddedDim; i++)
         for (int j=0; j<embeddedDim; j++)
             for (size_t k=0; k<coefficients_.size(); k++) {
-                dqdwF_times_dvq.data_[i][j][k] = 0;
+                dqdwF_times_dvq(i, j, k) = 0;
                 for (int l=0; l<embeddedDim; l++)
-                    dqdwF_times_dvq.data_[i][j][k] += dqdwF.data_[l][j][k] * dvq[l][i];
+                    dqdwF_times_dvq(i, j, k) += dqdwF(l, j, k) * dvq[l][i];
             }
 
     Tensor3<double, embeddedDim,embeddedDim,dim> foo;
@@ -482,7 +482,7 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
     for (int i=0; i<embeddedDim; i++)
         for (int j=0; j<embeddedDim; j++)
             for (int k=0; k<dim; k++)
-                foo[i][j][k] -= bar.data_[i][j][k];
+                foo[i][j][k] -= bar(i, j, k);
     
     result = 0;
     for (int i=0; i<embeddedDim; i++)
diff --git a/dune/gfe/tensorssd.hh b/dune/gfe/tensorssd.hh
index 17f724482e53096ce541c4a5bac6112945c7c6a6..876b9b09c1fe5587af4585e36d4644a1ca79c76e 100644
--- a/dune/gfe/tensorssd.hh
+++ b/dune/gfe/tensorssd.hh
@@ -41,6 +41,18 @@ public:
                 assert(false);
         }
     }
+
+    /** \brief Direct access to individual entries */
+    T& operator()(size_t i, size_t j, size_t k)
+    {
+        return data_[i][j][k];
+    }
+        
+    /** \brief Direct const access to individual entries */
+    const T& operator()(size_t i, size_t j, size_t k) const
+    {
+        return data_[i][j][k];
+    }
         
     /** \brief Assignment from scalar */
     TensorSSD<T,N1,N2>& operator=(const T& scalar)
@@ -84,6 +96,8 @@ public:
                     
         return result;
     }
+    
+private:
 
     // having the dynamic data type on the inside is kind of a stupid data layout
     Dune::array<Dune::array<std::vector<T>, N2>, N1> data_;