Handle irregularity

dune/tectonic/globalnonlinearity.hh

@@ -33,6 +33,11 @@ class GlobalNonlinearity {
       res->addGradient(v[i], gradient[i]);
     }
   }
+
+  double regularity(int i, typename VectorType::block_type const &x) const {
+    auto res = restriction(i);
+    return res->regularity(x);
+  }
 };
 }
 #endif

dune/tectonic/localnonlinearity.hh

@@ -26,6 +26,13 @@ template <int dimension> class LocalNonlinearity {
     return ret;
   }
 
+  double regularity(VectorType const &x) const {
+    if (x.two_norm() < 1e-8) // TODO
+      return std::numeric_limits<double>::infinity();
+
+    return func_->regularity(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 &u, VectorType const &v,
@@ -83,11 +90,12 @@ template <int dimension> class LocalNonlinearity {
 
   // TODO: do not evaluate at zero
   void addGradient(VectorType const &x, VectorType &y) const {
-    if (x == VectorType(0))
+    if (x.two_norm() <
+        1e-8 or std::isinf(func_->regularity(x.two_norm()))) // TODO
       return;
 
     VectorType tmp;
-    upperGradient(x, tmp); // TODO
+    upperGradient(x, tmp); // upper and lower gradient coincide in this case
     y += tmp;
   }
 

dune/tectonic/myblockproblem.hh

@@ -83,22 +83,9 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
     // determine truncation pattern
     linearization.truncation.resize(u.size());
     linearization.truncation.unsetAll();
-    Interval<double> ab;
-    for (int i = 0; i < u.size(); ++i) {
-      for (int j = 0; j < block_size; ++j) {
-        // TODO: handle smooth domain
-        // problem.phi.smoothDomain(i, u[i][j], regularity_tol, ab, j);
-
-        // // if (phi is not regular at u) or (u is not in smooth domain) or
-        // (component should be ignored)
-        // // truncate this component
-        // if ((problem.phi.regularity(i, u[i][j], j) > regularity_tol)
-        //     or (u[i][j] < ab[0])
-        //     or (u[i][j] > ab[1])
-        //     or ignore[i][j])
-        //   linearization.truncation[i][j] = true;
-      }
-    }
+    for (int i = 0; i < u.size(); ++i)
+      if (problem.phi.regularity(i, u[i]) > 1e8) // TODO
+        linearization.truncation[i] = true;
 
     // construct sparsity pattern for linearization
     Dune::MatrixIndexSet indices(problem.A.N(), problem.A.M());

dune/tectonic/nicefunction.hh

@@ -18,6 +18,10 @@ class NiceFunction : public VirtualFunction<double, double> {
   double virtual second_deriv(double x) const {
     DUNE_THROW(NotImplemented, "second derivative not implemented");
   }
+
+  double virtual regularity(double s) const {
+    DUNE_THROW(NotImplemented, "regularity not implemented");
+  }
 };
 
 class RuinaFunction : public NiceFunction {
@@ -82,14 +86,20 @@ class RuinaFunction : public NiceFunction {
     = a/x
   */
   double virtual second_deriv(double s) const {
-    if (eta * s < rho)
+    // includes the case eta * s = rho for which there is no second derivative
+    if (eta * s <= rho)
       return 0;
-    else if (eta * s == rho)
-      return 37; // TODO: not differentiable;
     else
       return coefficient * normalStress * (a / s);
   }
 
+  double virtual regularity(double s) const {
+    if (eta * s == rho)
+      return std::numeric_limits<double>::infinity();
+
+    return std::abs(second_deriv(s));
+  }
+
 private:
   double coefficient;
   double a;
@@ -116,6 +126,8 @@ class LinearFunction : public NiceFunction {
 
   double virtual second_deriv(double) const { return 0; }
 
+  double virtual regularity(double s) const { return 0; }
+
 private:
   double coefficient;
 };
@@ -153,6 +165,8 @@ class TrivialFunction : public NiceFunction {
   double virtual rightDifferential(double) const { return 0; }
 
   double virtual second_deriv(double) const { return 0; }
+
+  double virtual regularity(double) const { return 0; }
 };
 
 // slope in [n-1,n] is n

src/one-body-sample.parset

@@ -49,7 +49,7 @@ normalstress = 0.1 # laursen depends a lot more on this
 #    http://earthquake.usgs.gov/research/physics/lab/prediction.pdf
 mu = 0.5
 eta = 1
-model = Ruina # TODO: Laursen doesn't work with TNNMG
+model = Ruina
 
 [boundary.friction.ruina]
 # "For rocks, typical values of A and B range from 0.005 to 0.015"