new friction models: Tresca, None; Tresca not tested

dune/tectonic/data-structures/enumparser.cc

@@ -32,6 +32,12 @@ Config::FrictionModel StringToEnum<Config::FrictionModel>::convert(
   if (s == "Regularised")
     return Config::Regularised;
 
+  if (s == "Tresca")
+    return Config::Tresca;
+
+  if (s == "None")
+    return Config::None;
+
   DUNE_THROW(Dune::Exception, "failed to parse enum");
 }
 

dune/tectonic/data-structures/enums.hh

@@ -2,7 +2,7 @@
 #define SRC_ENUMS_HH
 
 struct Config {
-  enum FrictionModel { Truncated, Regularised };
+  enum FrictionModel { Truncated, Regularised, Tresca, None };
   enum stateModel { AgeingLaw, SlipLaw };
   enum scheme { Newmark, BackwardEuler };
   enum PatchType { Rectangular, Trapezoidal };

dune/tectonic/data-structures/friction/frictionpotential.hh

@@ -34,10 +34,17 @@ class TruncatedRateState : public FrictionPotential {
       : fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}
 
   double coefficientOfFriction(double V) const override {
-    if (V <= Vmin)
+    if (V <= Vmin or Vmin == 0.0)
       return 0.0;
 
-    return fd.a * std::log(V / Vmin);
+    //std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;
+
+    auto res = fd.a * std::log(V / Vmin);
+    if (std::isinf(res)) {
+        std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;
+    }
+
+    return res;
   }
 
   double differential(double V) const override {
@@ -48,7 +55,7 @@ class TruncatedRateState : public FrictionPotential {
   double second_deriv(double V) const override {
     //std::cout << "second deriv: " << ((V<=Vmin) ? 0.0 : -weightedNormalStress * (fd.a / V)) << std::endl;
 
-    if (V <= Vmin)
+    if (V <= Vmin or Vmin == 0.0)
       return 0.0;
 
     return -weightedNormalStress * (fd.a / V);
@@ -67,7 +74,7 @@ class TruncatedRateState : public FrictionPotential {
 
     //std::cout << "second deriv: " << ((V<=Vmin) ? 0.0 : weight * fd.C * V - weightedNormalStress * fd.a * V * (std::log(V / Vmin) - 1)) << std::endl;
 
-    if (V <= Vmin)
+    if (V <= Vmin or Vmin == 0.0)
         return 0.0;
 
       return weight * fd.C * V - weightedNormalStress * fd.a * V * (std::log(V / Vmin) - 1);
@@ -76,7 +83,9 @@ class TruncatedRateState : public FrictionPotential {
   void updateAlpha(double alpha) override {
     double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
     Vmin = fd.V0 / std::exp(logrest);
-    //std::cout << "Vmin: " << Vmin << std::endl;
+
+    if (Vmin == 0.0)
+        std::cout << " logrest: " << logrest << " alpha: " << alpha << std::endl;
   }
 
 private:
@@ -124,6 +133,38 @@ class RegularisedRateState : public FrictionPotential {
   double Vmin;
 };
 
+class Tresca : public FrictionPotential {
+public:
+  Tresca(double _weight, double _weightedNormalStress, FrictionData _fd)
+      : fd(_fd), weightedNormalStress(_weightedNormalStress) {}
+
+  double coefficientOfFriction(double V) const override {
+    return fd.mu0;
+  }
+
+  double differential(double V) const override {
+    return - weightedNormalStress * coefficientOfFriction(V);
+  }
+
+  double second_deriv(double V) const override {
+    return 0;
+  }
+
+  double regularity(double V) const override {
+    return std::abs(second_deriv(V));
+  }
+
+  double evaluate(double V) const override {
+      return 0;
+  }
+
+  void updateAlpha(double alpha) override {}
+
+private:
+  FrictionData const fd;
+  double const weightedNormalStress;
+};
+
 class ZeroFunction : public FrictionPotential {
 public:
   template <typename... Args>

dune/tectonic/data-structures/network/contactnetwork.cc

@@ -162,6 +162,16 @@ void ContactNetwork<HostGridType, VectorType>::assembleFriction(
               Matrix, VectorType, RegularisedRateState, GridType>>(
               nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
           break;
+        case Config::Tresca:
+          globalFriction_ = std::make_shared<GlobalRateStateFriction<
+            Matrix, VectorType, Tresca, GridType>>(
+            nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+          break;
+        case Config::None:
+          globalFriction_ = std::make_shared<GlobalRateStateFriction<
+            Matrix, VectorType, ZeroFunction, GridType>>(
+            nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+          break;
         default:
           assert(false);
           break;