diff --git a/dune/gfe/averageinterface.hh b/dune/gfe/averageinterface.hh
index 3adddd736ed2b5a6b8fe65dc37800adffd747a8d..af88f81a150dc14fc21db3bb4bd89508efbf32aa 100644
--- a/dune/gfe/averageinterface.hh
+++ b/dune/gfe/averageinterface.hh
@@ -10,7 +10,7 @@
 #include <dune/istl/preconditioners.hh>
 
 #include <dune/fufem/dgindexset.hh>
-#include <dune/fufem/crossproduct.hh>
+#include <dune/fufem/arithmetic.hh>
 #include <dune/fufem/surfmassmatrix.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 
@@ -500,7 +500,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, crossProduct(worldPos-center,value));
+            totalTorque.axpy(quad[ip].weight() * integrationElement, Arithmetic::crossProduct(worldPos-center,value));
 
         }
 
@@ -607,7 +607,7 @@ void computeAveragePressure(const typename RigidBodyMotion<double,GridView::dime
                         phi_i[j] = shapeFunctionValues[indexInFace];
 
                         mu_tilde[i][j].axpy(quad[qp].weight() * integrationElement,
-                                            crossProduct(Dune::FieldVector<double,dim>(worldPos-centerOfTorque), phi_i));
+                                            Arithmetic::crossProduct(Dune::FieldVector<double,dim>(worldPos-centerOfTorque), phi_i));
 
                     }
 
diff --git a/dune/gfe/rodfactory.hh b/dune/gfe/rodfactory.hh
index c48488c563003b32feee3a6853e30766dc15fe34..24ad9f86a8704cb69d3a65bfba1da85e14ff2ca9 100644
--- a/dune/gfe/rodfactory.hh
+++ b/dune/gfe/rodfactory.hh
@@ -3,7 +3,7 @@
 
 #include <vector>
 #include <dune/common/fvector.hh>
-#include <dune/fufem/crossproduct.hh>
+#include <dune/fufem/arithmetic.hh>
 
 #include <dune/gfe/rigidbodymotion.hh>
 #include <dune/gfe/localgeodesicfefunction.hh>
@@ -39,7 +39,7 @@ template <int dim>
 
     Dune::FieldVector<double,3> zAxis(0);
     zAxis[2] = 1;
-    Dune::FieldVector<double,3> axis = crossProduct(Dune::FieldVector<double,3>(end-beginning), zAxis);
+    Dune::FieldVector<double,3> axis = Arithmetic::crossProduct(Dune::FieldVector<double,3>(end-beginning), zAxis);
     if (axis.two_norm() != 0)
         axis /= -axis.two_norm();