diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index 81be4536caf79567e61f0eeab7df9994518f183e..521c513b30c767170052a6c611cdbfe90ea19276 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -52,28 +52,22 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
 template <class Factory, class Updaters, class ErrorNorm>
 FixedPointIterationCounter
 FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
-    Updaters updaters, const std::vector<const Matrix*>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
+    Updaters updaters, const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
     std::vector<Vector>& velocityIterates) {
 
-  EnergyNorm<Matrix, Vector> energyNorm(velocityMatrices[0]);
-  LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
-                                           velocityTolerance_, &energyNorm,
-                                           verbosity_, false); // absolute error
-
-  // assemble full global contact problem
   auto contactAssembler = factory_.getNBodyAssembler();
 
-  Matrix bilinearForm;
-  contactAssembler.assembleJacobian(velocityMatrices, bilinearForm);
-
-  Vector totalRhs;
-  contactAssembler.assembleRightHandSide(velocityRHSs, totalRhs);
-
-  Vector totalVelocityIterate;
-  contactAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
+  Matrix mat;
+  std::vector<const Matrix*> matrices_ptr(velocityMatrices.size());
+  for (size_t i=0; i<matrices_ptr.size(); i++) {
+      matrices_ptr[i] = &velocityMatrices[i];
+  }
+  contactAssembler.assembleJacobian(matrices_ptr, mat);
 
-  // contribution from nonlinearity
-  // add gradient to rhs, hessian to matrix!?
+  EnergyNorm<Matrix, Vector> energyNorm(mat);
+  LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
+                                           velocityTolerance_, &energyNorm,
+                                           verbosity_, false); // absolute error
 
   size_t fixedPointIteration;
   size_t multigridIterations = 0;
@@ -82,15 +76,43 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
   for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
        ++fixedPointIteration) {
 
-    // solve a velocity problem
+    // compute relative velocities
+    std::vector<Vector> v_rel;
+    relativeVelocities(velocityIterates, v_rel);
+
+    // contribution from nonlinearity
     for (size_t i=0; i<alpha.size(); i++) {
       globalFriction_[i]->updateAlpha(alpha[i]);
     }
 
+    std::vector<Matrix> matrices(velocityMatrices.size());
+    std::vector<Vector> rhs(velocityRHSs.size());
+    for (size_t i=0; i<globalFriction_.size(); i++) {
+      matrices[i] = velocityMatrices[i];
+      rhs[i] = velocityRHSs[i];
+
+      globalFriction_[i]->addHessian(v_rel[i], matrices[i]);
+      globalFriction_[i]->addGradient(v_rel[i], rhs[i]);
+
+      matrices_ptr[i] = &matrices[i];
+    }
+
+    // assemble full global contact problem
+    Matrix bilinearForm;
+    contactAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+    Vector totalRhs;
+    contactAssembler.assembleRightHandSide(rhs, totalRhs);
+
+    Vector totalVelocityIterate;
+    contactAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
+
+    // solve a velocity problem
     step_->setProblem(bilinearForm, totalVelocityIterate, totalRhs);
 
     velocityProblemSolver.preprocess();
     velocityProblemSolver.solve();
+    contactAssembler.postprocess(totalVelocityIterate, velocityIterates);
 
     multigridIterations += velocityProblemSolver.getResult().iterations;
 
@@ -103,19 +125,26 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
     }
 
     // compute relative velocities on contact boundaries
-    relativeVelocities(v_m);
+    relativeVelocities(v_m, v_rel);
 
     // solve a state problem
-    updaters.state_->solve(v_m);
-    ScalarVector newAlpha;
-   /* updaters.state_->extractAlpha(newAlpha);
+    updaters.state_->solve(v_rel);
+    std::vector<ScalarVector> newAlpha;
+    updaters.state_->extractAlpha(newAlpha);
 
-    if (lambda_ < 1e-12 or
-        errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
+    bool breakCriterion = true;
+    for (size_t i=0; i<alpha.size(); i++) {
+        if (errorNorm_.diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
+            breakCriterion = false;
+            break;
+        }
+    }
+
+    if (lambda_ < 1e-12 or breakCriterion) {
       fixedPointIteration++;
       break;
     }
-    alpha = newAlpha;*/
+    alpha = newAlpha;
   }
   if (fixedPointIteration == fixedPointMaxIterations_)
     DUNE_THROW(Dune::Exception, "FPI failed to converge");
@@ -137,7 +166,14 @@ std::ostream &operator<<(std::ostream &stream,
 }
 
 template <class Factory, class Updaters, class ErrorNorm>
