diff --git a/dune/gfe/averageinterface.hh b/dune/gfe/averageinterface.hh
index cdcd8e83af660086ccc9b462021df575ef7ef89d..c57c98a5ab8749cae92a86a018d8377a1eaa2696 100644
--- a/dune/gfe/averageinterface.hh
+++ b/dune/gfe/averageinterface.hh
@@ -9,6 +9,8 @@
 #include <dune/istl/solvers.hh>
 #include <dune/istl/preconditioners.hh>
 
+#include dune/matrix-vector/crossproduct.hh>
+
 #include <dune/fufem/dgindexset.hh>
 #include <dune/fufem/arithmetic.hh>
 #include <dune/fufem/surfmassmatrix.hh>
@@ -500,7 +502,7 @@ void computeTotalForceAndTorque(const BoundaryPatch<GridView>& interface,
             neumannFunction.evaluateLocal(*it->inside(), quadPos, value);
 
             totalForce.axpy(quad[ip].weight() * integrationElement, value);
-            totalTorque.axpy(quad[ip].weight() * integrationElement, Arithmetic::crossProduct(worldPos-center,value));
+            totalTorque.axpy(quad[ip].weight() * integrationElement, MatrixVector::crossProduct(worldPos-center,value));
 
         }
 
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index b995ccaa9c677e98a15d9cb6c704ecd50bff6e17..d2bdf188dd0af22222b5b030d28f61b49014a7d4 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -5,6 +5,8 @@
 #include <dune/common/parametertree.hh>
 #include <dune/geometry/quadraturerules.hh>
 
+#include <dune/matrix-vector/crossproduct.hh>
+
 #include <dune/fufem/functions/virtualgridfunction.hh>
 #include <dune/fufem/boundarypatch.hh>
 
@@ -351,7 +353,7 @@ energy(const typename Basis::LocalView& localView,
         aCovariant[i][j] = 0.0;
     }
 
-    aCovariant[2] = Arithmetic::crossProduct(aCovariant[0], aCovariant[1]);
+    aCovariant[2] = Dune::MatrixVector::crossProduct(aCovariant[0], aCovariant[1]);
     aCovariant[2] /= aCovariant[2].two_norm();
 
     auto aContravariant = aCovariant;
diff --git a/dune/gfe/rodfactory.hh b/dune/gfe/rodfactory.hh
index edfd47705fd4132f10c7fb19867fce4fdce975e8..94e06e451f012a8fe396f6124877f578a59d0624 100644
--- a/dune/gfe/rodfactory.hh
+++ b/dune/gfe/rodfactory.hh
@@ -5,6 +5,8 @@
 #include <dune/common/fvector.hh>
 #include <dune/fufem/arithmetic.hh>
 
+#include <dune/matrix-vector/crossproduct.hh>
+
 #include <dune/gfe/rigidbodymotion.hh>
 #include <dune/gfe/localgeodesicfefunction.hh>
 
@@ -39,7 +41,7 @@ template <int dim>
 
     Dune::FieldVector<double,3> zAxis(0);
     zAxis[2] = 1;
-    Dune::FieldVector<double,3> axis = Arithmetic::crossProduct(Dune::FieldVector<double,3>(end-beginning), zAxis);
+    Dune::FieldVector<double,3> axis = MatrixVector::crossProduct(Dune::FieldVector<double,3>(end-beginning), zAxis);
     if (axis.two_norm() != 0)
         axis /= -axis.two_norm();
 
diff --git a/dune/gfe/vertexnormals.hh b/dune/gfe/vertexnormals.hh
index fb899e0ee3af64a1e4ca608727db615999d0c113..49896216b810989f5c39ffb3b0b890d77300cc72 100644
--- a/dune/gfe/vertexnormals.hh
+++ b/dune/gfe/vertexnormals.hh
@@ -9,6 +9,8 @@
 
 #include <dune/geometry/referenceelements.hh>
 
+#include <dune/matrix-vector/crossproduct.hh>
+
 #include <dune/gfe/unitvector.hh>
 
 /** \brief Compute averaged vertex normals for a 2d-in-3d grid
@@ -29,7 +31,7 @@ std::vector<UnitVector<typename GridView::ctype,3> > computeVertexNormals(const
     {
       auto cornerPos = Dune::ReferenceElements<double,2>::general(element.type()).position(i,2);
       auto tangent = element.geometry().jacobianTransposed(cornerPos);
-      auto cornerNormal = Arithmetic::crossProduct(tangent[0], tangent[1]);
+      auto cornerNormal = Dune::MatrixVector::crossProduct(tangent[0], tangent[1]);
       cornerNormal /= cornerNormal.two_norm();
 
       unscaledNormals[indexSet.subIndex(element,i,2)] += cornerNormal;
@@ -50,4 +52,4 @@ std::vector<UnitVector<typename GridView::ctype,3> > computeVertexNormals(const
 }
 
 
-#endif
\ No newline at end of file
+#endif