Pass computed normalstress to nonlinearity asm

src/assemblers.cc

@@ -57,7 +57,8 @@ Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType> const>
 assemble_nonlinearity(
     Dune::ParameterTree const &parset,
     Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state) {
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state,
+    double scalarNormalStress) {
   auto const size = nodalIntegrals.size();
 
   typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
@@ -66,7 +67,7 @@ assemble_nonlinearity(
   mu = parset.get<double>("mu");
 
   SingletonVectorType normalStress(size);
-  normalStress = parset.get<double>("normalstress");
+  normalStress = scalarNormalStress;
 
   switch (parset.get<Config::model>("model")) {
     case Config::Exponential: {

src/assemblers.hh

@@ -27,6 +27,7 @@ Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType> const>
 assemble_nonlinearity(
     Dune::ParameterTree const &parset,
     Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state);
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state,
+    double scalarNormalStress);
 
 #endif

src/assemblers_tmpl.cc

@@ -32,4 +32,5 @@ template Dune::shared_ptr<
 assemble_nonlinearity<MatrixType, VectorType>(
     Dune::ParameterTree const &parset,
     Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state);
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state,
+    double scalarNormalStress);

src/one-body-sample.cc

@@ -167,6 +167,7 @@ int main(int argc, char *argv[]) {
     P0Basis const p0Basis(leafView);
     P1Basis const p1Basis(leafView);
 
+    double normalStress;
     MatrixType massMatrix;
     VectorType gravityFunctional;
     {
@@ -190,8 +191,7 @@ int main(int argc, char *argv[]) {
           area *= (upperRight[i] - lowerLeft[i]);
 
       double const gravity = 9.81;
-      double const normalStress =
-          parset.get<double>("boundary.friction.normalstress");
+      normalStress = parset.get<double>("boundary.friction.normalstress");
 
       // rho    = sigma     * A       / V   / g
       // kg/m^d = N/m^(d-1) * m^(d-1) / m^d / (N/kg)
@@ -344,7 +344,8 @@ int main(int argc, char *argv[]) {
         for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) {
           auto myGlobalNonlinearity =
               assemble_nonlinearity<MatrixType, VectorType>(
-                  parset.sub("boundary.friction"), *nodalIntegrals, alpha);
+                  parset.sub("boundary.friction"), *nodalIntegrals, alpha,
+                  normalStress);
 
           MyConvexProblemType const myConvexProblem(
               problem_A, *myGlobalNonlinearity, problem_rhs);