diff --git a/src/Makefile.am b/src/Makefile.am
index 170446c7d675b9659e1484f5f83f2fc7d54a8d66..46010a7ace993dce98a84add41d30bdab34ccb06 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -6,7 +6,9 @@ bin_PROGRAMS = \
SOURCES = \
assemblers.cc \
- compute_state_dieterich.cc \
+ compute_state_dieterich_newmark.cc \
+ compute_state_dieterich_euler.cc \
+ compute_state_dieterich_common.cc \
compute_state_ruina.cc \
mysolver.cc \
one-body-sample.cc \
diff --git a/src/compute_state_dieterich.cc b/src/compute_state_dieterich.cc
deleted file mode 100644
index 3b8441cfe58d489177e436f6bed89803ab0a0a77..0000000000000000000000000000000000000000
--- a/src/compute_state_dieterich.cc
+++ /dev/null
@@ -1,63 +0,0 @@
-#include <cassert>
-
-#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/fufem/interval.hh>
-
-#include <dune/tnnmg/problem-classes/bisection.hh>
-
-#include <dune/tectonic/mydirectionalconvexfunction.hh>
-
-#include "compute_state_dieterich.hh"
-
-// psi(beta) = V/L beta + e^(-beta)
-class DieterichNonlinearity {
-public:
- typedef Dune::FieldVector<double, 1> VectorType;
- typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
-
- DieterichNonlinearity(double _VoL) : VoL(_VoL) {}
-
- void directionalSubDiff(VectorType const &u, VectorType const &v,
- Interval<double> &D) const {
- D[0] = D[1] = v[0] * (VoL - std::exp(-u[0]));
- }
-
- void directionalDomain(VectorType const &, VectorType const &,
- Interval<double> &dom) const {
- dom[0] = -std::numeric_limits<double>::max();
- dom[1] = std::numeric_limits<double>::max();
- }
-
-private:
- double const VoL;
-};
-
-double state_update_dieterich_bisection(double tau, double VoL,
- double old_state) {
- DieterichNonlinearity::VectorType const start(0);
- DieterichNonlinearity::VectorType const direction(1);
- DieterichNonlinearity const phi(VoL);
- MyDirectionalConvexFunction<DieterichNonlinearity> const J(
- 1.0 / tau, old_state / tau, 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(double tau, double VoL, double old_state) {
- double const ret = state_update_dieterich_bisection(tau, VoL, old_state);
- /* We have
-
- ret - old_state + tau * VoL = tau*exp(-ret)
-
- or
-
- log((ret - old_state)/tau + VoL) = -ret
- */
- assert(std::min(std::abs(ret - old_state + tau * (VoL - std::exp(-ret))),
- std::abs(std::log((ret - old_state) / tau + VoL) + ret)) <
- 1e-12);
- return ret;
-}
diff --git a/src/compute_state_dieterich.hh b/src/compute_state_dieterich.hh
deleted file mode 100644
index 49fe3dcc3d7bc112b31d0d8d2193a12811b5121a..0000000000000000000000000000000000000000
--- a/src/compute_state_dieterich.hh
+++ /dev/null
@@ -1,5 +0,0 @@
-#ifndef COMPUTE_STATE_HH
-#define COMPUTE_STATE_HH
-
-double state_update_dieterich(double h, double VoL, double old_state);
-#endif
diff --git a/src/compute_state_dieterich_common.cc b/src/compute_state_dieterich_common.cc
new file mode 100644
index 0000000000000000000000000000000000000000..449bcbe356cfda2636cd21b24d7e0ffaa0470ad0
--- /dev/null
+++ b/src/compute_state_dieterich_common.cc
@@ -0,0 +1,16 @@
+#include "compute_state_dieterich_common.hh"
+
+DieterichNonlinearity::DieterichNonlinearity(double _VoL) : VoL(_VoL) {}
+
+void DieterichNonlinearity::directionalSubDiff(VectorType const &u,
+ VectorType const &v,
+ Interval<double> &D) const {
+ D[0] = D[1] = v[0] * (VoL - std::exp(-u[0]));
+}
+
+void DieterichNonlinearity::directionalDomain(VectorType const &,
+ VectorType const &,
+ Interval<double> &dom) const {
+ dom[0] = -std::numeric_limits<double>::max();
+ dom[1] = std::numeric_limits<double>::max();
+}
diff --git a/src/compute_state_dieterich_common.hh b/src/compute_state_dieterich_common.hh
new file mode 100644
index 0000000000000000000000000000000000000000..39f5f62745a7efb022868737a1bef52f4177239c
--- /dev/null
+++ b/src/compute_state_dieterich_common.hh
@@ -0,0 +1,18 @@
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/fufem/interval.hh>
+
+class DieterichNonlinearity {
+public:
+ typedef Dune::FieldVector<double, 1> VectorType;
+ typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
+
+ DieterichNonlinearity(double _VoL);
+ void directionalSubDiff(VectorType const &u, VectorType const &v,
+ Interval<double> &D) const;
+ void directionalDomain(VectorType const &, VectorType const &,
+ Interval<double> &dom) const;
+
+private:
+ double const VoL;
+};
diff --git a/src/compute_state_dieterich_euler.cc b/src/compute_state_dieterich_euler.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ca51487a4e892f8c9f8d15ed4fca21de844eb2ec
--- /dev/null
+++ b/src/compute_state_dieterich_euler.cc
@@ -0,0 +1,38 @@
+#include <cassert>
+
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+#include <dune/tectonic/mydirectionalconvexfunction.hh>
+
+#include "compute_state_dieterich_common.hh"
+#include "compute_state_dieterich_euler.hh"
+
+double state_update_dieterich_euler_bisection(double tau, double VoL,
+ double old_state) {
+ DieterichNonlinearity::VectorType const start(0);
+ DieterichNonlinearity::VectorType const direction(1);
+ DieterichNonlinearity const phi(VoL);
+ MyDirectionalConvexFunction<DieterichNonlinearity> const J(
+ 1.0 / tau, old_state / tau, 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_euler(double tau, double VoL, double old_state) {
+ double const ret =
+ state_update_dieterich_euler_bisection(tau, VoL, old_state);
+ /* We have
+
+ ret - old_state + tau * VoL = tau*exp(-ret)
+
+ or
+
+ log((ret - old_state)/tau + VoL) = -ret
+ */
+ assert(std::min(std::abs(ret - old_state + tau * (VoL - std::exp(-ret))),
+ std::abs(std::log((ret - old_state) / tau + VoL) + ret)) <
+ 1e-12);
+ return ret;
+}
diff --git a/src/compute_state_dieterich_euler.hh b/src/compute_state_dieterich_euler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5c6c5f9162cb9413b0b38dc1a5e13ade0c4fb6e7
--- /dev/null
+++ b/src/compute_state_dieterich_euler.hh
@@ -0,0 +1,5 @@
+#ifndef COMPUTE_STATE_DIETERICH_EULER_HH
+#define COMPUTE_STATE_DIETERICH_EULER_HH
+
+double state_update_dieterich_euler(double h, double VoL, double old_state);
+#endif
diff --git a/src/compute_state_dieterich_newmark.cc b/src/compute_state_dieterich_newmark.cc
new file mode 100644
index 0000000000000000000000000000000000000000..18704f488bd3eb96e2466857a763c9227d59cd20
--- /dev/null
+++ b/src/compute_state_dieterich_newmark.cc
@@ -0,0 +1,35 @@
+#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;
+}
diff --git a/src/compute_state_dieterich_newmark.hh b/src/compute_state_dieterich_newmark.hh
new file mode 100644
index 0000000000000000000000000000000000000000..755dbf03fbef3ca2de041e441c2351ea680d66d7
--- /dev/null
+++ b/src/compute_state_dieterich_newmark.hh
@@ -0,0 +1,6 @@
+#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
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index b4f376aec869b57f7719caea47f167cb69556d4b..8551d62ed936e46d6e9fa7439bea7365cae5db37 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -65,7 +65,8 @@
#include <dune/tectonic/myconvexproblem.hh>
#include "assemblers.hh"
-#include "compute_state_dieterich.hh"
+#include "compute_state_dieterich_euler.hh"
+#include "compute_state_dieterich_newmark.hh"
#include "compute_state_ruina.hh"
#include "mysolver.hh"
#include "vtk.hh"
@@ -255,6 +256,10 @@ 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;
// }}}
@@ -375,7 +380,21 @@ int main(int argc, char *argv[]) {
switch (parset.get<Config::state_model>(
"boundary.friction.state.model")) {
case Config::Dieterich:
- alpha[i] = state_update_dieterich(tau, V / L, alpha_old[i]);
+ 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
+ }
break;
case Config::Ruina:
alpha[i] = state_update_ruina(tau, V * tau / L, alpha_old[i]);
@@ -430,6 +449,7 @@ 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)
diff --git a/src/one-body-sample.parset b/src/one-body-sample.parset
index 64fa3a2b1df348f56e9be3ca418651030c437619..a6e3d7e76b5b4ba7c3a604d60a53530d12dc125f 100644
--- a/src/one-body-sample.parset
+++ b/src/one-body-sample.parset
@@ -11,6 +11,7 @@ printVelocitySteppingComparison = false
enableTimer = false
timeSteppingScheme = implicitEuler
+stateTimeSteppingScheme = implicitEuler
[grid]
refinements = 4