diff --git a/dune/tectonic/frictiondata.hh b/dune/tectonic/frictiondata.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8cc724698a9c163e7a5c157dbef4a2179e1f3817
--- /dev/null
+++ b/dune/tectonic/frictiondata.hh
@@ -0,0 +1,21 @@
+#ifndef FRICTIONDATA_HH
+#define FRICTIONDATA_HH
+#include <dune/common/parametertree.hh>
+
+struct FrictionData {
+ FrictionData(Dune::ParameterTree const &parset, double normalStress)
+ : 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) {}
+
+ double const L;
+ double const V0;
+ double const a;
+ double const b;
+ double const mu0;
+ double const normalStress;
+};
+#endif
diff --git a/dune/tectonic/globalruinanonlinearity.hh b/dune/tectonic/globalruinanonlinearity.hh
index b64d8bba57a970eb44749181e31d6a321d1674e2..bf7565f02c7d6acacda20fb35394d5afafb10694 100644
--- a/dune/tectonic/globalruinanonlinearity.hh
+++ b/dune/tectonic/globalruinanonlinearity.hh
@@ -23,9 +23,8 @@ class GlobalRuinaNonlinearity
public:
using GlobalNonlinearity<MatrixType, VectorType>::dim;
- GlobalRuinaNonlinearity(dataref nodalIntegrals, double a, double mu0,
- double V0, double normalStress, double b,
- dataref state, double L)
+ GlobalRuinaNonlinearity(dataref nodalIntegrals, FrictionData const &fd,
+ dataref state)
: restrictions(nodalIntegrals.size()) {
auto trivialNonlinearity = make_shared<LocalNonlinearity<dim> const>(
make_shared<TrivialFunction const>());
@@ -34,8 +33,7 @@ class GlobalRuinaNonlinearity
? trivialNonlinearity
: make_shared<LocalNonlinearity<dim> const>(
make_shared<RuinaFunction const>(
- nodalIntegrals[i], a, mu0, V0,
- normalStress, b, state[i], L));
+ nodalIntegrals[i], fd, state[i]));
}
}
diff --git a/dune/tectonic/nicefunction.hh b/dune/tectonic/nicefunction.hh
index 9a2446da6f99808bd06f480a7d0c42f563de29a7..7943a2c71e454bf3a5671042e6ca2e8134dcd2fe 100644
--- a/dune/tectonic/nicefunction.hh
+++ b/dune/tectonic/nicefunction.hh
@@ -9,6 +9,8 @@
#include <dune/common/exceptions.hh>
#include <dune/common/function.hh>
+#include "frictiondata.hh"
+
namespace Dune {
class NiceFunction {
protected:
@@ -45,12 +47,13 @@ class NiceFunction {
// i.e. K = -a log(V_m / V_0) = mu_0 + b log(V_0 / L) + b alpha
class RuinaFunction : public NiceFunction {
public:
- RuinaFunction(double coefficient, double a, double mu_0, double V_0,
- double normalStress, double b, double state, double L)
+ RuinaFunction(double coefficient, FrictionData const &fd, double state)
: NiceFunction(),
- coefficientProduct(coefficient * a * normalStress),
+ coefficientProduct(coefficient * fd.a * fd.normalStress),
// state is assumed to be logarithmic
- V_m(V_0 * std::exp(-(mu_0 + b * (state + std::log(V_0 / L))) / a))
+ V_m(fd.V0 *
+ std::exp(-(fd.mu0 + fd.b * (state + std::log(fd.V0 / fd.L))) /
+ fd.a))
// We could also compute V_m as
// V_0 * std::exp(-(mu_0 + b * state)/a)
// * std::pow(V_0 / L, -b/a)
diff --git a/src/assemblers.cc b/src/assemblers.cc
index aa19a12201815ad5613da8a6d870b6783e8a4dbd..49539136161e606acaa472548de3f9d13a335ae1 100644
--- a/src/assemblers.cc
+++ b/src/assemblers.cc
@@ -52,11 +52,11 @@ Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType> const>
assemble_nonlinearity(
Dune::ParameterTree const &parset,
Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
- double a, double mu0, double V0, double normalStress, double b,
- Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double L) {
+ FrictionData const &fd,
+ Dune::BlockVector<Dune::FieldVector<double, 1>> const &state) {
return Dune::make_shared<
Dune::GlobalRuinaNonlinearity<MatrixType, VectorType> const>(
- nodalIntegrals, a, mu0, V0, normalStress, b, state, L);
+ nodalIntegrals, fd, state);
}
#include "assemblers_tmpl.cc"
diff --git a/src/assemblers.hh b/src/assemblers.hh
index f29d6fad176a6e9bb52ead930bae684ded676204..5f5f27ea1f05f8eaa9fb042979eac9f13ffcb17f 100644
--- a/src/assemblers.hh
+++ b/src/assemblers.hh
@@ -29,6 +29,6 @@ Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType> const>
assemble_nonlinearity(
Dune::ParameterTree const &parset,
Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
- double a, double mu0, double V0, double normalStress, double b,
- Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double L);
+ FrictionData const &fd,
+ Dune::BlockVector<Dune::FieldVector<double, 1>> const &state);
#endif
diff --git a/src/assemblers_tmpl.cc b/src/assemblers_tmpl.cc
index 9a7fcc68a1f1eaec2bfb29e3ad284ae8bfe2417b..9faff8ebcd6b47d0ec43d7c0eae0c51be17fd11d 100644
--- a/src/assemblers_tmpl.cc
+++ b/src/assemblers_tmpl.cc
@@ -36,5 +36,5 @@ template Dune::shared_ptr<
assemble_nonlinearity<MatrixType, VectorType>(
Dune::ParameterTree const &parset,
Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
- double a, double mu0, double V0, double normalStress, double b,
- Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double L);
+ FrictionData const &fd,
+ Dune::BlockVector<Dune::FieldVector<double, 1>> const &state);
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 9d0deb709013c7146877bdb50949dfa3ae1af00b..ae0a0e4e7af9235be7643e7afae5923547dfee8d 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -116,16 +116,17 @@ template <class SingletonVectorType, class VectorType>
Dune::shared_ptr<StateUpdater<SingletonVectorType, VectorType>>
initStateUpdater(Config::state_model model,
SingletonVectorType const &alpha_initial,
- Dune::BitSetVector<1> const &frictionalNodes, double L) {
+ Dune::BitSetVector<1> const &frictionalNodes,
+ FrictionData const &fd) {
switch (model) {
case Config::Dieterich:
return Dune::make_shared<
DieterichStateUpdater<SingletonVectorType, VectorType>>(
- alpha_initial, frictionalNodes, L);
+ alpha_initial, frictionalNodes, fd.L);
case Config::Ruina:
return Dune::make_shared<
RuinaStateUpdater<SingletonVectorType, VectorType>>(
- alpha_initial, frictionalNodes, L);
+ alpha_initial, frictionalNodes, fd.L);
default:
assert(false);
}
@@ -142,10 +143,9 @@ template <class FunctionMap> void initPython(FunctionMap &functions) {
.toC<typename FunctionMap::Base>(functions);
}
-double computeCOF(double mu0, double a, double b, double V0, double L, double V,
- double log_state) {
- double const mu =
- mu0 + a * std::log(V / V0) + b * (log_state + std::log(V0 / L));
+double computeCOF(FrictionData const &fd, double V, double log_state) {
+ double const mu = fd.mu0 + fd.a * std::log(V / fd.V0) +
+ fd.b * (log_state + std::log(fd.V0 / fd.L));
return (mu > 0) ? mu : 0;
}
@@ -168,11 +168,6 @@ int main(int argc, char *argv[]) {
auto const E = parset.get<double>("body.E");
auto const nu = parset.get<double>("body.nu");
- auto const mu0 = parset.get<double>("boundary.friction.mu0");
- auto const a = parset.get<double>("boundary.friction.a");
- auto const b = parset.get<double>("boundary.friction.b");
- auto const V0 = parset.get<double>("boundary.friction.V0");
- auto const L = parset.get<double>("boundary.friction.L");
// {{{ Set up grid
using GridType =
@@ -294,6 +289,8 @@ int main(int argc, char *argv[]) {
gravityFunctional);
}
}
+ FrictionData const frictionData(parset.sub("boundary.friction"),
+ normalStress);
MatrixType stiffnessMatrix;
EnergyNorm<MatrixType, VectorType> const stiffnessMatrixNorm(
@@ -368,7 +365,7 @@ int main(int argc, char *argv[]) {
stiffnessMatrix, u_initial, ud_initial, udd_initial);
auto stateUpdater = initStateUpdater<SingletonVectorType, VectorType>(
parset.get<Config::state_model>("boundary.friction.state_model"),
- alpha_initial, frictionalNodes, L);
+ alpha_initial, frictionalNodes, frictionData);
auto const timesteps = parset.get<size_t>("timeSteps");
auto const tau = parset.get<double>("endOfTime") / timesteps;
@@ -412,8 +409,8 @@ int main(int argc, char *argv[]) {
SingletonVectorType const &_alpha) {
auto myGlobalNonlinearity =
assemble_nonlinearity<MatrixType, VectorType>(
- parset.sub("boundary.friction"), *nodalIntegrals, a, mu0, V0,
- normalStress, b, _alpha, L);
+ parset.sub("boundary.friction"), *nodalIntegrals, frictionData,
+ _alpha);
using MyConvexProblemType = MyConvexProblem<MatrixType, VectorType>;
MyConvexProblemType const myConvexProblem(
@@ -487,7 +484,7 @@ int main(int argc, char *argv[]) {
stateWriter << alpha[i][0] << " ";
displacementWriter << u[i][0] << " ";
velocityWriter << ud[i][0] << " ";
- coefficientWriter << computeCOF(mu0, a, b, V0, L, ud[i].two_norm(),
+ coefficientWriter << computeCOF(frictionData, ud[i].two_norm(),
alpha[i]) << " ";
}