manolas
/
MySources
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Added tbb for parallel drm. Refactoring. DRM settings save and load functions

This commit is contained in:
iasonmanolas 2021-07-22 18:06:46 +03:00
parent 90a00c5415
commit 8aff310a83
6 changed files with 323 additions and 172 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@ -50,6 +50,7 @@ download_project(PROJ TBB
PREFIX ${EXTERNAL_DEPS_DIR} PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE} ${UPDATE_DISCONNECTED_IF_AVAILABLE}
) )
option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF)
add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR}) add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR})
link_directories(${TBB_BINARY_DIR}) link_directories(${TBB_BINARY_DIR})
message(${TBB_BINARY_DIR}) message(${TBB_BINARY_DIR})

424
drmsimulationmodel.cpp
View File

@ -1,4 +1,5 @@
#include "drmsimulationmodel.hpp" #include "drmsimulationmodel.hpp"
#include "linearsimulationmodel.hpp"
#include "simulationhistoryplotter.hpp" #include "simulationhistoryplotter.hpp"
#include <Eigen/Geometry> #include <Eigen/Geometry>
#include <algorithm> #include <algorithm>
@ -810,7 +811,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
//#pragma omp parallel for schedule(static) num_threads(5) //#pragma omp parallel for schedule(static) num_threads(5)
//#endif //#endif
std::for_each( std::for_each(
#ifdef ENABLE_PARALLEL_DRM
std::execution::par_unseq, std::execution::par_unseq,
#endif
pMesh->edge.begin(), pMesh->edge.begin(),
pMesh->edge.end(), pMesh->edge.end(),
[&](const EdgeType &e) [&](const EdgeType &e)
@ -1049,6 +1052,7 @@ void DRMSimulationModel::updateNodalExternalForces(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices) const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices)
{ {
externalMomentsNorm = 0; externalMomentsNorm = 0;
double totalExternalForcesNorm = 0;
for (const std::pair<VertexIndex, Vector6d> &nodalForce : nodalForces) { for (const std::pair<VertexIndex, Vector6d> &nodalForce : nodalForces) {
const VertexIndex nodeIndex = nodalForce.first; const VertexIndex nodeIndex = nodalForce.first;
const bool isNodeConstrained = fixedVertices.contains(nodeIndex); const bool isNodeConstrained = fixedVertices.contains(nodeIndex);
@ -1086,7 +1090,10 @@ void DRMSimulationModel::updateNodalExternalForces(
} }
node.force.external = nodalExternalForce; node.force.external = nodalExternalForce;
totalExternalForcesNorm += node.force.external.norm();
} }
pMesh->totalExternalForcesNorm = totalExternalForcesNorm;
} }
void DRMSimulationModel::updateResidualForces() void DRMSimulationModel::updateResidualForces()
@ -1242,6 +1249,13 @@ void DRMSimulationModel::updateNodalMasses()
* rotationalSumSk_I3; * rotationalSumSk_I3;
pMesh->nodes[v].mass.rotationalJ = gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_J; pMesh->nodes[v].mass.rotationalJ = gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_J;
pMesh->nodes[v].mass_6d[DoF::Ux] = pMesh->nodes[v].mass.translational;
pMesh->nodes[v].mass_6d[DoF::Uy] = pMesh->nodes[v].mass.translational;
pMesh->nodes[v].mass_6d[DoF::Uz] = pMesh->nodes[v].mass.translational;
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;
assert(std::pow(mSettings.Dtini, 2.0) * translationalSumSk assert(std::pow(mSettings.Dtini, 2.0) * translationalSumSk
/ pMesh->nodes[v].mass.translational / pMesh->nodes[v].mass.translational
< 2); < 2);
@ -1295,17 +1309,17 @@ void DRMSimulationModel::updateNodalVelocities()
void DRMSimulationModel::updateNodalDisplacements() void DRMSimulationModel::updateNodalDisplacements()
{ {
const bool shouldCapDisplacements = mSettings.displacementCap.has_value(); // const bool shouldCapDisplacements = mSettings.displacementCap.has_value();
for (VertexType &v : pMesh->vert) { for (VertexType &v : pMesh->vert) {
Node &node = pMesh->nodes[v]; Node &node = pMesh->nodes[v];
Vector6d stepDisplacement = node.velocity * Dt; Vector6d stepDisplacement = node.velocity * Dt;
if (shouldCapDisplacements // if (shouldCapDisplacements
&& stepDisplacement.getTranslation().norm() > mSettings.displacementCap) { // && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) {
stepDisplacement = stepDisplacement // stepDisplacement = stepDisplacement
* (*mSettings.displacementCap // * (*mSettings.displacementCap
/ stepDisplacement.getTranslation().norm()); // / stepDisplacement.getTranslation().norm());
std::cout << "Displacement capped" << std::endl; // std::cout << "Displacement capped" << std::endl;
} // }
node.displacements = node.displacements + stepDisplacement; node.displacements = node.displacements + stepDisplacement;
// if (mSettings.isDebugMode && mSettings.beVerbose && pMesh->getIndex(v) == 40) { // if (mSettings.isDebugMode && mSettings.beVerbose && pMesh->getIndex(v) == 40) {
// std::cout << "Node " << node.vi << ":" << endl; // std::cout << "Node " << node.vi << ":" << endl;
@ -1397,7 +1411,7 @@ void DRMSimulationModel::updateNodeNormal(
* */ * */
const bool viHasMoments = node.force.external[3] != 0 || node.force.external[4] != 0; const bool viHasMoments = node.force.external[3] != 0 || node.force.external[4] != 0;
if (!checkedForMaximumMoment && viHasMoments) { if (!checkedForMaximumMoment && viHasMoments) {
mSettings.shouldUseTranslationalKineticEnergyThreshold = true; mSettings.totalTranslationalKineticEnergyThreshold = 1e-8;
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
std::cout << "Maximum moment reached.The Kinetic energy of the system will " std::cout << "Maximum moment reached.The Kinetic energy of the system will "
"be used as a convergence criterion" "be used as a convergence criterion"
@ -1830,9 +1844,14 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
// instead of full width. Must have // instead of full width. Must have
// matching PopItemWidth() below. // matching PopItemWidth() below.
ImGui::InputInt("Simulation debug step", static int debugModeStep = 0;
&mSettings.debugModeStep); // set a int variable if (ImGui::InputInt("Simulation debug step",
mSettings.drawingStep = mSettings.debugModeStep; &debugModeStep)) { // set a int variable
if (debugModeStep != 0) {
*mSettings.debugModeStep = debugModeStep;
}
}
// mSettings.drawingStep = mSettings.debugModeStep;
ImGui::Checkbox("Enable drawing", ImGui::Checkbox("Enable drawing",
&mSettings.shouldDraw); // set a int variable &mSettings.shouldDraw); // set a int variable
ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep); ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep);
@ -2024,6 +2043,42 @@ void DRMSimulationModel::applySolutionGuess(const SimulationResults &solutionGue
// //} // //}
// polyscope::show(); // polyscope::show();
} }
void DRMSimulationModel::applyKineticDamping(const std::shared_ptr<SimulationJob> &pJob)
{
// 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;
}
applyDisplacements(constrainedVertices);
if (!pJob->nodalForcedDisplacements.empty()) {
applyForcedDisplacements(
pJob->nodalForcedDisplacements);
}
updateElementalLengths();
// const double triggerPercentage = 0.01;
// const double xi_min = 0.55;
// const double xi_init = 0.9969;
// if (mSettings.totalResidualForcesNormThreshold / pMesh->totalResidualForcesNorm
// > triggerPercentage) {
// mSettings.xi = ((xi_min - xi_init)
// * (mSettings.totalResidualForcesNormThreshold
// / pMesh->totalResidualForcesNorm)
// + xi_init - triggerPercentage * xi_min)
// / (1 - triggerPercentage);
// }
resetVelocities();
++numOfDampings;
}
SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob, SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const Settings &settings, const Settings &settings,
@ -2072,7 +2127,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
} }
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw || mSettings.isDebugMode) { if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) {
PolyscopeInterface::init(); PolyscopeInterface::init();
polyscope::registerCurveNetwork(meshPolyscopeLabel, polyscope::registerCurveNetwork(meshPolyscopeLabel,
pMesh->getEigenVertices(), pMesh->getEigenVertices(),
@ -2090,35 +2145,44 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
} }
updateElementalFrames(); updateElementalFrames();
if (!solutionGuess.displacements.empty()) { if (!solutionGuess.displacements.empty()) {
//#ifdef POLYSCOPE_DEFINED assert(!mSettings.linearGuessForceScaleFactor.has_value());
// if (mSettings.shouldDraw || mSettings.isDebugMode) {
// solutionGuess.registerForDrawing();
// polyscope::show();
// solutionGuess.unregister();
// }
//#endif
applySolutionGuess(solutionGuess, pJob); applySolutionGuess(solutionGuess, pJob);
shouldTemporarilyDampForces = true; shouldTemporarilyDampForces = true;
// Dt = mSettings.Dtini * 0.9; }
if (mSettings.linearGuessForceScaleFactor.has_value()) {
const double forceScaleFactor = 1;
for (auto &externalForce : pJob->nodalExternalForces) {
externalForce.second = externalForce.second * forceScaleFactor;
}
LinearSimulationModel linearSimulationModel;
SimulationResults simulationResults_fullPatternLinearModel = linearSimulationModel
.executeSimulation(pJob);
// simulationResults_fullPatternLinearModel.save(fullResultsFolderPath);
for (auto &externalForce : pJob->nodalExternalForces) {
externalForce.second = externalForce.second / forceScaleFactor;
}
applySolutionGuess(simulationResults_fullPatternLinearModel, pJob);
shouldTemporarilyDampForces = true;
} }
updateNodalMasses(); updateNodalMasses();
std::unordered_map<VertexIndex, Vector6d> nodalExternalForces = pJob->nodalExternalForces; // std::unordered_map<VertexIndex, Vector6d> nodalExternalForces = pJob->nodalExternalForces;
double totalExternalForcesNorm = 0; // double totalExternalForcesNorm = 0;
// Vector6d sumOfExternalForces(0); // Vector6d sumOfExternalForces(0);
for (auto &nodalForce : nodalExternalForces) { // for (auto &nodalForce : nodalExternalForces) {
const double percentageOfExternalLoads = double(externalLoadStep) // const double percentageOfExternalLoads = double(externalLoadStep)
/ mSettings.desiredGradualExternalLoadsSteps; // / mSettings.desiredGradualExternalLoadsSteps;
nodalForce.second = nodalForce.second * percentageOfExternalLoads; // nodalForce.second = nodalForce.second * percentageOfExternalLoads;
totalExternalForcesNorm += nodalForce.second.norm(); // totalExternalForcesNorm += nodalForce.second.norm();
// sumOfExternalForces = sumOfExternalForces + nodalForce.second; // // sumOfExternalForces = sumOfExternalForces + nodalForce.second;
} // }
pMesh->totalExternalForcesNorm = totalExternalForcesNorm; updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
updateNodalExternalForces(nodalExternalForces, constrainedVertices); if (!pJob->nodalExternalForces.empty()) {
if (!nodalExternalForces.empty()) {
mSettings.totalResidualForcesNormThreshold mSettings.totalResidualForcesNormThreshold
= settings.totalExternalForcesNormPercentageTermination * totalExternalForcesNorm; = settings.totalExternalForcesNormPercentageTermination
* pMesh->totalExternalForcesNorm;
} else { } else {
mSettings.totalResidualForcesNormThreshold = 1e-3; mSettings.totalResidualForcesNormThreshold = 1e-3;
std::cout << "No forces setted default residual forces norm threshold" << std::endl; std::cout << "No forces setted default residual forces norm threshold" << std::endl;
@ -2126,9 +2190,9 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
if (mSettings.beVerbose) { if (mSettings.beVerbose) {
std::cout << "totalResidualForcesNormThreshold:" std::cout << "totalResidualForcesNormThreshold:"
<< mSettings.totalResidualForcesNormThreshold << std::endl; << mSettings.totalResidualForcesNormThreshold << std::endl;
if (mSettings.useAverage) { if (mSettings.averageResidualForcesCriterionThreshold.has_value()) {
std::cout << "average/extNorm threshold:" std::cout << "average/extNorm threshold:"
<< mSettings.averageResidualForcesCriterionThreshold << std::endl; << *mSettings.averageResidualForcesCriterionThreshold << std::endl;
} }
} }
const size_t movingAverageSampleSize = 200; const size_t movingAverageSampleSize = 200;
@ -2136,8 +2200,9 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
std::vector<double> percentageOfConvergence; std::vector<double> percentageOfConvergence;
// double residualForcesMovingAverageDerivativeMax = 0; // double residualForcesMovingAverageDerivativeMax = 0;
while (settings.maxDRMIterations == 0 || mCurrentSimulationStep < settings.maxDRMIterations) { while (!settings.maxDRMIterations.has_value()
if ((mSettings.isDebugMode && mCurrentSimulationStep == 50000)) { || mCurrentSimulationStep < settings.maxDRMIterations.value()) {
if ((mSettings.debugModeStep.has_value() && mCurrentSimulationStep == 50000)) {
// std::filesystem::create_directory("./PatternOptimizationNonConv"); // std::filesystem::create_directory("./PatternOptimizationNonConv");
// pJob->save("./PatternOptimizationNonConv"); // pJob->save("./PatternOptimizationNonConv");
// Dt = mSettings.Dtini; // Dt = mSettings.Dtini;
@ -2221,25 +2286,26 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// pMesh->currentTotalPotentialEnergykN; // pMesh->currentTotalPotentialEnergykN;
// pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000; // pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000;
// timePerNodePerIterationHistor.push_back(timePerNodePerIteration); // timePerNodePerIterationHistor.push_back(timePerNodePerIteration);
if (mSettings.beVerbose && mSettings.isDebugMode if (mSettings.beVerbose && mSettings.debugModeStep.has_value()
&& mCurrentSimulationStep % mSettings.debugModeStep == 0) { && mCurrentSimulationStep % mSettings.debugModeStep.value() == 0) {
printCurrentState(); printCurrentState();
// auto t2 = std::chrono::high_resolution_clock::now(); // auto t2 = std::chrono::high_resolution_clock::now();
// const double executionTime_min // const double executionTime_min
// = std::chrono::duration_cast<std::chrono::minutes>(t2 - beginTime).count(); // = std::chrono::duration_cast<std::chrono::minutes>(t2 - beginTime).count();
// std::cout << "Execution time(min):" << executionTime_min << std::endl; // std::cout << "Execution time(min):" << executionTime_min << std::endl;
if (mSettings.useAverage) { if (mSettings.averageResidualForcesCriterionThreshold.has_value()) {
std::cout << "Best percentage of target (average):" std::cout << "Best percentage of target (average):"
<< 100 * mSettings.averageResidualForcesCriterionThreshold << 100 * mSettings.averageResidualForcesCriterionThreshold.value()
* totalExternalForcesNorm * pMesh->totalExternalForcesNorm * pMesh->VN()
/ (minTotalResidualForcesNorm / pMesh->VN()) / minTotalResidualForcesNorm
<< "%" << std::endl; << "%" << std::endl;
} } else {
std::cout << "Best percentage of target:" std::cout << "Best percentage of target:"
<< 100 * mSettings.totalExternalForcesNormPercentageTermination << 100 * mSettings.totalExternalForcesNormPercentageTermination
* totalExternalForcesNorm / minTotalResidualForcesNorm * pMesh->totalExternalForcesNorm / minTotalResidualForcesNorm
<< "%" << std::endl; << "%" << std::endl;
}
// SimulationResultsReporter::createPlot("Number of Steps", // SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver", // "Residual Forces mov aver",
// history.residualForcesMovingAverage); // history.residualForcesMovingAverage);
@ -2252,14 +2318,15 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// timePerNodePerIterationHistory); // timePerNodePerIterationHistory);
} }
if ((mSettings.shouldCreatePlots || mSettings.isDebugMode) && mCurrentSimulationStep != 0) { if ((mSettings.shouldCreatePlots || mSettings.debugModeStep.has_value())
&& mCurrentSimulationStep != 0) {
history.stepPulse(*pMesh); history.stepPulse(*pMesh);
percentageOfConvergence.push_back(100 * mSettings.totalResidualForcesNormThreshold percentageOfConvergence.push_back(100 * mSettings.totalResidualForcesNormThreshold
/ pMesh->totalResidualForcesNorm); / pMesh->totalResidualForcesNorm);
} }
if (mSettings.shouldCreatePlots && mSettings.isDebugMode if (mSettings.shouldCreatePlots && mSettings.debugModeStep.has_value()
&& mCurrentSimulationStep % mSettings.debugModeStep == 0) { && mCurrentSimulationStep % mSettings.debugModeStep.value() == 0) {
// SimulationResultsReporter::createPlot("Number of Steps", // SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver deriv", // "Residual Forces mov aver deriv",
// movingAveragesDerivatives_norm); // movingAveragesDerivatives_norm);
@ -2295,8 +2362,8 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
} }
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw && mSettings.isDebugMode if (mSettings.shouldDraw && mSettings.debugModeStep.has_value()
&& mCurrentSimulationStep % mSettings.drawingStep == 0) /* && && mCurrentSimulationStep % mSettings.debugModeStep.value() == 0) /* &&
currentSimulationStep > maxDRMIterations*/ currentSimulationStep > maxDRMIterations*/
{ {
@ -2320,20 +2387,35 @@ currentSimulationStep > maxDRMIterations*/
// std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm
// << std::endl; // << std::endl;
const bool fullfillsResidualForcesNormTerminationCriterion const bool fullfillsResidualForcesNormTerminationCriterion
= !mSettings.useAverage = !mSettings.averageResidualForcesCriterionThreshold.has_value()
&& pMesh->totalResidualForcesNorm / totalExternalForcesNorm && pMesh->totalResidualForcesNorm / pMesh->totalExternalForcesNorm
< mSettings.totalExternalForcesNormPercentageTermination; < mSettings.totalExternalForcesNormPercentageTermination;
const bool fullfillsAverageResidualForcesNormTerminationCriterion const bool fullfillsAverageResidualForcesNormTerminationCriterion
= mSettings.useAverage = mSettings.averageResidualForcesCriterionThreshold.has_value()
&& (pMesh->totalResidualForcesNorm / pMesh->VN()) / totalExternalForcesNorm && (pMesh->totalResidualForcesNorm / pMesh->VN()) / pMesh->totalExternalForcesNorm
< mSettings.averageResidualForcesCriterionThreshold; < mSettings.averageResidualForcesCriterionThreshold.value()
// Residual forces norm convergence && numOfDampings > 2;
if (((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy if ((fullfillsAverageResidualForcesNormTerminationCriterion
|| fullfillsAverageResidualForcesNormTerminationCriterion
|| fullfillsResidualForcesNormTerminationCriterion) || fullfillsResidualForcesNormTerminationCriterion)
// && mCurrentSimulationStep > movingAverageSampleSize
&& (pJob->nodalForcedDisplacements.empty() && (pJob->nodalForcedDisplacements.empty()
|| mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)) || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)) {
if (mSettings.beVerbose /*&& !mSettings.isDebugMode*/) {
std::cout << "Simulation converged." << std::endl;
printCurrentState();
if (fullfillsResidualForcesNormTerminationCriterion) {
std::cout << "Converged using residual forces norm threshold criterion"
<< std::endl;
} else if (fullfillsAverageResidualForcesNormTerminationCriterion) {
std::cout << "Converged using average residual forces norm threshold criterion"
<< std::endl;
}
}
break;
}
// Residual forces norm convergence
if ((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy
// && mCurrentSimulationStep > movingAverageSampleSize
)
/* || (pMesh->totalResidualForcesNorm / mSettings.totalResidualForcesNormThreshold <= 1 /* || (pMesh->totalResidualForcesNorm / mSettings.totalResidualForcesNormThreshold <= 1
&& mCurrentSimulationStep > 1)*/ && mCurrentSimulationStep > 1)*/
/*|| /*||
@ -2342,103 +2424,61 @@ mesh->currentTotalPotentialEnergykN*/
/*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) { /*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) {
// if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) { // if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) {
const bool fullfillsKineticEnergyTerminationCriterion const bool fullfillsKineticEnergyTerminationCriterion
= mSettings.shouldUseTranslationalKineticEnergyThreshold = mSettings.totalTranslationalKineticEnergyThreshold.has_value()
&& pMesh->currentTotalTranslationalKineticEnergy && pMesh->currentTotalTranslationalKineticEnergy
< mSettings.totalTranslationalKineticEnergyThreshold < mSettings.totalTranslationalKineticEnergyThreshold.value()
&& mCurrentSimulationStep > 20 && numOfDampings > 0; && mCurrentSimulationStep > 20 && numOfDampings > 0;
const bool fullfillsMovingAverageNormTerminationCriterion // const bool fullfillsMovingAverageNormTerminationCriterion
= pMesh->residualForcesMovingAverage // = pMesh->residualForcesMovingAverage
< mSettings.residualForcesMovingAverageNormThreshold; // < mSettings.residualForcesMovingAverageNormThreshold;
/*pMesh->residualForcesMovingAverage < totalResidualForcesNormThreshold*/ /*pMesh->residualForcesMovingAverage < totalResidualForcesNormThreshold*/
const bool fullfillsMovingAverageDerivativeNormTerminationCriterion // const bool fullfillsMovingAverageDerivativeNormTerminationCriterion
= pMesh->residualForcesMovingAverageDerivativeNorm < 1e-4; // = pMesh->residualForcesMovingAverageDerivativeNorm < 1e-4;
const bool shouldTerminate const bool shouldTerminate = fullfillsKineticEnergyTerminationCriterion
= fullfillsKineticEnergyTerminationCriterion
// || fullfillsMovingAverageNormTerminationCriterion // || fullfillsMovingAverageNormTerminationCriterion
// || fullfillsMovingAverageDerivativeNormTerminationCriterion // || fullfillsMovingAverageDerivativeNormTerminationCriterion
|| fullfillsAverageResidualForcesNormTerminationCriterion ;
|| fullfillsResidualForcesNormTerminationCriterion; if (shouldTerminate) {
if (shouldTerminate && mCurrentSimulationStep > 100) {
if (mSettings.beVerbose /*&& !mSettings.isDebugMode*/) { if (mSettings.beVerbose /*&& !mSettings.isDebugMode*/) {
std::cout << "Simulation converged." << std::endl; std::cout << "Simulation converged." << std::endl;
printCurrentState(); printCurrentState();
if (fullfillsResidualForcesNormTerminationCriterion) { if (fullfillsKineticEnergyTerminationCriterion) {
std::cout << "Converged using residual forces norm threshold criterion"
<< std::endl;
} else if (fullfillsKineticEnergyTerminationCriterion) {
std::cout << "The kinetic energy of the system was " std::cout << "The kinetic energy of the system was "
" used as a convergence criterion" " used as a convergence criterion"
<< std::endl; << std::endl;
} else if (fullfillsMovingAverageNormTerminationCriterion) {
std::cout
<< "Converged using residual forces moving average derivative norm "
"threshold criterion"
<< std::endl;
} }
} }
if (mSettings.desiredGradualExternalLoadsSteps == externalLoadStep) { // if (mSettings.desiredGradualExternalLoadsSteps == externalLoadStep) {
break; break;
} else { // } else {
externalLoadStep++; // externalLoadStep++;
std::unordered_map<VertexIndex, Vector6d> nodalExternalForces // std::unordered_map<VertexIndex, Vector6d> nodalExternalForces
= pJob->nodalExternalForces; // = pJob->nodalExternalForces;
double totalExternalForcesNorm = 0; // double totalExternalForcesNorm = 0;
for (auto &nodalForce : nodalExternalForces) { // for (auto &nodalForce : nodalExternalForces) {
const double percentageOfExternalLoads // const double percentageOfExternalLoads
= double(externalLoadStep) / mSettings.desiredGradualExternalLoadsSteps; // = double(externalLoadStep) / mSettings.desiredGradualExternalLoadsSteps;
nodalForce.second = nodalForce.second * percentageOfExternalLoads; // nodalForce.second = nodalForce.second * percentageOfExternalLoads;
totalExternalForcesNorm += nodalForce.second.norm(); // totalExternalForcesNorm += nodalForce.second.norm();
} // }
updateNodalExternalForces(nodalExternalForces, constrainedVertices); // updateNodalExternalForces(nodalExternalForces, constrainedVertices);
if (!nodalExternalForces.empty()) { // if (!nodalExternalForces.empty()) {
mSettings.totalResidualForcesNormThreshold = 1e-2 * totalExternalForcesNorm; // mSettings.totalResidualForcesNormThreshold = 1e-2 * totalExternalForcesNorm;
} // }
if (mSettings.beVerbose) { // if (mSettings.beVerbose) {
std::cout << "totalResidualForcesNormThreshold:" // std::cout << "totalResidualForcesNormThreshold:"
<< mSettings.totalResidualForcesNormThreshold << std::endl; // << mSettings.totalResidualForcesNormThreshold << std::endl;
} // }
} // }
// } // }
} }
const bool shouldCapDisplacements = mSettings.displacementCap.has_value(); applyKineticDamping(pJob);
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;
}
applyDisplacements(constrainedVertices);
if (!pJob->nodalForcedDisplacements.empty()) {
applyForcedDisplacements(
pJob->nodalForcedDisplacements);
}
updateElementalLengths();
// const double triggerPercentage = 0.01;
// const double xi_min = 0.55;
// const double xi_init = 0.9969;
// if (mSettings.totalResidualForcesNormThreshold / pMesh->totalResidualForcesNorm
// > triggerPercentage) {
// mSettings.xi = ((xi_min - xi_init)
// * (mSettings.totalResidualForcesNormThreshold
// / pMesh->totalResidualForcesNorm)
// + xi_init - triggerPercentage * xi_min)
// / (1 - triggerPercentage);
// }
resetVelocities();
Dt *= mSettings.xi; Dt *= mSettings.xi;
// if (mSettings.isDebugMode) { // if (mSettings.isDebugMode) {
// std::cout << Dt << std::endl; // std::cout << Dt << std::endl;
// } // }
++numOfDampings;
if (mSettings.shouldCreatePlots) { if (mSettings.shouldCreatePlots) {
history.redMarks.push_back(mCurrentSimulationStep); history.redMarks.push_back(mCurrentSimulationStep);
} }
@ -2452,11 +2492,12 @@ mesh->currentTotalPotentialEnergykN*/
results.executionTime = std::chrono::duration_cast<std::chrono::seconds>(endTime - beginTime) results.executionTime = std::chrono::duration_cast<std::chrono::seconds>(endTime - beginTime)
.count(); .count();
if (mCurrentSimulationStep == settings.maxDRMIterations && mCurrentSimulationStep != 0) { if (settings.maxDRMIterations.has_value() && mCurrentSimulationStep == settings.maxDRMIterations
&& mCurrentSimulationStep != 0) {
if (mSettings.beVerbose) { if (mSettings.beVerbose) {
std::cout << "Did not reach equilibrium before reaching the maximum number " std::cout << "Did not reach equilibrium before reaching the maximum number "
"of DRM steps (" "of DRM steps ("
<< settings.maxDRMIterations << "). Breaking simulation" << std::endl; << settings.maxDRMIterations.value() << "). Breaking simulation" << std::endl;
} }
results.converged = false; results.converged = false;
} else if (std::isnan(pMesh->currentTotalKineticEnergy)) { } else if (std::isnan(pMesh->currentTotalKineticEnergy)) {
@ -2467,7 +2508,7 @@ mesh->currentTotalPotentialEnergykN*/
} }
// mesh.printVertexCoordinates(mesh.VN() / 2); // mesh.printVertexCoordinates(mesh.VN() / 2);
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw && !mSettings.isDebugMode) { if (mSettings.shouldDraw && !mSettings.debugModeStep.has_value()) {
draw(); draw();
} }
#endif #endif
@ -2568,13 +2609,13 @@ mesh->currentTotalPotentialEnergykN*/
} }
} }
if (!mSettings.isDebugMode && mSettings.shouldCreatePlots) { if (!mSettings.debugModeStep.has_value() && mSettings.shouldCreatePlots) {
SimulationResultsReporter reporter; SimulationResultsReporter reporter;
reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); reporter.reportResults({results}, "Results", pJob->pMesh->getLabel());
} }
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw || mSettings.isDebugMode) { if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) {
polyscope::removeCurveNetwork(meshPolyscopeLabel); polyscope::removeCurveNetwork(meshPolyscopeLabel);
polyscope::removeCurveNetwork("Initial_" + meshPolyscopeLabel); polyscope::removeCurveNetwork("Initial_" + meshPolyscopeLabel);
} }
@ -2583,18 +2624,103 @@ mesh->currentTotalPotentialEnergykN*/
return results; return results;
} }
bool DRMSimulationModel::Settings::save(const filesystem::__cxx11::path &folderPath) const void DRMSimulationModel::Settings::save(const std::filesystem::path &folderPath) const
{ {
bool returnValue = true;
std::filesystem::create_directories(folderPath); std::filesystem::create_directories(folderPath);
nlohmann::json json; nlohmann::json json;
json[jsonLabels.meshFilename]
= std::filesystem::relative(std::filesystem::path(meshFilename), json[jsonLabels.gamma] = gamma;
std::filesystem::path(jsonFilename).parent_path()) json[jsonLabels.shouldDraw] = shouldDraw ? "true" : "false";
.string(); json[jsonLabels.shouldCreatePlots] = shouldCreatePlots ? "true" : "false";
json[jsonLabels.beVerbose] = beVerbose ? "true" : "false";
json[jsonLabels.Dtini] = Dtini;
json[jsonLabels.xi] = xi;
if (maxDRMIterations.has_value()) {
json[jsonLabels.maxDRMIterations] = maxDRMIterations.value();
}
if (debugModeStep.has_value()) {
json[jsonLabels.debugModeStep] = debugModeStep.value();
}
if (totalTranslationalKineticEnergyThreshold.has_value()) {
json[jsonLabels.totalTranslationalKineticEnergyThreshold]
= totalTranslationalKineticEnergyThreshold.value();
}
if (averageResidualForcesCriterionThreshold.has_value()) {
json[jsonLabels.averageResidualForcesCriterionThreshold]
= averageResidualForcesCriterionThreshold.value();
}
if (averageResidualForcesCriterionThreshold.has_value()) {
json[jsonLabels.averageResidualForcesCriterionThreshold]
= averageResidualForcesCriterionThreshold.value();
}
if (linearGuessForceScaleFactor.has_value()) {
json[jsonLabels.linearGuessForceScaleFactor] = linearGuessForceScaleFactor.value();
}
const std::string jsonFilename = "drmSettings.json"; const std::string jsonFilename = "drmSettings.json";
std::ofstream jsonFile(std::filesystem::path(folderPath).append(jsonFilename)); std::ofstream jsonFile(std::filesystem::path(folderPath).append(jsonFilename));
jsonFile << json; jsonFile << json;
} }
bool DRMSimulationModel::Settings::load(const filesystem::__cxx11::path &filePath) const {} bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePath)
{
if (!std::filesystem::exists(std::filesystem::path(jsonFilePath))) {
std::cerr << "The json file does not exist. Json file provided:" << jsonFilePath.string()
<< std::endl;
assert(false);
return false;
}
if (std::filesystem::path(jsonFilePath).extension() != ".json") {
std::cerr << "A json file is expected as input. The given file has the "
"following extension:"
<< std::filesystem::path(jsonFilePath).extension() << std::endl;
assert(false);
return false;
}
nlohmann::json json;
std::ifstream ifs(jsonFilePath.string());
ifs >> json;
if (json.contains(jsonLabels.shouldDraw)) {
shouldDraw = json.at(jsonLabels.shouldDraw) == "true" ? true : false;
}
if (json.contains(jsonLabels.beVerbose)) {
beVerbose = json.at(jsonLabels.beVerbose) == "true" ? true : false;
}
if (json.contains(jsonLabels.shouldCreatePlots)) {
shouldCreatePlots = json.at(jsonLabels.shouldCreatePlots) == "true" ? true : false;
}
if (json.contains(jsonLabels.Dtini)) {
Dtini = json.at(jsonLabels.Dtini);
}
if (json.contains(jsonLabels.xi)) {
xi = json.at(jsonLabels.xi);
}
if (json.contains(jsonLabels.maxDRMIterations)) {
maxDRMIterations = json.at(jsonLabels.maxDRMIterations);
}
if (json.contains(jsonLabels.debugModeStep)) {
debugModeStep = json.at(jsonLabels.debugModeStep);
}
if (json.contains(jsonLabels.totalTranslationalKineticEnergyThreshold)) {
totalTranslationalKineticEnergyThreshold = json.at(
jsonLabels.totalTranslationalKineticEnergyThreshold);
}
if (json.contains(jsonLabels.gamma)) {
gamma = json.at(jsonLabels.gamma);
}
if (json.contains(jsonLabels.averageResidualForcesCriterionThreshold)) {
averageResidualForcesCriterionThreshold = json.at(
jsonLabels.averageResidualForcesCriterionThreshold);
}
if (json.contains(jsonLabels.linearGuessForceScaleFactor)) {
linearGuessForceScaleFactor = json.at(jsonLabels.linearGuessForceScaleFactor);
}
return true;
}

