diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp
index e1139fe..d36eeff 100755
--- a/drmsimulationmodel.cpp
+++ b/drmsimulationmodel.cpp
@@ -222,6 +222,7 @@ void DRMSimulationModel::reset()
     mCurrentSimulationStep = 0;
     history.clear();
     constrainedVertices.clear();
+    rigidSupports.clear();
     pMesh.reset();
     plotYValues.clear();
     plotHandle.reset();
@@ -692,8 +693,8 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e,
         derivativeF1 = derivativeT1_k - (n * (derivativeT1_k * n));
 
         const double f1Norm = f1.Norm();
-        const VectorType derivativeF1Normalized = -f1 * (f1 * derivativeF1)
-                                                      / (f1Norm * f1Norm * f1Norm)
+        const double derivativeF1Norm = f1 * derivativeF1 / f1Norm;
+        const VectorType derivativeF1Normalized = -f1 * derivativeF1Norm / (f1Norm * f1Norm)
                                                   + derivativeF1 / f1Norm;
 
         derivativeF2 = derivativeF1Normalized * cosRotationAngle
@@ -713,8 +714,8 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e,
             const VectorType &derivativeOfNormal = node.derivativeOfNormal[dofi];
             derivativeF1 = -(n * (t1_k * derivativeOfNormal) + derivativeOfNormal * (t1_k * n));
             const double f1Norm = f1.Norm();
-            const VectorType derivativeF1Normalized = -f1 * (f1 * derivativeF1)
-                                                          / (f1Norm * f1Norm * f1Norm)
+            const double derivativeF1Norm = f1 * derivativeF1 / f1Norm;
+            const VectorType derivativeF1Normalized = -f1 * derivativeF1Norm / (f1Norm * f1Norm)
                                                       + derivativeF1 / f1Norm;
 
             derivativeF2 = derivativeF1Normalized * cosRotationAngle +
@@ -805,9 +806,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
         pMesh->EN(), std::vector<std::pair<int, Vector6d>>(4, {-1, Vector6d()}));
     //  omp_lock_t writelock;
     //  omp_init_lock(&writelock);
