diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index e9db4febf7871c3cfdbb9aae21e5c3209ccd4f83..6fbc95c173549a5a73470dae08e0a95b5c8e0583 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -25,6 +25,12 @@
#include "../spatial-solving/tnnmg/zerononlinearity.hh"
#include "../utils/debugutils.hh"
+
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+
+#include "../spatial-solving/preconditioners/nbodycontacttransfer.hh"
+
template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class BodyState {
public:
using Vector = VectorTEMPLATE;
@@ -98,8 +104,17 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
const auto& preconditionerParset = parset.sub("solver.tnnmg.linear.preconditioner");
+ std::cout << "-- setupInitialConditions --" << std::endl;
+ std::cout << "----------------------------" << std::endl;
+
+ /*
+ std::cout << "Building preconditioner..." << std::endl;
using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
+ /*Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), false);
+ activeLevels.back() = true;
+ activeLevels[activeLevels.size()-2] = true;*/
+ /*print(activeLevels, "activeLevels:");
Preconditioner preconditioner(preconditionerParset, contactNetwork, activeLevels);
preconditioner.build();
@@ -110,7 +125,31 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
linearSolverStep.setPreconditioner(preconditioner);
EnergyNorm<Matrix, Vector> energyNorm(linearSolverStep);
- LinearSolver linearSolver(linearSolverStep, parset.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset.get<double>("solver.tnnmg.linear.tolerance"), energyNorm, Solver::QUIET);
+ LinearSolver linearSolver(linearSolverStep, parset.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset.get<double>("solver.tnnmg.linear.tolerance"), energyNorm, Solver::FULL);
+ */
+ using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+ using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+ using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+ TransferOperators transfer(contactNetwork.nLevels()-1);
+ for (size_t i=0; i<transfer.size(); i++) {
+ transfer[i] = std::make_shared<TransferOperator>();
+ transfer[i]->setup(contactNetwork, i, i+1);
+ }
+
+ // Remove any recompute filed so that initially the full transferoperator is assembled
+ for (size_t i=0; i<transfer.size(); i++)
+ std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+ auto smoother = TruncatedBlockGSStep<Matrix, Vector>{};
+ auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+ linearMultigridStep->setMGType(1, 3, 3);
+ 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::FULL);
+
// Solving a linear problem with a multigrid solver
auto const solveLinearProblem = [&](
@@ -177,7 +216,8 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
// set up TNMMG solver
using Factory = SolverFactory<Functional, BitVector>;
- Factory factory(parset.sub("solver.tnnmg"), J, linearSolver, _dirichletNodes);
+ //Factory factory(parset.sub("solver.tnnmg"), J, linearSolver, _dirichletNodes);
+ Factory factory(parset.sub("solver.tnnmg"), J, mgSolver, _dirichletNodes);
/* std::vector<BitVector> bodyDirichletNodes;
nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
@@ -270,10 +310,12 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
}
}*/
+ std::cout << "solving linear problem for u..." << std::endl;
+
solveLinearProblem(dirichletNodes, contactNetwork.matrices().elasticity, ell0, u,
parset.sub("u0.solver"));
- print(u, "initial u:");
+ //print(u, "initial u:");
// Initial acceleration: Computed in agreement with Ma = ell0 - Au
// (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
@@ -300,11 +342,13 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, u[i]);
}
+ std::cout << "solving linear problem for a..." << std::endl;
+
BitVector noNodes(dirichletNodes.size(), false);
solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, a,
parset.sub("a0.solver"));
- print(a, "initial a:");
+ //print(a, "initial a:");
}
private:
diff --git a/src/multi-body-problem-2D.cfg b/src/multi-body-problem-2D.cfg
index c0b758a8abf73e4ff14785138db64bf7afd2a99b..4fc63ee9679dee491c37c99fd18ea8694f476fa1 100644
--- a/src/multi-body-problem-2D.cfg
+++ b/src/multi-body-problem-2D.cfg
@@ -21,4 +21,4 @@ tolerance = 1e-5
tolerance = 1e-10
[solver.tnnmg.linear.preconditioner]
-tolerance = 1e-10
+tolerance = 1e-8
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 4215f278ae91870ffba03d7c21cd2eddf3a26e90..f7ef93e5d7857bff2473da67f5a0ff87385107d7 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -190,7 +190,7 @@ int main(int argc, char *argv[]) {
}
contactNetwork.setDeformation(def);*/
- const auto & coarseContactNetwork = *contactNetwork.level(3);
+ /* const auto & coarseContactNetwork = *contactNetwork.level(3);
const auto & fineContactNetwork = *contactNetwork.level(4);
// SupportPatchFactory<typename ContactNetwork::LevelContactNetwork> supportPatchFactory(coarseContactNetwork, fineContactNetwork);
size_t bodyID = 2;
@@ -222,7 +222,7 @@ int main(int argc, char *argv[]) {
dofs[j] = j;
}
writeToVTK(gv, dofs, "", "body_" + std::to_string(i) + "_fine");
- }
+ } */
/*
using Patch = typename SupportPatchFactory<typename ContactNetwork::LevelContactNetwork>::Patch;
@@ -267,9 +267,6 @@ int main(int argc, char *argv[]) {
}
} */
- std::cout << "Done! :)" << std::endl;
- //return 1;
-
//printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
// -----------------
@@ -417,7 +414,7 @@ int main(int argc, char *argv[]) {
}
}*/
- print(totalDirichletNodes, "totalDirichletNodes:");
+ //print(totalDirichletNodes, "totalDirichletNodes:");
using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
using NonlinearFactory = SolverFactory<Functional, BitVector>;
@@ -433,18 +430,18 @@ int main(int argc, char *argv[]) {
BoundaryNodes dirichletNodes;
contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
- for (size_t i=0; i<dirichletNodes.size(); i++) {
+ /*for (size_t i=0; i<dirichletNodes.size(); i++) {
for (size_t j=0; j<dirichletNodes[i].size(); j++) {
print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
}
- }
+ }*/
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(
@@ -508,7 +505,7 @@ int main(int argc, char *argv[]) {
std::signal(SIGINT, handleSignal);
std::signal(SIGTERM, handleSignal);
- /*
+/*
// set patch preconditioner for linear correction in TNNMG method
using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, PatchSolver, Matrix, Vector>;
@@ -527,21 +524,21 @@ int main(int argc, char *argv[]) {
Preconditioner preconditioner(contactNetwork, activeLevels, preconditionerParset.get<Preconditioner::Mode>("mode"));
preconditioner.setPatchSolver(patchSolver);
preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
-
+*/
// set adaptive time stepper
typename ContactNetwork::ExternalForces externalForces;
contactNetwork.externalForces(externalForces);
- AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(nBodyAssembler)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
- adaptiveTimeStepper(parset, nBodyAssembler, totalDirichletNodes, globalFriction, frictionNodes, current,
+ AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+ adaptiveTimeStepper(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes, current,
programState.relativeTime, programState.relativeTau,
externalForces, stateEnergyNorms, mustRefine);
while (!adaptiveTimeStepper.reachedEnd()) {
programState.timeStep++;
- preconditioner.build();
- iterationCount = adaptiveTimeStepper.advance(preconditioner);
+ //preconditioner.build();
+ iterationCount = adaptiveTimeStepper.advance();
programState.relativeTime = adaptiveTimeStepper.relativeTime_;
programState.relativeTau = adaptiveTimeStepper.relativeTau_;
@@ -559,6 +556,8 @@ 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
@@ -569,7 +568,7 @@ int main(int argc, char *argv[]) {
break;
}
}
- */
+
std::cout.rdbuf(coutbuf); //reset to standard output again
diff --git a/src/multi-body-problem.cfg b/src/multi-body-problem.cfg
index 8678a5ddf21e17c9da8fabd9253dee59868ae1d2..7188eed1b5851c795b2903d8000846c4a89c6450 100644
--- a/src/multi-body-problem.cfg
+++ b/src/multi-body-problem.cfg
@@ -11,7 +11,7 @@ vtk.write = true
[problem]
finalTime = 100 # [s] #1000
-bodyCount = 4
+bodyCount = 2
[body]
bulkModulus = 0.5e5 # [Pa]
@@ -49,11 +49,11 @@ relativeTau = 1e-4 # 1e-6
scheme = newmark
[u0.solver]
-maximumIterations = 10000
+maximumIterations = 100
verbosity = full
[a0.solver]
-maximumIterations = 10000
+maximumIterations = 100
verbosity = full
[v.solver]
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index 99456a3be28ad70481fda8cd4b59703215eae5f5..633eb1aa03df1c5bf50630d0a948b6ebbe2a0d5e 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -30,6 +30,11 @@
#include "../utils/tobool.hh"
#include "../utils/debugutils.hh"
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+
+#include "../spatial-solving/preconditioners/nbodycontacttransfer.hh"
void FixedPointIterationCounter::operator+=(
FixedPointIterationCounter const &other) {
@@ -37,16 +42,16 @@ void FixedPointIterationCounter::operator+=(
multigridIterations += other.multigridIterations;
}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::FixedPointIterator(
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::FixedPointIterator(
Dune::ParameterTree const &parset,
- NBodyAssembler& nBodyAssembler,
+ const ContactNetwork& contactNetwork,
const IgnoreVector& ignoreNodes,
GlobalFriction& globalFriction,
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
const ErrorNorms& errorNorms)
: parset_(parset),
- nBodyAssembler_(nBodyAssembler),
+ contactNetwork_(contactNetwork),
ignoreNodes_(ignoreNodes),
globalFriction_(globalFriction),
bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
@@ -58,17 +63,18 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::FixedPointIte
verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
errorNorms_(errorNorms) {}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-template <class Preconditioner>
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
FixedPointIterationCounter
-FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
- Updaters updaters, Preconditioner&& preconditioner,
+FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
+ Updaters updaters,
const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates) {
std::cout << "FixedPointIterator::run()" << std::endl;
- const auto nBodies = nBodyAssembler_.nGrids();
+ const auto& nBodyAssembler = contactNetwork_.nBodyAssembler();
+
+ const auto nBodies = nBodyAssembler.nGrids();
std::vector<const Matrix*> matrices_ptr(nBodies);
for (int i=0; i<nBodies; i++) {
@@ -77,17 +83,17 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
// assemble full global contact problem
Matrix bilinearForm;
- nBodyAssembler_.assembleJacobian(matrices_ptr, bilinearForm);
+ nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
//print(bilinearForm, "bilinearForm:");
Vector totalRhs;
- nBodyAssembler_.assembleRightHandSide(velocityRHSs, totalRhs);
+ nBodyAssembler.assembleRightHandSide(velocityRHSs, totalRhs);
//print(totalRhs, "totalRhs:");
// get lower and upper obstacles
- const auto& totalObstacles = nBodyAssembler_.totalObstacles_;
+ const auto& totalObstacles = nBodyAssembler.totalObstacles_;
Vector lower(totalObstacles.size());
Vector upper(totalObstacles.size());
@@ -101,10 +107,10 @@ FixedPointIterator<Factory, NBodyAssembler, 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
Vector vLower, vUpper;
@@ -113,20 +119,44 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
updaters.rate_->extractOldDisplacement(u0);
Vector totalu0, totalv0;
- nBodyAssembler_.concatenateVectors(u0, totalu0);
- nBodyAssembler_.concatenateVectors(v0, totalv0);
+ nBodyAssembler.concatenateVectors(u0, totalu0);
+ nBodyAssembler.concatenateVectors(v0, totalv0);
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;
+ using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, IgnoreVector>;
+ using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+ using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+ TransferOperators transfer(contactNetwork_.nLevels()-1);
+ for (size_t i=0; i<transfer.size(); i++) {
+ transfer[i] = std::make_shared<TransferOperator>();
+ transfer[i]->setup(contactNetwork_, i, i+1);
+ }
+
+ // Remove any recompute filed so that initially the full transferoperator is assembled
+ for (size_t i=0; i<transfer.size(); i++)
+ std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+ auto smoother = TruncatedBlockGSStep<Matrix, Vector>{};
+ auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+ linearMultigridStep->setMGType(1, 3, 3);
+ 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);
+
+
// set up functional and TNMMG solver
Functional J(bilinearForm, totalRhs, globalFriction_, vLower, vUpper);
- Factory solverFactory(parset_.sub("solver.tnnmg"), J, std::forward<Preconditioner>(preconditioner), ignoreNodes_);
+ Factory solverFactory(parset_.sub("solver.tnnmg"), J, mgSolver, ignoreNodes_);
auto step = solverFactory.step();
std::cout << "- Functional and TNNMG step setup: success" << std::endl;
@@ -147,7 +177,7 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
globalFriction_.updateAlpha(alpha);
Vector totalVelocityIterate;
- nBodyAssembler_.nodalToTransformed(velocityIterates, totalVelocityIterate);
+ nBodyAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
//print(velocityIterates, "velocityIterates:");
//print(totalVelocityIterate, "totalVelocityIterate:");
@@ -164,8 +194,8 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
std::cout << "-- Solved" << std::endl;
const auto& tnnmgSol = step->getSol();
- nBodyAssembler_.postprocess(tnnmgSol, velocityIterates);
- //nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
+ nBodyAssembler.postprocess(tnnmgSol, velocityIterates);
+ //nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
//print(totalVelocityIterate, "totalVelocityIterate:");
//print(velocityIterates, "velocityIterates:");
@@ -239,13 +269,14 @@ std::ostream &operator<<(std::ostream &stream,
<< ")";
}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::relativeVelocities(
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+void FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::relativeVelocities(
const Vector& v,
std::vector<Vector>& v_rel) const {
- const size_t nBodies = nBodyAssembler_.nGrids();
- // const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
+ const auto& nBodyAssembler = contactNetwork_.nBodyAssembler();
+ const size_t nBodies = nBodyAssembler.nGrids();
+ // const auto& contactCouplings = nBodyAssembler.getContactCouplings();
size_t globalIdx = 0;
v_rel.resize(nBodies);
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index c7c74189131f7cb053e85c223bdafe73639a4c92..c4d59c3dfdc1ff88ca850dba527077be7a65dbcf 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -23,7 +23,7 @@ struct FixedPointIterationCounter {
std::ostream &operator<<(std::ostream &stream,
FixedPointIterationCounter const &fpic);
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
class FixedPointIterator {
using Functional = typename Factory::Functional;
using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
@@ -46,22 +46,20 @@ class FixedPointIterator {
public:
FixedPointIterator(const Dune::ParameterTree& parset,
- NBodyAssembler& nBodyAssembler,
+ const ContactNetwork& contactNetwork,
const IgnoreVector& ignoreNodes,
GlobalFriction& globalFriction,
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
const ErrorNorms& errorNorms);
- template <class Preconditioner>
FixedPointIterationCounter run(Updaters updaters,
- Preconditioner&& preconditioner,
const std::vector<Matrix>& velocityMatrices,
const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates);
private:
const Dune::ParameterTree& parset_;
- NBodyAssembler& nBodyAssembler_;
+ const ContactNetwork& contactNetwork_;
const IgnoreVector& ignoreNodes_;
GlobalFriction& globalFriction_;
diff --git a/src/spatial-solving/fixedpointiterator_tmpl.cc b/src/spatial-solving/fixedpointiterator_tmpl.cc
index 33acfce6a619c1be6a75b7c590503fb0b7beae44..5f4b30cff6d376cb6550cb2388b1b176a61df764 100644
--- a/src/spatial-solving/fixedpointiterator_tmpl.cc
+++ b/src/spatial-solving/fixedpointiterator_tmpl.cc
@@ -14,7 +14,6 @@
#include "../time-stepping/state/stateupdater.hh"
#include "../time-stepping/updaters.hh"
-using NBodyAssembler = typename MyContactNetwork::NBodyAssembler;
using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
@@ -25,4 +24,4 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
-template class FixedPointIterator<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
+template class FixedPointIterator<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;
diff --git a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
index bd5ff1a28e970d3c38ac1c2565fe601cba7813cc..93612ee520b8a60f0efd1d309b94c11f93c75133 100644
--- a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
@@ -13,6 +13,7 @@
#include "../../data-structures/levelcontactnetwork.hh"
#include "supportpatchfactory.hh"
+#include "../tnnmg/localproblem.hh"
#include <dune/localfunctions/lagrange/pqkfactory.hh>
@@ -108,6 +109,8 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
for (size_t bodyIdx=0; bodyIdx<levelContactNetwork_.nBodies(); bodyIdx++) {
const auto& gridView = levelContactNetwork_.body(bodyIdx)->gridView();
+ //printDofLocation(gridView);
+
for (const auto& e : elements(gridView)) {
const auto& refElement = Dune::ReferenceElements<double, dim>::general(e.type());
@@ -117,6 +120,27 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
if (!vertexVisited[globalIdx][0]) {
vertexVisited[globalIdx][0] = true;
patchFactory_.build(bodyIdx, e, i, patches_[globalIdx], patchDepth_);
+
+ /*print(patches_[globalIdx], "patch:");
+
+ size_t c = 0;
+ for (size_t j=0; j<levelContactNetwork_.nBodies(); j++) {
+ const auto& gv = fineContactNetwork_.body(j)->gridView();
+
+ printDofLocation(gv);
+
+ Dune::BlockVector<Dune::FieldVector<ctype, 1>> patchVec(gv.size(dim));
+ for (size_t l=0; l<patchVec.size(); l++) {
+ if (patches_[globalIdx][c++][0]) {
+ patchVec[l][0] = 1;
+ }
+ }
+
+ //print(patchVec, "patchVec");
+
+ // output patch
+ writeToVTK(gv, patchVec, "", "level_" + std::to_string(level_) + "_patch_" + std::to_string(globalIdx) + "_body_" + std::to_string(j));
+ }*/
}
}
}
@@ -160,6 +184,7 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
for (const auto& p : patches_) {
x = 0;
+ /*
auto& step = patchSolver_->getIterationStep();
dynamic_cast<LinearIterationStep<MatrixType, VectorType>&>(step).setProblem(*this->mat_, x, *this->rhs_);
step.setIgnore(p);
@@ -168,7 +193,29 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
patchSolver_->preprocess();
patchSolver_->solve();
- *(this->x_) += step.getSol();
+ *(this->x_) += step.getSol();*/
+
+ auto ignore = this->ignore();
+ for (size_t i=0; i<ignore.size(); i++) {
+ for (size_t d=0; d<dim; d++) {
+ if (p[i][d])
+ ignore[i][d] = true;
+ }
+ }
+
+ LocalProblem<MatrixType, VectorType> localProblem(*this->mat_, *this->rhs_, ignore);
+ Vector newR;
+ localProblem.getLocalRhs(x, newR);
+
+ /* print(ignore, "ignore:");
+ print(*this->mat_, "Mat:");
+ print(localProblem.getMat(), "localMat:");*/
+
+ patchSolver_->setProblem(localProblem.getMat(), x, newR);
+ patchSolver_->preprocess();
+ patchSolver_->solve();
+
+ *(this->x_) += x;
}
}
diff --git a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
index d16c2a652dffd03247141e5a308c7ad85ef639ab..68fa9b43f98bafdd3c6aa293f81981874b5f26b4 100644
--- a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
@@ -10,20 +10,27 @@
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/solvers/umfpacksolver.hh>
#include <dune/solvers/iterationsteps/blockgssteps.hh>
+//#include <dune/solvers/iterationsteps/blockgsstep.hh>
#include <dune/solvers/iterationsteps/cgstep.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
#include "nbodycontacttransfer.hh"
#include "levelpatchpreconditioner.hh"
+#include "../../utils/debugutils.hh"
+
template <class ContactNetwork, class MatrixType, class VectorType>
class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
public:
//struct Mode { enum ModeType {additive, multiplicative}; };
private:
- using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+ using Base = LinearIterationStep<MatrixType, VectorType>;
+
+ //using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+ using PatchSolver = Dune::Solvers::UMFPackSolver<MatrixType, VectorType>;
using LevelContactNetwork = typename ContactNetwork::LevelContactNetwork;
using LevelPatchPreconditioner = LevelPatchPreconditioner<LevelContactNetwork, PatchSolver, MatrixType, VectorType>;
@@ -44,10 +51,14 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
std::vector<VectorType> levelX_;
std::vector<VectorType> levelRhs_;
+ std::vector<BitVector> ignoreHierarchy_;
+
+ std::vector<std::shared_ptr<EnergyNorm<MatrixType, VectorType>>> levelErrorNorms_;
+ std::vector<std::shared_ptr<LinearIterationStep<MatrixType, VectorType>>> levelItSteps_;
std::vector<std::shared_ptr<PatchSolver>> levelSolvers_;
std::vector<BitVector> recompute_;
- std::vector<BitVector> truncation_;
+ //std::vector<BitVector> truncation_;
std::vector<std::shared_ptr<MGTransfer>> mgTransfer_; //TODO
public:
@@ -84,29 +95,36 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
}
}
+ levelMat_.resize(size()-1);
levelX_.resize(size());
levelRhs_.resize(size());
// init patch solvers
levelSolvers_.resize(size());
- auto coarseStep = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
- EnergyNorm<MatrixType, VectorType> coarseEnergyNorm(coarseStep);
- levelSolvers_[0] = std::make_shared<PatchSolver>(coarseStep,
+ levelItSteps_.resize(size());
+ levelErrorNorms_.resize(size());
+ //auto coarseStep =
+ levelItSteps_[0] = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::createPtr(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+ levelErrorNorms_[0] = std::make_shared<EnergyNorm<MatrixType, VectorType>>(*levelItSteps_[0].get());
+ levelSolvers_[0] = /*std::make_shared<PatchSolver>(*levelItSteps_[0].get(),
parset.get<size_t>("maximumIterations"),
parset.get<double>("tolerance"),
- coarseEnergyNorm,
+ *levelErrorNorms_[0].get(),
parset.get<Solver::VerbosityMode>("verbosity"),
- false); // absolute error
+ false); // absolute error */
+ std::make_shared<PatchSolver>();
for (size_t i=1; i<levelSolvers_.size(); i++) {
- auto gsStep = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
- EnergyNorm<MatrixType, VectorType> energyNorm(gsStep);
- levelSolvers_[i] = std::make_shared<PatchSolver>(gsStep,
+ //auto gsStep = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+ levelItSteps_[i] = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::createPtr(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+ levelErrorNorms_[i] = std::make_shared<EnergyNorm<MatrixType, VectorType>>(*levelItSteps_[i].get());
+ levelSolvers_[i] = /*std::make_shared<PatchSolver>(*levelItSteps_[i].get(),
parset.get<size_t>("maximumIterations"),
parset.get<double>("tolerance"),
- energyNorm,
+ *levelErrorNorms_[i].get(),
parset.get<Solver::VerbosityMode>("verbosity"),
- false); // absolute error
+ false); // absolute error */
+ std::make_shared<PatchSolver>();
levelPatchPreconditioners_[i]->setPatchSolver(levelSolvers_[i]);
}
@@ -118,37 +136,90 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
mgTransfer_[i]->setup(contactNetwork_, levelPatchPreconditioners_[i+1]->level(), levelPatchPreconditioners_[i+1]->fineLevel());
}
- //TODO
- //mgTransfer_.push_back(std::make_shared<MGTransfer>(coarseGlobalPreconditioner_->basis(), *allFaultBasis_));
-
// Remove any recompute filed so that initially the full transferoperator is assembled
for (size_t i=0; i<mgTransfer_.size(); i++)
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setRecomputeBitField(nullptr);
// specify the subset of entries to be reassembled at each iteration
- recompute_.resize(size());
+ /* recompute_.resize(size());
recompute_[recompute_.size()-1] = contactNetwork_.nBodyAssembler().totalHasObstacle_;
for (int i=mgTransfer_.size()-1; i>=0; i--) {
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrict(recompute_[i+1], recompute_[i]);
// std::dynamic_pointer_cast<DenseMultigridTransfer<VectorType, BitVector, MatrixType> >(mgTransfer_[i])->restrict(recompute_[i+1], recompute_[i]);
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setRecomputeBitField(&recompute_[i]);
- }
+ print(recompute_[i], "recompute: ");
+ }*/
}
+ void setMatrix(const MatrixType& mat) override {
+ Base::setMatrix(mat);
+
+ ignoreHierarchy_.resize(size());
+ ignoreHierarchy_.back() = Base::ignore();
+ for (int i=mgTransfer_.size()-1; i>=0; i--) {
+ std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrictToFathers(ignoreHierarchy_[i+1], ignoreHierarchy_[i]);
+ }
+
+ /* truncation_.resize(size());
+ truncation_[size()-1] = Base::ignore();
+ print(truncation_.back(), "truncation: ");
+ for (int i=mgTransfer_.size()-1; i>=0; i--) {
+ std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrict(truncation_[i+1], truncation_[i]);
+ std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setCriticalBitField(&truncation_[i+1]);
+ print(truncation_[i], "truncation: ");
+ }*/
+
+ size_t maxLevel = levelPatchPreconditioners_.size()-1;
+ levelPatchPreconditioners_[maxLevel]->setMatrix(mat);
+
+ //print(mat, "levelMat_" + std::to_string(maxLevel) + ":");
+
+ for (int i=levelPatchPreconditioners_.size()-2; i>0; i--) {
+ const auto& transfer = *mgTransfer_[i];
+
+ //print(*levelPatchPreconditioners_[i+1]->getMatrix(), "levelMat:");
+ //print(levelMat_[i], "coarseMat:");
- void setIgnore(const BitVector& ignore) {
- this->setIgnore(ignore);
+ transfer.galerkinRestrictSetOccupation(*levelPatchPreconditioners_[i+1]->getMatrix(), levelMat_[i]);
+ transfer.galerkinRestrict(*levelPatchPreconditioners_[i+1]->getMatrix(), levelMat_[i]);
+ levelPatchPreconditioners_[i]->setMatrix(levelMat_[i]);
+
+ //print(levelMat_[i], "levelMat_" + std::to_string(i) + ":");
+
+ // Set solution vector sizes for the lower levels
+ Dune::MatrixVector::resize(levelX_[i], levelMat_[i]);
+ //Dune::MatrixVector::resize(levelRhs_[i], levelMat_[i]);
+ }
+
+ // set coarse global problem
+ const auto& coarseTransfer = *mgTransfer_[0];
+ const auto& fineMat = *levelPatchPreconditioners_[1]->getMatrix();
+ coarseTransfer.galerkinRestrictSetOccupation(fineMat, levelMat_[0]);
+ coarseTransfer.galerkinRestrict(fineMat, levelMat_[0]);
+
+ //print(levelMat_[0], "levelMat_0:");
+
+ /*dynamic_cast<Base&>(levelSolvers_[0]->getIterationStep()).setMatrix(levelMat_[0]); */
+
+ Dune::MatrixVector::resize(levelX_[0], levelMat_[0]);
+ //Dune::MatrixVector::resize(levelRhs_[0], levelMat_[0]);
+ }
+
+ void setTruncation(const BitVector& truncation) {
+ std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_.back())->setCriticalBitField(&truncation);
+ }
+
+ /*void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) override {
+ LinearIterationStep<MatrixType, VectorType>::setProblem(mat, x, rhs);
truncation_.resize(size());
- truncation_[size()-1] = ignore;
+ truncation_[size()-1] = this->ignore();
+ print(truncation_.back(), "truncation: ");
for (int i=mgTransfer_.size()-1; i>=0; i--) {
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrict(truncation_[i+1], truncation_[i]);
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setCriticalBitField(&truncation_[i]);
+ print(truncation_[i], "truncation: ");
}
- }
-
- void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) override {
- this->setProblem(mat, x, rhs);
size_t maxLevel = levelPatchPreconditioners_.size()-1;
levelX_[maxLevel] = x;
@@ -176,10 +247,31 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
// EnergyNorm<MatrixType, VectorType> coarseErrorNorm(levelSolvers_[0]->getIterationStep());
// levelSolvers_[0]->setErrorNorm(coarseErrorNorm);
- }
-
+ } */
void iterate() {
+ size_t maxLevel = levelPatchPreconditioners_.size()-1;
+ levelX_[maxLevel] = *this->getIterate();
+ levelRhs_[maxLevel] = *Base::rhs_;
+
+ for (int i=levelPatchPreconditioners_.size()-2; i>0; i--) {
+ const auto& transfer = *mgTransfer_[i];
+
+ /*print(levelX_[i+1], "levelX_:");
+ print(levelX_[i], "levelX_:");
+ print(levelRhs_[i+1], "levelRhs_:");
+ print(levelRhs_[i], "levelRhs_:");*/
+
+ transfer.restrict(levelX_[i+1], levelX_[i]);
+ transfer.restrict(levelRhs_[i+1], levelRhs_[i]);
+ }
+
+ // set coarse global problem
+ const auto& coarseTransfer = *mgTransfer_[0];
+ coarseTransfer.restrict(levelX_[1], levelX_[0]);
+ coarseTransfer.restrict(levelRhs_[1], levelRhs_[0]);
+
+
if (mode_ == additive)
iterateAdd();
else
@@ -187,26 +279,46 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
}
void iterateAdd() {
+ //print(this->ignore(), "ignore: ");
+
*(this->x_) = 0;
VectorType x;
+ /*
// solve coarse global problem
auto& coarseSolver = levelSolvers_[0];
- coarseSolver->getIterationStep().setIgnore(truncation_[0]);
+ //coarseSolver->getIterationStep().setIgnore(ignoreHierarchy_[0]);
+
+ dynamic_cast<LinearIterationStep<MatrixType, VectorType>&>(levelSolvers_[0]->getIterationStep()).setProblem(levelMat_[0], levelX_[0], levelRhs_[0]);
coarseSolver->check();
coarseSolver->preprocess();
coarseSolver->solve();
- mgTransfer_[0]->prolong(coarseSolver->getIterationStep().getSol(), x);
+ mgTransfer_[0]->prolong(coarseSolver->getIterationStep().getSol(), x); */
- for (size_t i=1; i<levelPatchPreconditioners_.size(); i++) {
+
+ BitVector emptyIgnore0(levelX_[0].size(), false);
+ LocalProblem<MatrixType, VectorType> localProblem(levelMat_[0], levelRhs_[0], ignoreHierarchy_[0]);
+ Vector newR;
+ localProblem.getLocalRhs(levelX_[0], newR);
+
+ auto& coarseSolver = levelSolvers_[0];
+ coarseSolver->setProblem(localProblem.getMat(), levelX_[0], newR);
+ coarseSolver->preprocess();
+ coarseSolver->solve();
+
+ mgTransfer_[0]->prolong(levelX_[0], x);
+
+ //print(x, "MultilevelPatchPreconditioner: x0");
+
+ for (size_t i=1; i<levelPatchPreconditioners_.size()-1; i++) {
const auto& transfer = *mgTransfer_[i];
auto& preconditioner = *levelPatchPreconditioners_[i];
- preconditioner.setIgnore(truncation_[i]);
- preconditioner.iterate();
+ preconditioner.setIgnore(ignoreHierarchy_[i]);
+ preconditioner.apply(levelX_[i], levelRhs_[i]);
x += preconditioner.getSol();
@@ -214,8 +326,30 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
transfer.prolong(x, newX);
x = newX;
+ //print(x, "MultilevelPatchPreconditioner: x" + std::to_string(i));
}
+ auto& preconditioner = *levelPatchPreconditioners_.back();
+
+ preconditioner.setIgnore(ignoreHierarchy_.back());
+ preconditioner.apply(levelX_.back(), levelRhs_.back());
+
+ x += preconditioner.getSol();
+
+ //print(x, "MultilevelPatchPreconditioner: xmax");
+
+
+ BitVector emptyIgnore(this->x_->size(), false);
+ LocalProblem<MatrixType, VectorType> globalProblem(*this->mat_.get(), *this->rhs_, emptyIgnore);
+ globalProblem.getLocalRhs(*this->x_, newR);
+ VectorType totalX = *(this->x_);
+ auto& globalSolver = levelSolvers_.back();
+ globalSolver->setProblem(globalProblem.getMat(), totalX, newR);
+ globalSolver->preprocess();
+ globalSolver->solve();
+
+ //print(totalX, "MultilevelPatchPreconditioner: totalX");
+
*(this->x_) = x;
}
diff --git a/src/spatial-solving/preconditioners/supportpatchfactory.hh b/src/spatial-solving/preconditioners/supportpatchfactory.hh
index 153480970d3ad79f6336703c6ea365205d3d05a4..761e1eae193e8d14a2c01fb3eac9f25d2f8179bb 100644
--- a/src/spatial-solving/preconditioners/supportpatchfactory.hh
+++ b/src/spatial-solving/preconditioners/supportpatchfactory.hh
@@ -538,7 +538,7 @@ class SupportPatchFactory
std::set<size_t> isDofs;
intersectionDofs(fineIndices_, is, bodyID, isDofs);
- addToPatchBoundary(isDofs, patchBoundary, patchDofs);
+ addToPatch(isDofs, patchBoundary, patchDofs);
}
}
}
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index 77c31c1ce766266a8d833fbc8976037df6c90403..6b9a3119d186cdddee89dcb83f469a9a5b06fd55 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -20,6 +20,9 @@ SolverFactory<Functional, BitVector>::SolverFactory(
nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
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(), LineSearchSolver());
tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));
tnnmgStep_->setIgnore(ignoreNodes);
diff --git a/src/spatial-solving/tnnmg/functional.hh b/src/spatial-solving/tnnmg/functional.hh
index 761359448bfcc5fd192fd4b6c1ff57f2545c52af..4bad3e3feb7f37f5ea4d9073f7ce6b3269e070c0 100644
--- a/src/spatial-solving/tnnmg/functional.hh
+++ b/src/spatial-solving/tnnmg/functional.hh
@@ -97,8 +97,9 @@ class DirectionalRestriction :
{
phi_.directionalDomain(origin_, direction_, domain_);
- std::cout << domain_[0] << " " << domain_[1] << "Phi domain:" << std::endl;
- std::cout << this->defectLower_ << " " << this->defectUpper_ << "defect obstacles:" << std::endl;
+ //std::cout << domain_[0] << " " << domain_[1] << "Phi domain:" << std::endl;
+ //std::cout << this->defectLower_ << " " << this->defectUpper_ << "defect obstacles:" << std::endl;
+
if (domain_[0] < this->defectLower_) {
domain_[0] = this->defectLower_;
}
@@ -401,7 +402,7 @@ 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 Base::Base::operator()(v) + phi_.get().operator()(v);
+ return Dune::TNNMG::QuadraticFunctional<M,V,R>::operator()(v) + phi_.get().operator()(v);
return std::numeric_limits<Range>::max();
}
diff --git a/src/spatial-solving/tnnmg/linearcorrection.hh b/src/spatial-solving/tnnmg/linearcorrection.hh
index 2f99d1adfc6fb1a6d80613add41d39659facfd09..2b1d575232a2646a0f98d49583f7f4e16d8d97f2 100644
--- a/src/spatial-solving/tnnmg/linearcorrection.hh
+++ b/src/spatial-solving/tnnmg/linearcorrection.hh
@@ -9,6 +9,10 @@
#include <dune/solvers/solvers/linearsolver.hh>
#include <dune/solvers/common/canignore.hh>
+#include <dune/solvers/common/resize.hh>
+#include <dune/solvers/solvers/umfpacksolver.hh>
+#include "localproblem.hh"
+
/**
* \brief linear correction step for use by \ref TNNMGStep
*
@@ -44,7 +48,25 @@ LinearCorrection<Matrix, Vector>
makeLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > iterativeSolver)
{
return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
- using LinearIterationStep = Dune::Solvers::LinearIterationStep<Matrix, Vector>;
+
+ // compute reference solution directly
+ LocalProblem<Matrix, Vector> localProblem(A, b, ignore);
+ Vector newR, directX;
+ Dune::Solvers::resizeInitializeZero(directX, b);
+ localProblem.getLocalRhs(directX, newR);
+
+ /*print(*this->ignoreNodes_, "ignoreNodes:");
+ print(A, "A:");
+ print(localProblem.getMat(), "localMat:");*/
+
+ auto directSolver = std::make_shared<Dune::Solvers::UMFPackSolver<Matrix, Vector>>();
+
+ directSolver->setProblem(localProblem.getMat(), directX, newR);
+ directSolver->preprocess();
+ directSolver->solve();
+
+ x = directX;
+ /* using LinearIterationStep = Dune::Solvers::LinearIterationStep<Matrix, Vector>;
auto linearIterationStep = dynamic_cast<LinearIterationStep*>(&iterativeSolver->getIterationStep());
if (not linearIterationStep)
@@ -53,7 +75,13 @@ makeLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > i
linearIterationStep->setIgnore(ignore);
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;*/
};
}
diff --git a/src/spatial-solving/tnnmg/linearization.hh b/src/spatial-solving/tnnmg/linearization.hh
index 8dcecd1513f664bcb0ce8e7a68d7267eb92fc614..bc6095a8cb03395cfa236332a269f419962ea203 100644
--- a/src/spatial-solving/tnnmg/linearization.hh
+++ b/src/spatial-solving/tnnmg/linearization.hh
@@ -145,10 +145,8 @@ class Linearization {
negativeGradient_ *= -1;
// truncate gradient and hessian
- // truncateVector(negativeGradient_, truncationFlags_);
- // truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
-
- // the actual truncation is done by the solver of linearized problem
+ truncateVector(negativeGradient_, truncationFlags_);
+ truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
}
diff --git a/src/spatial-solving/tnnmg/tnnmgstep.hh b/src/spatial-solving/tnnmg/tnnmgstep.hh
index dc2ca30d6d8da342a664bc76ec141f0675ac6bb9..84d7e9147aca29bee78359c1558b34a245ade778 100644
--- a/src/spatial-solving/tnnmg/tnnmgstep.hh
+++ b/src/spatial-solving/tnnmg/tnnmgstep.hh
@@ -125,6 +125,15 @@ class TNNMGStep :
preSmoothingSteps_ = i;
}
+ auto energy(const Vector& v) const {
+ Vector temp;
+ Dune::Solvers::resizeInitializeZero(temp, v);
+ f_->quadraticPart().umv(v, temp);
+ temp *= 0.5;
+ temp -= f_->linearPart();
+ return temp*v;
+ }
+
/**
* \brief Do one TNNMG step
*/
@@ -141,12 +150,14 @@ class TNNMGStep :
//print(f.linearPart(), "f.linearPart():");
+ //std::cout << "-- energy before smoothing: " << energy(x) << std::endl;
+
// 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:");
// Compute constraint/truncated linearization
if (not linearization_)
@@ -157,48 +168,46 @@ class TNNMGStep :
auto&& A = linearization_->hessian();
auto&& r = linearization_->negativeGradient();
- //print(A, "TNNMG Linear problem matrix:");
- //print(r, "TNNMG Linear problem rhs:");
+ /*print(A, "TNNMG Linear problem matrix:");
+ print(r, "TNNMG Linear problem rhs:");
+ print(ignore, "ignore:");
+ print(linearization_->truncated(), "truncation:");*/
// Compute inexact solution of the linearized problem
Solvers::resizeInitializeZero(correction_, x);
Solvers::resizeInitializeZero(constrainedCorrection_, r);
- // std::cout << "- computing linear correction..." << std::endl;
+ //std::cout << "-- energy after presmoothing: " << energy(x) << std::endl;
- linearCorrection_(A, constrainedCorrection_, r, linearization_->truncated());
+ //print(constrainedCorrection_compare, "direct solver solution:");
- LocalProblem<ConstrainedMatrix, ConstrainedVector> localProblem(A, r, linearization_->truncated());
- Vector newR, constrainedCorrection_compare;
- Solvers::resizeInitializeZero(constrainedCorrection_compare, r);
- localProblem.getLocalRhs(constrainedCorrection_compare, newR);
+ //std::cout << "- computing linear correction..." << std::endl;
- /*print(*this->ignoreNodes_, "ignoreNodes:");
- print(A, "A:");
- print(localProblem.getMat(), "localMat:");*/
- auto directSolver = std::make_shared<Dune::Solvers::UMFPackSolver<ConstrainedMatrix, ConstrainedVector>>();
+ linearCorrection_(A, constrainedCorrection_, r, ignore);
- directSolver->setProblem(localProblem.getMat(), constrainedCorrection_compare, newR);
- directSolver->preprocess();
- directSolver->solve();
+ //DUNE_THROW(Dune::Exception, "TNNMGStep: Just need to terminate here!");
// linearCorrection_(A, constrainedCorrection_, r);
- std::cout << "- linear correction: success" << std::endl;
+ //std::cout << "- linear correction: success" << std::endl;
//print(constrainedCorrection_, "contrained correction:");
linearization_->extendCorrection(constrainedCorrection_, correction_);
- std::cout << "- extended correction: success" << std::endl;
- print(correction_, "correction:");
+ /*Vector h = x;
+ h += correction_;
+ std::cout << "-- energy after linear correction: " << energy(h) << std::endl;*/
+
+ //std::cout << "- extended correction: success" << std::endl;
+ //print(correction_, "correction:");
// Project onto admissible set
projection_(f, x, correction_);
- std::cout << "- projection onto admissible set: success" << std::endl;
+ //std::cout << "- projection onto admissible set: success" << std::endl;
- print(correction_, "correction:");
+ //print(correction_, "correction:");
double const correctionNorm = correction_.two_norm();
correction_ /= correctionNorm;
@@ -221,7 +230,11 @@ class TNNMGStep :
dampingFactor_ = 0;
correction_ *= dampingFactor_;
+ //std::cout << "damping factor: " << dampingFactor_ << std::endl;
+
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 fae06ed6b14302b5f4edbf2e41cca0dda6e4e337..c1869aa83726cbddc3f22a4fa27e52e61d873154 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -17,10 +17,10 @@ void IterationRegister::reset() {
finalCount = FixedPointIterationCounter();
}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::AdaptiveTimeStepper(
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::AdaptiveTimeStepper(
Dune::ParameterTree const &parset,
- NBodyAssembler& nBodyAssembler,
+ const ContactNetwork& contactNetwork,
const IgnoreVector& ignoreNodes,
GlobalFriction& globalFriction,
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
@@ -34,7 +34,7 @@ AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::AdaptiveTime
relativeTau_(relativeTau),
finalTime_(parset.get<double>("problem.finalTime")),
parset_(parset),
- nBodyAssembler_(nBodyAssembler),
+ contactNetwork_(contactNetwork),
ignoreNodes_(ignoreNodes),
globalFriction_(globalFriction),
bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
@@ -44,14 +44,13 @@ AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::AdaptiveTime
mustRefine_(mustRefine),
errorNorms_(errorNorms) {}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-bool AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::reachedEnd() {
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+bool AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::reachedEnd() {
return relativeTime_ >= 1.0;
}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-template <class Preconditioner>
-IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::advance(Preconditioner&& preconditioner) {
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::advance() {
/*
| C | We check here if making the step R1 of size tau is a
| R1 | R2 | good idea. To check if we can coarsen, we compare
@@ -64,7 +63,7 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
std::cout << "AdaptiveTimeStepper::advance()" << std::endl;
if (R1_.updaters == Updaters())
- R1_ = step(current_, relativeTime_, relativeTau_, std::forward<Preconditioner>(preconditioner));
+ R1_ = step(current_, relativeTime_, relativeTau_);
std::cout << "AdaptiveTimeStepper Step 1" << std::endl;
@@ -74,9 +73,9 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
UpdatersWithCount R2;
UpdatersWithCount C;
while (relativeTime_ + relativeTau_ <= 1.0) {
- R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_, std::forward<Preconditioner>(preconditioner));
+ R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_);
std::cout << "AdaptiveTimeStepper R2 computed!" << std::endl << std::endl;
- C = step(current_, relativeTime_, 2 * relativeTau_, std::forward<Preconditioner>(preconditioner));
+ C = step(current_, relativeTime_, 2 * relativeTau_);
std::cout << "AdaptiveTimeStepper C computed!" << std::endl << std::endl;
if (mustRefine_(R2.updaters, C.updaters))
break;
@@ -91,10 +90,10 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
UpdatersWithCount F2;
if (!didCoarsen) {
while (true) {
- F1 = step(current_, relativeTime_, relativeTau_ / 2.0, std::forward<Preconditioner>(preconditioner));
+ F1 = step(current_, relativeTime_, relativeTau_ / 2.0);
std::cout << "AdaptiveTimeStepper F1 computed!" << std::endl << std::endl;
F2 = step(F1.updaters, relativeTime_ + relativeTau_ / 2.0,
- relativeTau_ / 2.0, std::forward<Preconditioner>(preconditioner));
+ relativeTau_ / 2.0);
std::cout << "AdaptiveTimeStepper F2 computed!" << std::endl << std::endl;
if (!mustRefine_(F2.updaters, R1_.updaters)) {
std::cout << "Sufficiently refined!" << std::endl;
@@ -122,17 +121,16 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
return iterationRegister_;
}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-template <class Preconditioner>
-typename AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::UpdatersWithCount
-AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::step(
- Updaters const &oldUpdaters, double rTime, double rTau, Preconditioner&& preconditioner) {
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+typename AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::UpdatersWithCount
+AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::step(
+ Updaters const &oldUpdaters, double rTime, double rTau) {
UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
newUpdatersAndCount.count =
- MyCoupledTimeStepper(finalTime_, parset_, nBodyAssembler_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_,
+ MyCoupledTimeStepper(finalTime_, parset_, contactNetwork_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_,
newUpdatersAndCount.updaters, errorNorms_,
externalForces_)
- .step(rTime, rTau, std::forward<Preconditioner>(preconditioner));
+ .step(rTime, rTau);
iterationRegister_.registerCount(newUpdatersAndCount.count);
return newUpdatersAndCount;
}
diff --git a/src/time-stepping/adaptivetimestepper.hh b/src/time-stepping/adaptivetimestepper.hh
index 3e48b2218a06a954903e9c2d4d63d0fe0642c18c..ce090c012658d23ff9b34dd9f182ffe2421df013 100644
--- a/src/time-stepping/adaptivetimestepper.hh
+++ b/src/time-stepping/adaptivetimestepper.hh
@@ -15,7 +15,7 @@ struct IterationRegister {
FixedPointIterationCounter finalCount;
};
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
class AdaptiveTimeStepper {
struct UpdatersWithCount {
Updaters updaters;
@@ -27,7 +27,7 @@ class AdaptiveTimeStepper {
//using ConvexProblem = typename Factory::ConvexProblem;
//using Nonlinearity = typename Factory::Nonlinearity;
- using MyCoupledTimeStepper = CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>;
+ using MyCoupledTimeStepper = CoupledTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>;
using GlobalFriction = typename MyCoupledTimeStepper::GlobalFriction;
using BitVector = typename MyCoupledTimeStepper::BitVector;
@@ -36,7 +36,7 @@ class AdaptiveTimeStepper {
public:
AdaptiveTimeStepper(
Dune::ParameterTree const &parset,
- NBodyAssembler& nBodyAssembler,
+ const ContactNetwork& contactNetwork,
const IgnoreVector& ignoreNodes,
GlobalFriction& globalFriction,
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
@@ -49,20 +49,18 @@ class AdaptiveTimeStepper {
bool reachedEnd();
- template <class Preconditioner>
- IterationRegister advance(Preconditioner&& preconditioner);
+ IterationRegister advance();
double relativeTime_;
double relativeTau_;
private:
- template <class Preconditioner>
UpdatersWithCount step(Updaters const &oldUpdaters, double rTime,
- double rTau, Preconditioner&& preconditioner);
+ double rTau);
double finalTime_;
Dune::ParameterTree const &parset_;
- NBodyAssembler& nBodyAssembler_;
+ const ContactNetwork& contactNetwork_;
const IgnoreVector& ignoreNodes_;
GlobalFriction& globalFriction_;
diff --git a/src/time-stepping/adaptivetimestepper_tmpl.cc b/src/time-stepping/adaptivetimestepper_tmpl.cc
index bda3c4c774b1db0d13f12e4004465a6df0fa00f7..ff76be34499f936cd42839472ecc7207d81ebea5 100644
--- a/src/time-stepping/adaptivetimestepper_tmpl.cc
+++ b/src/time-stepping/adaptivetimestepper_tmpl.cc
@@ -14,7 +14,6 @@
#include "state/stateupdater.hh"
#include "updaters.hh"
-using NBodyAssembler = typename MyContactNetwork::NBodyAssembler;
using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
@@ -24,6 +23,6 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
-using MyAdaptiveTimeStepper = AdaptiveTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
+using MyAdaptiveTimeStepper = AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;
-template class AdaptiveTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
+template class AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;
diff --git a/src/time-stepping/coupledtimestepper.cc b/src/time-stepping/coupledtimestepper.cc
index 6b9daac7bef9458d91fff57ca8da2e68c775b9bb..1df547741bee8de1cfb4e2b339f77b72badc8308 100644
--- a/src/time-stepping/coupledtimestepper.cc
+++ b/src/time-stepping/coupledtimestepper.cc
@@ -4,10 +4,10 @@
#include "coupledtimestepper.hh"
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::CoupledTimeStepper(
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+CoupledTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::CoupledTimeStepper(
double finalTime, Dune::ParameterTree const &parset,
- NBodyAssembler& nBodyAssembler,
+ const ContactNetwork& contactNetwork,
const IgnoreVector& ignoreNodes,
GlobalFriction& globalFriction,
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
@@ -16,7 +16,7 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::CoupledTimeSt
ExternalForces& externalForces)
: finalTime_(finalTime),
parset_(parset),
- nBodyAssembler_(nBodyAssembler),
+ contactNetwork_(contactNetwork),
ignoreNodes_(ignoreNodes),
globalFriction_(globalFriction),
bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
@@ -24,11 +24,10 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::CoupledTimeSt
externalForces_(externalForces),
errorNorms_(errorNorms) {}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-template <class Preconditioner>
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
FixedPointIterationCounter
-CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::step(double relativeTime,
- double relativeTau, Preconditioner&& preconditioner) {
+CoupledTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::step(double relativeTime,
+ double relativeTau) {
std::cout << "CoupledTimeStepper::step()" << std::endl;
@@ -50,8 +49,8 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::step(double r
updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS, velocityIterate, velocityMatrix);
FixedPointIterator fixedPointIterator(
- parset_, nBodyAssembler_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_, errorNorms_);
- auto const iterations = fixedPointIterator.run(updaters_, std::forward<Preconditioner>(preconditioner),
+ parset_, contactNetwork_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_, errorNorms_);
+ auto const iterations = fixedPointIterator.run(updaters_,
velocityMatrix, velocityRHS, velocityIterate);
return iterations;
}
diff --git a/src/time-stepping/coupledtimestepper.hh b/src/time-stepping/coupledtimestepper.hh
index 0b2ac4689e334de1d14ff78edde65d3186c23e05..fc619036f26c54137f0d7a2c48d53f4fb593d6d4 100644
--- a/src/time-stepping/coupledtimestepper.hh
+++ b/src/time-stepping/coupledtimestepper.hh
@@ -8,12 +8,12 @@
#include "../spatial-solving/fixedpointiterator.hh"
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
class CoupledTimeStepper {
using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix;
using IgnoreVector = typename Factory::BitVector;
- using FixedPointIterator = FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>;
+ using FixedPointIterator = FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>;
public:
using GlobalFriction = typename FixedPointIterator::GlobalFriction;
@@ -23,7 +23,7 @@ class CoupledTimeStepper {
public:
CoupledTimeStepper(double finalTime,
Dune::ParameterTree const &parset,
- NBodyAssembler& nBodyAssembler,
+ const ContactNetwork& contactNetwork,
const IgnoreVector& ignoreNodes,
GlobalFriction& globalFriction,
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
@@ -31,13 +31,12 @@ class CoupledTimeStepper {
const ErrorNorms& errorNorms,
ExternalForces& externalForces);
- template <class Preconditioner>
- FixedPointIterationCounter step(double relativeTime, double relativeTau, Preconditioner&& preconditioner);
+ FixedPointIterationCounter step(double relativeTime, double relativeTau);
private:
double finalTime_;
Dune::ParameterTree const &parset_;
- NBodyAssembler& nBodyAssembler_;
+ const ContactNetwork& contactNetwork_;
const IgnoreVector& ignoreNodes_;
GlobalFriction& globalFriction_;
diff --git a/src/time-stepping/coupledtimestepper_tmpl.cc b/src/time-stepping/coupledtimestepper_tmpl.cc
index 4d9f7f8eb062de3f2731ad3c7205631fa9799ef7..3df79591219c4b017965e642c088cc922b676025 100644
--- a/src/time-stepping/coupledtimestepper_tmpl.cc
+++ b/src/time-stepping/coupledtimestepper_tmpl.cc
@@ -14,7 +14,6 @@
#include "state/stateupdater.hh"
#include "updaters.hh"
-using NBodyAssembler = typename MyContactNetwork::NBodyAssembler;
using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
@@ -25,4 +24,4 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
-template class CoupledTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
+template class CoupledTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;
diff --git a/src/time-stepping/rate/newmark.cc b/src/time-stepping/rate/newmark.cc
index a3b3638529c8a75a09f633a5ea6a26b9ef6eb867..dec66b98d732c2589b7bc79c4c371a47de336714 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>