From 0863db658d43968688a2ab8d50fda9e80168c177 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Wed, 30 Jun 2021 10:24:54 +0300 Subject: [PATCH] Added max iterations on all drm simulations.using global optimizer for dome magnitude.Added bool variable for success in force magnitude search. --- src/reducedmodeloptimizer.cpp | 208 ++++++++++++++++++++-------------- src/reducedmodeloptimizer.hpp | 3 +- 2 files changed, 127 insertions(+), 84 deletions(-) diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 3d56299..8ee30e3 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -692,9 +692,8 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza const Settings &settings, ReducedPatternOptimization::Results &results) { - results.baseTriangle = global.baseTriangle; //Number of crashes - results.numberOfSimulationCrashes = global.numOfSimulationCrashes; + // results.numberOfSimulationCrashes = global.numOfSimulationCrashes; //Value of optimal objective Y results.objectiveValue.total = optimizationResult_dlib.y; //Optimal X values @@ -932,7 +931,7 @@ void ReducedModelOptimizer::constructDomeSimulationScenario( = Eigen::Vector3d(-interfaceVector[0], -interfaceVector[1], 0) - * 0.001 + * 0.01 * std::abs( forceMagnitude); //NOTE:Should the forced displacement change relatively to the magnitude? // * std::abs(forceMagnitude / maxForceMagnitude_dome); @@ -1000,21 +999,31 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni DRMSimulationModel simulator; DRMSimulationModel::Settings settings; - settings.totalExternalForcesNormPercentageTermination = 1e-2; - // settings.shouldUseTranslationalKineticEnergyThreshold = true; - // settings.totalTranslationalKineticEnergyThreshold = 1e-7; + // settings.totalResidualForcesNormThreshold = 1e-3; + // settings.totalTranslationalKineticEnergyThreshold = 1e-6; + // settings.shouldUseTranslationalKineticEnergyThreshold = true; + settings.shouldDraw = true; + // settings.isDebugMode = true; + // settings.drawingStep = 500; + // settings.beVerbose = true; + // settings.debugModeStep = 200000; + // settings.shouldCreatePlots = true; + settings.maxDRMIterations = 100000; SimulationResults results = simulator.executeSimulation(std::make_shared(job), settings); const double &desiredRotationAngle = global.desiredMaxRotationAngle; + if (!results.converged) { + return std::numeric_limits::max(); + } const double 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 +#ifdef POLYSCOPE_DEFINED + std::cout << "Force:" << forceMagnitude << " Error is:" << vcg::math::ToDeg(error) << std::endl; +#endif return error; } @@ -1029,36 +1038,47 @@ double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMa DRMSimulationModel simulator; DRMSimulationModel::Settings settings; - settings.totalExternalForcesNormPercentageTermination = 1e-2; - settings.shouldUseTranslationalKineticEnergyThreshold = true; - settings.totalTranslationalKineticEnergyThreshold = 1e-14; - settings.shouldUseTranslationalKineticEnergyThreshold = true; + settings.totalResidualForcesNormThreshold = 1e-3; + // settings.totalTranslationalKineticEnergyThreshold = 1e-10; + // settings.shouldUseTranslationalKineticEnergyThreshold = true; + // settings.shouldDraw = true; + // settings.isDebugMode = true; + // settings.drawingStep = 200000; + // settings.beVerbose = true; + // settings.debugModeStep = 200000; + settings.maxDRMIterations = 50000; #ifdef POLYSCOPE_DEFINED // settings.shouldDraw = true; // settings.isDebugMode = true; // settings.drawingStep = 100000; +// settings.debugModeStep = 10000; // settings.shouldCreatePlots = true; #endif - SimulationResults results = simulator.executeSimulation(std::make_shared(job), - settings); + SimulationResults results = simulator.executeSimulation(std::make_shared(job)); const double &desiredDisplacementValue = global.desiredMaxDisplacementValue; + if (!results.converged) { + return std::numeric_limits::max(); + } const double error = std::abs( // results.displacements[global.fullPatternInterfaceViPairs[1].first].getTranslation().norm() results.displacements[global.interfaceViForComputingScenarioError].getTranslation().norm() - desiredDisplacementValue); +#ifdef POLYSCOPE_DEFINED + std::cout << "Force:" << forceMagnitude << " Error is:" << error << std::endl; +#endif + return error; } -double ReducedModelOptimizer::getFullPatternMaxSimulationForce(const BaseSimulationScenario &scenario) +double ReducedModelOptimizer::getFullPatternMaxSimulationForce( + const BaseSimulationScenario &scenario, bool &wasSuccessful) { double forceMagnitude = 1; - global.desiredMaxDisplacementValue = 0.1 - * vcg::Distance(global.baseTriangle.cP(0), - (global.baseTriangle.cP(1) - + global.baseTriangle.cP(2)) - / 2); - const double optimizationEpsilon = 1e-1; + const double forceMagnitudeEpsilon = 1e-2; + double minimumError; + double translationalOptimizationEpsilon; + dlib::function_evaluation result; switch (scenario) { case Axial: global.desiredMaxDisplacementValue = 0.03 @@ -1068,11 +1088,14 @@ double ReducedModelOptimizer::getFullPatternMaxSimulationForce(const BaseSimulat / 2); global.constructScenarioFunction = &ReducedModelOptimizer::constructAxialSimulationScenario; global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first; - dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective, - forceMagnitude, - 1e-8, - 1e8, - optimizationEpsilon); + minimumError + = dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective, + forceMagnitude, + 1e-2, + 1e2, + forceMagnitudeEpsilon); + translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue; + wasSuccessful = minimumError < translationalOptimizationEpsilon; break; case Shear: global.desiredMaxDisplacementValue = 0.04 @@ -1082,11 +1105,14 @@ double ReducedModelOptimizer::getFullPatternMaxSimulationForce(const BaseSimulat / 2); global.constructScenarioFunction = &ReducedModelOptimizer::constructShearSimulationScenario; global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first; - dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective, - forceMagnitude, - 1e-8, - 1e8, - optimizationEpsilon); + minimumError + = dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective, + forceMagnitude, + 1e-2, + 1e2, + forceMagnitudeEpsilon); + translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue; + wasSuccessful = minimumError < translationalOptimizationEpsilon; break; case Bending: global.desiredMaxDisplacementValue = 0.1 @@ -1096,28 +1122,33 @@ double ReducedModelOptimizer::getFullPatternMaxSimulationForce(const BaseSimulat / 2); global.constructScenarioFunction = &ReducedModelOptimizer::constructBendingSimulationScenario; global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first; - dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective, - forceMagnitude, - 1e-8, - 1e8, - optimizationEpsilon); + minimumError + = dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective, + forceMagnitude, + 1e-2, + 1e2, + forceMagnitudeEpsilon); + translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue; + wasSuccessful = minimumError < translationalOptimizationEpsilon; break; case Dome: - global.desiredMaxRotationAngle = vcg::math::ToRad(35.0); + global.desiredMaxRotationAngle = vcg::math::ToRad(25.0); global.constructScenarioFunction = &ReducedModelOptimizer::constructDomeSimulationScenario; global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first; - dlib::find_min_single_variable( - &fullPatternMaxSimulationForceRotationalObjective, - forceMagnitude, - 1e-8, - 1e8, - vcg::math::ToRad(1.0) - // global.desiredMaxRotationAngle * 0.5, - // optimizationEpsilon, - ); + // minimumError + // = dlib::find_min_single_variable(&fullPatternMaxSimulationForceRotationalObjective, + // forceMagnitude, + // 1e-2, + // 1e2, + // 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); break; case Saddle: - // global.desiredMaxDisplacementValue *= 2; global.desiredMaxDisplacementValue = 0.1 * vcg::Distance(global.baseTriangle.cP(0), (global.baseTriangle.cP(1) @@ -1125,14 +1156,30 @@ double ReducedModelOptimizer::getFullPatternMaxSimulationForce(const BaseSimulat / 2); global.constructScenarioFunction = &ReducedModelOptimizer::constructSaddleSimulationScenario; global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first; - dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective, - forceMagnitude, - 1e-8, - 1e8, - 1e-2); + minimumError + = dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective, + forceMagnitude, + 1e-2, + 1e2, + forceMagnitudeEpsilon); + translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue; + wasSuccessful = minimumError < translationalOptimizationEpsilon; break; } +#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); + } + return forceMagnitude; } @@ -1148,11 +1195,9 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( job.pMesh = pMesh; //// Axial -#ifdef POLYSCOPE_DEFINED - std::cout << "Computing Axial scenarios.." << std::endl; -#endif - const double maxForceMagnitude_axial = getFullPatternMaxSimulationForce( - BaseSimulationScenario::Axial); + bool wasSuccessful; + const double maxForceMagnitude_axial + = getFullPatternMaxSimulationForce(BaseSimulationScenario::Axial, wasSuccessful); const double minForceMagnitude_axial = -maxForceMagnitude_axial; const int numberOfSimulationScenarios_axial = simulationScenariosResolution[BaseSimulationScenario::Axial]; @@ -1167,9 +1212,6 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( for (int axialSimulationScenarioIndex = 0; axialSimulationScenarioIndex < numberOfSimulationScenarios_axial; axialSimulationScenarioIndex++) { -#ifdef POLYSCOPE_DEFINED - std::cout << "Computing Axial scenario:" << axialSimulationScenarioIndex << std::endl; -#endif job.nodalExternalForces.clear(); job.constrainedVertices.clear(); job.nodalForcedDisplacements.clear(); @@ -1185,8 +1227,8 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( } //// Shear - const double maxForceMagnitude_shear = getFullPatternMaxSimulationForce( - BaseSimulationScenario::Shear); + const double maxForceMagnitude_shear + = getFullPatternMaxSimulationForce(BaseSimulationScenario::Shear, wasSuccessful); const double minForceMagnitude_shear = -maxForceMagnitude_shear; const int numberOfSimulationScenarios_shear = simulationScenariosResolution[BaseSimulationScenario::Shear]; @@ -1201,9 +1243,6 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( for (int shearSimulationScenarioIndex = 0; shearSimulationScenarioIndex < numberOfSimulationScenarios_shear; shearSimulationScenarioIndex++) { -#ifdef POLYSCOPE_DEFINED - std::cout << "Computing shear scenario:" << shearSimulationScenarioIndex << std::endl; -#endif job.constrainedVertices.clear(); job.nodalExternalForces.clear(); job.label = baseSimulationScenarioNames[BaseSimulationScenario::Shear] + "_" @@ -1216,8 +1255,8 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( } //// Bending - const double maxForceMagnitude_bending = getFullPatternMaxSimulationForce( - BaseSimulationScenario::Bending); + const double maxForceMagnitude_bending + = getFullPatternMaxSimulationForce(BaseSimulationScenario::Bending, wasSuccessful); const double minForceMagnitude_bending = 0; const int numberOfSimulationScenarios_bending = simulationScenariosResolution[BaseSimulationScenario::Bending]; @@ -1230,9 +1269,6 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( for (int bendingSimulationScenarioIndex = 0; bendingSimulationScenarioIndex < numberOfSimulationScenarios_bending; bendingSimulationScenarioIndex++) { -#ifdef POLYSCOPE_DEFINED - std::cout << "Computing bending scenario:" << bendingSimulationScenarioIndex << std::endl; -#endif job.nodalExternalForces.clear(); job.constrainedVertices.clear(); job.label = baseSimulationScenarioNames[BaseSimulationScenario::Bending] + "_" @@ -1248,8 +1284,8 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( } //// Dome - const double maxForceMagnitude_dome = getFullPatternMaxSimulationForce( - BaseSimulationScenario::Dome); + const double maxForceMagnitude_dome + = getFullPatternMaxSimulationForce(BaseSimulationScenario::Dome, wasSuccessful); const double minForceMagnitude_dome = 0; const int numberOfSimulationScenarios_dome = simulationScenariosResolution[BaseSimulationScenario::Dome]; @@ -1262,9 +1298,6 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( for (int domeSimulationScenarioIndex = 0; domeSimulationScenarioIndex < numberOfSimulationScenarios_dome; domeSimulationScenarioIndex++) { -#ifdef POLYSCOPE_DEFINED - std::cout << "Computing dome scenario:" << domeSimulationScenarioIndex << std::endl; -#endif job.constrainedVertices.clear(); job.nodalExternalForces.clear(); job.nodalForcedDisplacements.clear(); @@ -1278,8 +1311,8 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( } //// Saddle - const double maxForceMagnitude_saddle = getFullPatternMaxSimulationForce( - BaseSimulationScenario::Saddle); + const double maxForceMagnitude_saddle + = getFullPatternMaxSimulationForce(BaseSimulationScenario::Saddle, wasSuccessful); const double minForceMagnitude_saddle = 0; const int numberOfSimulationScenarios_saddle = simulationScenariosResolution[BaseSimulationScenario::Saddle]; @@ -1292,9 +1325,6 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( for (int saddleSimulationScenarioIndex = 0; saddleSimulationScenarioIndex < numberOfSimulationScenarios_saddle; saddleSimulationScenarioIndex++) { -#ifdef POLYSCOPE_DEFINED - std::cout << "Computing saddle scenario:" << saddleSimulationScenarioIndex << std::endl; -#endif job.constrainedVertices.clear(); job.nodalExternalForces.clear(); job.nodalForcedDisplacements.clear(); @@ -1426,10 +1456,15 @@ void ReducedModelOptimizer::optimize( m_pFullPatternSimulationMesh); // polyscope::removeAllStructures(); + results.baseTriangle = global.baseTriangle; + DRMSimulationModel::Settings simulationSettings; - simulationSettings.shouldDraw = false; + simulationSettings.maxDRMIterations = 200000; +// simulationSettings.shouldDraw = true; +// simulationSettings.isDebugMode = true; +// simulationSettings.debugModeStep = 100000; #ifdef POLYSCOPE_DEFINED - const bool drawFullPatternSimulationResults = true; + const bool drawFullPatternSimulationResults = false; if (drawFullPatternSimulationResults) { global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( ReducedPatternOptimization::Colors::fullInitial); @@ -1439,8 +1474,14 @@ void ReducedModelOptimizer::optimize( for (int simulationScenarioIndex : global.simulationScenarioIndices) { const std::shared_ptr &pFullPatternSimulationJob = global.fullPatternSimulationJobs[simulationScenarioIndex]; + SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, simulationSettings); + if (!fullPatternResults.converged) { + results.wasSuccessful = false; + return; + } + results.wasSuccessful = true; #ifdef POLYSCOPE_DEFINED if (drawFullPatternSimulationResults) { // SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation( @@ -1481,4 +1522,5 @@ void ReducedModelOptimizer::optimize( std::cout << "Running reduced model optimization" << std::endl; #endif runOptimization(optimizationSettings, results); + results.notes = optimizationNotes; } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index f9411ea..c50e902 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -28,6 +28,7 @@ class ReducedModelOptimizer m_fullPatternOppositeInterfaceViPairs; std::unordered_map nodeToSlot; std::unordered_map> slotToNode; + std::string optimizationNotes; public: constexpr static std::array simulationScenariosResolution = {10, 10, 20, 20, 20}; @@ -186,7 +187,7 @@ private: ReducedPatternOptimization::Results &results); std::vector getFullPatternMaxSimulationForces(); double getFullPatternMaxSimulationForce( - const ReducedPatternOptimization::BaseSimulationScenario &scenario); + const ReducedPatternOptimization::BaseSimulationScenario &scenario, bool &wasSuccessfull); }; void updateMesh(long n, const double *x); #endif // REDUCEDMODELOPTIMIZER_HPP