Allow the normal stress to be adjusted dynamically

src/one-body-sample.cc

@@ -51,6 +51,8 @@
 #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/assemblers/operatorassembler.hh>
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/computestress.hh>
 #include <dune/fufem/dunepython.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
@@ -241,6 +243,7 @@ int main(int argc, char *argv[]) {
     auto const nodalIntegrals =
         assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
             leafView, p1Basis, frictionalNodes);
+    BoundaryPatch<GridView> const frictionalBoundary(leafView, frictionalNodes);
 
     // {{{ Initialise vectors
     VectorType u(finestSize);
@@ -367,6 +370,21 @@ int main(int argc, char *argv[]) {
           timeSteppingScheme->extractDisplacement(u);
           timeSteppingScheme->extractVelocity(ud);
 
+          // Update the normal stress
+          if (parset.get<bool>("updateNormalStress")) {
+            VectorType tmp(finestSize);
+            Stress<GridView>::getElasticSurfaceNormalStress(frictionalBoundary,
+                                                            u, tmp, E, nu);
+            for (size_t i = 0; i < frictionalNodes.size(); ++i)
+              if (frictionalNodes[i][0]) {
+                if (tmp[i][1] < 0) {
+                  DUNE_THROW(Dune::Exception,
+                             "Encountered negative normal stress");
+                } else
+                  surfaceNormalStress[i] = tmp[i][1];
+              }
+          }
+
           // Update the state
           for (size_t i = 0; i < frictionalNodes.size(); ++i) {
             if (frictionalNodes[i][0]) {

src/one-body-sample.parset

@@ -10,6 +10,8 @@ enableTimer = false
 
 timeSteppingScheme = implicitEuler
 
+updateNormalStress = false
+
 [grid]
 refinements = 4