diff --git a/examples/benchmark_FrustrumDC.cpp b/examples/benchmark_FrustrumDC.cpp
index c9e74c0..8e0b031 100644
--- a/examples/benchmark_FrustrumDC.cpp
+++ b/examples/benchmark_FrustrumDC.cpp
@@ -46,10 +46,11 @@ int main (int argc, char** argv)
     index_t impl = 1; // 1= analytical, 2= generalized, 3= spectral
 
     int result        = 0;
-    index_t maxit     = 50;
+    index_t maxit     = 20;
     // Arc length method options
-    real_t dL         = 2e-2; // General arc length
-    real_t tol        = 1e-6;
+    real_t dL         = -1; // General arc length
+    real_t tolU        = 1e-3;
+    real_t tolF        = 1e-3;
 
     std::string wn("data.csv");
 
@@ -112,6 +113,7 @@ int main (int argc, char** argv)
     */
     if (testCase == 0)
     {
+        if (dL==-1) dL = 5e-2;
         displ.setValue(-dL,3);
         BCs.addCondition(boundary::north, condition_type::dirichlet, 0, 0, false, 0 ); // unknown 2 - z
         BCs.addCondition(boundary::north, condition_type::dirichlet, 0, 0, false, 1 ); // unknown 2 - z
@@ -136,6 +138,7 @@ int main (int argc, char** argv)
     }
     else if (testCase == 1)
     {
+        if (dL==-1) dL = 4e-2;
         displ.setValue(-dL,3);
 
         BCs.addCondition(boundary::north, condition_type::dirichlet, &displ, 0, false, 2 ); // unknown 1 - y
@@ -298,40 +301,44 @@ int main (int argc, char** argv)
       return assembler->rhs();
     };
 
-    gsSparseMatrix<> matrix;
-    gsVector<> vector;
+    assembler = new gsThinShellAssembler<3, real_t, true >(mp_def,dbasis,BCs,force,materialMatrix);
+    assembler->assemble();
+    gsSparseMatrix<> matrix = assembler->matrix();
+    gsVector<> vector = assembler->rhs();
 
     gsStaticSolver<real_t> staticSolver(matrix,vector,Jacobian,Residual);
     gsOptionList solverOptions = staticSolver.options();
     solverOptions.setInt("Verbose",true);
     solverOptions.setInt("MaxIterations",maxit);
-    solverOptions.setReal("Tolerance",tol);
+    solverOptions.setReal("ToleranceU",tolU);
+    solverOptions.setReal("ToleranceF",tolF);
     staticSolver.setOptions(solverOptions);
 
     real_t dL0 = dL;
     int reset = 0;
     gsMultiPatch<> mp_def0 = mp_def;
     real_t indicator;
+
+    real_t D = 0;
+
     for (index_t k=0; k<step; k++)
     {
       gsInfo<<"Load step "<< k<<"\n";
 
-      assembler = new gsThinShellAssembler<3, real_t, true >(mp_def,dbasis,BCs,force,materialMatrix);
-
       stopwatch.restart();
       stopwatch2.restart();
-      assembler->assemble();
+      displ.setValue(D - dL,3);
+      assembler->updateBCs(BCs);
+
       time += stopwatch.stop();
 
-      matrix = assembler->matrix();
-      vector = assembler->rhs();
       solVector = staticSolver.solveNonlinear();
       totaltime += stopwatch2.stop();
 
       if (!staticSolver.converged())
       {
         dL = dL/2;
-        displ.setValue(-dL,3);
+        displ.setValue(D -dL,3);
         reset = 1;
         mp_def = mp_def0;
         gsInfo<<"Iterations did not converge\n";
@@ -382,11 +389,12 @@ int main (int argc, char** argv)
       if (reset!=1)
       {
         dL = dL0;
-        displ.setValue(-dL,3);
+        displ.setValue(D - dL,3);
       }
       reset = 0;
 
       mp_def0 = mp_def;
+      D -= dL;
 
       gsInfo<<"--------------------------------------------------------------------------------------------------------------\n";
     }
@@ -396,6 +404,9 @@ int main (int argc, char** argv)
     gsInfo<<"Total ellapsed assembly time: \t\t"<<time<<" s\n";
     gsInfo<<"Total ellapsed solution time (incl. assembly): \t"<<totaltime<<" s\n";
 
+  delete assembler;
+  delete materialMatrix;
+
   return result;
 }
 