-#ifdef ENABLE_OPENMP
-#pragma omp parallel for schedule(static) num_threads(5)
-#endif
+    //#ifdef ENABLE_OPENMP
+    //#pragma omp parallel for schedule(static) num_threads(8)
+    //#endif
     for (int ei = 0; ei < pMesh->EN(); ei++) {
         const EdgeType &e = pMesh->edge[ei];
         const SimulationMesh::VertexType &ev_j = *e.cV(0);
@@ -1007,7 +1008,6 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
         sumOfResidualForces = sumOfResidualForces + nodeResidualForce;
         const double residualForceNorm = nodeResidualForce.norm();
         //        const double residualForceNorm = nodeResidualForce.getTranslation().norm();
-        const bool shouldTrigger = std::isnan(residualForceNorm);
         totalResidualForcesNorm += residualForceNorm;
 #ifdef POLYSCOPE_DEFINED
         if (std::isnan(force.residual.norm())) {
@@ -1192,10 +1192,10 @@ void DRMSimulationModel::updateNodalMasses()
     if (shouldTemporarilyDampForces && mCurrentSimulationStep < untilStep) {
         gamma *= 1e6 * (1 - static_cast<double>(mCurrentSimulationStep) / untilStep);
     }
-    if (mCurrentSimulationStep == static_cast<size_t>(1.2 * untilStep)
-        && shouldTemporarilyDampForces) {
-        Dt = mSettings.Dtini;
-    }
+    //    if (mCurrentSimulationStep == static_cast<size_t>(1.4 * untilStep)
+    //        && shouldTemporarilyDampForces) {
+    //        Dt = mSettings.Dtini;
+    //    }
     for (VertexType &v : pMesh->vert) {
         const size_t vi = pMesh->getIndex(v);
         double translationalSumSk = 0;
@@ -1777,7 +1777,7 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
         accelerationX[pMesh->getIndex(v)] = pMesh->nodes[v].acceleration[0];
     }
     meshPolyscopeHandle->addNodeScalarQuantity("Kinetic Energy", kineticEnergies)->setEnabled(false);
-    meshPolyscopeHandle->addNodeVectorQuantity("Node normals", nodeNormals)->setEnabled(false);
+    meshPolyscopeHandle->addNodeVectorQuantity("Node normals", nodeNormals)->setEnabled(true);
     meshPolyscopeHandle->addNodeVectorQuantity("Internal force", internalForces)->setEnabled(false);
     meshPolyscopeHandle->addNodeVectorQuantity("External force", externalForces)->setEnabled(false);
     meshPolyscopeHandle->addNodeVectorQuantity("External Moments", externalMoments)->setEnabled(true);
@@ -2089,7 +2089,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
     vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh);
     computeRigidSupports();
     isVertexConstrained.resize(pMesh->VN(), false);
-    for (auto fixedVertex : pJob->constrainedVertices) {
+    for (auto fixedVertex : constrainedVertices) {
         isVertexConstrained[fixedVertex.first] = true;
     }
 
@@ -2285,19 +2285,19 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
             //            SimulationResultsReporter::createPlot("Number of Steps",
             //                                                  "Residual Forces mov aver deriv",
             //                                                  movingAveragesDerivatives_norm);
-            //            SimulationResultsReporter::createPlot("Number of Steps",
-            //                                                  "Residual Forces mov aver",
-            //                                                  history.residualForcesMovingAverage);
+            SimulationResultsReporter::createPlot("Number of Steps",
+                                                  "Residual Forces mov aver",
+                                                  history.residualForcesMovingAverage);
             //            SimulationResultsReporter::createPlot("Number of Steps",
             //                                                  "Log of Residual Forces",
             //                                                  history.logResidualForces,
             //                                                  {},
             //                                                  history.redMarks);
-            SimulationResultsReporter::createPlot("Number of Steps",
-                                                  "Log of Kinetic energy",
-                                                  history.kineticEnergy,
-                                                  {},
-                                                  history.redMarks);
+            //            SimulationResultsReporter::createPlot("Number of Steps",
+            //                                                  "Log of Kinetic energy",
+            //                                                  history.kineticEnergy,
+            //                                                  {},
+            //                                                  history.redMarks);
             //            SimulationResultsReporter reporter;
             //            reporter.reportHistory(history, "IntermediateResults", pJob->pMesh->getLabel());
             //            SimulationResultsReporter::createPlot("Number of Iterations",
@@ -2349,21 +2349,14 @@ currentSimulationStep > maxDRMIterations*/
             = mSettings.useAverage
               && (pMesh->totalResidualForcesNorm / pMesh->VN()) / totalExternalForcesNorm
                      < mSettings.averageResidualForcesCriterionThreshold;
-        if (fullfillsResidualForcesNormTerminationCriterion
-            || fullfillsResidualForcesNormTerminationCriterion) {
-            if (mSettings.beVerbose /*&& !mSettings.isDebugMode*/) {
-                std::cout << "Simulation converged." << std::endl;
-                printCurrentState();
-            }
-            break;
-        }
-        // Residual forces norm convergence
-        if ((mSettings.useKineticDamping)
-            && (pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy
 
-                //          && mCurrentSimulationStep > movingAverageSampleSize
-                && (pJob->nodalForcedDisplacements.empty()
-                    || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps))
+        // Residual forces norm convergence
+        if (((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy
+              || fullfillsAverageResidualForcesNormTerminationCriterion
+              || fullfillsResidualForcesNormTerminationCriterion)
+             //          && mCurrentSimulationStep > movingAverageSampleSize
+             && (pJob->nodalForcedDisplacements.empty()
+                 || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps))
             /* || (pMesh->totalResidualForcesNorm / mSettings.totalResidualForcesNormThreshold <= 1
               && mCurrentSimulationStep > 1)*/
             /*||
@@ -2388,7 +2381,7 @@ mesh->currentTotalPotentialEnergykN*/
                   //                || fullfillsMovingAverageDerivativeNormTerminationCriterion
                   || fullfillsAverageResidualForcesNormTerminationCriterion
                   || fullfillsResidualForcesNormTerminationCriterion;
-            if (shouldTerminate) {
+            if (shouldTerminate && mCurrentSimulationStep > 100) {
                 if (mSettings.beVerbose /*&& !mSettings.isDebugMode*/) {
                     std::cout << "Simulation converged." << std::endl;
                     printCurrentState();
diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp
index 680cc97..7e1d6d8 100755
--- a/drmsimulationmodel.hpp
+++ b/drmsimulationmodel.hpp
@@ -161,9 +161,9 @@ private:
 
   void updateElementalFrames();
 
-  VectorType computeDerivativeOfR(const EdgeType &e, const DifferentiateWithRespectTo &dui) const;
+  VectorType computeDerivativeOfR(const EdgeType &e,
+                                  const DifferentiateWithRespectTo &dui) const;
 
-  //  bool isRigidSupport(const VertexType &v) const;
 
   static double computeDerivativeOfNorm(const VectorType &x,
                                         const VectorType &derivativeOfX);