diff --git a/dune/tectonic/data-structures/enumparser.cc b/dune/tectonic/data-structures/enumparser.cc index 46c2e24b8a8aae44660ec16629adea51a0686f42..eafb98e16f6b8eb9125ca50dd7b8bf80397afa4a 100644 --- a/dune/tectonic/data-structures/enumparser.cc +++ b/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"); } diff --git a/dune/tectonic/data-structures/enums.hh b/dune/tectonic/data-structures/enums.hh index 0ec8540dfe171fb667a18dabff0570140e7c9969..d8393c7874d84a71bbb4c0fde065ff1bc74f6cb6 100644 --- a/dune/tectonic/data-structures/enums.hh +++ b/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 }; diff --git a/dune/tectonic/data-structures/friction/frictionpotential.hh b/dune/tectonic/data-structures/friction/frictionpotential.hh index 8b98b22b11c46637389fb6875d7e8fac3c804069..becfb81d31207cc737cbc6964666c69957db70eb 100644 --- a/dune/tectonic/data-structures/friction/frictionpotential.hh +++ b/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> diff --git a/dune/tectonic/data-structures/network/contactnetwork.cc b/dune/tectonic/data-structures/network/contactnetwork.cc index 7e03c5496f83e6226a2a61d67828e2393a46b801..1e55592d0c5fa95d2d0f412dee71bbd98c31235c 100644 --- a/dune/tectonic/data-structures/network/contactnetwork.cc +++ b/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;