From 7bfbe7e8137476172639366a1578f66c4c7d3dfc Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Fri, 26 Feb 2021 15:27:04 +0100
Subject: [PATCH] relative velocities in mortar basis yield small but non-zero
normal contributions, remove them here
---
.../data-structures/friction/localfriction.hh | 33 +++++++++++++------
1 file changed, 23 insertions(+), 10 deletions(-)
diff --git a/dune/tectonic/data-structures/friction/localfriction.hh b/dune/tectonic/data-structures/friction/localfriction.hh
index 461dcea9..bbd78ce3 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 &,
--
GitLab