From 42ccaba01ed77ca0720bc96a29d631cdd65e4749 Mon Sep 17 00:00:00 2001 From: Iason Date: Mon, 4 Jan 2021 14:12:25 +0200 Subject: [PATCH] Exporting simulation scenario when number of simulations is greater than some treshold. --- src/main.cpp | 37 ++- src/reducedmodeloptimizer.cpp | 467 ++++++++++++++++++++++------------ src/reducedmodeloptimizer.hpp | 35 ++- 3 files changed, 350 insertions(+), 189 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 631be58..f5ac455 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -38,10 +38,19 @@ //} 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/12:24_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(); const std::vector numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; @@ -86,15 +95,18 @@ int main(int argc, char *argv[]) { // 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(); + std::filesystem::path( + "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" + "1v_0v_2e_1e_1c_6fan/3/Valid/624.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(); // Use only the base triangle version const std::string tiledSuffix = "_tiled.ply"; @@ -119,9 +131,10 @@ int main(int argc, char *argv[]) { // optimizationExcludedEi.insert(0); // } std::cout << "Full pattern:" << filepathString << std::endl; - optimizer.initializePatterns(pattern, *reducedModels[1], - optimizationExcludedEi); - optimizer.optimize(); + FlatPattern cp; + cp.copy(*reducedModels[0]); + optimizer.initializePatterns(cp, *reducedModels[0], optimizationExcludedEi); + optimizer.optimize({ReducedModelOptimizer::SimulationScenario::Bending}); // } } diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 8b963c0..3495fe1 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -9,18 +9,25 @@ const bool gShouldDraw = true; FormFinder simulator; std::vector g_optimalReducedModelDisplacements; -std::vector g_reducedPatternSimulationJob; +std::vector g_fullPatternSimulationJob; +std::vector> g_reducedPatternSimulationJob; std::unordered_map - g_reducedToFullInterfaceViMap; + reducedToFullInterfaceViMap; +std::unordered_map + g_reducedToFullViMap; matplot::line_handle gPlotHandle; std::vector gObjectiveValueHistory; Eigen::Vector4d g_initialX; std::unordered_set g_reducedPatternExludedEdges; // double g_initialParameters; Eigen::VectorXd g_initialParameters; -ReducedModelOptimizer::SimulationScenario g_chosenSimulationScenarioName; +std::vector + g_simulationScenarioIndices; std::vector g_innerHexagonVectors{6, VectorType(0, 0, 0)}; double g_innerHexagonInitialPos = 0; +bool g_optimizeInnerHexagonSize{false}; +std::vector firstOptimizationRoundResults; +int g_firstRoundIterationIndex{0}; // struct OptimizationCallback { // double operator()(const size_t &iterations, const Eigen::VectorXd &x, @@ -103,7 +110,7 @@ double ReducedModelOptimizer::computeError( const SimulationResults &reducedPatternResults, const Eigen::MatrixX3d &optimalReducedPatternDisplacements) { double error = 0; - for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) { + for (const auto reducedFullViPair : g_reducedToFullViMap) { VertexIndex reducedModelVi = reducedFullViPair.first; // const auto pos = // g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP(); @@ -124,16 +131,21 @@ double ReducedModelOptimizer::computeError( } void updateMesh(long n, const double *x) { - std::shared_ptr pReducedPatternSimulationMesh = - g_reducedPatternSimulationJob[0].mesh; + std::shared_ptr &pReducedPatternSimulationMesh = + g_reducedPatternSimulationJob[0]->pMesh; + // const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0]; + // std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << " + // " + // << elem.firstBendingConstFactor << " " + // << elem.secondBendingConstFactor << std::endl; for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; // if (g_reducedPatternExludedEdges.contains(ei)) { // continue; // } // e.properties.E = g_initialParameters * x[ei]; - // e.properties.E = g_initialParameters(ei, 0) * x[ei]; - // e.properties.G = g_initialParameters(1) * x[1]; + // 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]; @@ -147,36 +159,51 @@ void updateMesh(long n, const double *x) { 2 * e.properties.E * e.properties.I3 / e.initialLength; } - if (n == 5) { + // std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << " + // " + // << elem.firstBendingConstFactor << " " + // << elem.secondBendingConstFactor << std::endl; + // const Element &e = pReducedPatternSimulationMesh->elements[0]; + // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " + // << e.firstBendingConstFactor << " " << + // e.secondBendingConstFactor + // << std::endl; + + if (g_optimizeInnerHexagonSize) { assert(pReducedPatternSimulationMesh->EN() == 12); for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN(); vi += 2) { pReducedPatternSimulationMesh->vert[vi].P() = - g_innerHexagonVectors[vi / 2] * x[4]; + g_innerHexagonVectors[vi / 2] * x[n - 1]; } - pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); pReducedPatternSimulationMesh->reset(); + pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); // pReducedPatternSimulationMesh->registerForDrawing("Optimized // hexagon"); polyscope::show(); } } -double ReducedModelOptimizer::objectiveAllScenarios(long n, const double *x) { +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; } + const Element &e = g_reducedPatternSimulationJob[0]->pMesh->elements[0]; + // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " + // << e.firstBendingConstFactor << " " << + // e.secondBendingConstFactor + // << std::endl; updateMesh(n, x); + std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " + << e.firstBendingConstFactor << " " << e.secondBendingConstFactor + << std::endl; // run simulations double error = 0; - for (int simulationScenarioIndex = SimulationScenario::Axial; - simulationScenarioIndex != - SimulationScenario::NumberOfSimulationScenarios; - simulationScenarioIndex++) { + for (const int simulationScenarioIndex : g_simulationScenarioIndices) { SimulationResults reducedModelResults = simulator.executeSimulation( g_reducedPatternSimulationJob[simulationScenarioIndex], false, false); error += computeError( @@ -188,32 +215,13 @@ double ReducedModelOptimizer::objectiveAllScenarios(long n, const double *x) { 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; + std::vector colors(gObjectiveValueHistory.size(), 2); + if (g_firstRoundIterationIndex != 0) { + for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(), + [](double &c) { c = 0.7; }); } - - 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()); - gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory); + gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory, 6, colors); + std::cout << std::endl; return error; } @@ -257,7 +265,7 @@ void ReducedModelOptimizer::computeMaps( } // Construct reduced->full and full->reduced interface vi map - g_reducedToFullInterfaceViMap.clear(); + reducedToFullInterfaceViMap.clear(); vcg::tri::Allocator::PointerUpdater pu_reducedModel; reducedPattern.deleteDanglingVertices(pu_reducedModel); @@ -267,13 +275,13 @@ void ReducedModelOptimizer::computeMaps( reducedPattern.VN() - 1 /*- reducedModelBaseTriangleInterfaceVi*/; reducedPattern.createFan(); for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { - g_reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex + - reducedModelBaseTriangleInterfaceVi] = + reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex + + reducedModelBaseTriangleInterfaceVi] = fullModelBaseTriangleInterfaceVi + fanIndex * fullPatternInterfaceVertexOffset; } m_fullToReducedInterfaceViMap.clear(); - constructInverseMap(g_reducedToFullInterfaceViMap, + constructInverseMap(reducedToFullInterfaceViMap, m_fullToReducedInterfaceViMap); // fullPattern.setLabel("FullPattern"); @@ -288,6 +296,8 @@ void ReducedModelOptimizer::computeMaps( m_fullPatternOppositeInterfaceViMap[vi0] = vi1; } + g_reducedToFullViMap = reducedToFullInterfaceViMap; + const bool debugMapping = false; if (debugMapping) { reducedPattern.registerForDrawing(); @@ -322,11 +332,11 @@ void ReducedModelOptimizer::computeMaps( std::vector nodeColorsReducedToFull_full(fullPattern.VN(), glm::vec3(0, 0, 0)); for (size_t vi = 0; vi < reducedPattern.VN(); vi++) { - if (g_reducedToFullInterfaceViMap.contains(vi)) { + if (reducedToFullInterfaceViMap.contains(vi)) { auto color = polyscope::getNextUniqueColor(); nodeColorsReducedToFull_reduced[vi] = color; - nodeColorsReducedToFull_full[g_reducedToFullInterfaceViMap[vi]] = color; + nodeColorsReducedToFull_full[reducedToFullInterfaceViMap[vi]] = color; } } polyscope::getCurveNetwork(reducedPattern.getLabel()) @@ -352,10 +362,10 @@ void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel, m_pReducedPatternSimulationMesh->setBeamCrossSection( CrossSectionType{0.002, 0.002}); m_pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1); - m_pFullModelSimulationMesh = std::make_shared(fullModel); - m_pFullModelSimulationMesh->setBeamCrossSection( + m_pFullPatternSimulationMesh = std::make_shared(fullModel); + m_pFullPatternSimulationMesh->setBeamCrossSection( CrossSectionType{0.002, 0.002}); - m_pFullModelSimulationMesh->setBeamMaterial(0.3, 1); + m_pFullPatternSimulationMesh->setBeamMaterial(0.3, 1); } ReducedModelOptimizer::ReducedModelOptimizer( @@ -367,6 +377,8 @@ ReducedModelOptimizer::ReducedModelOptimizer( void ReducedModelOptimizer::initializePatterns( FlatPattern &fullPattern, FlatPattern &reducedPattern, const std::unordered_set &reducedModelExcludedEdges) { + fullPattern.setLabel("full_pattern_" + fullPattern.getLabel()); + reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel()); assert(fullPattern.VN() == reducedPattern.VN() && fullPattern.EN() >= reducedPattern.EN()); polyscope::removeAllStructures(); @@ -375,7 +387,8 @@ void ReducedModelOptimizer::initializePatterns( FlatPattern copyReducedPattern; copyFullPattern.copy(fullPattern); copyReducedPattern.copy(reducedPattern); - if (copyReducedPattern.EN() == 2) { + g_optimizeInnerHexagonSize = copyReducedPattern.EN() == 2; + if (g_optimizeInnerHexagonSize) { const double h = copyReducedPattern.getBaseTriangleHeight(); double baseTriangle_bottomEdgeSize = 2 * h / tan(vcg::math::ToRad(60.0)); VectorType baseTriangle_leftBottomNode(-baseTriangle_bottomEdgeSize / 2, -h, @@ -391,19 +404,20 @@ void ReducedModelOptimizer::initializePatterns( baseTriangle_leftBottomNode; g_innerHexagonVectors[rotationCounter] = rotatedVector; } + const double innerHexagonInitialPos_x = + copyReducedPattern.vert[0].cP()[0] / g_innerHexagonVectors[0][0]; + const double innerHexagonInitialPos_y = + copyReducedPattern.vert[0].cP()[1] / g_innerHexagonVectors[0][1]; + g_innerHexagonInitialPos = innerHexagonInitialPos_x; } computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges); - const double innerHexagonInitialPos_x = - copyReducedPattern.vert[0].cP()[0] / g_innerHexagonVectors[0][0]; - const double innerHexagonInitialPos_y = - copyReducedPattern.vert[0].cP()[1] / g_innerHexagonVectors[0][1]; - g_innerHexagonInitialPos = innerHexagonInitialPos_x; createSimulationMeshes(copyFullPattern, copyReducedPattern); - initializeStiffnesses(); + initializeStiffnesses(m_pReducedPatternSimulationMesh); } -void ReducedModelOptimizer::initializeStiffnesses() { - g_initialParameters.resize(4); +void ReducedModelOptimizer::initializeStiffnesses( + const std::shared_ptr &mesh) { + g_initialParameters.resize(g_optimizeInnerHexagonSize ? 5 : 4); // Save save the beam stiffnesses // for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) { // Element &e = pReducedModelElementalMesh->elements[ei]; @@ -417,73 +431,62 @@ void ReducedModelOptimizer::initializeStiffnesses() { // g_initialParameters = // m_pReducedPatternSimulationMesh->elements[0].properties.E; // for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) { - // g_initialParameters(ei, 0) = - // m_pReducedPatternSimulationMesh->elements[ei].properties.E; // } - // g_initialParameters(1) = - // pReducedModelElementalMesh->elements[0].properties.G; - g_initialParameters(0) = - m_pReducedPatternSimulationMesh->elements[0].properties.A; - g_initialParameters(1) = - m_pReducedPatternSimulationMesh->elements[0].properties.J; - g_initialParameters(2) = - m_pReducedPatternSimulationMesh->elements[0].properties.I2; - g_initialParameters(3) = - m_pReducedPatternSimulationMesh->elements[0].properties.I3; + // 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; // } } void ReducedModelOptimizer::computeReducedModelSimulationJob( const SimulationJob &simulationJobOfFullModel, + const std::unordered_map &simulationJobFullToReducedMap, SimulationJob &simulationJobOfReducedModel) { + assert(simulationJobOfReducedModel.pMesh->VN() != 0); std::unordered_map> reducedModelFixedVertices; - for (auto fullModelFixedVertex : simulationJobOfFullModel.fixedVertices) { - reducedModelFixedVertices[m_fullToReducedInterfaceViMap.at( + for (auto fullModelFixedVertex : + simulationJobOfFullModel.constrainedVertices) { + reducedModelFixedVertices[simulationJobFullToReducedMap.at( fullModelFixedVertex.first)] = fullModelFixedVertex.second; } std::unordered_map reducedModelNodalForces; for (auto fullModelNodalForce : simulationJobOfFullModel.nodalExternalForces) { - reducedModelNodalForces[m_fullToReducedInterfaceViMap.at( + reducedModelNodalForces[simulationJobFullToReducedMap.at( fullModelNodalForce.first)] = fullModelNodalForce.second; } - std::unordered_map reducedModelNodalForcedNormals; - for (auto fullModelNodalForcedRotation : - simulationJobOfFullModel.nodalForcedNormals) { - reducedModelNodalForcedNormals[m_fullToReducedInterfaceViMap.at( - fullModelNodalForcedRotation.first)] = - fullModelNodalForcedRotation.second; - } - - simulationJobOfReducedModel = SimulationJob{m_pReducedPatternSimulationMesh, - reducedModelFixedVertices, - reducedModelNodalForces, - {}, - reducedModelNodalForcedNormals}; -} - -SimulationJob ReducedModelOptimizer::getReducedSimulationJob( - const SimulationJob &fullModelSimulationJob) { - SimulationJob reducedModelSimulationJob; - computeReducedModelSimulationJob(fullModelSimulationJob, - reducedModelSimulationJob); - return reducedModelSimulationJob; + // std::unordered_map + // reducedModelNodalForcedNormals; for (auto fullModelNodalForcedRotation : + // simulationJobOfFullModel.nodalForcedNormals) { + // reducedModelNodalForcedNormals[simulationJobFullToReducedMap.at( + // fullModelNodalForcedRotation.first)] = + // fullModelNodalForcedRotation.second; + // } + simulationJobOfReducedModel.constrainedVertices = reducedModelFixedVertices; + simulationJobOfReducedModel.nodalExternalForces = reducedModelNodalForces; + simulationJobOfReducedModel.label = simulationJobOfFullModel.getLabel(); + // simulationJobOfReducedModel.nodalForcedNormals = + // reducedModelNodalForcedNormals; } void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, + const std::unordered_map &displacementsReducedToFullMap, Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) { - optimalDisplacementsOfReducedModel.resize( - m_pReducedPatternSimulationMesh->VN(), 3); + assert(optimalDisplacementsOfReducedModel.rows() != 0 && + optimalDisplacementsOfReducedModel.cols() == 3); optimalDisplacementsOfReducedModel.setZero( optimalDisplacementsOfReducedModel.rows(), optimalDisplacementsOfReducedModel.cols()); - for (auto reducedFullViPair : g_reducedToFullInterfaceViMap) { + for (auto reducedFullViPair : displacementsReducedToFullMap) { const VertexIndex fullModelVi = reducedFullViPair.second; const Vector6d fullModelViDisplacements = fullModelResults.displacements[fullModelVi]; @@ -499,7 +502,7 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization( gObjectiveValueHistory.clear(); - const size_t n = m_pReducedPatternSimulationMesh->EN() == 12 ? 5 : 4; + const size_t n = g_optimizeInnerHexagonSize ? 5 : 4; const size_t npt = 2 * n; // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2); @@ -507,39 +510,58 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization( "conditions is not recommended."); // Set initial guess of solution - std::vector x(n, 100); - x[4] = g_innerHexagonInitialPos; - if (!initialGuess.empty()) { - x = initialGuess; - } // {0.10000000000000 001, 2, 1.9999999971613847, 6.9560343643347764}; + std::vector x(n, 50); + if (g_optimizeInnerHexagonSize) { + x[n - 1] = g_innerHexagonInitialPos; + } + /*if (!initialGuess.empty()) { + x = g_optimizationInitialGuess; + }*/ // {0.10000000000000 001, + // 2, 1.9999999971613847, 6.9560343643347764}; // {1, 5.9277}; // {0.0001, 2, 2.000000005047502, 1.3055270196964464}; // {initialGuess(0), initialGuess(1), initialGuess(2), // initialGuess(3)}; const double xMin = 1e-2; - const double xMax = 5000; + const double xMax = 500; // assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) { // return d >= xMax || d <= xMin; // })); std::vector xLow(x.size(), xMin); std::vector xUpper(x.size(), xMax); - if (n == 5) { - xLow[4] = 0.1; - xUpper[4] = 1; + if (g_optimizeInnerHexagonSize) { + xLow[n - 1] = 0.1; + xUpper[n - 1] = 0.9; } const double maxX = *std::max_element( x.begin(), x.end(), [](const double &a, const double &b) { return abs(a) < abs(b); }); // const double rhobeg = std::min(0.95, 0.2 * maxX); - const double rhobeg = 0.01; + const double rhobeg = 1; // const double rhobeg = 10; const double rhoend = rhobeg * 1e-6; const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2; std::vector w(wSize); // const size_t maxFun = 10 * (x.size() ^ 2); - const size_t maxFun = 120; + const size_t maxFun = 200; bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(), rhobeg, rhoend, maxFun, w.data()); + std::cout << "Finished first optimization round" << std::endl; + firstOptimizationRoundResults.resize(6); + for (int simulationScenarioIndex = SimulationScenario::Axial; + simulationScenarioIndex != + SimulationScenario::NumberOfSimulationScenarios; + simulationScenarioIndex++) { + SimulationResults reducedModelResults = simulator.executeSimulation( + g_reducedPatternSimulationJob[simulationScenarioIndex], false, false); + reducedModelResults.setLabelPrefix("FirstRound"); + firstOptimizationRoundResults[simulationScenarioIndex] = + std::move(reducedModelResults); + } + g_firstRoundIterationIndex = gObjectiveValueHistory.size(); + bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(), + rhobeg * 1e-3, 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++) { @@ -552,12 +574,13 @@ void ReducedModelOptimizer::setInitialGuess(std::vector v) { initialGuess = v; } -std::vector ReducedModelOptimizer::createScenarios( +std::vector> +ReducedModelOptimizer::createScenarios( const std::shared_ptr &pMesh) { - std::vector scenarios; + std::vector> scenarios; + scenarios.resize(SimulationScenario::NumberOfSimulationScenarios); std::unordered_map> fixedVertices; std::unordered_map nodalForces; - // NOTE: Assuming that the first interface node lays on the y axis const double forceMagnitude = 1; // Assuming the patterns lays on the x-y plane const CoordType patternPlaneNormal(0, 0, 1); @@ -594,6 +617,7 @@ std::vector ReducedModelOptimizer::createScenarios( // m_pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); //// Axial + SimulationScenario scenarioName = SimulationScenario::Axial; for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { CoordType forceDirection = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) @@ -604,9 +628,12 @@ std::vector ReducedModelOptimizer::createScenarios( fixedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; } - scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + scenarios[scenarioName] = std::make_shared( + SimulationJob(pMesh, simulationScenarioStrings[scenarioName], + fixedVertices, nodalForces, {})); - //// In-plane Bending + //// Shear + scenarioName = SimulationScenario::Shear; fixedVertices.clear(); nodalForces.clear(); for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { @@ -620,7 +647,9 @@ std::vector ReducedModelOptimizer::createScenarios( fixedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; } - scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + scenarios[scenarioName] = std::make_shared( + SimulationJob(pMesh, simulationScenarioStrings[scenarioName], + fixedVertices, nodalForces, {})); // //// Torsion // fixedVertices.clear(); @@ -645,7 +674,8 @@ std::vector ReducedModelOptimizer::createScenarios( // } // scenarios.push_back({pMesh, fixedVertices, nodalForces}); - //// Out - of - plane bending.Pull towards Z + //// Bending + scenarioName = SimulationScenario::Bending; fixedVertices.clear(); nodalForces.clear(); for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { @@ -653,9 +683,12 @@ std::vector ReducedModelOptimizer::createScenarios( fixedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; } - scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + scenarios[scenarioName] = std::make_shared( + SimulationJob(pMesh, simulationScenarioStrings[scenarioName], + fixedVertices, nodalForces, {})); //// Double using moments + scenarioName = SimulationScenario::Dome; fixedVertices.clear(); nodalForces.clear(); for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); @@ -676,9 +709,12 @@ std::vector ReducedModelOptimizer::createScenarios( nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude * 0.1; } - scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + scenarios[scenarioName] = std::make_shared( + SimulationJob(pMesh, simulationScenarioStrings[scenarioName], + fixedVertices, nodalForces, {})); //// Saddle + scenarioName = SimulationScenario::Saddle; fixedVertices.clear(); nodalForces.clear(); for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); @@ -702,76 +738,177 @@ std::vector ReducedModelOptimizer::createScenarios( Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.01 * forceMagnitude; } } - scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + scenarios[scenarioName] = std::make_shared( + SimulationJob(pMesh, simulationScenarioStrings[scenarioName], + fixedVertices, nodalForces, {})); return scenarios; } -void ReducedModelOptimizer::optimize(const int &simulationScenario) { +void ReducedModelOptimizer::runBeamOptimization() { + // load beams + VCGEdgeMesh fullBeam; + fullBeam.loadFromPly("/home/iason/Models/simple_beam_model_10elem_1m.ply"); + VCGEdgeMesh reducedBeam; + reducedBeam.loadFromPly("/home/iason/Models/simple_beam_model_4elem_1m.ply"); + fullBeam.registerForDrawing(); + reducedBeam.registerForDrawing(); + // polyscope::show(); + // maps + std::unordered_map displacementReducedToFullMap; + displacementReducedToFullMap[reducedBeam.VN() / 2] = fullBeam.VN() / 2; + g_reducedToFullViMap = displacementReducedToFullMap; + std::unordered_map jobFullToReducedMap; + jobFullToReducedMap[0] = 0; + jobFullToReducedMap[fullBeam.VN() - 1] = reducedBeam.VN() - 1; - std::vector simulationJobs = - createScenarios(m_pFullModelSimulationMesh); + // full model simuilation job + auto pFullPatternSimulationMesh = std::make_shared(fullBeam); + pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02}); + pFullPatternSimulationMesh->setBeamMaterial(0.3, 1); + std::unordered_map> fixedVertices; + fixedVertices[0] = ::unordered_set({0, 1, 2, 3, 4, 5}); + std::unordered_map forces; + forces[fullBeam.VN() - 1] = Vector6d({0, 0, 10, 0, 0, 0}); + const std::string fullBeamSimulationJobLabel = "Pull_Z"; + std::shared_ptr pFullModelSimulationJob = + make_shared(SimulationJob(pFullPatternSimulationMesh, + fullBeamSimulationJobLabel, + fixedVertices, forces)); + FormFinder formFinder; + auto fullModelResults = formFinder.executeSimulation(pFullModelSimulationJob); + + // Optimal reduced model displacements + const size_t numberOfSimulationScenarios = 1; + g_optimalReducedModelDisplacements.resize(numberOfSimulationScenarios); + g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1].resize( + reducedBeam.VN(), 3); + computeDesiredReducedModelDisplacements( + fullModelResults, displacementReducedToFullMap, + g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]); + + // reduced model simuilation job + auto reducedSimulationMesh = std::make_shared(reducedBeam); + reducedSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02}); + reducedSimulationMesh->setBeamMaterial(0.3, 1); + g_reducedPatternSimulationJob.resize(numberOfSimulationScenarios); + SimulationJob reducedSimJob; + computeReducedModelSimulationJob(*pFullModelSimulationJob, + jobFullToReducedMap, reducedSimJob); + reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] = + reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] * 0.1; + g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1] = + make_shared( + reducedSimulationMesh, fullBeamSimulationJobLabel, + reducedSimJob.constrainedVertices, reducedSimJob.nodalExternalForces); + initializeStiffnesses(reducedSimulationMesh); + + // const std::string simulationJobsPath = "SimulationJobs"; + // std::filesystem::create_directory(simulationJobsPath); + // g_reducedPatternSimulationJob[0].save(simulationJobsPath); + // g_reducedPatternSimulationJob[0].load( + // std::filesystem::path(simulationJobsPath) + // .append(g_reducedPatternSimulationJob[0].mesh->getLabel() + + // "_simScenario.json")); + + runOptimization(&objective); + + fullModelResults.registerForDrawing(); + SimulationResults reducedModelResults = simulator.executeSimulation( + g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1], false, + false); + double error = computeError( + reducedModelResults, + g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]); + reducedModelResults.registerForDrawing(); + std::cout << "Error between beams:" << error << endl; + // registerWorldAxes(); + polyscope::show(); + fullModelResults.unregister(); + reducedModelResults.unregister(); +} + +void ReducedModelOptimizer::visualizeResults( + const std::vector> + &fullPatternSimulationJobs, + const std::vector &simulationScenarios) { + m_pFullPatternSimulationMesh->registerForDrawing(); + double error = 0; + for (const int simulationScenarioIndex : simulationScenarios) { + const std::shared_ptr &pFullPatternSimulationJob = + fullPatternSimulationJobs[simulationScenarioIndex]; + pFullPatternSimulationJob->registerForDrawing( + m_pFullPatternSimulationMesh->getLabel()); + SimulationResults fullModelResults = + simulator.executeSimulation(pFullPatternSimulationJob, false); + fullModelResults.registerForDrawing(); + fullModelResults.saveDeformedModel(); + const std::shared_ptr &pReducedPatternSimulationJob = + g_reducedPatternSimulationJob[simulationScenarioIndex]; + SimulationResults reducedModelResults = + simulator.executeSimulation(pReducedPatternSimulationJob, false, false); + error += computeError( + reducedModelResults, + g_optimalReducedModelDisplacements[simulationScenarioIndex]); + std::cout << "Error of simulation scenario " + << simulationScenarioStrings[simulationScenarioIndex] << " is " + << error << std::endl; + reducedModelResults.registerForDrawing(); + firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing(); + reducedModelResults.saveDeformedModel(); + // registerWorldAxes(); + polyscope::show(); + fullModelResults.unregister(); + reducedModelResults.unregister(); + firstOptimizationRoundResults[simulationScenarioIndex].unregister(); + } +} + +void ReducedModelOptimizer::optimize( + const std::vector &simulationScenariosIndices) { + + std::vector> simulationJobs = + createScenarios(m_pFullPatternSimulationMesh); g_optimalReducedModelDisplacements.resize(6); g_reducedPatternSimulationJob.resize(6); + g_firstRoundIterationIndex = 0; // polyscope::removeAllStructures(); for (int simulationScenarioIndex = SimulationScenario::Axial; simulationScenarioIndex != SimulationScenario::NumberOfSimulationScenarios; simulationScenarioIndex++) { - const SimulationJob &fullPatternSimulationJob = + const std::shared_ptr &pFullPatternSimulationJob = simulationJobs[simulationScenarioIndex]; // fullPatternSimulationJob.mesh->savePly( // "Fanned_" + m_pFullModelSimulationMesh->getLabel() + ".ply"); SimulationResults fullModelResults = - simulator.executeSimulation(fullPatternSimulationJob, false); + simulator.executeSimulation(pFullPatternSimulationJob); // fullModelResults.label = // "fullModel_" + SimulationScenarioStrings[simulationScenarioIndex]; // fullModelResults.registerForDrawing(fullPatternSimulationJob); + g_optimalReducedModelDisplacements[simulationScenarioIndex].resize( + m_pReducedPatternSimulationMesh->VN(), 3); computeDesiredReducedModelDisplacements( - fullModelResults, + fullModelResults, g_reducedToFullViMap, g_optimalReducedModelDisplacements[simulationScenarioIndex]); - computeReducedModelSimulationJob( - fullPatternSimulationJob, - g_reducedPatternSimulationJob[simulationScenarioIndex]); + SimulationJob reducedPatternSimulationJob; + reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; + computeReducedModelSimulationJob(*pFullPatternSimulationJob, + m_fullToReducedInterfaceViMap, + reducedPatternSimulationJob); + g_reducedPatternSimulationJob[simulationScenarioIndex] = + std::make_shared(reducedPatternSimulationJob); } - if (simulationScenario == -1) { - double (*pObjectiveFunction)(long, const double *) = &objectiveAllScenarios; - runOptimization(pObjectiveFunction); - } else { // run chosen - g_chosenSimulationScenarioName = SimulationScenario(simulationScenario); - double (*pObjectiveFunction)(long, const double *) = - &objectiveSingleScenario; - runOptimization(pObjectiveFunction); + g_simulationScenarioIndices = simulationScenariosIndices; + if (simulationScenariosIndices.empty()) { + g_simulationScenarioIndices = { + SimulationScenario::Axial, SimulationScenario::Shear, + SimulationScenario::Bending, SimulationScenario::Dome, + SimulationScenario::Saddle}; } - 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(); - polyscope::show(); - fullModelResults.unregister(); - reducedModelResults.unregister(); - } + runOptimization(&objective); + visualizeResults(simulationJobs, g_simulationScenarioIndices); } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index a687f6d..eb7a52c 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -12,7 +12,7 @@ using ReducedPatternVertexIndex = VertexIndex; class ReducedModelOptimizer { std::shared_ptr m_pReducedPatternSimulationMesh; - std::shared_ptr m_pFullModelSimulationMesh; + std::shared_ptr m_pFullPatternSimulationMesh; std::unordered_map m_fullToReducedInterfaceViMap; std::unordered_map @@ -26,18 +26,19 @@ public: Axial, Shear, Bending, - Double, + Dome, Saddle, NumberOfSimulationScenarios }; - inline static const std::string SimulationScenarioStrings[] = { + inline static const std::string simulationScenarioStrings[] = { "Axial", "Shear", "Bending", "Double", "Saddle"}; - void optimize(const int &simulationScenario = -1); + void optimize(const std::vector &simulationScenarios); double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot); - void computeReducedModelSimulationJob( + static void computeReducedModelSimulationJob( const SimulationJob &simulationJobOfFullModel, + const std::unordered_map &simulationJobFullToReducedMap, SimulationJob &simulationJobOfReducedModel); SimulationJob @@ -49,25 +50,35 @@ public: void setInitialGuess(std::vector v); + static void runBeamOptimization(); + + static void runSimulation(const std::string &filename, + std::vector &x); + private: - void computeDesiredReducedModelDisplacements( + void + visualizeResults(const std::vector> + &fullPatternSimulationJobs, + const std::vector &simulationScenarios); + static void computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, + const std::unordered_map &displacementsReducedToFullMap, Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); - Eigen::VectorXd runOptimization(double (*pObjectiveFunction)(long, - const double *)); - std::vector + static Eigen::VectorXd + runOptimization(double (*pObjectiveFunction)(long, const double *)); + std::vector> createScenarios(const std::shared_ptr &pMesh); void computeMaps(FlatPattern &fullModel, FlatPattern &reducedPattern, const std::unordered_set &reducedModelExcludedEges); void createSimulationMeshes(FlatPattern &fullModel, FlatPattern &reducedModel); - void initializeStiffnesses(); + static void + initializeStiffnesses(const std::shared_ptr &mesh); 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); + static double objective(long n, const double *x); }; #endif // REDUCEDMODELOPTIMIZER_HPP