diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index e7b1074c99d1510fdc7df80c0ad7b0c0d67a2ccc..c33bb2c94465d92c6492d164f7d121f81c51ae3e 100644 --- a/src/one-body-sample.cc +++ b/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]) { diff --git a/src/one-body-sample.parset b/src/one-body-sample.parset index e24da9d6aed27791bedb36d0de25ac8f7791ba2a..d19c9673cf0a51683cfbafefeef7b117c686014c 100644 --- a/src/one-body-sample.parset +++ b/src/one-body-sample.parset @@ -10,6 +10,8 @@ enableTimer = false timeSteppingScheme = implicitEuler +updateNormalStress = false + [grid] refinements = 4