From a8b361f8597423138bfee5ef13f31179ab866f18 Mon Sep 17 00:00:00 2001
From: iasonmanolas <iasonmanolasece@gmail.com>
Date: Thu, 1 Jul 2021 02:45:03 +0300
Subject: [PATCH] Reusing initial version. reset residual force convergence to
 percentage

---
 drmsimulationmodel.cpp | 28 ++++++++++++++--------------
 drmsimulationmodel.hpp |  1 +
 2 files changed, 15 insertions(+), 14 deletions(-)

diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp
index 0f868bc..7fd187a 100755
--- a/drmsimulationmodel.cpp
+++ b/drmsimulationmodel.cpp
@@ -977,8 +977,8 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
         Node::Forces &force = pMesh->nodes[vi].force;
         const Vector6d &nodeResidualForce = force.residual;
         sumOfResidualForces = sumOfResidualForces + nodeResidualForce;
-        //        const double residualForceNorm = nodeResidualForce.norm();
-        const double residualForceNorm = nodeResidualForce.getTranslation().norm();
+        const double residualForceNorm = nodeResidualForce.norm();
+        //        const double residualForceNorm = nodeResidualForce.getTranslation().norm();
         const bool shouldTrigger = std::isnan(residualForceNorm);
         totalResidualForcesNorm += residualForceNorm;
     }
@@ -1461,7 +1461,7 @@ void DRMSimulationModel::updateKineticEnergy()
                  + node.mass.rotationalI3 * pow(node.velocity[4], 2)
                  + node.mass.rotationalI2 * pow(node.velocity[5], 2));
 
-        node.kineticEnergy += nodeTranslationalKineticEnergy /*+ nodeRotationalKineticEnergy*/;
+        node.kineticEnergy += nodeTranslationalKineticEnergy + nodeRotationalKineticEnergy;
         assert(node.kineticEnergy < 1e15);
 
         pMesh->currentTotalKineticEnergy += node.kineticEnergy;
@@ -1962,18 +1962,18 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
     auto t1 = std::chrono::high_resolution_clock::now();
     mSettings = settings;
     reset();
-    //    double totalExternalForcesNorm = 0;
-    //    for (const std::pair<VertexIndex, Vector6d> &externalForce : pJob->nodalExternalForces) {
-    //        totalExternalForcesNorm += externalForce.second.norm();
-    //    }
+    double totalExternalForcesNorm = 0;
+    for (const std::pair<VertexIndex, Vector6d> &externalForce : pJob->nodalExternalForces) {
+        totalExternalForcesNorm += externalForce.second.norm();
+    }
 
-    //    if (!pJob->nodalExternalForces.empty()) {
-    //        //        totalResidualForcesNormThreshold = settings.totalExternalForcesNormPercentageTermination
-    //        //                                           * totalExternalForcesNorm;
-    //    } else {
-    //        totalResidualForcesNormThreshold = 1;
-    //        std::cout << "No forces setted default residual forces norm threshold" << std::endl;
-    //    }
+    if (!pJob->nodalExternalForces.empty()) {
+        mSettings.totalResidualForcesNormThreshold
+            = settings.totalExternalForcesNormPercentageTermination * totalExternalForcesNorm;
+    } else {
+        mSettings.totalResidualForcesNormThreshold = 1e-3;
+        std::cout << "No forces setted default residual forces norm threshold" << std::endl;
+    }
 
     if (mSettings.beVerbose) {
         std::cout << "totalResidualForcesNormThreshold:"
diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp
index 4c62eb3..af43ad7 100755
--- a/drmsimulationmodel.hpp
+++ b/drmsimulationmodel.hpp
@@ -42,6 +42,7 @@ public:
         double gamma{0.8};
         std::optional<double> displacementCap;
         double totalResidualForcesNormThreshold{1e-3};
+        double totalExternalForcesNormPercentageTermination{1e-3};
         Settings() {}
     };