TODOs / comments

dune/tectonic/localnonlinearity.hh

@@ -81,6 +81,7 @@ template <int dimension> class LocalNonlinearity {
       A[k][k] += (h2 - h1ox) * x[k] * x[k] / normX2 + h1ox;
   }
 
+  // TODO: do not evaluate at zero
   void addGradient(VectorType const &x, VectorType &y) const {
     if (x == VectorType(0))
       return;
@@ -90,11 +91,13 @@ template <int dimension> class LocalNonlinearity {
     y += tmp;
   }
 
+  // TODO: do not evaluate at zero
   void upperGradient(VectorType const &x, VectorType &ret) const {
     ret = x;
     ret *= func_->rightDifferential(x.two_norm()) / x.two_norm();
   }
 
+  // TODO: do not evaluate at zero
   void lowerGradient(VectorType const &x, VectorType &ret) const {
     ret = x;
     ret *= func_->leftDifferential(x.two_norm()) / x.two_norm();

dune/tectonic/myblockproblem.hh

@@ -86,6 +86,7 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
     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
@@ -102,6 +103,7 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
     // construct sparsity pattern for linearization
     Dune::MatrixIndexSet indices(problem.A.N(), problem.A.M());
     indices.import(problem.A);
+    // TODO: handle smooth domain
     // problem.phi.addHessianIndices(indices);
 
     // construct matrix from pattern and initialize it

src/TODO.org

+* compute damping parameter
+* investigate why laursen fails
+* handle smooth domain properly
+* fix up octave bindings (low-pri)

src/one-body-sample.cc

@@ -244,13 +244,13 @@ int main(int argc, char *argv[]) {
     nonlinearGSStep.ignoreNodes_ = &ignoreNodes;
 
     VectorType u1(grid.size(grid.maxLevel(), dim));
-    u1 = 0;
+    u1 = 0.0; // Has to be zero!
     VectorType u2 = u1;
     VectorType u3 = u1;
     VectorType u4 = u1;
 
     VectorType u1_diff_new(grid.size(grid.maxLevel(), dim));
-    u1_diff_new = 0;
+    u1_diff_new = 0.0; // Has to be zero!
     VectorType u2_diff_new = u1_diff_new;
     VectorType u3_diff_new = u1_diff_new;
     VectorType u4_diff_new = u1_diff_new;

src/one-body-sample.parset

@@ -39,7 +39,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
+model = Ruina # TODO: Laursen doesn't work with TNNMG
 
 [boundary.friction.ruina]
 # "For rocks, typical values of A and B range from 0.005 to 0.015"