diff --git a/src/main.cpp b/src/main.cpp index 105fee1..837af09 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -114,21 +114,9 @@ int main(int argc, char *argv[]) { auto elapsed = std::chrono::duration_cast(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 diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index c11c6c2..1a9c4d3 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -4,6 +4,7 @@ #include "trianglepatterngeometry.hpp" #include "trianglepattterntopology.hpp" #include +#include 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 &desiredBaseSimulationScenarioIndices) { std::vector> 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>>( 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::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(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(job), settings); const double &desiredDisplacementValue = global.desiredMaxDisplacementValue; @@ -1138,7 +1152,6 @@ double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMa error = std::numeric_limits::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 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; } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index d7d6a16..54d54fd 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -35,6 +35,7 @@ class ReducedModelOptimizer SimulationJob &)>> constructBaseScenarioFunctions; std::vector scenarioIsSymmetrical; + int fullPatternNumberOfEdges; public: constexpr static std::array simulationScenariosResolution = {10, 10, 20, 20, 20};