Do not distinguish between left- and right gradient

dune/tectonic/ellipticenergy.hh

@@ -52,9 +52,9 @@ template <int dim> class EllipticEnergy {
   // returns false if the direction is tangential
   void descentDirectionNew(SmallVector const &x, SmallVector &ret) const {
     SmallVector pg;
-    upperGradient(x, pg);
+    gradient(x, pg);
     SmallVector mg;
-    lowerGradient(x, mg);
+    gradient(x, mg);
     double const pgx = pg * x;
     double const mgx = mg * x;
     if (pgx >= 0 && mgx >= 0) {
@@ -78,9 +78,9 @@ template <int dim> class EllipticEnergy {
     }
 
     SmallVector pg;
-    upperGradient(x, pg);
+    gradient(x, pg);
     SmallVector mg;
-    lowerGradient(x, mg);
+    gradient(x, mg);
     double const pgx = pg * x;
     double const mgx = mg * x;
     if (pgx >= 0 && mgx >= 0) {
@@ -118,17 +118,8 @@ template <int dim> class EllipticEnergy {
       y[ignore] = 0;
   }
 
-  void upperGradient(SmallVector const &x, SmallVector &y) const {
-    phi->upperGradient(x, y);
-    SmallVector z;
-    smoothGradient(x, z);
-    y += z;
-    if (ignore != dim)
-      y[ignore] = 0;
-  }
-
-  void lowerGradient(SmallVector const &x, SmallVector &y) const {
-    phi->lowerGradient(x, y);
+  void gradient(SmallVector const &x, SmallVector &y) const {
+    phi->gradient(x, y);
     SmallVector z;
     smoothGradient(x, z);
     y += z;

dune/tectonic/frictionpotential.hh

@@ -16,9 +16,7 @@ class FrictionPotentialWrapper {
 public:
   virtual ~FrictionPotentialWrapper() {}
 
-  double virtual leftDifferential(double s) const = 0;
-  double virtual rightDifferential(double s) const = 0;
-
+  double virtual differential(double s) const = 0;
   double virtual second_deriv(double x) const = 0;
 
   double virtual regularity(double s) const = 0;
@@ -61,7 +59,7 @@ class FrictionPotential : public FrictionPotentialWrapper {
 
   // log(V/V_m)  if V >= V_0
   // 0           otherwise
-  double virtual leftDifferential(double V) const {
+  double virtual differential(double V) const {
     assert(V >= 0);
     if (V <= V_m)
       return 0;
@@ -69,10 +67,6 @@ class FrictionPotential : public FrictionPotentialWrapper {
     return coefficientProduct * std::log(V / V_m);
   }
 
-  double virtual rightDifferential(double V) const {
-    return leftDifferential(V);
-  }
-
   // 1/V        if V >  V_0
   // undefined  if V == V_0
   // 0          if V <  V_0
@@ -104,9 +98,7 @@ class TrivialFunction : public FrictionPotentialWrapper {
 public:
   double virtual evaluate(double) const { return 0; }
 
-  double virtual leftDifferential(double) const { return 0; }
-
-  double virtual rightDifferential(double) const { return 0; }
+  double virtual differential(double) const { return 0; }
 
   double virtual second_deriv(double) const { return 0; }
 

dune/tectonic/localfriction.hh

@@ -35,21 +35,14 @@ template <int dimension> class LocalFriction {
 
   // directional subdifferential: at u on the line u + t*v
   // u and v are assumed to be non-zero
-  void directionalSubDiff(VectorType const &u, VectorType const &v,
+  void directionalSubDiff(VectorType const &x, VectorType const &v,
                           Interval<double> &D) const {
-    if (u.two_norm() == 0) {
-      D[0] = D[1] = func->rightDifferential(0) * v.two_norm();
+    if (x.two_norm() == 0) {
+      D[0] = D[1] = func->differential(0) * v.two_norm();
       return;
     }
-    double const un = u.two_norm();
-    double const ndotp = (u * v) / un;
-    if (ndotp > 0) {
-      D[1] = ndotp * func->rightDifferential(un);
-      D[0] = ndotp * func->leftDifferential(un);
-    } else {
-      D[1] = ndotp * func->leftDifferential(un);
-      D[0] = ndotp * func->rightDifferential(un);
-    }
+    double const un = x.two_norm();
+    D[0] = D[1] = func->differential(un) * (x * v) / un;
   }
 
   /** Formula for the derivative:
@@ -73,7 +66,7 @@ template <int dimension> class LocalFriction {
     double const xnorm = std::sqrt(xnorm2);
     double const xnorm3 = xnorm * xnorm2;
 
-    double const H1 = func->rightDifferential(xnorm);
+    double const H1 = func->differential(xnorm);
     double const H2 = func->second_deriv(xnorm);
 
     // TODO: potential optimisation: factor out (H1 / xnorm), get rid of xnorm3
@@ -120,24 +113,12 @@ template <int dimension> class LocalFriction {
     if (xnorm == 0 or std::isinf(func->regularity(xnorm)))
       return;
 
-    // left and right differential coincide in this case
-    y.axpy(func->rightDifferential(xnorm) / xnorm, x);
+    y.axpy(func->differential(xnorm) / xnorm, x);
   }
 
-  void upperGradient(VectorType const &x, VectorType &ret) const {
-    double const xnorm = x.two_norm();
-    assert(xnorm != 0);
-
-    ret = x;
-    ret *= func->rightDifferential(xnorm) / xnorm;
-  }
-
-  void lowerGradient(VectorType const &x, VectorType &ret) const {
-    double const xnorm = x.two_norm();
-    assert(xnorm != 0);
-
-    ret = x;
-    ret *= func->leftDifferential(xnorm) / xnorm;
+  void gradient(VectorType const &x, VectorType &ret) const {
+    ret = 0;
+    addGradient(x, ret);
   }
 
   void directionalDomain(VectorType const &, VectorType const &,