From 6b0114c7d4cabe9fae414ee3e2de99165d2195e4 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Mon, 15 Mar 2021 23:30:08 +0100
Subject: [PATCH] coordinateRestriction now yields ShiftedFunctional
---
.../spatial-solving/tnnmg/functional.hh | 20 ++++++++-
.../tnnmg/localbisectionsolver.hh | 41 ++++++++++---------
2 files changed, 39 insertions(+), 22 deletions(-)
diff --git a/dune/tectonic/spatial-solving/tnnmg/functional.hh b/dune/tectonic/spatial-solving/tnnmg/functional.hh
index 9d08eca2..2cb37df8 100755
--- a/dune/tectonic/spatial-solving/tnnmg/functional.hh
+++ b/dune/tectonic/spatial-solving/tnnmg/functional.hh
@@ -116,6 +116,22 @@ class ShiftedFunctional {
return phi_;
}
+ friend auto directionalRestriction(const ShiftedFunctional& f, const Vector& origin, const Vector& direction)
+ -> DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>
+ {
+ // rescale direction for numerical stability
+ auto dir = direction;
+ auto dirNorm = dir.two_norm();
+ if (dirNorm > 0.0)
+ dir /= dirNorm;
+ else {
+ dirNorm = 0.0;
+ dir = 0.0;
+ }
+
+ return DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>(f.quadraticPart(), f.originalLinearPart(), f.phi(), f.originalLowerObstacle(), f.originalUpperObstacle(), origin, dir, dirNorm);
+ }
+
protected:
Dune::Solvers::ConstCopyOrReference<M> quadraticPart_;
@@ -267,9 +283,9 @@ auto coordinateRestriction(const GlobalShiftedFunctional& f, const Index& i)
dli -= origini;
dui -= origini;
- return Functional<decltype(Aii), LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, Range>(Aii, std::move(ri), phii, std::move(dli), std::move(dui));
+ //return Functional<decltype(Aii), LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, Range>(Aii, std::move(ri), phii, std::move(dli), std::move(dui));
- //return ShiftedFunctional<decltype(Aii), LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, LocalVector, Range>(Aii, std::move(ri), phii, std::move(dli), std::move(dui), origini);
+ return ShiftedFunctional<decltype(Aii), LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, LocalVector, Range>(Aii, std::move(ri), phii, std::move(dli), std::move(dui), origini);
}
diff --git a/dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh b/dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh
index a9327ad7..ea694ea6 100644
--- a/dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh
+++ b/dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh
@@ -26,9 +26,8 @@ class LocalBisectionSolver
return;
}
- auto old_x = x;
- auto origin = x;
- auto linearPart = f.linearPart();
+ auto origin = f.origin();
+ auto linearPart = f.originalLinearPart();
//std::cout << x << std::endl;
@@ -36,11 +35,11 @@ class LocalBisectionSolver
//std::cout << "ignore: " << ignore << std::endl;
//std::cout << "orig linear part: " << linearPart << std::endl;
- auto lower = f.lowerObstacle();
- auto upper = f.upperObstacle();
+ auto lower = f.originalLowerObstacle();
+ auto upper = f.originalUpperObstacle();
for (size_t j = 0; j < lower.size(); ++j) {
if (std::abs(upper[j]-lower[j])<= safety) {
- //linearPart[j] = 0;
+ linearPart[j] = 0;
}
}
@@ -54,13 +53,14 @@ 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);
@@ -70,10 +70,10 @@ class LocalBisectionSolver
//std::cout << "origin: " << origin << std::endl;
//std::cout << " linearPart: " << linearPart << std::endl;
- std::cout << "maxEigenvalue " << f.quadraticPart() << std::endl;
- print(linearPart, "linearPart:");
- std::cout << "linearPartNorm " << linearPartNorm << std::endl;
- print(origin, "origin:");
+ //std::cout << "maxEigenvalue " << f.quadraticPart() << std::endl;
+ //print(linearPart, "linearPart:");
+ //std::cout << "linearPartNorm " << linearPartNorm << std::endl;
+ //print(origin, "origin:");
@@ -115,11 +115,11 @@ class LocalBisectionSolver
if (fv.scaling()<=0.0) {
alpha = 0.0;
- } /* else {
+ } else {
LineSearchSolver lineSearch;
lineSearch(alpha, fv, false);
- }*/
-
+ }
+/*
const auto& domain = fv.domain();
if (std::abs(domain[0]-domain[1])>safety) {
int bisectionsteps = 0;
@@ -128,15 +128,16 @@ class LocalBisectionSolver
alpha = globalBisection.minimize(fv, 0.0, 0.0, bisectionsteps);
} else {
alpha = domain[0];
- }
+ }*/
- std::cout << "alpha: " << alpha << std::endl;
+ //std::cout << "alpha: " << alpha << std::endl;
linearPart *= alpha;
- std::cout << linearPart << std::endl;
+ //std::cout << linearPart << std::endl;
+ //std::cout << "old_x: " << old_x << std::endl;
#ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
if (std::abs(alpha)> safety) {
x = origin;
@@ -145,14 +146,14 @@ class LocalBisectionSolver
x = origin;
}
#else
- std::cout << "use old" << std::endl;
+ //std::cout << "use old" << std::endl;
if (std::abs(alpha)> safety) {
x = origin;
x += linearPart;
- x -= old_x;//f.origin();
+ x -= f.origin();
}
#endif
- std::cout << "new x: "<< x << std::endl << std::endl;
+ //std::cout << "new x: "<< x << std::endl << std::endl;
}
};
--
GitLab