Store normal stress in a vector

src/assemblers.cc

@@ -58,7 +58,7 @@ assemble_nonlinearity(
     Dune::ParameterTree const &parset,
     Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
     Dune::BlockVector<Dune::FieldVector<double, 1>> const &state,
-    double scalarNormalStress) {
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &normalStress) {
   auto const size = nodalIntegrals.size();
 
   typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
@@ -66,9 +66,6 @@ assemble_nonlinearity(
   SingletonVectorType mu(size);
   mu = parset.get<double>("mu");
 
-  SingletonVectorType normalStress(size);
-  normalStress = scalarNormalStress;
-
   switch (parset.get<Config::model>("model")) {
     case Config::Exponential: {
       SingletonVectorType a(size);

src/assemblers.hh

@@ -28,6 +28,5 @@ assemble_nonlinearity(
     Dune::ParameterTree const &parset,
     Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
     Dune::BlockVector<Dune::FieldVector<double, 1>> const &state,
-    double scalarNormalStress);
-
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &normalStress);
 #endif

src/assemblers_tmpl.cc

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

src/one-body-sample.cc

@@ -256,6 +256,12 @@ int main(int argc, char *argv[]) {
     SingletonVectorType alpha(alpha_old);
 
     SingletonVectorType vonMisesStress;
+
+    SingletonVectorType surfaceNormalStress(finestSize);
+    surfaceNormalStress = 0.0;
+    for (size_t i = 0; i < frictionalNodes.size(); ++i)
+      if (frictionalNodes[i][0])
+        surfaceNormalStress[i] = normalStress;
     // }}}
 
     typedef MyConvexProblem<MatrixType, VectorType> MyConvexProblemType;
@@ -345,7 +351,7 @@ int main(int argc, char *argv[]) {
           auto myGlobalNonlinearity =
               assemble_nonlinearity<MatrixType, VectorType>(
                   parset.sub("boundary.friction"), *nodalIntegrals, alpha,
-                  normalStress);
+                  surfaceNormalStress);
 
           MyConvexProblemType const myConvexProblem(
               problem_A, *myGlobalNonlinearity, problem_rhs);