Skip to content
Snippets Groups Projects
Commit afbcb797 authored by podlesny's avatar podlesny
Browse files

.

parent b0c3b714
No related branches found
No related tags found
No related merge requests found
# -*- mode:conf -*-
[boundary.friction]
smallestDiameter = 0.1 # 2e-3 [m]
smallestDiameter = 0.5 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5 # 1e-5
......
......@@ -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]);
......
......@@ -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;
};
}
......
......@@ -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;
}
......
......@@ -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;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment