Add frictionData struct

dune/tectonic/frictiondata.hh

+#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

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]));
     }
   }
 

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)

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"

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

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);

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]) << " ";
         }