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