diff --git a/dune/tectonic/mydirectionalconvexfunction.hh b/dune/tectonic/mydirectionalconvexfunction.hh
index 0b2016d4230b7ffe34d62090c9079c8db304118b..9f147ae43b0e87f2b11af9cecc4576473530ef91 100644
--- a/dune/tectonic/mydirectionalconvexfunction.hh
+++ b/dune/tectonic/mydirectionalconvexfunction.hh
@@ -55,9 +55,9 @@ template <class NonlinearityType> class MyDirectionalConvexFunction {
   double linearPart() const { return b; }
 
   void subDiff(double x, Interval<double> &D) const {
-    VectorType tmp = u;
-    tmp.axpy(x, v);
-    phi.directionalSubDiff(tmp, v, D);
+    VectorType uxv = u;
+    Arithmetic::addProduct(uxv, x, v);
+    phi.directionalSubDiff(uxv, v, D);
     D[0] += A * x - b;
     D[1] += A * x - b;
   }
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index a7b0719d70eb24f3c6a94ac867650b02603358b3..0c01069ee8e2cdfa4afce641b63bc9800441edf0 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -447,7 +447,7 @@ int main(int argc, char *argv[]) {
             dampingWriter << "N ";
           } else { // alpha is still the old time step here
             alpha *= damping;
-            alpha.axpy(1.0 - damping, computed_state);
+            Arithmetic::addProduct(alpha, 1.0 - damping, computed_state);
             dampingWriter << "Y ";
           }
           lastCorrection = correction;