From cfdbe60f370dcfbb31325c167da2dd56902cdc13 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Tue, 9 Feb 2021 20:55:44 +0200 Subject: [PATCH] Introduced CW reduced model --- src/main.cpp | 7 +-- src/reducedmodeloptimizer.cpp | 80 ++++++++++++++++++++--------------- 2 files changed, 51 insertions(+), 36 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 4344e9c..429c0f2 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -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); diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index be9ee97..b557a7f 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -21,10 +21,10 @@ struct GlobalOptimizationVariables { std::unordered_set g_reducedPatternExludedEdges; Eigen::VectorXd g_initialParameters; std::vector - g_simulationScenarioIndices; + simulationScenarioIndices; std::vector g_innerHexagonVectors{6, VectorType(0, 0, 0)}; double g_innerHexagonInitialPos{0}; - bool g_optimizeInnerHexagonSize{false}; + bool optimizeInnerHexagonSize{false}; std::vector firstOptimizationRoundResults; int g_firstRoundIterationIndex{0}; double minY{DBL_MAX}; @@ -32,6 +32,7 @@ struct GlobalOptimizationVariables { std::vector> 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 &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 x{x0, x1, x2, x3}; +double ReducedModelOptimizer::objective(double b, double h, double E, + double innerHexagonSize) { + std::vector 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 &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 x(n, initialGuess); - if (global.g_optimizeInnerHexagonSize) { - x[n - 1] = global.g_innerHexagonInitialPos; - } + //const size_t initialGuess = 1; + //std::vector 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 w(wSize); // const size_t maxFun = std::min(100.0 * (x.size() + 1), 1000.0); - dlib::matrix xMin(settings.xRanges.size()); - dlib::matrix xMax(settings.xRanges.size()); - for (int i = 0; i < settings.xRanges.size(); i++) { + dlib::matrix xMin(global.numberOfOptimizationParameters); + dlib::matrix 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(end - start); Results results{global.numOfSimulationCrashes, global.minX, global.minY}; @@ -1031,9 +1045,9 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( const std::vector &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 &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; }