-void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(std::vector<Vector>& v_m) const {
+void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(const std::vector<Vector>& v, std::vector<Vector>& v_rel) const {
+  // init result
+  v_rel.resize(v.size());
+  for (size_t i=0; i<v_rel.size(); i++) {
+      v_rel[i].resize(v[i].size());
+      v_rel[i] = 0;
+  }
+
   // needs assemblers to obtain basis
     /*
   std::vector<std::shared_ptr<MyAssembler>> assemblers(bodyCount);
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index ed1fe9b4957fddfb26291f99042dedc708fac353..e579645bbab839c82a714dd098560f283e3ed59f 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -32,7 +32,7 @@ class FixedPointIterator {
   using DeformedGrid = typename Factory::DeformedGrid;
 
 private:
-  void relativeVelocities(std::vector<Vector>& v_m) const;
+  void relativeVelocities(const std::vector<Vector>& v, std::vector<Vector>& v_rel) const;
 
 public:
   FixedPointIterator(Factory& factory, const Dune::ParameterTree& parset,
@@ -40,7 +40,7 @@ class FixedPointIterator {
                      const ErrorNorm& errorNorm);
 
   FixedPointIterationCounter run(Updaters updaters,
-                                 const std::vector<const Matrix*>& velocityMatrices,
+                                 const std::vector<Matrix>& velocityMatrices,
                                  const std::vector<Vector>& velocityRHSs,
                                  std::vector<Vector>& velocityIterates);
 
diff --git a/src/time-stepping/adaptivetimestepper.cc b/src/time-stepping/adaptivetimestepper.cc
index 2598ecaf14e377f1bcf6b8dab695e83a1403c894..0808750ff923408f6039721d84ec9193c3d2bbb1 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -22,7 +22,7 @@ AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
     Factory &factory, Dune::ParameterTree const &parset,
     std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, Updaters &current,
     double relativeTime, double relativeTau,
-    std::function<void(double, Vector &)> externalForces,
+    std::vector<std::function<void(double, Vector &)>>& externalForces,
     ErrorNorm const &errorNorm,
     std::function<bool(Updaters &, Updaters &)> mustRefine)
     : relativeTime_(relativeTime),
diff --git a/src/time-stepping/adaptivetimestepper.hh b/src/time-stepping/adaptivetimestepper.hh
index 5d559a8b097758f4f0ea1633eef1cf8b2a11545f..a0bf611397ef6a25663127515f9d50439cbe0fb1 100644
--- a/src/time-stepping/adaptivetimestepper.hh
+++ b/src/time-stepping/adaptivetimestepper.hh
@@ -33,7 +33,7 @@ class AdaptiveTimeStepper {
                       std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
                       Updaters &current, double relativeTime,
                       double relativeTau,
-                      std::function<void(double, Vector &)> externalForces,
+                      std::vector<std::function<void(double, Vector &)>>& externalForces,
                       ErrorNorm const &errorNorm,
                       std::function<bool(Updaters &, Updaters &)> mustRefine);
 
@@ -53,7 +53,7 @@ class AdaptiveTimeStepper {
   std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
   Updaters &current_;
   UpdatersWithCount R1_;
-  std::function<void(double, Vector &)> externalForces_;
+  std::vector<std::function<void(double, Vector &)>>& externalForces_;
   std::function<bool(Updaters &, Updaters &)> mustRefine_;
   ErrorNorm const &errorNorm_;
 
diff --git a/src/time-stepping/coupledtimestepper.cc b/src/time-stepping/coupledtimestepper.cc
index 410a48f8d65812124e3edbeceb0b5700a015ddae..7501f68a18d0eedfeef3c70313314dcdb109363a 100644
--- a/src/time-stepping/coupledtimestepper.cc
+++ b/src/time-stepping/coupledtimestepper.cc
@@ -9,7 +9,7 @@ CoupledTimeStepper<Factory, Updaters, ErrorNorm>::CoupledTimeStepper(
     double finalTime, Factory &factory, Dune::ParameterTree const &parset,
     std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, Updaters updaters,
     ErrorNorm const &errorNorm,
-    std::function<void(double, Vector &)> externalForces)
+    std::vector<std::function<void(double, Vector &)>>& externalForces)
     : finalTime_(finalTime),
       factory_(factory),
       parset_(parset),
diff --git a/src/time-stepping/coupledtimestepper.hh b/src/time-stepping/coupledtimestepper.hh
index 52543d1e8dd2b9204dd6fd978ba26aa2f833be0b..92e3cd4aa4c94a6ea39eb65fefbe5cf3b08e8553 100644
--- a/src/time-stepping/coupledtimestepper.hh
+++ b/src/time-stepping/coupledtimestepper.hh
@@ -19,7 +19,7 @@ class CoupledTimeStepper {
                      Dune::ParameterTree const &parset,
                      std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
                      Updaters updaters, ErrorNorm const &errorNorm,
-                     std::function<void(double, Vector &)> externalForces);
+                     std::vector<std::function<void(double, Vector &)>>& externalForces);
 
   FixedPointIterationCounter step(double relativeTime, double relativeTau);
 
@@ -29,7 +29,7 @@ class CoupledTimeStepper {
   Dune::ParameterTree const &parset_;
   std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
   Updaters updaters_;
-  std::function<void(double, Vector &)> externalForces_;
+  std::vector<std::function<void(double, Vector &)>>& externalForces_;
   ErrorNorm const &errorNorm_;
 };
 #endif