[Extend] Allow for cohesion

dune/tectonic/frictiondata.hh

@@ -3,9 +3,11 @@
 #include <dune/common/parametertree.hh>
 
 struct FrictionData {
-  explicit FrictionData(double _L, double _V0, double _a, double _b, double _mu0)
-      : L(_L), V0(_V0), a(_a), b(_b), mu0(_mu0) {}
+  explicit FrictionData(double _C, double _L, double _V0, double _a, double _b,
+                        double _mu0)
+      : C(_C), L(_L), V0(_V0), a(_a), b(_b), mu0(_mu0) {}
 
+  double const C;
   double const L;
   double const V0;
   double const a;

dune/tectonic/frictionpotential.hh

@@ -41,7 +41,7 @@ class FrictionPotential : public FrictionPotentialWrapper {
   }
 
   double differential(double V) const {
-    return weight * (-normalStress) * coefficientOfFriction(V);
+    return weight * (fd.C - normalStress * coefficientOfFriction(V));
   }
 
   double second_deriv(double V) const {

dune/tectonic/globalfrictiondata.hh

@@ -18,7 +18,7 @@ double evaluateScalarFunction(
 template <int dimension> class GlobalFrictionData {
 public:
   FrictionData operator()(Dune::FieldVector<double, dimension> const &x) const {
-    return FrictionData(L(), V0(), evaluateScalarFunction(a(), x),
+    return FrictionData(C(), L(), V0(), evaluateScalarFunction(a(), x),
                         evaluateScalarFunction(b(), x), mu0());
   }
 
@@ -26,6 +26,7 @@ template <int dimension> class GlobalFrictionData {
   using VirtualFunction = Dune::VirtualFunction<
       Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
 
+  double virtual const &C() const = 0;
   double virtual const &L() const = 0;
   double virtual const &V0() const = 0;
   VirtualFunction virtual const &a() const = 0;

src/myglobalfrictiondata.hh

@@ -17,12 +17,14 @@ class MyGlobalFrictionData : public GlobalFrictionData<dimension> {
 
 public:
   MyGlobalFrictionData(Dune::ParameterTree const &parset, MyGeometry const &tri)
-      : L_(parset.get<double>("L")),
+      : C_(parset.get<double>("C")),
+        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")) {}
 
+  double virtual const &C() const override { return C_; }
   double virtual const &L() const override { return L_; }
   double virtual const &V0() const override { return V0_; }
   VirtualFunction virtual const &a() const override { return a_; }
@@ -33,6 +35,7 @@ class MyGlobalFrictionData : public GlobalFrictionData<dimension> {
   using MyConstantFunction = ConstantFunction<
       Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
 
+  double const C_;
   double const L_;
   double const V0_;
   MyConstantFunction const a_;

src/one-body-sample.parset

@@ -16,6 +16,7 @@ shearViscosity = 0
 bulkViscosity  = 0
 
 [boundary.friction]
+C            = 0.0
 mu0          = 0.6
 a            = 0.010
 b            = 0.015