Update the state

As a consequence, the nonlinearity is only assembled once

dune/tectonic/frictionpotential.hh

@@ -27,26 +27,15 @@ class FrictionPotentialWrapper {
   // Between 0 and this point, the function is constantly zero (and
   // thus so are all derivatives)
   double virtual smallestPositivePoint() const = 0;
+  void virtual updateState(double) = 0;
 };
 
 // phi(V) = V log(V/V_m) - V + V_m  if V >= V_m
 //        = 0                       otherwise
-// with V_m = V_0 exp(-K/a),
-// i.e.   K = -a log(V_m / V_0) = mu_0 + b log(V_0 / L) + b alpha
 class FrictionPotential : public FrictionPotentialWrapper {
 public:
-  FrictionPotential(double coefficient, FrictionData const &fd, double state)
-      : FrictionPotentialWrapper(),
-        coefficientProduct(coefficient * fd.a * fd.normalStress),
-        // state is assumed to be logarithmic
-        V_m(fd.V0 *
-            std::exp(-(fd.mu0 + fd.b * (state + std::log(fd.V0 / fd.L))) /
-                     fd.a))
-        // We could also compute V_m as
-        // V_0 * std::exp(-(mu_0 + b * state)/a)
-        //     * std::pow(V_0 / L, -b/a)
-        // which would avoid the std::exp(std::log())
-  {}
+  FrictionPotential(double coefficient, FrictionData const &fd)
+      : fd(fd), coefficientProduct(coefficient * fd.a * fd.normalStress) {}
 
   double evaluate(double V) const {
     assert(V >= 0);
@@ -89,7 +78,16 @@ class FrictionPotential : public FrictionPotentialWrapper {
 
   double smallestPositivePoint() const { return V_m; }
 
+  //    V_m = V_0 exp(-K/a),
+  // with K = -a log(V_m / V_0) = mu_0 + b [ alpha + log(V_0 / L) ]
+  void updateState(double state) {
+    // state is assumed to be logarithmic
+    V_m = fd.V0 *
+          std::exp(-(fd.mu0 + fd.b * (state + std::log(fd.V0 / fd.L))) / fd.a);
+  }
+
 private:
+  FrictionData const &fd;
   double const coefficientProduct;
   double V_m;
 };
@@ -107,6 +105,8 @@ class TrivialFunction : public FrictionPotentialWrapper {
   double smallestPositivePoint() const {
     return std::numeric_limits<double>::infinity();
   }
+
+  void updateState(double) {}
 };
 }
 #endif

dune/tectonic/globalnonlinearity.hh

@@ -14,6 +14,9 @@
 namespace Dune {
 template <class MatrixTypeTEMPLATE, class VectorTypeTEMPLATE>
 class GlobalNonlinearity {
+private:
+  using dataref = BlockVector<FieldVector<double, 1>> const &;
+
 public:
   using MatrixType = MatrixTypeTEMPLATE;
   using VectorType = VectorTypeTEMPLATE;
@@ -74,6 +77,8 @@ class GlobalNonlinearity {
     auto const res = restriction(i);
     return res->regularity(x);
   }
+
+  virtual void updateState(dataref state) = 0;
 };
 }
 #endif

dune/tectonic/globalruinanonlinearity.hh

@@ -24,20 +24,24 @@ class GlobalRuinaNonlinearity
   using GlobalNonlinearity<MatrixType, VectorType>::dim;
 
   GlobalRuinaNonlinearity(Dune::BitSetVector<1> const &frictionalNodes,
-                          dataref nodalIntegrals, FrictionData const &fd,
-                          dataref state)
+                          dataref nodalIntegrals, FrictionData const &fd)
       : restrictions(nodalIntegrals.size()) {
     auto trivialNonlinearity =
         make_shared<LocalFriction<dim>>(make_shared<TrivialFunction>());
     for (size_t i = 0; i < restrictions.size(); ++i) {
       restrictions[i] =
           frictionalNodes[i][0]
-              ? make_shared<LocalFriction<dim>>(make_shared<FrictionPotential>(
-                    nodalIntegrals[i], fd, state[i]))
+              ? make_shared<LocalFriction<dim>>(
+                    make_shared<FrictionPotential>(nodalIntegrals[i], fd))
               : trivialNonlinearity;
     }
   }
 
+  void updateState(dataref state) override {
+    for (size_t i = 0; i < restrictions.size(); ++i)
+      restrictions[i]->updateState(state[i]);
+  }
+
   /*
     Return a restriction of the outer function to the i'th node.
   */

dune/tectonic/localfriction.hh

@@ -25,6 +25,11 @@ template <int dimension> class LocalFriction {
     return func->evaluate(x.two_norm());
   }
 
+  void updateState(double state) {
+    func->updateState(state);
+    smp = func->smallestPositivePoint();
+  }
+
   double regularity(VectorType const &x) const {
     double const xnorm = x.two_norm();
     if (xnorm < 1e-14 && xnorm >= smp)

src/assemblers.cc

@@ -52,11 +52,10 @@ Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType>>
 assemble_nonlinearity(
     Dune::BitSetVector<1> const &frictionalNodes,
     Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    FrictionData const &fd,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state) {
+    FrictionData const &fd) {
   return Dune::make_shared<
       Dune::GlobalRuinaNonlinearity<MatrixType, VectorType>>(
-      frictionalNodes, nodalIntegrals, fd, state);
+      frictionalNodes, nodalIntegrals, fd);
 }
 
 #include "assemblers_tmpl.cc"

src/assemblers.hh

@@ -28,6 +28,5 @@ Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType>>
 assemble_nonlinearity(
     Dune::BitSetVector<1> const &frictionalNodes,
     Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    FrictionData const &fd,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state);
+    FrictionData const &fd);
 #endif

src/assemblers_tmpl.cc

@@ -35,5 +35,4 @@ template Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType>>
 assemble_nonlinearity<MatrixType, VectorType>(
     Dune::BitSetVector<1> const &frictionalNodes,
     Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    FrictionData const &fd,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state);
+    FrictionData const &fd);

src/one-body-sample.cc

@@ -321,6 +321,8 @@ int main(int argc, char *argv[]) {
 
     auto const nodalIntegrals = assemble_frictional<GridView, SmallVector>(
         leafView, p1Assembler, frictionalNodes);
+    auto myGlobalNonlinearity = assemble_nonlinearity<MatrixType, VectorType>(
+        frictionalNodes, *nodalIntegrals, frictionData);
 
     // {{{ Initial conditions
     VectorType u_initial(finestSize);
@@ -407,9 +409,7 @@ int main(int argc, char *argv[]) {
       size_t iterationCounter;
       auto solveVelocityProblem = [&](VectorType &_problem_iterate,
                                       SingletonVectorType const &_alpha) {
-        auto myGlobalNonlinearity =
-            assemble_nonlinearity<MatrixType, VectorType>(
-                frictionalNodes, *nodalIntegrals, frictionData, _alpha);
+        myGlobalNonlinearity->updateState(_alpha);
 
         using MyConvexProblemType = MyConvexProblem<MatrixType, VectorType>;
         MyConvexProblemType const myConvexProblem(