diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 9f2b6c6778f5c109dc5060ac8d9bdf9f83c9a5b7..14eb31c3b5f878521ecb1710752907678984b8ee 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -69,7 +69,7 @@
//#include "time-stepping/adaptivetimestepper.hh"
#include "time-stepping/rate.hh"
#include "time-stepping/state.hh"
-//#include "time-stepping/updaters.hh"
+#include "time-stepping/updaters.hh"
#include "vtk.hh"
// for getcwd
@@ -320,6 +320,7 @@ int main(int argc, char *argv[]) {
std::vector<Vector> gravityFunctionals(bodyCount);
std::vector<std::function<void(double, Vector&)>> externalForces(bodyCount);
+ std::vector<EnergyNorm<ScalarMatrix, ScalarVector>* > stateEnergyNorms(bodyCount);
for (size_t i=0; i<assemblers.size(); i++) {
auto& assembler = assemblers[i];
@@ -332,7 +333,7 @@ int main(int argc, char *argv[]) {
ScalarMatrix relativeFrictionalBoundaryMass;
assembler->assembleFrictionalBoundaryMass(frictionalBoundaries[i], relativeFrictionalBoundaryMass);
relativeFrictionalBoundaryMass /= frictionalBoundaries[i].area();
- EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(relativeFrictionalBoundaryMass);
+ stateEnergyNorms[i] = new EnergyNorm<ScalarMatrix, ScalarVector>(relativeFrictionalBoundaryMass);
// assemble forces
assembler->assembleBodyForce(body.getGravityField(), gravityFunctionals[i]);
@@ -477,13 +478,10 @@ int main(int argc, char *argv[]) {
report(true);
*/
-/*
+
// Set up TNNMG solver
- using NonlinearFactory = SolverFactory<
- dims,
- MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
- Grid>;
- NonlinearFactory factory(parset.sub("solver.tnnmg"), *grid, dirichletNodes);
+ using NonlinearFactory = SolverFactory<DeformedGrid, Matrix, Vector>;
+ NonlinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, dirichletNodes);
using MyUpdater = Updaters<RateUpdater<Vector, Matrix, Function, dims>,
StateUpdater<ScalarVector, Vector>>;
@@ -501,15 +499,20 @@ int main(int argc, char *argv[]) {
parset.get<double>("timeSteps.refinementTolerance");
auto const mustRefine = [&](MyUpdater &coarseUpdater,
MyUpdater &fineUpdater) {
- ScalarVector coarseAlpha;
+ std::vector<ScalarVector> coarseAlpha;
coarseUpdater.state_->extractAlpha(coarseAlpha);
- ScalarVector fineAlpha;
+ std::vector<ScalarVector> fineAlpha;
fineUpdater.state_->extractAlpha(fineAlpha);
- return stateEnergyNorm.diff(fineAlpha, coarseAlpha) > refinementTolerance;
+ ScalarVector::field_type energyNorm = 0;
+ for (size_t i=0; i<stateEnergyNorms.size(); i++) {
+ energyNorm += stateEnergyNorms[i]->diff(fineAlpha[i], coarseAlpha[i]);
+ }
+ return energyNorm > refinementTolerance;
};
+ /*
std::signal(SIGXCPU, handleSignal);
std::signal(SIGINT, handleSignal);
std::signal(SIGTERM, handleSignal);
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index 1d4756840593d9c537e373eb43f57a9a112b66d2..5fc74bfc06b80fb1aecfa4edd8b42d7879b2b472 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -50,8 +50,9 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
template <class Factory, class Updaters, class ErrorNorm>
FixedPointIterationCounter
FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
- Updaters updaters, Matrix const &velocityMatrix, Vector const &velocityRHS,
- Vector &velocityIterate) {
+ Updaters updaters, const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
+ std::vector<Vector>& velocityIterates) {
+
EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix);
LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
velocityTolerance_, &energyNorm,
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index e0fa8ec1968ab6ed398baaa9d521e18e54a93188..c68ba188e06f05eecdae9b1d639dc13527c38914 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -15,8 +15,8 @@
template <class DeformedGridTEMPLATE, class MatrixType, class VectorType>
SolverFactory<DeformedGrid, Matrix, Vector>::SolverFactory(
- Dune::ParameterTree const &parset, const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
- Dune::BitSetVector<dim> const &ignoreNodes)
+ const Dune::ParameterTree& parset, const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
+ const std::vector<Dune::BitSetVector<DeformedGrid::dimension>>& ignoreNodes)
: nBodyAssembler_(nBodyAssembler),
baseEnergyNorm_(baseSolverStep_),
baseSolver_(&baseSolverStep_,
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index ad6dbde284eda07e0d85364dfaf30beca03c257e..2c4bc74ed0c5ec9c1e63dcf50f3a0befe1a7377c 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -14,7 +14,7 @@
#include <dune/contact/assemblers/nbodyassembler.hh>
-#include <dune/contact/estimators/hierarchiccontactestimator.hh>
+//#include <dune/contact/estimators/hierarchiccontactestimator.hh>
#include <dune/contact/solvers/nsnewtonmgstep.hh>
#include <dune/contact/solvers/contacttransfer.hh>
#include <dune/contact/solvers/contactobsrestrict.hh>
@@ -25,8 +25,8 @@
template <class DeformedGridTEMPLATE, class MatrixType, class VectorType>
class SolverFactory {
-protected:
- const size_t dim = DeformedGrid::dimension;
+//protected:
+// const size_t dim = DeformedGrid::dimension;
public:
using Vector = VectorType;
@@ -38,7 +38,7 @@ class SolverFactory {
using Step = Dune::Contact::NonSmoothNewtonMGStep<Matrix, Vector>;
SolverFactory(Dune::ParameterTree const &parset, const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
- Dune::BitSetVector<dim> const &ignoreNodes);
+ const std::vector<Dune::BitSetVector<DeformedGrid::dimension>>& ignoreNodes);
~SolverFactory();
diff --git a/src/time-stepping/adaptivetimestepper_tmpl.cc b/src/time-stepping/adaptivetimestepper_tmpl.cc
index d3e8a9c082c3baa4bf46abaf176847096197ef3b..e369a8f5e2a77dfd50ff5ef85ca000fd84637832 100644
--- a/src/time-stepping/adaptivetimestepper_tmpl.cc
+++ b/src/time-stepping/adaptivetimestepper_tmpl.cc
@@ -19,10 +19,7 @@
#include "updaters.hh"
using Function = Dune::VirtualFunction<double, double>;
-using Factory = SolverFactory<
- MY_DIM,
- MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
- Grid>;
+using Factory = SolverFactory<DeformedGrid, Matrix, Vector>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
diff --git a/src/time-stepping/coupledtimestepper.cc b/src/time-stepping/coupledtimestepper.cc
index 3be5824daf2f02c27e89fad45a2c76fa90fd3dbd..b8a4129a891806ee3130814ceb3670998a7c5fa1 100644
--- a/src/time-stepping/coupledtimestepper.cc
+++ b/src/time-stepping/coupledtimestepper.cc
@@ -26,17 +26,17 @@ CoupledTimeStepper<Factory, Updaters, ErrorNorm>::step(double relativeTime,
updaters_.rate_->nextTimeStep();
auto const newRelativeTime = relativeTime + relativeTau;
- Vector ell;
+ std::vector<Vector> ell;
externalForces_(newRelativeTime, ell);
- Matrix velocityMatrix;
- Vector velocityRHS;
- Vector velocityIterate;
+ std::vector<Matrix> velocityMatrix;
+ std::vector<Vector> velocityRHS;
+ std::vector<Vector> velocityIterate;
auto const tau = relativeTau * finalTime_;
- updaters_.state_->setup(tau);
- updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS,
- velocityIterate, velocityMatrix);
+ updaters_.state_->setup(tau);
+ updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS, velocityIterate, velocityMatrix);
+
FixedPointIterator<Factory, Updaters, ErrorNorm> fixedPointIterator(
factory_, parset_, globalFriction_, errorNorm_);
auto const iterations = fixedPointIterator.run(updaters_, velocityMatrix,