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

Exporting of result to excel file

This commit is contained in:
Iason 2021-01-22 16:39:36 +02:00
parent 5eab2511e9
commit a61856211b
3 changed files with 208 additions and 274 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

303
src/main.cpp
View File

@ -1,4 +1,5 @@
#include "beamformfinder.hpp"
#include "csvfile.hpp"
#include "edgemesh.hpp"
#include "flatpattern.hpp"
#include "polyscope/curve_network.h"
@ -14,42 +15,7 @@
#include <string>
#include <vcg/complex/algorithms/update/position.h>
// void scale(const std::vector<size_t>& numberOfNodesPerSlot,
// FlatPattern &pattern) {
// const double desiredSize = 0.0025; // center to boundary
// const size_t interfaceSlotIndex = 4; // bottom edge
// std::unordered_map<size_t, size_t> nodeToSlot;
// std::unordered_map<size_t, std::unordered_set<size_t>> 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<FlatPattern>::Scale(full, scaleFactor);
//}
int main(int argc, char *argv[]) {
// ReducedModelOptimizer::runBeamOptimization();
// FlatPattern pattern("/home/iason/Models/valid_6777.ply");
// FlatPattern pattern("/home/iason/Models/simple_beam_paper_example.ply");
// pattern.savePly("fannedValid.ply");
// SimulationJob job;
// job.load("../"
// "ProblematicSimulationJobs/17:16_4_1_2021/"
// "reduced_pattern_Single bar reduced model_simScenario.json");
// FormFinder ff;
// ff.executeSimulation(std::make_shared<SimulationJob>(job), false, false,
// true);
// Create reduced models
// FormFinder::runUnitTests();
@ -76,179 +42,106 @@ int main(int argc, char *argv[]) {
&CWReducedModel, &CCWReducedModel};
ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
for (double dimRationMax = 0.75; dimRationMax < 2; dimRationMax += 0.05) {
ReducedModelOptimizer::xRange beamWidth{"B", 0.5, 1.5};
ReducedModelOptimizer::xRange beamDimensionsRatio{"bOverh", 0.7,
dimRationMax};
ReducedModelOptimizer::xRange beamE{"E", 0.07, 1.5};
std::string xRangesString = beamWidth.toString() + " " +
beamDimensionsRatio.toString() + " " +
beamE.toString();
std::cout << xRangesString << std::endl;
ReducedModelOptimizer::Settings settings;
settings.xRanges = {beamWidth, beamDimensionsRatio, beamE};
// FlatPattern pattern("/home/iason/Models/TestSet_validPatterns/59.ply");
// optimizer.initialize(pattern, *reducedModels[1], {});
// optimizer.setInitialGuess
// ({0.050000000000000003, 4.0001308096448325, 4.2377893536538567,
// 5.6122373437520334});
// Eigen::VectorXd optimalParameters = optimizer.optimize(0);
std::filesystem::path thisOptimizationDirectory(
std::filesystem::path("../OptimizationResults").append(xRangesString));
std::filesystem::create_directory(thisOptimizationDirectory);
csvfile thisOptimizationStatistics(
std::filesystem::path(thisOptimizationDirectory)
.append("statistics.csv")
.string(),
true);
std::string fullPatternsTestSetDirectory =
// "/home/iason/Models/TestSet_validPatterns";
"/home/iason/Documents/PhD/Research/Approximating shapes with flat "
"patterns/Pattern_enumerator/Results/1v_0v_2e_1e_1c_6fan/3/Valid";
for (const auto &entry :
filesystem::directory_iterator(fullPatternsTestSetDirectory)) {
const auto filepath =
// std::filesystem::path(fullPatternsTestSetDirectory).append("305.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;
double totalError = 0;
int totalNumberOfSimulationCrashes = 0;
std::vector<double> errors;
std::string fullPatternsTestSetDirectory = "../TestSet";
// "/home/iason/Documents/PhD/Research/Approximating shapes with flat "
// "patterns/Pattern_enumerator/Results/1v_0v_2e_1e_1c_6fan/3/Valid";
for (const auto &entry :
filesystem::directory_iterator(fullPatternsTestSetDirectory)) {
const auto filepath =
// std::filesystem::path(fullPatternsTestSetDirectory).append("305.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;
}
// std::cout << "Full pattern:" << filepathString << std::endl;
FlatPattern pattern(filepathString);
pattern.setLabel(filepath.stem().string());
pattern.scale(0.03);
// for (int reducedPatternIndex = 0;
// reducedPatternIndex < reducedModels.size();
// reducedPatternIndex++) {
// FlatPattern *pReducedModel =
// reducedModels[reducedPatternIndex];
std::unordered_set<size_t> optimizationExcludedEi;
// if (pReducedModel !=
// reducedModels[0]) { // assumes that the singleBar reduced model
// // is
// // the
// // first in the reducedModels vector
// optimizationExcludedEi.insert(0);
// }
// FlatPattern cp;
// cp.copy(*reducedModels[0]);
optimizer.initializePatterns(pattern, *reducedModels[0],
optimizationExcludedEi);
// optimizer.optimize({ReducedModelOptimizer::Axial});
ReducedModelOptimizer::Results optimizationResults =
optimizer.optimize(settings);
errors.push_back(optimizationResults.objectiveValue);
SimulationResultsReporter::createPlot(
"", "Objective value", errors,
std::filesystem::path(thisOptimizationDirectory)
.append("ObjectiveValues.png")
.string());
thisOptimizationStatistics << filepath.stem().string()
<< optimizationResults.objectiveValue;
if (optimizationResults.numberOfSimulationCrashes == 0) {
thisOptimizationStatistics << "No crashes";
} else {
thisOptimizationStatistics
<< optimizationResults.numberOfSimulationCrashes;
}
thisOptimizationStatistics << endrow;
totalError += optimizationResults.objectiveValue;
totalNumberOfSimulationCrashes +=
optimizationResults.numberOfSimulationCrashes;
// }
}
std::cout << "Full pattern:" << filepathString << std::endl;
FlatPattern pattern(filepathString);
pattern.setLabel(filepath.stem().string());
pattern.scale(0.03);
// for (int reducedPatternIndex = 0;
// reducedPatternIndex < reducedModels.size();
// reducedPatternIndex++) {
// FlatPattern *pReducedModel =
// reducedModels[reducedPatternIndex];
std::unordered_set<size_t> optimizationExcludedEi;
// if (pReducedModel !=
// reducedModels[0]) { // assumes that the singleBar reduced model
// // is
// // the
// // first in the reducedModels vector
// optimizationExcludedEi.insert(0);
// }
// FlatPattern cp;
// cp.copy(*reducedModels[0]);
optimizer.initializePatterns(pattern, *reducedModels[0],
optimizationExcludedEi);
optimizer.optimize({ReducedModelOptimizer::SimulationScenario::Bending});
// }
csvfile statistics(std::filesystem::path("../OptimizationResults")
.append("statistics.csv")
.string(),
false);
for (const auto &range : settings.xRanges) {
statistics << range.min << range.max;
}
statistics << settings.maxSimulations;
if (totalNumberOfSimulationCrashes == 0) {
statistics << "No crashes";
} else {
statistics << totalNumberOfSimulationCrashes;
}
statistics << totalError << endrow;
}
// // Full model simulation
// std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
// fixedVertices;
// // fixedVertices[0] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// fixedVertices[3] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// // fixedVertices[7] = std::unordered_set<DoFType>{0, 1, 2};
// std::unordered_map<VertexIndex, Vector6d> nodalForces{
// {15, {0, 0, 2000, 0, 0, 0}}};
// SimulationJob fullModelSimulationJob{
// std::make_shared<ElementalMesh>(pattern), fixedVertices, nodalForces,
// {}};
// Simulator formFinder;
// SimulationResults fullModelResults =
// formFinder.executeSimulation(fullModelSimulationJob);
// fullModelResults.simulationLabel = "Full Model";
// fullModelResults.draw(fullModelSimulationJob);
// fullModelSimulationJob.draw();
// double stiffnessFactor = 1;
// while (true) {
// Reduced model simulation
// SimulationJob reducedModelSimulationJob =
// optimizer.getReducedSimulationJob(fullModelSimulationJob);
// const std::vector<double> stiffnessVector{
// fullModelSimulationJob.mesh->elements[0].properties.A,
// stiffnessFactor *
// fullModelSimulationJob.mesh->elements[0].properties.J, stiffnessFactor
// * fullModelSimulationJob.mesh->elements[0].properties.I2,
// stiffnessFactor *
// fullModelSimulationJob.mesh->elements[0].properties.I3};
// for (EdgeIndex ei = 0; ei < reducedModelSimulationJob.mesh->EN(); ei++) {
// BeamFormFinder::Element &e =
// reducedModelSimulationJob.mesh->elements[ei]; e.properties.A =
// 0.00035185827018667374;
// // stifnessVector[0];
// e.properties.J = // stiffnessVector[1];
// 7.709325104874406e-08;
// e.properties.I2 = -0.0000015661453308127776;
// // stiffnessVector[2];
// e.properties.I3 = 3.7099813776947167e-07; // stiffnessVector[3];
// e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength;
// e.torsionConstFactor = e.properties.G * e.properties.J /
// e.initialLength; e.firstBendingConstFactor =
// 2 * e.properties.E * e.properties.I2 / e.initialLength;
// e.secondBendingConstFactor =
// 2 * e.properties.E * e.properties.I3 / e.initialLength;
// }
// SimulationResults reducedModelResults =
// formFinder.executeSimulation(reducedModelSimulationJob, true);
// reducedModelResults.simulationLabel =
// "Reduced Model_" + std::to_string(stiffnessFactor);
// reducedModelResults.draw(reducedModelSimulationJob);
// SimulationResultsReporter resultsReporter;
// resultsReporter.reportResults(
// {reducedModelResults},
// std::filesystem::current_path().append("Results"),
// "_" + std::to_string(stiffnessFactor));
// if (reducedModelResults.history.numberOfSteps ==
// BeamFormFinder::Simulator::maxDRMIterations) {
// break;
// }
// stiffnessFactor *= 1.5;
// }
// // Beam
// // registerWorldAxes();
// VCGEdgeMesh mesh;
// // mesh.loadFromPly("/home/iason/Models/simple_beam_model_10elem_1m.ply");
// mesh.loadFromPly("/home/iason/Coding/Projects/Approximating shapes with
// flat "
// "patterns/RodModelOptimizationForPatterns/build/Debug/"
// "Fanned_Single bar reduced model.ply");
// // vcg::Matrix44d R;
// // R.SetRotateDeg(45, {0, 0, 1});
// // vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(mesh, R);
// // mesh.updateEigenEdgeAndVertices();
// FormFinder formFinder;
// formFinder.runUnitTests();
// std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
// fixedVertices; fixedVertices[1] = std::unordered_set<DoFType>{2};
// fixedVertices[2] = std::unordered_set<DoFType>{2};
// fixedVertices[3] = std::unordered_set<DoFType>{1, 2};
// fixedVertices[4] = std::unordered_set<DoFType>{2};
// fixedVertices[5] = std::unordered_set<DoFType>{2};
// fixedVertices[6] = std::unordered_set<DoFType>{0, 1, 2};
// // Forced displacements
// std::unordered_map<size_t, Eigen::Vector3d> nodalForcedDisplacements;
// // nodalForcedDisplacements[10] = Eigen::Vector3d(0, 0.1, 0.1);
// std::unordered_map<size_t, VectorType> nodalForcedNormal;
// CoordType v14 = (mesh.vert[4].cP() - mesh.vert[0].cP()).Normalize();
// CoordType v25 = (mesh.vert[5].cP() - mesh.vert[0].cP()).Normalize();
// CoordType v36 = (mesh.vert[6].cP() - mesh.vert[0].cP()).Normalize();
// const double forceMagnitude = 0.4;
// std::unordered_map<VertexIndex, Vector6d> nodalForces{
// // {4, Vector6d({0, 0, 0, v14[0], v14[1], 0}) * forceMagnitude},
// // {1, Vector6d({0, 0, 0, -v14[0], -v14[1], 0}) *
// forceMagnitude},
// // {5, Vector6d({0, 0, 0, v25[0], v25[1], 0}) * forceMagnitude},
// // {2, Vector6d({0, 0, 0, -v25[0], -v25[1], 0}) *
// forceMagnitude},
// // {6, Vector6d({0, 0, 0, v36[0], v36[1], 0}) * forceMagnitude},
// {3, Vector6d({0, 0, 0, -v36[0], -v36[1], 0}) * forceMagnitude}};
// // nodalForcedNormal[10] = v;
// // nodalForcedNormal[0] = -v;
// // fixedVertices[10] = {0, 1};
// SimulationJob beamSimulationJob{std::make_shared<SimulationMesh>(mesh),
// fixedVertices,
// nodalForces,
// {},
// {}};
// beamSimulationJob.mesh->setBeamMaterial(0.3, 1);
// beamSimulationJob.mesh->setBeamCrossSection(
// RectangularBeamDimensions{0.002, 0.002});
// beamSimulationJob.registerForDrawing();
// registerWorldAxes();
// polyscope::show();
// SimulationResults beamSimulationResults =
// formFinder.executeSimulation(beamSimulationJob, true, true);
// beamSimulationResults.label = "Beam";
// beamSimulationResults.registerForDrawing(beamSimulationJob);
// polyscope::show();
return 0;
}

