From 1f470d23cede0963f5ac354acdcf2aff260b8082 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Fri, 19 Jul 2019 11:16:24 +0200
Subject: [PATCH] .
---
src/data-structures/contactnetwork.cc | 2 +-
src/multi-body-problem.cc | 24 ++++-----
src/spatial-solving/fixedpointiterator.cc | 50 +++++++++++++------
src/spatial-solving/fixedpointiterator.hh | 9 +++-
src/spatial-solving/solverfactory.cc | 2 +
src/spatial-solving/solverfactory.hh | 2 +-
src/spatial-solving/solverfactory_tmpl.cc | 16 +++---
src/spatial-solving/tnnmg/tnnmgstep.hh | 18 +++----
src/time-stepping/adaptivetimestepper_tmpl.cc | 1 -
src/time-stepping/rate/newmark.cc | 6 +--
src/utils/debugutils.hh | 8 +--
11 files changed, 79 insertions(+), 59 deletions(-)
diff --git a/src/data-structures/contactnetwork.cc b/src/data-structures/contactnetwork.cc
index 4e66cf76..1a54dda9 100644
--- a/src/data-structures/contactnetwork.cc
+++ b/src/data-structures/contactnetwork.cc
@@ -26,7 +26,7 @@ ContactNetwork<HostGridType, VectorType>::ContactNetwork(
couplings_(nCouplings),
frictionBoundaries_(nBodies),
stateEnergyNorms_(nBodies),
- nBodyAssembler_(nBodies, nCouplings)
+ nBodyAssembler_(nBodies, nCouplings, false, 0.0)
{}
template <class HostGridType, class VectorType>
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index f986549e..4582a652 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -136,22 +136,23 @@ int main(int argc, char *argv[]) {
const size_t bodyCount = contactNetwork.nBodies();
- /*
- for (size_t i=0; i<bodyCount; i++) {
+ for (size_t i=0; i<contactNetwork.nLevels(); i++) {
// printDofLocation(contactNetwork.body(i)->gridView());
//Vector def(contactNetwork.deformedGrids()[i]->size(dims));
//def = 1;
//deformedGridComplex.setDeformation(def, i);
- auto& levelViews = contactNetwork.levelViews(i);
+ const auto& level = *contactNetwork.level(i);
- for (size_t j=0; j<levelViews.size(); j++) {
- writeToVTK(*levelViews[j], "", "body_" + std::to_string(i) + "_level_" + std::to_string(j));
+ for (size_t j=0; j<level.nBodies(); j++) {
+ writeToVTK(level.body(j)->gridView(), "debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
}
+ }
- writeToVTK(contactNetwork.leafView(i), "", "body_" + std::to_string(i) + "_leaf");
- } */
+ for (size_t i=0; i<bodyCount; i++) {
+ writeToVTK(contactNetwork.body(i)->gridView(), "debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
+ }
// ----------------------------
// assemble contactNetwork
@@ -307,8 +308,9 @@ int main(int argc, char *argv[]) {
}
}*/
- //print(totalDirichletNodes, "totalDirichletNodes:");
+ print(totalDirichletNodes, "totalDirichletNodes:");
+ //using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, field_type>;
using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
using NonlinearFactory = SolverFactory<Functional, BitVector>;
@@ -332,9 +334,9 @@ int main(int argc, char *argv[]) {
std::vector<const Dune::BitSetVector<1>*> frictionNodes;
contactNetwork.frictionNodes(frictionNodes);
- /*for (size_t i=0; i<frictionNodes.size(); i++) {
+ for (size_t i=0; i<frictionNodes.size(); i++) {
print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
- }*/
+ }
Updaters current(
initRateUpdater(
@@ -449,8 +451,6 @@ int main(int argc, char *argv[]) {
report();
- break;
-
if (programState.timeStep==50) {
std::cout << "limit of timeSteps reached!" << std::endl;
break; // TODO remove after debugging
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index 633eb1aa..d84f5f09 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -36,6 +36,10 @@
#include "../spatial-solving/preconditioners/nbodycontacttransfer.hh"
+#include "tnnmg/functional.hh"
+#include "tnnmg/zerononlinearity.hh"
+#include "solverfactory.hh"
+
void FixedPointIterationCounter::operator+=(
FixedPointIterationCounter const &other) {
iterations += other.iterations;
@@ -74,6 +78,12 @@ FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
const auto& nBodyAssembler = contactNetwork_.nBodyAssembler();
+ // debugging
+ const auto& contactCouplings = nBodyAssembler.getContactCouplings();
+ for (size_t i=0; i<contactCouplings.size(); i++) {
+ print(*contactCouplings[i]->nonmortarBoundary().getVertices(), "nonmortarBoundaries:");
+ }
+
const auto nBodies = nBodyAssembler.nGrids();
std::vector<const Matrix*> matrices_ptr(nBodies);
@@ -85,12 +95,12 @@ FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
Matrix bilinearForm;
nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
- //print(bilinearForm, "bilinearForm:");
+ print(bilinearForm, "bilinearForm:");
Vector totalRhs;
nBodyAssembler.assembleRightHandSide(velocityRHSs, totalRhs);
- //print(totalRhs, "totalRhs:");
+ print(totalRhs, "totalRhs:");
// get lower and upper obstacles
const auto& totalObstacles = nBodyAssembler.totalObstacles_;
@@ -107,12 +117,12 @@ FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
}
}
- //print(totalObstacles, "totalObstacles:");
+ print(totalObstacles, "totalObstacles:");
- //print(lower, "lower obstacles:");
- //print(upper, "upper obstacles:");
+ print(lower, "lower obstacles:");
+ print(upper, "upper obstacles:");
- // comput velocity obstacles
+ // compute velocity obstacles
Vector vLower, vUpper;
std::vector<Vector> u0, v0;
updaters.rate_->extractOldVelocity(v0);
@@ -125,8 +135,8 @@ FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
updaters.rate_->velocityObstacles(totalu0, lower, totalv0, vLower);
updaters.rate_->velocityObstacles(totalu0, upper, totalv0, vUpper);
- //print(vLower, "vLower obstacles:");
- //print(vUpper, "vUpper obstacles:");
+ print(vLower, "vLower obstacles:");
+ print(vUpper, "vUpper obstacles:");
std::cout << "- Problem assembled: success" << std::endl;
@@ -147,16 +157,20 @@ FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
auto smoother = TruncatedBlockGSStep<Matrix, Vector>{};
auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
linearMultigridStep->setMGType(1, 3, 3);
- linearMultigridStep->setSmoother(&smoother);
+ linearMultigridStep->setSmoother(smoother);
linearMultigridStep->setTransferOperators(transfer);
EnergyNorm<Matrix, Vector> mgNorm(*linearMultigridStep);
LinearSolver mgSolver(linearMultigridStep, parset_.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset_.get<double>("solver.tnnmg.linear.tolerance"), mgNorm, Solver::QUIET);
+ print(ignoreNodes_, "ignoreNodes:");
// set up functional and TNMMG solver
- Functional J(bilinearForm, totalRhs, globalFriction_, vLower, vUpper);
- Factory solverFactory(parset_.sub("solver.tnnmg"), J, mgSolver, ignoreNodes_);
+ using ZeroSolverFactory = SolverFactory<Functional, IgnoreVector>;
+ Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), vLower, vUpper);
+ ZeroSolverFactory solverFactory(parset_.sub("solver.tnnmg"), J, mgSolver, ignoreNodes_);
+ /*Functional J(bilinearForm, totalRhs, globalFriction_, vLower, vUpper);
+ Factory solverFactory(parset_.sub("solver.tnnmg"), J, mgSolver, ignoreNodes_);*/
auto step = solverFactory.step();
std::cout << "- Functional and TNNMG step setup: success" << std::endl;
@@ -173,6 +187,8 @@ FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
++fixedPointIteration) {
+ print(alpha, "alpha:");
+
// contribution from nonlinearity
globalFriction_.updateAlpha(alpha);
@@ -197,8 +213,10 @@ FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
nBodyAssembler.postprocess(tnnmgSol, velocityIterates);
//nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
- //print(totalVelocityIterate, "totalVelocityIterate:");
- //print(velocityIterates, "velocityIterates:");
+ print(totalVelocityIterate, "totalVelocityIterate:");
+ print(velocityIterates, "velocityIterates:");
+
+ //DUNE_THROW(Dune::Exception, "Just need to stop here!");
multigridIterations += velocityProblemSolver.getResult().iterations;
@@ -227,12 +245,12 @@ FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
updaters.state_->extractAlpha(newAlpha);
bool breakCriterion = true;
- for (size_t i=0; i<nBodies; i++) {
+ for (int i=0; i<nBodies; i++) {
if (alpha[i].size()==0 || newAlpha[i].size()==0)
continue;
- //print(alpha[i], "alpha i:");
- //print(newAlpha[i], "new alpha i:");
+ print(alpha[i], "alpha i:");
+ print(newAlpha[i], "new alpha i:");
if (errorNorms_[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
breakCriterion = false;
std::cout << "fixedPoint error: " << errorNorms_[i]->diff(alpha[i], newAlpha[i]) << std::endl;
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index c4d59c3d..c4839a79 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -13,6 +13,10 @@
#include <dune/contact/assemblers/nbodyassembler.hh>
+#include "tnnmg/functional.hh"
+#include "tnnmg/zerononlinearity.hh"
+#include "solverfactory.hh"
+
struct FixedPointIterationCounter {
size_t iterations = 0;
size_t multigridIterations = 0;
@@ -25,11 +29,12 @@ std::ostream &operator<<(std::ostream &stream,
template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
class FixedPointIterator {
- using Functional = typename Factory::Functional;
using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix;
- using Nonlinearity = typename Factory::Functional::Nonlinearity;
+
+ using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, double>; //typename Factory::Functional;
+ using Nonlinearity = typename Factory::Functional::Nonlinearity;
const static int dims = Vector::block_type::dimension;
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index 6b9a3119..cfe36305 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -6,6 +6,8 @@
#include <dune/solvers/iterationsteps/blockgssteps.hh>
#include <dune/solvers/solvers/umfpacksolver.hh>
+//#include "solverfactory.hh"
+
#include "../utils/debugutils.hh"
template <class Functional, class BitVector>
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index 48e4a63c..5d4a5aa2 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -36,7 +36,7 @@ class SolverFactory {
template <class LinearSolver>
- SolverFactory(Dune::ParameterTree const &,
+ SolverFactory(const Dune::ParameterTree&,
Functional&,
LinearSolver&&,
const BitVector&);
diff --git a/src/spatial-solving/solverfactory_tmpl.cc b/src/spatial-solving/solverfactory_tmpl.cc
index 488c5191..a136230f 100644
--- a/src/spatial-solving/solverfactory_tmpl.cc
+++ b/src/spatial-solving/solverfactory_tmpl.cc
@@ -10,27 +10,23 @@
#include "../../dune/tectonic/globalfriction.hh"
#include "tnnmg/functional.hh"
#include "tnnmg/zerononlinearity.hh"
+
#include "solverfactory.hh"
using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
+
using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
using MyZeroFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, double>;
-using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
-using MyLinearSolver = Dune::Solvers::LoopSolver<Vector, BitVector>;
+using MyLinearSolver = Dune::Solvers::LoopSolver<Vector>;
using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
template class SolverFactory<MyFunctional, BitVector>;
+template<> template<> SolverFactory<MyFunctional, BitVector>::SolverFactory(const Dune::ParameterTree&, MyFunctional&, MyLinearSolver&&, const BitVector&);
-/*template<> SolverFactory<MyFunctional, BitVector>::SolverFactory(Dune::ParameterTree const &,
- MyFunctional&,
- MyLinearSolver&&,
- const BitVector&);*/
-
+using MyZeroSolverFactory = SolverFactory<MyZeroFunctional, BitVector>;
template class SolverFactory<MyZeroFunctional, BitVector>;
-
-/*template<>
-template<> SolverFactory<MyZeroFunctional, BitVector>::SolverFactory(Dune::ParameterTree const &,
+/*template<> SolverFactory<MyZeroFunctional, BitVector>::SolverFactory(Dune::ParameterTree const &,
MyZeroFunctional&,
MyLinearSolver&&,
const BitVector&);*/
diff --git a/src/spatial-solving/tnnmg/tnnmgstep.hh b/src/spatial-solving/tnnmg/tnnmgstep.hh
index 41780c07..844f0137 100644
--- a/src/spatial-solving/tnnmg/tnnmgstep.hh
+++ b/src/spatial-solving/tnnmg/tnnmgstep.hh
@@ -153,16 +153,17 @@ class TNNMGStep :
std::cout << "-- energy before smoothing: " << f(x) << std::endl;
+ print(x, "TNNMG iterate after smoothing:");
+
// Nonlinear presmoothing
for (std::size_t i=0; i<preSmoothingSteps_; ++i)
nonlinearSmoother_->iterate();
- std::cout << "- nonlinear presmoothing: success" << std::endl;
- //print(x, "TNNMG iterate after smoothing:");
+ //std::cout << "- nonlinear presmoothing: success" << std::endl;
+ print(x, "TNNMG iterate after smoothing:");
std::cout << "-- energy after presmoothing: " << f(x) << std::endl;
- /*
// Compute constraint/truncated linearization
if (not linearization_)
linearization_ = std::make_shared<Linearization>(f, ignore);
@@ -197,7 +198,7 @@ class TNNMGStep :
//print(constrainedCorrection_, "contrained correction:");
- linearization_->extendCorrection(constrainedCorrection_, correction_); */
+ linearization_->extendCorrection(constrainedCorrection_, correction_);
/*Vector h = x;
h += correction_;
@@ -206,7 +207,6 @@ class TNNMGStep :
//std::cout << "- extended correction: success" << std::endl;
//print(correction_, "correction:");
- /*
// Project onto admissible set
projection_(f, x, correction_);
@@ -219,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;
@@ -227,7 +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);
@@ -238,9 +238,9 @@ class TNNMGStep :
//std::cout << "damping factor: " << dampingFactor_ << std::endl;
- x += correction_;*/
+ x += correction_;
- //std::cout << "-- energy after correction: " << energy(x) << std::endl;
+ std::cout << "-- energy after correction: " << energy(x) << std::endl;
}
/**
diff --git a/src/time-stepping/adaptivetimestepper_tmpl.cc b/src/time-stepping/adaptivetimestepper_tmpl.cc
index ff76be34..fcb100f4 100644
--- a/src/time-stepping/adaptivetimestepper_tmpl.cc
+++ b/src/time-stepping/adaptivetimestepper_tmpl.cc
@@ -24,5 +24,4 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
using MyAdaptiveTimeStepper = AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;
-
template class AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;
diff --git a/src/time-stepping/rate/newmark.cc b/src/time-stepping/rate/newmark.cc
index dec66b98..a3b36385 100644
--- a/src/time-stepping/rate/newmark.cc
+++ b/src/time-stepping/rate/newmark.cc
@@ -95,20 +95,20 @@ void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
double dirichletValue;
bodyDirichletFunctions[bc]->evaluate(relativeTime, dirichletValue);
- //std::cout << "dirichletValue: " << dirichletValue << std::endl;
+ std::cout << "dirichletValue: " << dirichletValue << std::endl;
for (size_t k=0; k<bcDirichletNodes.size() ; ++k) {
for (size_t j=0; j<dim; ++j) {
if (bcDirichletNodes[k][j]) {
- //std::cout << k << " " << j << std::endl;
+ std::cout << k << " " << j << std::endl;
bodyIterate[k][j] = dirichletValue;
}
}
}
}
}
- //print(iterate, "iterate:");
+ print(iterate, "iterate:");
}
template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
diff --git a/src/utils/debugutils.hh b/src/utils/debugutils.hh
index c14b16f1..e80dafee 100644
--- a/src/utils/debugutils.hh
+++ b/src/utils/debugutils.hh
@@ -118,8 +118,8 @@
void print(const std::vector<T>& x, std::string message){
std::cout << message << std::endl;
for (size_t i=0; i<x.size(); i++) {
- //std::cout << x[i] << std::endl;
- print(x[i], "");
+ std::cout << x[i] << std::endl;
+ //print(x[i], "");
}
std::cout << std::endl << std::endl;
}
@@ -165,7 +165,7 @@
std::cout.rdbuf(lStream.rdbuf());
vtkwriter.addVertexData(x,"data");
- vtkwriter.pwrite(name, path, path);
+ vtkwriter.pwrite(name, path, "");
std::cout.rdbuf( lBufferOld );
}
@@ -178,7 +178,7 @@
std::streambuf* lBufferOld = std::cout.rdbuf();
std::cout.rdbuf(lStream.rdbuf());
- vtkwriter.pwrite(name, path, path);
+ vtkwriter.pwrite(name, path, "");
std::cout.rdbuf( lBufferOld );
}
--
GitLab