diff --git a/src/main.cpp b/src/main.cpp index 3df92e3..e8f1828 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,4 +1,5 @@ #include "beamformfinder.hpp" +#include "csvfile.hpp" #include "edgemesh.hpp" #include "flatpattern.hpp" #include "polyscope/curve_network.h" @@ -14,42 +15,7 @@ #include #include -// void scale(const std::vector& numberOfNodesPerSlot, -// FlatPattern &pattern) { -// const double desiredSize = 0.0025; // center to boundary -// const size_t interfaceSlotIndex = 4; // bottom edge -// std::unordered_map nodeToSlot; -// std::unordered_map> slotToNode; -// FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, -// nodeToSlot); FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, -// slotToNode); assert(slotToNode.find(interfaceSlotIndex) != slotToNode.end() -// && -// slotToNode.find(interfaceSlotIndex)->second.size() == 1); -// // Assuming that in the bottom edge there is only one vertex which is also -// the -// // interface -// const size_t baseTriangleInterfaceVi = -// *(slotToNode.find(interfaceSlotIndex)->second.begin()); -// const double currentSize = -// (pattern.vert[baseTriangleInterfaceVi].cP() - pattern.vert[0].cP()) -// .Norm(); -// const double scaleFactor = desiredSize / currentSize; -// vcg::tri::UpdatePosition::Scale(full, scaleFactor); -//} - int main(int argc, char *argv[]) { - // ReducedModelOptimizer::runBeamOptimization(); - // FlatPattern pattern("/home/iason/Models/valid_6777.ply"); - // FlatPattern pattern("/home/iason/Models/simple_beam_paper_example.ply"); - // pattern.savePly("fannedValid.ply"); - - // SimulationJob job; - // job.load("../" - // "ProblematicSimulationJobs/17:16_4_1_2021/" - // "reduced_pattern_Single bar reduced model_simScenario.json"); - // FormFinder ff; - // ff.executeSimulation(std::make_shared(job), false, false, - // true); // Create reduced models // FormFinder::runUnitTests(); @@ -76,179 +42,106 @@ int main(int argc, char *argv[]) { &CWReducedModel, &CCWReducedModel}; ReducedModelOptimizer optimizer(numberOfNodesPerSlot); + for (double dimRationMax = 0.75; dimRationMax < 2; dimRationMax += 0.05) { + ReducedModelOptimizer::xRange beamWidth{"B", 0.5, 1.5}; + ReducedModelOptimizer::xRange beamDimensionsRatio{"bOverh", 0.7, + dimRationMax}; + ReducedModelOptimizer::xRange beamE{"E", 0.07, 1.5}; + std::string xRangesString = beamWidth.toString() + " " + + beamDimensionsRatio.toString() + " " + + beamE.toString(); + std::cout << xRangesString << std::endl; + ReducedModelOptimizer::Settings settings; + settings.xRanges = {beamWidth, beamDimensionsRatio, beamE}; - // FlatPattern pattern("/home/iason/Models/TestSet_validPatterns/59.ply"); - // optimizer.initialize(pattern, *reducedModels[1], {}); - // optimizer.setInitialGuess - // ({0.050000000000000003, 4.0001308096448325, 4.2377893536538567, - // 5.6122373437520334}); - // Eigen::VectorXd optimalParameters = optimizer.optimize(0); + std::filesystem::path thisOptimizationDirectory( + std::filesystem::path("../OptimizationResults").append(xRangesString)); + std::filesystem::create_directory(thisOptimizationDirectory); + csvfile thisOptimizationStatistics( + std::filesystem::path(thisOptimizationDirectory) + .append("statistics.csv") + .string(), + true); - std::string fullPatternsTestSetDirectory = - // "/home/iason/Models/TestSet_validPatterns"; - "/home/iason/Documents/PhD/Research/Approximating shapes with flat " - "patterns/Pattern_enumerator/Results/1v_0v_2e_1e_1c_6fan/3/Valid"; - for (const auto &entry : - filesystem::directory_iterator(fullPatternsTestSetDirectory)) { - const auto filepath = - // std::filesystem::path(fullPatternsTestSetDirectory).append("305.ply"); - entry.path(); - const auto filepathString = filepath.string(); - // Use only the base triangle version - const std::string tiledSuffix = "_tiled.ply"; - if (filepathString.compare(filepathString.size() - tiledSuffix.size(), - tiledSuffix.size(), tiledSuffix) == 0) { - continue; + double totalError = 0; + int totalNumberOfSimulationCrashes = 0; + std::vector errors; + std::string fullPatternsTestSetDirectory = "../TestSet"; + // "/home/iason/Documents/PhD/Research/Approximating shapes with flat " + // "patterns/Pattern_enumerator/Results/1v_0v_2e_1e_1c_6fan/3/Valid"; + for (const auto &entry : + filesystem::directory_iterator(fullPatternsTestSetDirectory)) { + const auto filepath = + // std::filesystem::path(fullPatternsTestSetDirectory).append("305.ply"); + entry.path(); + const auto filepathString = filepath.string(); + // Use only the base triangle version + const std::string tiledSuffix = "_tiled.ply"; + if (filepathString.compare(filepathString.size() - tiledSuffix.size(), + tiledSuffix.size(), tiledSuffix) == 0) { + continue; + } + // std::cout << "Full pattern:" << filepathString << std::endl; + FlatPattern pattern(filepathString); + pattern.setLabel(filepath.stem().string()); + pattern.scale(0.03); + // for (int reducedPatternIndex = 0; + // reducedPatternIndex < reducedModels.size(); + // reducedPatternIndex++) { + // FlatPattern *pReducedModel = + // reducedModels[reducedPatternIndex]; + std::unordered_set optimizationExcludedEi; + // if (pReducedModel != + // reducedModels[0]) { // assumes that the singleBar reduced model + // // is + // // the + // // first in the reducedModels vector + // optimizationExcludedEi.insert(0); + // } + // FlatPattern cp; + // cp.copy(*reducedModels[0]); + optimizer.initializePatterns(pattern, *reducedModels[0], + optimizationExcludedEi); + // optimizer.optimize({ReducedModelOptimizer::Axial}); + ReducedModelOptimizer::Results optimizationResults = + optimizer.optimize(settings); + errors.push_back(optimizationResults.objectiveValue); + SimulationResultsReporter::createPlot( + "", "Objective value", errors, + std::filesystem::path(thisOptimizationDirectory) + .append("ObjectiveValues.png") + .string()); + thisOptimizationStatistics << filepath.stem().string() + << optimizationResults.objectiveValue; + if (optimizationResults.numberOfSimulationCrashes == 0) { + thisOptimizationStatistics << "No crashes"; + } else { + thisOptimizationStatistics + << optimizationResults.numberOfSimulationCrashes; + } + thisOptimizationStatistics << endrow; + + totalError += optimizationResults.objectiveValue; + totalNumberOfSimulationCrashes += + optimizationResults.numberOfSimulationCrashes; + // } } - std::cout << "Full pattern:" << filepathString << std::endl; - FlatPattern pattern(filepathString); - pattern.setLabel(filepath.stem().string()); - pattern.scale(0.03); - // for (int reducedPatternIndex = 0; - // reducedPatternIndex < reducedModels.size(); - // reducedPatternIndex++) { - // FlatPattern *pReducedModel = - // reducedModels[reducedPatternIndex]; - std::unordered_set optimizationExcludedEi; - // if (pReducedModel != - // reducedModels[0]) { // assumes that the singleBar reduced model - // // is - // // the - // // first in the reducedModels vector - // optimizationExcludedEi.insert(0); - // } - // FlatPattern cp; - // cp.copy(*reducedModels[0]); - optimizer.initializePatterns(pattern, *reducedModels[0], - optimizationExcludedEi); - optimizer.optimize({ReducedModelOptimizer::SimulationScenario::Bending}); - // } + + csvfile statistics(std::filesystem::path("../OptimizationResults") + .append("statistics.csv") + .string(), + false); + for (const auto &range : settings.xRanges) { + statistics << range.min << range.max; + } + statistics << settings.maxSimulations; + if (totalNumberOfSimulationCrashes == 0) { + statistics << "No crashes"; + } else { + statistics << totalNumberOfSimulationCrashes; + } + + statistics << totalError << endrow; } - - // // Full model simulation - // std::unordered_map> - // fixedVertices; - // // fixedVertices[0] = std::unordered_set{0, 1, 2, 3, 4, 5}; - // fixedVertices[3] = std::unordered_set{0, 1, 2, 3, 4, 5}; - // // fixedVertices[7] = std::unordered_set{0, 1, 2}; - // std::unordered_map nodalForces{ - // {15, {0, 0, 2000, 0, 0, 0}}}; - // SimulationJob fullModelSimulationJob{ - // std::make_shared(pattern), fixedVertices, nodalForces, - // {}}; - // Simulator formFinder; - // SimulationResults fullModelResults = - // formFinder.executeSimulation(fullModelSimulationJob); - // fullModelResults.simulationLabel = "Full Model"; - // fullModelResults.draw(fullModelSimulationJob); - // fullModelSimulationJob.draw(); - // double stiffnessFactor = 1; - - // while (true) { - // Reduced model simulation - // SimulationJob reducedModelSimulationJob = - // optimizer.getReducedSimulationJob(fullModelSimulationJob); - // const std::vector stiffnessVector{ - // fullModelSimulationJob.mesh->elements[0].properties.A, - // stiffnessFactor * - // fullModelSimulationJob.mesh->elements[0].properties.J, stiffnessFactor - // * fullModelSimulationJob.mesh->elements[0].properties.I2, - // stiffnessFactor * - // fullModelSimulationJob.mesh->elements[0].properties.I3}; - // for (EdgeIndex ei = 0; ei < reducedModelSimulationJob.mesh->EN(); ei++) { - // BeamFormFinder::Element &e = - // reducedModelSimulationJob.mesh->elements[ei]; e.properties.A = - // 0.00035185827018667374; - // // stifnessVector[0]; - // e.properties.J = // stiffnessVector[1]; - // 7.709325104874406e-08; - // e.properties.I2 = -0.0000015661453308127776; - // // stiffnessVector[2]; - // e.properties.I3 = 3.7099813776947167e-07; // stiffnessVector[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; - // } - - // SimulationResults reducedModelResults = - // formFinder.executeSimulation(reducedModelSimulationJob, true); - // reducedModelResults.simulationLabel = - // "Reduced Model_" + std::to_string(stiffnessFactor); - // reducedModelResults.draw(reducedModelSimulationJob); - // SimulationResultsReporter resultsReporter; - // resultsReporter.reportResults( - // {reducedModelResults}, - // std::filesystem::current_path().append("Results"), - // "_" + std::to_string(stiffnessFactor)); - // if (reducedModelResults.history.numberOfSteps == - // BeamFormFinder::Simulator::maxDRMIterations) { - // break; - // } - // stiffnessFactor *= 1.5; - // } - - // // Beam - // // registerWorldAxes(); - // VCGEdgeMesh mesh; - // // mesh.loadFromPly("/home/iason/Models/simple_beam_model_10elem_1m.ply"); - // mesh.loadFromPly("/home/iason/Coding/Projects/Approximating shapes with - // flat " - // "patterns/RodModelOptimizationForPatterns/build/Debug/" - // "Fanned_Single bar reduced model.ply"); - // // vcg::Matrix44d R; - // // R.SetRotateDeg(45, {0, 0, 1}); - // // vcg::tri::UpdatePosition::Matrix(mesh, R); - // // mesh.updateEigenEdgeAndVertices(); - // FormFinder formFinder; - // formFinder.runUnitTests(); - // std::unordered_map> - // fixedVertices; fixedVertices[1] = std::unordered_set{2}; - // fixedVertices[2] = std::unordered_set{2}; - // fixedVertices[3] = std::unordered_set{1, 2}; - // fixedVertices[4] = std::unordered_set{2}; - // fixedVertices[5] = std::unordered_set{2}; - // fixedVertices[6] = std::unordered_set{0, 1, 2}; - // // Forced displacements - // std::unordered_map nodalForcedDisplacements; - // // nodalForcedDisplacements[10] = Eigen::Vector3d(0, 0.1, 0.1); - // std::unordered_map nodalForcedNormal; - - // CoordType v14 = (mesh.vert[4].cP() - mesh.vert[0].cP()).Normalize(); - // CoordType v25 = (mesh.vert[5].cP() - mesh.vert[0].cP()).Normalize(); - // CoordType v36 = (mesh.vert[6].cP() - mesh.vert[0].cP()).Normalize(); - // const double forceMagnitude = 0.4; - // std::unordered_map nodalForces{ - // // {4, Vector6d({0, 0, 0, v14[0], v14[1], 0}) * forceMagnitude}, - // // {1, Vector6d({0, 0, 0, -v14[0], -v14[1], 0}) * - // forceMagnitude}, - // // {5, Vector6d({0, 0, 0, v25[0], v25[1], 0}) * forceMagnitude}, - // // {2, Vector6d({0, 0, 0, -v25[0], -v25[1], 0}) * - // forceMagnitude}, - // // {6, Vector6d({0, 0, 0, v36[0], v36[1], 0}) * forceMagnitude}, - // {3, Vector6d({0, 0, 0, -v36[0], -v36[1], 0}) * forceMagnitude}}; - // // nodalForcedNormal[10] = v; - // // nodalForcedNormal[0] = -v; - // // fixedVertices[10] = {0, 1}; - - // SimulationJob beamSimulationJob{std::make_shared(mesh), - // fixedVertices, - // nodalForces, - // {}, - // {}}; - // beamSimulationJob.mesh->setBeamMaterial(0.3, 1); - // beamSimulationJob.mesh->setBeamCrossSection( - // RectangularBeamDimensions{0.002, 0.002}); - // beamSimulationJob.registerForDrawing(); - // registerWorldAxes(); - // polyscope::show(); - // SimulationResults beamSimulationResults = - // formFinder.executeSimulation(beamSimulationJob, true, true); - // beamSimulationResults.label = "Beam"; - // beamSimulationResults.registerForDrawing(beamSimulationJob); - // polyscope::show(); - return 0; } diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index bd8082b..dbc12a9 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -33,7 +33,8 @@ std::vector firstOptimizationRoundResults; int g_firstRoundIterationIndex{0}; double minY{std::numeric_limits::max()}; std::vector minX; -std::vector failedSimulationsXRatio; +std::vector> failedSimulationsXRatio; +int numOfSimulationCrashes{false}; // struct OptimizationCallback { // double operator()(const size_t &iterations, const Eigen::VectorXd &x, @@ -157,8 +158,11 @@ 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.setDimensions(RectangularBeamDimensions(g_initialParameters(0) * x[0], - g_initialParameters(1) * x[1])); + e.setDimensions(RectangularBeamDimensions( + g_initialParameters(0) * x[0], + g_initialParameters(0) * x[0] / (g_initialParameters(1) * x[1]))); + e.setMaterial(ElementMaterial(e.material.poissonsRatio, + g_initialParameters(2) * x[2])); // e.properties.A = g_initialParameters(0) * x[0]; // e.properties.J = g_initialParameters(1) * x[1]; // e.properties.I2 = g_initialParameters(2) * x[2]; @@ -197,8 +201,8 @@ void updateMesh(long n, const double *x) { } } -double ReducedModelOptimizer::objective(double b, double h) { - std::vector x{b, h}; +double ReducedModelOptimizer::objective(double b, double h, double E) { + std::vector x{b, h, E}; return ReducedModelOptimizer::objective(x.size(), x.data()); } @@ -209,12 +213,12 @@ double ReducedModelOptimizer::objective(double x0, double x1, double x2, } double ReducedModelOptimizer::objective(long n, const double *x) { - std::cout.precision(17); + // std::cout.precision(17); - for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { - std::cout << "x[" + std::to_string(parameterIndex) + "]=" - << x[parameterIndex] << std::endl; - } + // for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { + // std::cout << "x[" + std::to_string(parameterIndex) + "]=" + // << x[parameterIndex] << std::endl; + // } const Element &e = g_reducedPatternSimulationJob[0]->pMesh->elements[0]; // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " // << e.firstBendingConstFactor << " " << @@ -235,21 +239,35 @@ double ReducedModelOptimizer::objective(long n, const double *x) { g_reducedPatternSimulationJob[simulationScenarioIndex], simulationSettings); std::string filename; - if (!reducedModelResults.converged) { + if (!reducedModelResults.converged /*&& + g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]] + ->pMesh->elements[0] + .A > 1e-8 & + x[0] / x[1] < 60*/) { 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]); + // if (failedSimulationsXRatio.empty()) { + // failedSimulationsXRatio.resize(2); + // } + // failedSimulationsXRatio[0].push_back(std::log(x[0] / x[1])); + // failedSimulationsXRatio[1].push_back( + // std::log(g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]] + // ->pMesh->elements[0] + // .A)); - SimulationResultsReporter::createPlot("", "b/h", failedSimulationsXRatio); + // SimulationResultsReporter::createPlot( + // "log(b/h)", "log(A)", failedSimulationsXRatio[0], + // failedSimulationsXRatio[1], "ratioToAPlot.png"); // std::cout << "Failed simulation" << std::endl; // simulationSettings.shouldDraw = true; // simulationSettings.debugMessages = true; // simulator.executeSimulation( // g_reducedPatternSimulationJob[simulationScenarioIndex], // simulationSettings); + numOfSimulationCrashes++; } else { error += computeError( reducedModelResults, @@ -259,19 +277,23 @@ double ReducedModelOptimizer::objective(long n, const double *x) { "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"; + auto pMesh = + g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]->pMesh; + + for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { + out << "x[" + std::to_string(parameterIndex) + "]=" << x[parameterIndex] + << std::endl; + } + out << pMesh->elements[0].dimensions.toString() + "\n" + + pMesh->elements[0].material.toString() + " \nA=" + << pMesh->elements[0].A << " \nratio=" + << pMesh->elements[0].dimensions.b / pMesh->elements[0].dimensions.h + << " \naxialRig:" << pMesh->elements[0].rigidity.axial + << " \ntorsionalRig:" << pMesh->elements[0].rigidity.torsional + << " \nfirstBendingRig:" << pMesh->elements[0].rigidity.firstBending + << " \nsecondBendingRig:" << pMesh->elements[0].rigidity.secondBending + << " \nscenario:" + simulationScenarioStrings[simulationScenarioIndex] + + "\n\n"; out.close(); } if (error < minY) { @@ -289,7 +311,6 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // [](double &c) { c = 0.7; }); // } // gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory, 6, colors); - std::cout << std::endl; // SimulationResultsReporter::createPlot("Number of Steps", "Objective // value", // gObjectiveValueHistory); @@ -488,7 +509,7 @@ void ReducedModelOptimizer::initializePatterns( void ReducedModelOptimizer::initializeOptimizationParameters( const std::shared_ptr &mesh) { - const int numberOfOptimizationParameters = 2; + const int numberOfOptimizationParameters = 3; g_initialParameters.resize(g_optimizeInnerHexagonSize ? numberOfOptimizationParameters + 1 : numberOfOptimizationParameters); @@ -503,9 +524,11 @@ void ReducedModelOptimizer::initializeOptimizationParameters( // e.secondBendingConstFactor *= stiffnessFactor; // } const double initialB = std::sqrt(mesh->elements[0].A); - const double initialH = initialB; + const double initialRatio = 1; + ; g_initialParameters(0) = initialB; - g_initialParameters(1) = initialH; + g_initialParameters(1) = initialRatio; + g_initialParameters(2) = mesh->elements[0].material.youngsModulus; // g_initialParameters = // m_pReducedPatternSimulationMesh->elements[0].properties.E; // for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) { @@ -575,13 +598,14 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( } } -Eigen::VectorXd ReducedModelOptimizer::runOptimization( +ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization( + const Settings &settings, double (*pObjectiveFunction)(long, const double *)) { gObjectiveValueHistory.clear(); const size_t n = g_initialParameters.rows(); - assert(n == 2); + assert(n == 3); // g_optimizeInnerHexagonSize ? 5: 4; // const size_t npt = (n + 1) * (n + 2) / 2; // // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; @@ -603,8 +627,6 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization( // {0.0001, 2, 2.000000005047502, 1.3055270196964464}; // {initialGuess(0), initialGuess(1), initialGuess(2), // initialGuess(3)}; - const double xMin = 1e-2; - const double xMax = 1e1; // assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) { // return d >= xMax || d <= xMin; // })); @@ -623,25 +645,29 @@ Eigen::VectorXd 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); - const size_t maxFun = 150; + dlib::matrix xMin(settings.xRanges.size()); + dlib::matrix xMax(settings.xRanges.size()); + for (int i = 0; i < settings.xRanges.size(); i++) { + xMin(i) = settings.xRanges[i].min; + xMax(i) = settings.xRanges[i].max; + } - double (*objF)(double, double) = &objective; + double (*objF)(double, double, double) = &objective; auto start = std::chrono::system_clock::now(); dlib::function_evaluation result = dlib::find_min_global( - objF, {xMin, xMin}, {xMax, xMax}, dlib::max_function_calls(maxFun)); + objF, xMin, xMax, dlib::max_function_calls(settings.maxSimulations)); auto end = std::chrono::system_clock::now(); auto elapsed = std::chrono::duration_cast(end - start); - std::cout << "Finished optimization with dlib" << endl; - std::cout << "Solution x:" << endl; - std::cout << result.x << endl; - std::cout << "Solution y:" << endl; - std::cout << result.y << endl; - std::cout << minY << endl; - std::cout << "Time(sec):" << elapsed.count() << std::endl; - std::cout << "Max function evaluations:" << maxFun << std::endl; - std::cout << "[" + std::to_string(xMin) + "," + std::to_string(xMax) + "]" - << std::endl; - std::cout << "Initial guess:" << initialGuess << std::endl; + Results results{numOfSimulationCrashes, minX, minY}; + std::cout << "Finished optimizing." << endl; + // std::cout << "Solution x:" << endl; + // std::cout << result.x << endl; + std::cout << "Objective value:" << minY << endl; + // std::cout << result.y << endl; + // std::cout << minY << endl; + // std::cout << "Time(sec):" << elapsed.count() << std::endl; + // std::cout << "Max function evaluations:" << maxFun << std::endl; + // std::cout << "Initial guess:" << initialGuess << std::endl; // const size_t maxFun = 200; // bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(), // rhobeg, rhoend, maxFun, w.data()); @@ -665,15 +691,7 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization( // rhobeg, rhoend, maxFun, w.data()); // std::cout << "Finished second optimization round" << std::endl; - Eigen::VectorXd eigenX(x.size(), 1); - for (size_t xi = 0; xi < x.size(); xi++) { - eigenX(xi) = minX[xi]; - // eigenX(xi) = x[xi]; - // eigenX(xi) = result.x(xi); - } - // for (size_t xi = 0; xi < x.size(); xi++) { - // } - return eigenX; + return results; } void ReducedModelOptimizer::setInitialGuess(std::vector v) { @@ -917,7 +935,7 @@ void ReducedModelOptimizer::runBeamOptimization() { // .append(g_reducedPatternSimulationJob[0].mesh->getLabel() + // "_simScenario.json")); - runOptimization(&objective); + runOptimization({}, &objective); fullModelResults.registerForDrawing(); SimulationResults reducedModelResults = simulator.executeSimulation( @@ -938,7 +956,6 @@ void ReducedModelOptimizer::visualizeResults( &fullPatternSimulationJobs, const std::vector &simulationScenarios) { m_pFullPatternSimulationMesh->registerForDrawing(); - double error = 0; for (const int simulationScenarioIndex : simulationScenarios) { const std::shared_ptr &pFullPatternSimulationJob = fullPatternSimulationJobs[simulationScenarioIndex]; @@ -952,7 +969,7 @@ void ReducedModelOptimizer::visualizeResults( g_reducedPatternSimulationJob[simulationScenarioIndex]; SimulationResults reducedModelResults = simulator.executeSimulation(pReducedPatternSimulationJob); - error += computeError( + double error = computeError( reducedModelResults, g_optimalReducedModelDisplacements[simulationScenarioIndex]); std::cout << "Error of simulation scenario " @@ -975,7 +992,8 @@ void ReducedModelOptimizer::visualizeResults( } } -void ReducedModelOptimizer::optimize( +ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( + const Settings &xRanges, const std::vector &simulationScenarios) { g_simulationScenarioIndices = simulationScenarios; @@ -992,6 +1010,7 @@ void ReducedModelOptimizer::optimize( g_reducedPatternSimulationJob.resize(6); g_firstRoundIterationIndex = 0; minY = std::numeric_limits::max(); + numOfSimulationCrashes = 0; // polyscope::removeAllStructures(); FormFinder::Settings settings; @@ -1014,8 +1033,8 @@ void ReducedModelOptimizer::optimize( g_reducedPatternSimulationJob[simulationScenarioIndex] = std::make_shared(reducedPatternSimulationJob); } - - Eigen::VectorXd optimalParameters = runOptimization(&objective); - updateMesh(optimalParameters.rows(), optimalParameters.data()); + Results optResults = runOptimization(xRanges, &objective); + updateMesh(optResults.x.size(), optResults.x.data()); // visualizeResults(simulationJobs, g_simulationScenarioIndices); + return optResults; } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 9adf71b..1f7a761 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -30,10 +30,31 @@ public: Saddle, NumberOfSimulationScenarios }; + struct Results { + int numberOfSimulationCrashes{0}; + std::vector x; + double objectiveValue; + }; + struct xRange { + std::string label; + double min; + double max; + std::string toString() const { + return label + "=[" + std::to_string(min) + "," + std::to_string(max) + + "]"; + } + }; + + struct Settings { + std::vector xRanges; + int maxSimulations{300}; + }; + inline static const std::string simulationScenarioStrings[] = { "Axial", "Shear", "Bending", "Double", "Saddle"}; - void optimize(const std::vector &simulationScenarios = - std::vector()); + Results optimize(const Settings &xRanges, + const std::vector &simulationScenarios = + std::vector()); double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot); @@ -57,7 +78,7 @@ public: std::vector &x); static double objective(double x0, double x1, double x2, double x3); - static double objective(double b, double h); + static double objective(double b, double h, double E); private: void @@ -68,8 +89,9 @@ private: const SimulationResults &fullModelResults, const std::unordered_map &displacementsReducedToFullMap, Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); - static Eigen::VectorXd - runOptimization(double (*pObjectiveFunction)(long, const double *)); + static Results runOptimization(const Settings &settings, + double (*pObjectiveFunction)(long, + const double *)); std::vector> createScenarios(const std::shared_ptr &pMesh); void computeMaps(FlatPattern &fullModel, FlatPattern &reducedPattern,