diff --git a/dune/tectonic/ellipticenergy.hh b/dune/tectonic/ellipticenergy.hh
index c63c7a76eee0ba1a9b683c5284f48786c0e18b22..31551ad5fe646c3b6ebf1c0e83758ac69585e391 100644
--- a/dune/tectonic/ellipticenergy.hh
+++ b/dune/tectonic/ellipticenergy.hh
@@ -6,6 +6,7 @@
 #include <dune/common/stdstreams.hh>
 
 #include <dune/fufem/interval.hh>
+#include <dune/fufem/arithmetic.hh>
 
 #include "localnonlinearity.hh"
 
@@ -25,7 +26,7 @@ template <int dim> class EllipticEnergy {
 
   double operator()(SmallVector const &v) const {
     SmallVector y(0);
-    A.usmv(0.5, v, y);        //  1/2 Av
+    Arithmetic::addProduct(y, 0.5, A, v);
     y -= b;                   //  1/2 Av - b
     return y * v + (*phi)(v); // <1/2 Av - b,v> + H(|v|)
   }
diff --git a/src/timestepping.cc b/src/timestepping.cc
index fd15efdd2e482c145d39e21b84fa3455ffdb3b08..1922e0f2dd5a77c59169b1342b212e0f00fe53f9 100644
--- a/src/timestepping.cc
+++ b/src/timestepping.cc
@@ -129,10 +129,10 @@ void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
   tau = _tau;
 
   problem_rhs = ell;
-  B.usmv(2.0 / tau, ud_old, problem_rhs);
-  B.usmv(1.0, udd_old, problem_rhs);
-  A.usmv(-1.0, u_old, problem_rhs);
-  A.usmv(-tau / 2.0, ud_old, problem_rhs);
+  Arithmetic::addProduct(problem_rhs, 2.0 / tau, B, ud_old);
+  Arithmetic::addProduct(problem_rhs, B, udd_old);
+  Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
+  Arithmetic::addProduct(problem_rhs, -tau / 2.0, A, ud_old);
 
   // For fixed tau, we'd only really have to do this once
   Dune::MatrixIndexSet indices(A.N(), A.M());
@@ -237,8 +237,8 @@ void EulerPair<VectorType, MatrixType, FunctionType, dim>::setup(
   tau = _tau;
 
   problem_rhs = ell;
-  B.usmv(1.0 / tau, ud_old, problem_rhs);
-  A.usmv(-1.0, u_old, problem_rhs);
+  Arithmetic::addProduct(problem_rhs, 1.0 / tau, B, ud_old);
+  Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
 
   // For fixed tau, we'd only really have to do this once
   Dune::MatrixIndexSet indices(A.N(), A.M());