diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp
index 8b7ffaa..5960aa0 100755
--- a/drmsimulationmodel.cpp
+++ b/drmsimulationmodel.cpp
@@ -217,14 +217,22 @@ void DRMSimulationModel::runUnitTests()
 
     //  polyscope::show();
 }
+
 void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob)
 {
-    //#ifdef USE_ENSMALLEN
-    //    this->pJob = pJob;
-    //#endif
-    pMesh.reset();
-    pMesh = std::make_unique<SimulationMesh>(*pJob->pMesh);
-    vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh);
+    mCurrentSimulationStep = 0;
+    history.clear();
+    history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
+    plotYValues.clear();
+    plotHandle.reset();
+    checkedForMaximumMoment = false;
+    externalMomentsNorm = 0;
+    Dt = mSettings.Dtini;
+    numOfDampings = 0;
+    shouldTemporarilyDampForces = false;
+    externalLoadStep = 1;
+    isVertexConstrained.clear();
+    minTotalResidualForcesNorm = std::numeric_limits<double>::max();
 
     constrainedVertices.clear();
     constrainedVertices = pJob->constrainedVertices;
@@ -239,26 +247,6 @@ void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob)
     for (auto fixedVertex : constrainedVertices) {
         isVertexConstrained[fixedVertex.first] = true;
     }
-}
-
-void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings)
-{
-    mSettings = settings;
-    mCurrentSimulationStep = 0;
-    history.clear();
-    history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
-    plotYValues.clear();
-    plotHandle.reset();
-    checkedForMaximumMoment = false;
-    externalMomentsNorm = 0;
-    Dt = settings.Dtini;
-    numOfDampings = 0;
-    shouldTemporarilyDampForces = false;
-    externalLoadStep = 1;
-    isVertexConstrained.clear();
-    minTotalResidualForcesNorm = std::numeric_limits<double>::max();
-
-    reset(pJob);
 
 #ifdef POLYSCOPE_DEFINED
     if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) {
@@ -269,7 +257,7 @@ void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob, const
         polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel + "_" + pMesh->getLabel(),
                                         pMesh->getEigenVertices(),
                                         pMesh->getEigenEdges())
-            ->setRadius(0.002);
+            ->setRadius(pMesh->elements[0].dimensions.getDrawingRadius());
         //    registerWorldAxes();
     }
 #endif
@@ -1426,6 +1414,7 @@ void DRMSimulationModel::updateResidualForces()
 
 void DRMSimulationModel::computeRigidSupports()
 {
+    assert(pMesh != nullptr);
     isRigidSupport.clear();
     isRigidSupport.resize(pMesh->VN(), false);
     for (const VertexType &v : pMesh->vert) {
@@ -2611,6 +2600,8 @@ void DRMSimulationModel::applyKineticDamping(const std::shared_ptr<SimulationJob
 
 SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob)
 {
+    assert(pMesh != nullptr && pMesh->VN() == pJob->pMesh->VN() && pMesh->EN() == pJob->pMesh->EN());
+    reset(pJob);
     auto beginTime = std::chrono::high_resolution_clock::now();
 
     updateNodalMasses();
@@ -2975,11 +2966,11 @@ currentSimulationStep > maxDRMIterations*/
 
             if (mSettings.useTranslationalKineticEnergyForKineticDamping) {
                 applyKineticDamping(pJob);
-                Dt *= mSettings.xi;
                 if (mSettings.shouldCreatePlots) {
                     history.redMarks.push_back(mCurrentSimulationStep);
                 }
             }
+            Dt *= mSettings.xi;
             //            if (mSettings.isDebugMode) {
             //                std::cout << Dt << std::endl;
             //            }
@@ -3011,9 +3002,11 @@ currentSimulationStep > maxDRMIterations*/
 
 void DRMSimulationModel::setStructure(const std::shared_ptr<SimulationMesh> &pMesh)
 {
-    std::cout << "This function is currently not implemented" << std::endl;
-    assert(false);
-    std::terminate();
+    this->pMesh.reset();
+    this->pMesh = std::make_unique<SimulationMesh>(*pMesh);
+    vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh);
+    //    assert(false);
+    //    std::terminate();
 }
 
 SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
@@ -3021,7 +3014,8 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
                                                         const SimulationResults &solutionGuess)
 {
     std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
-    reset(pJob, settings);
+    mSettings = settings;
+    setStructure(pJob->pMesh);
     assert(pJob->pMesh != nullptr);
 
     if (!solutionGuess.displacements.empty()) {
diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp
index e420b76..0267a6f 100755
--- a/drmsimulationmodel.hpp
+++ b/drmsimulationmodel.hpp
@@ -48,7 +48,7 @@ public:
         //        int desiredGradualExternalLoadsSteps{1};
         double gamma{0.8};
         double totalResidualForcesNormThreshold{1e-20};
-        double totalExternalForcesNormPercentageTermination{1e-5};
+        double totalExternalForcesNormPercentageTermination{1e-3};
         std::optional<int> maxDRMIterations;
         std::optional<int> debugModeStep;
         std::optional<double> totalTranslationalKineticEnergyThreshold;
@@ -108,7 +108,6 @@ private:
     //#ifdef USE_ENSMALLEN
     //    std::shared_ptr<SimulationJob> pJob;
     //#endif
-    void reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings);
     void updateNodalInternalForces(
         const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
     void updateNodalExternalForces(