diff --git a/src/data-structures/contactnetwork.cc b/src/data-structures/contactnetwork.cc
index a25f7e49997bbe4a46a8bf3cccef9ce46c3dc259..84517f203a370cca87a3ecf6c58b2a5a2d39f29f 100644
--- a/src/data-structures/contactnetwork.cc
+++ b/src/data-structures/contactnetwork.cc
@@ -26,7 +26,7 @@ ContactNetwork<HostGridType, VectorType>::ContactNetwork(
     couplings_(nCouplings),
     frictionBoundaries_(nBodies),
     stateEnergyNorms_(nBodies),
-    nBodyAssembler_(nBodies, nCouplings) //, false, 0.0)
+    nBodyAssembler_(nBodies, nCouplings, false, 0.0)
 {}
 
 template <class HostGridType, class VectorType>
@@ -324,13 +324,17 @@ void ContactNetwork<HostGridType, VectorType>::totalNodes(
         std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> boundaryConditions;
         body->boundaryConditions(tag, boundaryConditions);
 
-        const int idxBackup = idx;
-        for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
-            const auto& nodes = boundaryConditions[bc]->boundaryNodes();
-            for (size_t j=0; j<nodes->size(); j++, idx++)
-                for (int k=0; k<dim; k++)
-                    totalNodes[idx][k] = (*nodes)[j][k];
-            idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
+        if (boundaryConditions.size()>0) {
+            const int idxBackup = idx;
+            for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
+                const auto& nodes = boundaryConditions[bc]->boundaryNodes();
+                for (size_t j=0; j<nodes->size(); j++, idx++)
+                    for (int k=0; k<dim; k++)
+                        totalNodes[idx][k] = (*nodes)[j][k];
+                idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
+            }
+        } else {
+            idx += body->nVertices();
         }
     }
 }
diff --git a/src/multi-body-problem-2D.cfg b/src/multi-body-problem-2D.cfg
index a4a778085c7e6d7b05a32c241c308874c61027de..a308c80b63b9cb2b48f0cc096ade7d9cdf464e57 100644
--- a/src/multi-body-problem-2D.cfg
+++ b/src/multi-body-problem-2D.cfg
@@ -1,6 +1,6 @@
 # -*- mode:conf -*-
 [boundary.friction]
-smallestDiameter = 0.25  # 2e-3 [m]
+smallestDiameter = 0.01  # 2e-3 [m]
 
 [timeSteps]
 refinementTolerance = 1e-5 # 1e-5
diff --git a/src/multi-body-problem-data/bc.hh b/src/multi-body-problem-data/bc.hh
index c7c919ae225c007b16fcd9c3242e88bccd0e86f3..b9b7ecdad223151dee490f24dae46e312c7c649c 100644
--- a/src/multi-body-problem-data/bc.hh
+++ b/src/multi-body-problem-data/bc.hh
@@ -9,12 +9,12 @@ class VelocityDirichletCondition
     // Assumed to vanish at time zero
       double const finalVelocity = -5e-5;
     
-    std::cout << "VelocityDirichletCondition::evaluate()" << std::endl;
+    /*std::cout << "VelocityDirichletCondition::evaluate()" << std::endl;
     
     if (relativeTime <= 0.1)
         std::cout << "- loading phase" << std::endl;
     else
-        std::cout << "- final velocity reached" << std::endl;
+        std::cout << "- final velocity reached" << std::endl;*/
     
     y = (relativeTime <= 0.1)
             ? finalVelocity * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0
diff --git a/src/multi-body-problem.cfg b/src/multi-body-problem.cfg
index 7188eed1b5851c795b2903d8000846c4a89c6450..e8a3e97bb00a8c6044364082c4800956d6d2f085 100644
--- a/src/multi-body-problem.cfg
+++ b/src/multi-body-problem.cfg
@@ -11,7 +11,7 @@ vtk.write       = true
 
 [problem]
 finalTime       = 100  # [s] #1000
-bodyCount       = 2
+bodyCount       = 4
 
 [body]
 bulkModulus     = 0.5e5 # [Pa]
diff --git a/src/solverfactorytest.cc b/src/solverfactorytest.cc
index 19b6896559823b0c748b103fa599b82c2d135b94..1d0f983b4fc2ad19af9cd16d0581ef44ea04ab57 100644
--- a/src/solverfactorytest.cc
+++ b/src/solverfactorytest.cc
@@ -57,7 +57,7 @@
 
 const int dim = MY_DIM;
 
