[Problem] Compute normalStress

src/assemblers.cc

@@ -12,6 +12,7 @@
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh>
 #include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/computestress.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
 
@@ -90,6 +91,24 @@ void MyAssembler<GridView, dimension>::assembleNeumann(
                                              neumannBoundary);
 }
 
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleNormalStress(
+    BoundaryPatch<GridView> const &frictionalBoundary,
+    ScalarVector &normalStress, double youngModulus, double poissonRatio,
+    Vector const &displacement) {
+
+  Vector traction;
+  Stress<GridView>::getElasticSurfaceNormalStress // misnomer(!)
+      (frictionalBoundary, displacement, traction, youngModulus, poissonRatio);
+
+  std::vector<typename Vector::block_type> normals;
+  frictionalBoundary.getNormals(normals);
+  for (size_t i = 0; i < traction.size(); ++i) {
+    normalStress[i] = normals[i] * traction[i];
+    assert(normalStress[i] <= 0.0);
+  }
+}
+
 template <class GridView, int dimension>
 auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
     BoundaryPatch<GridView> const &frictionalBoundary,

src/assemblers.hh

@@ -67,6 +67,10 @@ template <class GridView, int dimension> class MyAssembler {
                        Dune::VirtualFunction<double, double> const &neumann,
                        double relativeTime);
 
+  void assembleNormalStress(BoundaryPatch<GridView> const &frictionalBoundary,
+                            ScalarVector &normalStress, double youngModulus,
+                            double poissonRatio, Vector const &displacement);
+
   std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector, GridView>>
   assembleFrictionNonlinearity(
       BoundaryPatch<GridView> const &frictionalBoundary,

src/one-body-sample.cc

@@ -190,21 +190,7 @@ int main(int argc, char *argv[]) {
     Vector ell(fineVertexCount);
     computeExternalForces(0.0, ell);
 
-    ScalarVector normalStress(fineVertexCount);
-    normalStress = (myGeometry.A[1] - myGeometry.C[1]) *
-                   parset.get<double>("gravity") *
-                   parset.get<double>("body.density");
-    MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"),
-                                            myGeometry);
-
     // {{{ Initial conditions
-    ScalarVector alpha_initial(fineVertexCount);
-    alpha_initial =
-        std::log(parset.get<double>("boundary.friction.initialState"));
-    auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
-        frictionalBoundary, frictionInfo, normalStress);
-    myGlobalNonlinearity->updateLogState(alpha_initial);
-
     using LinearFactory = SolverFactory<
         dims, BlockNonlinearTNNMGProblem<ConvexProblem<
                   ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
@@ -240,6 +226,20 @@ int main(int argc, char *argv[]) {
     u_initial = 0.0;
     solveLinearProblem(dirichletNodes, A, ell, u_initial, ANorm,
                        parset.sub("u0.solver"));
+
+    ScalarVector alpha_initial(fineVertexCount);
+    alpha_initial =
+        std::log(parset.get<double>("boundary.friction.initialState"));
+    ScalarVector normalStress(fineVertexCount);
+    myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
+                                     body.getYoungModulus(),
+                                     body.getPoissonRatio(), u_initial);
+    MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"),
+                                            myGeometry);
+    auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
+        frictionalBoundary, frictionInfo, normalStress);
+    myGlobalNonlinearity->updateLogState(alpha_initial);
+
     Vector v_initial(fineVertexCount);
     v_initial = 0.0;
     {