From d0e40d822d52443b71420d41e46dcadddb3a836f Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Tue, 16 Jul 2013 16:56:55 +0200
Subject: [PATCH] [Cleanup] Rearrange Newmark algorithm
---
src/timestepping/newmark.cc | 51 ++++++++++++++++++++++++++++---------
1 file changed, 39 insertions(+), 12 deletions(-)
diff --git a/src/timestepping/newmark.cc b/src/timestepping/newmark.cc
index a08dad6b..7ccb02ad 100644
--- a/src/timestepping/newmark.cc
+++ b/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;
--
GitLab