From bc9e179dbc72d8dcbd18976fa6410b4597072dd5 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Thu, 19 Jun 2014 17:00:54 +0200
Subject: [PATCH] [Extend] Add regularised rate/state law
---
dune/tectonic/frictionpotential.hh | 32 ++++++++++++++++++++++++++++++
src/assemblers.cc | 16 ++++++++++++---
src/assemblers.hh | 3 +++
src/enum_friction_model.cc | 15 ++++++++++++++
src/enums.hh | 4 ++++
src/sand-wedge-data/parset.cfg | 1 +
src/sand-wedge.cc | 2 ++
7 files changed, 70 insertions(+), 3 deletions(-)
create mode 100644 src/enum_friction_model.cc
diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh
index e3f1d27c..3ac28d2a 100644
--- a/dune/tectonic/frictionpotential.hh
+++ b/dune/tectonic/frictionpotential.hh
@@ -69,6 +69,38 @@ class TruncatedRateState : public FrictionPotential {
double Vmin;
};
+class RegularisedRateState : public FrictionPotential {
+public:
+ RegularisedRateState(double coefficient, double _normalStress,
+ FrictionData _fd)
+ : fd(_fd), weight(coefficient), normalStress(_normalStress) {}
+
+ double coefficientOfFriction(double V) const {
+ return fd.a * std::asinh(0.5 * V / Vmin);
+ }
+
+ double differential(double V) const {
+ return weight * (fd.C - normalStress * coefficientOfFriction(V));
+ }
+
+ double second_deriv(double V) const {
+ return weight * (-normalStress) * fd.a / std::hypot(2.0 * Vmin, V);
+ }
+
+ double regularity(double V) const { return std::abs(second_deriv(V)); }
+
+ void updateAlpha(double alpha) {
+ double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
+ Vmin = fd.V0 / std::exp(logrest);
+ }
+
+private:
+ FrictionData const fd;
+ double const weight;
+ double const normalStress;
+ double Vmin;
+};
+
class TrivialFunction : public FrictionPotential {
public:
double evaluate(double) const { return 0; }
diff --git a/src/assemblers.cc b/src/assemblers.cc
index 23602583..c2a9fae9 100644
--- a/src/assemblers.cc
+++ b/src/assemblers.cc
@@ -122,6 +122,7 @@ void MyAssembler<GridView, dimension>::assembleNormalStress(
template <class GridView, int dimension>
auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
+ Config::FrictionModel frictionModel,
BoundaryPatch<GridView> const &frictionalBoundary,
GlobalFrictionData<dimension> const &frictionInfo,
ScalarVector const &normalStress)
@@ -136,9 +137,18 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
weights, frictionalBoundary);
}
- return std::make_shared<
- GlobalRateStateFriction<Matrix, Vector, TruncatedRateState, GridView>>(
- frictionalBoundary, gridView, frictionInfo, weights, normalStress);
+ switch (frictionModel) {
+ case Config::Truncated:
+ return std::make_shared<GlobalRateStateFriction<
+ Matrix, Vector, TruncatedRateState, GridView>>(
+ frictionalBoundary, gridView, frictionInfo, weights, normalStress);
+ case Config::Regularised:
+ return std::make_shared<GlobalRateStateFriction<
+ Matrix, Vector, RegularisedRateState, GridView>>(
+ frictionalBoundary, gridView, frictionInfo, weights, normalStress);
+ default:
+ assert(false);
+ }
}
template <class GridView, int dimension>
diff --git a/src/assemblers.hh b/src/assemblers.hh
index 7b0c26c9..207bbde2 100644
--- a/src/assemblers.hh
+++ b/src/assemblers.hh
@@ -16,6 +16,8 @@
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/globalfrictiondata.hh>
+#include "enums.hh"
+
template <class GridView, int dimension> class MyAssembler {
public:
using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
@@ -76,6 +78,7 @@ template <class GridView, int dimension> class MyAssembler {
double poissonRatio, Vector const &displacement);
std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
+ Config::FrictionModel frictionModel,
BoundaryPatch<GridView> const &frictionalBoundary,
GlobalFrictionData<dimension> const &frictionInfo,
ScalarVector const &normalStress);
diff --git a/src/enum_friction_model.cc b/src/enum_friction_model.cc
new file mode 100644
index 00000000..e8a7d066
--- /dev/null
+++ b/src/enum_friction_model.cc
@@ -0,0 +1,15 @@
+#include <dune/common/exceptions.hh>
+
+#include "enums.hh"
+
+template <> struct StringToEnum<Config::FrictionModel> {
+ static Config::FrictionModel convert(std::string const &s) {
+ if (s == "Truncated")
+ return Config::Truncated;
+
+ if (s == "Regularised")
+ return Config::Regularised;
+
+ DUNE_THROW(Dune::Exception, "failed to parse enum");
+ }
+};
diff --git a/src/enums.hh b/src/enums.hh
index 7dbd9b95..5c04e4a4 100644
--- a/src/enums.hh
+++ b/src/enums.hh
@@ -2,6 +2,10 @@
#define SRC_ENUMS_HH
struct Config {
+ enum FrictionModel {
+ Truncated,
+ Regularised
+ };
enum stateModel {
AgeingLaw,
SlipLaw
diff --git a/src/sand-wedge-data/parset.cfg b/src/sand-wedge-data/parset.cfg
index 8a2b3bd9..9425327f 100644
--- a/src/sand-wedge-data/parset.cfg
+++ b/src/sand-wedge-data/parset.cfg
@@ -27,6 +27,7 @@ V0 = 5e-5 # [m/s]
L = 2e-5 # [m] ?
initialAlpha = 0.916290731874155 # [ ] ?
stateModel = AgeingLaw
+frictionModel = Truncated
[boundary.friction.weakening]
a = 0.015 # [1] ?
b = 0.030 # [1] ?
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index 8c562f87..d4053ad2 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -68,6 +68,7 @@
#include "assemblers.hh"
#include "tobool.hh"
#include "enum_parser.cc"
+#include "enum_friction_model.cc"
#include "enum_scheme.cc"
#include "enum_state_model.cc"
#include "enum_verbosity.cc"
@@ -248,6 +249,7 @@ int main(int argc, char *argv[]) {
body.getPoissonRatio(), u_initial);
MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"));
auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
+ parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
frictionalBoundary, frictionInfo, normalStress);
myGlobalFriction->updateAlpha(alpha_initial);
--
GitLab