44
drmsimulationmodel.hpp
View File

@ -25,32 +25,40 @@ class DRMSimulationModel
public: public:
struct Settings struct Settings
{ {
// bool isDebugMode{false};
std::optional<int> debugModeStep{0};
bool shouldDraw{false}; bool shouldDraw{false};
bool beVerbose{false}; bool beVerbose{false};
bool shouldCreatePlots{false}; bool shouldCreatePlots{false};
// int drawingStep{0};
// double residualForcesMovingAverageDerivativeNormThreshold{1e-8};
// double residualForcesMovingAverageNormThreshold{1e-8};
double Dtini{0.1}; double Dtini{0.1};
double xi{0.9969}; double xi{0.9969};
int maxDRMIterations{0}; int gradualForcedDisplacementSteps{50};
std::optional<double> shouldUseTranslationalKineticEnergyThreshold;
// int gradualForcedDisplacementSteps{50};
// int desiredGradualExternalLoadsSteps{1};
double gamma{0.8}; double gamma{0.8};
std::optional<double> displacementCap;
double totalResidualForcesNormThreshold{1e-3}; double totalResidualForcesNormThreshold{1e-3};
double totalExternalForcesNormPercentageTermination{1e-3}; double totalExternalForcesNormPercentageTermination{1e-3};
std::optional<double> averageResidualForcesCriterionThreshold{1e-5}; std::optional<int> maxDRMIterations;
std::optional<int> debugModeStep;
std::optional<double> totalTranslationalKineticEnergyThreshold;
std::optional<double> averageResidualForcesCriterionThreshold;
std::optional<double> linearGuessForceScaleFactor;
Settings() {} Settings() {}
bool save(const std::filesystem::path &folderPath) const; void save(const std::filesystem::path &folderPath) const;
bool load(const std::filesystem::path &filePath) const; bool load(const std::filesystem::path &filePath);
bool useViscousDamping{false};
bool useKineticDamping{true};
struct JsonLabels struct JsonLabels
{}; {
const std::string shouldDraw{"shouldDraw"};
const std::string beVerbose{"beVerbose"};
const std::string shouldCreatePlots{"shouldCreatePlots"};
const std::string Dtini{"DtIni"};
const std::string xi{"xi"};
const std::string maxDRMIterations{"maxIterations"};
const std::string debugModeStep{"debugModeStep"};
const std::string totalTranslationalKineticEnergyThreshold{
"totalTranslationaKineticEnergyThreshold"};
const std::string gamma{"gamma"};
const std::string averageResidualForcesCriterionThreshold{
"averageResidualForcesThreshold"};
const std::string linearGuessForceScaleFactor{"linearGuessForceScaleFactor"};
};
static JsonLabels jsonLabels;
}; };
private: private:
@ -219,6 +227,8 @@ private:
void updateNodeNr(VertexType &v); void updateNodeNr(VertexType &v);
void applyKineticDamping(const std::shared_ptr<SimulationJob> &pJob);
public: public:
DRMSimulationModel(); DRMSimulationModel();
SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob, SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
@ -228,6 +238,8 @@ private:
static void runUnitTests(); static void runUnitTests();
}; };
inline DRMSimulationModel::Settings::JsonLabels DRMSimulationModel::Settings::jsonLabels;
template <typename PointType> PointType Cross(PointType p1, PointType p2) { template <typename PointType> PointType Cross(PointType p1, PointType p2) {
return p1 ^ p2; return p1 ^ p2;
} }

