diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh
index 3bf6a11debb3173bbd97c5ff53e0795a2dff76ec..60b474ed767bb14fbbdbbd18565fddd816195db7 100644
--- a/dune/tectonic/frictionpotential.hh
+++ b/dune/tectonic/frictionpotential.hh
@@ -44,7 +44,7 @@ class TruncatedRateState : public FrictionPotential {
return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
}
- double second_deriv(double V) const override {
+ double second_deriv(double V) const override {
if (V <= Vmin)
return 0;
diff --git a/dune/tectonic/localfriction.hh b/dune/tectonic/localfriction.hh
index 227b8ff5c3af92f90a9c1c83d6bad6a163c9a079..12b98e5606f3424e477929fead99d77f66ab31ee 100644
--- a/dune/tectonic/localfriction.hh
+++ b/dune/tectonic/localfriction.hh
@@ -53,7 +53,7 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
double regularity(VectorType const &x) const override {
double const xnorm = x.two_norm();
- if (xnorm <= 0.0)
+ if (xnorm < 0.0)
return std::numeric_limits<double>::infinity();
return func_.regularity(xnorm);
diff --git a/src/multi-body-problem.cfg b/src/multi-body-problem.cfg
index 77de9404eafc006ed026f47ebfd43b58a4db1fe1..83476cb76ec11d3bd5d1e9f0375e731d578692ba 100644
--- a/src/multi-body-problem.cfg
+++ b/src/multi-body-problem.cfg
@@ -32,7 +32,7 @@ V0 = 5e-5 # [m/s]
L = 2.25e-5 # [m]
initialAlpha = 0 # [ ]
stateModel = AgeingLaw
-frictionModel = Regularised
+frictionModel = Truncated #Regularised
[boundary.friction.weakening]
a = 0.002 # [ ]
b = 0.017 # [ ]
@@ -47,10 +47,10 @@ relativeTau = 1e-4 # 1e-6
[timeSteps]
scheme = newmark
-timeSteps = 1
+timeSteps = 100
[u0.solver]
-maximumIterations = 20
+maximumIterations = 100
verbosity = full
[a0.solver]
diff --git a/src/solverfactorytest.cc b/src/solverfactorytest.cc
index 3a4ceebe7d28710f6197627839ffefe651442271..a2c95ae0e3596cbe4fe6f128f870bdc4c872ab6a 100644
--- a/src/solverfactorytest.cc
+++ b/src/solverfactorytest.cc
@@ -90,6 +90,9 @@ Dune::Solvers::Criterion reductionFactorCriterion(
double normOfCorrection = norm.diff(*lastIterate, *iterationStep.getIterate());
double convRate = (normOfOldCorrection > 0) ? normOfCorrection / normOfOldCorrection : 0.0;
+ if (convRate>1.0)
+ DUNE_THROW(Dune::Exception, "Solver convergence rate of " + std::to_string(convRate));
+
normOfOldCorrection = normOfCorrection;
*lastIterate = *iterationStep.getIterate();
@@ -202,18 +205,30 @@ void solveProblem(const ContactNetwork& contactNetwork,
refSolver.preprocess();
refSolver.solve();
- /*if (initial) {
+ //print(refX, "refX: ");
+
+ if (initial) {
x = refX;
return;
- }*/
+ }
// set up solver factory solver
// set up functional
- /*auto& globalFriction = contactNetwork.globalFriction();
+ auto& globalFriction = contactNetwork.globalFriction();
+
+ /* std::vector<const Dune::BitSetVector<1>*> nodes;
+ contactNetwork.frictionNodes(nodes);
+ print(*nodes[0], "frictionNodes: ");
+ print(*nodes[1], "frictionNodes: ");
+
+ print(ignore, "ignore: ");*/
+
using MyFunctional = Functional<Matrix&, Vector&, std::decay_t<decltype(globalFriction)>&, Vector&, Vector&, typename Matrix::field_type>;
- MyFunctional J(mat, rhs, globalFriction, lower, upper);*/
- using MyFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
- MyFunctional J(mat, rhs, ZeroNonlinearity(), lower, upper);
+ MyFunctional J(mat, rhs, globalFriction, lower, upper);
+ //using MyFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity, Vector&, Vector&, typename Matrix::field_type>;
+ //MyFunctional J(mat, rhs, ZeroNonlinearity(), lower, upper);
+
+ //std::cout << "ref energy: " << J(refX) << std::endl;
// set up TNMMG solver
// dummy solver, uses direct solver for linear correction no matter what is set here
@@ -239,7 +254,7 @@ void solveProblem(const ContactNetwork& contactNetwork,
LoopSolver<Vector> solver(
tnnmgStep.get(), parset.get<size_t>("u0.solver.maximumIterations"),
parset.get<double>("u0.solver.tolerance"), &norm,
- Solver::FULL, true, &refX); // absolute error
+ Solver::FULL); //, true, &refX); // absolute error
solver.addCriterion(
[&](){
@@ -260,7 +275,7 @@ void solveProblem(const ContactNetwork& contactNetwork,
solver.addCriterion(
[&](){
- return Dune::formatString(" % 12.5e", step.lastDampingFactor());
+ return Dune::formatString(" % 12.5e", tnnmgStep->lastDampingFactor());
},
" damping ");
@@ -270,6 +285,10 @@ void solveProblem(const ContactNetwork& contactNetwork,
},
" truncated ");
+
+ std::vector<double> factors;
+ solver.addCriterion(reductionFactorCriterion(*tnnmgStep, norm, factors));
+
solver.preprocess();
solver.solve();
@@ -527,10 +546,26 @@ void run(Updaters& updaters, ContactNetwork& contactNetwork,
Vector totalVelocityIterate;
nBodyAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
+ // project in onto admissible set
+ const size_t blocksize = Vector::block_type::dimension;
+ for (size_t i=0; i<totalVelocityIterate.size(); i++) {
+ for (size_t j=0; j<blocksize; j++) {
+ if (totalVelocityIterate[i][j] < lower[i][j]) {
+ totalVelocityIterate[i][j] = lower[i][j];
+ }
+
+ if (totalVelocityIterate[i][j] > upper[i][j]) {
+ totalVelocityIterate[i][j] = upper[i][j];
+ }
+ }
+ }
+
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations;
++fixedPointIteration) {
// contribution from nonlinearity
+ //print(alpha, "alpha: ");
+
globalFriction.updateAlpha(alpha);
solveProblem(contactNetwork, bilinearForm, totalRhs, totalVelocityIterate, totalDirichletNodes, lower, upper, false);
@@ -549,6 +584,8 @@ void run(Updaters& updaters, ContactNetwork& contactNetwork,
std::vector<Vector> v_rel;
relativeVelocities(contactNetwork, totalVelocityIterate, v_rel);
+ //print(v_rel, "v_rel");
+
std::cout << "- State problem set" << std::endl;
// solve a state problem
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index d682cd0594f4065759e0ed3960a4ce2f9f84e289..bca8290aab62771a1ff759ca6b09256e8b6aba4d 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, DefectProjection(), Dune::TNNMG::ScalarObstacleSolver());
+ //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 2e1855e05ad051cec6e0315db54d2815adc543e8..1613dd45ba77e1984d69b405d9445eb935982145 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -33,8 +33,8 @@ class SolverFactory {
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, Dune::TNNMG::ScalarObstacleSolver>;
+ 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
old mode 100644
new mode 100755
index 97dab68876812aefab43d4fe30549bdecf20429e..eff8e22d34a8572122b67b503ed573463117d2d5
--- a/src/spatial-solving/tnnmg/functional.hh
+++ b/src/spatial-solving/tnnmg/functional.hh
@@ -21,7 +21,7 @@
template<class M, class V, class N, class L, class U, class R>
class DirectionalRestriction;
-template<class M, class E, class V, class N, class R>
+template<class M, class E, class V, class N, class L, class U, class O, class R>
class ShiftedFunctional;
template <class M, class V, class N, class L, class U, class R>
@@ -35,32 +35,82 @@ class Functional;
* \tparam V nonlinearity type
* \tparam R Range type
*/
-template<class M, class E, class V, class N, class R>
-class ShiftedFunctional : public Dune::TNNMG::ShiftedBoxConstrainedQuadraticFunctional<M,V,R> {
- using Base = Dune::TNNMG::ShiftedBoxConstrainedQuadraticFunctional<M,V,R>;
+template<class M, class E, class V, class N, class L, class U, class O, class R>
+class ShiftedFunctional {
public:
- using Matrix = M;
+ using Matrix = std::decay_t<M>;
using Eigenvalues = E;
- using Vector = V;
- using Nonlinearity = N;
- using LowerObstacle = Vector;
- using UpperObstacle = Vector;
+ using Vector = std::decay_t<V>;
+ using Nonlinearity = std::decay_t<N>;
+ using LowerObstacle = std::decay_t<L>;
+ using UpperObstacle = std::decay_t<U>;
+ using Origin = std::decay_t<O>;
using Range = R;
- ShiftedFunctional(const Matrix& quadraticPart, const Eigenvalues& maxEigenvalues, const Vector& linearPart, const Nonlinearity& phi, const LowerObstacle& lowerObstacle, const UpperObstacle& upperObstacle, const Vector& origin) :
- Base(quadraticPart, linearPart, lowerObstacle, upperObstacle, origin),
- phi_(phi),
- maxEigenvalues_(maxEigenvalues)
+ template <class MM, class VV, class LL, class UU, class OO>
+ ShiftedFunctional(MM&& matrix, const Eigenvalues& maxEigenvalues, VV&& linearPart, const Nonlinearity& phi,
+ LL&& lower, UU&& upper, OO&& origin) :
+ quadraticPart_(std::forward<MM>(matrix)),
+ maxEigenvalues_(maxEigenvalues),
+ originalLinearPart_(std::forward<VV>(linearPart)),
+ phi_(phi),
+ originalLowerObstacle_(std::forward<LL>(lower)),
+ originalUpperObstacle_(std::forward<UU>(upper)),
+ origin_(std::forward<OO>(origin))
{}
Range operator()(const Vector& v) const
{
- auto temp = Base::Base::origin();
- temp += v;
- if (checkInsideBox(temp, this->originalLowerObstacle(), this->originalUpperObstacle()))
- return Base::Base::operator()(v) + phi_.operator()(temp);
+ auto w = origin();
+ w += v;
+ if (Dune::TNNMG::checkInsideBox(w, originalLowerObstacle(), originalUpperObstacle())) {
+ Vector temp;
+ Dune::Solvers::resizeInitializeZero(temp,v);
+ quadraticPart().umv(w, temp);
+ temp *= 0.5;
+ temp -= originalLinearPart();
+ return temp*w + phi()(w);
+ }
return std::numeric_limits<Range>::max();
+
+ /*Vector temp;
+ Dune::Solvers::resizeInitializeZero(temp,v);
+ quadraticPart().umv(w, temp);
+ temp *= 0.5;
+ temp -= originalLinearPart();
+ return temp*w + phi()(w);*/
+ }
+
+ const Matrix& quadraticPart() const
+ {
+ return quadraticPart_.get();
+ }
+
+ const Vector& originalLinearPart() const
+ {
+ return originalLinearPart_.get();
+ }
+
+ const Origin& origin() const
+ {
+ return origin_.get();
+ }
+
+ void updateOrigin()
+ {}
+
+ void updateOrigin(std::size_t i)
+ {}
+
+ const LowerObstacle& originalLowerObstacle() const
+ {
+ return originalLowerObstacle_.get();
+ }
+
+ const UpperObstacle& originalUpperObstacle() const
+ {
+ return originalUpperObstacle_.get();
}
const auto& phi() const {
@@ -72,13 +122,17 @@ class ShiftedFunctional : public Dune::TNNMG::ShiftedBoxConstrainedQuadraticFunc
}
protected:
- const Nonlinearity& phi_;
+ Dune::Solvers::ConstCopyOrReference<M> quadraticPart_;
const Eigenvalues& maxEigenvalues_;
-};
-
+ Dune::Solvers::ConstCopyOrReference<V> originalLinearPart_;
+ const Nonlinearity& phi_;
+ Dune::Solvers::ConstCopyOrReference<L> originalLowerObstacle_;
+ Dune::Solvers::ConstCopyOrReference<U> originalUpperObstacle_;
+ Dune::Solvers::ConstCopyOrReference<O> origin_;
+};
/** \brief Coordinate restriction of box constrained quadratic functional with nonlinearity;
* mainly used for line search step in TNNMG algorithm
@@ -183,11 +237,11 @@ class FirstOrderFunctional {
using field_type = typename V::field_type;
public:
- template <class NN, class LL, class UU, class OO, class DD>
+ template <class LL, class UU, class OO, class DD>
FirstOrderFunctional(
const Range& maxEigenvalue,
const Range& linearPart,
- NN&& phi,
+ const Nonlinearity& phi,
LL&& lower,
UU&& upper,
OO&& origin,
@@ -198,7 +252,7 @@ class FirstOrderFunctional {
upper_(std::forward<UU>(upper)),
origin_(std::forward<OO>(origin)),
direction_(std::forward<DD>(direction)),
- phi_(std::forward<NN>(phi)) {
+ phi_(phi) {
// set defect obstacles
defectLower_ = -std::numeric_limits<field_type>::max();
@@ -206,15 +260,15 @@ class FirstOrderFunctional {
Dune::TNNMG::directionalObstacles(origin_.get(), direction_.get(), lower_.get(), upper_.get(), defectLower_, defectUpper_);
// set domain
- phi_.get().directionalDomain(origin_.get(), direction_.get(), domain_);
+ phi_.directionalDomain(origin_.get(), direction_.get(), domain_);
- /*if (domain_[0] < defectLower_) {
+ if (domain_[0] < defectLower_) {
domain_[0] = defectLower_;
}
if (domain_[1] > defectUpper_) {
domain_[1] = defectUpper_;
- }*/
+ }
}
Range operator()(const Vector& v) const
@@ -227,7 +281,7 @@ class FirstOrderFunctional {
Vector uxv = origin_.get();
Dune::MatrixVector::addProduct(uxv, x, direction_.get());
- phi_.get().directionalSubDiff(uxv, direction_.get(), Di);
+ phi_.directionalSubDiff(uxv, direction_.get(), Di);
const auto Axmb = quadraticPart_ * x - linearPart_;
Di[0] += Axmb;
@@ -273,7 +327,7 @@ class FirstOrderFunctional {
Dune::Solvers::ConstCopyOrReference<Vector> origin_;
Dune::Solvers::ConstCopyOrReference<Vector> direction_;
- Dune::Solvers::ConstCopyOrReference<N> phi_;
+ const Nonlinearity& phi_;
Interval domain_;
@@ -285,10 +339,10 @@ class FirstOrderFunctional {
// but gcc-4.9.2 shipped with debian jessie does not accept
// inline friends with auto return type due to bug-59766.
// Notice, that this is fixed in gcc-4.9.3.
-template<class Index, class M, class E, class V, class Nonlinearity, class R>
-auto coordinateRestriction(const ShiftedFunctional<M, E, V, Nonlinearity, R>& f, const Index& i)
+template<class GlobalShiftedFunctional, class Index>
+auto coordinateRestriction(const GlobalShiftedFunctional& f, const Index& i)
{
- using Range = R;
+ using Range = typename GlobalShiftedFunctional::Range;
using LocalMatrix = std::decay_t<decltype(f.quadraticPart()[i][i])>;
using LocalVector = std::decay_t<decltype(f.originalLinearPart()[i])>;
using LocalLowerObstacle = std::decay_t<decltype(f.originalLowerObstacle()[i])>;
@@ -299,7 +353,8 @@ auto coordinateRestriction(const ShiftedFunctional<M, E, V, Nonlinearity, R>& f,
const LocalMatrix* Aii_p = nullptr;
- //print(f.origin(), "f.origin: ");
+ //print(f.originalLinearPart(), "f.linearPart: ");
+ //print(f.quadraticPart(), "f.quadraticPart: ");
LocalVector ri = f.originalLinearPart()[i];
const auto& Ai = f.quadraticPart()[i];
@@ -311,9 +366,15 @@ auto coordinateRestriction(const ShiftedFunctional<M, E, V, Nonlinearity, R>& f,
Dune::TNNMG::Imp::mmv(Aij, f.origin()[j], ri, Dune::PriorityTag<1>());
});
- auto&& phii = f.phi().restriction(i);
+ //print(*Aii_p, "Aii_p:");
+ //print(ri, "ri:");
+
+ //print(f.originalLowerObstacle()[i], "lower:");
+ //print(f.originalUpperObstacle()[i], "upper:");
+
+ auto& phii = f.phi().restriction(i);
- return ShiftedFunctional<LocalMatrix, Range, LocalVector, decltype(phii), Range>(*Aii_p, f.maxEigenvalues()[i], std::move(ri), std::move(phii), f.originalLowerObstacle()[i], f.originalUpperObstacle()[i], f.origin()[i]);
+ return ShiftedFunctional<LocalMatrix&, Range, LocalVector, std::decay_t<decltype(phii)>, LocalLowerObstacle&, LocalUpperObstacle&, LocalVector&, Range>(*Aii_p, f.maxEigenvalues()[i], std::move(ri), phii, f.originalLowerObstacle()[i], f.originalUpperObstacle()[i], f.origin()[i]);
}
@@ -378,8 +439,12 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
Range operator()(const Vector& v) const
{
- if (Dune::TNNMG::checkInsideBox(v, this->lower_.get(), this->upper_.get()))
- return Dune::TNNMG::QuadraticFunctional<M,V,R>::operator()(v) + phi_.get().operator()(v);
+ //print(v, "v:");
+ //print(this->lowerObstacle(), "lower: ");
+ //print(this->upperObstacle(), "upper: ");
+ //std::cout << Dune::TNNMG::checkInsideBox(v, this->lowerObstacle(), this->upperObstacle()) << " " << Dune::TNNMG::QuadraticFunctional<M,V,R>::operator()(v) << " " << phi_.get().operator()(v) << std::endl;
+ if (Dune::TNNMG::checkInsideBox(v, this->lowerObstacle(), this->upperObstacle()))
+ return Dune::TNNMG::QuadraticFunctional<M,V,R>::operator()(v) + phi_.get()(v);
return std::numeric_limits<Range>::max();
}
@@ -390,9 +455,9 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
}
friend auto shift(const Functional& f, const Vector& origin)
- -> ShiftedFunctional<Matrix, Eigenvalues, Vector, Nonlinearity, Range>
+ -> ShiftedFunctional<Matrix&, Eigenvalues, Vector&, Nonlinearity&, LowerObstacle&, UpperObstacle&, Vector&, Range>
{
- return ShiftedFunctional<Matrix, Eigenvalues, Vector, Nonlinearity, Range>(f.quadraticPart(), f.maxEigenvalues(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin);
+ return ShiftedFunctional<Matrix&, Eigenvalues, Vector&, Nonlinearity&, LowerObstacle&, UpperObstacle&, Vector&, Range>(f.quadraticPart(), f.maxEigenvalues(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin);
}
};
diff --git a/src/spatial-solving/tnnmg/linearization.hh b/src/spatial-solving/tnnmg/linearization.hh
index bc6095a8cb03395cfa236332a269f419962ea203..362dcec47a4daa8a6b43432339368980ed3b3586 100644
--- a/src/spatial-solving/tnnmg/linearization.hh
+++ b/src/spatial-solving/tnnmg/linearization.hh
@@ -33,7 +33,14 @@ class Linearization {
template <class T>
void determineRegularityTruncation(const Vector& x, BitVector& truncationFlags, const T& truncationTolerance)
{
+ //std::cout << "x: " << x << std::endl;
+
for (size_t i = 0; i < x.size(); ++i) {
+ //if (truncationFlags[i].all())
+ // continue;
+
+ //std::cout << f_.phi().regularity(i, x[i]) << " xnorm: " << x[i].two_norm() << std::endl;
+
if (f_.phi().regularity(i, x[i]) > truncationTolerance) {
truncationFlags[i] = true;
}
@@ -90,7 +97,8 @@ class Linearization {
Linearization(const F& f, const BitVector& ignore) :
f_(f),
ignore_(ignore),
- truncationTolerance_(1e-10)
+ obstacleTruncationTolerance_(1e-10),
+ regularityTruncationTolerance_(1e8)
{}
void bind(const Vector& x) {
@@ -101,11 +109,28 @@ class Linearization {
// determine which components to truncate
// ----------------
- truncationFlags_ = ignore_;
+ BitVector obstacleTruncationFlags = ignore_;
+ BitVector regularityTruncationFlags = ignore_;
+
+ //std::cout << "ignore truncation: " << truncationFlags_.count();
+
+ determineTruncation(x, f_.lowerObstacle(), f_.upperObstacle(), obstacleTruncationFlags, obstacleTruncationTolerance_); // obstacle truncation
+
+ //std::cout << " obstacle truncation: " << truncationFlags_.count();
+ //std::cout << " " << obstacleTruncationTolerance_;
- determineTruncation(x, f_.lowerObstacle(), f_.upperObstacle(), truncationFlags_, truncationTolerance_); // obstacle truncation
- determineRegularityTruncation(x, truncationFlags_, truncationTolerance_); // truncation due to regularity deficit of nonlinearity
+ determineRegularityTruncation(x, regularityTruncationFlags, regularityTruncationTolerance_); // truncation due to regularity deficit of nonlinearity
+ //std::cout << " regularity truncation: " << truncationFlags_.count() << std::endl;
+ //std::cout << " " << regularityTruncationTolerance_;
+
+ truncationFlags_ = obstacleTruncationFlags;
+ size_t blocksize = truncationFlags_[0].size();
+ for (size_t i=0; i<truncationFlags_.size(); i++) {
+ for (size_t j=0; j<blocksize; j++) {
+ //truncationFlags_[i][j] = truncationFlags_[i][j] and regularityTruncationFlags[i][j];
+ }
+ }
// compute hessian
// ---------------
@@ -126,11 +151,15 @@ class Linearization {
hessian_[i][it.index()] += *it;
}
+ truncateMatrix(hessian_, obstacleTruncationFlags, obstacleTruncationFlags);
+
//print(hessian_, "Hessian:");
//std::cout << "--------" << (double) hessian_[21][21] << "------------" << std::endl;
// nonlinearity part
phi.addHessian(x, hessian_);
+ truncateMatrix(hessian_, regularityTruncationFlags, regularityTruncationFlags);
+
// compute gradient
// ----------------
@@ -144,9 +173,11 @@ class Linearization {
// -grad is needed for Newton step
negativeGradient_ *= -1;
- // truncate gradient and hessian
+ // truncate gradient
truncateVector(negativeGradient_, truncationFlags_);
- truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
+
+ //print(hessian_, "hessian: ");
+ //print(negativeGradient_, "negativeGradient: ");
}
@@ -171,7 +202,8 @@ class Linearization {
const F& f_;
const BitVector& ignore_;
- double truncationTolerance_;
+ double obstacleTruncationTolerance_;
+ double regularityTruncationTolerance_;
Vector negativeGradient_;
Matrix hessian_;
diff --git a/src/spatial-solving/tnnmg/linesearchsolver.hh b/src/spatial-solving/tnnmg/linesearchsolver.hh
index 93471396dc10a5e89299f59e5d3a6928ddd719fb..587c81567d57ba7f49e90243018ba692b18c02e2 100644
--- a/src/spatial-solving/tnnmg/linesearchsolver.hh
+++ b/src/spatial-solving/tnnmg/linesearchsolver.hh
@@ -18,28 +18,24 @@ class LineSearchSolver
void operator()(Vector& x, const Functional& f, const BitVector& ignore) const {
x = 0;
- Dune::Solvers::Interval<double> D;
+ if (ignore)
+ return;
+
+ /*Dune::Solvers::Interval<double> D;
D = f.subDifferential(0);
if (D[1] > 0) // NOTE: Numerical instability can actually get us here
return;
-
+*/
if (almost_equal(f.domain()[0], f.domain()[1], 2)) {
x = f.domain()[0];
return;
}
- if (not ignore) {
- int bisectionsteps = 0;
- const Bisection globalBisection;
- // std::cout << "LineSearchSolver: starting bisection" << std::endl;
- x = globalBisection.minimize(f, 1.0, 0.0, bisectionsteps);
+ int bisectionsteps = 0;
+ const Bisection globalBisection(0.0, 1.0, 0.0, 0.0);;
- x = f.domain().projectIn(x);
-
-
- //std::cout << "LineSearchSolver: bisection terminated" << std::endl;
- // x = globalBisection.minimize(f, vNorm, 0.0, bisectionsteps) / vNorm; // used rescaled v in f?
- }
+ x = globalBisection.minimize(f, 0.0, 0.0, bisectionsteps);
+ //x = f.domain().projectIn(x);
}
};
diff --git a/src/spatial-solving/tnnmg/localbisectionsolver.hh b/src/spatial-solving/tnnmg/localbisectionsolver.hh
index b24eb8ee6f84da9652d11924ea1a44056a1836c8..8df8265c5adf5fc8ee2ed7c5134f6b3524c1323a 100644
--- a/src/spatial-solving/tnnmg/localbisectionsolver.hh
+++ b/src/spatial-solving/tnnmg/localbisectionsolver.hh
@@ -30,6 +30,7 @@ class LocalBisectionSolver
auto origin = f.origin();
auto linearPart = f.originalLinearPart();
+
Dune::MatrixVector::addProduct(linearPart, maxEigenvalue, origin);
for (size_t j = 0; j < ignore.size(); ++j)
@@ -43,16 +44,31 @@ class LocalBisectionSolver
return;
linearPart /= linearPartNorm;
+ //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>;
- FirstOrderFunctional firstOrderFunctional(maxEigenvalue, linearPartNorm, f.phi(), f.originalLowerObstacle(), f.originalUpperObstacle(), origin, linearPart);
+ auto lower = f.originalLowerObstacle();
+ auto upper = f.originalUpperObstacle();
+
+ //print(lower, "lower:");
+ //print(upper, "upper:");
+
+ FirstOrderFunctional firstOrderFunctional(maxEigenvalue, linearPartNorm, f.phi(),
+ lower, upper,
+ 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);*/
+ //Dune::TNNMG::ScalarObstacleSolver obstacleSolver;
+ //obstacleSolver(alpha, firstOrderFunctional, false);
//direction *= alpha;
@@ -71,31 +87,36 @@ class LocalBisectionSolver
auto D = directionalF.subDifferential(0);
std::cout << "subDiff at x: " << D[0] << " " << D[1] << std::endl;*/
+ //std::cout << "domain: " << firstOrderFunctional.domain()[0] << " " << firstOrderFunctional.domain()[1] << std::endl;
+
+
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);
-
- if (alpha < firstOrderFunctional.lowerObstacle() ) {
- alpha = firstOrderFunctional.lowerObstacle();
- }
-
- if (alpha > firstOrderFunctional.upperObstacle() ) {
- alpha = firstOrderFunctional.upperObstacle();
- }
+ } else {
+ alpha = domain[0];
}
linearPart *= alpha;
+ //std::cout << linearPart << std::endl;
+
#ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
- x = origin;
- x += linearPart;
+ if (std::abs(alpha)> safety) {
+ x = origin;
+ x += linearPart;
+ } else {
+ x = f.origin();
+ }
#else
- //x = origin;
- //x -= f.origin();
- x = linearPart;
+ if (std::abs(alpha)> safety) {
+ x = origin;
+ x += linearPart;
+ x -= f.origin();
+ }
#endif
}
};
diff --git a/src/time-stepping/rate/newmark.cc b/src/time-stepping/rate/newmark.cc
index 0ec590126ea5fd19755e6699324a259b5b5daeb8..dec66b98d732c2589b7bc79c4c371a47de336714 100644
--- a/src/time-stepping/rate/newmark.cc
+++ b/src/time-stepping/rate/newmark.cc
@@ -108,7 +108,7 @@ void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
}
}
}
- //sprint(iterate, "iterate:");
+ //print(iterate, "iterate:");
}
template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>