From 81ac6c20b8b3249a1fa21f3b6fa88ccce86458cc Mon Sep 17 00:00:00 2001 From: Iason Date: Mon, 21 Dec 2020 17:56:21 +0200 Subject: [PATCH] Optimizes the reduced model for a set of simulation scenarios instead of a single scenario --- .gitignore | 4 + CMakeLists.txt | 1 - src/main.cpp | 11 +- src/main.cpp.autosave | 251 ---------------------------------- src/reducedmodeloptimizer.cpp | 247 ++++++++++++++++++++------------- src/reducedmodeloptimizer.hpp | 27 +++- 6 files changed, 181 insertions(+), 360 deletions(-) create mode 100644 .gitignore delete mode 100644 src/main.cpp.autosave diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..abffbab --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +/build +/TestSet +*.user + diff --git a/CMakeLists.txt b/CMakeLists.txt index 0c80a90..dcb4f0b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -75,7 +75,6 @@ target_include_directories(${PROJECT_NAME} PRIVATE ${BOBYQA_SOURCE_DIR}/include/ ) -#Use C++17 if(MSVC) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /std:c++20") else(MSVC) diff --git a/src/main.cpp b/src/main.cpp index 13fa179..7e427d9 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -93,7 +93,6 @@ int main(int argc, char *argv[]) { // "/home/iason/Models/TestSet_validPatterns/59.ply"); entry.path(); const auto filepathString = filepath.string(); - std::cout << "Full pattern:" << filepathString << std::endl; // Use only the base triangle version const std::string tiledSuffix = "_tiled.ply"; if (filepathString.compare(filepathString.size() - tiledSuffix.size(), @@ -102,11 +101,11 @@ int main(int argc, char *argv[]) { } FlatPattern pattern(filepathString); pattern.setLabel(filepath.stem().string()); - std::cout << "Testing Pattern:" << filepathString << std::endl; // for (int reducedPatternIndex = 0; // reducedPatternIndex < reducedModels.size(); // reducedPatternIndex++) { - // FlatPattern *pReducedModel = reducedModels[reducedPatternIndex]; + // FlatPattern *pReducedModel = + // reducedModels[reducedPatternIndex]; std::unordered_set optimizationExcludedEi; // if (pReducedModel != // reducedModels[0]) { // assumes that the singleBar reduced model @@ -116,9 +115,9 @@ int main(int argc, char *argv[]) { // optimizationExcludedEi.insert(0); // } std::cout << "Full pattern:" << filepathString << std::endl; - std::cout << "Reduced pattern:" << 1 << std::endl; - optimizer.initialize(pattern, *reducedModels[1], optimizationExcludedEi); - Eigen::VectorXd optimalParameters = optimizer.optimize(); + optimizer.initializePatterns(*reducedModels[2], *reducedModels[0], + optimizationExcludedEi); + optimizer.optimize(); // } } diff --git a/src/main.cpp.autosave b/src/main.cpp.autosave deleted file mode 100644 index 13fa179..0000000 --- a/src/main.cpp.autosave +++ /dev/null @@ -1,251 +0,0 @@ -#include "beamformfinder.hpp" -#include "edgemesh.hpp" -#include "flatpattern.hpp" -#include "polyscope/curve_network.h" -#include "polyscope/point_cloud.h" -#include "polyscope/polyscope.h" -#include "reducedmodeloptimizer.hpp" -#include "simulationhistoryplotter.hpp" -#include "trianglepattterntopology.hpp" -#include -#include -#include -#include -#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[]) { - // FlatPattern pattern("/home/iason/Models/valid_6777.ply"); - // FlatPattern pattern("/home/iason/Models/simple_beam_paper_example.ply"); - // pattern.savePly("fannedValid.ply"); - - // Create reduced models - // FormFinder::runUnitTests(); - const std::vector numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; - std::vector singleBarReducedModelEdges{vcg::Point2i(0, 3)}; - FlatPattern singleBarReducedModel(numberOfNodesPerSlot, - singleBarReducedModelEdges); - singleBarReducedModel.setLabel("Single bar reduced model"); - - std::vector CCWreducedModelEdges{vcg::Point2i(1, 5), - vcg::Point2i(3, 5)}; - FlatPattern CWReducedModel(numberOfNodesPerSlot, CCWreducedModelEdges); - CWReducedModel.setLabel("CW reduced model"); - - std::vector CWreducedModelEdges{vcg::Point2i(1, 5), - vcg::Point2i(3, 1)}; - FlatPattern CCWReducedModel(numberOfNodesPerSlot, CWreducedModelEdges); - CCWReducedModel.setLabel("CCW reduced model"); - - std::vector reducedModels{&singleBarReducedModel, - &CWReducedModel, &CCWReducedModel}; - - ReducedModelOptimizer optimizer(numberOfNodesPerSlot); - - // 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::string fullPatternsTestSetDirectory = - // "/home/iason/Models/TestSet_validPatterns"; - "/home/iason/Documents/PhD/Research/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("/home/iason/Models/valid_6777.ply"); - // std::filesystem::path( - // "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" - // "1v_0v_2e_1e_1c_6fan/3/Valid/222.ply"); - // std::filesystem::path( - // "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" - // "1v_0v_2e_1e_1c_6fan/2/Valid/33.ply"); - // std::filesystem::path( - // "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" - // "1v_0v_2e_1e_1c_6fan/3/Valid/431.ply"); - // std::filesystem::path( - // "/home/iason/Models/TestSet_validPatterns/59.ply"); - entry.path(); - const auto filepathString = filepath.string(); - std::cout << "Full pattern:" << filepathString << std::endl; - // Use only the base triangle version - const std::string tiledSuffix = "_tiled.ply"; - if (filepathString.compare(filepathString.size() - tiledSuffix.size(), - tiledSuffix.size(), tiledSuffix) == 0) { - continue; - } - FlatPattern pattern(filepathString); - pattern.setLabel(filepath.stem().string()); - std::cout << "Testing Pattern:" << filepathString << std::endl; - // 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); - // } - std::cout << "Full pattern:" << filepathString << std::endl; - std::cout << "Reduced pattern:" << 1 << std::endl; - optimizer.initialize(pattern, *reducedModels[1], optimizationExcludedEi); - Eigen::VectorXd optimalParameters = optimizer.optimize(); - // } - } - - // // 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 8cb18f5..14f5bf4 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -8,8 +8,8 @@ const bool gShouldDraw = true; FormFinder simulator; -Eigen::MatrixX3d g_optimalReducedModelDisplacements; -SimulationJob g_reducedPatternSimulationJob; +std::vector g_optimalReducedModelDisplacements; +std::vector g_reducedPatternSimulationJob; std::unordered_map g_reducedToFullInterfaceViMap; matplot::line_handle gPlotHandle; @@ -18,6 +18,7 @@ Eigen::Vector4d g_initialX; std::unordered_set g_reducedPatternExludedEdges; // double g_initialParameters; Eigen::VectorXd g_initialParameters; +ReducedModelOptimizer::SimulationScenario g_chosenSimulationScenarioName; // struct OptimizationCallback { // double operator()(const size_t &iterations, const Eigen::VectorXd &x, @@ -96,16 +97,35 @@ Eigen::VectorXd g_initialParameters; // } //}; -double ReducedModelOptimizer::objective(long n, const double *x) { - std::cout.precision(17); - - for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { - std::cout << "x[" + std::to_string(parameterIndex) + "]=" - << x[parameterIndex] << std::endl; +double ReducedModelOptimizer::computeError( + const SimulationResults &reducedPatternResults, + const Eigen::MatrixX3d &optimalReducedPatternDisplacements) { + double error = 0; + for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) { + VertexIndex reducedModelVi = reducedFullViPair.first; + // const auto pos = + // g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP(); + // std::cout << "Interface vi " << reducedModelVi << " is at position " + // << pos[0] << " " << pos[1] << " " << pos[2] << std::endl; + Eigen::Vector3d vertexDisplacement( + reducedPatternResults.displacements[reducedModelVi][0], + reducedPatternResults.displacements[reducedModelVi][1], + reducedPatternResults.displacements[reducedModelVi][2]); + Eigen::Vector3d errorVector = + Eigen::Vector3d( + optimalReducedPatternDisplacements.row(reducedModelVi)) - + vertexDisplacement; + // error += errorVector.squaredNorm(); + error += errorVector.norm(); } + return error; +} - for (EdgeIndex ei = 0; ei < g_reducedPatternSimulationJob.mesh->EN(); ei++) { - Element &e = g_reducedPatternSimulationJob.mesh->elements[ei]; +void updateMesh(long n, const double *x) { + std::shared_ptr pReducedPatternSimulationMesh = + g_reducedPatternSimulationJob[0].mesh; + for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { + Element &e = pReducedPatternSimulationMesh->elements[ei]; // if (g_reducedPatternExludedEdges.contains(ei)) { // continue; // } @@ -124,30 +144,56 @@ double ReducedModelOptimizer::objective(long n, const double *x) { e.secondBendingConstFactor = 2 * e.properties.E * e.properties.I3 / e.initialLength; } - // run simulation - SimulationResults reducedModelResults = - simulator.executeSimulation(g_reducedPatternSimulationJob, false, false); +} + +double ReducedModelOptimizer::objectiveAllScenarios(long n, const double *x) { + std::cout.precision(17); + + for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { + std::cout << "x[" + std::to_string(parameterIndex) + "]=" + << x[parameterIndex] << std::endl; + } + updateMesh(n, x); + + // run simulations + double error = 0; + for (int simulationScenarioIndex = SimulationScenario::Axial; + simulationScenarioIndex != + SimulationScenario::NumberOfSimulationScenarios; + simulationScenarioIndex++) { + SimulationResults reducedModelResults = simulator.executeSimulation( + g_reducedPatternSimulationJob[simulationScenarioIndex], false, false); + error += computeError( + reducedModelResults, + g_optimalReducedModelDisplacements[simulationScenarioIndex]); + } // compute error and return it - double error = 0; - for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) { - VertexIndex reducedModelVi = reducedFullViPair.first; - // const auto pos = - // g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP(); - // std::cout << "Interface vi " << reducedModelVi << " is at position " - // << pos[0] << " " << pos[1] << " " << pos[2] << std::endl; - Eigen::Vector3d vertexDisplacement( - reducedModelResults.displacements[reducedModelVi][0], - reducedModelResults.displacements[reducedModelVi][1], - reducedModelResults.displacements[reducedModelVi][2]); - Eigen::Vector3d errorVector = - Eigen::Vector3d( - g_optimalReducedModelDisplacements.row(reducedModelVi)) - - vertexDisplacement; - // error += errorVector.squaredNorm(); - error += errorVector.norm(); + gObjectiveValueHistory.push_back(error); + auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), + gObjectiveValueHistory.size()); + gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory); + + return error; +} + +double ReducedModelOptimizer::objectiveSingleScenario(long n, const double *x) { + std::cout.precision(17); + + for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { + std::cout << "x[" + std::to_string(parameterIndex) + "]=" + << x[parameterIndex] << std::endl; } + updateMesh(n, x); + + SimulationResults reducedModelResults = simulator.executeSimulation( + g_reducedPatternSimulationJob[g_chosenSimulationScenarioName], false, + false); + const double error = computeError( + reducedModelResults, + g_optimalReducedModelDisplacements[g_chosenSimulationScenarioName]); + gObjectiveValueHistory.push_back(error); auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), gObjectiveValueHistory.size()); @@ -302,7 +348,7 @@ ReducedModelOptimizer::ReducedModelOptimizer( FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode); } -void ReducedModelOptimizer::initialize( +void ReducedModelOptimizer::initializePatterns( FlatPattern &fullPattern, FlatPattern &reducedPattern, const std::unordered_set &reducedModelExcludedEdges) { assert(fullPattern.VN() == reducedPattern.VN() && @@ -316,6 +362,8 @@ void ReducedModelOptimizer::initialize( computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges); createSimulationMeshes(copyFullPattern, copyReducedPattern); initializeStiffnesses(); + copyFullPattern.registerForDrawing(); + polyscope::show(); } void ReducedModelOptimizer::initializeStiffnesses() { @@ -393,7 +441,6 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) { - optimalDisplacementsOfReducedModel.resize(0, 3); optimalDisplacementsOfReducedModel.resize( m_pReducedPatternSimulationMesh->VN(), 3); optimalDisplacementsOfReducedModel.setZero( @@ -411,28 +458,11 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( } } -Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( - const SimulationJob &fullPatternSimulationJob) { +Eigen::VectorXd ReducedModelOptimizer::runOptimization( + double (*pObjectiveFunction)(long, const double *)) { - // fullPatternSimulationJob.registerForDrawing(); - // polyscope::show(); gObjectiveValueHistory.clear(); - fullPatternSimulationJob.mesh->savePly( - "Fanned_" + m_pFullModelSimulationMesh->getLabel() + ".ply"); - SimulationResults fullModelResults = - simulator.executeSimulation(fullPatternSimulationJob, false); - fullModelResults.label = "fullModel"; - fullModelResults.registerForDrawing(fullPatternSimulationJob); - registerWorldAxes(); - polyscope::show(); - computeDesiredReducedModelDisplacements(fullModelResults, - g_optimalReducedModelDisplacements); - computeReducedModelSimulationJob(fullPatternSimulationJob, - g_reducedPatternSimulationJob); - // gReducedPatternSimulationJob.registerForDrawing(); - // polyscope::show(); - double (*pObjectiveFunction)(long, const double *) = &objective; // const size_t n = m_pReducedPatternSimulationMesh->EN(); const size_t n = 4; const size_t npt = 2 * n; @@ -442,7 +472,7 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( "conditions is not recommended."); // Set initial guess of solution - std::vector x(n, 4); + std::vector x(n, 2); if (!initialGuess.empty()) { x = initialGuess; } // {0.10000000000000 001, 2, 1.9999999971613847, 6.9560343643347764}; @@ -450,8 +480,8 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( // {0.0001, 2, 2.000000005047502, 1.3055270196964464}; // {initialGuess(0), initialGuess(1), initialGuess(2), // initialGuess(3)}; - const double xMin = 0.12; - const double xMax = 10; + const double xMin = 1e-2; + const double xMax = 5000; // assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) { // return d >= xMax || d <= xMin; // })); @@ -470,29 +500,6 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(), rhobeg, rhoend, maxFun, w.data()); - SimulationResults reducedModelOptimizedResults = - simulator.executeSimulation(g_reducedPatternSimulationJob); - - double error = 0; - for (const auto reducedToFullViPair : g_reducedToFullInterfaceViMap) { - const size_t reducedInterfaceVi = reducedToFullViPair.first; - error += - (Eigen::Vector3d( - g_optimalReducedModelDisplacements.row(reducedInterfaceVi)) - - Eigen::Vector3d( - reducedModelOptimizedResults.displacements[reducedInterfaceVi][0], - reducedModelOptimizedResults.displacements[reducedInterfaceVi][1], - reducedModelOptimizedResults.displacements[reducedInterfaceVi][2])) - .norm(); - } - - std::cout << "Final objective value:" << error << std::endl; - - reducedModelOptimizedResults.label = "reducedModel"; - reducedModelOptimizedResults.registerForDrawing( - g_reducedPatternSimulationJob); - polyscope::show(); - Eigen::VectorXd eigenX(x.size(), 1); for (size_t xi = 0; xi < x.size(); xi++) { eigenX(xi) = x[xi]; @@ -659,29 +666,77 @@ std::vector ReducedModelOptimizer::createScenarios( return scenarios; } -Eigen::VectorXd ReducedModelOptimizer::optimize(const int &simulationScenario) { +void ReducedModelOptimizer::optimize(const int &simulationScenario) { + std::vector simulationJobs = createScenarios(m_pFullModelSimulationMesh); - std::vector results; - if (simulationScenario == -1) { // run all scenarios - for (int i = 0; i < simulationJobs.size(); i++) { - std::cout << "Running simulation job:" << i << std::endl; - const SimulationJob &job = simulationJobs[i]; - polyscope::removeAllStructures(); - auto result = optimizeForSimulationJob(job); - results.push_back(result); - } + g_optimalReducedModelDisplacements.resize(6); + g_reducedPatternSimulationJob.resize(6); + polyscope::removeAllStructures(); + + for (int simulationScenarioIndex = SimulationScenario::Axial; + simulationScenarioIndex != + SimulationScenario::NumberOfSimulationScenarios; + simulationScenarioIndex++) { + const SimulationJob &fullPatternSimulationJob = + simulationJobs[simulationScenarioIndex]; + // fullPatternSimulationJob.mesh->savePly( + // "Fanned_" + m_pFullModelSimulationMesh->getLabel() + ".ply"); + SimulationResults fullModelResults = + simulator.executeSimulation(fullPatternSimulationJob, false); + // fullModelResults.label = + // "fullModel_" + SimulationScenarioStrings[simulationScenarioIndex]; + // fullModelResults.registerForDrawing(fullPatternSimulationJob); + computeDesiredReducedModelDisplacements( + fullModelResults, + g_optimalReducedModelDisplacements[simulationScenarioIndex]); + computeReducedModelSimulationJob( + fullPatternSimulationJob, + g_reducedPatternSimulationJob[simulationScenarioIndex]); + } + + if (simulationScenario == -1) { + double (*pObjectiveFunction)(long, const double *) = &objectiveAllScenarios; + runOptimization(pObjectiveFunction); } else { // run chosen - std::cout << "Running simulation job:" << simulationScenario << std::endl; - const SimulationJob &job = simulationJobs[simulationScenario]; - polyscope::removeAllStructures(); - auto result = optimizeForSimulationJob(job); - results.push_back(result); + g_chosenSimulationScenarioName = SimulationScenario(simulationScenario); + double (*pObjectiveFunction)(long, const double *) = + &objectiveSingleScenario; + runOptimization(pObjectiveFunction); } - if (results.empty()) { - return Eigen::VectorXd(); + m_pFullModelSimulationMesh->registerForDrawing(); + double error = 0; + for (int simulationScenarioIndex = SimulationScenario::Axial; + simulationScenarioIndex != + SimulationScenario::NumberOfSimulationScenarios; + simulationScenarioIndex++) { + const SimulationJob &fullPatternSimulationJob = + simulationJobs[simulationScenarioIndex]; + SimulationResults fullModelResults = + simulator.executeSimulation(fullPatternSimulationJob, false); + fullModelResults.label = + "fullModel_" + SimulationScenarioStrings[simulationScenarioIndex]; + fullModelResults.registerForDrawing(fullPatternSimulationJob); + const SimulationJob &reducedPatternSimulationJob = + g_reducedPatternSimulationJob[simulationScenarioIndex]; + SimulationResults reducedModelResults = + simulator.executeSimulation(reducedPatternSimulationJob, false, false); + error += computeError( + reducedModelResults, + g_optimalReducedModelDisplacements[simulationScenarioIndex]); + reducedModelResults.label = + "reducedModel_" + SimulationScenarioStrings[simulationScenarioIndex]; + reducedModelResults.registerForDrawing(reducedPatternSimulationJob); + // registerWorldAxes(); + std::cout << "A full:" + << m_pFullModelSimulationMesh->elements[0].properties.A + << std::endl; + std::cout << "A reduced:" + << m_pReducedPatternSimulationMesh->elements[0].properties.A + << std::endl; + polyscope::show(); + fullModelResults.unregister(); + reducedModelResults.unregister(); } - - return results[0]; } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 286ab8f..a687f6d 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -22,7 +22,17 @@ class ReducedModelOptimizer { std::vector initialGuess; public: - Eigen::VectorXd optimize(const int &simulationScenario = -1); + enum SimulationScenario { + Axial, + Shear, + Bending, + Double, + Saddle, + NumberOfSimulationScenarios + }; + inline static const std::string SimulationScenarioStrings[] = { + "Axial", "Shear", "Bending", "Double", "Saddle"}; + void optimize(const int &simulationScenario = -1); double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot); @@ -33,8 +43,9 @@ public: SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob); - void initialize(FlatPattern &fullPattern, FlatPattern &reducedPatterm, - const std::unordered_set &reducedModelExcludedEges); + void initializePatterns( + FlatPattern &fullPattern, FlatPattern &reducedPatterm, + const std::unordered_set &reducedModelExcludedEges); void setInitialGuess(std::vector v); @@ -42,8 +53,8 @@ private: void computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); - Eigen::VectorXd - optimizeForSimulationJob(const SimulationJob &fullModelSimulationJob); + Eigen::VectorXd runOptimization(double (*pObjectiveFunction)(long, + const double *)); std::vector createScenarios(const std::shared_ptr &pMesh); void computeMaps(FlatPattern &fullModel, FlatPattern &reducedPattern, @@ -52,7 +63,11 @@ private: FlatPattern &reducedModel); void initializeStiffnesses(); - static double objective(long n, const double *x); + static double + computeError(const SimulationResults &reducedPatternResults, + const Eigen::MatrixX3d &optimalReducedPatternDisplacements); + static double objectiveSingleScenario(long n, const double *x); + static double objectiveAllScenarios(long n, const double *x); }; #endif // REDUCEDMODELOPTIMIZER_HPP