diff --git a/dune/gfe/tensor3.hh b/dune/gfe/tensor3.hh
index 10144b7926e7bde06945864f74c68a4440c0ed72..fa18c10c04aa5106ae68681615a0a86b60d6108a 100644
--- a/dune/gfe/tensor3.hh
+++ b/dune/gfe/tensor3.hh
@@ -35,6 +35,13 @@ class Tensor3
             return norm;
         }
         
+        Tensor3<T,N1,N2,N3>& axpy(const T& alpha, const Tensor3<T,N1,N2,N3>& other)
+        {
+            for (int i=0; i<N1; i++)
+                (*this)[i].axpy(alpha,other[i]);
+            return *this;
+        }
+        
         static Tensor3<T,N1,N2,N3> product(const Dune::FieldVector<T,N1>& a, const Dune::FieldVector<T,N2>& b, const Dune::FieldVector<T,N3>& c)
         {
             Tensor3<T,N1,N2,N3> result;
@@ -87,6 +94,22 @@ class Tensor3
             return result;
         }
 
+        template <int N4>
+        friend Tensor3<T,N1,N3,N4> operator*(const Dune::FieldMatrix<T,N1,N2>& a, const Tensor3<T,N2,N3,N4>& b)
+        {
+            Tensor3<T,N1,N3,N4> result;
+            
+            for (int i=0; i<N1; i++)
+                for (int j=0; j<N3; j++)
+                    for (int k=0; k<N4; k++) {
+                        result[i][j][k] = 0;
+                        for (int l=0; l<N2; l++)
+                            result[i][j][k] += a[i][l]*b[l][j][k];
+                    }
+                    
+            return result;
+        }
+
     friend Tensor3<T,N1,N2,N3> operator+(const Tensor3<T,N1,N2,N3>& a, const Tensor3<T,N1,N2,N3>& b)
     {
             Tensor3<T,N1,N2,N3> result;