diff --git a/gsStaticSolver.h b/gsStaticSolver.h
index 00c95ab..5c2e005 100644
--- a/gsStaticSolver.h
+++ b/gsStaticSolver.h
@@ -91,7 +91,7 @@ class gsStaticSolver
     const std::function < gsSparseMatrix<T> ( gsVector<T> const & ) > m_nonlinear;
     const std::function < gsVector<T> ( gsVector<T> const & ) > m_residual;
 
-    gsVector<T> m_solVec;
+    gsVector<T> m_solVec, m_DeltaU, m_deltaU;
 
     mutable index_t m_verbose;
     bool m_NL;
diff --git a/gsStaticSolver.hpp b/gsStaticSolver.hpp
index c8c9ce3..d63c8b7 100644
--- a/gsStaticSolver.hpp
+++ b/gsStaticSolver.hpp
@@ -45,8 +45,8 @@ gsVector<T> gsStaticSolver<T>::solveNonlinear()
     if (residual==0) residual=1;
     T residual0 = residual;
     T residualOld = residual;
-    gsVector<T> DeltaU = gsVector<T>::Zero(m_solVec.rows());
-    gsVector<T> deltaU = gsVector<T>::Zero(m_solVec.rows());
+    m_DeltaU = gsVector<T>::Zero(m_solVec.rows());
+    m_deltaU = gsVector<T>::Zero(m_solVec.rows());
     gsSparseMatrix<T> jacMat;
 
 
@@ -66,17 +66,17 @@ gsVector<T> gsStaticSolver<T>::solveNonlinear()
 
     for (m_iterations = 0; m_iterations != m_maxIterations; ++m_iterations)
     {
-        jacMat = m_nonlinear(m_solVec+DeltaU);
+        jacMat = m_nonlinear(m_solVec+m_DeltaU);
         if (m_verbose==2)
         {
             gsInfo<<"Matrix: \n"<<jacMat.toDense()<<"\n";
             gsInfo<<"Vector: \n"<<resVec<<"\n";
         }
         m_solver.compute(jacMat);
-        deltaU = m_solver.solve(resVec); // this is the UPDATE
-        DeltaU += m_relax * deltaU;
+        m_deltaU = m_solver.solve(resVec); // this is the UPDATE
+        m_DeltaU += m_relax * m_deltaU;
 
-        resVec = m_residual(m_solVec+DeltaU);
+        resVec = m_residual(m_solVec+m_DeltaU);
         residual = resVec.norm();
 
         if (m_verbose>0)
@@ -85,9 +85,9 @@ gsVector<T> gsStaticSolver<T>::solveNonlinear()
             gsInfo<<std::setw(4)<<std::left<<m_iterations;
             gsInfo<<std::setw(17)<<std::left<<residual;
             gsInfo<<std::setw(17)<<std::left<<residual/residual0;
-            gsInfo<<std::setw(17)<<std::left<<m_relax * deltaU.norm();
-            gsInfo<<std::setw(17)<<std::left<<m_relax * deltaU.norm()/DeltaU.norm();
-            gsInfo<<std::setw(17)<<std::left<<m_relax * deltaU.norm()/m_solVec.norm();
+            gsInfo<<std::setw(17)<<std::left<<m_relax * m_deltaU.norm();
+            gsInfo<<std::setw(17)<<std::left<<m_relax * m_deltaU.norm()/m_DeltaU.norm();
+            gsInfo<<std::setw(17)<<std::left<<m_relax * m_deltaU.norm()/m_solVec.norm();
             gsInfo<<std::setw(17)<<std::left<<math::log10(residualOld/residual0);
             gsInfo<<std::setw(17)<<std::left<<math::log10(residual/residual0);
             gsInfo<<"\n";
@@ -95,10 +95,10 @@ gsVector<T> gsStaticSolver<T>::solveNonlinear()
 
         residualOld = residual;
 
-        if (m_relax * deltaU.norm()/m_solVec.norm()  < m_toleranceU && residual/residual0 < m_toleranceF)
+        if (m_relax * m_deltaU.norm()/m_solVec.norm()  < m_toleranceU && residual/residual0 < m_toleranceF)
         {
             m_converged = true;
-            m_solVec+=DeltaU;
+            m_solVec+=m_DeltaU;
             break;
         }
         else if (m_iterations+1 == m_maxIterations)