diff --git a/CMakeLists.txt b/CMakeLists.txt index ed623b3..c733ada 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,19 +1,19 @@ cmake_minimum_required(VERSION 3.0) project(ReducedModelOptimization) +#message(STATUS "The compiler ${CMAKE_CXX_COMPILER}") set(CMAKE_CXX_STANDARD 20) set(CMAKE_CXX_STANDARD_REQUIRED ON) -message(STATUS "The compiler ${CMAKE_CXX_COMPILER}") #Add the project cmake scripts to the module path list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake) #Download external dependencies NOTE: If the user has one of these libs it shouldn't be downloaded again. include(${CMAKE_MODULE_PATH}/DownloadProject.cmake) -if (CMAKE_VERSION VERSION_LESS 3.2) - set(UPDATE_DISCONNECTED_IF_AVAILABLE "") -else() - set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1") -endif() +#if (CMAKE_VERSION VERSION_LESS 3.2) +# set(UPDATE_DISCONNECTED_IF_AVAILABLE "") +#else() +# set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1") +#endif() if(NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE "RelWithDebInfo") @@ -25,7 +25,6 @@ endif() ##Create directory for the external libraries file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR}) - #Add the project sources file(GLOB SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/src/*.hpp ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cpp) @@ -42,17 +41,22 @@ if(${CMAKE_BUILD_TYPE} STREQUAL "Release") else() set_target_properties(${PROJECT_NAME} PROPERTIES POSITION_INDEPENDENT_CODE TRUE) set(USE_POLYSCOPE true) - add_compile_definitions(POLYSCOPE_DEFINED) + if(${USE_POLYSCOPE}) + add_compile_definitions(POLYSCOPE_DEFINED) + endif() set(MYSOURCES_STATIC_LINK false) set(USE_ENSMALLEN true) if(${USE_ENSMALLEN}) add_compile_definitions(USE_ENSMALLEN) endif() +# if(NOT ${MYSOURCES_STATIC_LINK}) +# set( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fPIC" ) +# set( CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fPIC" ) +# endif() endif() #dlib -set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin) -file(MAKE_DIRECTORY ${DLIB_BIN_DIR}) +set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib) download_project(PROJ DLIB GIT_REPOSITORY https://github.com/davisking/dlib.git GIT_TAG master @@ -60,15 +64,16 @@ download_project(PROJ DLIB PREFIX ${EXTERNAL_DEPS_DIR} ${UPDATE_DISCONNECTED_IF_AVAILABLE} ) -add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BINARY_DIR}) +add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BIN_DIR}) if(${MYSOURCES_STATIC_LINK}) - target_link_libraries(${PROJECT_NAME} "-static" dlib::dlib) + target_link_libraries(${PROJECT_NAME} PUBLIC -static dlib::dlib) else() - target_link_libraries(${PROJECT_NAME} dlib::dlib) + target_link_libraries(${PROJECT_NAME} PUBLIC dlib::dlib) endif() add_compile_definitions(DLIB_DEFINED) set(MYSOURCES_SOURCE_DIR "/home/iason/Coding/Libraries/MySources") +set(MYSOURCES_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/MySourcesBinDir) if (EXISTS ${MYSOURCES_SOURCE_DIR}) else() ##MySources @@ -76,21 +81,20 @@ download_project(PROJ MYSOURCES GIT_REPOSITORY https://gitea-s2i2s.isti.cnr.it/manolas/MySources.git GIT_TAG master PREFIX ${EXTERNAL_DEPS_DIR} + BINARY_DIR ${MYSOURCES_BIN_DIR} ${UPDATE_DISCONNECTED_IF_AVAILABLE} ) endif() -add_subdirectory(${MYSOURCES_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}/MySourcesBinDir) +add_subdirectory(${MYSOURCES_SOURCE_DIR} ${MYSOURCES_BIN_DIR}) -target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_20) target_include_directories(${PROJECT_NAME} PUBLIC ${MYSOURCES_SOURCE_DIR} ) + if(${MYSOURCES_STATIC_LINK}) - target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen MySources) + target_link_libraries(${PROJECT_NAME} PUBLIC #[[-static]] MySources) else() - target_link_libraries(${PROJECT_NAME} Eigen3::Eigen MySources) + target_link_libraries(${PROJECT_NAME} PUBLIC MySources) endif() -##Eigen 3 NOTE: Eigen is required on the system the code is ran -find_package(Eigen3 3.3 REQUIRED) - +target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_20) diff --git a/main.py b/main.py index b60599a..2f52a6c 100644 --- a/main.py +++ b/main.py @@ -14,6 +14,7 @@ import shutil import sys from datetime import datetime from subprocess import check_output +from shutil import copyfile numberOfOptimizedPatterns=0 numberOfSkippedPatterns=0 @@ -51,8 +52,8 @@ def listener(q): def optimize(fullPatternFilepath, reducedPatternFilepath,translationalObjectiveWeight): """Call run(), catch exceptions.""" - # dirname = os.path.abspath(os.path.dirname(__file__)) - dirname="/home/iason/Coding/build/ReducedModelOptimization/Release" + dirname = os.path.abspath(os.path.dirname(__file__)) + # dirname="/home/iason/Coding/build/ReducedModelOptimization/Release" executableFilepath = os.path.join(dirname, 'ReducedModelOptimization') #putResultsTo=os.path.join(resultsDir,str(translationalObjectiveWeight)) #putResultsTo=os.path.join(putResultsTo,os.path.basename(os.path.dirname(os.path.dirname(fullPatternFilepath)))) @@ -83,10 +84,11 @@ def optimize(fullPatternFilepath, reducedPatternFilepath,translationalObjectiveW # Press the green button in the gutter to run the script. if __name__ == '__main__': dirOfThisFile = os.path.abspath(os.path.dirname(__file__)) + copyFrom="/home/iason/Coding/build/ReducedModelOptimization/Release/ReducedModelOptimization" + copyfile(copyFrom, os.path.join(dirOfThisFile,"ReducedModelOptimization")) fullPatternDirectory= "/home/iason/Coding/Projects/Approximating shapes with flat patterns/ReducedModelOptimization/TestSet/FullPatterns/selectionOfPatterns" - #resultsDir=os.path.join(dirOfThisFile,os.path.join('Results/OptimizationResults/',os.path.basename(fullPatternDirectory))) numberOfFunctionCalls=100000 - optimizationBatchName='variableComparison_allVars_joint_forConfirmation'+'_'+str(int(numberOfFunctionCalls/1000))+'k' + optimizationBatchName='variableComparison_AllVars_ensmallen'+'_'+str(int(numberOfFunctionCalls/1000))+'k' resultsDir=os.path.join(dirOfThisFile,os.path.join('Results/OptimizationResults/',optimizationBatchName)) #print(resultsDir) if not os.path.exists(resultsDir): @@ -119,7 +121,7 @@ if __name__ == '__main__': print("Number of function calls:",numberOfFunctionCalls) print("Full pattern test set directory:"+fullPatternDirectory) print("Total number of optimization jobs:"+ str(totalNumberOfPatterns)) - + print("Results directory:"+str(resultsDir)) print("Start time:", start_time) pool.starmap(optimize,jobs) print("Completed") diff --git a/src/main.cpp b/src/main.cpp index e68b54b..3af526a 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -55,13 +55,14 @@ int main(int argc, char *argv[]) { : 100; settings_optimization.normalizationStrategy = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon; - settings_optimization.splitGeometryMaterialOptimization = false; +#ifdef POLYSCOPE_DEFINED + settings_optimization.translationNormalizationParameter = 0; + settings_optimization.rotationNormalizationParameter = 0; +#else settings_optimization.translationNormalizationParameter = 3e-4; settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0); - // settings_optimization.translationNormalizationParameter = 1e-15; - // settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(1e-15); - settings_optimization.solverAccuracy = 1e-3; - // settings_optimization.solverAccuracy = 1e-1; +#endif + settings_optimization.solverAccuracy = 1e-1; settings_optimization.objectiveWeights.translational = std::atof(argv[4]); settings_optimization.objectiveWeights.rotational = 2 - std::atof(argv[4]); @@ -150,6 +151,45 @@ int main(int argc, char *argv[]) { } #ifdef POLYSCOPE_DEFINED + std::vector scenarioLabels( + optimizationResults.objectiveValue.perSimulationScenario_total.size()); + const double colorAxial = 1; + const double colorShear = 3; + const double colorBending = 5; + const double colorDome = 7; + const double colorSaddle = 9; + std::vector colors(optimizationResults.objectiveValue.perSimulationScenario_total.size()); + for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) { + scenarioLabels[scenarioIndex] + = optimizationResults.reducedPatternSimulationJobs[scenarioIndex]->getLabel(); + if (scenarioLabels[scenarioIndex].rfind("Axial", 0) == 0) { + colors[scenarioIndex] = colorAxial; + + } else if (scenarioLabels[scenarioIndex].rfind("Shear", 0) == 0) { + colors[scenarioIndex] = colorShear; + + } else if (scenarioLabels[scenarioIndex].rfind("Bending", 0) == 0) { + colors[scenarioIndex] = colorBending; + + } else if (scenarioLabels[scenarioIndex].rfind("Dome", 0) == 0) { + colors[scenarioIndex] = colorDome; + } else if (scenarioLabels[scenarioIndex].rfind("Saddle", 0) == 0) { + colors[scenarioIndex] = colorSaddle; + } else { + std::cerr << "Label could not be identified" << std::endl; + } + } + std::vector y(optimizationResults.objectiveValue.perSimulationScenario_total.size()); + for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) { + y[scenarioIndex] + = optimizationResults.objectiveValue.perSimulationScenario_rawTranslational[scenarioIndex] + + optimizationResults.objectiveValue.perSimulationScenario_rawRotational[scenarioIndex]; + // y[scenarioIndex] = optimizationResults.objectiveValue + // .perSimulationScenario_total[scenarioIndex]; + } + std::vector x = matplot::linspace(0, y.size() - 1, y.size()); + std::vector markerSizes(y.size(), 5); + SimulationResultsReporter::createPlot("scenario index", "error", x, y, markerSizes, colors); // optimizationResults.saveMeshFiles(); optimizationResults.draw(); #endif diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 5e2d933..d867e65 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -4,15 +4,19 @@ #include "trianglepatterngeometry.hpp" #include "trianglepattterntopology.hpp" #include +#include #include #include +#include #include +#include + +#define USE_SCENARIO_WEIGHTS using namespace ReducedPatternOptimization; struct GlobalOptimizationVariables { - // std::vector> fullPatternDisplacements; std::vector fullPatternResults; std::vector translationalDisplacementNormalizationValues; std::vector rotationalDisplacementNormalizationValues; @@ -25,14 +29,18 @@ struct GlobalOptimizationVariables std::vector objectiveValueHistory_iteration; std::vector objectiveValueHistory; std::vector plotColors; - std::array parametersInitialValue; - std::array optimizationInitialValue; + std::array + parametersInitialValue; + std::array + optimizationInitialValue; std::vector simulationScenarioIndices; double minY{DBL_MAX}; std::vector minX; int numOfSimulationCrashes{false}; int numberOfFunctionCalls{0}; - int numberOfOptimizationParameters{5}; + int numberOfOptimizationParameters{5}; ReducedPatternOptimization::Settings optimizationSettings; vcg::Triangle3 baseTriangle; //Variables for finding the full pattern simulation forces @@ -51,8 +59,7 @@ struct GlobalOptimizationVariables functions_updateReducedPatternParameter; std::vector xMin; std::vector xMax; - std::vector scenarioEqualizationWeight; - std::vector scenarioUserWeights; + std::vector scenarioWeights; } global; double ReducedModelOptimizer::computeDisplacementError( @@ -128,23 +135,35 @@ double ReducedModelOptimizer::computeError( const std::unordered_map &reducedToFullInterfaceViMap, const double &normalizationFactor_translationalDisplacement, - const double &normalizationFactor_rotationalDisplacement) + const double &normalizationFactor_rotationalDisplacement, + const double &scenarioWeight) { const double translationalError = computeDisplacementError(simulationResults_fullPattern.displacements, simulationResults_reducedPattern.displacements, reducedToFullInterfaceViMap, normalizationFactor_translationalDisplacement); + // const double translationalError + // = computeRawTranslationalError(simulationResults_fullPattern.displacements, + // simulationResults_reducedPattern.displacements, + // reducedToFullInterfaceViMap); // std::cout << "normalization factor:" << normalizationFactor_rotationalDisplacement << std::endl; - // std::cout << "trans error:" << translationalError << std::endl; +// std::cout << "trans error:" << translationalError << std::endl; const double rotationalError = computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, simulationResults_reducedPattern.rotationalDisplacementQuaternion, reducedToFullInterfaceViMap, normalizationFactor_rotationalDisplacement); - return global.optimizationSettings.objectiveWeights.translational * translationalError - + global.optimizationSettings.objectiveWeights.rotational * rotationalError; + // computeRawRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, + // simulationResults_reducedPattern.rotationalDisplacementQuaternion, + // reducedToFullInterfaceViMap); + // std::cout << "rot error:" << rotationalError<< std::endl; + const double simulationError = global.optimizationSettings.objectiveWeights.translational + * translationalError + + global.optimizationSettings.objectiveWeights.rotational + * rotationalError; + return simulationError * simulationError * scenarioWeight; } //double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3, @@ -162,13 +181,17 @@ struct EnsmallenOptimizationObjective std::vector x(x_arma.begin(), x_arma.end()); for (int xi = 0; xi < x.size(); xi++) { if (x[xi] > global.xMax[xi]) { + //#ifdef POLYSCOPE_DEFINED + // std::cout << "Out of range" << std::endl; + //#endif return std::numeric_limits::max(); - x[xi] = global.xMax[xi]; - // std::cout << "Out of range" << std::endl; + // x[xi] = global.xMax[xi]; } else if (x[xi] < global.xMin[xi]) { + //#ifdef POLYSCOPE_DEFINED + // std::cout << "Out of range" << std::endl; + //#endif return std::numeric_limits::max(); - x[xi] = global.xMin[xi]; - // std::cout << "Out of range" << std::endl; + // x[xi] = global.xMin[xi]; } } return ReducedModelOptimizer::objective(x); @@ -185,7 +208,7 @@ double ReducedModelOptimizer::objective(const dlib::matrix &x) double ReducedModelOptimizer::objective(const double &xValue) { - return objective({xValue}); + return objective(std::vector({xValue})); } double ReducedModelOptimizer::objective(const std::vector &x) @@ -212,13 +235,8 @@ double ReducedModelOptimizer::objective(const std::vector &x) // polyscope::show(); // global.reducedPatternSimulationJobs[0]->pMesh->unregister(); - // std::cout << e.axialConstFact|A,I2,I3,J,r,thetaor << " " << e.torsionConstFactor << " " - // << e.firstBendingConstFactor << " " << - // e.secondBendingConstFactor - // << std::endl; - // run simulations - double totalError = 0; + std::vector simulationErrorsPerScenario(global.simulationScenarioIndices.size()); LinearSimulationModel simulator; simulator.setStructure(pReducedPatternSimulationMesh); // simulator.initialize(); @@ -238,73 +256,60 @@ double ReducedModelOptimizer::objective(const std::vector &x) SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob); // std::string filename; if (!reducedModelResults.converged) { - totalError += std::numeric_limits::max(); + simulationErrorsPerScenario[simulationScenarioIndex] + = std::numeric_limits::max(); global.numOfSimulationCrashes++; #ifdef POLYSCOPE_DEFINED std::cout << "Failed simulation" << std::endl; #endif } else { - // double simulationScenarioError; - // constexpr bool usePotentialEnergy = false; - // if (usePotentialEnergy) { - // simulationScenarioError = std::abs( - // reducedModelResults.internalPotentialEnergy - // - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy); - // } else { +// double simulationScenarioError = 0; +// constexpr bool usePotentialEnergy = false; +// if (usePotentialEnergy) { +// simulationScenarioError += std::abs( +// reducedModelResults.internalPotentialEnergy +// - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy +// / global.fullPatternResults[simulationScenarioIndex] +// .internalPotentialEnergy); +// } else { +#ifdef POLYSCOPE_DEFINED +#ifdef USE_SCENARIO_WEIGHTS + const double scenarioWeight = global.scenarioWeights[simulationScenarioIndex]; +#else + const double scenarioWeight = 1; +#endif +#else + const double scenarioWeight = 1; +#endif const double simulationScenarioError = computeError( global.fullPatternResults[simulationScenarioIndex], reducedModelResults, global.reducedToFullInterfaceViMap, global.translationalDisplacementNormalizationValues[simulationScenarioIndex], - global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); + global.rotationalDisplacementNormalizationValues[simulationScenarioIndex], + scenarioWeight); + + simulationErrorsPerScenario[simulationScenarioIndex] = simulationScenarioError; // } - //#ifdef POLYSCOPE_DEFINED - // reducedJob->pMesh->registerForDrawing(Colors::reducedInitial); - // reducedModelResults.registerForDrawing(Colors::reducedDeformed); - // global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); - // global.fullPatternResults[simulationScenarioIndex].registerForDrawing( - // Colors::fullDeformed); - // polyscope::show(); - // reducedModelResults.unregister(); - // global.pFullPatternSimulationMesh->unregister(); - // global.fullPatternResults[simulationScenarioIndex].unregister(); - //#endif - // if (global.optimizationSettings.normalizationStrategy != - // NormalizationStrategy::Epsilon && - // simulationScenarioError > 1) { - // std::cout << "Simulation scenario " - // << - // simulationScenarioStrings[simulationScenarioIndex] - // << " results in an error bigger than one." << - // std::endl; - // for (size_t parameterIndex = 0; parameterIndex < n; - // parameterIndex++) { - // std::cout << "x[" + std::to_string(parameterIndex) + "]=" - // << x[parameterIndex] << std::endl; - // } - // } - - // #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 - totalError += /*global.scenarioUserWeights[simulationScenarioIndex] - **/ - /*global.scenarioEqualizationWeight[simulationScenarioIndex] - **/ - simulationScenarioError * simulationScenarioError; +// #ifdef POLYSCOPE_DEFINED +// reducedJob->pMesh->registerForDrawing(Colors::reducedInitial); +// reducedModelResults.registerForDrawing(Colors::reducedDeformed); +// global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); +// global.fullPatternResults[simulationScenarioIndex].registerForDrawing( +// Colors::fullDeformed); +// polyscope::show(); +// reducedModelResults.unregister(); +// global.pFullPatternSimulationMesh->unregister(); +// global.fullPatternResults[simulationScenarioIndex].unregister(); +// #endif } }); + // double totalError = std::accumulate(simulationErrorsPerScenario.begin(), + // simulationErrorsPerScenario.end(), + // 0); + double totalError = std::reduce(std::execution::par_unseq, + simulationErrorsPerScenario.begin(), + simulationErrorsPerScenario.end()); // std::cout << totalError << std::endl; // global.objectiveValueHistory.push_back(totalError); // global.plotColors.push_back(10); @@ -313,13 +318,15 @@ double ReducedModelOptimizer::objective(const std::vector &x) global.minY = totalError; // global.objectiveValueHistory.push_back(totalError); // global.objectiveValueHistory_iteration.push_back(global.numberOfFunctionCalls); - // std::cout << "New best:" << totalError << std::endl; +#ifdef POLYSCOPE_DEFINED + std::cout << "New best:" << totalError << std::endl; // // global.minX.assign(x.begin(), x.begin() + n); - // std::cout.precision(17); - // for (int i = 0; i < x.size(); i++) { - // std::cout << x(i) << " "; - // } - // std::cout << std::endl; + std::cout.precision(17); + for (int i = 0; i < x.size(); i++) { + std::cout << x[i] << " "; + } + std::cout << std::endl; +#endif // global.objectiveValueHistoryY.push_back(std::log(totalError)); // global.objectiveValueHistoryX.push_back(global.numberOfFunctionCalls); // global.plotColors.push_back(0.1); @@ -360,11 +367,13 @@ void ReducedModelOptimizer::createSimulationMeshes( pFullPatternSimulationMesh = std::make_shared(fullModel); pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002}); pFullPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus); + pFullPatternSimulationMesh->reset(); // Reduced pattern pReducedPatternSimulationMesh = std::make_shared(reducedModel); pReducedPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002}); pReducedPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus); + pReducedPatternSimulationMesh->reset(); } void ReducedModelOptimizer::createSimulationMeshes(PatternGeometry &fullModel, @@ -616,7 +625,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() [](const double &newI2, std::shared_ptr &pReducedPatternSimulationMesh) { for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; - e.inertia.I2 = newI2; + e.dimensions.inertia.I2 = newI2; e.updateRigidity(); } }; @@ -624,7 +633,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() [](const double &newI3, std::shared_ptr &pReducedPatternSimulationMesh) { for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; - e.inertia.I3 = newI3; + e.dimensions.inertia.I3 = newI3; e.updateRigidity(); } }; @@ -632,7 +641,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() [](const double &newJ, std::shared_ptr &pReducedPatternSimulationMesh) { for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; - e.inertia.J = newJ; + e.dimensions.inertia.J = newJ; e.updateRigidity(); } }; @@ -657,22 +666,23 @@ void ReducedModelOptimizer::initializeOptimizationParameters( global.optimizationInitialValue[optimizationParameterIndex] = 1; break; case A: - global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].A; + global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0] + .dimensions.A; global.optimizationInitialValue[optimizationParameterIndex] = 1; break; case I2: - global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0] - .inertia.I2; + global.parametersInitialValue[optimizationParameterIndex] + = mesh->elements[0].dimensions.inertia.I2; global.optimizationInitialValue[optimizationParameterIndex] = 1; break; case I3: - global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0] - .inertia.I3; + global.parametersInitialValue[optimizationParameterIndex] + = mesh->elements[0].dimensions.inertia.I3; global.optimizationInitialValue[optimizationParameterIndex] = 1; break; case J: - global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0] - .inertia.J; + global.parametersInitialValue[optimizationParameterIndex] + = mesh->elements[0].dimensions.inertia.J; global.optimizationInitialValue[optimizationParameterIndex] = 1; break; case R: @@ -742,6 +752,8 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv const Settings &settings, ReducedPatternOptimization::Results &results) { + std::shared_ptr &pReducedPatternSimulationMesh + = global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh; //Number of crashes // results.numberOfSimulationCrashes = global.numOfSimulationCrashes; //Value of optimal objective Y @@ -752,8 +764,7 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv << std::endl; } //Optimal X values - results.optimalXNameValuePairs.resize(settings.parameterRanges.size()); - std::vector optimalX(settings.parameterRanges.size()); + results.optimalXNameValuePairs.resize(NumberOfOptimizationParameters); for (int optimizationParameterIndex = E; optimizationParameterIndex != NumberOfOptimizationParameters; optimizationParameterIndex++) { @@ -762,7 +773,7 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv < 1e-10) { continue; } - if (optimizationParameterIndex < settings.parameterRanges.size() - 2) { + if (optimizationParameterIndex != OptimizationParameterIndex::R && optimizationParameterIndex != OptimizationParameterIndex::Theta ) { results.optimalXNameValuePairs[optimizationParameterIndex] = std::make_pair(settings.parameterRanges[optimizationParameterIndex].label, optimalObjective.x[optimizationParameterIndex] @@ -774,24 +785,22 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv optimalObjective.x[optimizationParameterIndex]); } - optimalX[optimizationParameterIndex] = optimalObjective.x[optimizationParameterIndex]; + global.functions_updateReducedPatternParameter[optimizationParameterIndex]( + results.optimalXNameValuePairs[optimizationParameterIndex].second, + pReducedPatternSimulationMesh); #ifdef POLYSCOPE_DEFINED std::cout << results.optimalXNameValuePairs[optimizationParameterIndex].first << ":" - << optimalX[optimizationParameterIndex] << " "; + << results.optimalXNameValuePairs[optimizationParameterIndex].second << " "; #endif } #ifdef POLYSCOPE_DEFINED std::cout << std::endl; #endif + pReducedPatternSimulationMesh->reset(); results.fullPatternYoungsModulus = youngsModulus; // Compute obj value per simulation scenario and the raw objective value // updateMeshFunction(optimalX); - std::shared_ptr &pReducedPatternSimulationMesh - = global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh; - function_updateReducedPattern(std::vector(optimalObjective.x.begin(), - optimalObjective.x.end()), - pReducedPatternSimulationMesh); // results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios); //TODO:use push_back it will make the code more readable LinearSimulationModel simulator; @@ -816,12 +825,22 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv SimulationResults reducedModelResults = simulator.executeSimulation( global.reducedPatternSimulationJobs[simulationScenarioIndex]); +#ifdef POLYSCOPE_DEFINED +#ifdef USE_SCENARIO_WEIGHTS + const double scenarioWeight = global.scenarioWeights[simulationScenarioIndex]; +#else + const double scenarioWeight = 1; +#endif +#else + const double scenarioWeight = 1; +#endif results.objectiveValue.perSimulationScenario_total[i] = computeError( global.fullPatternResults[simulationScenarioIndex], reducedModelResults, global.reducedToFullInterfaceViMap, global.translationalDisplacementNormalizationValues[simulationScenarioIndex], - global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); + global.rotationalDisplacementNormalizationValues[simulationScenarioIndex], + scenarioWeight); // results.objectiveValue.total += results.objectiveValue.perSimulationScenario_total[i]; //Raw translational @@ -911,30 +930,28 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv // results.draw(); } -std::vector> +std::array ReducedModelOptimizer::computeFullPatternMaxSimulationForces( - const std::vector &desiredBaseSimulationScenario) + const std::vector &desiredBaseSimulationScenario) const { - std::vector> fullPatternMaxSimulationForces; - fullPatternMaxSimulationForces.reserve(desiredBaseSimulationScenario.size()); + std::array fullPatternMaxSimulationForces; for (const BaseSimulationScenario &scenario : desiredBaseSimulationScenario) { const double maxForce = computeFullPatternMaxSimulationForce(scenario); - fullPatternMaxSimulationForces.emplace_back(std::make_pair(scenario, maxForce)); + fullPatternMaxSimulationForces[scenario] = maxForce; } return fullPatternMaxSimulationForces; } -std::vector> +std::array ReducedModelOptimizer::getFullPatternMaxSimulationForces( - const std::vector &desiredBaseSimulationScenarioIndices) + const std::vector &desiredBaseSimulationScenarioIndices, + const std::filesystem::path &intermediateResultsDirectoryPath) { - std::vector> fullPatternSimulationScenarioMaxMagnitudes; + std::array fullPatternSimulationScenarioMaxMagnitudes; #ifdef POLYSCOPE_DEFINED - const std::filesystem::path forceMagnitudesDirectoryPath(std::filesystem::current_path() - .parent_path() - .append("IntermediateResults") - .append("ForceMagnitudes")); + const std::filesystem::path forceMagnitudesDirectoryPath( + std::filesystem::path(intermediateResultsDirectoryPath).append("ForceMagnitudes")); std::filesystem::path patternMaxForceMagnitudesFilePath( std::filesystem::path(forceMagnitudesDirectoryPath) .append(m_pFullPatternSimulationMesh->getLabel() + ".json")); @@ -945,13 +962,8 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces( std::ifstream ifs(patternMaxForceMagnitudesFilePath.string()); ifs >> json; fullPatternSimulationScenarioMaxMagnitudes - = static_cast>>( - json.at("maxMagn")); - const bool shouldRecompute = fullPatternSimulationScenarioMaxMagnitudes.size() - != desiredBaseSimulationScenarioIndices.size(); - if (!shouldRecompute) { - return fullPatternSimulationScenarioMaxMagnitudes; - } + = static_cast>(json.at("maxMagn")); + return fullPatternSimulationScenarioMaxMagnitudes; } #endif fullPatternSimulationScenarioMaxMagnitudes = computeFullPatternMaxSimulationForces( @@ -985,16 +997,18 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, #endif #ifdef USE_ENSMALLEN #else - double (*objF)(const dlib::matrix &) = &objective; #endif enum OptimizationParameterComparisonScenarioIndex { - AllVar, - YMAndGeo, - noYM, - SplittedMatGeo, - SplittedYM_before, - SplittedYM_after, - SplittedGeoYMMat, + AllVar, //ok + GeoYM, //ok + noYM, //ok + YMMat_Geo, //ok + YM_MatGeo, //executing + MatGeo_YM, //ok + Geo_YM_Mat, //ok + YM_Geo_Mat, //ok + Geo_Mat, //ok + YMGeo_Mat, NumberOfScenarios }; const std::vector>> scenarioParameters = @@ -1002,38 +1016,57 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, std::vector>> scenarioParameters( NumberOfScenarios); scenarioParameters[AllVar] = {{E, A, I2, I3, J, R, Theta}}; - scenarioParameters[YMAndGeo] = {{R, Theta, E}}; + scenarioParameters[GeoYM] = {{R, Theta, E}}; scenarioParameters[noYM] = {{A, I2, I3, J, R, Theta}}; - scenarioParameters[SplittedMatGeo] = {{E, A, I2, I3, J}, {R, Theta}}; - scenarioParameters[SplittedYM_before] = {{E}, {A, I2, I3, J, R, Theta}}; - scenarioParameters[SplittedYM_after] = {{A, I2, I3, J, R, Theta}, {E}}; - scenarioParameters[SplittedGeoYMMat] = {{R, Theta}, {E}, {A, I2, I3, J}}; + scenarioParameters[YMMat_Geo] = {{E, A, I2, I3, J}, {R, Theta}}; + scenarioParameters[YM_MatGeo] = {{E}, {A, I2, I3, J, R, Theta}}; + scenarioParameters[MatGeo_YM] = {{A, I2, I3, J, R, Theta}, {E}}; + scenarioParameters[Geo_YM_Mat] = {{R, Theta}, {E}, {A, I2, I3, J}}; + scenarioParameters[YM_Geo_Mat] = {{E},{R, Theta}, {A, I2, I3, J}}; + scenarioParameters[Geo_Mat] = {{R, Theta}, {A, I2, I3, J}}; + scenarioParameters[YMGeo_Mat] = {{E, R, Theta}, {A, I2, I3, J}}; return scenarioParameters; + }(); - constexpr OptimizationParameterComparisonScenarioIndex scenario = noYM; + constexpr OptimizationParameterComparisonScenarioIndex scenario = YM_MatGeo; const std::vector> scenarioParameterGroups = scenarioParameters[scenario]; + const int totalNumberOfOptimizationParameters + = std::accumulate(scenarioParameterGroups.begin(), + scenarioParameterGroups.end(), + 0, + [](const int &sum, + const std::vector ¶meterGroup) { + return sum + parameterGroup.size(); + }); - //TODO:Ensure that the reduced pattern mesh has the initial parameter values before starting the optimization - FunctionEvaluation optimalResult; + FunctionEvaluation optimization_optimalResult; + optimization_optimalResult.x.resize(NumberOfOptimizationParameters,0); + for (int optimizationParameterIndex = E; + optimizationParameterIndex != NumberOfOptimizationParameters; + optimizationParameterIndex++) { + optimization_optimalResult.x[optimizationParameterIndex]=global.parametersInitialValue[optimizationParameterIndex]; + } for (const std::vector ¶meterGroup : scenarioParameterGroups) { + FunctionEvaluation parameterGroup_optimalResult; //Set update function. TODO: Make this function immutable by defining it once and using the global variable to set parameterGroup function_updateReducedPattern = [&](const std::vector &x, std::shared_ptr &pMesh) { - for (int xIndex = 0; xIndex < x.size(); xIndex++) { + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; // const double parameterInitialValue=optimizationSettings.parameterRanges[parameterIndex].initialValue; - const double parameterInitialValue = global.parametersInitialValue[parameterIndex]; const double parameterNewValue = [&]() { - if (parameterIndex == R - || parameterIndex - == Theta) { //NOTE:Here I explore the geometry parameters linearly + if (parameterIndex == R || parameterIndex == Theta) { return x[xIndex] /*+ parameterInitialValue*/; } //and the material parameters exponentially(?).TODO: Check what happens if I make all linear + const double parameterInitialValue + = global.parametersInitialValue[parameterIndex]; return x[xIndex] * parameterInitialValue; }(); + // std::cout << "Optimization parameter:" << parameterIndex << std::endl; + // std::cout << "New value:" << parameterNewValue << std::endl; global.functions_updateReducedPatternParameter[parameterIndex](parameterNewValue, pMesh); } @@ -1072,28 +1105,29 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, // ens::SPSA optimizer; // arma::mat xMin_arma(global.xMin); // arma::mat xMax_arma(global.xMax); - // ens::LBestPSO optimizer(64, xMin_arma, xMax_arma, 1000); - ens::LBestPSO optimizer; + // ens::LBestPSO optimizer(64, xMin_arma, xMax_arma, 1000); + ens::LBestPSO optimizer(64, 1, 1, 3000, 350, 1e-1, 2.05, 2.15); + // ens::LBestPSO optimizer; const double minima = optimizer.Optimize(optimizationFunction, x); - for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { - if (x[xIndex] > global.xMax[xIndex]) { - x[xIndex] = global.xMax[xIndex]; - } else if (x[xIndex] < global.xMin[xIndex]) { - x[xIndex] = global.xMin[xIndex]; - } - } + // for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + // if (x[xIndex] > global.xMax[xIndex]) { + // x[xIndex] = global.xMax[xIndex]; + // } else if (x[xIndex] < global.xMin[xIndex]) { + // x[xIndex] = global.xMin[xIndex]; + // } + // } #ifdef POLYSCOPE_DEFINED std::cout << "optimal x:" << "\n" << x << std::endl; std::cout << "Minima:" << minima << std::endl; + std::cout << "optimal objective:" << optimizationFunction.Evaluate(x) << std::endl; #endif - optimalResult.x.clear(); - optimalResult.x.resize(parameterGroup.size()); - std::copy(x.begin(), x.end(), optimalResult.x.begin()); - optimalResult.y = minima; - + parameterGroup_optimalResult.x.clear(); + parameterGroup_optimalResult.x.resize(parameterGroup.size()); + std::copy(x.begin(), x.end(), parameterGroup_optimalResult.x.begin()); + parameterGroup_optimalResult.y = minima; #else //Set min max values dlib::matrix xMin_dlib(parameterGroup.size()); @@ -1104,14 +1138,39 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, xMin_dlib(xIndex) = settings.parameterRanges[parameterIndex].min; xMax_dlib(xIndex) = settings.parameterRanges[parameterIndex].max; } - const dlib::function_evaluation optimalResult_dlib - = dlib::find_min_global(objF, - xMin_dlib, - xMax_dlib, - dlib::max_function_calls(settings.numberOfFunctionCalls), - std::chrono::hours(24 * 365 * 290), - settings.solverAccuracy); + const int optimizationNumberOfFunctionCalls = static_cast( + (static_cast(parameterGroup.size()) / totalNumberOfOptimizationParameters) + * settings.numberOfFunctionCalls); + const dlib::function_evaluation optimalResult_dlib = [&]() { + if (parameterGroup.size() == 1) { + dlib::function_evaluation result; + double optimalX = global.optimizationInitialValue[parameterGroup[0]]; + double (*singleVariableObjective)(const double &) = &objective; + try { + result.y = dlib::find_min_single_variable(singleVariableObjective, + optimalX, + xMin_dlib(0), + xMax_dlib(0), + 1e-5, + optimizationNumberOfFunctionCalls); + } catch (const std::exception &e) { // reference to the base of a polymorphic object +#ifdef POLYSCOPE_DEFINED + std::cout << e.what() << std::endl; +#endif + } + result.x.set_size(1); + result.x(0) = optimalX; + return result; + } + double (*objF)(const dlib::matrix &) = &objective; + return dlib::find_min_global(objF, + xMin_dlib, + xMax_dlib, + dlib::max_function_calls(optimizationNumberOfFunctionCalls), + std::chrono::hours(24 * 365 * 290), + settings.solverAccuracy); + }(); // constexpr bool useBOBYQA = false; // if (useBOBYQA) { // const size_t npt = 2 * global.numberOfOptimizationParameters; @@ -1141,20 +1200,40 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, // rhoend, // bobyqa_maxFunctionCalls); // } - optimalResult.x.clear(); - optimalResult.x.resize(parameterGroup.size()); - std::copy(optimalResult_dlib.x.begin(), optimalResult_dlib.x.end(), optimalResult.x.begin()); - optimalResult.y = optimalResult_dlib.y; + + parameterGroup_optimalResult.x.clear(); + parameterGroup_optimalResult.x.resize(parameterGroup.size()); + std::copy(optimalResult_dlib.x.begin(), + optimalResult_dlib.x.end(), + parameterGroup_optimalResult.x.begin()); + parameterGroup_optimalResult.y = optimalResult_dlib.y; #endif + function_updateReducedPattern( - std::vector(optimalResult.x.begin(), optimalResult.x.end()), + std::vector(parameterGroup_optimalResult.x.begin(), parameterGroup_optimalResult.x.end()), global.reducedPatternSimulationJobs[0] ->pMesh); //TODO: Check if its ok to update only the changed parameters - // std::cout << "GLOBAL MIN:" << global.minY << std::endl; - // std::cout << "opt res y:" << optimalResult.y << std::endl; +#ifdef POLYSCOPE_DEFINED + std::cout << "Optimal x:" + << "\n"; + for (const auto &x : parameterGroup_optimalResult.x) { + std::cout << x << std::endl; + } + std::cout << "optimal objective:" << objective(parameterGroup_optimalResult.x) << std::endl; + // std::cout << "Minima:" << minima << std::endl; + std::cout << "GLOBAL MIN:" << global.minY << std::endl; + std::cout << "opt res y:" << parameterGroup_optimalResult.y << std::endl; +#endif + +optimization_optimalResult.y=parameterGroup_optimalResult.y; + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; +optimization_optimalResult.x[parameterIndex]=parameterGroup_optimalResult.x[xIndex]; +} + } - getResults(optimalResult, settings, results); + getResults(optimization_optimalResult, settings, results); } void ReducedModelOptimizer::constructAxialSimulationScenario( @@ -1443,7 +1522,7 @@ struct ForceMagnitudeOptimization #endif double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( - const BaseSimulationScenario &scenario) + const BaseSimulationScenario &scenario) const { double forceMagnitude = 1; double minimumError; @@ -1537,9 +1616,61 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( return forceMagnitude; } +void ReducedModelOptimizer::computeScenarioWeights( + const std::vector &baseSimulationScenarios, + const std::array &maxForceMagnitudePerBaseScenario) +{ + assert(maxForceMagnitudePerBaseScenario.size() == NumberOfBaseSimulationScenarios); + const double userWeightsSum = std::accumulate(baseScenarioWeights.begin(), + baseScenarioWeights.end(), + 0); + const int numberOfBaseScenarios = baseSimulationScenarios.size(); + + Eigen::VectorXd baseScenario_equalizationWeights(numberOfBaseScenarios); + Eigen::VectorXd baseScenario_userWeights(numberOfBaseScenarios); + Eigen::Matrix baseScenario_maxMagnitudeWeights( + maxForceMagnitudePerBaseScenario.data()); + for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) { + const double baseScenarioEqualizationWeight + = totalNumberOfSimulationScenarios + / static_cast(simulationScenariosResolution[baseScenario]); + baseScenario_equalizationWeights(baseScenario) = baseScenarioEqualizationWeight; + const double baseScenarioUserWeight = static_cast(baseScenarioWeights[baseScenario]) + / userWeightsSum; + baseScenario_userWeights(baseScenario) = baseScenarioUserWeight; + } + + baseScenario_maxMagnitudeWeights.normalize(); + baseScenario_equalizationWeights.normalize(); + baseScenario_userWeights.normalize(); + + global.scenarioWeights.resize( + totalNumberOfSimulationScenarios); //This essentially holds the base scenario weights but I use totalNumberOfSimulationScenarios elements instead in order to have O(1) access time during the optimization + for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) { + const int baseSimulationScenarioIndexOffset + = std::accumulate(simulationScenariosResolution.begin(), + simulationScenariosResolution.begin() + baseScenario, + 0); + const double baseScenarioEqualizationWeight = baseScenario_equalizationWeights(baseScenario); + // const double baseScenarioUserWeight = baseScenario_userWeights(baseScenario); + const double baseScenarioUserWeight = baseScenario_userWeights(baseScenario); + + const double baseScenarioMaxMagnitudeWeight = baseScenario_maxMagnitudeWeights(baseScenario); + for (int simulationScenarioIndex = 0; + simulationScenarioIndex < simulationScenariosResolution[baseScenario]; + simulationScenarioIndex++) { + const int globalScenarioIndex = baseSimulationScenarioIndexOffset + + simulationScenarioIndex; + const double scenarioWeight = /*baseScenarioUserWeight +*/ baseScenarioEqualizationWeight + + baseScenarioMaxMagnitudeWeight; + global.scenarioWeights[globalScenarioIndex] = scenarioWeight; + } + } +} + std::vector> ReducedModelOptimizer::createFullPatternSimulationJobs( const std::shared_ptr &pMesh, - const std::vector> &scenarioMaxForceMagnitudePairs) + const std::array &baseScenarioMaxForceMagnitudes) const { std::vector> scenarios; scenarios.resize(totalNumberOfSimulationScenarios); @@ -1547,18 +1678,10 @@ std::vector> ReducedModelOptimizer::createFullPat SimulationJob job; job.pMesh = pMesh; - global.scenarioEqualizationWeight.resize( - totalNumberOfSimulationScenarios); //This essentially holds the base scenario weights but I use totalNumberOfSimulationScenarios elements instead in order to have O(1) access time during the optimization - global.scenarioUserWeights.resize(totalNumberOfSimulationScenarios); - const double userWeightsSum = std::accumulate(baseScenarioWeights.begin(), - baseScenarioWeights.end(), - 0); - for (std::pair scenarioMaxForceMagnitudePair : - scenarioMaxForceMagnitudePairs) { - const BaseSimulationScenario scenario = scenarioMaxForceMagnitudePair.first; - const double maxForceMagnitude = scenarioMaxForceMagnitudePair.second; - const int numberOfSimulationScenarios = simulationScenariosResolution[scenario]; - const double minForceMagnitude = scenarioIsSymmetrical[scenario] + for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios; baseScenario++) { + const double maxForceMagnitude = baseScenarioMaxForceMagnitudes[baseScenario]; + const int numberOfSimulationScenarios = simulationScenariosResolution[baseScenario]; + const double minForceMagnitude = scenarioIsSymmetrical[baseScenario] ? maxForceMagnitude / numberOfSimulationScenarios : -maxForceMagnitude; const double forceMagnitudeStep = numberOfSimulationScenarios == 1 @@ -1567,32 +1690,28 @@ std::vector> ReducedModelOptimizer::createFullPat / (numberOfSimulationScenarios); const int baseSimulationScenarioIndexOffset = std::accumulate(simulationScenariosResolution.begin(), - simulationScenariosResolution.begin() + scenario, + simulationScenariosResolution.begin() + baseScenario, 0); - const double baseScenarioEqualizationWeight = static_cast( - numberOfSimulationScenarios) - / totalNumberOfSimulationScenarios; - const double baseScenarioUserWeight = static_cast(baseScenarioWeights[scenario]) - / userWeightsSum; for (int simulationScenarioIndex = 0; simulationScenarioIndex < numberOfSimulationScenarios; simulationScenarioIndex++) { + const int scenarioIndex = baseSimulationScenarioIndexOffset + simulationScenarioIndex; + if (baseScenarioMaxForceMagnitudes[baseScenario] == -1) { + scenarios[scenarioIndex] = nullptr; + continue; + } job.nodalExternalForces.clear(); job.constrainedVertices.clear(); job.nodalForcedDisplacements.clear(); - job.label = baseSimulationScenarioNames[scenario] + "_" + job.label = baseSimulationScenarioNames[baseScenario] + "_" + std::to_string(simulationScenarioIndex); const double forceMagnitude = (forceMagnitudeStep * simulationScenarioIndex + minForceMagnitude); - constructBaseScenarioFunctions[scenario](forceMagnitude, - m_fullPatternOppositeInterfaceViPairs, - job); + constructBaseScenarioFunctions[baseScenario](forceMagnitude, + m_fullPatternOppositeInterfaceViPairs, + job); - scenarios[baseSimulationScenarioIndexOffset + simulationScenarioIndex] - = std::make_shared(job); - global.scenarioEqualizationWeight[simulationScenarioIndex] - = baseScenarioEqualizationWeight; - global.scenarioUserWeights[simulationScenarioIndex] = baseScenarioUserWeight; + scenarios[scenarioIndex] = std::make_shared(job); } } @@ -1712,13 +1831,25 @@ void ReducedModelOptimizer::optimize( global.optimizationSettings = optimizationSettings; global.pFullPatternSimulationMesh = m_pFullPatternSimulationMesh; - std::vector> fullPatternSimulationScenarioMaxMagnitudes - = getFullPatternMaxSimulationForces(desiredBaseSimulationScenarioIndices); +#ifdef POLYSCOPE_DEFINED + const std::filesystem::path intermediateResultsDirectoryPath( + std::filesystem::current_path().parent_path().append("IntermediateResults")); +#else + const std::filesystem::path intermediateResultsDirectoryPath( + std::filesystem::current_path().append("IntermediateResults")); + +#endif + std::array fullPatternSimulationScenarioMaxMagnitudes + = getFullPatternMaxSimulationForces(desiredBaseSimulationScenarioIndices, + intermediateResultsDirectoryPath); global.fullPatternSimulationJobs = createFullPatternSimulationJobs(m_pFullPatternSimulationMesh, fullPatternSimulationScenarioMaxMagnitudes); // polyscope::removeAllStructures(); + computeScenarioWeights(desiredBaseSimulationScenarioIndices, + fullPatternSimulationScenarioMaxMagnitudes); + results.baseTriangle = global.baseTriangle; DRMSimulationModel::Settings simulationSettings; @@ -1745,33 +1876,43 @@ void ReducedModelOptimizer::optimize( const std::shared_ptr &pFullPatternSimulationJob = global.fullPatternSimulationJobs[simulationScenarioIndex]; - // std::filesystem::path patternMaxForceMagnitudesFilePath( - // std::filesystem::path(forceMagnitudesDirectoryPath) - // .append(pFullPatternSimulationJob->getLabel() + ".json")); - // const bool fullPatternScenarioMagnitudesExist = std::filesystem::exists( - // patternMaxForceMagnitudesFilePath); - // if (fullPatternScenarioMagnitudesExist) { - // nlohmann::json json; - // std::ifstream ifs(patternMaxForceMagnitudesFilePath.string()); - // ifs >> json; - // fullPatternSimulationScenarioMaxMagnitudes - // = static_cast>>( - // json.at("maxMagn")); - // const bool shouldRecompute = fullPatternSimulationScenarioMaxMagnitudes.size() - // != desiredBaseSimulationScenarioIndices.size(); - // if (!shouldRecompute) { - // return fullPatternSimulationScenarioMaxMagnitudes; - // } - // } -#ifdef POLYSCOPE_DEFINED - LinearSimulationModel linearSimulator; - SimulationResults fullPatternResults = linearSimulator.executeSimulation( - pFullPatternSimulationJob); + std::filesystem::path jobResultsDirectoryPath( + std::filesystem::path(intermediateResultsDirectoryPath) + .append("FullPatternResults") + .append(m_pFullPatternSimulationMesh->getLabel()) + .append(pFullPatternSimulationJob->getLabel())); + // .append(pFullPatternSimulationJob->getLabel() + ".json") + SimulationResults fullPatternResults; + if (std::filesystem::exists(jobResultsDirectoryPath)) { + fullPatternResults.load(std::filesystem::path(jobResultsDirectoryPath).append("Results"), + pFullPatternSimulationJob); + } else { + // const bool fullPatternScenarioMagnitudesExist = std::filesystem::exists( + // patternMaxForceMagnitudesFilePath); + // if (fullPatternScenarioMagnitudesExist) { + // nlohmann::json json; + // std::ifstream ifs(patternMaxForceMagnitudesFilePath.string()); + // ifs >> json; + // fullPatternSimulationScenarioMaxMagnitudes + // = static_cast>>( + // json.at("maxMagn")); + // const bool shouldRecompute = fullPatternSimulationScenarioMaxMagnitudes.size() + // != desiredBaseSimulationScenarioIndices.size(); + // if (!shouldRecompute) { + // return fullPatternSimulationScenarioMaxMagnitudes; + // } + // } + //#ifdef POLYSCOPE_DEFINED + // LinearSimulationModel linearSimulator; + // SimulationResults fullPatternResults = linearSimulator.executeSimulation( + // pFullPatternSimulationJob); -#else - SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, - simulationSettings); -#endif + //#else + fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, + simulationSettings); + fullPatternResults.save(jobResultsDirectoryPath); + } +//#endif // if (!fullPatternResults.converged) { // DRMSimulationModel::Settings simulationSettings_secondRound; // simulationSettings_secondRound.viscousDampingFactor = 2e-3; @@ -1816,9 +1957,9 @@ void ReducedModelOptimizer::optimize( std::filesystem::create_directories(outputPath); pFullPatternSimulationJob->save(outputPath); simulationSettings.save(outputPath); -#endif std::terminate(); return; +#endif // } } #ifdef POLYSCOPE_DEFINED diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 638ef1e..bb9ce98 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -5,7 +5,9 @@ #include "drmsimulationmodel.hpp" #include "edgemesh.hpp" #include "linearsimulationmodel.hpp" +#ifdef POLYSCOPE_DEFINED #include "matplot/matplot.h" +#endif #include "reducedmodeloptimizer_structs.hpp" #include "simulationmesh.hpp" #include @@ -64,8 +66,10 @@ public: enum OptimizationParameterIndex { E, A, I2, I3, J, R, Theta, NumberOfOptimizationParameters }; - constexpr static std::array simulationScenariosResolution = {11, 11, 20, 20, 20}; - constexpr static std::array baseScenarioWeights = {1, 1, 5, 5, 5}; + constexpr static std::array + simulationScenariosResolution = {11, 11, 20, 20, 20}; + constexpr static std::array + baseScenarioWeights = {1, 1, 100, 50, 50}; // constexpr static std::array simulationScenariosResolution = {3, 3, 3, 3, 3}; inline static int totalNumberOfSimulationScenarios = std::accumulate(simulationScenariosResolution.begin(), @@ -149,14 +153,13 @@ public: double innerHexagonSize, double innerHexagonRotationAngle); static double computeRawRotationalError( - const std::vector> &rotatedQuaternion_fullPattern, - const std::vector> &rotatedQuaternion_reducedPattern, + const std::vector &rotatedQuaternion_fullPattern, + const std::vector &rotatedQuaternion_reducedPattern, const std::unordered_map &reducedToFullInterfaceViMap); - static double computeRotationalError( - const std::vector> &rotatedQuaternion_fullPattern, - const std::vector> &rotatedQuaternion_reducedPattern, + static double computeRotationalError(const std::vector &rotatedQuaternion_fullPattern, + const std::vector &rotatedQuaternion_reducedPattern, const std::unordered_map &reducedToFullInterfaceViMap, const double &normalizationFactor); @@ -166,7 +169,8 @@ public: const std::unordered_map &reducedToFullInterfaceViMap, const double &normalizationFactor_translationalDisplacement, - const double &normalizationFactor_rotationalDisplacement); + const double &normalizationFactor_rotationalDisplacement, + const double &scenarioWeight); static void constructAxialSimulationScenario( const double &forceMagnitude, const std::vector> @@ -228,10 +232,6 @@ private: Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); static void runOptimization(const ReducedPatternOptimization::Settings &settings, ReducedPatternOptimization::Results &results); - std::vector> createFullPatternSimulationJobs( - const std::shared_ptr &pMesh, - const std::vector> - &maxForceMagnitudes); void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern); void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); static void initializeOptimizationParameters( @@ -245,20 +245,27 @@ private: const ReducedPatternOptimization::Settings &settings, ReducedPatternOptimization::Results &results); double computeFullPatternMaxSimulationForce( - const ReducedPatternOptimization::BaseSimulationScenario &scenario); - std::vector> - computeFullPatternMaxSimulationForces( - const std::vector - &desiredBaseSimulationScenarioIndices); - - std::vector> - getFullPatternMaxSimulationForces( - const std::vector - &desiredBaseSimulationScenarioIndices); + const ReducedPatternOptimization::BaseSimulationScenario &scenario) const; #ifdef DLIB_DEFINED static double objective(const dlib::matrix &x); #endif + std::array + getFullPatternMaxSimulationForces( + const std::vector + &desiredBaseSimulationScenarioIndices, + const std::filesystem::path &intermediateResultsDirectoryPath); + void computeScenarioWeights(const std::vector + &baseSimulationScenarios, + const std::array &maxForceMagnitudePerBaseScenario); + std::array + computeFullPatternMaxSimulationForces( + const std::vector + &desiredBaseSimulationScenario) const; + std::vector> createFullPatternSimulationJobs( + const std::shared_ptr &pMesh, + const std::array + &baseScenarioMaxForceMagnitudes) const; }; inline std::function &pReducedPatternSimulationMesh)>