From afbcb7977e6ea3e0177662c19960ba2cbb46dc6c Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Tue, 16 Jul 2019 16:53:06 +0200
Subject: [PATCH] .
---
src/multi-body-problem-2D.cfg | 2 +-
src/multi-body-problem.cc | 2 ++
src/spatial-solving/tnnmg/linearcorrection.hh | 25 ++++++++++++----
src/spatial-solving/tnnmg/tnnmgstep.hh | 30 +++++++++++--------
src/time-stepping/adaptivetimestepper.cc | 10 +++----
5 files changed, 45 insertions(+), 24 deletions(-)
diff --git a/src/multi-body-problem-2D.cfg b/src/multi-body-problem-2D.cfg
index 9112827c..676d1a40 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.1 # 2e-3 [m]
+smallestDiameter = 0.5 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5 # 1e-5
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 6e32d93b..5e4c505e 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -199,6 +199,8 @@ int main(int argc, char *argv[]) {
programState.setupInitialConditions(parset, contactNetwork);
+ DUNE_THROW(Dune::Exception, "Just need to stop here!");
+
auto& nBodyAssembler = contactNetwork.nBodyAssembler();
for (size_t i=0; i<bodyCount; i++) {
contactNetwork.body(i)->setDeformation(programState.u[i]);
diff --git a/src/spatial-solving/tnnmg/linearcorrection.hh b/src/spatial-solving/tnnmg/linearcorrection.hh
index 2b1d5752..d06155e5 100644
--- a/src/spatial-solving/tnnmg/linearcorrection.hh
+++ b/src/spatial-solving/tnnmg/linearcorrection.hh
@@ -25,6 +25,17 @@
template<typename Matrix, typename Vector>
using LinearCorrection = std::function<void(const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore)>;
+template<typename Vector>
+Dune::Solvers::DefaultBitVector_t<Vector>
+emptyIgnore(const Vector& v)
+{
+ // TNNMGStep assumes that the linearization and the solver for the
+ // linearized problem will not use the ignoreNodes field
+ Dune::Solvers::DefaultBitVector_t<Vector> ignore;
+ Dune::Solvers::resizeInitialize(ignore, v, false);
+ return ignore;
+}
+
template<typename Matrix, typename Vector>
LinearCorrection<Matrix, Vector>
makeLinearCorrection(std::shared_ptr< Dune::Solvers::LinearSolver<Matrix, Vector> > linearSolver)
@@ -65,23 +76,25 @@ makeLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > i
directSolver->preprocess();
directSolver->solve();
- x = directX;
- /* using LinearIterationStep = Dune::Solvers::LinearIterationStep<Matrix, Vector>;
+ //x = directX;
+ using LinearIterationStep = Dune::Solvers::LinearIterationStep<Matrix, Vector>;
auto linearIterationStep = dynamic_cast<LinearIterationStep*>(&iterativeSolver->getIterationStep());
if (not linearIterationStep)
DUNE_THROW(Dune::Exception, "iterative solver must use a linear iteration step");
- linearIterationStep->setIgnore(ignore);
+ auto empty = emptyIgnore(x);
+
+ linearIterationStep->setIgnore(empty);
linearIterationStep->setProblem(A, x, b);
iterativeSolver->preprocess();
- iterativeSolver->solve();*/
+ iterativeSolver->solve();
+
- /*
const auto& norm = iterativeSolver->getErrorNorm();
auto error = norm.diff(linearIterationStep->getSol(), directX);
- std::cout << "CG Solver iterations: " << iterativeSolver->getResult().iterations << " Error: " << error << std::endl;*/
+ std::cout << "Linear solver iterations: " << iterativeSolver->getResult().iterations << " Error: " << error << std::endl;
};
}
diff --git a/src/spatial-solving/tnnmg/tnnmgstep.hh b/src/spatial-solving/tnnmg/tnnmgstep.hh
index 84d7e914..d09f3fd1 100644
--- a/src/spatial-solving/tnnmg/tnnmgstep.hh
+++ b/src/spatial-solving/tnnmg/tnnmgstep.hh
@@ -19,7 +19,8 @@
#include "localproblem.hh"
-#include "linearcorrection.hh"
+//#include "linearcorrection.hh"
+#include <dune/tnnmg/iterationsteps/linearcorrection.hh>
#include <dune/tectonic/../../src/utils/debugutils.hh>
@@ -144,7 +145,7 @@ class TNNMGStep :
const auto& ignore = (*this->ignoreNodes_);
auto& x = *getIterate();
- //std::cout << "TNNMGStep::iterate " << std::endl;
+ std::cout << "TNNMGStep::iterate " << std::endl;
//print(f.quadraticPart(), "f.quadraticPart():");
@@ -156,9 +157,12 @@ class TNNMGStep :
for (std::size_t i=0; i<preSmoothingSteps_; ++i)
nonlinearSmoother_->iterate();
- //std::cout << "- nonlinear presmoothing: success" << std::endl;
+ std::cout << "- nonlinear presmoothing: success" << std::endl;
//print(x, "TNNMG iterate after smoothing:");
+ std::cout << "-- energy after presmoothing: " << energy(x) << std::endl;
+
+ /*
// Compute constraint/truncated linearization
if (not linearization_)
linearization_ = std::make_shared<Linearization>(f, ignore);
@@ -168,32 +172,32 @@ class TNNMGStep :
auto&& A = linearization_->hessian();
auto&& r = linearization_->negativeGradient();
- /*print(A, "TNNMG Linear problem matrix:");
+ print(A, "TNNMG Linear problem matrix:");
print(r, "TNNMG Linear problem rhs:");
print(ignore, "ignore:");
- print(linearization_->truncated(), "truncation:");*/
+ print(linearization_->truncated(), "truncation:");
// Compute inexact solution of the linearized problem
Solvers::resizeInitializeZero(correction_, x);
Solvers::resizeInitializeZero(constrainedCorrection_, r);
- //std::cout << "-- energy after presmoothing: " << energy(x) << std::endl;
+
//print(constrainedCorrection_compare, "direct solver solution:");
- //std::cout << "- computing linear correction..." << std::endl;
+ std::cout << "- computing linear correction..." << std::endl;
- linearCorrection_(A, constrainedCorrection_, r, ignore);
+ //linearCorrection_(A, constrainedCorrection_, r, ignore);
//DUNE_THROW(Dune::Exception, "TNNMGStep: Just need to terminate here!");
- // linearCorrection_(A, constrainedCorrection_, r);
+ linearCorrection_(A, constrainedCorrection_, r);
//std::cout << "- linear correction: success" << std::endl;
//print(constrainedCorrection_, "contrained correction:");
- linearization_->extendCorrection(constrainedCorrection_, correction_);
+ linearization_->extendCorrection(constrainedCorrection_, correction_); */
/*Vector h = x;
h += correction_;
@@ -202,6 +206,7 @@ class TNNMGStep :
//std::cout << "- extended correction: success" << std::endl;
//print(correction_, "correction:");
+ /*
// Project onto admissible set
projection_(f, x, correction_);
@@ -214,7 +219,7 @@ class TNNMGStep :
// Line search
- auto fv = directionalRestriction(f, x, correction_);
+ auto fv = directionalRestriction(f, x, correction_); */
/* std::cout << fv.quadraticPart() << " f.quadraticPart():" << std::endl;
std::cout << fv.linearPart() << " f.linearPart():" << std::endl;
@@ -222,6 +227,7 @@ class TNNMGStep :
std::cout << fv.domain()[0] << " " << fv.domain()[0] << " f.domain():" << std::endl;
std::cout << "- setup directional restriction: success" << std::endl; */
+ /*
dampingFactor_ = 0;
lineSolver_(dampingFactor_, fv, false);
@@ -232,7 +238,7 @@ class TNNMGStep :
//std::cout << "damping factor: " << dampingFactor_ << std::endl;
- x += correction_;
+ x += correction_;*/
//std::cout << "-- energy after correction: " << energy(x) << std::endl;
}
diff --git a/src/time-stepping/adaptivetimestepper.cc b/src/time-stepping/adaptivetimestepper.cc
index c1869aa8..3a349a7a 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -67,7 +67,7 @@ IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNo
std::cout << "AdaptiveTimeStepper Step 1" << std::endl;
-
+/*
bool didCoarsen = false;
iterationRegister_.reset();
UpdatersWithCount R2;
@@ -105,16 +105,16 @@ IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNo
R2 = F2;
}
}
-
+ */
std::cout << "AdaptiveTimeStepper::advance() ...";
iterationRegister_.registerFinalCount(R1_.count);
relativeTime_ += relativeTau_;
current_ = R1_.updaters;
- //UpdatersWithCount emptyR1;
- //R1_ = emptyR1;
- R1_ = R2;
+ UpdatersWithCount emptyR1;
+ R1_ = emptyR1;
+ //R1_ = R2;
std::cout << " done" << std::endl;
--
GitLab