From b58ae0e13bc3879fff1d399dfac5f72041a8b5eb Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Thu, 14 Oct 2021 14:58:34 +0300 Subject: [PATCH] CG first Submission state --- src/main.cpp | 10 +- src/reducedmodeloptimizer.cpp | 205 ++++++++++++++++++++-------------- src/reducedmodeloptimizer.hpp | 3 +- 3 files changed, 129 insertions(+), 89 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 837af09..5ae434d 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -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; } diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 1a9c4d3..8f10819 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -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( - // reducedPatternMaximumDisplacementSimulationJobs - // [simulationScenarioIndex]), - // global.reducedToFullInterfaceViMap, true); - // polyscope::removeAllStructures(); - //#endif // POLYSCOPE_DEFINED + // ReducedModelOptimizer::visualizeResults( + // global.fullPatternSimulationJobs[simulationScenarioIndex], + // std::make_shared( + // reducedPatternMaximumDisplacementSimulationJobs + // [simulationScenarioIndex]), + // global.reducedToFullInterfaceViMap, true); + // polyscope::removeAllStructures(); + // #endif // POLYSCOPE_DEFINED totalError += simulationScenarioError; } } @@ -272,14 +275,14 @@ void ReducedModelOptimizer::createSimulationMeshes( pFullPatternSimulationMesh = std::make_shared(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(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 pFullPatternSimulationMesh = -// fullPatternSimulationJobs[0]->pMesh; +// std::shared_ptr 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 &pFullPatternSimulationJob = -// fullPatternSimulationJobs[simulationScenarioIndex]; -// pFullPatternSimulationJob->registerForDrawing( -// pFullPatternSimulationMesh->getLabel()); -// SimulationResults fullModelResults = -// simulator.executeSimulation(pFullPatternSimulationJob); +// const std::shared_ptr &pFullPatternSimulationJob +// = fullPatternSimulationJobs[simulationScenarioIndex]; +// pFullPatternSimulationJob->registerForDrawing(pFullPatternSimulationMesh->getLabel()); +// SimulationResults fullModelResults = simulator.executeSimulation(pFullPatternSimulationJob); // fullModelResults.registerForDrawing(); // // fullModelResults.saveDeformedModel(); -// const std::shared_ptr &pReducedPatternSimulationJob = -// reducedPatternSimulationJobs[simulationScenarioIndex]; -// SimulationResults reducedModelResults = -// simulator.executeSimulation(pReducedPatternSimulationJob); +// const std::shared_ptr &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 " -// (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::max(); @@ -1108,12 +1108,21 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true; // auto debugResults = simulator.executeSimulation(std::make_shared(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(job), settings); const double &desiredDisplacementValue = global.desiredMaxDisplacementValue; double error; if (!results.converged) { error = std::numeric_limits::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> 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 &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, diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 54d54fd..a8280e3 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -38,7 +38,8 @@ class ReducedModelOptimizer int fullPatternNumberOfEdges; public: - constexpr static std::array simulationScenariosResolution = {10, 10, 20, 20, 20}; + constexpr static std::array simulationScenariosResolution = {11, 11, 20, 20, 20}; + // constexpr static std::array simulationScenariosResolution = {3, 3, 3, 3, 3}; inline static int totalNumberOfSimulationScenarios = std::accumulate(simulationScenariosResolution.begin(), simulationScenariosResolution.end(),