From 513529f0a4e333aa84a2e613061f23bd323082c8 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Fri, 16 Jul 2021 09:48:14 +0300 Subject: [PATCH] Refactoring --- drmsimulationmodel.cpp | 109 ++++++++++++++++++++++++++--------------- drmsimulationmodel.hpp | 3 +- simulationmesh.hpp | 1 + 3 files changed, 73 insertions(+), 40 deletions(-) diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index b2d6645..e1139fe 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -1235,7 +1235,24 @@ void DRMSimulationModel::updateNodalMasses() pMesh->nodes[v].mass_6d[DoF::Nx] = pMesh->nodes[v].mass.rotationalJ; pMesh->nodes[v].mass_6d[DoF::Ny] = pMesh->nodes[v].mass.rotationalI3; pMesh->nodes[v].mass_6d[DoF::Nr] = pMesh->nodes[v].mass.rotationalI2; - + if (mSettings.useViscousDamping) { + //fill 6d damping vector + const double translationalDampingFactor + = 2 * std::sqrt(translationalSumSk * pMesh->nodes[v].mass.translational); + pMesh->nodes[v].damping_6d[DoF::Ux] = translationalDampingFactor; + pMesh->nodes[v].damping_6d[DoF::Uy] = translationalDampingFactor; + pMesh->nodes[v].damping_6d[DoF::Uz] = translationalDampingFactor; + pMesh->nodes[v].damping_6d[DoF::Nx] = 2 + * std::sqrt(rotationalSumSk_J + * pMesh->nodes[v].mass_6d[DoF::Nx]); + pMesh->nodes[v].damping_6d[DoF::Ny] = 2 + * std::sqrt(rotationalSumSk_I3 + * pMesh->nodes[v].mass_6d[DoF::Ny]); + pMesh->nodes[v].damping_6d[DoF::Nr] = 2 + * std::sqrt(rotationalSumSk_I2 + * pMesh->nodes[v].mass_6d[DoF::Nr]); + pMesh->nodes[v].damping_6d = pMesh->nodes[v].damping_6d * 1e-2; + } assert(std::pow(mSettings.Dtini, 2.0) * translationalSumSk / pMesh->nodes[v].mass.translational < 2); @@ -1277,21 +1294,23 @@ void DRMSimulationModel::updateNodalVelocities() for (VertexType &v : pMesh->vert) { const VertexIndex vi = pMesh->getIndex(v); Node &node = pMesh->nodes[v]; -#ifdef POLYSCOPE_DEFINED - if (std::isnan(node.velocity.norm())) { - std::cout << "Velocity " << vi << ":" << node.velocity.toString() << std::endl; - } -#endif if (mSettings.useViscousDamping) { const Vector6d massOverDt = node.mass_6d / Dt; - const Vector6d denominator = massOverDt + Vector6d(viscuousDampingConstant / 2); - const Vector6d firstTermNominator = massOverDt - Vector6d(viscuousDampingConstant / 2); + const Vector6d visciousDampingFactor(viscuousDampingConstant / 2); + // const Vector6d &visciousDampingFactor = node.damping_6d; + const Vector6d denominator = massOverDt + visciousDampingFactor / 2; + const Vector6d firstTermNominator = massOverDt - visciousDampingFactor / 2; const Vector6d firstTerm = node.velocity * firstTermNominator / denominator; const Vector6d secondTerm = node.force.residual / denominator; node.velocity = firstTerm + secondTerm; } else { node.velocity = node.velocity + node.acceleration * Dt; } +#ifdef POLYSCOPE_DEFINED + if (std::isnan(node.velocity.norm())) { + std::cout << "Velocity " << vi << ":" << node.velocity.toString() << std::endl; + } +#endif } updateKineticEnergy(); } @@ -2266,19 +2285,19 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrpMesh->getLabel()); // SimulationResultsReporter::createPlot("Number of Iterations", @@ -2330,13 +2349,21 @@ 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 (((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy - || fullfillsAverageResidualForcesNormTerminationCriterion - || fullfillsResidualForcesNormTerminationCriterion) - // && mCurrentSimulationStep > movingAverageSampleSize - && (pJob->nodalForcedDisplacements.empty() - || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)) + if ((mSettings.useKineticDamping) + && (pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy + + // && mCurrentSimulationStep > movingAverageSampleSize + && (pJob->nodalForcedDisplacements.empty() + || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)) /* || (pMesh->totalResidualForcesNorm / mSettings.totalResidualForcesNormThreshold <= 1 && mCurrentSimulationStep > 1)*/ /*|| @@ -2405,25 +2432,27 @@ mesh->currentTotalPotentialEnergykN*/ // } } - const bool shouldCapDisplacements = mSettings.displacementCap.has_value(); - for (VertexType &v : pMesh->vert) { - Node &node = pMesh->nodes[v]; - Vector6d stepDisplacement = node.velocity * 0.5 * Dt; - if (shouldCapDisplacements - && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) { - stepDisplacement = stepDisplacement - * (*mSettings.displacementCap - / stepDisplacement.getTranslation().norm()); + if (!mSettings.useViscousDamping) { + const bool shouldCapDisplacements = mSettings.displacementCap.has_value(); + for (VertexType &v : pMesh->vert) { + Node &node = pMesh->nodes[v]; + Vector6d stepDisplacement = node.velocity * 0.5 * Dt; + if (shouldCapDisplacements + && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) { + stepDisplacement = stepDisplacement + * (*mSettings.displacementCap + / stepDisplacement.getTranslation().norm()); + } + node.displacements = node.displacements - stepDisplacement; } - node.displacements = node.displacements - stepDisplacement; - } - applyDisplacements(constrainedVertices); - if (!pJob->nodalForcedDisplacements.empty()) { - applyForcedDisplacements( + applyDisplacements(constrainedVertices); + if (!pJob->nodalForcedDisplacements.empty()) { + applyForcedDisplacements( - pJob->nodalForcedDisplacements); + pJob->nodalForcedDisplacements); + } + updateElementalLengths(); } - updateElementalLengths(); // const double triggerPercentage = 0.01; // const double xi_min = 0.55; @@ -2436,7 +2465,9 @@ mesh->currentTotalPotentialEnergykN*/ // + xi_init - triggerPercentage * xi_min) // / (1 - triggerPercentage); // } - resetVelocities(); + if (mSettings.useKineticDamping) { + resetVelocities(); + } Dt *= mSettings.xi; // if (mSettings.isDebugMode) { // std::cout << Dt << std::endl; diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index e2afd98..680cc97 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -48,6 +48,7 @@ public: double averageResidualForcesCriterionThreshold{1e-5}; Settings() {} bool useViscousDamping{false}; + bool useKineticDamping{true}; }; private: @@ -62,7 +63,7 @@ private: std::vector plotYValues; size_t numOfDampings{0}; int externalLoadStep{1}; - const double viscuousDampingConstant{0.01}; + const double viscuousDampingConstant{100}; std::vector isVertexConstrained; std::vector isRigidSupport; double minTotalResidualForcesNorm{std::numeric_limits::max()}; diff --git a/simulationmesh.hpp b/simulationmesh.hpp index 234a535..a7b85c9 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -149,6 +149,7 @@ struct Node { Mass mass; Vector6d mass_6d; + Vector6d damping_6d; VertexIndex vi; CoordType initialLocation; CoordType initialNormal;