relative velocities in mortar basis yield small but non-zero normal contributions, remove them here

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 &,