Introduced CW reduced model

This commit is contained in:
iasonmanolas 2021-02-09 20:55:44 +02:00
parent a0466fd1cb
commit cfdbe60f37
2 changed files with 51 additions and 36 deletions

View File

@ -46,8 +46,9 @@ int main(int argc, char *argv[]) {
ReducedModelOptimizer::xRange beamWidth{"B", 0.5, 1.5};
ReducedModelOptimizer::xRange beamDimensionsRatio{"bOverh", 0.7, 1.3};
ReducedModelOptimizer::xRange beamE{"E", 0.1, 1.9};
ReducedModelOptimizer::xRange innerHexagonSize{"HexagonSize", 0.1, 0.9};
// Test set of full patterns
std::string fullPatternsTestSetDirectory = "../TestSet";
std::string fullPatternsTestSetDirectory = "TestSet";
if (!std::filesystem::exists(
std::filesystem::path(fullPatternsTestSetDirectory))) {
std::cerr << "Full pattern directory does not exist: "
@ -73,7 +74,7 @@ int main(int argc, char *argv[]) {
pFullPattern->setLabel(filepath.stem().string());
pFullPattern->scale(0.03);
FlatPattern *pReducedPattern = new FlatPattern();
pReducedPattern->copy(*reducedModels[0]);
pReducedPattern->copy(*reducedModels[1]);
patternPairs.push_back(std::make_pair(pFullPattern, pReducedPattern));
}
@ -86,7 +87,7 @@ int main(int argc, char *argv[]) {
beamDimensionsRatio.toString() + " " +
beamE.toString();
std::cout << xRangesString << std::endl;
settings.xRanges = {beamWidth, beamDimensionsRatio, beamE};
settings.xRanges = {beamWidth, beamDimensionsRatio, beamE,innerHexagonSize};
// std::filesystem::path thisOptimizationDirectory(
// std::filesystem::path("../OptimizationResults").append(xRangesString));
// std::filesystem::create_directories(thisOptimizationDirectory);

View File

@ -21,10 +21,10 @@ struct GlobalOptimizationVariables {
std::unordered_set<size_t> g_reducedPatternExludedEdges;
Eigen::VectorXd g_initialParameters;
std::vector<ReducedModelOptimizer::SimulationScenario>
g_simulationScenarioIndices;
simulationScenarioIndices;
std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)};
double g_innerHexagonInitialPos{0};
bool g_optimizeInnerHexagonSize{false};
bool optimizeInnerHexagonSize{false};
std::vector<SimulationResults> firstOptimizationRoundResults;
int g_firstRoundIterationIndex{0};
double minY{DBL_MAX};
@ -32,6 +32,7 @@ struct GlobalOptimizationVariables {
std::vector<std::vector<double>> failedSimulationsXRatio;
int numOfSimulationCrashes{false};
int numberOfFunctionCalls{0};
int numberOfOptimizationParameters{ 3 };
};
// static GlobalOptimizationVariables global;
@ -155,7 +156,7 @@ void updateMesh(long n, const double *x) {
auto &global = tls[omp_get_thread_num()];
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
global
.g_reducedPatternSimulationJob[global.g_simulationScenarioIndices[0]]
.g_reducedPatternSimulationJob[global.simulationScenarioIndices[0]]
->pMesh;
// const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0];
// std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << "
@ -199,7 +200,7 @@ void updateMesh(long n, const double *x) {
// e.secondBendingConstFactor
// << std::endl;
if (global.g_optimizeInnerHexagonSize) {
if (global.optimizeInnerHexagonSize) {
assert(pReducedPatternSimulationMesh->EN() == 12);
for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN();
vi += 2) {
@ -219,9 +220,9 @@ double ReducedModelOptimizer::objective(double b, double h, double E) {
return ReducedModelOptimizer::objective(x.size(), x.data());
}
double ReducedModelOptimizer::objective(double x0, double x1, double x2,
double x3) {
std::vector<double> x{x0, x1, x2, x3};
double ReducedModelOptimizer::objective(double b, double h, double E,
double innerHexagonSize) {
std::vector<double> x{b, h, E,innerHexagonSize};
return ReducedModelOptimizer::objective(x.size(), x.data());
}
@ -249,8 +250,8 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
double error = 0;
FormFinder simulator;
FormFinder::Settings simulationSettings;
// simulationSettings.shouldDraw = true;
for (const int simulationScenarioIndex : global.g_simulationScenarioIndices) {
simulationSettings.shouldDraw = true;
for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
SimulationResults reducedModelResults = simulator.executeSimulation(
global.g_reducedPatternSimulationJob[simulationScenarioIndex],
simulationSettings);
@ -296,7 +297,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
auto pMesh =
global
.g_reducedPatternSimulationJob[global
.g_simulationScenarioIndices[0]]
.simulationScenarioIndices[0]]
->pMesh;
for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
@ -319,10 +320,11 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
global.minY = error;
global.minX.assign(x, x + n);
}
if (++global.numberOfFunctionCalls % 50 == 0) {
global.numberOfFunctionCalls++;
//if (++global.numberOfFunctionCalls % 50 == 0) {
std::cout << "Number of function calls:" << global.numberOfFunctionCalls
<< std::endl;
}
//}
// compute error and return it
global.gObjectiveValueHistory.push_back(error);
@ -506,8 +508,8 @@ void ReducedModelOptimizer::initializePatterns(
copyFullPattern.copy(fullPattern);
copyReducedPattern.copy(reducedPattern);
auto &global = tls[omp_get_thread_num()];
global.g_optimizeInnerHexagonSize = copyReducedPattern.EN() == 2;
if (global.g_optimizeInnerHexagonSize) {
global.optimizeInnerHexagonSize = copyReducedPattern.EN() == 2;
if (global.optimizeInnerHexagonSize) {
const double h = copyReducedPattern.getBaseTriangleHeight();
double baseTriangle_bottomEdgeSize = 2 * h / tan(vcg::math::ToRad(60.0));
VectorType baseTriangle_leftBottomNode(-baseTriangle_bottomEdgeSize / 2, -h,
@ -537,10 +539,10 @@ void ReducedModelOptimizer::initializePatterns(
void ReducedModelOptimizer::initializeOptimizationParameters(
const std::shared_ptr<SimulationMesh> &mesh) {
auto &global = tls[omp_get_thread_num()];
const int numberOfOptimizationParameters = 3;
global.g_initialParameters.resize(global.g_optimizeInnerHexagonSize
? numberOfOptimizationParameters + 1
: numberOfOptimizationParameters);
global.numberOfOptimizationParameters = 3;
global.g_initialParameters.resize(global.optimizeInnerHexagonSize
? ++global.numberOfOptimizationParameters
: global.numberOfOptimizationParameters);
// Save save the beam stiffnesses
// for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) {
// Element &e = pReducedModelElementalMesh->elements[ei];
@ -557,6 +559,9 @@ void ReducedModelOptimizer::initializeOptimizationParameters(
global.g_initialParameters(0) = initialB;
global.g_initialParameters(1) = initialRatio;
global.g_initialParameters(2) = mesh->elements[0].material.youngsModulus;
if (global.optimizeInnerHexagonSize) {
global.g_initialParameters(3) =global.g_innerHexagonInitialPos;
}
// g_initialParameters =
// m_pReducedPatternSimulationMesh->elements[0].properties.E;
// for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) {
@ -643,11 +648,11 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(
// "conditions is not recommended.");
// Set initial guess of solution
const size_t initialGuess = 1;
std::vector<double> x(n, initialGuess);
if (global.g_optimizeInnerHexagonSize) {
x[n - 1] = global.g_innerHexagonInitialPos;
}
//const size_t initialGuess = 1;
//std::vector<double> x(n, initialGuess);
//if (global.optimizeInnerHexagonSize) {
// x[n - 1] = global.g_innerHexagonInitialPos;
//}
/*if (!initialGuess.empty()) {
x = g_optimizationInitialGuess;
}*/ // {0.10000000000000 001,
@ -674,19 +679,28 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(
// const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2;
// std::vector<double> w(wSize);
// const size_t maxFun = std::min(100.0 * (x.size() + 1), 1000.0);
dlib::matrix<double, 0, 1> xMin(settings.xRanges.size());
dlib::matrix<double, 0, 1> xMax(settings.xRanges.size());
for (int i = 0; i < settings.xRanges.size(); i++) {
dlib::matrix<double, 0, 1> xMin(global.numberOfOptimizationParameters);
dlib::matrix<double, 0, 1> xMax(global.numberOfOptimizationParameters);
for (int i = 0; i < global.numberOfOptimizationParameters; i++) {
xMin(i) = settings.xRanges[i].min;
xMax(i) = settings.xRanges[i].max;
}
global.numberOfFunctionCalls = 0;
double (*objF)(double, double, double) = &objective;
auto start = std::chrono::system_clock::now();
dlib::function_evaluation result = dlib::find_min_global(
dlib::function_evaluation result;
if (global.optimizeInnerHexagonSize) {
double (*objF)(double, double, double,double) = &objective;
result = dlib::find_min_global(
objF, xMin, xMax, dlib::max_function_calls(settings.maxSimulations),
std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy);
}
else {
double (*objF)(double, double, double) = &objective;
result = dlib::find_min_global(
objF, xMin, xMax, dlib::max_function_calls(settings.maxSimulations),
std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy);
}
auto end = std::chrono::system_clock::now();
auto elapsed = std::chrono::duration_cast<std::chrono::seconds>(end - start);
Results results{global.numOfSimulationCrashes, global.minX, global.minY};
@ -1031,9 +1045,9 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
const std::vector<SimulationScenario> &simulationScenarios) {
auto &global = tls[omp_get_thread_num()];
global.g_simulationScenarioIndices = simulationScenarios;
if (global.g_simulationScenarioIndices.empty()) {
global.g_simulationScenarioIndices = {
global.simulationScenarioIndices = simulationScenarios;
if (global.simulationScenarioIndices.empty()) {
global.simulationScenarioIndices = {
SimulationScenario::Axial, SimulationScenario::Shear,
SimulationScenario::Bending, SimulationScenario::Dome,
SimulationScenario::Saddle};
@ -1050,7 +1064,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
FormFinder::Settings settings;
// settings.shouldDraw = true;
for (int simulationScenarioIndex : global.g_simulationScenarioIndices) {
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
simulationJobs[simulationScenarioIndex];
SimulationResults fullModelResults =
@ -1070,6 +1084,6 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
}
Results optResults = runOptimization(xRanges, &objective);
updateMesh(optResults.x.size(), optResults.x.data());
// visualizeResults(simulationJobs, g_simulationScenarioIndices);
visualizeResults(simulationJobs, global.simulationScenarioIndices);
return optResults;
}