From 7e679c0252d9c42eac19ecb6bd27b57d4275dc76 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Sat, 17 Mar 2018 17:02:05 +0100
Subject: [PATCH] .

---
 src/multi-body-problem.cc                     | 25 +++++++++++--------
 src/spatial-solving/fixedpointiterator.cc     |  5 ++--
 src/spatial-solving/solverfactory.cc          |  4 +--
 src/spatial-solving/solverfactory.hh          |  8 +++---
 src/time-stepping/adaptivetimestepper_tmpl.cc |  5 +---
 src/time-stepping/coupledtimestepper.cc       | 14 +++++------
 6 files changed, 31 insertions(+), 30 deletions(-)

diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 9f2b6c67..14eb31c3 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 1d475684..5fc74bfc 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 e0fa8ec1..c68ba188 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 ad6dbde2..2c4bc74e 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 d3e8a9c0..e369a8f5 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 3be5824d..b8a4129a 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,
-- 
GitLab