adapt to removal of FirstOrderFunctional, use LineSearchSolver, remove contact...

adapt to removal of FirstOrderFunctional, use LineSearchSolver, remove contact obstacle dofs from direction

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

@@ -26,47 +26,55 @@ class LocalBisectionSolver
           return;
       }
 
-      auto maxEigenvalue = f.maxEigenvalues();
+      auto origin = x;
+      auto linearPart = f.linearPart();
 
-      auto origin = f.origin();
-      auto linearPart = f.originalLinearPart();
 
-      Dune::MatrixVector::addProduct(linearPart, maxEigenvalue, origin);
+      /*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;*/
 
-      for (size_t j = 0; j < ignore.size(); ++j)
-        if (ignore[j])
+      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;
+        }
+      }
+      std::cout << "lower: " << lower << " upper: " << upper << std::endl;
+
+      for (size_t j = 0; j < ignore.size(); ++j) {
+        if (ignore[j]) {
           linearPart[j] = 0; // makes sure result remains in subspace after correction
-        else
+        } else {
           origin[j] = 0; // shift affine offset
+        }
+      }
+
+      auto fv = directionalRestriction(f, x, linearPart);
 
-      double const linearPartNorm = linearPart.two_norm();
-      if (linearPartNorm <= 0.0)
-        return;
-      linearPart /= linearPartNorm;
+      //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:");
 
-      using Nonlinearity = std::decay_t<decltype(f.phi())>;
-      using Range = std::decay_t<decltype(f.maxEigenvalues())>;
-      using FirstOrderFunctional = FirstOrderFunctional<Nonlinearity, Vector, Range>;
-
-      auto lower = f.originalLowerObstacle();
-      auto upper = f.originalUpperObstacle();
 
       //print(lower, "lower:");
       //print(upper, "upper:");
 
-      FirstOrderFunctional firstOrderFunctional(maxEigenvalue, linearPartNorm, f.phi(),
+
+
+      /*FirstOrderFunctional firstOrderFunctional(maxEigenvalue, linearPartNorm, f.phi(),
                                                 lower, upper,
-                                                origin, linearPart);
+                                                origin, linearPart);*/
 
       //std::cout << "lower: " << firstOrderFunctional.lowerObstacle() << " upper: " << firstOrderFunctional.upperObstacle() << std::endl;
 
       // scalar obstacle solver without nonlinearity
-      typename Functional::Range alpha;
       //Dune::TNNMG::ScalarObstacleSolver obstacleSolver;
       //obstacleSolver(alpha, firstOrderFunctional, false);
 
@@ -89,16 +97,11 @@ class LocalBisectionSolver
 
       //std::cout << "domain: " << firstOrderFunctional.domain()[0] << " " << firstOrderFunctional.domain()[1] << std::endl;
 
+      typename Functional::Range alpha;
+      LineSearchSolver lineSearch;
+      lineSearch(alpha, fv, false);
 
-      const auto& domain = firstOrderFunctional.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(firstOrderFunctional, 0.0, 0.0, bisectionsteps);
-      } else {
-          alpha = domain[0];
-      }
+      //std::cout << "alpha: " << alpha << std::endl;
 
       linearPart *= alpha;
 
@@ -109,15 +112,16 @@ class LocalBisectionSolver
         x = origin;
         x += linearPart;
       } else {
-        x = f.origin();
+        x = origin;
       }
 #else
       if (std::abs(alpha)> safety) {
         x = origin;
         x += linearPart;
-        x -= f.origin();
+        //x -= f.origin();
       }
 #endif
+      //std::cout << x << std::endl << std::endl;
   }
 };