.

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

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];
       };

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]

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,

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);