Implement evolution of state

src/one-body-sample.cc

@@ -15,6 +15,8 @@
 
 #include "boost/format.hpp"
 
+#include <cmath>
+
 #include "LambertW.h"
 
 #include <dune/common/bitsetvector.hh>
@@ -58,6 +60,7 @@
 #include <dune/tectonic/globalruinanonlinearity.hh>
 #include <dune/tectonic/myconvexproblem.hh>
 #include <dune/tectonic/myblockproblem.hh>
+#include <dune/tectonic/mydirectionalconvexfunction.hh>
 #include <dune/grid/io/file/vtk/vtkwriter.hh>
 
 int const dim = 3;
@@ -180,6 +183,41 @@ assemble_nonlinearity(
   }
 }
 
+// The nonlinearity exp(-x)
+class DecayingExponential {
+public:
+  typedef Dune::FieldVector<double, 1> VectorType;
+  typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
+
+  double operator()(VectorType const &u) const { return exp(-u[0]); }
+
+  void directionalSubDiff(VectorType const &u, VectorType const &v,
+                          Interval<double> &D) const {
+    D[0] = D[1] = v[0] * (-exp(-u[0]));
+  }
+
+  void directionalDomain(VectorType const &, VectorType const &,
+                         Interval<double> &dom) const {
+    dom[0] = -std::numeric_limits<double>::max();
+    dom[1] = std::numeric_limits<double>::max();
+  }
+};
+
+double compute_state_update_bisection(double h, double unorm, double L,
+                                      double old_state) {
+  MyDirectionalConvexFunction<DecayingExponential> const Jplus(
+      1.0 / h, (old_state - unorm / L) / h, DecayingExponential(), 0, 1);
+  int bisectionsteps = 0;
+  Bisection const bisection(0.0, 1.0, 1e-12, true, 0); // TODO
+  return bisection.minimize(Jplus, 0.0, 0.0, bisectionsteps);
+}
+
+double compute_state_update_lambert(double h, double unorm, double L,
+                                    double old_state) {
+  double const rhs = unorm / L - old_state;
+  return LambertW(0, h * exp(rhs)) - rhs;
+}
+
 int main(int argc, char *argv[]) {
   try {
     Dune::ParameterTree parset;
@@ -258,12 +296,21 @@ int main(int argc, char *argv[]) {
     VectorType u3 = u1;
     VectorType u4 = u1;
 
+    Dune::BlockVector<Dune::FieldVector<double, 1>> s4_old(
+        grid->size(grid->maxLevel(), dim));
+    s4_old = 50; // FIXME: magic value (-500 is still workable; -1000 is not)
+
     VectorType u1_diff(grid->size(grid->maxLevel(), dim));
     u1_diff = 0.0; // Has to be zero!
     VectorType u2_diff = u1_diff;
     VectorType u3_diff = u1_diff;
     VectorType u4_diff = u1_diff;
 
+    auto s4_new =
+        Dune::make_shared<Dune::BlockVector<Dune::FieldVector<double, 1>>>(
+            grid->size(grid->maxLevel(), dim));
+    *s4_new = s4_old;
+
     CellVectorType vonMisesStress;
 
     VectorType b1;
@@ -360,27 +407,47 @@ int main(int argc, char *argv[]) {
       u1 += u1_diff;
 
       if (parset.get<bool>("solver.tnnmg.use")) {
-        auto state =
-            Dune::make_shared<Dune::BlockVector<Dune::FieldVector<double, 1>>>(
-                grid->size(grid->maxLevel(), dim));
-        *state = 0.0;
-        auto myGlobalNonlinearity =
-            assemble_nonlinearity<VectorType, OperatorType>(
-                grid->size(grid->maxLevel(), dim), parset, nodalIntegrals,
-                state, h);
-        MyConvexProblemType myConvexProblem(stiffnessMatrix,
-                                            *myGlobalNonlinearity, b4);
-
-        MyBlockProblemType myBlockProblem(parset, myConvexProblem);
-        multigridStep.setProblem(u4_diff, myBlockProblem);
-
-        LoopSolver<VectorType> overallSolver(
-            &multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
-            solver_tolerance, &energyNorm, verbosity);
-        overallSolver.solve();
+        for (int state_fpi = 0; state_fpi < parset.get<int>("state.iterations");
+             ++state_fpi) {
+          auto myGlobalNonlinearity =
+              assemble_nonlinearity<VectorType, OperatorType>(
+                  grid->size(grid->maxLevel(), dim), parset, nodalIntegrals,
+                  s4_new, h);
+          MyConvexProblemType myConvexProblem(stiffnessMatrix,
+                                              *myGlobalNonlinearity, b4);
+
+          MyBlockProblemType myBlockProblem(parset, myConvexProblem);
+          multigridStep.setProblem(u4_diff, myBlockProblem);
+
+          LoopSolver<VectorType> overallSolver(
+              &multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
+              solver_tolerance, &energyNorm, verbosity);
+          overallSolver.solve();
+
+          if (parset.get<bool>("state.enable")) {
+            for (size_t i = 0; i < frictionalNodes.size(); ++i)
+              if (frictionalNodes[i][0]) {
+                double const L = 1e-4; // FIXME: magic value
+                double const unorm = u4_diff[i].two_norm();
+
+                double ret1 =
+                    compute_state_update_bisection(h, unorm, L, s4_old[i][0]);
+                assert(abs(1.0 / h * ret1 - (s4_old[i] - unorm / L) / h -
+                           exp(-ret1)) < 1e-12);
+                double ret2 =
+                    compute_state_update_lambert(h, unorm, L, s4_old[i][0]);
+                assert(abs(1.0 / h * ret2 - (s4_old[i] - unorm / L) / h -
+                           exp(-ret2)) < 1e-12);
+                assert(abs(ret2 - ret1) < 1e-14);
+
+                (*s4_new)[i][0] = ret1;
+              }
+          }
+        }
       }
 
       u4 += u4_diff;
+      s4_old = *s4_new;
 
       // Compute von Mises stress and write everything to a file
       if (parset.get<bool>("writeVTK")) {

src/one-body-sample.parset

@@ -8,6 +8,10 @@ printDifference = false
 
 writeVTK = false
 
+[state]
+enable = true
+iterations = 6
+
 [grid]
 refinements = 3