manolas
/
ReducedModelOptimization
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Optimizes for b and h but some simulations explode.

This commit is contained in:
Iason 2021-01-17 13:46:33 +02:00
parent 50c110c08e
commit 5eab2511e9
3 changed files with 107 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/main.cpp
View File

@ -121,7 +121,7 @@ int main(int argc, char *argv[]) {
// cp.copy(*reducedModels[0]); // cp.copy(*reducedModels[0]);
optimizer.initializePatterns(pattern, *reducedModels[0], optimizer.initializePatterns(pattern, *reducedModels[0],
optimizationExcludedEi); optimizationExcludedEi);
optimizer.optimize(); optimizer.optimize({ReducedModelOptimizer::SimulationScenario::Bending});
// } // }
} }

149
src/reducedmodeloptimizer.cpp
View File

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

3
src/reducedmodeloptimizer.hpp
View File

@ -57,6 +57,7 @@ public:
std::vector<double> &x); std::vector<double> &x);
static double objective(double x0, double x1, double x2, double x3); static double objective(double x0, double x1, double x2, double x3);
static double objective(double b, double h);
private: private:
void void
@ -76,7 +77,7 @@ private:
void createSimulationMeshes(FlatPattern &fullModel, void createSimulationMeshes(FlatPattern &fullModel,
FlatPattern &reducedModel); FlatPattern &reducedModel);
static void static void
initializeStiffnesses(const std::shared_ptr<SimulationMesh> &mesh); initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh);
static double static double
computeError(const SimulationResults &reducedPatternResults, computeError(const SimulationResults &reducedPatternResults,