147
src/reducedmodeloptimizer.cpp
View File

@ -33,7 +33,8 @@ std::vector<SimulationResults> firstOptimizationRoundResults;
int g_firstRoundIterationIndex{0};
double minY{std::numeric_limits<double>::max()};
std::vector<double> minX;
std::vector<double> failedSimulationsXRatio;
std::vector<std::vector<double>> failedSimulationsXRatio;
int numOfSimulationCrashes{false};
// struct OptimizationCallback {
// double operator()(const size_t &iterations, const Eigen::VectorXd &x,
@ -157,8 +158,11 @@ void updateMesh(long n, const double *x) {
// e.properties.E = g_initialParameters * x[ei];
// e.properties.E = g_initialParameters(0) * x[0];
// e.properties.G = g_initialParameters(1) * x[1];
e.setDimensions(RectangularBeamDimensions(g_initialParameters(0) * x[0],
g_initialParameters(1) * x[1]));
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]));
// e.properties.A = g_initialParameters(0) * x[0];
// e.properties.J = g_initialParameters(1) * x[1];
// e.properties.I2 = g_initialParameters(2) * x[2];
@ -197,8 +201,8 @@ void updateMesh(long n, const double *x) {
}
}
double ReducedModelOptimizer::objective(double b, double h) {
std::vector<double> x{b, h};
double ReducedModelOptimizer::objective(double b, double h, double E) {
std::vector<double> x{b, h, E};
return ReducedModelOptimizer::objective(x.size(), x.data());
}
@ -209,12 +213,12 @@ double ReducedModelOptimizer::objective(double x0, double x1, double x2,
}
double ReducedModelOptimizer::objective(long n, const double *x) {
std::cout.precision(17);
// std::cout.precision(17);
for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
std::cout << "x[" + std::to_string(parameterIndex) + "]="
<< x[parameterIndex] << std::endl;
}
// for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
// std::cout << "x[" + std::to_string(parameterIndex) + "]="
// << x[parameterIndex] << std::endl;
// }
const Element &e = g_reducedPatternSimulationJob[0]->pMesh->elements[0];
// std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
// << e.firstBendingConstFactor << " " <<
@ -235,21 +239,35 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
g_reducedPatternSimulationJob[simulationScenarioIndex],
simulationSettings);
std::string filename;
if (!reducedModelResults.converged) {
if (!reducedModelResults.converged /*&&
g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]
->pMesh->elements[0]
.A > 1e-8 &
x[0] / x[1] < 60*/) {
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]);
// 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("", "b/h", failedSimulationsXRatio);
// SimulationResultsReporter::createPlot(
// "log(b/h)", "log(A)", failedSimulationsXRatio[0],
// failedSimulationsXRatio[1], "ratioToAPlot.png");
// std::cout << "Failed simulation" << std::endl;
// simulationSettings.shouldDraw = true;
// simulationSettings.debugMessages = true;
// simulator.executeSimulation(
// g_reducedPatternSimulationJob[simulationScenarioIndex],
// simulationSettings);
numOfSimulationCrashes++;
} else {
error += computeError(
reducedModelResults,
@ -259,19 +277,23 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
"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";
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";
out.close();
}
if (error < minY) {
@ -289,7 +311,6 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
// [](double &c) { c = 0.7; });
// }
// gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory, 6, colors);
std::cout << std::endl;
// SimulationResultsReporter::createPlot("Number of Steps", "Objective
// value",
// gObjectiveValueHistory);
@ -488,7 +509,7 @@ void ReducedModelOptimizer::initializePatterns(
void ReducedModelOptimizer::initializeOptimizationParameters(
const std::shared_ptr<SimulationMesh> &mesh) {
const int numberOfOptimizationParameters = 2;
const int numberOfOptimizationParameters = 3;
g_initialParameters.resize(g_optimizeInnerHexagonSize
? numberOfOptimizationParameters + 1
: numberOfOptimizationParameters);
@ -503,9 +524,11 @@ void ReducedModelOptimizer::initializeOptimizationParameters(
// e.secondBendingConstFactor *= stiffnessFactor;
// }
const double initialB = std::sqrt(mesh->elements[0].A);
const double initialH = initialB;
const double initialRatio = 1;
;
g_initialParameters(0) = initialB;
g_initialParameters(1) = initialH;
g_initialParameters(1) = initialRatio;
g_initialParameters(2) = mesh->elements[0].material.youngsModulus;
// g_initialParameters =
// m_pReducedPatternSimulationMesh->elements[0].properties.E;
// for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) {
@ -575,13 +598,14 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
}
}
Eigen::VectorXd ReducedModelOptimizer::runOptimization(
ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(
const Settings &settings,
double (*pObjectiveFunction)(long, const double *)) {
gObjectiveValueHistory.clear();
const size_t n = g_initialParameters.rows();
assert(n == 2);
assert(n == 3);
// g_optimizeInnerHexagonSize ? 5: 4;
// const size_t npt = (n + 1) * (n + 2) / 2;
// // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
@ -603,8 +627,6 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization(
// {0.0001, 2, 2.000000005047502, 1.3055270196964464};
// {initialGuess(0), initialGuess(1), initialGuess(2),
// initialGuess(3)};
const double xMin = 1e-2;
const double xMax = 1e1;
// assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) {
// return d >= xMax || d <= xMin;
// }));
@ -623,25 +645,29 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization(
// 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);
const size_t maxFun = 150;
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;
}
double (*objF)(double, double) = &objective;
double (*objF)(double, double, double) = &objective;
auto start = std::chrono::system_clock::now();
dlib::function_evaluation result = dlib::find_min_global(
objF, {xMin, xMin}, {xMax, xMax}, dlib::max_function_calls(maxFun));
objF, xMin, xMax, dlib::max_function_calls(settings.maxSimulations));
auto end = std::chrono::system_clock::now();
auto elapsed = std::chrono::duration_cast<std::chrono::seconds>(end - start);
std::cout << "Finished optimization with dlib" << endl;
std::cout << "Solution x:" << endl;
std::cout << result.x << endl;
std::cout << "Solution y:" << 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 << "[" + std::to_string(xMin) + "," + std::to_string(xMax) + "]"
<< std::endl;
std::cout << "Initial guess:" << initialGuess << std::endl;
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;
// const size_t maxFun = 200;
// bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(),
// rhobeg, rhoend, maxFun, w.data());
@ -665,15 +691,7 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization(
// rhobeg, rhoend, maxFun, w.data());
// std::cout << "Finished second optimization round" << std::endl;
Eigen::VectorXd eigenX(x.size(), 1);
for (size_t xi = 0; xi < x.size(); xi++) {
eigenX(xi) = minX[xi];
// eigenX(xi) = x[xi];
// eigenX(xi) = result.x(xi);
}
// for (size_t xi = 0; xi < x.size(); xi++) {
// }
return eigenX;
return results;
}
void ReducedModelOptimizer::setInitialGuess(std::vector<double> v) {
@ -917,7 +935,7 @@ void ReducedModelOptimizer::runBeamOptimization() {
// .append(g_reducedPatternSimulationJob[0].mesh->getLabel() +
// "_simScenario.json"));
runOptimization(&objective);
runOptimization({}, &objective);
fullModelResults.registerForDrawing();
SimulationResults reducedModelResults = simulator.executeSimulation(
@ -938,7 +956,6 @@ void ReducedModelOptimizer::visualizeResults(
&fullPatternSimulationJobs,
const std::vector<SimulationScenario> &simulationScenarios) {
m_pFullPatternSimulationMesh->registerForDrawing();
double error = 0;
for (const int simulationScenarioIndex : simulationScenarios) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
fullPatternSimulationJobs[simulationScenarioIndex];
@ -952,7 +969,7 @@ void ReducedModelOptimizer::visualizeResults(
g_reducedPatternSimulationJob[simulationScenarioIndex];
SimulationResults reducedModelResults =
simulator.executeSimulation(pReducedPatternSimulationJob);
error += computeError(
double error = computeError(
reducedModelResults,
g_optimalReducedModelDisplacements[simulationScenarioIndex]);
std::cout << "Error of simulation scenario "
@ -975,7 +992,8 @@ void ReducedModelOptimizer::visualizeResults(
}
}
void ReducedModelOptimizer::optimize(
ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
const Settings &xRanges,
const std::vector<SimulationScenario> &simulationScenarios) {
g_simulationScenarioIndices = simulationScenarios;
@ -992,6 +1010,7 @@ void ReducedModelOptimizer::optimize(
g_reducedPatternSimulationJob.resize(6);
g_firstRoundIterationIndex = 0;
minY = std::numeric_limits<double>::max();
numOfSimulationCrashes = 0;
// polyscope::removeAllStructures();
FormFinder::Settings settings;
@ -1014,8 +1033,8 @@ void ReducedModelOptimizer::optimize(
g_reducedPatternSimulationJob[simulationScenarioIndex] =
std::make_shared<SimulationJob>(reducedPatternSimulationJob);
}
Eigen::VectorXd optimalParameters = runOptimization(&objective);
updateMesh(optimalParameters.rows(), optimalParameters.data());
Results optResults = runOptimization(xRanges, &objective);
updateMesh(optResults.x.size(), optResults.x.data());
// visualizeResults(simulationJobs, g_simulationScenarioIndices);
return optResults;
}

32
src/reducedmodeloptimizer.hpp
View File

@ -30,10 +30,31 @@ public:
Saddle,
NumberOfSimulationScenarios
};
struct Results {
int numberOfSimulationCrashes{0};
std::vector<double> x;
double objectiveValue;
};
struct xRange {
std::string label;
double min;
double max;
std::string toString() const {
return label + "=[" + std::to_string(min) + "," + std::to_string(max) +
"]";
}
};
struct Settings {
std::vector<xRange> xRanges;
int maxSimulations{300};
};
inline static const std::string simulationScenarioStrings[] = {
"Axial", "Shear", "Bending", "Double", "Saddle"};
void optimize(const std::vector<SimulationScenario> &simulationScenarios =
std::vector<SimulationScenario>());
Results optimize(const Settings &xRanges,
const std::vector<SimulationScenario> &simulationScenarios =
std::vector<SimulationScenario>());
double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;
ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot);
@ -57,7 +78,7 @@ public:
std::vector<double> &x);
static double objective(double x0, double x1, double x2, double x3);
static double objective(double b, double h);
static double objective(double b, double h, double E);
private:
void
@ -68,8 +89,9 @@ private:
const SimulationResults &fullModelResults,
const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
static Eigen::VectorXd
runOptimization(double (*pObjectiveFunction)(long, const double *));
static Results runOptimization(const Settings &settings,
double (*pObjectiveFunction)(long,
const double *));
std::vector<std::shared_ptr<SimulationJob>>
createScenarios(const std::shared_ptr<SimulationMesh> &pMesh);
void computeMaps(FlatPattern &fullModel, FlatPattern &reducedPattern,