2020-11-27 11:45:20 +01:00
|
|
|
|
#include "reducedmodeloptimizer.hpp"
|
|
|
|
|
#include "flatpattern.hpp"
|
2020-11-23 10:06:45 +01:00
|
|
|
|
#include "gradientDescent.h"
|
|
|
|
|
#include "simulationhistoryplotter.hpp"
|
2020-11-27 11:45:20 +01:00
|
|
|
|
#include "trianglepattterntopology.hpp"
|
2021-01-12 13:41:40 +01:00
|
|
|
|
#include <chrono>
|
|
|
|
|
#include <dlib/global_optimization.h>
|
|
|
|
|
#include <dlib/optimization.h>
|
2020-11-23 10:06:45 +01:00
|
|
|
|
|
2020-11-27 11:45:20 +01:00
|
|
|
|
const bool gShouldDraw = true;
|
2020-12-09 16:58:48 +01:00
|
|
|
|
|
2020-11-27 11:45:20 +01:00
|
|
|
|
FormFinder simulator;
|
2020-12-21 16:56:21 +01:00
|
|
|
|
std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
|
2021-01-04 13:12:25 +01:00
|
|
|
|
std::vector<SimulationJob> g_fullPatternSimulationJob;
|
|
|
|
|
std::vector<std::shared_ptr<SimulationJob>> g_reducedPatternSimulationJob;
|
2020-12-14 10:07:43 +01:00
|
|
|
|
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
|
2021-01-04 13:12:25 +01:00
|
|
|
|
reducedToFullInterfaceViMap;
|
|
|
|
|
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
|
|
|
|
|
g_reducedToFullViMap;
|
2020-11-27 11:45:20 +01:00
|
|
|
|
matplot::line_handle gPlotHandle;
|
|
|
|
|
std::vector<double> gObjectiveValueHistory;
|
2021-01-17 12:46:33 +01:00
|
|
|
|
Eigen::Vector2d g_initialX;
|
2020-11-27 11:45:20 +01:00
|
|
|
|
std::unordered_set<size_t> g_reducedPatternExludedEdges;
|
2020-12-16 20:31:58 +01:00
|
|
|
|
// double g_initialParameters;
|
|
|
|
|
Eigen::VectorXd g_initialParameters;
|
2021-01-04 13:12:25 +01:00
|
|
|
|
std::vector<ReducedModelOptimizer::SimulationScenario>
|
|
|
|
|
g_simulationScenarioIndices;
|
2020-12-22 17:47:29 +01:00
|
|
|
|
std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)};
|
|
|
|
|
double g_innerHexagonInitialPos = 0;
|
2021-01-04 13:12:25 +01:00
|
|
|
|
bool g_optimizeInnerHexagonSize{false};
|
|
|
|
|
std::vector<SimulationResults> firstOptimizationRoundResults;
|
|
|
|
|
int g_firstRoundIterationIndex{0};
|
2021-01-12 13:41:40 +01:00
|
|
|
|
double minY{std::numeric_limits<double>::max()};
|
|
|
|
|
std::vector<double> minX;
|
2021-01-22 15:39:36 +01:00
|
|
|
|
std::vector<std::vector<double>> failedSimulationsXRatio;
|
|
|
|
|
int numOfSimulationCrashes{false};
|
2021-01-31 22:36:46 +01:00
|
|
|
|
int numberOfFunctionCalls{ 0 };
|
2020-12-09 16:58:48 +01:00
|
|
|
|
|
|
|
|
|
// 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<SimulationMesh> 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;
|
|
|
|
|
// }
|
|
|
|
|
//};
|
|
|
|
|
|
2020-12-21 16:56:21 +01:00
|
|
|
|
double ReducedModelOptimizer::computeError(
|
|
|
|
|
const SimulationResults &reducedPatternResults,
|
|
|
|
|
const Eigen::MatrixX3d &optimalReducedPatternDisplacements) {
|
|
|
|
|
double error = 0;
|
2021-01-04 13:12:25 +01:00
|
|
|
|
for (const auto reducedFullViPair : g_reducedToFullViMap) {
|
2020-12-21 16:56:21 +01:00
|
|
|
|
VertexIndex reducedModelVi = reducedFullViPair.first;
|
|
|
|
|
// const auto pos =
|
|
|
|
|
// g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP();
|
|
|
|
|
// std::cout << "Interface vi " << reducedModelVi << " is at position "
|
|
|
|
|
// << pos[0] << " " << pos[1] << " " << pos[2] << std::endl;
|
|
|
|
|
Eigen::Vector3d vertexDisplacement(
|
|
|
|
|
reducedPatternResults.displacements[reducedModelVi][0],
|
|
|
|
|
reducedPatternResults.displacements[reducedModelVi][1],
|
|
|
|
|
reducedPatternResults.displacements[reducedModelVi][2]);
|
2021-01-12 13:41:40 +01:00
|
|
|
|
if (!std::isfinite(vertexDisplacement[0]) ||
|
|
|
|
|
!std::isfinite(vertexDisplacement[1]) ||
|
|
|
|
|
!std::isfinite(vertexDisplacement[2])) {
|
|
|
|
|
return std::numeric_limits<double>::max();
|
|
|
|
|
}
|
2020-12-21 16:56:21 +01:00
|
|
|
|
Eigen::Vector3d errorVector =
|
|
|
|
|
Eigen::Vector3d(
|
|
|
|
|
optimalReducedPatternDisplacements.row(reducedModelVi)) -
|
|
|
|
|
vertexDisplacement;
|
|
|
|
|
// error += errorVector.squaredNorm();
|
|
|
|
|
error += errorVector.norm();
|
2020-11-23 10:06:45 +01:00
|
|
|
|
}
|
2020-12-21 16:56:21 +01:00
|
|
|
|
return error;
|
|
|
|
|
}
|
2020-11-23 10:06:45 +01:00
|
|
|
|
|
2020-12-21 16:56:21 +01:00
|
|
|
|
void updateMesh(long n, const double *x) {
|
2021-01-04 13:12:25 +01:00
|
|
|
|
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
|
2021-01-17 12:46:33 +01:00
|
|
|
|
g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]->pMesh;
|
2021-01-04 13:12:25 +01:00
|
|
|
|
// const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0];
|
|
|
|
|
// std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << "
|
|
|
|
|
// "
|
|
|
|
|
// << elem.firstBendingConstFactor << " "
|
|
|
|
|
// << elem.secondBendingConstFactor << std::endl;
|
2020-12-21 16:56:21 +01:00
|
|
|
|
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
|
|
|
|
|
Element &e = pReducedPatternSimulationMesh->elements[ei];
|
2020-12-09 16:58:48 +01:00
|
|
|
|
// if (g_reducedPatternExludedEdges.contains(ei)) {
|
|
|
|
|
// continue;
|
|
|
|
|
// }
|
2020-12-16 20:31:58 +01:00
|
|
|
|
// e.properties.E = g_initialParameters * x[ei];
|
2021-01-04 13:12:25 +01:00
|
|
|
|
// e.properties.E = g_initialParameters(0) * x[0];
|
|
|
|
|
// e.properties.G = g_initialParameters(1) * x[1];
|
2021-01-22 15:39:36 +01:00
|
|
|
|
e.setDimensions(RectangularBeamDimensions(
|
|
|
|
|
g_initialParameters(0) * x[0],
|
|
|
|
|
g_initialParameters(0) * x[0] / (g_initialParameters(1) * x[1])));
|
|
|
|
|
e.setMaterial(ElementMaterial(e.material.poissonsRatio,
|
|
|
|
|
g_initialParameters(2) * x[2]));
|
2021-01-17 12:46:33 +01:00
|
|
|
|
// e.properties.A = g_initialParameters(0) * x[0];
|
|
|
|
|
// e.properties.J = g_initialParameters(1) * x[1];
|
|
|
|
|
// e.properties.I2 = g_initialParameters(2) * x[2];
|
|
|
|
|
// e.properties.I3 = g_initialParameters(3) * x[3];
|
2020-12-16 20:31:58 +01:00
|
|
|
|
// e.properties.G = e.properties.E / (2 * (1 + 0.3));
|
2021-01-17 12:46:33 +01:00
|
|
|
|
// 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;
|
2020-11-23 10:06:45 +01:00
|
|
|
|
}
|
2020-12-22 17:47:29 +01:00
|
|
|
|
|
2021-01-04 13:12:25 +01:00
|
|
|
|
// 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) {
|
2020-12-22 17:47:29 +01:00
|
|
|
|
assert(pReducedPatternSimulationMesh->EN() == 12);
|
|
|
|
|
for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN();
|
|
|
|
|
vi += 2) {
|
|
|
|
|
pReducedPatternSimulationMesh->vert[vi].P() =
|
2021-01-04 13:12:25 +01:00
|
|
|
|
g_innerHexagonVectors[vi / 2] * x[n - 1];
|
2020-12-22 17:47:29 +01:00
|
|
|
|
}
|
|
|
|
|
|
2020-12-22 19:21:28 +01:00
|
|
|
|
pReducedPatternSimulationMesh->reset();
|
2021-01-04 13:12:25 +01:00
|
|
|
|
pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
|
2020-12-22 17:47:29 +01:00
|
|
|
|
// pReducedPatternSimulationMesh->registerForDrawing("Optimized
|
|
|
|
|
// hexagon"); polyscope::show();
|
|
|
|
|
}
|
2020-12-21 16:56:21 +01:00
|
|
|
|
}
|
2020-11-27 11:45:20 +01:00
|
|
|
|
|
2021-01-22 15:39:36 +01:00
|
|
|
|
double ReducedModelOptimizer::objective(double b, double h, double E) {
|
|
|
|
|
std::vector<double> x{b, h, E};
|
2021-01-17 12:46:33 +01:00
|
|
|
|
return ReducedModelOptimizer::objective(x.size(), x.data());
|
|
|
|
|
}
|
|
|
|
|
|
2021-01-12 13:41:40 +01:00
|
|
|
|
double ReducedModelOptimizer::objective(double x0, double x1, double x2,
|
|
|
|
|
double x3) {
|
|
|
|
|
std::vector<double> x{x0, x1, x2, x3};
|
2021-01-17 12:46:33 +01:00
|
|
|
|
return ReducedModelOptimizer::objective(x.size(), x.data());
|
2021-01-12 13:41:40 +01:00
|
|
|
|
}
|
|
|
|
|
|
2021-01-04 13:12:25 +01:00
|
|
|
|
double ReducedModelOptimizer::objective(long n, const double *x) {
|
2021-01-22 15:39:36 +01:00
|
|
|
|
// std::cout.precision(17);
|
2020-12-21 16:56:21 +01:00
|
|
|
|
|
2021-01-22 15:39:36 +01:00
|
|
|
|
// for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
|
|
|
|
|
// std::cout << "x[" + std::to_string(parameterIndex) + "]="
|
|
|
|
|
// << x[parameterIndex] << std::endl;
|
|
|
|
|
// }
|
2021-01-04 13:12:25 +01:00
|
|
|
|
const Element &e = g_reducedPatternSimulationJob[0]->pMesh->elements[0];
|
|
|
|
|
// std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
|
|
|
|
|
// << e.firstBendingConstFactor << " " <<
|
|
|
|
|
// e.secondBendingConstFactor
|
|
|
|
|
// << std::endl;
|
2020-12-21 16:56:21 +01:00
|
|
|
|
updateMesh(n, x);
|
2021-01-17 12:46:33 +01:00
|
|
|
|
// std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
|
|
|
|
|
// << e.firstBendingConstFactor << " " <<
|
|
|
|
|
// e.secondBendingConstFactor
|
|
|
|
|
// << std::endl;
|
2020-12-21 16:56:21 +01:00
|
|
|
|
|
|
|
|
|
// run simulations
|
2020-11-23 10:06:45 +01:00
|
|
|
|
double error = 0;
|
2021-01-17 12:46:33 +01:00
|
|
|
|
FormFinder::Settings simulationSettings;
|
|
|
|
|
// simulationSettings.shouldDraw = true;
|
2021-01-04 13:12:25 +01:00
|
|
|
|
for (const int simulationScenarioIndex : g_simulationScenarioIndices) {
|
2020-12-21 16:56:21 +01:00
|
|
|
|
SimulationResults reducedModelResults = simulator.executeSimulation(
|
2021-01-17 12:46:33 +01:00
|
|
|
|
g_reducedPatternSimulationJob[simulationScenarioIndex],
|
|
|
|
|
simulationSettings);
|
|
|
|
|
std::string filename;
|
2021-01-22 15:39:36 +01:00
|
|
|
|
if (!reducedModelResults.converged /*&&
|
|
|
|
|
g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]
|
|
|
|
|
->pMesh->elements[0]
|
|
|
|
|
.A > 1e-8 &
|
|
|
|
|
x[0] / x[1] < 60*/) {
|
2021-01-17 12:46:33 +01:00
|
|
|
|
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";
|
2021-01-22 15:39:36 +01:00
|
|
|
|
// if (failedSimulationsXRatio.empty()) {
|
|
|
|
|
// failedSimulationsXRatio.resize(2);
|
|
|
|
|
// }
|
|
|
|
|
// failedSimulationsXRatio[0].push_back(std::log(x[0] / x[1]));
|
|
|
|
|
// failedSimulationsXRatio[1].push_back(
|
|
|
|
|
// std::log(g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]
|
|
|
|
|
// ->pMesh->elements[0]
|
|
|
|
|
// .A));
|
|
|
|
|
|
|
|
|
|
// SimulationResultsReporter::createPlot(
|
|
|
|
|
// "log(b/h)", "log(A)", failedSimulationsXRatio[0],
|
|
|
|
|
// failedSimulationsXRatio[1], "ratioToAPlot.png");
|
2021-01-17 12:46:33 +01:00
|
|
|
|
// std::cout << "Failed simulation" << std::endl;
|
|
|
|
|
// simulationSettings.shouldDraw = true;
|
|
|
|
|
// simulationSettings.debugMessages = true;
|
|
|
|
|
// simulator.executeSimulation(
|
|
|
|
|
// g_reducedPatternSimulationJob[simulationScenarioIndex],
|
|
|
|
|
// simulationSettings);
|
2021-01-22 15:39:36 +01:00
|
|
|
|
numOfSimulationCrashes++;
|
2021-01-17 12:46:33 +01:00
|
|
|
|
} 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);
|
2021-01-22 15:39:36 +01:00
|
|
|
|
auto pMesh =
|
|
|
|
|
g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]->pMesh;
|
|
|
|
|
|
|
|
|
|
for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
|
|
|
|
|
out << "x[" + std::to_string(parameterIndex) + "]=" << x[parameterIndex]
|
|
|
|
|
<< std::endl;
|
|
|
|
|
}
|
|
|
|
|
out << pMesh->elements[0].dimensions.toString() + "\n" +
|
|
|
|
|
pMesh->elements[0].material.toString() + " \nA="
|
|
|
|
|
<< pMesh->elements[0].A << " \nratio="
|
|
|
|
|
<< pMesh->elements[0].dimensions.b / pMesh->elements[0].dimensions.h
|
|
|
|
|
<< " \naxialRig:" << pMesh->elements[0].rigidity.axial
|
|
|
|
|
<< " \ntorsionalRig:" << pMesh->elements[0].rigidity.torsional
|
|
|
|
|
<< " \nfirstBendingRig:" << pMesh->elements[0].rigidity.firstBending
|
|
|
|
|
<< " \nsecondBendingRig:" << pMesh->elements[0].rigidity.secondBending
|
|
|
|
|
<< " \nscenario:" + simulationScenarioStrings[simulationScenarioIndex] +
|
|
|
|
|
"\n\n";
|
2021-01-17 12:46:33 +01:00
|
|
|
|
out.close();
|
2020-12-21 16:56:21 +01:00
|
|
|
|
}
|
2021-01-12 13:41:40 +01:00
|
|
|
|
if (error < minY) {
|
|
|
|
|
minY = error;
|
2021-01-17 12:46:33 +01:00
|
|
|
|
minX.assign(x, x + n);
|
2021-01-12 13:41:40 +01:00
|
|
|
|
}
|
2021-01-31 22:36:46 +01:00
|
|
|
|
if(++numberOfFunctionCalls%50==0){
|
|
|
|
|
std::cout << "Number of function calls:"<<numberOfFunctionCalls << std::endl;
|
|
|
|
|
}
|
2020-12-21 16:56:21 +01:00
|
|
|
|
|
|
|
|
|
// compute error and return it
|
|
|
|
|
gObjectiveValueHistory.push_back(error);
|
2021-01-12 13:41:40 +01:00
|
|
|
|
// auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
|
|
|
|
|
// gObjectiveValueHistory.size());
|
|
|
|
|
// std::vector<double> colors(gObjectiveValueHistory.size(), 2);
|
|
|
|
|
// if (g_firstRoundIterationIndex != 0) {
|
|
|
|
|
// for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(),
|
|
|
|
|
// [](double &c) { c = 0.7; });
|
|
|
|
|
// }
|
|
|
|
|
// gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory, 6, colors);
|
2021-01-17 12:46:33 +01:00
|
|
|
|
// SimulationResultsReporter::createPlot("Number of Steps", "Objective
|
|
|
|
|
// value",
|
|
|
|
|
// gObjectiveValueHistory);
|
2020-11-27 11:45:20 +01:00
|
|
|
|
|
2020-11-23 10:06:45 +01:00
|
|
|
|
return error;
|
|
|
|
|
}
|
|
|
|
|
|
2020-11-27 11:45:20 +01:00
|
|
|
|
void ReducedModelOptimizer::computeMaps(
|
|
|
|
|
FlatPattern &fullPattern, FlatPattern &reducedPattern,
|
|
|
|
|
const std::unordered_set<size_t> &reducedModelExcludedEges) {
|
|
|
|
|
// Compute the offset between the interface nodes
|
|
|
|
|
const size_t interfaceSlotIndex = 4; // bottom edge
|
|
|
|
|
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());
|
|
|
|
|
|
|
|
|
|
vcg::tri::Allocator<FlatPattern>::PointerUpdater<FlatPattern::VertexPointer>
|
|
|
|
|
pu_fullModel;
|
|
|
|
|
fullPattern.deleteDanglingVertices(pu_fullModel);
|
|
|
|
|
const size_t fullModelBaseTriangleInterfaceVi =
|
|
|
|
|
pu_fullModel.remap.empty() ? baseTriangleInterfaceVi
|
|
|
|
|
: pu_fullModel.remap[baseTriangleInterfaceVi];
|
|
|
|
|
const size_t fullModelBaseTriangleVN = fullPattern.VN();
|
|
|
|
|
fullPattern.createFan();
|
|
|
|
|
const size_t duplicateVerticesPerFanPair =
|
|
|
|
|
fullModelBaseTriangleVN - fullPattern.VN() / 6;
|
|
|
|
|
const size_t fullPatternInterfaceVertexOffset =
|
|
|
|
|
fullModelBaseTriangleVN - duplicateVerticesPerFanPair;
|
|
|
|
|
// std::cout << "Dups in fan pair:" << duplicateVerticesPerFanPair <<
|
|
|
|
|
// std::endl;
|
|
|
|
|
|
|
|
|
|
// Save excluded edges
|
|
|
|
|
g_reducedPatternExludedEdges.clear();
|
|
|
|
|
const size_t fanSize = 6;
|
|
|
|
|
const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN();
|
|
|
|
|
for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
|
|
|
|
|
for (const size_t ei : reducedModelExcludedEges) {
|
|
|
|
|
g_reducedPatternExludedEdges.insert(
|
|
|
|
|
fanIndex * reducedBaseTriangleNumberOfEdges + ei);
|
|
|
|
|
}
|
2020-11-23 10:06:45 +01:00
|
|
|
|
}
|
|
|
|
|
|
2020-11-27 11:45:20 +01:00
|
|
|
|
// Construct reduced->full and full->reduced interface vi map
|
2021-01-04 13:12:25 +01:00
|
|
|
|
reducedToFullInterfaceViMap.clear();
|
2020-11-27 11:45:20 +01:00
|
|
|
|
vcg::tri::Allocator<FlatPattern>::PointerUpdater<FlatPattern::VertexPointer>
|
|
|
|
|
pu_reducedModel;
|
|
|
|
|
reducedPattern.deleteDanglingVertices(pu_reducedModel);
|
|
|
|
|
const size_t reducedModelBaseTriangleInterfaceVi =
|
|
|
|
|
pu_reducedModel.remap[baseTriangleInterfaceVi];
|
|
|
|
|
const size_t reducedModelInterfaceVertexOffset =
|
|
|
|
|
reducedPattern.VN() - 1 /*- reducedModelBaseTriangleInterfaceVi*/;
|
|
|
|
|
reducedPattern.createFan();
|
|
|
|
|
for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
|
2021-01-04 13:12:25 +01:00
|
|
|
|
reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex +
|
|
|
|
|
reducedModelBaseTriangleInterfaceVi] =
|
2020-11-27 11:45:20 +01:00
|
|
|
|
fullModelBaseTriangleInterfaceVi +
|
|
|
|
|
fanIndex * fullPatternInterfaceVertexOffset;
|
|
|
|
|
}
|
|
|
|
|
m_fullToReducedInterfaceViMap.clear();
|
2021-01-04 13:12:25 +01:00
|
|
|
|
constructInverseMap(reducedToFullInterfaceViMap,
|
2020-11-27 11:45:20 +01:00
|
|
|
|
m_fullToReducedInterfaceViMap);
|
|
|
|
|
|
|
|
|
|
// fullPattern.setLabel("FullPattern");
|
|
|
|
|
// reducedPattern.setLabel("ReducedPattern");
|
|
|
|
|
// Create opposite vertex map
|
|
|
|
|
m_fullPatternOppositeInterfaceViMap.clear();
|
2020-12-16 20:31:58 +01:00
|
|
|
|
for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) {
|
2020-11-27 11:45:20 +01:00
|
|
|
|
const size_t vi0 = fullModelBaseTriangleInterfaceVi +
|
|
|
|
|
fanIndex * fullPatternInterfaceVertexOffset;
|
|
|
|
|
const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset;
|
|
|
|
|
assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN());
|
|
|
|
|
m_fullPatternOppositeInterfaceViMap[vi0] = vi1;
|
|
|
|
|
}
|
|
|
|
|
|
2021-01-04 13:12:25 +01:00
|
|
|
|
g_reducedToFullViMap = reducedToFullInterfaceViMap;
|
|
|
|
|
|
2020-11-27 11:45:20 +01:00
|
|
|
|
const bool debugMapping = false;
|
|
|
|
|
if (debugMapping) {
|
|
|
|
|
reducedPattern.registerForDrawing();
|
|
|
|
|
std::vector<glm::vec3> colors_reducedPatternExcludedEdges(
|
|
|
|
|
reducedPattern.EN(), glm::vec3(0, 0, 0));
|
|
|
|
|
for (const size_t ei : g_reducedPatternExludedEdges) {
|
|
|
|
|
colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0);
|
|
|
|
|
}
|
|
|
|
|
const std::string label = reducedPattern.getLabel();
|
|
|
|
|
polyscope::getCurveNetwork(label)
|
|
|
|
|
->addEdgeColorQuantity("Excluded edges",
|
|
|
|
|
colors_reducedPatternExcludedEdges)
|
|
|
|
|
->setEnabled(true);
|
|
|
|
|
polyscope::show();
|
|
|
|
|
|
|
|
|
|
std::vector<glm::vec3> nodeColorsOpposite(fullPattern.VN(),
|
|
|
|
|
glm::vec3(0, 0, 0));
|
|
|
|
|
for (const std::pair<size_t, size_t> oppositeVerts :
|
|
|
|
|
m_fullPatternOppositeInterfaceViMap) {
|
|
|
|
|
auto color = polyscope::getNextUniqueColor();
|
|
|
|
|
nodeColorsOpposite[oppositeVerts.first] = color;
|
|
|
|
|
nodeColorsOpposite[oppositeVerts.second] = color;
|
|
|
|
|
}
|
|
|
|
|
fullPattern.registerForDrawing();
|
|
|
|
|
polyscope::getCurveNetwork(fullPattern.getLabel())
|
|
|
|
|
->addNodeColorQuantity("oppositeMap", nodeColorsOpposite)
|
|
|
|
|
->setEnabled(true);
|
|
|
|
|
polyscope::show();
|
|
|
|
|
|
|
|
|
|
std::vector<glm::vec3> nodeColorsReducedToFull_reduced(reducedPattern.VN(),
|
|
|
|
|
glm::vec3(0, 0, 0));
|
|
|
|
|
std::vector<glm::vec3> nodeColorsReducedToFull_full(fullPattern.VN(),
|
|
|
|
|
glm::vec3(0, 0, 0));
|
|
|
|
|
for (size_t vi = 0; vi < reducedPattern.VN(); vi++) {
|
2021-01-04 13:12:25 +01:00
|
|
|
|
if (reducedToFullInterfaceViMap.contains(vi)) {
|
2020-11-27 11:45:20 +01:00
|
|
|
|
|
|
|
|
|
auto color = polyscope::getNextUniqueColor();
|
|
|
|
|
nodeColorsReducedToFull_reduced[vi] = color;
|
2021-01-04 13:12:25 +01:00
|
|
|
|
nodeColorsReducedToFull_full[reducedToFullInterfaceViMap[vi]] = color;
|
2020-11-27 11:45:20 +01:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
polyscope::getCurveNetwork(reducedPattern.getLabel())
|
|
|
|
|
->addNodeColorQuantity("reducedToFull_reduced",
|
|
|
|
|
nodeColorsReducedToFull_reduced)
|
|
|
|
|
->setEnabled(true);
|
|
|
|
|
polyscope::getCurveNetwork(fullPattern.getLabel())
|
|
|
|
|
->addNodeColorQuantity("reducedToFull_full",
|
|
|
|
|
nodeColorsReducedToFull_full)
|
|
|
|
|
->setEnabled(true);
|
|
|
|
|
polyscope::show();
|
2020-11-23 10:06:45 +01:00
|
|
|
|
}
|
2020-11-27 11:45:20 +01:00
|
|
|
|
}
|
2020-11-23 10:06:45 +01:00
|
|
|
|
|
2020-11-27 11:45:20 +01:00
|
|
|
|
void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel,
|
|
|
|
|
FlatPattern &reducedModel) {
|
2020-12-16 20:31:58 +01:00
|
|
|
|
if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) {
|
2020-12-17 21:33:42 +01:00
|
|
|
|
std::cerr << "Error: A rectangular cross section is expected." << std::endl;
|
|
|
|
|
terminate();
|
2020-12-16 20:31:58 +01:00
|
|
|
|
}
|
|
|
|
|
m_pReducedPatternSimulationMesh =
|
2020-12-14 10:07:43 +01:00
|
|
|
|
std::make_shared<SimulationMesh>(reducedModel);
|
2020-12-16 20:31:58 +01:00
|
|
|
|
m_pReducedPatternSimulationMesh->setBeamCrossSection(
|
2020-12-09 16:58:48 +01:00
|
|
|
|
CrossSectionType{0.002, 0.002});
|
2021-01-17 12:46:33 +01:00
|
|
|
|
m_pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
m_pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
|
|
|
|
|
m_pFullPatternSimulationMesh->setBeamCrossSection(
|
2020-12-14 10:07:43 +01:00
|
|
|
|
CrossSectionType{0.002, 0.002});
|
2021-01-17 12:46:33 +01:00
|
|
|
|
m_pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
|
2020-11-27 11:45:20 +01:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
ReducedModelOptimizer::ReducedModelOptimizer(
|
|
|
|
|
const std::vector<size_t> &numberOfNodesPerSlot) {
|
|
|
|
|
FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlot);
|
|
|
|
|
FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode);
|
|
|
|
|
}
|
|
|
|
|
|
2020-12-21 16:56:21 +01:00
|
|
|
|
void ReducedModelOptimizer::initializePatterns(
|
2020-11-27 11:45:20 +01:00
|
|
|
|
FlatPattern &fullPattern, FlatPattern &reducedPattern,
|
|
|
|
|
const std::unordered_set<size_t> &reducedModelExcludedEdges) {
|
2021-01-04 13:12:25 +01:00
|
|
|
|
fullPattern.setLabel("full_pattern_" + fullPattern.getLabel());
|
|
|
|
|
reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel());
|
2020-11-27 11:45:20 +01:00
|
|
|
|
assert(fullPattern.VN() == reducedPattern.VN() &&
|
2020-12-14 10:07:43 +01:00
|
|
|
|
fullPattern.EN() >= reducedPattern.EN());
|
2020-11-27 11:45:20 +01:00
|
|
|
|
polyscope::removeAllStructures();
|
|
|
|
|
// Create copies of the input models
|
|
|
|
|
FlatPattern copyFullPattern;
|
|
|
|
|
FlatPattern copyReducedPattern;
|
|
|
|
|
copyFullPattern.copy(fullPattern);
|
|
|
|
|
copyReducedPattern.copy(reducedPattern);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
g_optimizeInnerHexagonSize = copyReducedPattern.EN() == 2;
|
|
|
|
|
if (g_optimizeInnerHexagonSize) {
|
2020-12-22 17:47:29 +01:00
|
|
|
|
const double h = copyReducedPattern.getBaseTriangleHeight();
|
|
|
|
|
double baseTriangle_bottomEdgeSize = 2 * h / tan(vcg::math::ToRad(60.0));
|
|
|
|
|
VectorType baseTriangle_leftBottomNode(-baseTriangle_bottomEdgeSize / 2, -h,
|
|
|
|
|
0);
|
|
|
|
|
|
|
|
|
|
const int fanSize = 6;
|
|
|
|
|
const CoordType rotationAxis(0, 0, 1);
|
|
|
|
|
for (int rotationCounter = 0; rotationCounter < fanSize;
|
|
|
|
|
rotationCounter++) {
|
|
|
|
|
VectorType rotatedVector =
|
|
|
|
|
vcg::RotationMatrix(rotationAxis,
|
|
|
|
|
vcg::math::ToRad(rotationCounter * 60.0)) *
|
|
|
|
|
baseTriangle_leftBottomNode;
|
|
|
|
|
g_innerHexagonVectors[rotationCounter] = rotatedVector;
|
|
|
|
|
}
|
2021-01-04 13:12:25 +01:00
|
|
|
|
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;
|
2020-12-22 17:47:29 +01:00
|
|
|
|
}
|
2020-11-27 11:45:20 +01:00
|
|
|
|
computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges);
|
|
|
|
|
createSimulationMeshes(copyFullPattern, copyReducedPattern);
|
2021-01-17 12:46:33 +01:00
|
|
|
|
initializeOptimizationParameters(m_pReducedPatternSimulationMesh);
|
2020-11-27 11:45:20 +01:00
|
|
|
|
}
|
|
|
|
|
|
2021-01-17 12:46:33 +01:00
|
|
|
|
void ReducedModelOptimizer::initializeOptimizationParameters(
|
2021-01-04 13:12:25 +01:00
|
|
|
|
const std::shared_ptr<SimulationMesh> &mesh) {
|
2021-01-22 15:39:36 +01:00
|
|
|
|
const int numberOfOptimizationParameters = 3;
|
2021-01-17 12:46:33 +01:00
|
|
|
|
g_initialParameters.resize(g_optimizeInnerHexagonSize
|
|
|
|
|
? numberOfOptimizationParameters + 1
|
|
|
|
|
: numberOfOptimizationParameters);
|
2020-11-27 11:45:20 +01:00
|
|
|
|
// Save save the beam stiffnesses
|
2020-12-09 16:58:48 +01:00
|
|
|
|
// 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;
|
|
|
|
|
// }
|
2021-01-17 12:46:33 +01:00
|
|
|
|
const double initialB = std::sqrt(mesh->elements[0].A);
|
2021-01-22 15:39:36 +01:00
|
|
|
|
const double initialRatio = 1;
|
|
|
|
|
;
|
2021-01-17 12:46:33 +01:00
|
|
|
|
g_initialParameters(0) = initialB;
|
2021-01-22 15:39:36 +01:00
|
|
|
|
g_initialParameters(1) = initialRatio;
|
|
|
|
|
g_initialParameters(2) = mesh->elements[0].material.youngsModulus;
|
2020-12-16 20:31:58 +01:00
|
|
|
|
// g_initialParameters =
|
|
|
|
|
// m_pReducedPatternSimulationMesh->elements[0].properties.E;
|
|
|
|
|
// for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) {
|
|
|
|
|
// }
|
2021-01-04 13:12:25 +01:00
|
|
|
|
// g_initialParameters(0) = mesh->elements[0].properties.E;
|
|
|
|
|
// g_initialParameters(1) = mesh->elements[0].properties.G;
|
2021-01-17 12:46:33 +01:00
|
|
|
|
// 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;
|
2020-12-09 16:58:48 +01:00
|
|
|
|
// }
|
2020-11-23 10:06:45 +01:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void ReducedModelOptimizer::computeReducedModelSimulationJob(
|
|
|
|
|
const SimulationJob &simulationJobOfFullModel,
|
2021-01-04 13:12:25 +01:00
|
|
|
|
const std::unordered_map<size_t, size_t> &simulationJobFullToReducedMap,
|
2020-11-23 10:06:45 +01:00
|
|
|
|
SimulationJob &simulationJobOfReducedModel) {
|
2021-01-04 13:12:25 +01:00
|
|
|
|
assert(simulationJobOfReducedModel.pMesh->VN() != 0);
|
2020-11-23 10:06:45 +01:00
|
|
|
|
std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
|
|
|
|
|
reducedModelFixedVertices;
|
2021-01-04 13:12:25 +01:00
|
|
|
|
for (auto fullModelFixedVertex :
|
|
|
|
|
simulationJobOfFullModel.constrainedVertices) {
|
|
|
|
|
reducedModelFixedVertices[simulationJobFullToReducedMap.at(
|
2020-11-23 10:06:45 +01:00
|
|
|
|
fullModelFixedVertex.first)] = fullModelFixedVertex.second;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
std::unordered_map<VertexIndex, Vector6d> reducedModelNodalForces;
|
|
|
|
|
for (auto fullModelNodalForce :
|
|
|
|
|
simulationJobOfFullModel.nodalExternalForces) {
|
2021-01-04 13:12:25 +01:00
|
|
|
|
reducedModelNodalForces[simulationJobFullToReducedMap.at(
|
2020-11-27 11:45:20 +01:00
|
|
|
|
fullModelNodalForce.first)] = fullModelNodalForce.second;
|
2020-11-23 10:06:45 +01:00
|
|
|
|
}
|
|
|
|
|
|
2021-01-04 13:12:25 +01:00
|
|
|
|
// std::unordered_map<VertexIndex, VectorType>
|
|
|
|
|
// 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;
|
2020-11-23 10:06:45 +01:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
|
|
|
|
|
const SimulationResults &fullModelResults,
|
2021-01-04 13:12:25 +01:00
|
|
|
|
const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
|
2020-11-27 11:45:20 +01:00
|
|
|
|
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) {
|
|
|
|
|
|
2021-01-04 13:12:25 +01:00
|
|
|
|
assert(optimalDisplacementsOfReducedModel.rows() != 0 &&
|
|
|
|
|
optimalDisplacementsOfReducedModel.cols() == 3);
|
2020-11-27 11:45:20 +01:00
|
|
|
|
optimalDisplacementsOfReducedModel.setZero(
|
|
|
|
|
optimalDisplacementsOfReducedModel.rows(),
|
|
|
|
|
optimalDisplacementsOfReducedModel.cols());
|
|
|
|
|
|
2021-01-04 13:12:25 +01:00
|
|
|
|
for (auto reducedFullViPair : displacementsReducedToFullMap) {
|
2020-11-23 10:06:45 +01:00
|
|
|
|
const VertexIndex fullModelVi = reducedFullViPair.second;
|
2020-11-27 11:45:20 +01:00
|
|
|
|
const Vector6d fullModelViDisplacements =
|
2020-11-23 10:06:45 +01:00
|
|
|
|
fullModelResults.displacements[fullModelVi];
|
2020-11-27 11:45:20 +01:00
|
|
|
|
optimalDisplacementsOfReducedModel.row(reducedFullViPair.first) =
|
2020-11-23 10:06:45 +01:00
|
|
|
|
Eigen::Vector3d(fullModelViDisplacements[0],
|
|
|
|
|
fullModelViDisplacements[1],
|
|
|
|
|
fullModelViDisplacements[2]);
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2021-01-22 15:39:36 +01:00
|
|
|
|
ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(
|
|
|
|
|
const Settings &settings,
|
2020-12-21 16:56:21 +01:00
|
|
|
|
double (*pObjectiveFunction)(long, const double *)) {
|
2020-11-27 11:45:20 +01:00
|
|
|
|
|
2020-12-14 10:07:43 +01:00
|
|
|
|
gObjectiveValueHistory.clear();
|
2020-12-21 16:56:21 +01:00
|
|
|
|
|
2021-01-17 12:46:33 +01:00
|
|
|
|
const size_t n = g_initialParameters.rows();
|
2021-01-22 15:39:36 +01:00
|
|
|
|
assert(n == 3);
|
2021-01-17 12:46:33 +01:00
|
|
|
|
// g_optimizeInnerHexagonSize ? 5: 4;
|
2021-01-12 13:41:40 +01:00
|
|
|
|
// const size_t npt = (n + 1) * (n + 2) / 2;
|
|
|
|
|
// // ((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.");
|
2020-11-23 10:06:45 +01:00
|
|
|
|
// Set initial guess of solution
|
2020-12-16 20:31:58 +01:00
|
|
|
|
|
2021-01-12 13:41:40 +01:00
|
|
|
|
const size_t initialGuess = 1;
|
|
|
|
|
std::vector<double> x(n, initialGuess);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
if (g_optimizeInnerHexagonSize) {
|
|
|
|
|
x[n - 1] = g_innerHexagonInitialPos;
|
|
|
|
|
}
|
|
|
|
|
/*if (!initialGuess.empty()) {
|
|
|
|
|
x = g_optimizationInitialGuess;
|
|
|
|
|
}*/ // {0.10000000000000 001,
|
|
|
|
|
// 2, 1.9999999971613847, 6.9560343643347764};
|
2020-12-09 16:58:48 +01:00
|
|
|
|
// {1, 5.9277};
|
2020-12-16 20:31:58 +01:00
|
|
|
|
// {0.0001, 2, 2.000000005047502, 1.3055270196964464};
|
2020-12-09 16:58:48 +01:00
|
|
|
|
// {initialGuess(0), initialGuess(1), initialGuess(2),
|
|
|
|
|
// initialGuess(3)};
|
|
|
|
|
// assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) {
|
|
|
|
|
// return d >= xMax || d <= xMin;
|
|
|
|
|
// }));
|
2021-01-12 13:41:40 +01:00
|
|
|
|
// std::vector<double> xLow(x.size(), xMin);
|
|
|
|
|
// std::vector<double> xUpper(x.size(), xMax);
|
|
|
|
|
// 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); });
|
2020-12-22 17:47:29 +01:00
|
|
|
|
// const double rhobeg = std::min(0.95, 0.2 * maxX);
|
2021-01-12 13:41:40 +01:00
|
|
|
|
// double rhobeg = 1;
|
|
|
|
|
// double rhoend = rhobeg * 1e-8;
|
|
|
|
|
// const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2;
|
|
|
|
|
// std::vector<double> w(wSize);
|
|
|
|
|
// const size_t maxFun = std::min(100.0 * (x.size() + 1), 1000.0);
|
2021-01-22 15:39:36 +01:00
|
|
|
|
dlib::matrix<double, 0, 1> xMin(settings.xRanges.size());
|
|
|
|
|
dlib::matrix<double, 0, 1> xMax(settings.xRanges.size());
|
|
|
|
|
for (int i = 0; i < settings.xRanges.size(); i++) {
|
|
|
|
|
xMin(i) = settings.xRanges[i].min;
|
|
|
|
|
xMax(i) = settings.xRanges[i].max;
|
|
|
|
|
}
|
2021-01-12 13:41:40 +01:00
|
|
|
|
|
2021-01-31 22:36:46 +01:00
|
|
|
|
numberOfFunctionCalls = 0;
|
2021-01-22 15:39:36 +01:00
|
|
|
|
double (*objF)(double, double, double) = &objective;
|
2021-01-12 13:41:40 +01:00
|
|
|
|
auto start = std::chrono::system_clock::now();
|
|
|
|
|
dlib::function_evaluation result = dlib::find_min_global(
|
2021-01-29 18:07:13 +01:00
|
|
|
|
objF, xMin, xMax, dlib::max_function_calls(settings.maxSimulations),
|
|
|
|
|
std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy);
|
2021-01-12 13:41:40 +01:00
|
|
|
|
auto end = std::chrono::system_clock::now();
|
|
|
|
|
auto elapsed = std::chrono::duration_cast<std::chrono::seconds>(end - start);
|
2021-01-22 15:39:36 +01:00
|
|
|
|
Results results{numOfSimulationCrashes, minX, minY};
|
|
|
|
|
std::cout << "Finished optimizing." << endl;
|
|
|
|
|
// std::cout << "Solution x:" << endl;
|
|
|
|
|
// std::cout << result.x << endl;
|
|
|
|
|
std::cout << "Objective value:" << minY << endl;
|
|
|
|
|
// std::cout << result.y << endl;
|
|
|
|
|
// std::cout << minY << endl;
|
|
|
|
|
// std::cout << "Time(sec):" << elapsed.count() << std::endl;
|
|
|
|
|
// std::cout << "Max function evaluations:" << maxFun << std::endl;
|
|
|
|
|
// std::cout << "Initial guess:" << initialGuess << std::endl;
|
2021-01-12 13:41:40 +01:00
|
|
|
|
// 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();
|
|
|
|
|
// rhobeg *= 1e1;
|
|
|
|
|
// // rhoend *= 1e2;
|
|
|
|
|
// bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(),
|
|
|
|
|
// rhobeg, rhoend, maxFun, w.data());
|
|
|
|
|
// std::cout << "Finished second optimization round" << std::endl;
|
2020-12-09 16:58:48 +01:00
|
|
|
|
|
2021-01-22 15:39:36 +01:00
|
|
|
|
return results;
|
2020-11-23 10:06:45 +01:00
|
|
|
|
}
|
|
|
|
|
|
2020-12-16 20:31:58 +01:00
|
|
|
|
void ReducedModelOptimizer::setInitialGuess(std::vector<double> v) {
|
|
|
|
|
initialGuess = v;
|
|
|
|
|
}
|
|
|
|
|
|
2021-01-04 13:12:25 +01:00
|
|
|
|
std::vector<std::shared_ptr<SimulationJob>>
|
|
|
|
|
ReducedModelOptimizer::createScenarios(
|
2020-11-27 11:45:20 +01:00
|
|
|
|
const std::shared_ptr<SimulationMesh> &pMesh) {
|
2021-01-04 13:12:25 +01:00
|
|
|
|
std::vector<std::shared_ptr<SimulationJob>> scenarios;
|
|
|
|
|
scenarios.resize(SimulationScenario::NumberOfSimulationScenarios);
|
2020-11-27 11:45:20 +01:00
|
|
|
|
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
|
|
|
|
|
std::unordered_map<VertexIndex, Vector6d> nodalForces;
|
2020-12-09 16:58:48 +01:00
|
|
|
|
const double forceMagnitude = 1;
|
|
|
|
|
// Assuming the patterns lays on the x-y plane
|
|
|
|
|
const CoordType patternPlaneNormal(0, 0, 1);
|
2020-12-14 10:07:43 +01:00
|
|
|
|
// Make the first interface node lay on the x axis
|
2020-12-16 20:31:58 +01:00
|
|
|
|
// const size_t fullPatternFirstInterfaceNodeIndex =
|
|
|
|
|
// m_fullPatternOppositeInterfaceViMap.begin()->second;
|
|
|
|
|
// CoordType fullPatternFirstInterfaceNodePosition =
|
|
|
|
|
// m_pFullModelSimulationMesh->vert[fullPatternFirstInterfaceNodeIndex].cP();
|
|
|
|
|
// CoordType centerOfMass(0, 0, 0);
|
|
|
|
|
// for (size_t vi = 0; vi < pMesh->VN(); vi++) {
|
|
|
|
|
// centerOfMass = centerOfMass + pMesh->vert[vi].P();
|
2020-11-27 11:45:20 +01:00
|
|
|
|
// }
|
2020-12-16 20:31:58 +01:00
|
|
|
|
// centerOfMass /= pMesh->VN();
|
|
|
|
|
// vcg::tri::UpdatePosition<SimulationMesh>::Translate(
|
|
|
|
|
// *m_pFullModelSimulationMesh, -centerOfMass);
|
|
|
|
|
// vcg::tri::UpdatePosition<SimulationMesh>::Translate(
|
|
|
|
|
// *m_pReducedPatternSimulationMesh, centerOfMass);
|
|
|
|
|
|
|
|
|
|
// const vcg::Matrix33d R = vcg::RotationMatrix(
|
|
|
|
|
// fullPatternFirstInterfaceNodePosition,
|
|
|
|
|
// CoordType(fullPatternFirstInterfaceNodePosition.Norm(), 0, 0), false);
|
|
|
|
|
// std::for_each(m_pFullModelSimulationMesh->vert.begin(),
|
|
|
|
|
// m_pFullModelSimulationMesh->vert.end(), [&](auto &v) {
|
|
|
|
|
// v.P() = R * v.P();
|
|
|
|
|
// v.N() = R * v.N();
|
|
|
|
|
// });
|
|
|
|
|
// std::for_each(m_pReducedPatternSimulationMesh->vert.begin(),
|
|
|
|
|
// m_pReducedPatternSimulationMesh->vert.end(), [&](auto &v) {
|
|
|
|
|
// v.P() = R * v.P();
|
|
|
|
|
// v.N() = R * v.N();
|
|
|
|
|
// });
|
|
|
|
|
|
|
|
|
|
// m_pFullModelSimulationMesh->updateEigenEdgeAndVertices();
|
|
|
|
|
// m_pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
|
|
|
|
|
|
|
|
|
|
//// Axial
|
2021-01-04 13:12:25 +01:00
|
|
|
|
SimulationScenario scenarioName = SimulationScenario::Axial;
|
2020-12-16 20:31:58 +01:00
|
|
|
|
for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
|
|
|
|
|
CoordType forceDirection =
|
|
|
|
|
(pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
|
|
|
|
|
.Normalize();
|
|
|
|
|
nodalForces[viPair.first] = Vector6d({forceDirection[0], forceDirection[1],
|
|
|
|
|
forceDirection[2], 0, 0, 0}) *
|
|
|
|
|
forceMagnitude * 10;
|
|
|
|
|
fixedVertices[viPair.second] =
|
|
|
|
|
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
|
|
|
|
|
}
|
2021-01-04 13:12:25 +01:00
|
|
|
|
scenarios[scenarioName] = std::make_shared<SimulationJob>(
|
|
|
|
|
SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
|
|
|
|
|
fixedVertices, nodalForces, {}));
|
2020-11-27 11:45:20 +01:00
|
|
|
|
|
2021-01-04 13:12:25 +01:00
|
|
|
|
//// Shear
|
|
|
|
|
scenarioName = SimulationScenario::Shear;
|
2020-12-16 20:31:58 +01:00
|
|
|
|
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 = (v ^ patternPlaneNormal).Normalize();
|
|
|
|
|
nodalForces[viPair.first] = Vector6d({forceDirection[0], forceDirection[1],
|
|
|
|
|
forceDirection[2], 0, 0, 0}) *
|
|
|
|
|
0.40 * forceMagnitude;
|
|
|
|
|
fixedVertices[viPair.second] =
|
|
|
|
|
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
|
|
|
|
|
}
|
2021-01-04 13:12:25 +01:00
|
|
|
|
scenarios[scenarioName] = std::make_shared<SimulationJob>(
|
|
|
|
|
SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
|
|
|
|
|
fixedVertices, nodalForces, {}));
|
2020-11-27 11:45:20 +01:00
|
|
|
|
|
2020-12-14 10:07:43 +01:00
|
|
|
|
// //// Torsion
|
2020-12-09 16:58:48 +01:00
|
|
|
|
// 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();
|
2020-12-14 10:07:43 +01:00
|
|
|
|
// nodalForces[viPair.first] = Vector6d{
|
|
|
|
|
// 0, 0, 0, normalVDerivativeDir[0], normalVDerivativeDir[1], 0};
|
|
|
|
|
// fixedVertices[viPair.second] =
|
|
|
|
|
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
|
2020-12-09 16:58:48 +01:00
|
|
|
|
// fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
|
|
|
|
|
// } else {
|
2020-12-14 10:07:43 +01:00
|
|
|
|
// fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
|
|
|
|
|
// fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
|
2020-11-27 11:45:20 +01:00
|
|
|
|
// }
|
|
|
|
|
// }
|
2020-12-14 10:07:43 +01:00
|
|
|
|
// scenarios.push_back({pMesh, fixedVertices, nodalForces});
|
2020-12-09 16:58:48 +01:00
|
|
|
|
|
2021-01-04 13:12:25 +01:00
|
|
|
|
//// Bending
|
|
|
|
|
scenarioName = SimulationScenario::Bending;
|
2020-12-16 20:31:58 +01:00
|
|
|
|
fixedVertices.clear();
|
|
|
|
|
nodalForces.clear();
|
|
|
|
|
for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
|
|
|
|
|
nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) * 1;
|
|
|
|
|
fixedVertices[viPair.second] =
|
|
|
|
|
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
|
|
|
|
|
}
|
2021-01-04 13:12:25 +01:00
|
|
|
|
scenarios[scenarioName] = std::make_shared<SimulationJob>(
|
|
|
|
|
SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
|
|
|
|
|
fixedVertices, nodalForces, {}));
|
2020-11-27 11:45:20 +01:00
|
|
|
|
|
2020-12-16 20:31:58 +01:00
|
|
|
|
//// Double using moments
|
2021-01-04 13:12:25 +01:00
|
|
|
|
scenarioName = SimulationScenario::Dome;
|
2020-12-16 20:31:58 +01:00
|
|
|
|
fixedVertices.clear();
|
|
|
|
|
nodalForces.clear();
|
|
|
|
|
for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
|
|
|
|
|
viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) {
|
|
|
|
|
const auto viPair = *viPairIt;
|
|
|
|
|
if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
|
|
|
|
|
fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
|
|
|
|
|
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 2};
|
|
|
|
|
} else {
|
|
|
|
|
fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
|
|
|
|
|
fixedVertices[viPair.second] = std::unordered_set<DoFType>{2};
|
|
|
|
|
}
|
|
|
|
|
CoordType v =
|
|
|
|
|
(pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
|
|
|
|
|
.Normalize();
|
|
|
|
|
nodalForces[viPair.first] =
|
|
|
|
|
Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude * 0.1;
|
|
|
|
|
nodalForces[viPair.second] =
|
|
|
|
|
Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude * 0.1;
|
|
|
|
|
}
|
2021-01-04 13:12:25 +01:00
|
|
|
|
scenarios[scenarioName] = std::make_shared<SimulationJob>(
|
|
|
|
|
SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
|
|
|
|
|
fixedVertices, nodalForces, {}));
|
2020-11-23 10:06:45 +01:00
|
|
|
|
|
2020-12-14 10:07:43 +01:00
|
|
|
|
//// Saddle
|
2021-01-04 13:12:25 +01:00
|
|
|
|
scenarioName = SimulationScenario::Saddle;
|
2020-12-09 16:58:48 +01:00
|
|
|
|
fixedVertices.clear();
|
|
|
|
|
nodalForces.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();
|
|
|
|
|
if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
|
2020-12-14 10:07:43 +01:00
|
|
|
|
nodalForces[viPair.first] =
|
|
|
|
|
Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.02 * forceMagnitude;
|
|
|
|
|
nodalForces[viPair.second] =
|
|
|
|
|
Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.02 * forceMagnitude;
|
2020-12-09 16:58:48 +01:00
|
|
|
|
} else {
|
2020-12-14 10:07:43 +01:00
|
|
|
|
fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
|
|
|
|
|
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
|
2020-12-09 16:58:48 +01:00
|
|
|
|
|
2020-12-14 10:07:43 +01:00
|
|
|
|
nodalForces[viPair.first] =
|
|
|
|
|
Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.01 * forceMagnitude;
|
|
|
|
|
nodalForces[viPair.second] =
|
|
|
|
|
Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.01 * forceMagnitude;
|
2020-12-09 16:58:48 +01:00
|
|
|
|
}
|
|
|
|
|
}
|
2021-01-04 13:12:25 +01:00
|
|
|
|
scenarios[scenarioName] = std::make_shared<SimulationJob>(
|
|
|
|
|
SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
|
|
|
|
|
fixedVertices, nodalForces, {}));
|
2020-11-23 10:06:45 +01:00
|
|
|
|
|
|
|
|
|
return scenarios;
|
|
|
|
|
}
|
|
|
|
|
|
2021-01-04 13:12:25 +01:00
|
|
|
|
void ReducedModelOptimizer::runBeamOptimization() {
|
|
|
|
|
// load beams
|
|
|
|
|
VCGEdgeMesh fullBeam;
|
2021-01-17 12:46:33 +01:00
|
|
|
|
fullBeam.loadPly("/home/iason/Models/simple_beam_model_10elem_1m.ply");
|
2021-01-04 13:12:25 +01:00
|
|
|
|
VCGEdgeMesh reducedBeam;
|
2021-01-17 12:46:33 +01:00
|
|
|
|
reducedBeam.loadPly("/home/iason/Models/simple_beam_model_4elem_1m.ply");
|
2021-01-04 13:12:25 +01:00
|
|
|
|
fullBeam.registerForDrawing();
|
|
|
|
|
reducedBeam.registerForDrawing();
|
|
|
|
|
// polyscope::show();
|
|
|
|
|
// maps
|
|
|
|
|
std::unordered_map<size_t, size_t> displacementReducedToFullMap;
|
|
|
|
|
displacementReducedToFullMap[reducedBeam.VN() / 2] = fullBeam.VN() / 2;
|
|
|
|
|
g_reducedToFullViMap = displacementReducedToFullMap;
|
|
|
|
|
std::unordered_map<size_t, size_t> jobFullToReducedMap;
|
|
|
|
|
jobFullToReducedMap[0] = 0;
|
|
|
|
|
jobFullToReducedMap[fullBeam.VN() - 1] = reducedBeam.VN() - 1;
|
|
|
|
|
|
|
|
|
|
// full model simuilation job
|
|
|
|
|
auto pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullBeam);
|
|
|
|
|
pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02});
|
2021-01-17 12:46:33 +01:00
|
|
|
|
pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
std::unordered_map<VertexIndex, std::unordered_set<int>> fixedVertices;
|
|
|
|
|
fixedVertices[0] = ::unordered_set<int>({0, 1, 2, 3, 4, 5});
|
|
|
|
|
std::unordered_map<VertexIndex, Vector6d> forces;
|
|
|
|
|
forces[fullBeam.VN() - 1] = Vector6d({0, 0, 10, 0, 0, 0});
|
|
|
|
|
const std::string fullBeamSimulationJobLabel = "Pull_Z";
|
|
|
|
|
std::shared_ptr<SimulationJob> pFullModelSimulationJob =
|
|
|
|
|
make_shared<SimulationJob>(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<SimulationMesh>(reducedBeam);
|
|
|
|
|
reducedSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02});
|
2021-01-17 12:46:33 +01:00
|
|
|
|
reducedSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
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<SimulationJob>(
|
|
|
|
|
reducedSimulationMesh, fullBeamSimulationJobLabel,
|
|
|
|
|
reducedSimJob.constrainedVertices, reducedSimJob.nodalExternalForces);
|
2021-01-17 12:46:33 +01:00
|
|
|
|
initializeOptimizationParameters(reducedSimulationMesh);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
|
|
|
|
|
// 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"));
|
|
|
|
|
|
2021-01-22 15:39:36 +01:00
|
|
|
|
runOptimization({}, &objective);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
|
|
|
|
|
fullModelResults.registerForDrawing();
|
|
|
|
|
SimulationResults reducedModelResults = simulator.executeSimulation(
|
2021-01-12 13:41:40 +01:00
|
|
|
|
g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1]);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
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<std::shared_ptr<SimulationJob>>
|
|
|
|
|
&fullPatternSimulationJobs,
|
|
|
|
|
const std::vector<SimulationScenario> &simulationScenarios) {
|
|
|
|
|
m_pFullPatternSimulationMesh->registerForDrawing();
|
|
|
|
|
for (const int simulationScenarioIndex : simulationScenarios) {
|
|
|
|
|
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
|
|
|
|
|
fullPatternSimulationJobs[simulationScenarioIndex];
|
|
|
|
|
pFullPatternSimulationJob->registerForDrawing(
|
|
|
|
|
m_pFullPatternSimulationMesh->getLabel());
|
|
|
|
|
SimulationResults fullModelResults =
|
2021-01-12 13:41:40 +01:00
|
|
|
|
simulator.executeSimulation(pFullPatternSimulationJob);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
fullModelResults.registerForDrawing();
|
|
|
|
|
fullModelResults.saveDeformedModel();
|
|
|
|
|
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
|
|
|
|
|
g_reducedPatternSimulationJob[simulationScenarioIndex];
|
|
|
|
|
SimulationResults reducedModelResults =
|
2021-01-12 13:41:40 +01:00
|
|
|
|
simulator.executeSimulation(pReducedPatternSimulationJob);
|
2021-01-22 15:39:36 +01:00
|
|
|
|
double error = computeError(
|
2021-01-04 13:12:25 +01:00
|
|
|
|
reducedModelResults,
|
|
|
|
|
g_optimalReducedModelDisplacements[simulationScenarioIndex]);
|
|
|
|
|
std::cout << "Error of simulation scenario "
|
|
|
|
|
<< simulationScenarioStrings[simulationScenarioIndex] << " is "
|
|
|
|
|
<< error << std::endl;
|
|
|
|
|
reducedModelResults.registerForDrawing();
|
2021-01-12 13:41:40 +01:00
|
|
|
|
// firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing();
|
|
|
|
|
// reducedModelResults.saveDeformedModel();
|
2021-01-04 13:12:25 +01:00
|
|
|
|
// registerWorldAxes();
|
2021-01-12 13:41:40 +01:00
|
|
|
|
const std::string screenshotFilename =
|
|
|
|
|
"/home/iason/Coding/Projects/Approximating shapes with flat "
|
|
|
|
|
"patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" +
|
|
|
|
|
m_pFullPatternSimulationMesh->getLabel() + "_" +
|
|
|
|
|
simulationScenarioStrings[simulationScenarioIndex];
|
2021-01-04 13:12:25 +01:00
|
|
|
|
polyscope::show();
|
2021-01-12 13:41:40 +01:00
|
|
|
|
polyscope::screenshot(screenshotFilename, false);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
fullModelResults.unregister();
|
|
|
|
|
reducedModelResults.unregister();
|
2021-01-12 13:41:40 +01:00
|
|
|
|
// firstOptimizationRoundResults[simulationScenarioIndex].unregister();
|
2021-01-04 13:12:25 +01:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2021-01-22 15:39:36 +01:00
|
|
|
|
ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
|
|
|
|
|
const Settings &xRanges,
|
2021-01-12 13:41:40 +01:00
|
|
|
|
const std::vector<SimulationScenario> &simulationScenarios) {
|
|
|
|
|
|
|
|
|
|
g_simulationScenarioIndices = simulationScenarios;
|
|
|
|
|
if (g_simulationScenarioIndices.empty()) {
|
|
|
|
|
g_simulationScenarioIndices = {
|
|
|
|
|
SimulationScenario::Axial, SimulationScenario::Shear,
|
|
|
|
|
SimulationScenario::Bending, SimulationScenario::Dome,
|
|
|
|
|
SimulationScenario::Saddle};
|
|
|
|
|
}
|
2020-12-21 16:56:21 +01:00
|
|
|
|
|
2021-01-04 13:12:25 +01:00
|
|
|
|
std::vector<std::shared_ptr<SimulationJob>> simulationJobs =
|
|
|
|
|
createScenarios(m_pFullPatternSimulationMesh);
|
2020-12-21 16:56:21 +01:00
|
|
|
|
g_optimalReducedModelDisplacements.resize(6);
|
|
|
|
|
g_reducedPatternSimulationJob.resize(6);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
g_firstRoundIterationIndex = 0;
|
2021-01-12 13:41:40 +01:00
|
|
|
|
minY = std::numeric_limits<double>::max();
|
2021-01-22 15:39:36 +01:00
|
|
|
|
numOfSimulationCrashes = 0;
|
2020-12-22 17:47:29 +01:00
|
|
|
|
// polyscope::removeAllStructures();
|
2020-12-21 16:56:21 +01:00
|
|
|
|
|
2021-01-17 12:46:33 +01:00
|
|
|
|
FormFinder::Settings settings;
|
|
|
|
|
// settings.shouldDraw = true;
|
2021-01-12 13:41:40 +01:00
|
|
|
|
for (int simulationScenarioIndex : g_simulationScenarioIndices) {
|
2021-01-04 13:12:25 +01:00
|
|
|
|
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
|
2020-12-21 16:56:21 +01:00
|
|
|
|
simulationJobs[simulationScenarioIndex];
|
|
|
|
|
SimulationResults fullModelResults =
|
2021-01-17 12:46:33 +01:00
|
|
|
|
simulator.executeSimulation(pFullPatternSimulationJob, settings);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
g_optimalReducedModelDisplacements[simulationScenarioIndex].resize(
|
|
|
|
|
m_pReducedPatternSimulationMesh->VN(), 3);
|
2020-12-21 16:56:21 +01:00
|
|
|
|
computeDesiredReducedModelDisplacements(
|
2021-01-04 13:12:25 +01:00
|
|
|
|
fullModelResults, g_reducedToFullViMap,
|
2020-12-21 16:56:21 +01:00
|
|
|
|
g_optimalReducedModelDisplacements[simulationScenarioIndex]);
|
2021-01-04 13:12:25 +01:00
|
|
|
|
SimulationJob reducedPatternSimulationJob;
|
|
|
|
|
reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
|
|
|
|
|
computeReducedModelSimulationJob(*pFullPatternSimulationJob,
|
|
|
|
|
m_fullToReducedInterfaceViMap,
|
|
|
|
|
reducedPatternSimulationJob);
|
|
|
|
|
g_reducedPatternSimulationJob[simulationScenarioIndex] =
|
|
|
|
|
std::make_shared<SimulationJob>(reducedPatternSimulationJob);
|
2020-11-23 10:06:45 +01:00
|
|
|
|
}
|
2021-01-22 15:39:36 +01:00
|
|
|
|
Results optResults = runOptimization(xRanges, &objective);
|
|
|
|
|
updateMesh(optResults.x.size(), optResults.x.data());
|
2021-01-17 12:46:33 +01:00
|
|
|
|
// visualizeResults(simulationJobs, g_simulationScenarioIndices);
|
2021-01-22 15:39:36 +01:00
|
|
|
|
return optResults;
|
2020-11-23 10:06:45 +01:00
|
|
|
|
}
|