diff --git a/src/matrices.hh b/src/matrices.hh
index 4fcaabfb452845a158793dc50739612e64c60988..01fc682dcbaff5b21b88a03f8bb958faa3c52be2 100644
--- a/src/matrices.hh
+++ b/src/matrices.hh
@@ -8,9 +8,15 @@ template <class Matrix> class Matrices {
   std::vector<std::shared_ptr<Matrix>> mass;
 
   Matrices(size_t n) {
-      elasticity.resize(n);
-      damping.resize(n);
-      mass.resize(n);
+    elasticity.resize(n);
+    damping.resize(n);
+    mass.resize(n);
+
+      for (size_t i=0; i<n; i++) {
+        elasticity[i] = std::make_shared<Matrix>();
+        damping[i] = std::make_shared<Matrix>();
+        mass[i] = std::make_shared<Matrix>();
+      }
   }
 };
 #endif
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index e1183e14dff399f7b52e25134b7024f79c5c01f3..5e3e54231ce4383fb6700c1fa16800cfd2128c8b 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -237,6 +237,8 @@ int main(int argc, char *argv[]) {
     std::vector<std::shared_ptr<CouplingPair>> couplings(bodyCount-1);
 
     for (size_t i=0; i<couplings.size(); i++) {
+      couplings[i] = std::make_shared<CouplingPair>();
+
       auto nonmortarPatch = std::make_shared<LevelBoundaryPatch>(levelFaces[i]->upper);
       auto mortarPatch = std::make_shared<LevelBoundaryPatch>(levelFaces[i+1]->lower);
 
@@ -285,8 +287,7 @@ int main(int argc, char *argv[]) {
       std::cout << "Grid" << i << " Number of DOFs: " << leafVertexCount << std::endl;
 
       // Neumann boundary
-      auto& neumannBoundary = neumannBoundaries[i];
-      neumannBoundary = BoundaryPatch<DefLeafGridView>(*leafViews[i]);
+      neumannBoundaries[i] = BoundaryPatch<DefLeafGridView>(*leafViews[i]);
 
       // friction boundary
       auto& frictionalBoundary = frictionalBoundaries[i];
@@ -348,7 +349,7 @@ int main(int argc, char *argv[]) {
       // Problem formulation: right-hand side
       externalForces[i] =
             [&](double _relativeTime, Vector &_ell) {
-              assembler->assembleNeumann(neumannBoundaries[i], _ell, neumannFunction,
+              assemblers[i]->assembleNeumann(neumannBoundaries[i], _ell, neumannFunction,
                                           _relativeTime);
               _ell += gravityFunctionals[i];
       };
diff --git a/src/multi-body-problem.cfg b/src/multi-body-problem.cfg
index b8358f558c5a1b5fbc853aeec8e7c6d2e17878c2..b850b5cc829734e50342d1e812254778f74ad97b 100644
--- a/src/multi-body-problem.cfg
+++ b/src/multi-body-problem.cfg
@@ -2,11 +2,11 @@
 gravity         = 9.81  # [m/s^2]
 
 [io]
-data.write      = true
+data.write      = false #true
 printProgress   = false
 restarts.first  = 0
 restarts.spacing= 20
-restarts.write  = true
+restarts.write  = false #true
 vtk.write       = false
 
 [problem]
diff --git a/src/program_state.hh b/src/program_state.hh
index 68c465b472a0cc5bee3d28f9b4b32b27270b8bfc..dfacce35338ca3ab701274a3ff8d2c80619d35d3 100644
--- a/src/program_state.hh
+++ b/src/program_state.hh
@@ -101,7 +101,6 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
         const std::vector<Vector>& _rhs, std::vector<Vector>& _x,
         Dune::ParameterTree const &_localParset) {
 
-
       std::vector<const Matrix*> matrices_ptr(_matrices.size());
       for (size_t i=0; i<matrices_ptr.size(); i++) {
             matrices_ptr[i] = _matrices[i].get();
@@ -159,7 +158,9 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       v[i] = 0.0;
 
       ell0[i].resize(u[i].size());
-      externalForces[i](relativeTime, ell0[i]);
+      ell0[i] = 0.0;
+      // TODO
+      //externalForces[i](relativeTime, ell0[i]);
     }
 
     // Initial displacement: Start from a situation of minimal stress,
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index faf14f1c281a5509fee6c4d824c9e2ca7183fc70..f9220a826a7e5e9a5aae422ea8def92933421458 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -25,6 +25,7 @@ SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::SolverFactory(
                        Solver::QUIET),
       transferOperators_(nBodyAssembler.getGrids().at(0)->maxLevel()) {
 
+  multigridStep_ = std::make_shared<Step>();
   // tnnmg iteration step
   multigridStep_->setMGType(parset.get<int>("main.multi"), parset.get<int>("main.pre"), parset.get<int>("main.post"));
   //multigridStep_->setIgnore(ignoreNodes);