[Cleanup] Rearrange Newmark algorithm

src/timestepping/newmark.cc

@@ -27,20 +27,47 @@ void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
 
   tau = _tau;
 
-  problem_rhs = ell;
-  Arithmetic::addProduct(problem_rhs, 2.0 / tau, B, v_o);
-  Arithmetic::addProduct(problem_rhs, B, a_o);
-  Arithmetic::subtractProduct(problem_rhs, A, u_o);
-  Arithmetic::subtractProduct(problem_rhs, tau / 2.0, A, v_o);
-
-  // For fixed tau, we'd only really have to do this once
-  Dune::MatrixIndexSet indices(A.N(), A.M());
-  indices.import(A);
-  indices.import(B);
-  indices.exportIdx(problem_AB);
+  /* We start out with the formulation
+
+       B a + A u = ell
+
+     Newmark means
+
+       a1 = 2/tau ( v1 - v0 ) - a0
+       u1 = tau/2 ( v1 + v0 ) + u0
+
+     in summary, we get at time t=1
+
+       B [2/tau ( u1 - u0 ) - a0]
+       + A [tau/2 ( v1 + v0 ) + u0] = ell
+
+     or
+
+       2/tau B v1 + tau/2 A v1
+       = [2/tau B + tau/2 A] v1
+       = ell + 2/tau B v0 + B a0
+       - tau/2 A v0 - A u0
+  */
+
+  // set up LHS (for fixed tau, we'd only really have to do this once)
+  {
+    Dune::MatrixIndexSet indices(A.N(), A.M());
+    indices.import(A);
+    indices.import(B);
+    indices.exportIdx(problem_AB);
+  }
   problem_AB = 0.0;
-  Arithmetic::addProduct(problem_AB, tau / 2.0, A);
   Arithmetic::addProduct(problem_AB, 2.0 / tau, B);
+  Arithmetic::addProduct(problem_AB, tau / 2.0, A);
+
+  // set up RHS
+  {
+    problem_rhs = ell;
+    Arithmetic::addProduct(problem_rhs, 2.0 / tau, B, v_o);
+    Arithmetic::addProduct(problem_rhs, B, a_o);
+    Arithmetic::subtractProduct(problem_rhs, tau / 2.0, A, v_o);
+    Arithmetic::subtractProduct(problem_rhs, A, u_o);
+  }
 
   // v_o makes a good initial iterate; we could use anything, though
   problem_iterate = 0.0;