diff --git a/src/main.cpp b/src/main.cpp index ec5c1ee..3df92e3 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -121,7 +121,7 @@ int main(int argc, char *argv[]) { // cp.copy(*reducedModels[0]); optimizer.initializePatterns(pattern, *reducedModels[0], optimizationExcludedEi); - optimizer.optimize(); + optimizer.optimize({ReducedModelOptimizer::SimulationScenario::Bending}); // } } diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index bac17b0..bd8082b 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -20,7 +20,7 @@ std::unordered_map g_reducedToFullViMap; matplot::line_handle gPlotHandle; std::vector gObjectiveValueHistory; -Eigen::Vector4d g_initialX; +Eigen::Vector2d g_initialX; std::unordered_set g_reducedPatternExludedEdges; // double g_initialParameters; Eigen::VectorXd g_initialParameters; @@ -33,6 +33,7 @@ std::vector firstOptimizationRoundResults; int g_firstRoundIterationIndex{0}; double minY{std::numeric_limits::max()}; std::vector minX; +std::vector failedSimulationsXRatio; // struct OptimizationCallback { // double operator()(const size_t &iterations, const Eigen::VectorXd &x, @@ -142,7 +143,7 @@ double ReducedModelOptimizer::computeError( void updateMesh(long n, const double *x) { std::shared_ptr &pReducedPatternSimulationMesh = - g_reducedPatternSimulationJob[0]->pMesh; + g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]->pMesh; // const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0]; // std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << " // " @@ -156,17 +157,19 @@ void updateMesh(long n, const double *x) { // e.properties.E = g_initialParameters * x[ei]; // e.properties.E = g_initialParameters(0) * x[0]; // e.properties.G = g_initialParameters(1) * x[1]; - e.properties.A = g_initialParameters(0) * x[0]; - e.properties.J = g_initialParameters(1) * x[1]; - e.properties.I2 = g_initialParameters(2) * x[2]; - e.properties.I3 = g_initialParameters(3) * x[3]; + e.setDimensions(RectangularBeamDimensions(g_initialParameters(0) * x[0], + g_initialParameters(1) * x[1])); + // e.properties.A = g_initialParameters(0) * x[0]; + // e.properties.J = g_initialParameters(1) * x[1]; + // e.properties.I2 = g_initialParameters(2) * x[2]; + // e.properties.I3 = g_initialParameters(3) * x[3]; // e.properties.G = e.properties.E / (2 * (1 + 0.3)); - e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength; - e.torsionConstFactor = e.properties.G * e.properties.J / e.initialLength; - e.firstBendingConstFactor = - 2 * e.properties.E * e.properties.I2 / e.initialLength; - e.secondBendingConstFactor = - 2 * e.properties.E * e.properties.I3 / e.initialLength; + // e.axialConstFactor = e.properties.E * e.properties.A / + // e.initialLength; e.torsionConstFactor = e.properties.G * + // e.properties.J / e.initialLength; e.firstBendingConstFactor = + // 2 * e.properties.E * e.properties.I2 / e.initialLength; + // e.secondBendingConstFactor = + // 2 * e.properties.E * e.properties.I3 / e.initialLength; } // std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << " @@ -194,10 +197,15 @@ void updateMesh(long n, const double *x) { } } +double ReducedModelOptimizer::objective(double b, double h) { + std::vector x{b, h}; + 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}; - return ReducedModelOptimizer::objective(4, x.data()); + return ReducedModelOptimizer::objective(x.size(), x.data()); } double ReducedModelOptimizer::objective(long n, const double *x) { @@ -213,22 +221,62 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // e.secondBendingConstFactor // << std::endl; updateMesh(n, x); - std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " - << e.firstBendingConstFactor << " " << e.secondBendingConstFactor - << std::endl; + // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " + // << e.firstBendingConstFactor << " " << + // e.secondBendingConstFactor + // << std::endl; // run simulations double error = 0; + FormFinder::Settings simulationSettings; + // simulationSettings.shouldDraw = true; for (const int simulationScenarioIndex : g_simulationScenarioIndices) { SimulationResults reducedModelResults = simulator.executeSimulation( - g_reducedPatternSimulationJob[simulationScenarioIndex]); - error += computeError( - reducedModelResults, - g_optimalReducedModelDisplacements[simulationScenarioIndex]); + g_reducedPatternSimulationJob[simulationScenarioIndex], + simulationSettings); + std::string filename; + if (!reducedModelResults.converged) { + std::cout << "Failed simulation" << std::endl; + error += std::numeric_limits::max(); + filename = "/home/iason/Coding/Projects/Approximating shapes with flat " + "patterns/RodModelOptimizationForPatterns/build/" + "ProblematicSimulationJobs/nonConv_dimensions.txt"; + failedSimulationsXRatio.push_back(x[0] / x[1]); + + SimulationResultsReporter::createPlot("", "b/h", failedSimulationsXRatio); + // std::cout << "Failed simulation" << std::endl; + // simulationSettings.shouldDraw = true; + // simulationSettings.debugMessages = true; + // simulator.executeSimulation( + // g_reducedPatternSimulationJob[simulationScenarioIndex], + // simulationSettings); + } else { + error += computeError( + reducedModelResults, + g_optimalReducedModelDisplacements[simulationScenarioIndex]); + filename = "/home/iason/Coding/Projects/Approximating shapes with flat " + "patterns/RodModelOptimizationForPatterns/build/" + "ProblematicSimulationJobs/conv_dimensions.txt"; + } + std::ofstream out(filename, std::ios_base::app); + out << g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]] + ->pMesh->elements[0] + .dimensions.toString() + + " ratio=" + + std::to_string( + g_reducedPatternSimulationJob[simulationScenarioIndex] + ->pMesh->elements[0] + .dimensions.b / + g_reducedPatternSimulationJob[simulationScenarioIndex] + ->pMesh->elements[0] + .dimensions.h) + + " scenario:" + + simulationScenarioStrings[simulationScenarioIndex] + "\n"; + out.close(); } if (error < minY) { minY = error; - minX = std::vector({x[0], x[1], x[2], x[3]}); + minX.assign(x, x + n); } // compute error and return it @@ -242,8 +290,9 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // } // gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory, 6, colors); std::cout << std::endl; - SimulationResultsReporter::createPlot("Number of Steps", "Objective value", - gObjectiveValueHistory); + // SimulationResultsReporter::createPlot("Number of Steps", "Objective + // value", + // gObjectiveValueHistory); return error; } @@ -383,11 +432,11 @@ void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel, std::make_shared(reducedModel); m_pReducedPatternSimulationMesh->setBeamCrossSection( CrossSectionType{0.002, 0.002}); - m_pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1); + m_pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); m_pFullPatternSimulationMesh = std::make_shared(fullModel); m_pFullPatternSimulationMesh->setBeamCrossSection( CrossSectionType{0.002, 0.002}); - m_pFullPatternSimulationMesh->setBeamMaterial(0.3, 1); + m_pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); } ReducedModelOptimizer::ReducedModelOptimizer( @@ -434,12 +483,15 @@ void ReducedModelOptimizer::initializePatterns( } computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges); createSimulationMeshes(copyFullPattern, copyReducedPattern); - initializeStiffnesses(m_pReducedPatternSimulationMesh); + initializeOptimizationParameters(m_pReducedPatternSimulationMesh); } -void ReducedModelOptimizer::initializeStiffnesses( +void ReducedModelOptimizer::initializeOptimizationParameters( const std::shared_ptr &mesh) { - g_initialParameters.resize(g_optimizeInnerHexagonSize ? 5 : 4); + const int numberOfOptimizationParameters = 2; + g_initialParameters.resize(g_optimizeInnerHexagonSize + ? numberOfOptimizationParameters + 1 + : numberOfOptimizationParameters); // Save save the beam stiffnesses // for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) { // Element &e = pReducedModelElementalMesh->elements[ei]; @@ -450,16 +502,20 @@ void ReducedModelOptimizer::initializeStiffnesses( // e.firstBendingConstFactor *= stiffnessFactor; // e.secondBendingConstFactor *= stiffnessFactor; // } + const double initialB = std::sqrt(mesh->elements[0].A); + const double initialH = initialB; + g_initialParameters(0) = initialB; + g_initialParameters(1) = initialH; // g_initialParameters = // m_pReducedPatternSimulationMesh->elements[0].properties.E; // for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) { // } // g_initialParameters(0) = mesh->elements[0].properties.E; // g_initialParameters(1) = mesh->elements[0].properties.G; - g_initialParameters(0) = mesh->elements[0].properties.A; - g_initialParameters(1) = mesh->elements[0].properties.J; - g_initialParameters(2) = mesh->elements[0].properties.I2; - g_initialParameters(3) = mesh->elements[0].properties.I3; + // g_initialParameters(0) = mesh->elements[0].properties.A; + // g_initialParameters(1) = mesh->elements[0].properties.J; + // g_initialParameters(2) = mesh->elements[0].properties.I2; + // g_initialParameters(3) = mesh->elements[0].properties.I3; // } } @@ -524,7 +580,9 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization( gObjectiveValueHistory.clear(); - const size_t n = g_optimizeInnerHexagonSize ? 5 : 4; + const size_t n = g_initialParameters.rows(); + assert(n == 2); + // g_optimizeInnerHexagonSize ? 5: 4; // const size_t npt = (n + 1) * (n + 2) / 2; // // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; // assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2); @@ -546,7 +604,7 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization( // {initialGuess(0), initialGuess(1), initialGuess(2), // initialGuess(3)}; const double xMin = 1e-2; - const double xMax = 1e2; + const double xMax = 1e1; // assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) { // return d >= xMax || d <= xMin; // })); @@ -567,11 +625,10 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization( // const size_t maxFun = std::min(100.0 * (x.size() + 1), 1000.0); const size_t maxFun = 150; - double (*objF)(double, double, double, double) = &objective; + double (*objF)(double, double) = &objective; auto start = std::chrono::system_clock::now(); dlib::function_evaluation result = dlib::find_min_global( - objF, {xMin, xMin, xMin, xMin}, {xMax, xMax, xMax, xMax}, - dlib::max_function_calls(maxFun)); + objF, {xMin, xMin}, {xMax, xMax}, dlib::max_function_calls(maxFun)); auto end = std::chrono::system_clock::now(); auto elapsed = std::chrono::duration_cast(end - start); std::cout << "Finished optimization with dlib" << endl; @@ -797,9 +854,9 @@ ReducedModelOptimizer::createScenarios( void ReducedModelOptimizer::runBeamOptimization() { // load beams VCGEdgeMesh fullBeam; - fullBeam.loadFromPly("/home/iason/Models/simple_beam_model_10elem_1m.ply"); + fullBeam.loadPly("/home/iason/Models/simple_beam_model_10elem_1m.ply"); VCGEdgeMesh reducedBeam; - reducedBeam.loadFromPly("/home/iason/Models/simple_beam_model_4elem_1m.ply"); + reducedBeam.loadPly("/home/iason/Models/simple_beam_model_4elem_1m.ply"); fullBeam.registerForDrawing(); reducedBeam.registerForDrawing(); // polyscope::show(); @@ -814,7 +871,7 @@ void ReducedModelOptimizer::runBeamOptimization() { // full model simuilation job auto pFullPatternSimulationMesh = std::make_shared(fullBeam); pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02}); - pFullPatternSimulationMesh->setBeamMaterial(0.3, 1); + pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); std::unordered_map> fixedVertices; fixedVertices[0] = ::unordered_set({0, 1, 2, 3, 4, 5}); std::unordered_map forces; @@ -839,7 +896,7 @@ void ReducedModelOptimizer::runBeamOptimization() { // reduced model simuilation job auto reducedSimulationMesh = std::make_shared(reducedBeam); reducedSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02}); - reducedSimulationMesh->setBeamMaterial(0.3, 1); + reducedSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); g_reducedPatternSimulationJob.resize(numberOfSimulationScenarios); SimulationJob reducedSimJob; computeReducedModelSimulationJob(*pFullModelSimulationJob, @@ -850,7 +907,7 @@ void ReducedModelOptimizer::runBeamOptimization() { make_shared( reducedSimulationMesh, fullBeamSimulationJobLabel, reducedSimJob.constrainedVertices, reducedSimJob.nodalExternalForces); - initializeStiffnesses(reducedSimulationMesh); + initializeOptimizationParameters(reducedSimulationMesh); // const std::string simulationJobsPath = "SimulationJobs"; // std::filesystem::create_directory(simulationJobsPath); @@ -937,11 +994,13 @@ void ReducedModelOptimizer::optimize( minY = std::numeric_limits::max(); // polyscope::removeAllStructures(); + FormFinder::Settings settings; + // settings.shouldDraw = true; for (int simulationScenarioIndex : g_simulationScenarioIndices) { const std::shared_ptr &pFullPatternSimulationJob = simulationJobs[simulationScenarioIndex]; SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob); + simulator.executeSimulation(pFullPatternSimulationJob, settings); g_optimalReducedModelDisplacements[simulationScenarioIndex].resize( m_pReducedPatternSimulationMesh->VN(), 3); computeDesiredReducedModelDisplacements( @@ -957,6 +1016,6 @@ void ReducedModelOptimizer::optimize( } Eigen::VectorXd optimalParameters = runOptimization(&objective); - updateMesh(4, optimalParameters.data()); - visualizeResults(simulationJobs, g_simulationScenarioIndices); + updateMesh(optimalParameters.rows(), optimalParameters.data()); + // visualizeResults(simulationJobs, g_simulationScenarioIndices); } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 78b0561..9adf71b 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -57,6 +57,7 @@ public: std::vector &x); static double objective(double x0, double x1, double x2, double x3); + static double objective(double b, double h); private: void @@ -76,7 +77,7 @@ private: void createSimulationMeshes(FlatPattern &fullModel, FlatPattern &reducedModel); static void - initializeStiffnesses(const std::shared_ptr &mesh); + initializeOptimizationParameters(const std::shared_ptr &mesh); static double computeError(const SimulationResults &reducedPatternResults,