From 476b58b6889c45e604a1c6704aa38b8734ff5086 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sun, 5 Jan 2014 20:31:29 +0100
Subject: [PATCH] [Extend] Allow normal stress to vary in space
---
dune/tectonic/frictiondata.hh | 6 ++----
dune/tectonic/frictionpotential.hh | 12 +++++++-----
dune/tectonic/globalruinanonlinearity.hh | 6 ++++--
src/assemblers.cc | 5 +++--
src/assemblers.hh | 4 ++--
src/one-body-sample.cc | 13 +++++++------
6 files changed, 25 insertions(+), 21 deletions(-)
diff --git a/dune/tectonic/frictiondata.hh b/dune/tectonic/frictiondata.hh
index 8cc72469..0262f2df 100644
--- a/dune/tectonic/frictiondata.hh
+++ b/dune/tectonic/frictiondata.hh
@@ -3,19 +3,17 @@
#include <dune/common/parametertree.hh>
struct FrictionData {
- FrictionData(Dune::ParameterTree const &parset, double normalStress)
+ FrictionData(Dune::ParameterTree const &parset)
: L(parset.get<double>("L")),
V0(parset.get<double>("V0")),
a(parset.get<double>("a")),
b(parset.get<double>("b")),
- mu0(parset.get<double>("mu0")),
- normalStress(normalStress) {}
+ mu0(parset.get<double>("mu0")) {}
double const L;
double const V0;
double const a;
double const b;
double const mu0;
- double const normalStress;
};
#endif
diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh
index f0d6e64a..eb447bfb 100644
--- a/dune/tectonic/frictionpotential.hh
+++ b/dune/tectonic/frictionpotential.hh
@@ -28,15 +28,16 @@ class FrictionPotentialWrapper {
class FrictionPotential : public FrictionPotentialWrapper {
public:
- FrictionPotential(double coefficient, FrictionData const &_fd)
- : fd(_fd), weightTimesNormalStress(coefficient * (-fd.normalStress)) {}
+ FrictionPotential(double coefficient, double _normalStress,
+ FrictionData const &_fd)
+ : fd(_fd), weight(coefficient), normalStress(_normalStress) {}
double differential(double V) const {
assert(V >= 0.0);
if (V <= V_cutoff)
return 0.0;
- return weightTimesNormalStress * fd.a * (std::log(V) - logV_m);
+ return weight * (-normalStress) * fd.a * (std::log(V) - logV_m);
}
double second_deriv(double V) const {
@@ -44,7 +45,7 @@ class FrictionPotential : public FrictionPotentialWrapper {
if (V <= V_cutoff)
return 0;
- return weightTimesNormalStress * fd.a / V;
+ return weight * (-normalStress) * (fd.a / V);
}
double regularity(double V) const {
@@ -64,7 +65,8 @@ class FrictionPotential : public FrictionPotentialWrapper {
private:
FrictionData const &fd;
- double const weightTimesNormalStress;
+ double const weight;
+ double const normalStress;
double logV_m;
double V_cutoff; // velocity at which mu falls short of mumin
};
diff --git a/dune/tectonic/globalruinanonlinearity.hh b/dune/tectonic/globalruinanonlinearity.hh
index 0b365b1a..c95ed1f7 100644
--- a/dune/tectonic/globalruinanonlinearity.hh
+++ b/dune/tectonic/globalruinanonlinearity.hh
@@ -23,7 +23,8 @@ class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
public:
GlobalRuinaNonlinearity(Dune::BitSetVector<1> const &frictionalNodes,
- ScalarVector const &weights, FrictionData const &fd)
+ ScalarVector const &weights, FrictionData const &fd,
+ ScalarVector const &normalStress)
: restrictions(frictionalNodes.size()) {
auto trivialNonlinearity =
std::make_shared<Friction>(std::make_shared<TrivialFunction>());
@@ -32,7 +33,8 @@ class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
restrictions[i] = trivialNonlinearity;
continue;
}
- auto const fp = std::make_shared<FrictionPotential>(weights[i], fd);
+ auto const fp =
+ std::make_shared<FrictionPotential>(weights[i], normalStress[i], fd);
restrictions[i] = std::make_shared<Friction>(fp);
}
}
diff --git a/src/assemblers.cc b/src/assemblers.cc
index 83565cce..205ef090 100644
--- a/src/assemblers.cc
+++ b/src/assemblers.cc
@@ -92,7 +92,8 @@ void MyAssembler<GridView, dimension>::assembleNeumann(
template <class GridView, int dimension>
auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
- BoundaryPatch<GridView> const &frictionalBoundary, FrictionData const &fd)
+ BoundaryPatch<GridView> const &frictionalBoundary, FrictionData const &fd,
+ ScalarVector const &normalStress)
-> std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector>> {
// Lump negative normal stress (kludge)
ScalarVector weights;
@@ -110,7 +111,7 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
Dune::BitSetVector<1> frictionalNodes;
frictionalBoundary.getVertices(frictionalNodes);
return std::make_shared<GlobalRuinaNonlinearity<Matrix, Vector>>(
- frictionalNodes, weights, fd);
+ frictionalNodes, weights, fd, normalStress);
}
template <class GridView, int dimension>
diff --git a/src/assemblers.hh b/src/assemblers.hh
index a0068158..06723500 100644
--- a/src/assemblers.hh
+++ b/src/assemblers.hh
@@ -68,8 +68,8 @@ template <class GridView, int dimension> class MyAssembler {
std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector>>
assembleFrictionNonlinearity(
- BoundaryPatch<GridView> const &frictionalBoundary,
- FrictionData const &fd);
+ BoundaryPatch<GridView> const &frictionalBoundary, FrictionData const &fd,
+ ScalarVector const &normalStress);
void assembleVonMisesStress(double youngModulus, double poissonRatio,
Vector const &u, ScalarVector &stress);
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 1d2f8b27..651553c4 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -188,18 +188,19 @@ int main(int argc, char *argv[]) {
};
Vector ell(fineVertexCount);
computeExternalForces(0.0, ell);
- double normalStress = (myGeometry.A[1] - myGeometry.C[1]) *
- parset.get<double>("gravity") *
- parset.get<double>("body.density");
- FrictionData const frictionData(parset.sub("boundary.friction"),
- normalStress);
+
+ ScalarVector normalStress(fineVertexCount);
+ normalStress = (myGeometry.A[1] - myGeometry.C[1]) *
+ parset.get<double>("gravity") *
+ parset.get<double>("body.density");
+ FrictionData const frictionData(parset.sub("boundary.friction"));
// {{{ Initial conditions
ScalarVector alpha_initial(fineVertexCount);
alpha_initial =
std::log(parset.get<double>("boundary.friction.initialState"));
auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
- frictionalBoundary, frictionData);
+ frictionalBoundary, frictionData, normalStress);
myGlobalNonlinearity->updateLogState(alpha_initial);
using LinearFactory = SolverFactory<
--
GitLab