From 4d4aa346b5b42de8bc92a6512b32adaf219fbe40 Mon Sep 17 00:00:00 2001 From: Iason Date: Wed, 9 Dec 2020 17:58:48 +0200 Subject: [PATCH] Created simulation scenarios. --- src/main.cpp | 181 ++++++----- src/reducedmodeloptimizer.cpp | 592 ++++++++++++++++++---------------- src/reducedmodeloptimizer.hpp | 3 + 3 files changed, 420 insertions(+), 356 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 7e38d96..9948fa9 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -6,18 +6,42 @@ #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"); - registerWorldAxes(); // Create reduced models const std::vector numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; std::vector singleBarReducedModelEdges{vcg::Point2i(0, 3)}; @@ -25,50 +49,56 @@ int main(int argc, char *argv[]) { 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 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 CWreducedModelEdges{vcg::Point2i(1, 5), + // vcg::Point2i(3, 1)}; + // FlatPattern CCWReducedModel(numberOfNodesPerSlot, CWreducedModelEdges); + // CCWReducedModel.setLabel("CCW reduced model"); - std::vector reducedModels{&singleBarReducedModel, - &CWReducedModel, &CCWReducedModel}; + // std::vector reducedModels{&singleBarReducedModel, + // &CWReducedModel, + // &CCWReducedModel}; - ReducedModelOptimizer optimizer(numberOfNodesPerSlot); - std::string fullPatternsTestSetDirectory = - "/home/iason/Models/TestSet_validPatterns"; - for (const auto &entry : - filesystem::directory_iterator(fullPatternsTestSetDirectory)) { - const auto filepath = - // entry.path(); - std::filesystem::path("/home/iason/Models/valid_6777.ply"); - 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; - } - FlatPattern pattern(filepathString); - pattern.setLabel(filepath.stem().string()); - std::cout << "Testing Pattern:" << filepathString << std::endl; - for (FlatPattern *pReducedModel : reducedModels) { - pReducedModel = reducedModels[0]; - std::unordered_set optimizationExcludedEi; - if (pReducedModel != - reducedModels[0]) { // assumes that the singleBar reduced model is the - // first in the reducedModels vector - optimizationExcludedEi.insert(0); - } - optimizer.initialize(pattern, *pReducedModel, optimizationExcludedEi); - optimizer.initialize(pattern, *pReducedModel, optimizationExcludedEi); - Eigen::VectorXd optimalParameters = optimizer.optimize(); - } - } + // ReducedModelOptimizer optimizer(numberOfNodesPerSlot); + // std::string fullPatternsTestSetDirectory = + // "/home/iason/Models/TestSet_validPatterns"; + // 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/431.ply"); + // // std::filesystem::path( + // // "/home/iason/Models/TestSet_validPatterns/865.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; + // } + // FlatPattern pattern(filepathString); + // pattern.setLabel(filepath.stem().string()); + // std::cout << "Testing Pattern:" << filepathString << std::endl; + // for (FlatPattern *pReducedModel : reducedModels) { + // // pReducedModel = reducedModels[1]; + // std::unordered_set optimizationExcludedEi; + // if (pReducedModel != + // reducedModels[0]) { // assumes that the singleBar reduced model is + // // the + // // first in the reducedModels vector + // optimizationExcludedEi.insert(0); + // } + // optimizer.initialize(pattern, *pReducedModel, optimizationExcludedEi); + // Eigen::VectorXd optimalParameters = optimizer.optimize(); + // } + // } // // Full model simulation // std::unordered_map> @@ -135,43 +165,34 @@ int main(int argc, char *argv[]) { // stiffnessFactor *= 1.5; // } - // Beam - // VCGEdgeMesh beamMesh; - // beamMesh.loadFromPly("/home/iason/Models/simple_beam_model_2elem_1m.ply"); - // const VertexIndex reducedModelOpposite_vi = - // std::ceil(reducedModelSimulationJob.mesh->VN() / 2.0); - // auto v0 = - // reducedModelSimulationJob.mesh->vert[reducedModelOpposite_vi].cP() - - // reducedModelSimulationJob.mesh->vert[1].cP(); - // auto v1 = beamMesh.vert[2].cP() - beamMesh.vert[0].cP(); - // vcg::Matrix44d R; - // const double rotationTheta = - // std::asin((v0 ^ v1).Norm() / (v0.Norm() * v1.Norm())); - // R.SetRotateRad(-rotationTheta, v0 ^ v1); - // vcg::tri::UpdatePosition::Matrix(beamMesh, R); - // vcg::tri::UpdatePosition::Scale(beamMesh, v0.Norm() / - // v1.Norm()); vcg::tri::UpdatePosition::Translate( - // beamMesh, reducedModelSimulationJob.mesh->vert[1].cP()); - // std::unordered_map> - // beamFixedVertices; - // beamFixedVertices[0] = std::unordered_set{0, 1, 2, 3, 4, 5}; - // std::unordered_map beamNodalForces{ - // {2, Eigen::Vector3d(0, 0, 1000)}}; - // SimulationJob beamSimulationJob{std::make_shared(beamMesh), - // beamFixedVertices, - // beamNodalForces, - // {}}; - // // for (EdgeIndex ei = 0; ei < beamSimulationJob.mesh->EN(); ei++) { - // // BeamFormFinder::Element &e = beamSimulationJob.mesh->elements[ei]; - // // e.properties.A = stifnessVector[0]; - // // e.properties.J = stifnessVector[1]; - // // e.properties.I2 = stifnessVector[2]; - // // e.properties.I3 = stifnessVector[3]; - // // } - // // beamSimulationJob.draw(); - // SimulationResults beamSimulationResults = - // formFinder.executeSimulation(beamSimulationJob); - // beamSimulationResults.simulationLabel = "Beam"; - // beamSimulationResults.draw(beamSimulationJob); + // Beam + VCGEdgeMesh mesh; + mesh.loadFromPly("/home/iason/Models/simple_beam_model_10elem_1m.ply"); + FormFinder formFinder; + std::unordered_map> fixedVertices; + fixedVertices[0] = std::unordered_set{0, 1, 2, 3, 4, 5}; + std::unordered_map nodalForces{ + {10, Vector6d({0, 0, 0, 0.00001, 0, 0})}}; + // Forced displacements + std::unordered_map nodalForcedDisplacements; + std::unordered_map nodalForcedNormal; + + // CoordType v = (mesh.vert[10].cP() - mesh.vert[0].cP()).Normalize(); + // nodalForcedNormal[10] = v; + fixedVertices[10] = {0, 1, 2}; + + SimulationJob beamSimulationJob{std::make_shared(mesh), + fixedVertices, nodalForces, + nodalForcedDisplacements, nodalForcedNormal}; + beamSimulationJob.mesh->setBeamMaterial(0.3, 1); + beamSimulationJob.mesh->setBeamCrossSection( + RectangularBeamDimensions{0.002, 0.002}); + // registerWorldAxes(); + SimulationResults beamSimulationResults = + formFinder.executeSimulation(beamSimulationJob, true, true); + beamSimulationResults.simulationLabel = "Beam"; + beamSimulationResults.registerForDrawing(beamSimulationJob); + polyscope::show(); + return 0; } diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index c026a9f..3dd77e0 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -2,126 +2,137 @@ #include "bobyqa.h" #include "flatpattern.hpp" #include "gradientDescent.h" -#include "matplot/matplot.h" #include "simulationhistoryplotter.hpp" #include "trianglepattterntopology.hpp" const bool gShouldDraw = true; -size_t g_numberOfOptimizationRounds{0}; + FormFinder simulator; Eigen::MatrixX3d g_optimalReducedModelDisplacements; -SimulationJob gReducedPatternSimulationJob; +SimulationJob g_reducedPatternSimulationJob; std::unordered_map g_reducedToFullInterfaceViMap; matplot::line_handle gPlotHandle; std::vector gObjectiveValueHistory; Eigen::Vector4d g_initialX; std::unordered_set g_reducedPatternExludedEdges; -Eigen::MatrixX4d g_initialStiffnessFactors; +double g_initialParameters; -struct OptimizationCallback { - double operator()(const size_t &iterations, const Eigen::VectorXd &x, - const double &fval, Eigen::VectorXd &gradient) const { - // run simulation - // SimulationResults reducedModelResults = - // simulator.executeSimulation(reducedModelSimulationJob); - // reducedModelResults.draw(reducedModelSimulationJob); - gObjectiveValueHistory.push_back(fval); - auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), - gObjectiveValueHistory.size()); - gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory); - // const std::string plotImageFilename = "objectivePlot.png"; - // matplot::save(plotImageFilename); - // if (numberOfOptimizationRounds % 30 == 0) { - // std::filesystem::copy_file( - // std::filesystem::path(plotImageFilename), - // std::filesystem::path("objectivePlot_copy.png")); +// struct OptimizationCallback { +// double operator()(const size_t &iterations, const Eigen::VectorXd &x, +// const double &fval, Eigen::VectorXd &gradient) const { +// // run simulation +// // SimulationResults reducedModelResults = +// // simulator.executeSimulation(reducedModelSimulationJob); +// // reducedModelResults.draw(reducedModelSimulationJob); +// gObjectiveValueHistory.push_back(fval); +// auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), +// gObjectiveValueHistory.size()); +// gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory); +// // const std::string plotImageFilename = "objectivePlot.png"; +// // matplot::save(plotImageFilename); +// // if (numberOfOptimizationRounds % 30 == 0) { +// // std::filesystem::copy_file( +// // std::filesystem::path(plotImageFilename), +// // std::filesystem::path("objectivePlot_copy.png")); +// // } +// // std::stringstream ss; +// // ss << x; +// // reducedModelResults.simulationLabel = ss.str(); +// // SimulationResultsReporter resultsReporter; +// // resultsReporter.reportResults( +// // {reducedModelResults}, +// // std::filesystem::current_path().append("Results")); + +// return true; +// } +//}; + +// struct Objective { +// double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const { +// assert(x.rows() == 4); + +// // drawSimulationJob(simulationJob); +// // Set mesh from x +// std::shared_ptr reducedModel = +// g_reducedPatternSimulationJob.mesh; +// for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) { +// if (g_reducedPatternExludedEdges.contains(ei)) { +// continue; +// } +// Element &e = reducedModel->elements[ei]; +// e.axialConstFactor = g_initialStiffnessFactors(ei, 0) * x(0); +// e.torsionConstFactor = g_initialStiffnessFactors(ei, 1) * x(1); +// e.firstBendingConstFactor = g_initialStiffnessFactors(ei, 2) * x(2); +// e.secondBendingConstFactor = g_initialStiffnessFactors(ei, 3) * x(3); +// } +// // run simulation +// SimulationResults reducedModelResults = +// simulator.executeSimulation(g_reducedPatternSimulationJob); +// // std::stringstream ss; +// // ss << x; +// // reducedModelResults.simulationLabel = ss.str(); +// // SimulationResultsReporter resultsReporter; +// // resultsReporter.reportResults( +// // {reducedModelResults}, +// // std::filesystem::current_path().append("Results")); +// // compute error and return it +// double error = 0; +// for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) { +// VertexIndex reducedModelVi = reducedFullViPair.first; +// 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.norm(); +// } + +// return error; +// } +//}; + +double ReducedModelOptimizer::objective(long n, const double *x) { + + for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { + std::cout << "x[" + std::to_string(parameterIndex) + "]=" + << x[parameterIndex] << std::endl; + } + + for (EdgeIndex ei = 0; ei < g_reducedPatternSimulationJob.mesh->EN(); ei++) { + Element &e = g_reducedPatternSimulationJob.mesh->elements[ei]; + // if (g_reducedPatternExludedEdges.contains(ei)) { + // continue; // } - // std::stringstream ss; - // ss << x; - // reducedModelResults.simulationLabel = ss.str(); - // SimulationResultsReporter resultsReporter; - // resultsReporter.reportResults( - // {reducedModelResults}, - // std::filesystem::current_path().append("Results")); - - return true; - } -}; - -struct Objective { - double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const { - assert(x.rows() == 4); - - // drawSimulationJob(simulationJob); - // Set mesh from x - std::shared_ptr reducedModel = - gReducedPatternSimulationJob.mesh; - for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) { - if (g_reducedPatternExludedEdges.contains(ei)) { - continue; - } - Element &e = reducedModel->elements[ei]; - e.axialConstFactor *= x(0); - e.torsionConstFactor *= x(1); - e.firstBendingConstFactor *= x(2); - e.secondBendingConstFactor *= x(3); - } - // run simulation - SimulationResults reducedModelResults = - simulator.executeSimulation(gReducedPatternSimulationJob); - // std::stringstream ss; - // ss << x; - // reducedModelResults.simulationLabel = ss.str(); - // SimulationResultsReporter resultsReporter; - // resultsReporter.reportResults( - // {reducedModelResults}, - // std::filesystem::current_path().append("Results")); - // compute error and return it - double error = 0; - for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) { - VertexIndex reducedModelVi = reducedFullViPair.first; - 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.norm(); - } - - return error; - } -}; - -double objective(long n, const double *x) { - Eigen::VectorXd eigenX(n, 1); - for (size_t xi = 0; xi < n; xi++) { - eigenX(xi) = x[xi]; - } - - std::shared_ptr reducedPattern = - gReducedPatternSimulationJob.mesh; - for (EdgeIndex ei = 0; ei < reducedPattern->EN(); ei++) { - Element &e = reducedPattern->elements[ei]; - if (g_reducedPatternExludedEdges.contains(ei)) { - continue; - } - e.axialConstFactor = g_initialStiffnessFactors(ei, 0) * eigenX(0); - e.torsionConstFactor = g_initialStiffnessFactors(ei, 1) * eigenX(1); - e.firstBendingConstFactor = g_initialStiffnessFactors(ei, 2) * eigenX(2); - e.secondBendingConstFactor = g_initialStiffnessFactors(ei, 3) * eigenX(3); + e.properties.E = g_initialParameters * x[0]; + // e.properties.G = g_initialParameters(1) * x[1]; + // e.properties.A = g_initialParameters(0, 0) * x[0]; + // e.properties.J = g_initialParameters(0, 1) * x[1]; + // e.properties.I2 = g_initialParameters(0, 2) * x[2]; + // e.properties.I3 = g_initialParameters(0, 3) * x[3]; + e.properties.G = e.properties.E / (2 * (1 + 0.3)); + e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength; + e.torsionConstFactor = e.properties.G * e.properties.J / e.initialLength; + e.firstBendingConstFactor = + 2 * e.properties.E * e.properties.I2 / e.initialLength; + e.secondBendingConstFactor = + 2 * e.properties.E * e.properties.I3 / e.initialLength; } // run simulation SimulationResults reducedModelResults = - simulator.executeSimulation(gReducedPatternSimulationJob, false, false); + simulator.executeSimulation(g_reducedPatternSimulationJob, false, false); // 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], @@ -130,6 +141,7 @@ double objective(long n, const double *x) { Eigen::Vector3d( g_optimalReducedModelDisplacements.row(reducedModelVi)) - vertexDisplacement; + // error += errorVector.squaredNorm(); error += errorVector.norm(); } @@ -267,7 +279,12 @@ void ReducedModelOptimizer::computeMaps( void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel, FlatPattern &reducedModel) { pReducedModelElementalMesh = std::make_shared(reducedModel); + pReducedModelElementalMesh->setBeamCrossSection( + CrossSectionType{0.002, 0.002}); + pReducedModelElementalMesh->setBeamMaterial(0.3, 1); pFullModelElementalMesh = std::make_shared(fullModel); + pFullModelElementalMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002}); + pFullModelElementalMesh->setBeamMaterial(0.3, 1); } ReducedModelOptimizer::ReducedModelOptimizer( @@ -296,22 +313,29 @@ void ReducedModelOptimizer::initialize( } void ReducedModelOptimizer::initializeStiffnesses() { - g_initialStiffnessFactors.resize(pReducedModelElementalMesh->EN(), 4); + // g_initialParameters.resize(1, 2); // Save save the beam stiffnesses - for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) { - Element &e = pReducedModelElementalMesh->elements[ei]; - // if (g_reducedPatternExludedEdges.contains(ei)) { - // const double stiffnessFactor = 1; - // e.axialConstFactor *= stiffnessFactor; - // e.torsionConstFactor *= stiffnessFactor; - // e.firstBendingConstFactor *= stiffnessFactor; - // e.secondBendingConstFactor *= stiffnessFactor; - // } - g_initialStiffnessFactors(ei, 0) = e.axialConstFactor; - g_initialStiffnessFactors(ei, 1) = e.torsionConstFactor; - g_initialStiffnessFactors(ei, 2) = e.firstBendingConstFactor; - g_initialStiffnessFactors(ei, 3) = e.secondBendingConstFactor; - } + // for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) { + // Element &e = pReducedModelElementalMesh->elements[ei]; + // if (g_reducedPatternExludedEdges.contains(ei)) { + // const double stiffnessFactor = 5; + // e.axialConstFactor *= stiffnessFactor; + // e.torsionConstFactor *= stiffnessFactor; + // e.firstBendingConstFactor *= stiffnessFactor; + // e.secondBendingConstFactor *= stiffnessFactor; + // } + g_initialParameters = pReducedModelElementalMesh->elements[0].properties.E; + // g_initialParameters(1) = + // pReducedModelElementalMesh->elements[0].properties.G; + // g_initialParameters(0, 0) = + // pReducedModelElementalMesh->elements[0].properties.A; + // g_initialParameters(0, 1) = + // pReducedModelElementalMesh->elements[0].properties.J; + // g_initialParameters(0, 2) = + // pReducedModelElementalMesh->elements[0].properties.I2; + // g_initialParameters(0, 3) = + // pReducedModelElementalMesh->elements[0].properties.I3; + // } } void ReducedModelOptimizer::computeReducedModelSimulationJob( @@ -331,10 +355,19 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob( fullModelNodalForce.first)] = fullModelNodalForce.second; } + std::unordered_map reducedModelNodalForcedNormals; + for (auto fullModelNodalForcedRotation : + simulationJobOfFullModel.nodalForcedNormals) { + reducedModelNodalForcedNormals[m_fullToReducedInterfaceViMap.at( + fullModelNodalForcedRotation.first)] = + fullModelNodalForcedRotation.second; + } + simulationJobOfReducedModel = SimulationJob{pReducedModelElementalMesh, reducedModelFixedVertices, reducedModelNodalForces, - {}}; + {}, + reducedModelNodalForcedNormals}; } SimulationJob ReducedModelOptimizer::getReducedSimulationJob( @@ -349,6 +382,7 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) { + optimalDisplacementsOfReducedModel.resize(0, 3); optimalDisplacementsOfReducedModel.resize(pReducedModelElementalMesh->VN(), 3); optimalDisplacementsOfReducedModel.setZero( @@ -367,115 +401,88 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( } Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( - const SimulationJob &fullModelSimulationJob) { + const SimulationJob &fullPatternSimulationJob) { gObjectiveValueHistory.clear(); + // fullPatternSimulationJob.registerForDrawing(); + // polyscope::show(); + fullPatternSimulationJob.mesh->savePly( + "Fanned_" + pFullModelElementalMesh->getLabel() + ".ply"); SimulationResults fullModelResults = - simulator.executeSimulation(fullModelSimulationJob, false, false); + simulator.executeSimulation(fullPatternSimulationJob, false, true); fullModelResults.simulationLabel = "fullModel"; computeDesiredReducedModelDisplacements(fullModelResults, g_optimalReducedModelDisplacements); - computeReducedModelSimulationJob(fullModelSimulationJob, - gReducedPatternSimulationJob); + computeReducedModelSimulationJob(fullPatternSimulationJob, + g_reducedPatternSimulationJob); // fullModelSimulationJob.registerForDrawing(); // polyscope::show(); // gReducedPatternSimulationJob.registerForDrawing(); // polyscope::show(); - fullModelResults.registerForDrawing(fullModelSimulationJob); + fullModelResults.registerForDrawing(fullPatternSimulationJob); polyscope::show(); + double (*pObjectiveFunction)(long, const double *) = &objective; + const size_t n = 1; + // g_initialParameters.rows(); + const size_t npt = ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; + assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2); + assert(npt <= 2 * n + 1 && "The choice of the number of interpolation " + "conditions is not recommended."); // Set initial guess of solution - Eigen::VectorXd initialGuess(4); - const double stifnessFactor = 1; - initialGuess(0) = stifnessFactor; - initialGuess(1) = stifnessFactor; - initialGuess(2) = stifnessFactor; - initialGuess(3) = stifnessFactor; + const double initialGuess = 1; + std::vector x(n, initialGuess); + // {1, 5.9277}; + // {initialGuess(0), initialGuess(1), initialGuess(2), + // initialGuess(3)}; + const double xMin = 1; + const double xMax = 100; + // 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); + 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 = 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; + bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(), + rhobeg, rhoend, maxFun, w.data()); - const bool useGradientDescent = false; - if (useGradientDescent) { + SimulationResults reducedModelOptimizedResults = + simulator.executeSimulation(g_reducedPatternSimulationJob); - // gdc::GradientDescent, - // OptimizationCallback> - gdc::GradientDescent, - OptimizationCallback> - // gdc::GradientDescent, - // OptimizationCallback> - // gdc::GradientDescent, - // OptimizationCallback> - optimizer; - - // Turn verbosity on, so the optimizer prints status updates after each - // iteration. - optimizer.setVerbosity(1); - - // Set initial guess. - matplot::xlabel("Optimization iterations"); - matplot::ylabel("Objective value"); - // matplot::figure(false); - matplot::grid(matplot::on); - // Start the optimization - auto result = optimizer.minimize(initialGuess); - - std::cout << "Done! Converged: " << (result.converged ? "true" : "false") - << " Iterations: " << result.iterations << std::endl; - - // do something with final function value - std::cout << "Final fval: " << result.fval << std::endl; - - // do something with final x-value - std::cout << "Final xval: " << result.xval.transpose() << std::endl; - - SimulationResults reducedModelOptimizedResults = - simulator.executeSimulation(gReducedPatternSimulationJob); - reducedModelOptimizedResults.simulationLabel = "reducedModel"; - reducedModelOptimizedResults.registerForDrawing( - gReducedPatternSimulationJob); - - return result.xval; - } else { // use bobyqa - double (*pObjectiveFunction)(long, const double *) = &objective; - const size_t n = 4; - const size_t npt = 8; - assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2); - assert(npt < 2 * n + 1 && "The choice of the number of interpolation " - "conditions is not recommended."); - std::vector x - // {1.03424, 0.998456, 0.619916, -0.202997}; - {initialGuess(0), initialGuess(1), initialGuess(2), initialGuess(3)}; - std::vector xLow(x.size(), -100); - std::vector xUpper(x.size(), 100); - 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 rhoend = rhobeg * 1e-6; - const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2; - std::vector w(wSize); - bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(), - rhobeg, rhoend, 100, w.data()); - - std::cout << "Final objective value:" << objective(n, x.data()) - << std::endl; - - Eigen::VectorXd eigenX(x.size(), 1); - - for (size_t xi = 0; xi < x.size(); xi++) { - eigenX(xi) = x[xi]; - } - - SimulationResults reducedModelOptimizedResults = - simulator.executeSimulation(gReducedPatternSimulationJob); - reducedModelOptimizedResults.simulationLabel = "reducedModel"; - reducedModelOptimizedResults.registerForDrawing( - gReducedPatternSimulationJob); - polyscope::show(); - - return eigenX; + 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.simulationLabel = "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]; + } + return eigenX; } std::vector ReducedModelOptimizer::createScenarios( @@ -483,7 +490,10 @@ std::vector ReducedModelOptimizer::createScenarios( std::vector scenarios; std::unordered_map> fixedVertices; std::unordered_map nodalForces; - const double forceMagnitude = 250; + std::unordered_map nodalForcedNormals; + const double forceMagnitude = 1; + // Assuming the patterns lays on the x-y plane + const CoordType patternPlaneNormal(0, 0, 1); // // Axial // for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { @@ -493,82 +503,100 @@ std::vector ReducedModelOptimizer::createScenarios( // nodalForces[viPair.first] = Vector6d({forceDirection[0], // forceDirection[1], // forceDirection[2], 0, 0, 0}) * - // forceMagnitude; + // forceMagnitude * 10; // fixedVertices[viPair.second] = // std::unordered_set{0, 1, 2, 3, 4, 5}; // } // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); // // In-plane Bending - // // Assuming the patterns lay on the x-y plane - // const CoordType patternPlane(0, 0, 1); // fixedVertices.clear(); // nodalForces.clear(); // for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { // CoordType v = // (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) // .Normalize(); - // CoordType forceDirection = patternPlane ^ v; + // CoordType forceDirection = (v ^ patternPlaneNormal).Normalize(); // nodalForces[viPair.first] = Vector6d({forceDirection[0], // forceDirection[1], // forceDirection[2], 0, 0, 0}) * - // forceMagnitude; + // 0.40 * forceMagnitude; // fixedVertices[viPair.second] = // std::unordered_set{0, 1, 2, 3, 4, 5}; // } // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); // // Torsion - // { - // fixedVertices.clear(); - // nodalForces.clear(); - // for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); - // viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) - // { - // const auto &viPair = *viPairIt; - // if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) { - // CoordType v = - // (pMesh->vert[viPair.first].cP() - - // pMesh->vert[viPair.second].cP()) - // .Normalize(); - // nodalForces[viPair.first] = - // Vector6d({0, 0, 0, v[1], v[0], 0}) * forceMagnitude; - // } else { - // fixedVertices[viPair.first] = - // std::unordered_set{0, 1, 2, 3, 4, 5}; - // } - // fixedVertices[viPair.second] = - // std::unordered_set{0, 1, 2, 3, 4, 5}; - // } - // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); - // } - // Out-of-plane bending . Pull towards Z - fixedVertices.clear(); - nodalForces.clear(); - for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { - nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}); - fixedVertices[viPair.second] = - std::unordered_set{0, 1, 2, 3, 4, 5}; - } - scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); - - // // Dou?? // fixedVertices.clear(); // nodalForces.clear(); + // for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); + // viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { + // const auto &viPair = *viPairIt; + // if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) { + // CoordType v = + // (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) + // .Normalize(); + // CoordType normalVDerivativeDir = (v ^ patternPlaneNormal).Normalize(); + // fixedVertices[viPair.first] = std::unordered_set{0, 1, 2}; + // nodalForcedNormals[viPair.first] = normalVDerivativeDir; + // } else { + // fixedVertices[viPair.first] = + // std::unordered_set{0, 1, 2, 3, 4, 5}; + // } + // fixedVertices[viPair.second] = + // std::unordered_set{0, 1, 2, 3, 4, 5}; + // } + // scenarios.push_back( + // {pMesh, fixedVertices, nodalForces, {}, nodalForcedNormals}); + + // // Out - of - plane bending.Pull towards Z + // fixedVertices.clear(); + // nodalForces.clear(); + // for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { + // nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) * + // 2; fixedVertices[viPair.second] = + // std::unordered_set{0, 1, 2, 3, 4, 5}; + // } + // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + + // Double + // fixedVertices.clear(); + // nodalForces.clear(); + // nodalForcedNormals.clear(); + // int counter = 0; // for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { // CoordType v = // (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) // .Normalize(); - // CoordType momentDirection = patternPlane ^ v; - // nodalForces[viPair.first] = - // Vector6d({0, 0, 0, momentDirection[0], momentDirection[1], 0}) * - // forceMagnitude; - // fixedVertices[viPair.first] = std::unordered_set{2}; + // nodalForcedNormals[viPair.first] = -v; + // if (counter++ == 1) { + // fixedVertices[viPair.first] = std::unordered_set{0}; + + // } else { + // fixedVertices[viPair.first] = std::unordered_set{0, 1, 2}; + // } // fixedVertices[viPair.second] = std::unordered_set{0, 1, 2}; // } - // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + // scenarios.push_back( + // {pMesh, fixedVertices, nodalForces, {}, nodalForcedNormals}); - // // Saddle + // Double + fixedVertices.clear(); + nodalForces.clear(); + for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { + CoordType v = + (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) + .Normalize(); + CoordType momentDirection = patternPlaneNormal ^ v; + nodalForces[viPair.first] = + Vector6d({0, 0, 0, momentDirection[1], momentDirection[0], 0}) * + forceMagnitude; + fixedVertices[viPair.first] = std::unordered_set{2}; + fixedVertices[viPair.second] = std::unordered_set{0, 1, 2}; + } + scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + + //// Saddle using moments // fixedVertices.clear(); // nodalForces.clear(); // for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); @@ -577,38 +605,50 @@ std::vector ReducedModelOptimizer::createScenarios( // CoordType v = // (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) // .Normalize(); - // CoordType momentDirection = patternPlane ^ v; + // // CoordType momentDirection = (patternPlane ^ v).Normalize(); // if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) { // nodalForces[viPair.first] = - // Vector6d({0, 0, 0, momentDirection[0], momentDirection[1], 0}) * 3 - // * forceMagnitude; + // Vector6d({0, 0, 0, v[0], v[1], 0}) * 5 * forceMagnitude; // nodalForces[viPair.second] = - // Vector6d({0, 0, 0, momentDirection[0], momentDirection[1], 0}) * 3 - // * - // (-forceMagnitude); + // Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 5 * forceMagnitude; // } else { // fixedVertices[viPair.first] = std::unordered_set{2}; + // // nodalForces[viPair.first] = + // // Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 2 * forceMagnitude; + // fixedVertices[viPair.second] = std::unordered_set{2}; + // nodalForces[viPair.first] = - // Vector6d({0, 0, 0, momentDirection[0], momentDirection[1], 0}) * - // (-forceMagnitude); - // fixedVertices[viPair.second] = std::unordered_set{0, 1, 2}; + // Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 2 * forceMagnitude; + // nodalForces[viPair.second] = + // Vector6d({0, 0, 0, v[0], v[1], 0}) * 2 * forceMagnitude; // } // } // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); - // std::unordered_map> - // saddleFixedVertices; - // // saddle_fixedVertices[3] = std::unordered_set{0, 1, 2}; - // saddleFixedVertices[7] = std::unordered_set{0, 1, 2}; - // saddleFixedVertices[11] = std::unordered_set{0, 1, 2}; - // // saddle_fixedVertices[15] = std::unordered_set{0, 1, 2}; - // // saddle_fixedVertices[19] = std::unordered_set{0, 1, 2}; - // // saddle_fixedVertices[23] = std::unordered_set{0, 1, 2}; - // std::unordered_map saddleNodalForces{ - // {15, {0, 0, 0, 0, -4 * 90, 0}}, {3, {0, 0, 0, 0, 4 * 90, 0}}, - // {7, {0, 0, 0, 4 * 70, 0, 0}}, {11, {0, 0, 0, 4 * 70, 0, 0}}, - // {19, {0, 0, 0, -4 * 70, 0, 0}}, {23, {0, 0, 0, -4 * 70, 0, 0}}}; - // scenarios.push_back({pMesh, saddleFixedVertices, saddleNodalForces, {}}); + //// Saddle using forced normals + fixedVertices.clear(); + nodalForces.clear(); + nodalForcedNormals.clear(); + for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); + viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { + const auto &viPair = *viPairIt; + CoordType v = + (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) + .Normalize(); + // CoordType momentDirection = (patternPlane ^ v).Normalize(); + if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) { + nodalForcedNormals[viPair.first] = v; + nodalForcedNormals[viPair.second] = -v; + } else { + + nodalForcedNormals[viPair.first] = -v; + nodalForcedNormals[viPair.second] = v; + } + fixedVertices[viPair.first] = std::unordered_set{0, 1, 2}; + fixedVertices[viPair.second] = std::unordered_set{0, 1, 2}; + } + scenarios.push_back( + {pMesh, fixedVertices, nodalForces, {}, nodalForcedNormals}); return scenarios; } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index c29d5e3..ec89164 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -4,6 +4,7 @@ #include "beamformfinder.hpp" #include "edgemesh.hpp" #include "elementalmesh.hpp" +#include "matplot/matplot.h" #include using FullModelVertexIndex = VertexIndex; @@ -47,6 +48,8 @@ private: void createSimulationMeshes(FlatPattern &fullModel, FlatPattern &reducedModel); void initializeStiffnesses(); + + static double objective(long n, const double *x); }; #endif // REDUCEDMODELOPTIMIZER_HPP