From 7c7282b649dac901b59e6887b98da12270f5716f Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Thu, 19 Dec 2013 03:11:56 +0100
Subject: [PATCH] [Problem] Compute normalStress
---
src/assemblers.cc | 19 +++++++++++++++++++
src/assemblers.hh | 4 ++++
src/one-body-sample.cc | 28 ++++++++++++++--------------
3 files changed, 37 insertions(+), 14 deletions(-)
diff --git a/src/assemblers.cc b/src/assemblers.cc
index b9a3e76b..d78bb999 100644
--- a/src/assemblers.cc
+++ b/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,
diff --git a/src/assemblers.hh b/src/assemblers.hh
index a281251e..c90235cc 100644
--- a/src/assemblers.hh
+++ b/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,
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 918f4994..eb4b38ac 100644
--- a/src/one-body-sample.cc
+++ b/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;
{
--
GitLab