Refactoring. Added kinetic energy as a convergence criterion for the full pattern simulations

This commit is contained in:
iasonmanolas 2021-07-20 16:30:41 +03:00
parent 19d9fe434a
commit f949c2d793
3 changed files with 100 additions and 135 deletions

View File

@ -114,21 +114,9 @@ int main(int argc, char *argv[]) {
auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
optimizationResults.time = elapsed.count() / 1000.0;
// Export results
if (!optimizationResults.wasSuccessful) {
resultsOutputDir = crashedJobsDirPath.string();
} else {
resultsOutputDir = convergedJobsDirPath.string();
}
optimizationResults.save(resultsOutputDir, true);
// Write results in csv
csvFile csv_results({}, false);
// csvFile csv_results(std::filesystem::path(dirPath_thisOptimization)
// .append("results.csv")
// .string(),
// false);
if (optimizationResults.wasSuccessful) {
resultsOutputDir = convergedJobsDirPath.string();
csvFile csv_results({}, false);
csv_results << "Name";
optimizationResults.writeHeaderTo(csv_results);
settings_optimization.writeHeaderTo(csv_results);
@ -137,7 +125,10 @@ int main(int argc, char *argv[]) {
optimizationResults.writeResultsTo(settings_optimization, csv_results);
settings_optimization.writeSettingsTo(csv_results);
csv_results << endrow;
} else {
resultsOutputDir = crashedJobsDirPath.string();
}
optimizationResults.save(resultsOutputDir, true);
}
//#ifdef POLYSCOPE_DEFINED

View File

@ -4,6 +4,7 @@
#include "trianglepatterngeometry.hpp"
#include "trianglepattterntopology.hpp"
#include <chrono>
#include <functional>
using namespace ReducedPatternOptimization;
@ -459,8 +460,8 @@ void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern,
PatternGeometry &reducedPattern,
const int &optimizationParameters)
{
assert(fullPattern.VN() == reducedPattern.VN() &&
fullPattern.EN() >= reducedPattern.EN());
assert(fullPattern.VN() == reducedPattern.VN() && fullPattern.EN() >= reducedPattern.EN());
fullPatternNumberOfEdges = fullPattern.EN();
#if POLYSCOPE_DEFINED
polyscope::removeAllStructures();
#endif
@ -873,6 +874,7 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces(
const std::vector<BaseSimulationScenario> &desiredBaseSimulationScenarioIndices)
{
std::vector<std::pair<BaseSimulationScenario, double>> fullPatternSimulationScenarioMaxMagnitudes;
#ifdef POLYSCOPE_DEFINED
const std::filesystem::path forceMagnitudesDirectoryPath(std::filesystem::current_path()
.parent_path()
.append("IntermediateResults")
@ -890,8 +892,11 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces(
= static_cast<std::vector<std::pair<BaseSimulationScenario, double>>>(
json.at("maxMagn"));
} else {
#endif
fullPatternSimulationScenarioMaxMagnitudes = computeFullPatternMaxSimulationForces(
desiredBaseSimulationScenarioIndices);
#ifdef POLYSCOPE_DEFINED
nlohmann::json json;
json["maxMagn"] = fullPatternSimulationScenarioMaxMagnitudes;
@ -899,6 +904,7 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces(
std::ofstream jsonFile(patternMaxForceMagnitudesFilePath.string());
jsonFile << json;
}
#endif
assert(fullPatternSimulationScenarioMaxMagnitudes.size()
== desiredBaseSimulationScenarioIndices.size());
@ -1069,10 +1075,10 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni
DRMSimulationModel simulator;
DRMSimulationModel::Settings settings;
settings.totalExternalForcesNormPercentageTermination = 1e-2;
settings.totalTranslationalKineticEnergyThreshold = 1e-11;
settings.totalTranslationalKineticEnergyThreshold = 1e-9;
settings.shouldUseTranslationalKineticEnergyThreshold = true;
// settings.shouldDraw = true;
settings.useAverage = true;
// settings.useAverage = true;
// settings.isDebugMode = true;
// settings.drawingStep = 500;
@ -1084,17 +1090,31 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni
settings);
const double &desiredRotationAngle = global.desiredMaxRotationAngle;
double error;
#ifdef POLYSCOPE_DEFINED
if (!results.converged) {
// std::cout << "Force used:" << forceMagnitude << std::endl;
error = std::numeric_limits<double>::max();
// DRMSimulationModel::Settings debugSimulationSettings;
// debugSimulationSettings.isDebugMode = true;
// debugSimulationSettings.debugModeStep = 1000;
// // debugSimulationSettings.maxDRMIterations = 100000;
// debugSimulationSettings.shouldDraw = true;
// debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep;
// debugSimulationSettings.shouldCreatePlots = true;
// // debugSimulationSettings.Dtini = 0.06;
// debugSimulationSettings.beVerbose = true;
// debugSimulationSettings.useAverage = true;
// // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3;
// debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
// auto debugResults = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
// debugSimulationSettings);
} else {
error = std::abs(
// results.displacements[global.fullPatternInterfaceViPairs[1].first].getTranslation().norm()
results.rotationalDisplacementQuaternion[global.interfaceViForComputingScenarioError]
.angularDistance(Eigen::Quaterniond::Identity())
- desiredRotationAngle);
}
#ifdef POLYSCOPE_DEFINED
std::cout << "Force:" << forceMagnitude << " Error is:" << vcg::math::ToDeg(error) << std::endl;
#endif
return error;
@ -1113,23 +1133,17 @@ double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMa
DRMSimulationModel::Settings settings;
// settings.totalResidualForcesNormThreshold = 1e-3;
settings.totalExternalForcesNormPercentageTermination = 1e-2;
settings.totalTranslationalKineticEnergyThreshold = 1e-10;
settings.totalTranslationalKineticEnergyThreshold = 1e-9;
settings.shouldUseTranslationalKineticEnergyThreshold = true;
// settings.totalTranslationalKineticEnergyThreshold = 1e-10;
// settings.useAverage = true;
// settings.totalTranslationalKineticEnergyThreshold = 1e-10;
// settings.shouldUseTranslationalKineticEnergyThreshold = true;
// settings.shouldDraw = true;
// settings.isDebugMode = true;
// settings.drawingStep = 200000;
// settings.beVerbose = true;
// settings.debugModeStep = 200000;
settings.maxDRMIterations = 250000;
#ifdef POLYSCOPE_DEFINED
// settings.shouldDraw = true;
// settings.isDebugMode = true;
// settings.drawingStep = 100000;
// settings.debugModeStep = 10000;
// settings.shouldCreatePlots = true;
#endif
settings.maxDRMIterations = 100000;
SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
settings);
const double &desiredDisplacementValue = global.desiredMaxDisplacementValue;
@ -1138,7 +1152,6 @@ double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMa
error = std::numeric_limits<double>::max();
} else {
error = std::abs(
// results.displacements[global.fullPatternInterfaceViPairs[1].first].getTranslation().norm()
results.displacements[global.interfaceViForComputingScenarioError].getTranslation().norm()
- desiredDisplacementValue);
}
@ -1154,114 +1167,68 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce(
const BaseSimulationScenario &scenario)
{
double forceMagnitude = 1;
const double forceMagnitudeEpsilon = 1e-2;
double minimumError;
double translationalOptimizationEpsilon;
dlib::function_evaluation result;
bool wasSuccessful = false;
global.constructScenarioFunction = constructBaseScenarioFunctions[scenario];
const double baseTriangleHeight = vcg::Distance(global.baseTriangle.cP(0),
(global.baseTriangle.cP(1)
+ global.baseTriangle.cP(2))
/ 2);
std::function<double(const double &)> objectiveFunction;
double objectiveEpsilon = translationalOptimizationEpsilon;
objectiveFunction = &fullPatternMaxSimulationForceTranslationalObjective;
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
global.desiredMaxDisplacementValue = 0.1 * baseTriangleHeight;
double forceMagnitudeEpsilon = 1e-4;
switch (scenario) {
case Axial:
global.desiredMaxDisplacementValue = 0.03
* vcg::Distance(global.baseTriangle.cP(0),
(global.baseTriangle.cP(1)
+ global.baseTriangle.cP(2))
/ 2);
global.constructScenarioFunction = &ReducedModelOptimizer::constructAxialSimulationScenario;
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
minimumError
= dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective,
forceMagnitude,
1e-2,
1e2,
forceMagnitudeEpsilon);
translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue;
wasSuccessful = minimumError < translationalOptimizationEpsilon;
global.desiredMaxDisplacementValue = 0.03 * baseTriangleHeight;
break;
case Shear:
global.desiredMaxDisplacementValue = 0.04
* vcg::Distance(global.baseTriangle.cP(0),
(global.baseTriangle.cP(1)
+ global.baseTriangle.cP(2))
/ 2);
global.constructScenarioFunction = &ReducedModelOptimizer::constructShearSimulationScenario;
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
minimumError
= dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective,
forceMagnitude,
1e-2,
1e2,
forceMagnitudeEpsilon);
translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue;
wasSuccessful = minimumError < translationalOptimizationEpsilon;
global.desiredMaxDisplacementValue = 0.04 * baseTriangleHeight;
break;
case Bending:
global.desiredMaxDisplacementValue = 0.1
* vcg::Distance(global.baseTriangle.cP(0),
(global.baseTriangle.cP(1)
+ global.baseTriangle.cP(2))
/ 2);
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first;
minimumError
= dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective,
forceMagnitude,
1e-2,
1e2,
forceMagnitudeEpsilon);
translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue;
wasSuccessful = minimumError < translationalOptimizationEpsilon;
break;
case Dome:
forceMagnitude = 0.005;
while (!wasSuccessful) {
global.desiredMaxRotationAngle = vcg::math::ToRad(35.0);
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
minimumError
= dlib::find_min_single_variable(&fullPatternMaxSimulationForceRotationalObjective,
forceMagnitude,
1e-3,
0.5,
forceMagnitudeEpsilon);
wasSuccessful = minimumError < vcg::math::ToRad(3.0);
}
// result = dlib::find_min_global(&fullPatternMaxSimulationForceRotationalObjective,
// 1e-2,
// 1,
// dlib::max_function_calls(50));
// wasSuccessful = result.y < vcg::math::ToRad(3.0);
global.desiredMaxRotationAngle = vcg::math::ToRad(35.0);
objectiveFunction = &fullPatternMaxSimulationForceRotationalObjective;
forceMagnitudeEpsilon *= 1e-2;
objectiveEpsilon = vcg::math::ToRad(3.0);
break;
case Saddle:
global.desiredMaxDisplacementValue = 0.1
* vcg::Distance(global.baseTriangle.cP(0),
(global.baseTriangle.cP(1)
+ global.baseTriangle.cP(2))
/ 2);
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first;
minimumError
= dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective,
forceMagnitude,
1e-2,
1e2,
forceMagnitudeEpsilon);
translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue;
wasSuccessful = minimumError < translationalOptimizationEpsilon;
break;
}
minimumError = dlib::find_min_single_variable(objectiveFunction,
forceMagnitude,
1e-4,
1e4,
forceMagnitudeEpsilon,
1000);
wasSuccessful = minimumError < objectiveEpsilon;
#ifdef POLYSCOPE_DEFINED
std::cout << "Max " << ReducedPatternOptimization::baseSimulationScenarioNames[scenario]
<< " magnitude:" << forceMagnitude << std::endl;
if (!wasSuccessful)
std::cout << "Was not successfull" << std::endl;
#endif
if (!wasSuccessful) {
const std::string debugMessage
= ReducedPatternOptimization::baseSimulationScenarioNames[scenario]
+ " max scenario magnitude was not succefully determined.";
optimizationNotes.append(debugMessage);
SimulationJob job;
job.pMesh = global.pFullPatternSimulationMesh;
global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job);
std::cout << ReducedPatternOptimization::baseSimulationScenarioNames[scenario]
+ " max scenario magnitude was not succefully determined."
<< std::endl;
std::filesystem::path outputPath(
std::filesystem::path("../nonConvergingJobs")
.append(m_pFullPatternSimulationMesh->getLabel())
.append("mag_" + ReducedPatternOptimization::baseSimulationScenarioNames[scenario]));
std::filesystem::create_directories(outputPath);
job.save(outputPath);
}
#endif
return forceMagnitude;
}
@ -1429,9 +1396,9 @@ void ReducedModelOptimizer::optimize(
results.baseTriangle = global.baseTriangle;
DRMSimulationModel::Settings simulationSettings;
simulationSettings.maxDRMIterations = 50000;
simulationSettings.maxDRMIterations = 100000;
simulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-10;
simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8;
// simulationSettings.beVerbose = true;
// simulationSettings.shouldDraw = true;
// simulationSettings.isDebugMode = true;
@ -1455,24 +1422,30 @@ void ReducedModelOptimizer::optimize(
#ifdef POLYSCOPE_DEFINED
std::cout << "Simulation job " << pFullPatternSimulationJob->getLabel()
<< " did not converge." << std::endl;
DRMSimulationModel::Settings debugSimulationSettings;
debugSimulationSettings.isDebugMode = true;
debugSimulationSettings.debugModeStep = 1000;
// debugSimulationSettings.maxDRMIterations = 100000;
debugSimulationSettings.shouldDraw = true;
debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep;
debugSimulationSettings.shouldCreatePlots = true;
// debugSimulationSettings.Dtini = 0.06;
debugSimulationSettings.beVerbose = true;
debugSimulationSettings.useAverage = true;
// debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3;
debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
auto debugResults = simulator.executeSimulation(pFullPatternSimulationJob,
debugSimulationSettings);
debugResults.setLabelPrefix("debugResults");
debugResults.registerForDrawing();
polyscope::show();
debugResults.unregister();
// DRMSimulationModel::Settings debugSimulationSettings;
// debugSimulationSettings.isDebugMode = true;
// debugSimulationSettings.debugModeStep = 1000;
// // debugSimulationSettings.maxDRMIterations = 100000;
// debugSimulationSettings.shouldDraw = true;
// debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep;
// debugSimulationSettings.shouldCreatePlots = true;
// // debugSimulationSettings.Dtini = 0.06;
// debugSimulationSettings.beVerbose = true;
// debugSimulationSettings.useAverage = true;
// // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3;
// // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
// auto debugResults = simulator.executeSimulation(pFullPatternSimulationJob,
// debugSimulationSettings);
// debugResults.setLabelPrefix("debugResults");
// debugResults.registerForDrawing();
// polyscope::show();
// debugResults.unregister();
std::filesystem::path outputPath(
std::filesystem::path("../nonConvergingJobs")
.append(m_pFullPatternSimulationMesh->getLabel())
.append("final_" + pFullPatternSimulationJob->getLabel()));
std::filesystem::create_directories(outputPath);
pFullPatternSimulationJob->save(outputPath);
#endif
return;
}

View File

@ -35,6 +35,7 @@ class ReducedModelOptimizer
SimulationJob &)>>
constructBaseScenarioFunctions;
std::vector<bool> scenarioIsSymmetrical;
int fullPatternNumberOfEdges;
public:
constexpr static std::array<int, 5> simulationScenariosResolution = {10, 10, 20, 20, 20};