Remove Newmark step in the time stepping schem

This partially reverts d84ddb712891e2fe8a4bd067589a31c86e8a3eea

src/Makefile.am

@@ -6,7 +6,6 @@ bin_PROGRAMS = \
 
 SOURCES = \
 	assemblers.cc \
-	compute_state_dieterich_newmark.cc \
 	compute_state_dieterich_euler.cc \
 	compute_state_dieterich_common.cc \
 	compute_state_ruina.cc \

src/compute_state_dieterich_newmark.cc

-#include <cassert>
-
-#include <dune/tnnmg/problem-classes/bisection.hh>
-
-#include <dune/tectonic/mydirectionalconvexfunction.hh>
-
-#include "compute_state_dieterich_common.hh"
-#include "compute_state_dieterich_newmark.hh"
-
-double state_update_dieterich_newmark_bisection(double tau, double VoL,
-                                                double old_state,
-                                                double old_stated) {
-  DieterichNonlinearity::VectorType const start(0);
-  DieterichNonlinearity::VectorType const direction(1);
-  DieterichNonlinearity const phi(VoL);
-  MyDirectionalConvexFunction<DieterichNonlinearity> const J(
-      2.0 / tau, old_stated + 2.0 / tau * old_state, phi, start, direction);
-
-  int bisectionsteps = 0;
-  Bisection const bisection(0.0, 1.0, 1e-12, true, 0);    // TODO
-  return bisection.minimize(J, 0.0, 0.0, bisectionsteps); // TODO
-}
-
-double state_update_dieterich_newmark(double tau, double VoL, double old_state,
-                                      double old_stated) {
-  double const new_state =
-      state_update_dieterich_newmark_bisection(tau, VoL, old_state, old_stated);
-  double const new_stated = 2.0 / tau * (new_state - old_state) - old_stated;
-  double const rhs = -new_stated;
-  double const lhs = VoL - std::exp(-new_state);
-
-  assert(std::abs(lhs - rhs) < 1e-10);
-
-  return new_state;
-}

src/compute_state_dieterich_newmark.hh

-#ifndef COMPUTE_STATE_DIETERICH_NEWMARK_HH
-#define COMPUTE_STATE_DIETERICH_NEWMARK_HH
-
-double state_update_dieterich_newmark(double h, double VoL, double old_state,
-                                      double old_stated);
-#endif

src/compute_state_dieterich_newmark.tex

-\documentclass{scrartcl}
-\usepackage{amsmath}
-\begin{document}
-\noindent The Newmark scheme in its classical form with $\gamma = 1/2$
-and $\beta = 1/4$ reads
-\begin{align}
-  \label{eq:1} \dot \alpha_1
-  &= \dot \alpha_0 + \frac \tau 2 (\ddot \alpha_0 + \ddot \alpha_1 )\\
-  \label{eq:2} \alpha_1
-  &= \alpha_0 + \tau \dot \alpha_0 + \frac {\tau^2}4 ( \ddot \alpha_0 + \ddot \alpha_1 )\text.
-  \intertext{We can also write \eqref{eq:2} as}
-  \nonumber \ddot \alpha_1 + \ddot \alpha_0
-  &= \frac 4{\tau^2} ( \alpha_1 - \alpha_0 - \tau \dot \alpha_0)
-  \intertext{so that it yields}
-  \label{eq:3} \dot \alpha_1
-  &= \dot \alpha_0 + \frac 2\tau ( \alpha_1 - \alpha_0 - \tau \dot \alpha_0) = \frac 2\tau ( \alpha_1 - \alpha_0) - \dot \alpha_0
-\end{align}
-in conjunction with \eqref{eq:1}. The problem
-\begin{align*}
-  -\dot \alpha_1 \in \partial \psi(\alpha_1)
-\end{align*}
-then becomes
-\begin{align*}
-  \dot \alpha_0 -\frac 2\tau ( \alpha_1 - \alpha_0)
-  &\in \partial \psi(\alpha_1)\\
-  \psi(\beta) - \psi(\alpha_1)
-  &\ge (\dot \alpha_0 -\frac 2\tau ( \alpha_1 - \alpha_0), \beta - \alpha_1)
-  \quad \forall \beta\\
-  \frac 2\tau ( \alpha_1, \beta - \alpha_1) + \psi(\beta) - \psi(\alpha_1)
-  &\ge (\dot \alpha_0 + \frac 2\tau \alpha_0, \beta - \alpha_1)
-  \quad \forall \beta
-\end{align*}
-After which $\dot \alpha_1$ can be computed according to \eqref{eq:3}.
-\end{document}

src/one-body-sample.cc

@@ -66,7 +66,6 @@
 
 #include "assemblers.hh"
 #include "compute_state_dieterich_euler.hh"
-#include "compute_state_dieterich_newmark.hh"
 #include "compute_state_ruina.hh"
 #include "mysolver.hh"
 #include "vtk.hh"
@@ -256,10 +255,6 @@ int main(int argc, char *argv[]) {
     alpha_old = parset.get<double>("boundary.friction.state.initial");
     SingletonVectorType alpha(alpha_old);
 
-    SingletonVectorType alphad_old(finestSize);
-    alphad_old = 0.0;
-    SingletonVectorType alphad(alphad_old);
-
     SingletonVectorType vonMisesStress;
     // }}}
 
@@ -380,21 +375,8 @@ int main(int argc, char *argv[]) {
               switch (parset.get<Config::state_model>(
                   "boundary.friction.state.model")) {
                 case Config::Dieterich:
-                  switch (
-                      parset.get<Config::scheme>("stateTimeSteppingScheme")) {
-                    case Config::ImplicitEuler:
-                      alpha[i] = state_update_dieterich_euler(tau, V / L,
-                                                              alpha_old[i]);
-                      break;
-                    case Config::Newmark:
-                      alpha[i] = state_update_dieterich_newmark(
-                          tau, V / L, alpha_old[i], alphad_old[i]);
-                      alphad[i] =
-                          2.0 / tau * (alpha[i] - alpha_old[i]) - alphad_old[i];
-                      break;
-                    default:
-                      assert(false); // Nothing else implemented
-                  }
+                  alpha[i] =
+                      state_update_dieterich_euler(tau, V / L, alpha_old[i]);
                   break;
                 case Config::Ruina:
                   alpha[i] = state_update_ruina(tau, V * tau / L, alpha_old[i]);
@@ -449,7 +431,6 @@ int main(int argc, char *argv[]) {
       }
 
       alpha_old = alpha;
-      alphad_old = alphad;
 
       { // Write all kinds of data
         for (size_t i = 0; i < frictionalNodes.size(); ++i)

src/one-body-sample.parset

@@ -11,7 +11,6 @@ printVelocitySteppingComparison = false
 enableTimer = false
 
 timeSteppingScheme = implicitEuler
-stateTimeSteppingScheme = implicitEuler
 
 [grid]
 refinements = 4