diff --git a/dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh b/dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh index 8df8265c5adf5fc8ee2ed7c5134f6b3524c1323a..c9f6bf9eb063d86cd1a20b707cb8df538510c7ba 100644 --- a/dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh +++ b/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; } };