debug print

dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh

@@ -26,23 +26,25 @@ class LocalBisectionSolver
           return;
       }
 
+      auto old_x = x;
       auto origin = x;
       auto linearPart = f.linearPart();
 
 
-      /*std::cout << x << std::endl;
-      std::cout << "eigenvalue: " << maxEigenvalue << std::endl;
-      std::cout << "ignore: " << ignore << std::endl;
-      std::cout << "orig linear part: " << linearPart << std::endl;*/
+      //std::cout << x << std::endl;
+      //std::cout << "eigenvalue: " << f.quadraticPart()[0][0] << std::endl;
+      //std::cout << "ignore: " << ignore << std::endl;
+      //std::cout << "orig linear part: " << linearPart << std::endl;
 
       auto lower = f.lowerObstacle();
       auto upper = f.upperObstacle();
       for (size_t j = 0; j < lower.size(); ++j) {
         if (std::abs(upper[j]-lower[j])<= safety) {
-          linearPart[j] = 0;
+          //linearPart[j] = 0;
         }
       }
-      std::cout << "lower: " << lower << " upper: " << upper << std::endl;
+
+      //std::cout << "lower: " << lower << " upper: " << upper << std::endl;
 
       for (size_t j = 0; j < ignore.size(); ++j) {
         if (ignore[j]) {
@@ -52,15 +54,27 @@ class LocalBisectionSolver
         }
       }
 
+      auto linearPartNorm = linearPart.two_norm();
+      if (linearPartNorm > 0.0)
+          linearPart /= linearPartNorm;
+      else {
+          linearPartNorm = 0.0;
+          linearPart = 0.0;
+      }
+
+      //using DirectionalRestriction = DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>;
+      //DirectionalRestriction fv(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), f.origin(), linearPart, linearPartNorm);
+
       auto fv = directionalRestriction(f, x, linearPart);
 
       //std::cout << "origin: " << origin << std::endl;
       //std::cout << " linearPart: " << linearPart << std::endl;
 
-      //std::cout << "maxEigenvalue " << maxEigenvalue << std::endl;
-      //std::cout << "linearPartNorm " << linearPartNorm << std::endl;
-      //print(origin, "origin:");
-      //print(linearPart, "linearPart:");
+      std::cout << "maxEigenvalue " << f.quadraticPart() << std::endl;
+      print(linearPart, "linearPart:");
+      std::cout << "linearPartNorm " << linearPartNorm << std::endl;
+      print(origin, "origin:");
+
 
 
       //print(lower, "lower:");
@@ -98,14 +112,30 @@ class LocalBisectionSolver
       //std::cout << "domain: " << firstOrderFunctional.domain()[0] << " " << firstOrderFunctional.domain()[1] << std::endl;
 
       typename Functional::Range alpha;
-      LineSearchSolver lineSearch;
-      lineSearch(alpha, fv, false);
 
-      //std::cout << "alpha: " << alpha << std::endl;
+      if (fv.scaling()<=0.0) {
+        alpha = 0.0;
+      } /* else {
+        LineSearchSolver lineSearch;
+        lineSearch(alpha, fv, false);
+      }*/
+
+      const auto& domain = fv.domain();
+      if (std::abs(domain[0]-domain[1])>safety) {
+        int bisectionsteps = 0;
+        const Bisection globalBisection(0.0, 1.0, 0.0, 0.0);
+
+        alpha = globalBisection.minimize(fv, 0.0, 0.0, bisectionsteps);
+      } else {
+          alpha = domain[0];
+      }
+
+
+      std::cout << "alpha: " << alpha << std::endl;
 
       linearPart *= alpha;
 
-      //std::cout << linearPart << std::endl;
+      std::cout << linearPart << std::endl;
 
 #ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
       if (std::abs(alpha)> safety) {
@@ -115,13 +145,14 @@ class LocalBisectionSolver
         x = origin;
       }
 #else
+      std::cout << "use old" << std::endl;
       if (std::abs(alpha)> safety) {
         x = origin;
         x += linearPart;
-        //x -= f.origin();
+        x -= old_x;//f.origin();
       }
 #endif
-      //std::cout << x << std::endl << std::endl;
+      std::cout << "new x: "<< x << std::endl << std::endl;
   }
 };