2
polymesh.hpp
View File

@ -192,10 +192,10 @@ public:
// getVertices(); // getVertices();
// } // }
vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this); vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::EdgeEdge(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this); vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this); vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::AllocateEdge(*this); vcg::tri::UpdateTopology<VCGPolyMesh>::AllocateEdge(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::EdgeEdge(*this);
// vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this); // vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
return true; return true;

8
topologyenumerator.cpp
View File

@ -170,7 +170,7 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
// std::filesystem::path(resultsPath).append("patterns.patt")); // std::filesystem::path(resultsPath).append("patterns.patt"));
// } // }
if (numberOfDesiredEdges == -1) { if (numberOfDesiredEdges == -1) {
for (size_t numberOfEdges = 2; numberOfEdges < validEdges.size(); numberOfEdges++) { for (size_t numberOfEdges = 2; numberOfEdges <= validEdges.size(); numberOfEdges++) {
std::cout << "Computing " + setupString << " with " << numberOfEdges << " edges." std::cout << "Computing " + setupString << " with " << numberOfEdges << " edges."
<< std::endl; << std::endl;
@ -478,6 +478,7 @@ void TopologyEnumerator::computeValidPatterns(
const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold
= patternGeometry.hasAngleSmallerThanThreshold(numberOfNodesPerSlot, 15); = patternGeometry.hasAngleSmallerThanThreshold(numberOfNodesPerSlot, 15);
if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) { if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) {
statistics.numberOfPatternsViolatingAngleThreshold++;
if (debugIsOn /*|| savePlyFiles*/) { if (debugIsOn /*|| savePlyFiles*/) {
if (savePlyFiles) { if (savePlyFiles) {
exportPattern(std::filesystem::path(resultsPath) exportPattern(std::filesystem::path(resultsPath)
@ -493,6 +494,7 @@ void TopologyEnumerator::computeValidPatterns(
const bool tiledPatternHasNodeWithValenceGreaterThanDesired const bool tiledPatternHasNodeWithValenceGreaterThanDesired
= patternGeometry.hasValenceGreaterThan(numberOfNodesPerSlot, 6); = patternGeometry.hasValenceGreaterThan(numberOfNodesPerSlot, 6);
if (tiledPatternHasNodeWithValenceGreaterThanDesired) { if (tiledPatternHasNodeWithValenceGreaterThanDesired) {
statistics.numberOfPatternsViolatingValenceThreshold++;
if (debugIsOn) { if (debugIsOn) {
if (savePlyFiles) { if (savePlyFiles) {
auto highValencePath = std::filesystem::path(resultsPath) auto highValencePath = std::filesystem::path(resultsPath)
@ -648,6 +650,9 @@ void TopologyEnumerator::computeValidPatterns(
PatternIO::save(compressedPatternsFilePath, patternSet); PatternIO::save(compressedPatternsFilePath, patternSet);
} }
std::cout << "patternIndex:" << patternIndex << std::endl;
std::cout << "numberOfPatterns:" << numberOfPatterns << std::endl;
if (!validPatternsExist) { if (!validPatternsExist) {
std::filesystem::remove_all(validPatternsPath); std::filesystem::remove_all(validPatternsPath);
if (!debugIsOn) { if (!debugIsOn) {
@ -655,6 +660,7 @@ void TopologyEnumerator::computeValidPatterns(
} }
} }
} }
std::vector<size_t> TopologyEnumerator::getCoincideEdges( std::vector<size_t> TopologyEnumerator::getCoincideEdges(
const std::vector<size_t> &edgeNodeIndices) const const std::vector<size_t> &edgeNodeIndices) const
{ {

6
topologyenumerator.hpp
View File

@ -62,6 +62,8 @@ class TopologyEnumerator {
size_t numberOfPatternsWithMoreThanASingleCC{0}; size_t numberOfPatternsWithMoreThanASingleCC{0};
size_t numberOfPatternsWithADanglingEdgeOrNode{0}; size_t numberOfPatternsWithADanglingEdgeOrNode{0};
size_t numberOfPatternsWithArticulationPoints{0}; size_t numberOfPatternsWithArticulationPoints{0};
size_t numberOfPatternsViolatingAngleThreshold{0};
size_t numberOfPatternsViolatingValenceThreshold{0};
size_t numberOfValidPatterns{0}; size_t numberOfValidPatterns{0};
nlohmann::json convertToJson() const nlohmann::json convertToJson() const
{ {
@ -79,6 +81,8 @@ class TopologyEnumerator {
json["numPatternsWithDangling"] = numberOfPatternsWithADanglingEdgeOrNode; json["numPatternsWithDangling"] = numberOfPatternsWithADanglingEdgeOrNode;
json["numPatternsWithArticulationPoints"] = numberOfPatternsWithArticulationPoints; json["numPatternsWithArticulationPoints"] = numberOfPatternsWithArticulationPoints;
json["numValidPatterns"] = numberOfValidPatterns; json["numValidPatterns"] = numberOfValidPatterns;
json["violatingAngle"] = numberOfPatternsViolatingAngleThreshold;
json["violatingValence"] = numberOfPatternsViolatingValenceThreshold;
return json; return json;
} }
@ -141,6 +145,8 @@ class TopologyEnumerator {
numberOfPatternsWithADanglingEdgeOrNode = 0; numberOfPatternsWithADanglingEdgeOrNode = 0;
numberOfPatternsWithArticulationPoints = 0; numberOfPatternsWithArticulationPoints = 0;
numberOfValidPatterns = 0; numberOfValidPatterns = 0;
numberOfPatternsViolatingAngleThreshold = 0;
numberOfPatternsViolatingValenceThreshold = 0;
} }
}; };