Introduce L and b

dune/tectonic/globalruinanonlinearity.hh

@@ -26,13 +26,17 @@ class GlobalRuinaNonlinearity
       shared_ptr<BlockVector<FieldVector<double, 1>> const> mu,
       shared_ptr<BlockVector<FieldVector<double, 1>> const> eta,
       shared_ptr<BlockVector<FieldVector<double, 1>> const> normalStress,
-      shared_ptr<BlockVector<FieldVector<double, 1>> const> state, double h)
+      shared_ptr<BlockVector<FieldVector<double, 1>> const> b,
+      shared_ptr<BlockVector<FieldVector<double, 1>> const> state,
+      shared_ptr<BlockVector<FieldVector<double, 1>> const> L, double h)
       : nodalIntegrals(nodalIntegrals),
         a(a),
         mu(mu),
         eta(eta),
         normalStress(normalStress),
+        b(b),
         state(state),
+        L(L),
         h(h),
         trivialNonlinearity(
             new LocalNonlinearity<dim>(make_shared<TrivialFunction const>())),
@@ -55,7 +59,7 @@ class GlobalRuinaNonlinearity
 
     auto const func = make_shared<RuinaFunction const>(
         (*nodalIntegrals)[i][0], (*a)[i][0], (*mu)[i][0], (*eta)[i][0],
-        (*normalStress)[i][0], (*state)[i][0], h);
+        (*normalStress)[i][0], (*b)[i][0], (*state)[i][0], (*L)[i][0], h);
     restrictions[i] = make_shared<LocalNonlinearity<dim> const>(func);
     return restrictions[i];
   }
@@ -68,7 +72,9 @@ class GlobalRuinaNonlinearity
   shared_ptr<BlockVector<FieldVector<double, 1>> const> mu;
   shared_ptr<BlockVector<FieldVector<double, 1>> const> eta;
   shared_ptr<BlockVector<FieldVector<double, 1>> const> normalStress;
+  shared_ptr<BlockVector<FieldVector<double, 1>> const> b;
   shared_ptr<BlockVector<FieldVector<double, 1>> const> state;
+  shared_ptr<BlockVector<FieldVector<double, 1>> const> L;
   double const h;
 
   std::vector<shared_ptr<LocalNonlinearity<dim> const>> mutable restrictions;

dune/tectonic/nicefunction.hh

@@ -31,14 +31,14 @@ class NiceFunction : public VirtualFunction<double, double> {
 class RuinaFunction : public NiceFunction {
 public:
   RuinaFunction(double coefficient, double a, double mu, double eta,
-                double normalStress, double state, double h)
+                double normalStress, double b, double state, double L, double h)
       : coefficient(coefficient),
         a(a),
         mu(mu),
         rho(exp(-mu / a)),
         eta(eta),
         normalStress(normalStress),
-        state(state),
+        compound_state(b * (state - log(eta * L))),
         h(h) {}
 
   /*
@@ -57,8 +57,8 @@ class RuinaFunction : public NiceFunction {
   */
   void virtual evaluate(double const &x, double &y) const {
     double const arg = std::max(eta * x / h, rho);
-    double const r = arg * (a * (std::log(arg) - 1) + mu + state);
-    y = 1 / eta * normalStress * r + a * rho + mu + state;
+    double const r = arg * (a * (std::log(arg) - 1) + mu + compound_state);
+    y = 1 / eta * normalStress * r + a * rho + mu + compound_state;
     y *= coefficient;
     y *= h;
   }
@@ -75,7 +75,8 @@ class RuinaFunction : public NiceFunction {
     if (arg <= rho)
       return 0;
 
-    return coefficient * normalStress * (a * std::log(arg) + mu + state);
+    return coefficient * normalStress *
+           (a * std::log(arg) + mu + compound_state);
   }
 
   /* see above */
@@ -84,7 +85,8 @@ class RuinaFunction : public NiceFunction {
     if (arg <= rho)
       return 0;
 
-    return coefficient * normalStress * (a * std::log(arg) + mu + state);
+    return coefficient * normalStress *
+           (a * std::log(arg) + mu + compound_state);
   }
 
   /*
@@ -113,7 +115,7 @@ class RuinaFunction : public NiceFunction {
   double const mu;
   double const eta;
   double const normalStress;
-  double const state;
+  double const compound_state;
   double const h;
 
   double const rho;

src/one-body-sample.cc

@@ -168,9 +168,15 @@ assemble_nonlinearity(
     auto eta = Dune::make_shared<SingletonVectorType>(size);
     *eta = parset.get<double>("boundary.friction.eta");
 
+    auto b = Dune::make_shared<SingletonVectorType>(size);
+    *b = 1; // FIXME
+
+    auto L = Dune::make_shared<SingletonVectorType>(size);
+    *L = 1.0 / parset.get<double>("boundary.friction.eta"); // FIXME
+
     return Dune::make_shared<
         Dune::GlobalRuinaNonlinearity<VectorType, MatrixType> const>(
-        nodalIntegrals, a, mu, eta, normalStress, state, h);
+        nodalIntegrals, a, mu, eta, normalStress, b, state, L, h);
   } else if (friction_model == std::string("Laursen")) {
     return
         // TODO: take state and h into account