From 19dcc1beace9f25d75bebecd6d3fe3cb2819660c Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Fri, 2 Aug 2019 10:38:07 +0200
Subject: [PATCH] .
---
src/CMakeLists.txt | 2 +
src/multi-body-problem-2D.cfg | 2 +-
src/multi-body-problem.cfg | 4 +-
src/solverfactorytest.cc | 96 +++++++++++--
src/spatial-solving/solverfactory.cc | 4 +-
src/spatial-solving/solverfactory.hh | 8 +-
src/spatial-solving/tnnmg/functional.hh | 126 +++++++++---------
.../tnnmg/localbisectionsolver.hh | 79 +++++++----
8 files changed, 209 insertions(+), 112 deletions(-)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index b9748ffc..6a92ecf9 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -106,6 +106,8 @@ foreach(_dim 2 3)
set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
set_property(TARGET ${_msw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+ #set_property(TARGET ${_msw_target} APPEND PROPERTY COMPILE_DEFINITIONS "NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY=1")
set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
set_property(TARGET ${_sft_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+ #set_property(TARGET ${_sft_target} APPEND PROPERTY COMPILE_DEFINITIONS "NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY=1")
endforeach()
diff --git a/src/multi-body-problem-2D.cfg b/src/multi-body-problem-2D.cfg
index 647e99cb..a308c80b 100644
--- a/src/multi-body-problem-2D.cfg
+++ b/src/multi-body-problem-2D.cfg
@@ -1,6 +1,6 @@
# -*- mode:conf -*-
[boundary.friction]
-smallestDiameter = 0.05 # 2e-3 [m]
+smallestDiameter = 0.01 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5 # 1e-5
diff --git a/src/multi-body-problem.cfg b/src/multi-body-problem.cfg
index b3a4ff9e..d58f80cf 100644
--- a/src/multi-body-problem.cfg
+++ b/src/multi-body-problem.cfg
@@ -47,10 +47,10 @@ relativeTau = 1e-4 # 1e-6
[timeSteps]
scheme = newmark
-timeSteps = 5
+timeSteps = 1
[u0.solver]
-maximumIterations = 100
+maximumIterations = 20
verbosity = full
[a0.solver]
diff --git a/src/solverfactorytest.cc b/src/solverfactorytest.cc
index 5365726f..3a4ceebe 100644
--- a/src/solverfactorytest.cc
+++ b/src/solverfactorytest.cc
@@ -74,6 +74,31 @@ size_t timeSteps = 0;
const std::string path = "";
const std::string outputFile = "solverfactorytest.log";
+std::vector<std::vector<double>> allReductionFactors;
+
+template<class IterationStepType, class NormType, class ReductionFactorContainer>
+Dune::Solvers::Criterion reductionFactorCriterion(
+ const IterationStepType& iterationStep,
+ const NormType& norm,
+ ReductionFactorContainer& reductionFactors)
+{
+ double normOfOldCorrection = 1;
+ auto lastIterate = std::make_shared<typename IterationStepType::Vector>(*iterationStep.getIterate());
+
+ return Dune::Solvers::Criterion(
+ [&, lastIterate, normOfOldCorrection] () mutable {
+ double normOfCorrection = norm.diff(*lastIterate, *iterationStep.getIterate());
+ double convRate = (normOfOldCorrection > 0) ? normOfCorrection / normOfOldCorrection : 0.0;
+
+ normOfOldCorrection = normOfCorrection;
+ *lastIterate = *iterationStep.getIterate();
+
+ reductionFactors.push_back(convRate);
+ return std::make_tuple(convRate < 0, Dune::formatString(" % '.5f", convRate));
+ },
+ " reductionFactor");
+}
+
template <class ContactNetwork>
void solveProblem(const ContactNetwork& contactNetwork,
const Matrix& mat, const Vector& rhs, Vector& x,
@@ -141,15 +166,15 @@ void solveProblem(const ContactNetwork& contactNetwork,
refSolver.addCriterion(
[&](){
- return Dune::formatString(" % 12.5e", I(x));
+ return Dune::formatString(" % 12.5e", I(refX));
},
" energy ");
- double initialEnergy = I(x);
+ double initialEnergy = I(refX);
refSolver.addCriterion(
[&](){
static double oldEnergy=initialEnergy;
- double currentEnergy = I(x);
+ double currentEnergy = I(refX);
double decrease = currentEnergy - oldEnergy;
oldEnergy = currentEnergy;
return Dune::formatString(" % 12.5e", decrease);
@@ -169,13 +194,13 @@ void solveProblem(const ContactNetwork& contactNetwork,
},
" truncated ");
-
- std::vector<double> correctionNorms;
- refSolver.addCriterion(Dune::Solvers::correctionNormCriterion(step, norm, 1e-10, correctionNorms));
+ if (timeStep>0 and initial) {
+ allReductionFactors[timeStep].resize(0);
+ refSolver.addCriterion(reductionFactorCriterion(step, norm, allReductionFactors[timeStep]));
+ }
refSolver.preprocess();
refSolver.solve();
- std::cout << correctionNorms.size() << std::endl;
/*if (initial) {
x = refX;
@@ -216,6 +241,35 @@ void solveProblem(const ContactNetwork& contactNetwork,
parset.get<double>("u0.solver.tolerance"), &norm,
Solver::FULL, true, &refX); // absolute error
+ solver.addCriterion(
+ [&](){
+ return Dune::formatString(" % 12.5e", J(x));
+ },
+ " energy ");
+
+ initialEnergy = J(x);
+ solver.addCriterion(
+ [&](){
+ static double oldEnergy=initialEnergy;
+ double currentEnergy = J(x);
+ double decrease = currentEnergy - oldEnergy;
+ oldEnergy = currentEnergy;
+ return Dune::formatString(" % 12.5e", decrease);
+ },
+ " decrease ");
+
+ solver.addCriterion(
+ [&](){
+ return Dune::formatString(" % 12.5e", step.lastDampingFactor());
+ },
+ " damping ");
+
+ solver.addCriterion(
+ [&](){
+ return Dune::formatString(" % 12d", tnnmgStep->linearization().truncated().count());
+ },
+ " truncated ");
+
solver.preprocess();
solver.solve();
@@ -479,7 +533,7 @@ void run(Updaters& updaters, ContactNetwork& contactNetwork,
// contribution from nonlinearity
globalFriction.updateAlpha(alpha);
- solveProblem(contactNetwork, bilinearForm, totalRhs, totalVelocityIterate, totalDirichletNodes, lower, upper);
+ solveProblem(contactNetwork, bilinearForm, totalRhs, totalVelocityIterate, totalDirichletNodes, lower, upper, false);
nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
@@ -619,6 +673,8 @@ int main(int argc, char *argv[]) { try {
// -----------------
// init input/output
// -----------------
+ timeSteps = parset.get<size_t>("timeSteps.timeSteps");
+ allReductionFactors.resize(timeSteps+1);
setupInitialConditions(contactNetwork);
@@ -739,7 +795,6 @@ int main(int argc, char *argv[]) { try {
//const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
- timeSteps = parset.get<size_t>("timeSteps.timeSteps");
for (timeStep=1; timeStep<=timeSteps; timeStep++) {
advanceStep(current, contactNetwork);
@@ -756,6 +811,29 @@ int main(int argc, char *argv[]) { try {
report();
}
+ // output reduction factors
+ size_t count = 0;
+ for (size_t i=0; i<allReductionFactors.size(); i++) {
+ count = std::max(count, allReductionFactors[i].size());
+ }
+ std::vector<double> avgReductionFactors(count);
+ for (size_t i=0; i<count; i++) {
+ avgReductionFactors[i] = 1;
+ size_t c = 0;
+
+ for (size_t j=0; j<allReductionFactors.size(); j++) {
+ if (!(i<allReductionFactors[j].size()))
+ continue;
+
+ avgReductionFactors[i] *= allReductionFactors[j][i];
+ c++;
+ }
+
+ avgReductionFactors[i] = std::pow(avgReductionFactors[i], 1.0/((double)c));
+ }
+
+ print(avgReductionFactors, "average reduction factors: ");
+
bool passed = true;
std::cout << "Overall, the test " << (passed ? "was successful!" : "failed!") << std::endl;
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index 7084672a..d682cd05 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -26,9 +26,9 @@ SolverFactory<Functional, BitVector>::SolverFactory(
auto linearSolver_ptr = Dune::Solvers::wrap_own_share<std::decay_t<LinearSolver>>(std::forward<LinearSolver>(linearSolver));
- //tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_ptr, 10, DefectProjection(), LineSearchSolver());
+ tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_ptr, DefectProjection(), Dune::TNNMG::ScalarObstacleSolver());
- tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_ptr, DefectProjection(), LineSearchSolver());
+ //tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_ptr, DefectProjection(), LineSearchSolver());
tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));
tnnmgStep_->setIgnore(ignoreNodes);
}
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index a68755ec..2e1855e0 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -1,7 +1,7 @@
#ifndef SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
#define SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
-#define NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
+//#define NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
@@ -28,13 +28,13 @@ class SolverFactory {
using Vector = typename Functional::Vector;
using BitVector = BitVectorType;
- using LocalSolver = LocalBisectionSolver;//Dune::TNNMG::ScalarObstacleSolver;//LocalBisectionSolver;
- using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, LocalBisectionSolver, BitVector>;//Dune::TNNMG::NonlinearGSStep<Functional, Dune::TNNMG::GaussSeidelLocalSolver<LocalSolver>, BitVector>;
+ using LocalSolver = LocalBisectionSolver; //Dune::TNNMG::ScalarObstacleSolver;//LocalBisectionSolver;
+ using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, LocalBisectionSolver, BitVector>; //Dune::TNNMG::NonlinearGSStep<Functional, Dune::TNNMG::GaussSeidelLocalSolver<LocalSolver>, BitVector>;
using Linearization = Linearization<Functional, BitVector>;
using DefectProjection = typename Dune::TNNMG::ObstacleDefectProjection;
//using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
- using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
+ using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, Dune::TNNMG::ScalarObstacleSolver>;
template <class LinearSolver>
SolverFactory(const Dune::ParameterTree&,
diff --git a/src/spatial-solving/tnnmg/functional.hh b/src/spatial-solving/tnnmg/functional.hh
index 8305b04b..5a1d3b37 100644
--- a/src/spatial-solving/tnnmg/functional.hh
+++ b/src/spatial-solving/tnnmg/functional.hh
@@ -16,6 +16,9 @@
#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
+template<class M, class V, class N, class L, class U, class R>
+class DirectionalRestriction;
+
template<class M, class V, class N, class R>
class ShiftedFunctional;
@@ -60,6 +63,12 @@ class ShiftedFunctional : public Dune::TNNMG::ShiftedBoxConstrainedQuadraticFunc
return phi_;
}
+ friend auto directionalRestriction(const ShiftedFunctional& f, const Vector& origin, const Vector& direction)
+ -> DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>
+ {
+ return DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>(f.quadraticPart(), f.originalLinearPart(), f.phi(), f.originalLowerObstacle(), f.originalUpperObstacle(), origin, direction);
+ }
+
protected:
const Nonlinearity& phi_;
};
@@ -148,6 +157,7 @@ class DirectionalRestriction :
Interval domain_;
};
+
/** \brief Box constrained quadratic functional with nonlinearity
* Note: call setIgnore() to set up functional in affine subspace with ignore information
*
@@ -157,7 +167,7 @@ class DirectionalRestriction :
* \tparam U Upper obstacle type
* \tparam R Range type
*/
-template<class V, class N, class L, class U, class O, class D, class R>
+/*template<class V, class N, class L, class U, class O, class D, class R>
class FirstOrderModelFunctional {
using Interval = typename Dune::Solvers::Interval<R>;
@@ -170,10 +180,30 @@ class FirstOrderModelFunctional {
using field_type = typename Vector::field_type;
-private:
- void init() {
- auto origin = origin_.get();
- auto direction = direction_.get();
+public:
+ template <class VV1, class NN, class LL, class UU, class VV2>
+ FirstOrderModelFunctional(
+ const Range& maxEigenvalue,
+ VV1&& linearTerm,
+ NN&& phi,
+ LL&& lower,
+ UU&& upper,
+ VV2&& origin) :
+ maxEigenvalue_(maxEigenvalue),
+ linearTerm_(std::forward<VV1>(linearTerm)),
+ lower_(std::forward<LL>(lower)),
+ upper_(std::forward<UU>(upper)),
+ origin_(std::forward<VV2>(origin)),
+ //direction_(linearTerm_),
+ phi_(std::forward<NN>(phi)) {
+
+ auto& origin = origin_.get();
+ auto& direction = direction_.get();
+
+ auto v = ri;
+ double const vnorm = v.two_norm();
+ if (vnorm > 1.0)
+ v /= vnorm;
// set defect obstacles
defectLower_ = -std::numeric_limits<field_type>::max();
@@ -190,61 +220,10 @@ class FirstOrderModelFunctional {
if (domain_[1] > defectUpper_) {
domain_[1] = defectUpper_;
}
-
- // set quadratic and linear parts of functional
- quadraticPart_ = direction*direction;
- quadraticPart_ *= maxEigenvalue_;
-
- linearPart_ = linearTerm_.get()*direction - maxEigenvalue_*(direction*origin);
- }
-
-public:
- template <class MM, class VV1, class NN, class LL, class UU, class VV2, class VV3>
- FirstOrderModelFunctional(
- const MM& matrix,
- VV1&& linearTerm,
- NN&& phi,
- LL&& lower,
- UU&& upper,
- VV2&& origin,
- VV3&& direction) :
- linearTerm_(std::forward<VV1>(linearTerm)),
- lower_(std::forward<LL>(lower)),
- upper_(std::forward<UU>(upper)),
- origin_(std::forward<VV2>(origin)),
- direction_(std::forward<VV3>(direction)),
- phi_(std::forward<NN>(phi)) {
-
- // set maxEigenvalue_ from matrix
- Vector eigenvalues;
- Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
- maxEigenvalue_ = *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
-
- init();
- }
-
- template <class VV1, class NN, class LL, class UU, class VV2, class VV3>
- FirstOrderModelFunctional(
- const Range& maxEigenvalue,
- VV1&& linearTerm,
- NN&& phi,
- LL&& lower,
- UU&& upper,
- VV2&& origin,
- VV3&& direction) :
- linearTerm_(std::forward<VV1>(linearTerm)),
- lower_(std::forward<LL>(lower)),
- upper_(std::forward<UU>(upper)),
- origin_(std::forward<VV2>(origin)),
- direction_(std::forward<VV3>(direction)),
- phi_(std::forward<NN>(phi)),
- maxEigenvalue_(maxEigenvalue) {
-
- init();
}
template <class BitVector>
- auto getIgnoreFunctional(const BitVector& ignore) const {
+ auto setIgnore(const BitVector& ignore) const {
Vector direction = direction_.get();
Vector origin = origin_.get();
// Dune::Solvers::resizeInitializeZero(direction, linearPart());
@@ -254,7 +233,9 @@ class FirstOrderModelFunctional {
else
origin[j] = 0; // shift affine offset
- return FirstOrderModelFunctional<Vector, Nonlinearity&, LowerObstacle, UpperObstacle, Vector, Vector, Range>(maxEigenvalue_, linearTerm_, phi_, lower_, upper_, std::move(origin), std::move(direction));
+ // set quadratic and linear parts of functional
+ quadraticPart_ = maxEigenvalue_ * direction.two_norm2();
+ linearPart_ = linearTerm_.get() * direction;
}
Range operator()(const Vector& v) const
@@ -315,7 +296,7 @@ class FirstOrderModelFunctional {
field_type defectLower_;
field_type defectUpper_;
-};
+};*/
// \ToDo This should be an inline friend of ShiftedBoxConstrainedQuadraticFunctional
// but gcc-4.9.2 shipped with debian jessie does not accept
@@ -345,6 +326,18 @@ auto coordinateRestriction(const ShiftedFunctional<M, V, Nonlinearity, R>& f, co
Dune::TNNMG::Imp::mmv(Aij, f.origin()[j], ri, Dune::PriorityTag<1>());
});
+ // set maxEigenvalue from matrix
+ LocalVector eigenvalues;
+ Dune::FMatrixHelp::eigenValues(*Aii_p, eigenvalues);
+ auto maxEigenvalue = *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
+
+ Dune::MatrixVector::addProduct(ri, maxEigenvalue, f.origin()[i]);
+
+ LocalMatrix Aii = 0;
+ for (size_t d=0; d<Aii.N(); d++) {
+ Aii[d][d] = maxEigenvalue;
+ }
+
LocalLowerObstacle dli = f.originalLowerObstacle()[i];
LocalUpperObstacle dui = f.originalUpperObstacle()[i];
dli -= f.origin()[i];
@@ -352,13 +345,16 @@ auto coordinateRestriction(const ShiftedFunctional<M, V, Nonlinearity, R>& f, co
auto&& phii = f.phi().restriction(i);
- auto v = ri;
- double const vnorm = v.two_norm();
- if (vnorm > 1.0)
- v /= vnorm;
+ /*std::cout << "maxEigenvalue matrix: " << std::endl;
+ for (size_t i=0; i<Aii.N(); i++) {
+ for (size_t j=0; j<Aii.M(); j++) {
+ std::cout << Aii[i][j] << " ";
+ }
+ std::cout << std::endl;
+ }*/
- return FirstOrderModelFunctional<LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, LocalVector, LocalVector, Range>(*Aii_p, std::move(ri), std::move(phii), std::move(dli), std::move(dui), std::move(f.origin()[i]), std::move(v));
- //return Functional<LocalMatrix&, LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, Range>(*Aii_p, std::move(ri), std::move(phii), std::move(dli), std::move(dui));
+ //return FirstOrderModelFunctional<LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, LocalVector, LocalVector, Range>(maxEigenvalue, std::move(ri), std::move(phii), std::move(dli), std::move(dui), std::move(f.origin()[i]));
+ return ShiftedFunctional<LocalMatrix, LocalVector, decltype(phii), Range>(std::move(Aii), std::move(ri), std::move(phii), std::move(dli), std::move(dui), f.origin()[i]);
}
diff --git a/src/spatial-solving/tnnmg/localbisectionsolver.hh b/src/spatial-solving/tnnmg/localbisectionsolver.hh
index 50618dd6..33f69e95 100644
--- a/src/spatial-solving/tnnmg/localbisectionsolver.hh
+++ b/src/spatial-solving/tnnmg/localbisectionsolver.hh
@@ -3,7 +3,9 @@
#ifndef DUNE_TECTONIC_LOCALBISECTIONSOLVER_HH
#define DUNE_TECTONIC_LOCALBISECTIONSOLVER_HH
+#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>
#include <dune/tnnmg/problem-classes/bisection.hh>
+#include "../../utils/debugutils.hh"
/**
* \brief A local solver for quadratic obstacle problems with nonlinearity
@@ -14,48 +16,67 @@ class LocalBisectionSolver
public:
template<class Vector, class Functional, class BitVector>
void operator()(Vector& x, const Functional& f, const BitVector& ignore) const {
+ double safety = 1e-14; //or (f.upperObstacle()-f.lowerObstacle()<safety)
if (ignore.all()) {
return;
}
- double safety = 1e-14;
- if (f.defectUpper() - f.defectLower() < safety) {
- x = f.origin();
- return;
- }
+ Vector origin = f.origin();
+ Vector direction = f.originalLinearPart();
- int bisectionsteps = 0;
- const Bisection globalBisection;
+ for (size_t j = 0; j < ignore.size(); ++j)
+ if (ignore[j])
+ direction[j] = 0; // makes sure result remains in subspace after correction
+ else
+ origin[j] = 0; // shift affine offset
- if (ignore.any()) {
- auto&& restrictedf = f.getIgnoreFunctional(ignore);
- const auto alpha = globalBisection.minimize(restrictedf, 0.0, 0.0, bisectionsteps);
+ double const directionNorm = direction.two_norm();
+ if (directionNorm <= 0.0)
+ return;
+ direction /= directionNorm;
-#ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
- x = restrictedf.origin();
- auto correction = restrictedf.direction();
- correction *= alpha;
- x += correction;
-#else
- auto x = restrictedf.direction();
- x *= alpha;
-#endif
+ auto&& directionalF = directionalRestriction(f, origin, direction);
+
+ // scalar obstacle solver without nonlinearity
+ typename Functional::Range alpha;
+ Dune::TNNMG::ScalarObstacleSolver obstacleSolver;
+ obstacleSolver(alpha, directionalF, false);
+ //direction *= alpha;
- } else { // ignore.none()
- const auto alpha = globalBisection.minimize(f, 0.0, 0.0, bisectionsteps);
+ /*const auto& A = f.quadraticPart();
+ std::cout << "f.quadratic(): " << std::endl;
+ for (size_t i=0; i<A.N(); i++) {
+ for (size_t j=0; j<A.M(); j++) {
+ std::cout << A[i][j] << " ";
+ }
+ std::cout << std::endl;
+ }*/
+ //std::cout << f.quadraticPart() << std::endl;
+ //std::cout << "A: " << directionalF.quadraticPart() << " " << (directionalF.quadraticPart()==f.quadraticPart()[0][0]) << std::endl;
+ /*std::cout << "b: " << directionalF.linearPart() << std::endl;
+ std::cout << "domain: " << directionalF.domain()[0] << " " << directionalF.domain()[1] << std::endl;
+ auto D = directionalF.subDifferential(0);
+ std::cout << "subDiff at x: " << D[0] << " " << D[1] << std::endl;*/
+
+ const auto& domain = directionalF.domain();
+ if (std::abs(domain[0]-domain[1])>safety) {
+ int bisectionsteps = 0;
+ const Bisection globalBisection;
+
+ alpha = globalBisection.minimize(directionalF, alpha, 0.0, bisectionsteps);
+ }
+
+ direction *= alpha;
#ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
- x = f.origin();
- auto correction = f.direction();
- correction *= alpha;
- x += correction;
+ x = origin;
+ x += direction;
#else
- auto x = f.direction();
- x *= alpha;
+ //x = origin;
+ x -= f.origin();
+ x += direction;
#endif
-
- }
}
};
--
GitLab