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

CG first Submission state

This commit is contained in:
iasonmanolas 2021-10-14 14:58:34 +03:00
parent f949c2d793
commit b58ae0e13b
3 changed files with 129 additions and 89 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/main.cpp
View File

@ -53,6 +53,7 @@ int main(int argc, char *argv[]) {
input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100;
settings_optimization.normalizationStrategy
= ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
settings_optimization.normalizationParameter = 0.0003;
settings_optimization.solverAccuracy = 0.001;
settings_optimization.objectiveWeights.translational = std::atof(argv[4]);
@ -117,6 +118,7 @@ int main(int argc, char *argv[]) {
if (optimizationResults.wasSuccessful) {
resultsOutputDir = convergedJobsDirPath.string();
csvFile csv_results({}, false);
// csvFile csv_results(std::filesystem::path(resultsOutputDir).append("optimizationDistances.csv"), false);
csv_results << "Name";
optimizationResults.writeHeaderTo(csv_results);
settings_optimization.writeHeaderTo(csv_results);
@ -131,10 +133,10 @@ int main(int argc, char *argv[]) {
optimizationResults.save(resultsOutputDir, true);
}
//#ifdef POLYSCOPE_DEFINED
// optimizationResults.saveMeshFiles();
// optimizationResults.draw();
//#endif
#ifdef POLYSCOPE_DEFINED
// optimizationResults.saveMeshFiles();
optimizationResults.draw();
#endif
return 0;
}

205
src/reducedmodeloptimizer.cpp
View File

@ -39,6 +39,7 @@ struct GlobalOptimizationVariables {
FullPatternVertexIndex interfaceViForComputingScenarioError;
double desiredMaxDisplacementValue;
double desiredMaxRotationAngle;
std::string currentScenarioName;
} global;
double ReducedModelOptimizer::computeDisplacementError(
@ -190,15 +191,17 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
}
//#ifdef POLYSCOPE_DEFINED
// std::cout << "sim error:" << simulationScenarioError << std::endl;
// global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->registerForDrawing(
// ReducedModelOptimization::Colors::reducedInitial);
// reducedModelResults.registerForDrawing(
// ReducedModelOptimization::Colors::reducedDeformed);
// polyscope::show();
// reducedModelResults.unregister();
//#endif
#ifdef POLYSCOPE_DEFINED
// std::cout << reducedJob->getLabel() << " sim error:" << simulationScenarioError
// << std::endl;
// reducedJob->pMesh->registerForDrawing(Colors::reducedInitial);
// reducedModelResults.registerForDrawing(Colors::reducedDeformed);
// global.fullPatternResults[simulationScenarioIndex].registerForDrawing(
// Colors::fullDeformed);
// polyscope::show();
// reducedModelResults.unregister();
// global.fullPatternResults[simulationScenarioIndex].unregister();
#endif
// if (global.optimizationSettings.normalizationStrategy !=
// NormalizationStrategy::Epsilon &&
// simulationScenarioError > 1) {
@ -214,20 +217,20 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
// }
// }
//#ifdef POLYSCOPE_DEFINED
// ReducedModelOptimizer::visualizeResults(
// global.fullPatternSimulationJobs[simulationScenarioIndex],
// global.reducedPatternSimulationJobs[simulationScenarioIndex],
// global.reducedToFullInterfaceViMap, false);
// #ifdef POLYSCOPE_DEFINED
// ReducedModelOptimizer::visualizeResults(
// global.fullPatternSimulationJobs[simulationScenarioIndex],
// global.reducedPatternSimulationJobs[simulationScenarioIndex],
// global.reducedToFullInterfaceViMap, false);
// ReducedModelOptimizer::visualizeResults(
// global.fullPatternSimulationJobs[simulationScenarioIndex],
// std::make_shared<SimulationJob>(
// reducedPatternMaximumDisplacementSimulationJobs
// [simulationScenarioIndex]),
// global.reducedToFullInterfaceViMap, true);
// polyscope::removeAllStructures();
//#endif // POLYSCOPE_DEFINED
// ReducedModelOptimizer::visualizeResults(
// global.fullPatternSimulationJobs[simulationScenarioIndex],
// std::make_shared<SimulationJob>(
// reducedPatternMaximumDisplacementSimulationJobs
// [simulationScenarioIndex]),
// global.reducedToFullInterfaceViMap, true);
// polyscope::removeAllStructures();
// #endif // POLYSCOPE_DEFINED
totalError += simulationScenarioError;
}
}
@ -272,14 +275,14 @@ void ReducedModelOptimizer::createSimulationMeshes(
pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
pFullPatternSimulationMesh->setBeamCrossSection(
CrossSectionType{0.002, 0.002});
pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
pFullPatternSimulationMesh->setBeamMaterial(0.3, 2.3465 * 1e9);
// Reduced pattern
pReducedPatternSimulationMesh =
std::make_shared<SimulationMesh>(reducedModel);
pReducedPatternSimulationMesh->setBeamCrossSection(
CrossSectionType{0.002, 0.002});
pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
pReducedPatternSimulationMesh->setBeamMaterial(0.3, 2.3465 * 1e9);
}
void ReducedModelOptimizer::createSimulationMeshes(
@ -474,7 +477,7 @@ void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern,
computeMaps(copyFullPattern, copyReducedPattern);
createSimulationMeshes(copyFullPattern, copyReducedPattern);
initializeOptimizationParameters(m_pReducedPatternSimulationMesh,optimizationParameters);
initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters);
}
void updateMesh(long n, const double *x) {
@ -488,7 +491,7 @@ void updateMesh(long n, const double *x) {
const double beamHeight=beamWidth;
const double J=global.initialParameters(2) * x[2];
const double I2=global.initialParameters(3) * x[3];
const double I2 = global.initialParameters(3) * x[3];
const double I3=global.initialParameters(4) * x[4];
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
@ -594,29 +597,26 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob(
// &reducedToFullInterfaceViMap)
//{
// DRMSimulationModel simulator;
// std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
// fullPatternSimulationJobs[0]->pMesh;
// std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh = fullPatternSimulationJobs[0]->pMesh;
// pFullPatternSimulationMesh->registerForDrawing();
// pFullPatternSimulationMesh->save(pFullPatternSimulationMesh->getLabel() + "_undeformed.ply");
// reducedPatternSimulationJobs[0]->pMesh->save(reducedPatternSimulationJobs[0]->pMesh->getLabel()
// + "_undeformed.ply");
// double totalError = 0;
// for (const int simulationScenarioIndex : simulationScenarios) {
// const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
// fullPatternSimulationJobs[simulationScenarioIndex];
// pFullPatternSimulationJob->registerForDrawing(
// pFullPatternSimulationMesh->getLabel());
// SimulationResults fullModelResults =
// simulator.executeSimulation(pFullPatternSimulationJob);
// const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
// = fullPatternSimulationJobs[simulationScenarioIndex];
// pFullPatternSimulationJob->registerForDrawing(pFullPatternSimulationMesh->getLabel());
// SimulationResults fullModelResults = simulator.executeSimulation(pFullPatternSimulationJob);
// fullModelResults.registerForDrawing();
// // fullModelResults.saveDeformedModel();
// const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
// reducedPatternSimulationJobs[simulationScenarioIndex];
// SimulationResults reducedModelResults =
// simulator.executeSimulation(pReducedPatternSimulationJob);
// const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob
// = reducedPatternSimulationJobs[simulationScenarioIndex];
// SimulationResults reducedModelResults = simulator.executeSimulation(
// pReducedPatternSimulationJob);
// double normalizationFactor = 1;
// if (global.optimizationSettings.normalizationStrategy !=
// ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) {
// if (global.optimizationSettings.normalizationStrategy
// != ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) {
// normalizationFactor
// = global.translationalDisplacementNormalizationValues[simulationScenarioIndex];
// }
@ -627,8 +627,7 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob(
// reducedToFullInterfaceViMap,
// normalizationFactor);
// std::cout << "Error of simulation scenario "
// <baseSimulationScenarioNames[simulationScenarioIndex]
// << " is " << error << std::endl;
// < baseSimulationScenarioNames[simulationScenarioIndex] << " is " << error << std::endl;
// totalError += error;
// reducedModelResults.registerForDrawing();
// // firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing();
@ -1075,22 +1074,23 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni
DRMSimulationModel simulator;
DRMSimulationModel::Settings settings;
settings.totalExternalForcesNormPercentageTermination = 1e-2;
settings.totalTranslationalKineticEnergyThreshold = 1e-9;
settings.shouldUseTranslationalKineticEnergyThreshold = true;
// settings.shouldDraw = true;
// settings.useAverage = true;
// settings.isDebugMode = true;
// settings.drawingStep = 500;
// settings.beVerbose = true;
// settings.debugModeStep = 200000;
// settings.shouldCreatePlots = true;
settings.maxDRMIterations = 100000;
settings.totalTranslationalKineticEnergyThreshold = 1e-8;
settings.viscousDampingFactor = 1e-2;
settings.useKineticDamping = true;
settings.shouldDraw = false;
settings.debugModeStep = 200;
// settings.averageResidualForcesCriterionThreshold = 1e-5;
settings.maxDRMIterations = 200000;
SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
settings);
const double &desiredRotationAngle = global.desiredMaxRotationAngle;
double error;
#ifdef POLYSCOPE_DEFINED
job.pMesh->setLabel("initial");
// job.pMesh->registerForDrawing();
// results.registerForDrawing();
// polyscope::show();
std::string saveJobToPath;
if (!results.converged) {
// std::cout << "Force used:" << forceMagnitude << std::endl;
error = std::numeric_limits<double>::max();
@ -1108,12 +1108,21 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni
// debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
// auto debugResults = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
// debugSimulationSettings);
// std::terminate();
saveJobToPath = "../nonConvergingJobs";
} else {
error = std::abs(
results.rotationalDisplacementQuaternion[global.interfaceViForComputingScenarioError]
.angularDistance(Eigen::Quaterniond::Identity())
- desiredRotationAngle);
saveJobToPath = "../convergingJobs";
}
std::filesystem::path outputPath(std::filesystem::path(saveJobToPath)
.append(job.pMesh->getLabel())
.append("mag_" + global.currentScenarioName));
std::filesystem::create_directories(outputPath);
job.save(outputPath);
settings.save(outputPath);
std::cout << "Force:" << forceMagnitude << " Error is:" << vcg::math::ToDeg(error) << std::endl;
#endif
@ -1133,23 +1142,32 @@ double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMa
DRMSimulationModel::Settings settings;
// settings.totalResidualForcesNormThreshold = 1e-3;
settings.totalExternalForcesNormPercentageTermination = 1e-2;
settings.totalTranslationalKineticEnergyThreshold = 1e-9;
settings.shouldUseTranslationalKineticEnergyThreshold = true;
// settings.useAverage = true;
settings.totalTranslationalKineticEnergyThreshold = 1e-8;
settings.viscousDampingFactor = 5e-3;
settings.useKineticDamping = true;
// settings.averageResidualForcesCriterionThreshold = 1e-5;
// 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 = 100000;
// settings.debugModeStep = 200;
settings.maxDRMIterations = 200000;
SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
settings);
const double &desiredDisplacementValue = global.desiredMaxDisplacementValue;
double error;
if (!results.converged) {
error = std::numeric_limits<double>::max();
std::filesystem::path outputPath(std::filesystem::path("../nonConvergingJobs")
.append(job.pMesh->getLabel())
.append("mag_" + global.currentScenarioName));
std::filesystem::create_directories(outputPath);
job.save(outputPath);
settings.save(outputPath);
// std::terminate();
} else {
error = std::abs(
results.displacements[global.interfaceViForComputingScenarioError].getTranslation().norm()
@ -1180,6 +1198,7 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce(
objectiveFunction = &fullPatternMaxSimulationForceTranslationalObjective;
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
global.desiredMaxDisplacementValue = 0.1 * baseTriangleHeight;
global.currentScenarioName = baseSimulationScenarioNames[scenario];
double forceMagnitudeEpsilon = 1e-4;
switch (scenario) {
@ -1187,16 +1206,16 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce(
global.desiredMaxDisplacementValue = 0.03 * baseTriangleHeight;
break;
case Shear:
global.desiredMaxDisplacementValue = 0.04 * baseTriangleHeight;
global.desiredMaxDisplacementValue = 0.03 * baseTriangleHeight;
break;
case Bending:
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first;
break;
case Dome:
global.desiredMaxRotationAngle = vcg::math::ToRad(35.0);
global.desiredMaxRotationAngle = vcg::math::ToRad(25.0);
objectiveFunction = &fullPatternMaxSimulationForceRotationalObjective;
forceMagnitudeEpsilon *= 1e-2;
objectiveEpsilon = vcg::math::ToRad(3.0);
forceMagnitude = 0.6;
break;
case Saddle:
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first;
@ -1224,9 +1243,11 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce(
std::filesystem::path outputPath(
std::filesystem::path("../nonConvergingJobs")
.append(m_pFullPatternSimulationMesh->getLabel())
.append("mag_" + ReducedPatternOptimization::baseSimulationScenarioNames[scenario]));
.append("magFinal_"
+ ReducedPatternOptimization::baseSimulationScenarioNames[scenario]));
std::filesystem::create_directories(outputPath);
job.save(outputPath);
std::terminate();
}
#endif
@ -1247,15 +1268,14 @@ std::vector<std::shared_ptr<SimulationJob>> ReducedModelOptimizer::createFullPat
scenarioMaxForceMagnitudePairs) {
const BaseSimulationScenario scenario = scenarioMaxForceMagnitudePair.first;
const double maxForceMagnitude = scenarioMaxForceMagnitudePair.second;
const double minForceMagnitude = scenarioIsSymmetrical[scenario] ? 0 : -maxForceMagnitude;
const int numberOfSimulationScenarios = simulationScenariosResolution[scenario];
const int forceMagnitudeSamples = scenarioIsSymmetrical[scenario]
? numberOfSimulationScenarios - 1
: numberOfSimulationScenarios;
const double minForceMagnitude = scenarioIsSymmetrical[scenario]
? maxForceMagnitude / numberOfSimulationScenarios
: -maxForceMagnitude;
const double forceMagnitudeStep = numberOfSimulationScenarios == 1
? maxForceMagnitude
: (maxForceMagnitude - minForceMagnitude)
/ (forceMagnitudeSamples);
/ (numberOfSimulationScenarios);
const int baseSimulationScenarioIndexOffset
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin() + scenario,
@ -1344,11 +1364,10 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
global.translationalDisplacementNormalizationValues[simulationScenarioIndex]
= std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex],
epsilon_translationalDisplacement);
// const double epsilon_rotationalDisplacement = vcg::math::ToRad(10.0);
const double epsilon_rotationalDisplacement = vcg::math::ToRad(3.0);
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
= /*std::max(*/ fullPatternAngularDistance[simulationScenarioIndex] /*,
epsilon_rotationalDisplacement)*/
;
= std::max(fullPatternAngularDistance[simulationScenarioIndex],
epsilon_rotationalDisplacement);
} else {
global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = 1;
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = 1;
@ -1396,9 +1415,12 @@ void ReducedModelOptimizer::optimize(
results.baseTriangle = global.baseTriangle;
DRMSimulationModel::Settings simulationSettings;
simulationSettings.maxDRMIterations = 100000;
simulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
simulationSettings.maxDRMIterations = 200000;
simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8;
simulationSettings.viscousDampingFactor = 5e-3;
simulationSettings.useKineticDamping = true;
// simulationSettings.averageResidualForcesCriterionThreshold = 1e-5;
// simulationSettings.viscousDampingFactor = 1e-3;
// simulationSettings.beVerbose = true;
// simulationSettings.shouldDraw = true;
// simulationSettings.isDebugMode = true;
@ -1410,30 +1432,44 @@ void ReducedModelOptimizer::optimize(
ReducedPatternOptimization::Colors::fullInitial);
}
#endif
// LinearSimulationModel linearSimulator;
results.wasSuccessful = true;
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
= global.fullPatternSimulationJobs[simulationScenarioIndex];
SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob,
simulationSettings);
// if (!fullPatternResults.converged) {
// DRMSimulationModel::Settings simulationSettings_secondRound;
// simulationSettings_secondRound.viscousDampingFactor = 2e-3;
// simulationSettings_secondRound.useKineticDamping = true;
// simulationSettings.maxDRMIterations = 200000;
// fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob,
// simulationSettings_secondRound);
#ifdef POLYSCOPE_DEFINED
// std::cout << "Simulation job " << pFullPatternSimulationJob->getLabel()
// << " used viscous damping." << std::endl;
#endif
if (!fullPatternResults.converged) {
results.wasSuccessful = false;
#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.debugModeStep = 50;
// // // debugSimulationSettings.maxDRMIterations = 100000;
// debugSimulationSettings.shouldDraw = true;
// debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep;
// // debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep;
// debugSimulationSettings.shouldCreatePlots = true;
// // debugSimulationSettings.Dtini = 0.06;
// // // debugSimulationSettings.Dtini = 0.06;
// debugSimulationSettings.beVerbose = true;
// debugSimulationSettings.useAverage = true;
// // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3;
// // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
// debugSimulationSettings.averageResidualForcesCriterionThreshold = 1e-5;
// debugSimulationSettings.maxDRMIterations = 100000;
// debugSimulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8;
// debugSimulationSettings.viscousDampingFactor = 1e-2;
// // // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3;
// // // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
// auto debugResults = simulator.executeSimulation(pFullPatternSimulationJob,
// debugSimulationSettings);
// debugResults.setLabelPrefix("debugResults");
@ -1446,16 +1482,17 @@ void ReducedModelOptimizer::optimize(
.append("final_" + pFullPatternSimulationJob->getLabel()));
std::filesystem::create_directories(outputPath);
pFullPatternSimulationJob->save(outputPath);
simulationSettings.save(outputPath);
#endif
std::terminate();
return;
// }
}
results.wasSuccessful = true;
#ifdef POLYSCOPE_DEFINED
if (drawFullPatternSimulationResults) {
// SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation(
// pFullPatternSimulationJob);
fullPatternResults.registerForDrawing(ReducedPatternOptimization::Colors::fullDeformed,
true,
true);
// fullPatternResults_linear.labelPrefix += "_linear";
// fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed,

3
src/reducedmodeloptimizer.hpp
View File

@ -38,7 +38,8 @@ class ReducedModelOptimizer
int fullPatternNumberOfEdges;
public:
constexpr static std::array<int, 5> simulationScenariosResolution = {10, 10, 20, 20, 20};
constexpr static std::array<int, 5> simulationScenariosResolution = {11, 11, 20, 20, 20};
// constexpr static std::array<int, 5> simulationScenariosResolution = {3, 3, 3, 3, 3};
inline static int totalNumberOfSimulationScenarios
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.end(),