.

src/multi-body-problem-2D.cfg

 # -*- mode:conf -*-
 [boundary.friction]
-smallestDiameter = 0.01  # 2e-3 [m]
+smallestDiameter = 0.05  # 2e-3 [m]
 
 [timeSteps]
 refinementTolerance = 1e-5 # 1e-5

src/multi-body-problem.cfg

@@ -47,6 +47,7 @@ relativeTau = 1e-4 # 1e-6
 
 [timeSteps]
 scheme = newmark
+timeSteps = 5
 
 [u0.solver]
 maximumIterations = 100
@@ -76,7 +77,7 @@ maximumIterations = 10000
 verbosity         = quiet
 
 [solver.tnnmg.main]
-pre   = 1
+pre   = 10
 multi = 5 # number of multigrid steps
 post  = 0
 

src/solverfactorytest.cc

@@ -57,8 +57,6 @@
 
 const int dim = MY_DIM;
 
-const int timeSteps = 5;
-
 Dune::ParameterTree parset;
 size_t bodyCount;
 
@@ -71,6 +69,7 @@ std::vector<ScalarVector> weightedNormalStress;
 double relativeTime;
 double relativeTau;
 size_t timeStep = 0;
+size_t timeSteps = 0;
 
 const std::string path = "";
 const std::string outputFile = "solverfactorytest.log";
@@ -78,7 +77,7 @@ const std::string outputFile = "solverfactorytest.log";
 template <class ContactNetwork>
 void solveProblem(const ContactNetwork& contactNetwork,
                   const Matrix& mat, const Vector& rhs, Vector& x,
-                  const BitVector& ignore, const Vector& lower, const Vector& upper) {
+                  const BitVector& ignore, const Vector& lower, const Vector& upper, bool initial = false) {
 
     using Solver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
     using Norm =  EnergyNorm<Matrix, Vector>;
@@ -96,6 +95,8 @@ void solveProblem(const ContactNetwork& contactNetwork,
     using ContactFunctional = Dune::TNNMG::BoxConstrainedQuadraticFunctional<Matrix&, Vector&, Vector&, Vector&, double>;
     auto I = ContactFunctional(mat, rhs, lower, upper);
 
+    std::cout << "Energy start iterate: " << I(x) << std::endl;
+
     using LocalSolver = Dune::TNNMG::ScalarObstacleSolver;
     auto localSolver = Dune::TNNMG::gaussSeidelLocalSolver(LocalSolver());
 
@@ -176,10 +177,16 @@ void solveProblem(const ContactNetwork& contactNetwork,
     refSolver.solve();
     std::cout << correctionNorms.size() << std::endl;
 
-
+    if (initial) {
+        x = refX;
+        return;
+    }
     // set up solver factory solver
 
     // set up functional
+    /*auto& globalFriction = contactNetwork.globalFriction();
+    using MyFunctional = Functional<Matrix&, Vector&, std::decay_t<decltype(globalFriction)>&, Vector&, Vector&, typename Matrix::field_type>;
+    MyFunctional J(mat, rhs, globalFriction, lower, upper);*/
     using MyFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
     MyFunctional J(mat, rhs, ZeroNonlinearity(), lower, upper);
 
@@ -216,7 +223,7 @@ void solveProblem(const ContactNetwork& contactNetwork,
     diff -= refX;
     std::cout << "Solver error in energy norm: " << norm(diff) << std::endl;
 
-    x = refX;
+    std::cout << "Energy end iterate: " << J(x) << std::endl;
 }
 
 
@@ -274,7 +281,7 @@ void setupInitialConditions(const ContactNetwork& contactNetwork) {
     print(lower, "lower");
     print(upper, "upper");*/
 
-    solveProblem(contactNetwork, bilinearForm, totalRhs, totalX, _dirichletNodes, lower, upper);
+    solveProblem(contactNetwork, bilinearForm, totalRhs, totalX, _dirichletNodes, lower, upper, true);
 
     nBodyAssembler.postprocess(totalX, _x);
   };
@@ -434,17 +441,17 @@ void run(Updaters& updaters, ContactNetwork& contactNetwork,
   }
 
   // compute velocity obstacles
-  Vector vLower, vUpper;
+  /*Vector vLower, vUpper;
   std::vector<Vector> u0, v0;
   updaters.rate_->extractOldVelocity(v0);
   updaters.rate_->extractOldDisplacement(u0);
 
   Vector totalu0, totalv0;
-  nBodyAssembler.concatenateVectors(u0, totalu0);
-  nBodyAssembler.concatenateVectors(v0, totalv0);
+  nBodyAssembler.nodalToTransformed(u0, totalu0);
+  nBodyAssembler.nodalToTransformed(v0, totalv0);
 
-  //updaters.rate_->velocityObstacles(totalu0, lower, totalv0, vLower);
-  //updaters.rate_->velocityObstacles(totalu0, upper, totalv0, vUpper);
+  updaters.rate_->velocityObstacles(totalu0, lower, totalv0, vLower);
+  updaters.rate_->velocityObstacles(totalu0, upper, totalv0, vUpper);*/
 
   const auto& errorNorms = contactNetwork.stateEnergyNorms();
 
@@ -459,22 +466,24 @@ void run(Updaters& updaters, ContactNetwork& contactNetwork,
 
   size_t fixedPointIteration;
   size_t multigridIterations = 0;
+
   std::vector<ScalarVector> alpha(nBodies);
   updaters.state_->extractAlpha(alpha);
+
+  Vector totalVelocityIterate;
+  nBodyAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
+
   for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations;
        ++fixedPointIteration) {
 
     // contribution from nonlinearity
     globalFriction.updateAlpha(alpha);
 
-    Vector totalVelocityIterate;
-    nBodyAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
-
     solveProblem(contactNetwork, bilinearForm, totalRhs, totalVelocityIterate, totalDirichletNodes, lower, upper);
 
     nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
 
-    std::vector<Vector> v_m;
+    std::vector<Vector> v_m; //TODO : wrong, isnt used atm;
     updaters.rate_->extractOldVelocity(v_m);
 
     for (size_t i=0; i<v_m.size(); i++) {
@@ -730,6 +739,7 @@ int main(int argc, char *argv[]) { try {
 
     //const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
 
+    timeSteps = parset.get<size_t>("timeSteps.timeSteps");
     for (timeStep=1; timeStep<=timeSteps; timeStep++) {
 
       advanceStep(current, contactNetwork);

src/spatial-solving/solverfactory.cc

@@ -20,14 +20,10 @@ SolverFactory<Functional, BitVector>::SolverFactory(
         J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J))) {
 
     auto localSolver = Dune::TNNMG::gaussSeidelLocalSolver(LocalSolver());
-    //nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, localSolver);
-
-    //using MyNonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, decltype(localSolver), BitVector>;
     nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, localSolver);
 
     //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());

src/spatial-solving/solverfactory.hh

@@ -33,7 +33,7 @@ class SolverFactory {
     using Linearization = Linearization<Functional, BitVector>;
     using DefectProjection = typename Dune::TNNMG::ObstacleDefectProjection;
 
-    //using Step = TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
+    //using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
     using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
 
   template <class LinearSolver>