Normalizing using an estimate of the maximum displacement. Removed the R and D command line argument.Refactoring

This commit is contained in:
Iason 2021-02-24 19:41:12 +02:00
parent e52b195f81
commit 94eeb0a9df
4 changed files with 423 additions and 339 deletions

View File

@ -37,7 +37,6 @@ download_project(PROJ MYSOURCES
endif()
#Polyscope
download_project(PROJ POLYSCOPE
GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git
@ -89,4 +88,5 @@ target_include_directories(${PROJECT_NAME}
)
link_directories(${MYSOURCES_SOURCE_DIR}/boost_graph/libs)
target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib)
target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen matplot dlib::dlib)
#target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib)

View File

@ -2,9 +2,6 @@
#include "csvfile.hpp"
#include "edgemesh.hpp"
#include "flatpattern.hpp"
#include "polyscope/curve_network.h"
#include "polyscope/point_cloud.h"
#include "polyscope/polyscope.h"
#include "reducedmodeloptimizer.hpp"
#include "simulationhistoryplotter.hpp"
#include "trianglepattterntopology.hpp"
@ -16,6 +13,11 @@
#include <string>
#include <vcg/complex/algorithms/update/position.h>
#ifdef POLYSCOPE_DEFINED
#include "polyscope/curve_network.h"
#include "polyscope/point_cloud.h"
#include "polyscope/polyscope.h"
#endif
int main(int argc, char *argv[]) {
if (argc < 3) {
std::cerr << "Specify at least the two pattern filepaths to be "
@ -46,10 +48,11 @@ int main(int argc, char *argv[]) {
ReducedModelOptimizer::Settings settings_optimization;
settings_optimization.xRanges = {beamWidth, beamDimensionsRatio, beamE,
innerHexagonSize};
const bool input_numberOfFunctionCallsDefined = argc >= 5;
const bool input_numberOfFunctionCallsDefined = argc >= 4;
settings_optimization.numberOfFunctionCalls =
input_numberOfFunctionCallsDefined ? std::atoi(argv[4]) : 100;
settings_optimization.normalizeObjectiveValue = true;
input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100;
settings_optimization.shouldNormalizeObjectiveValue = true;
settings_optimization.solutionAccuracy = 0.1;
// Optimize pair
const std::string pairName =
@ -62,9 +65,9 @@ int main(int argc, char *argv[]) {
optimizer.optimize(settings_optimization);
// Export results
const bool input_resultDirectoryDefined = argc >= 6;
const bool input_resultDirectoryDefined = argc >= 5;
std::string optimizationResultsDirectory =
input_resultDirectoryDefined ? argv[5] : "OptimizationResults";
input_resultDirectoryDefined ? argv[4] : "OptimizationResults";
if (optimizationResults.numberOfSimulationCrashes != 0) {
const auto crashedJobsDirPath =
std::filesystem::path(optimizationResultsDirectory)

View File

@ -9,7 +9,7 @@
struct GlobalOptimizationVariables {
std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
std::vector<std::vector<Vector6d>> fullPatternDisplacements;
std::vector<double> fullPatternDisplacementNormSum;
std::vector<double> reducedPatternMaximumDisplacementNormSum;
std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
@ -35,26 +35,49 @@ struct GlobalOptimizationVariables {
ReducedModelOptimizer::Settings optimizationSettings;
} global;
std::vector<SimulationJob> reducedPatternMaximumDisplacementSimulationJobs;
double ReducedModelOptimizer::computeError(
const std::vector<Vector6d> &reducedPatternDisplacements,
const std::vector<Vector6d> &fullPatternDisplacements,
const double &interfaceDisplacementsNormSum,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const double &normalizationFactor) {
const double rawError =
computeRawError(reducedPatternDisplacements, fullPatternDisplacements,
reducedToFullInterfaceViMap);
if (global.optimizationSettings.shouldNormalizeObjectiveValue) {
assert(rawError <= normalizationFactor);
if (rawError > normalizationFactor) {
std::cout << "Normalized error is bigger than one." << std::endl;
}
return rawError /
normalizationFactor; // std::max(interfaceDisplacementsNormSum,
// 0.00003);
}
return rawError;
}
double ReducedModelOptimizer::computeRawError(
const std::vector<Vector6d> &reducedPatternDisplacements,
const std::vector<Vector6d> &fullPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap) {
double error = 0;
for (const auto reducedFullViPair : reducedToFullInterfaceViMap) {
VertexIndex reducedModelVi = reducedFullViPair.first;
// const auto pos =
// g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP();
// std::cout << "Interface vi " << reducedModelVi << " is at position "
// << pos[0] << " " << pos[1] << " " << pos[2] << std::endl;
Eigen::Vector3d reducedVertexDisplacement(
reducedPatternDisplacements[reducedModelVi][0],
reducedPatternDisplacements[reducedModelVi][1],
reducedPatternDisplacements[reducedModelVi][2]);
reducedPatternDisplacements[reducedFullViPair.first][0],
reducedPatternDisplacements[reducedFullViPair.first][1],
reducedPatternDisplacements[reducedFullViPair.first][2]);
if (!std::isfinite(reducedVertexDisplacement[0]) ||
!std::isfinite(reducedVertexDisplacement[1]) ||
!std::isfinite(reducedVertexDisplacement[2])) {
std::cout << "Displacements are not finite" << std::endl;
std::terminate();
}
Eigen::Vector3d fullVertexDisplacement(
@ -67,9 +90,6 @@ double ReducedModelOptimizer::computeError(
error += errorVector.norm();
}
if (global.optimizationSettings.normalizeObjectiveValue) {
return error / std::max(interfaceDisplacementsNormSum, 0.00003);
}
return error;
}
@ -132,8 +152,8 @@ void updateMesh(long n, const double *x) {
pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
}
double ReducedModelOptimizer::objective(double b, double h, double E) {
std::vector<double> x{b, h, E};
double ReducedModelOptimizer::objective(double b, double r, double E) {
std::vector<double> x{b, r, E};
return ReducedModelOptimizer::objective(x.size(), x.data());
}
@ -146,12 +166,9 @@ double ReducedModelOptimizer::objective(double b, double h, double E,
double ReducedModelOptimizer::objective(long n, const double *x) {
// std::cout.precision(17);
// for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
// std::cout << "x[" + std::to_string(parameterIndex) + "]="
// << x[parameterIndex] << std::endl;
// }
const Element &e = global.reducedPatternSimulationJobs[0]->pMesh->elements[0];
// std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
// const Element &e =
// global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout <<
// e.axialConstFactor << " " << e.torsionConstFactor << " "
// << e.firstBendingConstFactor << " " <<
// e.secondBendingConstFactor
// << std::endl;
@ -167,51 +184,79 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
FormFinder::Settings simulationSettings;
// simulationSettings.shouldDraw = true;
for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
// std::cout << "Executing "
// << simulationScenarioStrings[simulationScenarioIndex]
// << std::endl;
SimulationResults reducedModelResults = simulator.executeSimulation(
global.reducedPatternSimulationJobs[simulationScenarioIndex],
simulationSettings);
std::string filename;
// std::string filename;
if (!reducedModelResults.converged) {
error += std::numeric_limits<double>::max();
global.numOfSimulationCrashes++;
if(beVerbose){
std::cout << "Failed simulation" << std::endl;
} }else {
error += computeError(
reducedModelResults.displacements,
global.fullPatternDisplacements[simulationScenarioIndex],
global.fullPatternDisplacementNormSum[simulationScenarioIndex],
global.reducedToFullInterfaceViMap);
}
std::ofstream out(filename, std::ios_base::app);
auto pMesh =
global
.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
->pMesh;
bool beVerbose = false; // TODO: use an extern or data member for that
if (beVerbose) {
std::cout << "Failed simulation" << std::endl;
}
} else {
const double thisSimulationScenarioError =
computeError(reducedModelResults.displacements,
global.fullPatternDisplacements[simulationScenarioIndex],
global.reducedToFullInterfaceViMap,
global.reducedPatternMaximumDisplacementNormSum
[simulationScenarioIndex]);
if (thisSimulationScenarioError > 1) {
std::cout << "Simulation scenario "
<< simulationScenarioStrings[simulationScenarioIndex]
<< " results in an error bigger than one." << std::endl;
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++) {
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();
#ifdef POLYSCOPE_DEFINED
ReducedModelOptimizer::visualizeResults(
global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs,
std::vector<SimulationScenario>{
static_cast<SimulationScenario>(simulationScenarioIndex)},
global.reducedToFullInterfaceViMap);
ReducedModelOptimizer::visualizeResults(
global.fullPatternSimulationJobs[simulationScenarioIndex],
global.reducedPatternSimulationJobs[simulationScenarioIndex],
global.reducedToFullInterfaceViMap, false);
ReducedModelOptimizer::visualizeResults(
global.fullPatternSimulationJobs[simulationScenarioIndex],
std::make_shared<SimulationJob>(
reducedPatternMaximumDisplacementSimulationJobs
[simulationScenarioIndex]),
global.reducedToFullInterfaceViMap, true);
polyscope::removeAllStructures();
#endif // POLYSCOPE_DEFINED
}
error += thisSimulationScenarioError;
}
// std::ofstream out(filename, std::ios_base::app);
// auto pMesh =
// global
// .reducedPatternSimulationJobs[global.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();
}
// std::cout << error << std::endl;
if (error < global.minY) {
@ -533,69 +578,50 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob(
#if POLYSCOPE_DEFINED
void ReducedModelOptimizer::visualizeResults(
const std::vector<std::shared_ptr<SimulationJob>>
&fullPatternSimulationJobs,
const std::vector<std::shared_ptr<SimulationJob>>
&reducedPatternSimulationJobs,
const std::vector<SimulationScenario> &simulationScenarios,
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap) {
&reducedToFullInterfaceViMap,
const bool draw = true) {
FormFinder simulator;
std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
fullPatternSimulationJobs[0]->pMesh;
pFullPatternSimulationJob->pMesh;
pFullPatternSimulationMesh->registerForDrawing();
pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() +
"_undeformed.ply");
reducedPatternSimulationJobs[0]->pMesh->savePly(
reducedPatternSimulationJobs[0]->pMesh->getLabel() + "_undeformed.ply");
double totalError = 0;
for (const int simulationScenarioIndex : simulationScenarios) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
fullPatternSimulationJobs[simulationScenarioIndex];
pFullPatternSimulationJob->registerForDrawing(
pFullPatternSimulationMesh->getLabel());
SimulationResults fullModelResults =
simulator.executeSimulation(pFullPatternSimulationJob);
fullModelResults.registerForDrawing();
// fullModelResults.saveDeformedModel();
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
reducedPatternSimulationJobs[simulationScenarioIndex];
SimulationResults reducedModelResults =
simulator.executeSimulation(pReducedPatternSimulationJob);
double interfaceDisplacementNormSum = 0;
for (const auto &interfaceViPair : reducedToFullInterfaceViMap) {
const int fullPatternInterfaceIndex = interfaceViPair.second;
Eigen::Vector3d fullPatternDisplacementVector(
fullModelResults.displacements[fullPatternInterfaceIndex][0],
fullModelResults.displacements[fullPatternInterfaceIndex][1],
fullModelResults.displacements[fullPatternInterfaceIndex][2]);
interfaceDisplacementNormSum += fullPatternDisplacementVector.norm();
}
reducedModelResults.saveDeformedModel();
fullModelResults.saveDeformedModel();
double error = computeError(
reducedModelResults.displacements, fullModelResults.displacements,
interfaceDisplacementNormSum, reducedToFullInterfaceViMap);
std::cout << "Error of simulation scenario "
<< simulationScenarioStrings[simulationScenarioIndex] << " is "
<< error << std::endl;
totalError += error;
reducedModelResults.registerForDrawing();
// firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing();
// registerWorldAxes();
// pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel()
// +
// "_undeformed.ply");
// reducedPatternSimulationJobs[0]->pMesh->savePly(
// reducedPatternSimulationJobs[0]->pMesh->getLabel() +
// "_undeformed.ply");
pFullPatternSimulationJob->registerForDrawing(
pFullPatternSimulationMesh->getLabel());
SimulationResults fullModelResults =
simulator.executeSimulation(pFullPatternSimulationJob);
fullModelResults.registerForDrawing();
// fullModelResults.saveDeformedModel();
SimulationResults reducedModelResults =
simulator.executeSimulation(pReducedPatternSimulationJob);
// reducedModelResults.saveDeformedModel();
// fullModelResults.saveDeformedModel();
double error = computeRawError(
reducedModelResults.displacements, fullModelResults.displacements,
reducedToFullInterfaceViMap/*,
global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/);
std::cout << "Error is " << error << std::endl;
reducedModelResults.registerForDrawing();
if (draw) {
polyscope::show();
const std::string screenshotFilename =
"/home/iason/Coding/Projects/Approximating shapes with flat "
"patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
"Images/" +
pFullPatternSimulationMesh->getLabel() + "_" +
simulationScenarioStrings[simulationScenarioIndex];
polyscope::show();
pFullPatternSimulationMesh->getLabel() + "_" + "noScenarioName";
// simulationScenarioStrings[simulationScenarioIndex];
polyscope::screenshot(screenshotFilename, false);
fullModelResults.unregister();
reducedModelResults.unregister();
// firstOptimizationRoundResults[simulationScenarioIndex].unregister();
pFullPatternSimulationMesh->unregister();
}
std::cout << "Total error:" << totalError << std::endl;
}
#endif // POLYSCOPE_DEFINED
@ -694,7 +720,19 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
results.numberOfSimulationCrashes = global.numOfSimulationCrashes;
results.x = global.minX;
results.objectiveValue = global.minY;
if (settings.shouldNormalizeObjectiveValue) {
global.optimizationSettings.shouldNormalizeObjectiveValue = false;
if (global.optimizeInnerHexagonSize) {
results.objectiveValue = objective(4, results.x.data());
} else {
results.objectiveValue =
objective(results.x[0], results.x[1], results.x[2]);
}
global.optimizationSettings.shouldNormalizeObjectiveValue = true;
}
results.time = elapsed.count() / 1000.0;
const bool printDebugInfo = false;
if (printDebugInfo) {
std::cout << "Finished optimizing." << endl;
// std::cout << "Solution x:" << endl;
@ -1006,7 +1044,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
global.g_optimalReducedModelDisplacements.resize(6);
global.reducedPatternSimulationJobs.resize(6);
global.fullPatternDisplacements.resize(6);
global.fullPatternDisplacementNormSum.resize(6);
global.reducedPatternMaximumDisplacementNormSum.resize(6);
global.g_firstRoundIterationIndex = 0;
global.minY = std::numeric_limits<double>::max();
global.numOfSimulationCrashes = 0;
@ -1014,34 +1052,18 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
global.optimizationSettings = optimizationSettings;
global.fullPatternSimulationJobs =
createScenarios(m_pFullPatternSimulationMesh);
reducedPatternMaximumDisplacementSimulationJobs.resize(6);
// polyscope::removeAllStructures();
FormFinder::Settings settings;
FormFinder::Settings simulationSettings;
// settings.shouldDraw = true;
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
global.fullPatternSimulationJobs[simulationScenarioIndex];
SimulationResults fullModelResults =
simulator.executeSimulation(pFullPatternSimulationJob, settings);
SimulationResults fullModelResults = simulator.executeSimulation(
pFullPatternSimulationJob, simulationSettings);
global.fullPatternDisplacements[simulationScenarioIndex] =
fullModelResults.displacements;
double interfaceDisplacementNormSum = 0;
for (const auto &interfaceViPair : global.reducedToFullInterfaceViMap) {
const int fullPatternInterfaceIndex = interfaceViPair.second;
Eigen::Vector3d fullPatternDisplacementVector(
fullModelResults.displacements[fullPatternInterfaceIndex][0],
fullModelResults.displacements[fullPatternInterfaceIndex][1],
fullModelResults.displacements[fullPatternInterfaceIndex][2]);
interfaceDisplacementNormSum += fullPatternDisplacementVector.norm();
}
global.fullPatternDisplacementNormSum[simulationScenarioIndex] =
interfaceDisplacementNormSum;
// global.g_optimalReducedModelDisplacements[simulationScenarioIndex].resize(
// m_pReducedPatternSimulationMesh->VN(), 3);
// computeDesiredReducedModelDisplacements(
// fullModelResults, global.reducedToFullInterfaceViMap,
// global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);
SimulationJob reducedPatternSimulationJob;
reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
computeReducedModelSimulationJob(*pFullPatternSimulationJob,
@ -1049,6 +1071,54 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
reducedPatternSimulationJob);
global.reducedPatternSimulationJobs[simulationScenarioIndex] =
std::make_shared<SimulationJob>(reducedPatternSimulationJob);
if (optimizationSettings.shouldNormalizeObjectiveValue) {
// Compute the quantities for normalizing the obj value
const double minB = optimizationSettings.xRanges[0].min;
const double maxRatio = optimizationSettings.xRanges[1].max;
const double minE = optimizationSettings.xRanges[2].min;
const double minHS = 0.3;
std::vector<double> mostFlexibleOptimizationParameters{minB, maxRatio,
minE, minHS};
if (global.optimizeInnerHexagonSize) {
updateMesh(4, mostFlexibleOptimizationParameters.data());
} else {
updateMesh(3, mostFlexibleOptimizationParameters.data());
}
reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex] =
global.reducedPatternSimulationJobs[simulationScenarioIndex]
->getCopy();
reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex]
.pMesh->setLabel("reduced_maxDisplacement");
SimulationResults reducedModelResults = simulator.executeSimulation(
global.reducedPatternSimulationJobs[simulationScenarioIndex],
simulationSettings);
const double errorOfMaxDisplacedReduced = computeRawError(
reducedModelResults.displacements, fullModelResults.displacements,
global.reducedToFullInterfaceViMap);
const double errorOfNonDisplacedReduced = computeRawError(
std::vector<Vector6d>(reducedModelResults.displacements.size(),
Vector6d(0)),
fullModelResults.displacements, global.reducedToFullInterfaceViMap);
double displacementMultiplier = 2;
// if (global.optimizeInnerHexagonSize) {
// displacementMultiplier = 2;
// } else {
// displacementMultiplier = 2;
// }
global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex] =
displacementMultiplier *
std::max(errorOfMaxDisplacedReduced, errorOfNonDisplacedReduced);
// if (errorOfMaxDisplacedReduced > errorOfNonDisplacedReduced) {
// std::cout << "Max disp results in a bigger error for scenario "
// << simulationScenarioStrings[simulationScenarioIndex]
// << std::endl;
// } else {
// std::cout << "Zero disp results in a bigger error for scenario "
// << simulationScenarioStrings[simulationScenarioIndex]
// << std::endl;
// }
}
}
Results optResults = runOptimization(optimizationSettings);
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
@ -1061,7 +1131,10 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
// optResults.draw();
// visualizeResults(simulationJobs, global.simulationScenarioIndices);
// visualizeResults(simulationJobs, global.reducedPatternSimulationJobs,
// global.simulationScenarioIndices,global.reducedToFullInterfaceViMap);
#ifdef POLYSCOPE_DEFINED
// visualizeResults(
// global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs,
// global.simulationScenarioIndices, global.reducedToFullInterfaceViMap);
#endif // POLYSCOPE_DEFINED
return optResults;
}

View File

@ -6,6 +6,7 @@
#include "edgemesh.hpp"
#include "elementalmesh.hpp"
#include "matplot/matplot.h"
#include "optimizationstructs.hpp"
#include <Eigen/Dense>
using FullPatternVertexIndex = VertexIndex;
@ -32,6 +33,7 @@ public:
NumberOfSimulationScenarios
};
struct Results;
struct xRange {
std::string label;
double min;
@ -41,13 +43,12 @@ public:
"]";
}
};
struct Results;
struct Settings {
std::vector<xRange> xRanges;
int numberOfFunctionCalls{100};
double solutionAccuracy{1e-2};
bool normalizeObjectiveValue{true};
bool shouldNormalizeObjectiveValue{true};
std::string toString() const {
std::string settingsString;
@ -58,9 +59,10 @@ public:
}
settingsString += xRangesString;
}
settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls) +
" Accuracy=" + std::to_string(solutionAccuracy) +
" Norm=" + (normalizeObjectiveValue ? "yes" : "no");
settingsString +=
"FuncCalls=" + std::to_string(numberOfFunctionCalls) +
" Accuracy=" + std::to_string(solutionAccuracy) +
" Norm=" + (shouldNormalizeObjectiveValue ? "yes" : "no");
return settingsString;
}
@ -88,9 +90,189 @@ public:
os << solutionAccuracy;
}
};
struct Results {
double time{-1};
int numberOfSimulationCrashes{0};
std::vector<double> x;
double objectiveValue;
std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
void save(const string &saveToPath) const {
assert(std::filesystem::is_directory(saveToPath));
const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
for (int simulationJobIndex = 0;
simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
fullPatternSimulationJobs[simulationJobIndex];
std::filesystem::path simulationJobFolderPath(
std::filesystem::path(saveToPath)
.append(pFullPatternSimulationJob->label));
std::filesystem::create_directory(simulationJobFolderPath);
const auto fullPatternDirectoryPath =
std::filesystem::path(simulationJobFolderPath).append("Full");
std::filesystem::create_directory(fullPatternDirectoryPath);
pFullPatternSimulationJob->save(fullPatternDirectoryPath.string());
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
reducedPatternSimulationJobs[simulationJobIndex];
const auto reducedPatternDirectoryPath =
std::filesystem::path(simulationJobFolderPath).append("Reduced");
if (!std::filesystem::exists(reducedPatternDirectoryPath)) {
std::filesystem::create_directory(reducedPatternDirectoryPath);
}
pReducedPatternSimulationJob->save(
reducedPatternDirectoryPath.string());
}
}
void load(const string &loadFromPath) {
assert(std::filesystem::is_directory(loadFromPath));
for (const auto &directoryEntry :
filesystem::directory_iterator(loadFromPath)) {
const auto simulationScenarioPath = directoryEntry.path();
if (!std::filesystem::is_directory(simulationScenarioPath)) {
continue;
}
// Load reduced pattern files
for (const auto &fileEntry : filesystem::directory_iterator(
std::filesystem::path(simulationScenarioPath)
.append("Full"))) {
const auto filepath = fileEntry.path();
if (filepath.extension() == ".json") {
SimulationJob job;
job.load(filepath.string());
fullPatternSimulationJobs.push_back(
std::make_shared<SimulationJob>(job));
}
}
// Load full pattern files
for (const auto &fileEntry : filesystem::directory_iterator(
std::filesystem::path(simulationScenarioPath)
.append("Reduced"))) {
const auto filepath = fileEntry.path();
if (filepath.extension() == ".json") {
SimulationJob job;
job.load(filepath.string());
reducedPatternSimulationJobs.push_back(
std::make_shared<SimulationJob>(job));
}
}
}
}
#if POLYSCOPE_DEFINED
void draw() const {
initPolyscope();
FormFinder simulator;
assert(fullPatternSimulationJobs.size() ==
reducedPatternSimulationJobs.size());
fullPatternSimulationJobs[0]->pMesh->registerForDrawing();
reducedPatternSimulationJobs[0]->pMesh->registerForDrawing();
const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
for (int simulationJobIndex = 0;
simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
// Drawing of full pattern results
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
fullPatternSimulationJobs[simulationJobIndex];
pFullPatternSimulationJob->registerForDrawing(
fullPatternSimulationJobs[0]->pMesh->getLabel());
SimulationResults fullModelResults =
simulator.executeSimulation(pFullPatternSimulationJob);
fullModelResults.registerForDrawing();
// Drawing of reduced pattern results
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
reducedPatternSimulationJobs[simulationJobIndex];
SimulationResults reducedModelResults =
simulator.executeSimulation(pReducedPatternSimulationJob);
reducedModelResults.registerForDrawing();
polyscope::show();
// Save a screensh
// const std::string screenshotFilename =
// "/home/iason/Coding/Projects/Approximating shapes with flat "
// "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
// + m_pFullPatternSimulationMesh->getLabel() + "_" +
// simulationScenarioStrings[simulationScenarioIndex];
// polyscope::screenshot(screenshotFilename, false);
fullModelResults.unregister();
reducedModelResults.unregister();
// double error = computeError(
// reducedModelResults,
// global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);
// std::cout << "Error of simulation scenario "
// << simulationScenarioStrings[simulationScenarioIndex] <<
// " is "
// << error << std::endl;
}
}
#endif // POLYSCOPE_DEFINED
void saveMeshFiles() const {
const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
assert(numberOfSimulationJobs != 0 &&
fullPatternSimulationJobs.size() ==
reducedPatternSimulationJobs.size());
fullPatternSimulationJobs[0]->pMesh->savePly(
"FullPattern_undeformed.ply");
reducedPatternSimulationJobs[0]->pMesh->savePly(
"ReducedPattern_undeformed.ply");
FormFinder simulator;
for (int simulationJobIndex = 0;
simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
// Drawing of full pattern results
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
fullPatternSimulationJobs[simulationJobIndex];
SimulationResults fullModelResults =
simulator.executeSimulation(pFullPatternSimulationJob);
fullModelResults.saveDeformedModel();
// Drawing of reduced pattern results
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
reducedPatternSimulationJobs[simulationJobIndex];
SimulationResults reducedModelResults =
simulator.executeSimulation(pReducedPatternSimulationJob);
reducedModelResults.saveDeformedModel();
}
}
void
writeHeaderTo(const ReducedModelOptimizer::Settings &settings_optimization,
csvFile &os) {
os << "Obj value";
for (const ReducedModelOptimizer::xRange &range :
settings_optimization.xRanges) {
os << range.label;
}
os << "Time(s)";
os << "#Crashes";
}
void
writeResultsTo(const ReducedModelOptimizer::Settings &settings_optimization,
csvFile &os) const {
os << objectiveValue;
for (const double &optimalX : x) {
os << optimalX;
}
for (int unusedXVarCounter = 0;
unusedXVarCounter < settings_optimization.xRanges.size() - x.size();
unusedXVarCounter++) {
os << "-";
}
os << time;
if (numberOfSimulationCrashes == 0) {
os << "-";
} else {
os << numberOfSimulationCrashes;
}
}
};
inline static const std::string simulationScenarioStrings[] = {
"Axial", "Shear", "Bending", "Double", "Saddle"};
"Axial", "Shear", "Bending", "Dome", "Saddle"};
Results optimize(const Settings &xRanges,
const std::vector<SimulationScenario> &simulationScenarios =
std::vector<SimulationScenario>());
@ -117,7 +299,7 @@ public:
std::vector<double> &x);
static double objective(double x0, double x1, double x2, double x3);
static double objective(double b, double h, double E);
static double objective(double b, double r, double E);
static std::vector<std::shared_ptr<SimulationJob>>
createScenarios(const std::shared_ptr<SimulationMesh> &pMesh,
@ -137,22 +319,28 @@ public:
&fullToReducedInterfaceViMap,
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
&fullPatternOppositeInterfaceViMap);
static void
visualizeResults(const std::vector<std::shared_ptr<SimulationJob>>
&fullPatternSimulationJobs,
const std::vector<std::shared_ptr<SimulationJob>>
&reducedPatternSimulationJobs,
const std::vector<SimulationScenario> &simulationScenarios,
const std::unordered_map<ReducedPatternVertexIndex,
FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
static void visualizeResults(
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
const std::unordered_map<ReducedPatternVertexIndex,
FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const bool dontShow);
static double computeError(const std::vector<Vector6d> &reducedPatternResults,
const std::vector<Vector6d> &fullPatternResults,
const double &interfaceDisplacementNormSum,
const std::unordered_map<ReducedPatternVertexIndex,
FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
static double
computeRawError(const std::vector<Vector6d> &reducedPatternDisplacements,
const std::vector<Vector6d> &fullPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex,
FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
static double
computeError(const std::vector<Vector6d> &reducedPatternDisplacements,
const std::vector<Vector6d> &fullPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex,
FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const double &normalizationFactor);
private:
void
@ -173,188 +361,8 @@ private:
static void
initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh);
static double
computeError(const SimulationResults &reducedPatternResults,
const Eigen::MatrixX3d &optimalReducedPatternDisplacements);
static double objective(long n, const double *x);
FormFinder simulator;
};
struct ReducedModelOptimizer::Results {
double time{-1};
int numberOfSimulationCrashes{0};
std::vector<double> x;
double objectiveValue;
std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
void save(const string &saveToPath) const {
assert(std::filesystem::is_directory(saveToPath));
const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
for (int simulationJobIndex = 0;
simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
fullPatternSimulationJobs[simulationJobIndex];
std::filesystem::path simulationJobFolderPath(
std::filesystem::path(saveToPath)
.append(pFullPatternSimulationJob->label));
std::filesystem::create_directory(simulationJobFolderPath);
const auto fullPatternDirectoryPath =
std::filesystem::path(simulationJobFolderPath).append("Full");
std::filesystem::create_directory(fullPatternDirectoryPath);
pFullPatternSimulationJob->save(fullPatternDirectoryPath.string());
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
reducedPatternSimulationJobs[simulationJobIndex];
const auto reducedPatternDirectoryPath =
std::filesystem::path(simulationJobFolderPath).append("Reduced");
if (!std::filesystem::exists(reducedPatternDirectoryPath)) {
std::filesystem::create_directory(reducedPatternDirectoryPath);
}
pReducedPatternSimulationJob->save(reducedPatternDirectoryPath.string());
}
}
void load(const string &loadFromPath) {
assert(std::filesystem::is_directory(loadFromPath));
for (const auto &directoryEntry :
filesystem::directory_iterator(loadFromPath)) {
const auto simulationScenarioPath = directoryEntry.path();
if (!std::filesystem::is_directory(simulationScenarioPath)) {
continue;
}
// Load reduced pattern files
for (const auto &fileEntry : filesystem::directory_iterator(
std::filesystem::path(simulationScenarioPath).append("Full"))) {
const auto filepath = fileEntry.path();
if (filepath.extension() == ".json") {
SimulationJob job;
job.load(filepath.string());
fullPatternSimulationJobs.push_back(
std::make_shared<SimulationJob>(job));
}
}
// Load full pattern files
for (const auto &fileEntry : filesystem::directory_iterator(
std::filesystem::path(simulationScenarioPath)
.append("Reduced"))) {
const auto filepath = fileEntry.path();
if (filepath.extension() == ".json") {
SimulationJob job;
job.load(filepath.string());
reducedPatternSimulationJobs.push_back(
std::make_shared<SimulationJob>(job));
}
}
}
}
#if POLYSCOPE_DEFINED
void draw() const {
initPolyscope();
FormFinder simulator;
assert(fullPatternSimulationJobs.size() ==
reducedPatternSimulationJobs.size());
fullPatternSimulationJobs[0]->pMesh->registerForDrawing();
reducedPatternSimulationJobs[0]->pMesh->registerForDrawing();
const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
for (int simulationJobIndex = 0;
simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
// Drawing of full pattern results
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
fullPatternSimulationJobs[simulationJobIndex];
pFullPatternSimulationJob->registerForDrawing(
fullPatternSimulationJobs[0]->pMesh->getLabel());
SimulationResults fullModelResults =
simulator.executeSimulation(pFullPatternSimulationJob);
fullModelResults.registerForDrawing();
// Drawing of reduced pattern results
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
reducedPatternSimulationJobs[simulationJobIndex];
SimulationResults reducedModelResults =
simulator.executeSimulation(pReducedPatternSimulationJob);
reducedModelResults.registerForDrawing();
polyscope::show();
// Save a screensh
// const std::string screenshotFilename =
// "/home/iason/Coding/Projects/Approximating shapes with flat "
// "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
// + m_pFullPatternSimulationMesh->getLabel() + "_" +
// simulationScenarioStrings[simulationScenarioIndex];
// polyscope::screenshot(screenshotFilename, false);
fullModelResults.unregister();
reducedModelResults.unregister();
// double error = computeError(
// reducedModelResults,
// global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);
// std::cout << "Error of simulation scenario "
// << simulationScenarioStrings[simulationScenarioIndex] << "
// is "
// << error << std::endl;
}
}
#endif // POLYSCOPE_DEFINED
void saveMeshFiles() const {
const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
assert(numberOfSimulationJobs != 0 &&
fullPatternSimulationJobs.size() ==
reducedPatternSimulationJobs.size());
fullPatternSimulationJobs[0]->pMesh->savePly("FullPattern_undeformed.ply");
reducedPatternSimulationJobs[0]->pMesh->savePly(
"ReducedPattern_undeformed.ply");
FormFinder simulator;
for (int simulationJobIndex = 0;
simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
// Drawing of full pattern results
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
fullPatternSimulationJobs[simulationJobIndex];
SimulationResults fullModelResults =
simulator.executeSimulation(pFullPatternSimulationJob);
fullModelResults.saveDeformedModel();
// Drawing of reduced pattern results
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
reducedPatternSimulationJobs[simulationJobIndex];
SimulationResults reducedModelResults =
simulator.executeSimulation(pReducedPatternSimulationJob);
reducedModelResults.saveDeformedModel();
}
}
void
writeHeaderTo(const ReducedModelOptimizer::Settings &settings_optimization,
csvFile &os) {
os << "Obj value";
for (const ReducedModelOptimizer::xRange &range :
settings_optimization.xRanges) {
os << range.label;
}
os << "Time(s)";
os << "#Crashes";
}
void
writeResultsTo(const ReducedModelOptimizer::Settings &settings_optimization,
csvFile &os) const {
os << objectiveValue;
for (const double &optimalX : x) {
os << optimalX;
}
for (int unusedXVarCounter = 0;
unusedXVarCounter < settings_optimization.xRanges.size() - x.size();
unusedXVarCounter++) {
os << "-";
}
os << time;
if (numberOfSimulationCrashes == 0) {
os << "-";
} else {
os << numberOfSimulationCrashes;
}
}
};
#endif // REDUCEDMODELOPTIMIZER_HPP