Add Ruina's approach

dune/tectonic/globalruinanonlinearity.hh

+/* -*- mode:c++; mode:semantic -*- */
+
+#ifndef DUNE_TECTONIC_GLOBAL_RUINA_NONLINEARITY_HH
+#define DUNE_TECTONIC_GLOBAL_RUINA_NONLINEARITY_HH
+
+#include <dune/common/fvector.hh>
+
+#include "localnonlinearity.hh"
+#include "globalnonlinearity.hh"
+
+namespace Dune {
+template <int dim>
+class GlobalRuinaNonlinearity
+    : public Dune::GlobalNonlinearity<dim, Dune::RuinaFunction> {
+public:
+  GlobalRuinaNonlinearity(
+      std::vector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
+      std::vector<double> const &a,
+      std::vector<double> const &coefficientOfFriction,
+      std::vector<double> const &eta, std::vector<double> const &normalStress)
+      : nodalIntegrals(nodalIntegrals),
+        a(a),
+        coefficientOfFriction(coefficientOfFriction),
+        eta(eta),
+        normalStress(normalStress) {}
+
+  /*
+    Return a restriction of the outer function to the i'th node.
+  */
+  RuinaFunction virtual restriction(int i) const {
+    return Dune::RuinaFunction(nodalIntegrals[i][0], a[i],
+                               coefficientOfFriction[i], eta[i],
+                               normalStress[i]);
+  }
+
+private:
+  std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
+  std::vector<double> a;
+  std::vector<double> coefficientOfFriction;
+  std::vector<double> eta;
+  std::vector<double> normalStress;
+};
+}
+#endif

dune/tectonic/nicefunction.hh

@@ -16,6 +16,64 @@ class NiceFunction : public VirtualFunction<double, double> {
   virtual double rightDifferential(double s) const = 0;
 };
 
+class RuinaFunction : public NiceFunction {
+public:
+  RuinaFunction() {}
+
+  RuinaFunction(double coefficient, double a, double mu, double eta,
+                double normalStress)
+      : coefficient(coefficient),
+        a(a),
+        mu(mu),
+        rho(exp(-mu / a)),
+        eta(eta),
+        normalStress(normalStress) {}
+
+  /*
+    If mu and sigma_n denote the coefficient of friction and the
+    normal stress, respectively, this function is given by
+
+    1/eta sigma_n h(max(eta id, rho)) + a rho + mu
+
+    with the constants a and b from Ruina's model and
+
+    h(beta) = beta (a (log beta - 1) + mu)
+
+    as well as
+
+    rho = exp(-mu/a)
+  */
+  void virtual evaluate(double const &x, double &y) const {
+    double const arg = std::max(eta * x, rho);
+    double const h = arg * (a * (std::log(arg) - 1) + mu);
+    y = 1 / eta * normalStress * h + a * rho + mu;
+    y *= coefficient;
+  }
+
+  double virtual leftDifferential(double s) const {
+    if (eta * s <= rho)
+      return 0;
+
+    return normalStress * (a * std::log(s) + mu);
+  }
+
+  double virtual rightDifferential(double s) const {
+    if (eta * s <= rho)
+      return 0;
+
+    return normalStress * (a * std::log(s) + mu);
+  }
+
+private:
+  double coefficient;
+  double a;
+  double mu;
+  double eta;
+  double normalStress;
+
+  double rho;
+};
+
 class LinearFunction : public NiceFunction {
 public:
   LinearFunction() {}

src/one-body-sample.cc

@@ -42,6 +42,7 @@
 #include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
 
 #include <dune/tectonic/globallaursennonlinearity.hh>
+#include <dune/tectonic/globalruinanonlinearity.hh>
 #include <dune/tectonic/myconvexproblem.hh>
 #include <dune/tectonic/myblockproblem.hh>
 
@@ -186,6 +187,8 @@ int main(int argc, char *argv[]) {
 
     typedef MyConvexProblem<OperatorType, VectorType, Dune::LinearFunction>
     MyConvexProblemType;
+    // typedef MyConvexProblem<OperatorType, VectorType, Dune::RuinaFunction>
+    // MyConvexProblemType;
     typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType;
     GenericNonlinearGS<MyBlockProblemType> nonlinearGSStep;
     nonlinearGSStep.ignoreNodes_ = &ignoreNodes;
@@ -216,18 +219,29 @@ int main(int argc, char *argv[]) {
       stiffnessMatrix.umv(u3, b3);
 
       // {{{ Assemble terms for the nonlinearity
-      std::vector<double> normalStress;
-      normalStress.resize(grid.size(grid.maxLevel(), dim));
-      std::fill(normalStress.begin(), normalStress.end(),
-                parset.get<double>("boundary.normalstress"));
+      std::vector<double> a;
+      a.resize(grid.size(grid.maxLevel(), dim));
+      std::fill(a.begin(), a.end(), parset.get<double>("boundary.a"));
 
       std::vector<double> coefficientOfFriction;
       coefficientOfFriction.resize(grid.size(grid.maxLevel(), dim));
       std::fill(coefficientOfFriction.begin(), coefficientOfFriction.end(),
                 parset.get<double>("boundary.mu"));
 
+      std::vector<double> eta;
+      eta.resize(grid.size(grid.maxLevel(), dim));
+      std::fill(eta.begin(), eta.end(), parset.get<double>("boundary.eta"));
+
+      std::vector<double> normalStress;
+      normalStress.resize(grid.size(grid.maxLevel(), dim));
+      std::fill(normalStress.begin(), normalStress.end(),
+                parset.get<double>("boundary.normalstress"));
+
       Dune::GlobalLaursenNonlinearity<dim, Dune::LinearFunction>
       myGlobalNonlinearity(coefficientOfFriction, normalStress, nodalIntegrals);
+      // Dune::GlobalRuinaNonlinearity<dim>
+      //   myGlobalNonlinearity(nodalIntegrals, a, coefficientOfFriction, eta,
+      // normalStress);
       // }}}
 
       {

src/one-body-sample.parset

@@ -12,5 +12,7 @@ maxiterations = 100000
 tolerance = 1e-6
 
 [boundary]
+a = 0.0015
 normalstress = 0.1
 mu = 0.75
+eta = 1
\ No newline at end of file