diff --git a/dune/tectonic/data-structures/friction/localfriction.hh b/dune/tectonic/data-structures/friction/localfriction.hh index 461dcea90abc643b3478db2370e15525eaa2f5dc..bbd78ce36f7c84172fbb4f817a28db5c6ca9f16c 100644 --- a/dune/tectonic/data-structures/friction/localfriction.hh +++ b/dune/tectonic/data-structures/friction/localfriction.hh @@ -40,19 +40,28 @@ class WrappedScalarFriction : public LocalFriction<dimension> { using VectorType = typename LocalFriction<dimension>::VectorType; using MatrixType = typename LocalFriction<dimension>::MatrixType; +private: + VectorType removeNormal(const VectorType& x) { + VectorType res = x; + res[0] = 0.0; + return res; + } + public: template <typename... Args> WrappedScalarFriction(Args... args) : func_(args...) {} double operator()(VectorType const &x) const override { - return func_.evaluate(x.two_norm()); + auto tangential_x = removeNormal(x); + return func_.evaluate(tangential_x.two_norm()); } void updateAlpha(double alpha) override { func_.updateAlpha(alpha); } double regularity(VectorType const &x) const override { - double const xnorm = x.two_norm(); + auto tangential_x = removeNormal(x); + double const xnorm = tangential_x.two_norm(); if (xnorm <= 0.0) return std::numeric_limits<double>::infinity(); @@ -60,18 +69,21 @@ class WrappedScalarFriction : public LocalFriction<dimension> { } double coefficientOfFriction(VectorType const &x) const override { - return func_.coefficientOfFriction(x.two_norm()); + auto tangential_x = removeNormal(x); + return func_.coefficientOfFriction(tangential_x.two_norm()); } // directional subdifferential: at u on the line u + t*v // u and v are assumed to be non-zero void directionalSubDiff(VectorType const &x, VectorType const &v, Dune::Solvers::Interval<double> &D) const override { - double const xnorm = x.two_norm(); + + auto tangential_x = removeNormal(x); + double const xnorm = tangential_x.two_norm(); if (xnorm <= 0.0) D[0] = D[1] = func_.differential(0.0) * v.two_norm(); else - D[0] = D[1] = func_.differential(xnorm) * (x * v) / xnorm; + D[0] = D[1] = func_.differential(xnorm) * (tangential_x * v) / xnorm; } /** Formula for the derivative: @@ -122,12 +134,12 @@ class WrappedScalarFriction : public LocalFriction<dimension> { double const entry = tensorweight * x[k] * x[k]; A[k][k] += entry + idweight; }*/ - - const double xnorm = x.two_norm(); + auto tangential_x = removeNormal(x); + const double xnorm = tangential_x.two_norm(); if (xnorm <= 0.0) return; - VectorType y = x; + VectorType y = tangential_x; y /= xnorm; double const H1 = func_.differential(xnorm); @@ -154,12 +166,13 @@ class WrappedScalarFriction : public LocalFriction<dimension> { } void addGradient(VectorType const &x, VectorType &y) const override { - double const xnorm = x.two_norm(); + auto tangential_x = removeNormal(x); + double const xnorm = tangential_x.two_norm(); if (std::isinf(func_.regularity(xnorm))) return; if (xnorm > 0.0) - Arithmetic::addProduct(y, func_.differential(xnorm) / xnorm, x); + Arithmetic::addProduct(y, func_.differential(xnorm) / xnorm, tangential_x); } void directionalDomain(VectorType const &, VectorType const &,