diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh index 3ac28d2ade2a9c4ced7213872daa665099a4ac8a..de2d1dc53fc5d924700e166f4de189478241a370 100644 --- a/dune/tectonic/frictionpotential.hh +++ b/dune/tectonic/frictionpotential.hh @@ -29,8 +29,9 @@ class FrictionPotential { class TruncatedRateState : public FrictionPotential { public: - TruncatedRateState(double coefficient, double _normalStress, FrictionData _fd) - : fd(_fd), weight(coefficient), normalStress(_normalStress) {} + TruncatedRateState(double _weight, double _weightedNormalStress, + FrictionData _fd) + : fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {} double coefficientOfFriction(double V) const { if (V <= Vmin) @@ -40,14 +41,14 @@ class TruncatedRateState : public FrictionPotential { } double differential(double V) const { - return weight * (fd.C - normalStress * coefficientOfFriction(V)); + return weight * fd.C - weightedNormalStress * coefficientOfFriction(V); } double second_deriv(double V) const { if (V <= Vmin) return 0; - return weight * (-normalStress) * (fd.a / V); + return -weightedNormalStress * (fd.a / V); } double regularity(double V) const { @@ -65,26 +66,26 @@ class TruncatedRateState : public FrictionPotential { private: FrictionData const fd; double const weight; - double const normalStress; + double const weightedNormalStress; double Vmin; }; class RegularisedRateState : public FrictionPotential { public: - RegularisedRateState(double coefficient, double _normalStress, + RegularisedRateState(double _weight, double _weightedNormalStress, FrictionData _fd) - : fd(_fd), weight(coefficient), normalStress(_normalStress) {} + : fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {} double coefficientOfFriction(double V) const { return fd.a * std::asinh(0.5 * V / Vmin); } double differential(double V) const { - return weight * (fd.C - normalStress * coefficientOfFriction(V)); + return weight * fd.C - weightedNormalStress * coefficientOfFriction(V); } double second_deriv(double V) const { - return weight * (-normalStress) * fd.a / std::hypot(2.0 * Vmin, V); + return -weightedNormalStress * fd.a / std::hypot(2.0 * Vmin, V); } double regularity(double V) const { return std::abs(second_deriv(V)); } @@ -97,7 +98,7 @@ class RegularisedRateState : public FrictionPotential { private: FrictionData const fd; double const weight; - double const normalStress; + double const weightedNormalStress; double Vmin; }; diff --git a/dune/tectonic/globalratestatefriction.hh b/dune/tectonic/globalratestatefriction.hh index 1e7e83445717119afaf860be3797ae8c1c279d43..a43619c1b59fd38b07242cc2cf0f4fba24ea5735 100644 --- a/dune/tectonic/globalratestatefriction.hh +++ b/dune/tectonic/globalratestatefriction.hh @@ -26,10 +26,9 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> { GlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary, GridView const &gridView, GlobalFrictionData<block_size> const &frictionInfo, - // Note: passing the following two makes no sense ScalarVector const &weights, - ScalarVector const &normalStress) - : restrictions(normalStress.size()) { + ScalarVector const &weightedNormalStress) + : restrictions(weightedNormalStress.size()) { auto trivialNonlinearity = std::make_shared<Friction>(std::make_shared<TrivialFunction>()); @@ -44,7 +43,7 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> { continue; } auto const fp = std::make_shared<ScalarFriction>( - weights[i], normalStress[i], frictionInfo(coordinate)); + weights[i], weightedNormalStress[i], frictionInfo(coordinate)); restrictions[i] = std::make_shared<Friction>(fp); } } diff --git a/src/assemblers.cc b/src/assemblers.cc index 29ceb758f9c7eecda7539426aee415102e9113db..1dc5959399de87d2177762b9e33c1a24e5b7e091 100644 --- a/src/assemblers.cc +++ b/src/assemblers.cc @@ -142,15 +142,21 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity( vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler, weights, frictionalBoundary); } + ScalarVector weightedNormalStress(weights.size()); + for (size_t i = 0; i < weights.size(); ++i) + weightedNormalStress[i] = weights[i] * normalStress[i]; + switch (frictionModel) { case Config::Truncated: return std::make_shared<GlobalRateStateFriction< Matrix, Vector, TruncatedRateState, GridView>>( - frictionalBoundary, gridView, frictionInfo, weights, normalStress); + frictionalBoundary, gridView, frictionInfo, weights, + weightedNormalStress); case Config::Regularised: return std::make_shared<GlobalRateStateFriction< Matrix, Vector, RegularisedRateState, GridView>>( - frictionalBoundary, gridView, frictionInfo, weights, normalStress); + frictionalBoundary, gridView, frictionInfo, weights, + weightedNormalStress); default: assert(false); }