Integrate a state variable

We now re-assemble the global nonlinearity for every step

dune/tectonic/globalruinanonlinearity.hh

@@ -25,12 +25,14 @@ class GlobalRuinaNonlinearity
       shared_ptr<BlockVector<FieldVector<double, 1>> const> a,
       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> normalStress,
+      shared_ptr<BlockVector<FieldVector<double, 1>> const> state)
       : nodalIntegrals(nodalIntegrals),
         a(a),
         mu(mu),
         eta(eta),
         normalStress(normalStress),
+        state(state),
         trivialNonlinearity(
             new LocalNonlinearity<dim>(make_shared<TrivialFunction const>())),
         restrictions(nodalIntegrals->size()) // TODO: can we get the size from
@@ -52,7 +54,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]);
+        (*normalStress)[i][0], (*state)[i][0]);
     restrictions[i] = make_shared<LocalNonlinearity<dim> const>(func);
     return restrictions[i];
   }
@@ -65,6 +67,7 @@ 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;
 
   std::vector<shared_ptr<LocalNonlinearity<dim> const>> mutable restrictions;
 };

dune/tectonic/nicefunction.hh

@@ -31,13 +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 normalStress, double state)
       : coefficient(coefficient),
         a(a),
         mu(mu),
         rho(exp(-mu / a)),
         eta(eta),
-        normalStress(normalStress) {}
+        normalStress(normalStress),
+        state(state) {}
 
   /*
     If mu and sigma_n denote the coefficient of friction and the
@@ -55,8 +56,8 @@ class RuinaFunction : public NiceFunction {
   */
   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;
+    double const h = arg * (a * (std::log(arg) - 1) + mu + state);
+    y = 1 / eta * normalStress * h + a * rho + mu + state;
     y *= coefficient;
   }
 
@@ -71,7 +72,7 @@ class RuinaFunction : public NiceFunction {
     if (eta * s <= rho)
       return 0;
 
-    return coefficient * normalStress * (a * std::log(eta * s) + mu);
+    return coefficient * normalStress * (a * std::log(eta * s) + mu + state);
   }
 
   /* see above */
@@ -79,7 +80,7 @@ class RuinaFunction : public NiceFunction {
     if (eta * s <= rho)
       return 0;
 
-    return coefficient * normalStress * (a * std::log(eta * s) + mu);
+    return coefficient * normalStress * (a * std::log(eta * s) + mu + state);
   }
 
   /*
@@ -108,6 +109,7 @@ class RuinaFunction : public NiceFunction {
   double const mu;
   double const eta;
   double const normalStress;
+  double const state;
 
   double const rho;
 };

src/one-body-sample.cc

@@ -144,7 +144,8 @@ Dune::shared_ptr<Dune::GlobalNonlinearity<VectorType, MatrixType> const>
 assemble_nonlinearity(
     int size, Dune::ParameterTree const &parset,
     Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-        nodalIntegrals) {
+        nodalIntegrals,
+    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state) {
   typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
 
   // {{{ Assemble terms for the nonlinearity
@@ -165,7 +166,7 @@ assemble_nonlinearity(
 
     return Dune::make_shared<
         Dune::GlobalRuinaNonlinearity<VectorType, MatrixType> const>(
-        nodalIntegrals, a, mu, eta, normalStress);
+        nodalIntegrals, a, mu, eta, normalStress, state);
   } else if (friction_model == std::string("Laursen")) {
     return Dune::make_shared<Dune::GlobalLaursenNonlinearity<
         Dune::LinearFunction, VectorType, MatrixType> const>(mu, normalStress,
@@ -269,8 +270,10 @@ int main(int argc, char *argv[]) {
     auto nodalIntegrals =
         assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
             leafView, p1Basis, frictionalNodes);
-    auto myGlobalNonlinearity = assemble_nonlinearity<VectorType, OperatorType>(
-        grid->size(grid->maxLevel(), dim), parset, nodalIntegrals);
+    auto state =
+        Dune::make_shared<Dune::BlockVector<Dune::FieldVector<double, 1>>>(
+            grid->size(grid->maxLevel(), dim));
+    *state = 0.0;
 
     // {{{ Set up TNNMG solver
     // linear iteration step components
@@ -334,6 +337,10 @@ int main(int argc, char *argv[]) {
       stiffnessMatrix.mmv(u4, b4);
 
       if (parset.get<bool>("solver.nonlineargs.use")) {
+        auto myGlobalNonlinearity =
+            assemble_nonlinearity<VectorType, OperatorType>(
+                grid->size(grid->maxLevel(), dim), parset, nodalIntegrals,
+                state);
         MyConvexProblemType myConvexProblem(stiffnessMatrix,
                                             *myGlobalNonlinearity, b1);
         MyBlockProblemType myBlockProblem(parset, myConvexProblem);
@@ -349,6 +356,10 @@ int main(int argc, char *argv[]) {
       u1 += u1_diff;
 
       if (parset.get<bool>("solver.tnnmg.use")) {
+        auto myGlobalNonlinearity =
+            assemble_nonlinearity<VectorType, OperatorType>(
+                grid->size(grid->maxLevel(), dim), parset, nodalIntegrals,
+                state);
         MyConvexProblemType myConvexProblem(stiffnessMatrix,
                                             *myGlobalNonlinearity, b4);