-const int timeSteps = 1;
+const int timeSteps = 5;
 
 Dune::ParameterTree parset;
 size_t bodyCount;
@@ -86,6 +86,11 @@ void solveProblem(const ContactNetwork& contactNetwork,
     //using LocalSolver = LocalBisectionSolver;
     //using Linearization = Linearization<Functional, BitVector>;
 
+    /*print(ignore, "ignore:");
+    for (size_t i=0; i<x.size(); i++) {
+        std::cout << x[i] << std::endl;
+    }*/
+
     // set up reference solver
     Vector refX = x;
     using ContactFunctional = Dune::TNNMG::BoxConstrainedQuadraticFunctional<Matrix&, Vector&, Vector&, Vector&, double>;
@@ -438,8 +443,8 @@ void run(Updaters& updaters, ContactNetwork& contactNetwork,
   nBodyAssembler.concatenateVectors(u0, totalu0);
   nBodyAssembler.concatenateVectors(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();
 
@@ -465,7 +470,7 @@ void run(Updaters& updaters, ContactNetwork& contactNetwork,
     Vector totalVelocityIterate;
     nBodyAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
 
-    solveProblem(contactNetwork, bilinearForm, totalRhs, totalVelocityIterate, totalDirichletNodes, vLower, vUpper);
+    solveProblem(contactNetwork, bilinearForm, totalRhs, totalVelocityIterate, totalDirichletNodes, lower, upper);
 
     nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
 
@@ -496,8 +501,8 @@ void run(Updaters& updaters, ContactNetwork& contactNetwork,
         if (alpha[i].size()==0 || newAlpha[i].size()==0)
             continue;
 
-        print(alpha[i], "alpha i:");
-        print(newAlpha[i], "new alpha i:");
+        //print(alpha[i], "alpha i:");
+        //print(newAlpha[i], "new alpha i:");
         if (errorNorms[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance) {
             breakCriterion = false;
             std::cout << "fixedPoint error: " << errorNorms[i]->diff(alpha[i], newAlpha[i]) << std::endl;
@@ -674,10 +679,10 @@ int main(int argc, char *argv[]) { try {
     // Set up TNNMG solver
     // -------------------
 
-    BitVector totalDirichletNodes;
+    /*BitVector totalDirichletNodes;
     contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
 
-    //print(totalDirichletNodes, "totalDirichletNodes:");
+    print(totalDirichletNodes, "totalDirichletNodes:");*/
 
     //using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, field_type>;
     //using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
diff --git a/src/time-stepping/rate/newmark.cc b/src/time-stepping/rate/newmark.cc
index a3b3638529c8a75a09f633a5ea6a26b9ef6eb867..0ec590126ea5fd19755e6699324a259b5b5daeb8 100644
--- a/src/time-stepping/rate/newmark.cc
+++ b/src/time-stepping/rate/newmark.cc
@@ -95,20 +95,20 @@ void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
             double dirichletValue;
             bodyDirichletFunctions[bc]->evaluate(relativeTime, dirichletValue);
 
-            std::cout << "dirichletValue: " << dirichletValue << std::endl;
+            //std::cout << "dirichletValue: " << dirichletValue << std::endl;
 
                   
             for (size_t k=0; k<bcDirichletNodes.size() ; ++k) {
                 for (size_t j=0; j<dim; ++j) {
                     if (bcDirichletNodes[k][j]) {
-                              std::cout << k << " " << j << std::endl;
+                              //std::cout << k << " " << j << std::endl;
                         bodyIterate[k][j] = dirichletValue;
                     }
                 }
             }
         }
     }
-    print(iterate, "iterate:");
+    //sprint(iterate, "iterate:");
 }
 
 template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
diff --git a/src/time-stepping/state/ageinglawstateupdater.cc b/src/time-stepping/state/ageinglawstateupdater.cc
index 796f8878618a8a46174f8dac5dbdb5e7acedc621..1a7dcea3d34302f93f0ca1cb9fc2f8d9e2868f87 100644
--- a/src/time-stepping/state/ageinglawstateupdater.cc
+++ b/src/time-stepping/state/ageinglawstateupdater.cc
@@ -66,7 +66,7 @@ template <class ScalarVector, class Vector>
 void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
         ScalarVector& alpha) {
 
-    std::cout << "alpha size: " << alpha.size() << " nodes_.size() " << nodes_.size() << std::endl;
+    //std::cout << "alpha size: " << alpha.size() << " nodes_.size() " << nodes_.size() << std::endl;
     if (alpha.size() != nodes_.size()) {
         alpha.resize(nodes_.size());
     }