Add Ruina state evolution law

src/Makefile.am

@@ -25,6 +25,7 @@ run-one-body-sample-gdb: one-body-sample
 one_body_sample_SOURCES = \
 	assemblers.cc \
 	compute_state.cc \
+	compute_state_ruina.cc \
 	mysolver.cc \
 	one-body-sample.cc \
 	print.cc \

src/compute_state_ruina.cc

+#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_ruina.hh"
+
+// Ruina's nonlinearity
+class Ruina {
+public:
+  typedef Dune::FieldVector<double, 1> VectorType;
+  typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
+
+  Ruina(double h, double uol) : h(h), uol(uol) {}
+
+  void directionalSubDiff(VectorType const &u, VectorType const &v,
+                          Interval<double> &D) const {
+    if (uol == 0)
+      D[0] = D[1] = v[0] * uol * (u[0]);
+    else
+      D[0] = D[1] = v[0] * uol * (u[0] + std::log(uol / h));
+  }
+
+  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 h;
+  double const uol;
+};
+
+double state_update_ruina_bisection(double h, double uol, double old_state) {
+  Ruina::VectorType const start(0);
+  Ruina::VectorType const direction(1);
+  Ruina const phi(h, uol);
+  MyDirectionalConvexFunction<Ruina> const J(1.0, 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_ruina(double h, double uol, // unorm/L
+                          double old_state) {
+  double const uol_log_uol = (uol == 0) ? 0.0 : uol * std::log(uol);
+
+  double ret1 = state_update_ruina_bisection(h, uol, old_state);
+  assert(std::abs(ret1 + uol * ret1 - old_state + uol_log_uol -
+                  uol * std::log(h)) < 1e-14);
+
+  return ret1;
+}

src/compute_state_ruina.hh

+#ifndef COMPUTE_STATE_RUINA_HH
+#define COMPUTE_STATE_RUINA_HH
+
+double state_update_ruina(double h, double uol, double old_state);
+
+#endif

src/one-body-sample.cc

@@ -63,6 +63,7 @@
 
 #include "assemblers.hh"
 #include "compute_state.hh"
+#include "compute_state_ruina.hh"
 #include "mysolver.hh"
 #include "print.hh"
 #include "vtk.hh"
@@ -277,7 +278,14 @@ int main(int argc, char *argv[]) {
               // velocity
               // std::cout << std::log(L/unorm * h) << std::endl;
 
-              (*s4_new)[i] = state_update_dieterich(h, unorm / L, s4_old[i]);
+              auto const model =
+                  parset.get<std::string>("boundary.friction.state.model");
+              if (model == std::string("Dieterich"))
+                (*s4_new)[i] = state_update_dieterich(h, unorm / L, s4_old[i]);
+              else if (model == std::string("Ruina"))
+                (*s4_new)[i] = state_update_ruina(h, unorm / L, s4_old[i]);
+              else
+                assert(false);
             }
           }
           if (parset.get<bool>("printProgress")) {

src/one-body-sample.parset

@@ -73,6 +73,7 @@ model = Ruina
 evolve = true
 # log(alpha(0))
 initial = -2.42037 # Dirichlet: -4.72295 (1e-4), -2.42037 (1e-3), -0.117783 (1e-2)
+model = Dieterich # Ruina
 
 [boundary.friction.ruina]
 # "For rocks, typical values of A and B range from 0.005 to 0.015"