Add indices before computing A+B

src/timestepping.cc

@@ -2,6 +2,7 @@
 #include "config.h"
 #endif
 
+#include <dune/istl/matrixindexset.hh>
 #include <dune/fufem/arithmetic.hh>
 #include "timestepping.hh"
 
@@ -262,8 +263,12 @@ void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
   A.usmv(-tau / 2.0, ud_old, problem_rhs);
 
   // For fixed tau, we'd only really have to do this once
-  problem_A = A;
-  problem_A *= tau / 2.0;
+  Dune::MatrixIndexSet indices(A.N(), A.M());
+  indices.import(A);
+  indices.import(B);
+  indices.exportIdx(problem_A);
+  problem_A = 0.0;
+  Arithmetic::addProduct(problem_A, tau / 2.0, A);
   Arithmetic::addProduct(problem_A, 2.0 / tau, B);
 
   // ud_old makes a good initial iterate; we could use anything, though