From 2509b0a7952151b5d6a4c4da50c7472dcb54f291 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Wed, 7 Jul 2021 22:43:44 +0300 Subject: [PATCH] Added code for exporting pe.gradual application of external forces. average residual forces threshold --- drmsimulationmodel.cpp | 250 ++++++++++++++++++------------ drmsimulationmodel.hpp | 4 + reducedmodeloptimizer_structs.hpp | 24 +++ simulation_structs.hpp | 4 +- simulationmesh.hpp | 1 + utilities.hpp | 10 ++ 6 files changed, 193 insertions(+), 100 deletions(-) diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 7fd187a..eb7d742 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -231,6 +231,7 @@ void DRMSimulationModel::reset() Dt = mSettings.Dtini; numOfDampings = 0; shouldTemporarilyDampForces = false; + externalLoadStep = 1; } VectorType DRMSimulationModel::computeDisplacementDifferenceDerivative( @@ -958,7 +959,6 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( force.residual = force.external; force.internal = 0; } - double totalResidualForcesNorm = 0; for (size_t ei = 0; ei < pMesh->EN(); ei++) { for (int i = 0; i < 4; i++) { std::pair internalForcePair @@ -972,6 +972,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( force.residual = force.residual + (internalForcePair.second * -1); } } + double totalResidualForcesNorm = 0; Vector6d sumOfResidualForces(0); for (size_t vi = 0; vi < pMesh->VN(); vi++) { Node::Forces &force = pMesh->nodes[vi].force; @@ -984,7 +985,8 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( } pMesh->previousTotalResidualForcesNorm = pMesh->totalResidualForcesNorm; pMesh->totalResidualForcesNorm = totalResidualForcesNorm; - // mesh->totalResidualForcesNorm = sumOfResidualForces.norm(); + pMesh->averageResidualForcesNorm = totalResidualForcesNorm / pMesh->VN(); + // pMesh->averageResidualForcesNorm = sumOfResidualForces.norm() / pMesh->VN(); // plotYValues.push_back(totalResidualForcesNorm); // auto xPlot = matplot::linspace(0, plotYValues.size(), plotYValues.size()); @@ -1147,7 +1149,7 @@ DRMSimulationModel::DRMSimulationModel() {} void DRMSimulationModel::updateNodalMasses() { double gamma = mSettings.gamma; - const size_t untilStep = 500; + const size_t untilStep = 4000; if (shouldTemporarilyDampForces && mCurrentSimulationStep < untilStep) { gamma *= 1e6 * (1 - static_cast(mCurrentSimulationStep) / untilStep); } @@ -1233,10 +1235,9 @@ void DRMSimulationModel::updateNodalVelocities() const VertexIndex vi = pMesh->getIndex(v); Node &node = pMesh->nodes[v]; node.velocity = node.velocity + node.acceleration * Dt; - // if (vi == 10) { - // std::cout << "Velocity:" << node.velocity[0] << " " << node.velocity[1] << " " - // << node.velocity[2] << std::endl; - // } + if (std::isnan(node.velocity.norm())) { + std::cout << "Velocity:" << node.velocity.toString() << std::endl; + } } updateKineticEnergy(); } @@ -1461,7 +1462,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; @@ -1959,26 +1960,10 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrpMesh.operator bool()); - auto t1 = std::chrono::high_resolution_clock::now(); + auto beginTime = std::chrono::high_resolution_clock::now(); mSettings = settings; reset(); - double totalExternalForcesNorm = 0; - for (const std::pair &externalForce : pJob->nodalExternalForces) { - totalExternalForcesNorm += externalForce.second.norm(); - } - 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:" - << mSettings.totalResidualForcesNormThreshold << std::endl; - } history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel(); // if (!pJob->nodalExternalForces.empty()) { @@ -2044,12 +2029,38 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrnodalExternalForces, constrainedVertices); - const size_t movingAverageSampleSize = 200; - std::queue residualForcesMovingAverageHistorySample; + std::unordered_map nodalExternalForces = pJob->nodalExternalForces; + double totalExternalForcesNorm = 0; + // Vector6d sumOfExternalForces(0); + for (auto &nodalForce : nodalExternalForces) { + const double percentageOfExternalLoads = double(externalLoadStep) + / mSettings.desiredGradualExternalLoadsSteps; + nodalForce.second = nodalForce.second * percentageOfExternalLoads; + totalExternalForcesNorm += nodalForce.second.norm(); + // sumOfExternalForces = sumOfExternalForces + nodalForce.second; + } + updateNodalExternalForces(nodalExternalForces, constrainedVertices); + double averageExternalForcesNorm = + // sumOfExternalForces.norm() / pMesh->VN(); + totalExternalForcesNorm / pMesh->VN(); + if (!nodalExternalForces.empty()) { + mSettings.totalResidualForcesNormThreshold + = std::min(settings.totalExternalForcesNormPercentageTermination + * totalExternalForcesNorm, + 1e-3); + } else { + mSettings.totalResidualForcesNormThreshold = 1e-3; + std::cout << "No forces setted default residual forces norm threshold" << std::endl; + } + if (mSettings.beVerbose) { + std::cout << "totalResidualForcesNormThreshold:" + << mSettings.totalResidualForcesNormThreshold << std::endl; + } + // const size_t movingAverageSampleSize = 200; + // std::queue residualForcesMovingAverageHistorySample; std::vector percentageOfConvergence; - double residualForcesMovingAverageDerivativeMax = 0; + // double residualForcesMovingAverageDerivativeMax = 0; while (settings.maxDRMIterations == 0 || mCurrentSimulationStep < settings.maxDRMIterations) { if ((mSettings.isDebugMode && mCurrentSimulationStep == 50000)) { // std::filesystem::create_directory("./PatternOptimizationNonConv"); @@ -2087,38 +2098,46 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrVN(); vi++) { - sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm(); + if (std::isnan(pMesh->currentTotalKineticEnergy)) { + std::cout << pMesh->currentTotalKineticEnergy << std::endl; + if (mSettings.beVerbose) { + std::cerr << "Infinite kinetic energy at step " << mCurrentSimulationStep + << ". Exiting.." << std::endl; + } + break; } - sumOfDisplacementsNorm /= pMesh->bbox.Diag(); - pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm; - // pMesh->residualForcesMovingAverage = (pMesh->totalResidualForcesNorm - // + mCurrentSimulationStep - // * pMesh->residualForcesMovingAverage) - // / (1 + mCurrentSimulationStep); - pMesh->residualForcesMovingAverage += pMesh->totalResidualForcesNorm - / movingAverageSampleSize; - residualForcesMovingAverageHistorySample.push(pMesh->residualForcesMovingAverage); - if (movingAverageSampleSize < residualForcesMovingAverageHistorySample.size()) { - const double firstElementValue = residualForcesMovingAverageHistorySample.front(); - pMesh->residualForcesMovingAverage -= firstElementValue / movingAverageSampleSize; - residualForcesMovingAverageHistorySample.pop(); - pMesh->residualForcesMovingAverageDerivativeNorm - = std::abs((residualForcesMovingAverageHistorySample.back() - - residualForcesMovingAverageHistorySample.front())) - / (movingAverageSampleSize); - residualForcesMovingAverageDerivativeMax - = std::max(pMesh->residualForcesMovingAverageDerivativeNorm, - residualForcesMovingAverageDerivativeMax); - pMesh->residualForcesMovingAverageDerivativeNorm - /= residualForcesMovingAverageDerivativeMax; - // std::cout << "Normalized derivative:" - // << residualForcesMovingAverageDerivativeNorm - // << std::endl; - } + // double sumOfDisplacementsNorm = 0; + // for (size_t vi = 0; vi < pMesh->VN(); vi++) { + // sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm(); + // } + // sumOfDisplacementsNorm /= pMesh->bbox.Diag(); + // pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm; + // // pMesh->residualForcesMovingAverage = (pMesh->totalResidualForcesNorm + // // + mCurrentSimulationStep + // // * pMesh->residualForcesMovingAverage) + // // / (1 + mCurrentSimulationStep); + // pMesh->residualForcesMovingAverage += pMesh->totalResidualForcesNorm + // / movingAverageSampleSize; + // residualForcesMovingAverageHistorySample.push(pMesh->residualForcesMovingAverage); + // if (movingAverageSampleSize < residualForcesMovingAverageHistorySample.size()) { + // const double firstElementValue = residualForcesMovingAverageHistorySample.front(); + // pMesh->residualForcesMovingAverage -= firstElementValue / movingAverageSampleSize; + // residualForcesMovingAverageHistorySample.pop(); + + // pMesh->residualForcesMovingAverageDerivativeNorm + // = std::abs((residualForcesMovingAverageHistorySample.back() + // - residualForcesMovingAverageHistorySample.front())) + // / (movingAverageSampleSize); + // residualForcesMovingAverageDerivativeMax + // = std::max(pMesh->residualForcesMovingAverageDerivativeNorm, + // residualForcesMovingAverageDerivativeMax); + // pMesh->residualForcesMovingAverageDerivativeNorm + // /= residualForcesMovingAverageDerivativeMax; + // // std::cout << "Normalized derivative:" + // // << residualForcesMovingAverageDerivativeNorm + // // << std::endl; + // } // pMesh->previousTotalPotentialEnergykN = // pMesh->currentTotalPotentialEnergykN; @@ -2127,10 +2146,21 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrtotalResidualForcesNorm - << "%" << std::endl; + + // auto t2 = std::chrono::high_resolution_clock::now(); + // const double executionTime_min + // = std::chrono::duration_cast(t2 - beginTime).count(); + // std::cout << "Execution time(min):" << executionTime_min << std::endl; + if (mSettings.useAverageResidualForcesNorm) { + const double percOfConv = averageExternalForcesNorm + / pMesh->averageResidualForcesNorm; + std::cout << "Percentage of target:" << percOfConv << "%" << std::endl; + } else { + std::cout << "Percentage of target:" + << 100 * mSettings.totalResidualForcesNormThreshold + / pMesh->totalResidualForcesNorm + << "%" << std::endl; + } // SimulationResultsReporter::createPlot("Number of Steps", // "Residual Forces mov aver", // movingAverages); @@ -2149,13 +2179,6 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrtotalResidualForcesNorm); } - if (std::isnan(pMesh->currentTotalKineticEnergy)) { - if (mSettings.beVerbose) { - std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl; - } - break; - } - if (mSettings.shouldCreatePlots && mSettings.isDebugMode && mCurrentSimulationStep % mSettings.debugModeStep == 0) { // SimulationResultsReporter::createPlot("Number of Steps", @@ -2164,18 +2187,18 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrpMesh->getLabel()); + // SimulationResultsReporter reporter; + // reporter.reportHistory(history, "IntermediateResults", pJob->pMesh->getLabel()); // SimulationResultsReporter::createPlot("Number of Iterations", // "Sum of normalized displacement norms", // history.sumOfNormalizedDisplacementNorms /*, @@ -2234,8 +2257,12 @@ mesh->currentTotalPotentialEnergykN*/ && pMesh->currentTotalTranslationalKineticEnergy < mSettings.totalTranslationalKineticEnergyThreshold && mCurrentSimulationStep > 20 && numOfDampings > 0; + const bool fullfillsAverageResidualForcesNormCriterion + = mSettings.useAverageResidualForcesNorm + && pMesh->averageResidualForcesNorm / averageExternalForcesNorm < 1e-2; const bool fullfillsResidualForcesNormTerminationCriterion - = pMesh->totalResidualForcesNorm < mSettings.totalResidualForcesNormThreshold; + = mSettings.useResidualForcesNorm + && pMesh->totalResidualForcesNorm < mSettings.totalResidualForcesNormThreshold; const bool fullfillsMovingAverageNormTerminationCriterion = pMesh->residualForcesMovingAverage < mSettings.residualForcesMovingAverageNormThreshold; @@ -2246,7 +2273,8 @@ mesh->currentTotalPotentialEnergykN*/ = fullfillsKineticEnergyTerminationCriterion // || fullfillsMovingAverageNormTerminationCriterion // || fullfillsMovingAverageDerivativeNormTerminationCriterion - || fullfillsResidualForcesNormTerminationCriterion; + || fullfillsResidualForcesNormTerminationCriterion + || fullfillsAverageResidualForcesNormCriterion; if (shouldTerminate) { if (mSettings.beVerbose && !mSettings.isDebugMode) { std::cout << "Simulation converged." << std::endl; @@ -2265,29 +2293,52 @@ mesh->currentTotalPotentialEnergykN*/ << std::endl; } } - break; + if (mSettings.desiredGradualExternalLoadsSteps == externalLoadStep) { + break; + } else { + externalLoadStep++; + std::unordered_map nodalExternalForces + = pJob->nodalExternalForces; + double totalExternalForcesNorm = 0; + for (auto &nodalForce : nodalExternalForces) { + const double percentageOfExternalLoads + = double(externalLoadStep) / mSettings.desiredGradualExternalLoadsSteps; + nodalForce.second = nodalForce.second * percentageOfExternalLoads; + totalExternalForcesNorm += nodalForce.second.norm(); + } + updateNodalExternalForces(nodalExternalForces, constrainedVertices); + + averageExternalForcesNorm = totalExternalForcesNorm / pMesh->VN(); + if (!nodalExternalForces.empty()) { + mSettings.totalResidualForcesNormThreshold = 1e-2 * totalExternalForcesNorm; + } + if (mSettings.beVerbose) { + std::cout << "totalResidualForcesNormThreshold:" + << mSettings.totalResidualForcesNormThreshold << std::endl; + } + } // } } - // const bool shouldCapDisplacements = mSettings.displacementCap.has_value(); - // for (VertexType &v : pMesh->vert) { - // Node &node = pMesh->nodes[v]; - // Vector6d stepDisplacement = node.velocity * Dt; - // if (shouldCapDisplacements - // && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) { - // stepDisplacement = stepDisplacement - // * (*mSettings.displacementCap - // / stepDisplacement.getTranslation().norm()); - // } - // node.displacements = node.displacements - stepDisplacement; - // } - // applyDisplacements(constrainedVertices); - // if (!pJob->nodalForcedDisplacements.empty()) { - // applyForcedDisplacements( + const bool shouldCapDisplacements = mSettings.displacementCap.has_value(); + for (VertexType &v : pMesh->vert) { + Node &node = pMesh->nodes[v]; + Vector6d stepDisplacement = node.velocity * Dt; + if (shouldCapDisplacements + && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) { + stepDisplacement = stepDisplacement + * (*mSettings.displacementCap + / stepDisplacement.getTranslation().norm()); + } + node.displacements = node.displacements - stepDisplacement; + } + applyDisplacements(constrainedVertices); + if (!pJob->nodalForcedDisplacements.empty()) { + applyForcedDisplacements( - // pJob->nodalForcedDisplacements); - // } - // updateElementalLengths(); + pJob->nodalForcedDisplacements); + } + updateElementalLengths(); // const double triggerPercentage = 0.01; // const double xi_min = 0.55; @@ -2313,8 +2364,11 @@ mesh->currentTotalPotentialEnergykN*/ // draw(); } } + auto endTime = std::chrono::high_resolution_clock::now(); SimulationResults results = computeResults(pJob); + results.executionTime = std::chrono::duration_cast(endTime - beginTime) + .count(); if (mCurrentSimulationStep == settings.maxDRMIterations && mCurrentSimulationStep != 0) { if (mSettings.beVerbose) { diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index af43ad7..96422bf 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -39,10 +39,13 @@ public: int maxDRMIterations{0}; bool shouldUseTranslationalKineticEnergyThreshold{false}; int gradualForcedDisplacementSteps{100}; + int desiredGradualExternalLoadsSteps{1}; double gamma{0.8}; std::optional displacementCap; double totalResidualForcesNormThreshold{1e-3}; double totalExternalForcesNormPercentageTermination{1e-3}; + bool useResidualForcesNorm{true}; + bool useAverageResidualForcesNorm{false}; Settings() {} }; @@ -57,6 +60,7 @@ private: matplot::line_handle plotHandle; std::vector plotYValues; size_t numOfDampings{0}; + int externalLoadStep{1}; const std::string meshPolyscopeLabel{"Simulation mesh"}; std::unique_ptr pMesh; diff --git a/reducedmodeloptimizer_structs.hpp b/reducedmodeloptimizer_structs.hpp index 40878ce..6d42852 100644 --- a/reducedmodeloptimizer_structs.hpp +++ b/reducedmodeloptimizer_structs.hpp @@ -267,6 +267,7 @@ struct Colors std::vector perSimulationScenario_rotational; std::vector perSimulationScenario_total; } objectiveValue; + std::vector perScenario_fullPatternPotentialEnergy; // std::vector> optimalXNameValuePairs; std::vector> optimalXNameValuePairs; std::vector> fullPatternSimulationJobs; @@ -284,6 +285,7 @@ struct Colors inline static std::string Label{"Label"}; inline static std::string FullPatternLabel{"FullPatternLabel"}; inline static std::string Settings{"OptimizationSettings"}; + inline static std::string FullPatternPotentialEnergies{"PE"}; }; void save(const string &saveToPath, const bool shouldExportDebugFiles = false) @@ -323,6 +325,16 @@ struct Colors baseTriangle.cP2(0)[2]}; baseTriangleFullPattern.save(std::filesystem::path(saveToPath).string()); json_optimizationResults[JsonKeys::FullPatternLabel] = baseTriangleFullPattern.getLabel(); + + //potential energies + const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); + std::vector fullPatternPE(numberOfSimulationJobs); + for (int simulationScenarioIndex = 0; simulationScenarioIndex < numberOfSimulationJobs; + simulationScenarioIndex++) { + fullPatternPE[simulationScenarioIndex] + = perScenario_fullPatternPotentialEnergy[simulationScenarioIndex]; + } + json_optimizationResults[JsonKeys::FullPatternPotentialEnergies] = fullPatternPE; ////Save to json file std::filesystem::path jsonFilePath( std::filesystem::path(saveToPath).append(JsonKeys::filename)); @@ -652,6 +664,13 @@ struct Colors os << "Obj value Rot " + simulationScenarioName; os << "Obj value Total " + simulationScenarioName; } + for (int simulationScenarioIndex = 0; + simulationScenarioIndex < fullPatternSimulationJobs.size(); + simulationScenarioIndex++) { + const std::string simulationScenarioName + = fullPatternSimulationJobs[simulationScenarioIndex]->getLabel(); + os << "PE " + simulationScenarioName; + } for (const auto &nameValuePair : optimalXNameValuePairs) { os << nameValuePair.first; } @@ -671,6 +690,11 @@ struct Colors os << objectiveValue.perSimulationScenario_rotational[simulationScenarioIndex]; os << objectiveValue.perSimulationScenario_total[simulationScenarioIndex]; } + for (int simulationScenarioIndex = 0; + simulationScenarioIndex < fullPatternSimulationJobs.size(); + simulationScenarioIndex++) { + os << perScenario_fullPatternPotentialEnergy[simulationScenarioIndex]; + } for (const auto &optimalXNameValuePair : optimalXNameValuePairs) { os << optimalXNameValuePair.second; } diff --git a/simulation_structs.hpp b/simulation_structs.hpp index cf515cd..f27e8be 100755 --- a/simulation_structs.hpp +++ b/simulation_structs.hpp @@ -492,7 +492,7 @@ struct SimulationResults return labelPrefix + deliminator + job->pMesh->getLabel() + deliminator + job->getLabel(); } - void saveDeformedModel(const std::string &outputFolder = std::string()) + bool saveDeformedModel(const std::string &outputFolder = std::string()) { VCGEdgeMesh m; vcg::tri::Append::MeshCopy(m, *job->pMesh); @@ -503,7 +503,7 @@ struct SimulationResults displacements[vi][1], displacements[vi][2]); } - m.save(std::filesystem::path(outputFolder).append(getLabel() + ".ply").string()); + return m.save(std::filesystem::path(outputFolder).append(getLabel() + ".ply").string()); } void save(const std::string &outputFolder = std::string()) diff --git a/simulationmesh.hpp b/simulationmesh.hpp index ca3dc38..150f8d3 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -43,6 +43,7 @@ public: double currentTotalKineticEnergy{0}; double currentTotalTranslationalKineticEnergy{0}; double totalResidualForcesNorm{0}; + double averageResidualForcesNorm{0}; double currentTotalPotentialEnergykN{0}; double previousTotalPotentialEnergykN{0}; double residualForcesMovingAverageDerivativeNorm{0}; diff --git a/utilities.hpp b/utilities.hpp index cf13769..863768a 100755 --- a/utilities.hpp +++ b/utilities.hpp @@ -116,6 +116,16 @@ struct Vector6d : public std::array { { return Eigen::Vector3d(this->operator[](3), this->operator[](4), this->operator[](5)); } + + std::string toString() const + { + std::string s; + for (int i = 0; i < 6; i++) { + s.append(std::to_string(this->operator[](i)) + ","); + } + s.pop_back(); + return s; + } }; namespace Utilities {