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

refactoring

This commit is contained in:
iasonmanolas 2021-12-21 21:24:45 +02:00
parent 535dd8a44d
commit 3eaef7a37c
5 changed files with 475 additions and 281 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
CMakeLists.txt
View File

@ -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)

12
main.py
View File

@ -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")

50
src/main.cpp
View File

@ -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<std::string> 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<double> 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<double> 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<double> x = matplot::linspace(0, y.size() - 1, y.size());
std::vector<double> markerSizes(y.size(), 5);
SimulationResultsReporter::createPlot("scenario index", "error", x, y, markerSizes, colors);
// optimizationResults.saveMeshFiles();
optimizationResults.draw();
#endif

599
src/reducedmodeloptimizer.cpp
View File

@ -4,15 +4,19 @@
#include "trianglepatterngeometry.hpp"
#include "trianglepattterntopology.hpp"
#include <armadillo>
#include <atomic>
#include <chrono>
#include <execution>
#include <experimental/numeric>
#include <functional>
#include <parallel/parallel.h>
#define USE_SCENARIO_WEIGHTS
using namespace ReducedPatternOptimization;
struct GlobalOptimizationVariables
{
// std::vector<std::vector<Vector6d>> fullPatternDisplacements;
std::vector<SimulationResults> fullPatternResults;
std::vector<double> translationalDisplacementNormalizationValues;
std::vector<double> rotationalDisplacementNormalizationValues;
@ -25,14 +29,18 @@ struct GlobalOptimizationVariables
std::vector<size_t> objectiveValueHistory_iteration;
std::vector<double> objectiveValueHistory;
std::vector<double> plotColors;
std::array<double, 7> parametersInitialValue;
std::array<double, 7> optimizationInitialValue;
std::array<double,
ReducedModelOptimizer::OptimizationParameterIndex::NumberOfOptimizationParameters>
parametersInitialValue;
std::array<double,
ReducedModelOptimizer::OptimizationParameterIndex::NumberOfOptimizationParameters>
optimizationInitialValue;
std::vector<int> simulationScenarioIndices;
double minY{DBL_MAX};
std::vector<double> minX;
int numOfSimulationCrashes{false};
int numberOfFunctionCalls{0};
int numberOfOptimizationParameters{5};
int numberOfOptimizationParameters{5};
ReducedPatternOptimization::Settings optimizationSettings;
vcg::Triangle3<double> baseTriangle;
//Variables for finding the full pattern simulation forces
@ -51,8 +59,7 @@ struct GlobalOptimizationVariables
functions_updateReducedPatternParameter;
std::vector<double> xMin;
std::vector<double> xMax;
std::vector<double> scenarioEqualizationWeight;
std::vector<double> scenarioUserWeights;
std::vector<double> scenarioWeights;
} global;
double ReducedModelOptimizer::computeDisplacementError(
@ -128,23 +135,35 @@ double ReducedModelOptimizer::computeError(
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&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<double> 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<double>::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<double>::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<double, 0, 1> &x)
double ReducedModelOptimizer::objective(const double &xValue)
{
return objective({xValue});
return objective(std::vector<double>({xValue}));
}
double ReducedModelOptimizer::objective(const std::vector<double> &x)
@ -212,13 +235,8 @@ double ReducedModelOptimizer::objective(const std::vector<double> &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<double> simulationErrorsPerScenario(global.simulationScenarioIndices.size());
LinearSimulationModel simulator;
simulator.setStructure(pReducedPatternSimulationMesh);
// simulator.initialize();
@ -238,73 +256,60 @@ double ReducedModelOptimizer::objective(const std::vector<double> &x)
SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob);
// std::string filename;
if (!reducedModelResults.converged) {
totalError += std::numeric_limits<double>::max();
simulationErrorsPerScenario[simulationScenarioIndex]
= std::numeric_limits<double>::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<SimulationJob>(
// 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<double> &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<SimulationMesh>(fullModel);
pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002});
pFullPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus);
pFullPatternSimulationMesh->reset();
// Reduced pattern
pReducedPatternSimulationMesh = std::make_shared<SimulationMesh>(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<SimulationMesh> &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<SimulationMesh> &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<SimulationMesh> &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<SimulationMesh> &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<double> 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<SimulationMesh> &pReducedPatternSimulationMesh
= global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh;
function_updateReducedPattern(std::vector<double>(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::pair<BaseSimulationScenario, double>>
std::array<double, NumberOfBaseSimulationScenarios>
ReducedModelOptimizer::computeFullPatternMaxSimulationForces(
const std::vector<BaseSimulationScenario> &desiredBaseSimulationScenario)
const std::vector<BaseSimulationScenario> &desiredBaseSimulationScenario) const
{
std::vector<std::pair<BaseSimulationScenario, double>> fullPatternMaxSimulationForces;
fullPatternMaxSimulationForces.reserve(desiredBaseSimulationScenario.size());
std::array<double, NumberOfBaseSimulationScenarios> 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::pair<BaseSimulationScenario, double>>
std::array<double, NumberOfBaseSimulationScenarios>
ReducedModelOptimizer::getFullPatternMaxSimulationForces(
const std::vector<BaseSimulationScenario> &desiredBaseSimulationScenarioIndices)
const std::vector<BaseSimulationScenario> &desiredBaseSimulationScenarioIndices,
const std::filesystem::path &intermediateResultsDirectoryPath)
{
std::vector<std::pair<BaseSimulationScenario, double>> fullPatternSimulationScenarioMaxMagnitudes;
std::array<double, NumberOfBaseSimulationScenarios> 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<std::vector<std::pair<BaseSimulationScenario, double>>>(
json.at("maxMagn"));
const bool shouldRecompute = fullPatternSimulationScenarioMaxMagnitudes.size()
!= desiredBaseSimulationScenarioIndices.size();
if (!shouldRecompute) {
return fullPatternSimulationScenarioMaxMagnitudes;
}
= static_cast<std::array<double, NumberOfBaseSimulationScenarios>>(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<double, 0, 1> &) = &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<std::vector<std::vector<OptimizationParameterIndex>>> scenarioParameters =
@ -1002,38 +1016,57 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings,
std::vector<std::vector<std::vector<OptimizationParameterIndex>>> 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<std::vector<OptimizationParameterIndex>> scenarioParameterGroups
= scenarioParameters[scenario];
const int totalNumberOfOptimizationParameters
= std::accumulate(scenarioParameterGroups.begin(),
scenarioParameterGroups.end(),
0,
[](const int &sum,
const std::vector<OptimizationParameterIndex> &parameterGroup) {
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<OptimizationParameterIndex> &parameterGroup : 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<double> &x,
std::shared_ptr<SimulationMesh> &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<double, 0, 1> 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<int>(
(static_cast<double>(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<double, 0, 1> &) = &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<double>(optimalResult.x.begin(), optimalResult.x.end()),
std::vector<double>(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<BaseSimulationScenario> &baseSimulationScenarios,
const std::array<double, 5> &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<double, NumberOfBaseSimulationScenarios, 1> baseScenario_maxMagnitudeWeights(
maxForceMagnitudePerBaseScenario.data());
for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) {
const double baseScenarioEqualizationWeight
= totalNumberOfSimulationScenarios
/ static_cast<double>(simulationScenariosResolution[baseScenario]);
baseScenario_equalizationWeights(baseScenario) = baseScenarioEqualizationWeight;
const double baseScenarioUserWeight = static_cast<double>(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<std::shared_ptr<SimulationJob>> ReducedModelOptimizer::createFullPatternSimulationJobs(
const std::shared_ptr<SimulationMesh> &pMesh,
const std::vector<std::pair<BaseSimulationScenario, double>> &scenarioMaxForceMagnitudePairs)
const std::array<double, NumberOfBaseSimulationScenarios> &baseScenarioMaxForceMagnitudes) const
{
std::vector<std::shared_ptr<SimulationJob>> scenarios;
scenarios.resize(totalNumberOfSimulationScenarios);
@ -1547,18 +1678,10 @@ std::vector<std::shared_ptr<SimulationJob>> 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<BaseSimulationScenario, double> 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<std::shared_ptr<SimulationJob>> ReducedModelOptimizer::createFullPat
/ (numberOfSimulationScenarios);
const int baseSimulationScenarioIndexOffset
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin() + scenario,
simulationScenariosResolution.begin() + baseScenario,
0);
const double baseScenarioEqualizationWeight = static_cast<double>(
numberOfSimulationScenarios)
/ totalNumberOfSimulationScenarios;
const double baseScenarioUserWeight = static_cast<double>(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<SimulationJob>(job);
global.scenarioEqualizationWeight[simulationScenarioIndex]
= baseScenarioEqualizationWeight;
global.scenarioUserWeights[simulationScenarioIndex] = baseScenarioUserWeight;
scenarios[scenarioIndex] = std::make_shared<SimulationJob>(job);
}
}
@ -1712,13 +1831,25 @@ void ReducedModelOptimizer::optimize(
global.optimizationSettings = optimizationSettings;
global.pFullPatternSimulationMesh = m_pFullPatternSimulationMesh;
std::vector<std::pair<BaseSimulationScenario, double>> 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<double, NumberOfBaseSimulationScenarios> 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<SimulationJob> &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<std::vector<std::pair<BaseSimulationScenario, double>>>(
// 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<std::vector<std::pair<BaseSimulationScenario, double>>>(
// 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

51
src/reducedmodeloptimizer.hpp
View File

@ -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 <Eigen/Dense>
@ -64,8 +66,10 @@ public:
enum OptimizationParameterIndex { E, A, I2, I3, J, R, Theta, NumberOfOptimizationParameters };
constexpr static std::array<int, 5> simulationScenariosResolution = {11, 11, 20, 20, 20};
constexpr static std::array<int, 5> baseScenarioWeights = {1, 1, 5, 5, 5};
constexpr static std::array<int, ReducedPatternOptimization::NumberOfBaseSimulationScenarios>
simulationScenariosResolution = {11, 11, 20, 20, 20};
constexpr static std::array<int, ReducedPatternOptimization::NumberOfBaseSimulationScenarios>
baseScenarioWeights = {1, 1, 100, 50, 50};
// constexpr static std::array<int, 5> 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<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
const std::vector<Eigen::Quaterniond> &rotatedQuaternion_fullPattern,
const std::vector<Eigen::Quaterniond> &rotatedQuaternion_reducedPattern,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
static double computeRotationalError(
const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
static double computeRotationalError(const std::vector<Eigen::Quaterniond> &rotatedQuaternion_fullPattern,
const std::vector<Eigen::Quaterniond> &rotatedQuaternion_reducedPattern,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const double &normalizationFactor);
@ -166,7 +169,8 @@ public:
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&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<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
@ -228,10 +232,6 @@ private:
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
static void runOptimization(const ReducedPatternOptimization::Settings &settings,
ReducedPatternOptimization::Results &results);
std::vector<std::shared_ptr<SimulationJob>> createFullPatternSimulationJobs(
const std::shared_ptr<SimulationMesh> &pMesh,
const std::vector<std::pair<ReducedPatternOptimization::BaseSimulationScenario, double>>
&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<std::pair<ReducedPatternOptimization::BaseSimulationScenario, double>>
computeFullPatternMaxSimulationForces(
const std::vector<ReducedPatternOptimization::BaseSimulationScenario>
&desiredBaseSimulationScenarioIndices);
std::vector<std::pair<ReducedPatternOptimization::BaseSimulationScenario, double>>
getFullPatternMaxSimulationForces(
const std::vector<ReducedPatternOptimization::BaseSimulationScenario>
&desiredBaseSimulationScenarioIndices);
const ReducedPatternOptimization::BaseSimulationScenario &scenario) const;
#ifdef DLIB_DEFINED
static double objective(const dlib::matrix<double, 0, 1> &x);
#endif
std::array<double, ReducedPatternOptimization::NumberOfBaseSimulationScenarios>
getFullPatternMaxSimulationForces(
const std::vector<ReducedPatternOptimization::BaseSimulationScenario>
&desiredBaseSimulationScenarioIndices,
const std::filesystem::path &intermediateResultsDirectoryPath);
void computeScenarioWeights(const std::vector<ReducedPatternOptimization::BaseSimulationScenario>
&baseSimulationScenarios,
const std::array<double, 5> &maxForceMagnitudePerBaseScenario);
std::array<double, ReducedPatternOptimization::NumberOfBaseSimulationScenarios>
computeFullPatternMaxSimulationForces(
const std::vector<ReducedPatternOptimization::BaseSimulationScenario>
&desiredBaseSimulationScenario) const;
std::vector<std::shared_ptr<SimulationJob>> createFullPatternSimulationJobs(
const std::shared_ptr<SimulationMesh> &pMesh,
const std::array<double, ReducedPatternOptimization::NumberOfBaseSimulationScenarios>
&baseScenarioMaxForceMagnitudes) const;
};
inline std::function<void(const double &newE,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>