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

Parallelization on the set of valid patterns.

This commit is contained in:
Iason 2021-02-01 16:10:24 +02:00
parent de83cb3d67
commit 858282c859
3 changed files with 232 additions and 205 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

85
src/main.cpp
View File

@ -11,10 +11,10 @@
#include <chrono> #include <chrono>
#include <filesystem> #include <filesystem>
#include <iostream> #include <iostream>
#include <iterator>
#include <stdexcept> #include <stdexcept>
#include <string> #include <string>
#include <vcg/complex/algorithms/update/position.h> #include <vcg/complex/algorithms/update/position.h>
#include <iterator>
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
@ -43,11 +43,35 @@ int main(int argc, char *argv[]) {
std::vector<FlatPattern *> reducedModels{&singleBarReducedModel, std::vector<FlatPattern *> reducedModels{&singleBarReducedModel,
&CWReducedModel, &CCWReducedModel}; &CWReducedModel, &CCWReducedModel};
ReducedModelOptimizer optimizer(numberOfNodesPerSlot); // Test set of full patterns
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";
std::vector<std::pair<FlatPattern *, FlatPattern *>> patternPairs;
for (const auto &entry :
filesystem::directory_iterator(fullPatternsTestSetDirectory)) {
const auto filepath =
// std::filesystem::path(fullPatternsTestSetDirectory).append("305.ply");
entry.path();
const std::string filepathString = filepath.string();
const std::string tiledSuffix = "_tiled.ply";
if (filepathString.compare(filepathString.size() - tiledSuffix.size(),
tiledSuffix.size(), tiledSuffix) == 0) {
continue;
}
FlatPattern *pFullPattern = new FlatPattern(filepathString);
pFullPattern->setLabel(filepath.stem().string());
pFullPattern->scale(0.03);
FlatPattern *pReducedPattern = new FlatPattern();
pReducedPattern->copy(*reducedModels[0]);
patternPairs.push_back(std::make_pair(pFullPattern, pReducedPattern));
}
// for (double rangeOffset = 0.15; rangeOffset <= 0.95; rangeOffset += 0.05) // for (double rangeOffset = 0.15; rangeOffset <= 0.95; rangeOffset += 0.05)
// { // {
ReducedModelOptimizer::Settings settings; ReducedModelOptimizer::Settings settings;
for (settings.maxSimulations = 2200; settings.maxSimulations < 5000; for (settings.maxSimulations = 2600; settings.maxSimulations < 5000;
settings.maxSimulations += 200) { settings.maxSimulations += 200) {
ReducedModelOptimizer::xRange beamWidth{"B", 0.5, 1.5}; ReducedModelOptimizer::xRange beamWidth{"B", 0.5, 1.5};
ReducedModelOptimizer::xRange beamDimensionsRatio{"bOverh", 0.7, 1.3}; ReducedModelOptimizer::xRange beamDimensionsRatio{"bOverh", 0.7, 1.3};
@ -68,40 +92,17 @@ int main(int argc, char *argv[]) {
double totalError = 0; double totalError = 0;
int totalNumberOfSimulationCrashes = 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";
std::vector<std::pair<std::string, ReducedModelOptimizer::Results>> std::vector<std::pair<std::string, ReducedModelOptimizer::Results>>
resultsPerPattern; resultsPerPattern(patternPairs.size());
auto start = std::chrono::high_resolution_clock::now(); auto start = std::chrono::high_resolution_clock::now();
int patternsOptimized = 0;
//extract paths #pragma omp parallel for // schedule(static) num_threads(8)
std::vector<std::filesystem::path> fullPatternSetFilenames; for (int patternPairIndex = 0; patternPairIndex < patternPairs.size();
for (const auto& entry : patternPairIndex++) {
filesystem::directory_iterator(fullPatternsTestSetDirectory)) { FlatPattern *pPattern = patternPairs[patternPairIndex].first;
const auto filepath = // const auto filepathString = filepath.string();
// std::filesystem::path(fullPatternsTestSetDirectory).append("305.ply");
entry.path();
fullPatternSetFilenames.push_back(filepath);
}
//#pragma omp parallel for //schedule(static) num_threads(8)
for (int fullPatternIndex = 0; fullPatternIndex < fullPatternSetFilenames.size();fullPatternIndex++) {
const auto& filepath = fullPatternSetFilenames[fullPatternIndex];
const std::string& filepathString = fullPatternSetFilenames[fullPatternIndex].string();
//const auto filepathString = filepath.string();
// Use only the base triangle version // 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; // std::cout << "Full pattern:" << filepathString << std::endl;
FlatPattern pattern(filepathString);
pattern.setLabel(filepath.stem().string());
pattern.scale(0.03);
// for (int reducedPatternIndex = 0; // for (int reducedPatternIndex = 0;
// reducedPatternIndex < reducedModels.size(); // reducedPatternIndex < reducedModels.size();
// reducedPatternIndex++) { // reducedPatternIndex++) {
@ -117,7 +118,10 @@ int main(int argc, char *argv[]) {
// } // }
// FlatPattern cp; // FlatPattern cp;
// cp.copy(*reducedModels[0]); // cp.copy(*reducedModels[0]);
optimizer.initializePatterns(pattern, *reducedModels[0], const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
optimizer.initializePatterns(*pPattern,
*patternPairs[patternPairIndex].second,
optimizationExcludedEi); optimizationExcludedEi);
// optimizer.optimize({ReducedModelOptimizer::Axial}); // optimizer.optimize({ReducedModelOptimizer::Axial});
ReducedModelOptimizer::Results optimizationResults = ReducedModelOptimizer::Results optimizationResults =
@ -139,11 +143,16 @@ int main(int argc, char *argv[]) {
// thisOptimizationStatistics << endrow; // thisOptimizationStatistics << endrow;
totalError += optimizationResults.objectiveValue; totalError += optimizationResults.objectiveValue;
resultsPerPattern.push_back( resultsPerPattern[patternPairIndex] =
std::make_pair(filepath.stem().string(), optimizationResults)); std::make_pair(pPattern->getLabel(), optimizationResults);
totalNumberOfSimulationCrashes += totalNumberOfSimulationCrashes +=
optimizationResults.numberOfSimulationCrashes; optimizationResults.numberOfSimulationCrashes;
std::cout << "Have optimized " <<++patternsOptimized<<"/"<< static_cast<int>(std::distance(std::filesystem::directory_iterator(fullPatternsTestSetDirectory),std::filesystem::directory_iterator()))<<" patterns." << std::endl; // std::cout << "Have optimized " << ++patternsOptimized << "/"
// << static_cast<int>(
// std::distance(std::filesystem::directory_iterator(
// fullPatternsTestSetDirectory),
// std::filesystem::directory_iterator()))
// << " patterns." << std::endl;
// } // }
} }
auto end = std::chrono::high_resolution_clock::now(); auto end = std::chrono::high_resolution_clock::now();
@ -177,5 +186,9 @@ int main(int argc, char *argv[]) {
statistics << endrow; statistics << endrow;
} }
for(auto patternPair:patternPairs){
delete patternPair.first;
delete patternPair.second;
}
return 0; return 0;
} }

351
src/reducedmodeloptimizer.cpp
View File

@ -7,35 +7,35 @@
#include <dlib/global_optimization.h> #include <dlib/global_optimization.h>
#include <dlib/optimization.h> #include <dlib/optimization.h>
const bool gShouldDraw = true; struct GlobalOptimizationVariables {
const bool gShouldDraw = true;
FormFinder simulator; std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements; std::vector<SimulationJob> g_fullPatternSimulationJob;
std::vector<SimulationJob> g_fullPatternSimulationJob; std::vector<std::shared_ptr<SimulationJob>> g_reducedPatternSimulationJob;
std::vector<std::shared_ptr<SimulationJob>> g_reducedPatternSimulationJob; std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> reducedToFullInterfaceViMap;
reducedToFullInterfaceViMap; std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> g_reducedToFullViMap;
g_reducedToFullViMap; matplot::line_handle gPlotHandle;
matplot::line_handle gPlotHandle; std::vector<double> gObjectiveValueHistory;
std::vector<double> gObjectiveValueHistory; Eigen::Vector2d g_initialX;
Eigen::Vector2d g_initialX; std::unordered_set<size_t> g_reducedPatternExludedEdges;
std::unordered_set<size_t> g_reducedPatternExludedEdges; // double g_initialParameters;
// double g_initialParameters; Eigen::VectorXd g_initialParameters;
Eigen::VectorXd g_initialParameters; std::vector<ReducedModelOptimizer::SimulationScenario>
std::vector<ReducedModelOptimizer::SimulationScenario> g_simulationScenarioIndices;
g_simulationScenarioIndices; std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)};
std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)}; double g_innerHexagonInitialPos = 0;
double g_innerHexagonInitialPos = 0; bool g_optimizeInnerHexagonSize{false};
bool g_optimizeInnerHexagonSize{false}; std::vector<SimulationResults> firstOptimizationRoundResults;
std::vector<SimulationResults> firstOptimizationRoundResults; int g_firstRoundIterationIndex{0};
int g_firstRoundIterationIndex{0}; double minY{std::numeric_limits<double>::max()};
double minY{std::numeric_limits<double>::max()}; std::vector<double> minX;
std::vector<double> minX; std::vector<std::vector<double>> failedSimulationsXRatio;
std::vector<std::vector<double>> failedSimulationsXRatio; int numOfSimulationCrashes{false};
int numOfSimulationCrashes{false}; int numberOfFunctionCalls{0};
int numberOfFunctionCalls{ 0 }; } global;
#pragma omp threadprivate(global)
// struct OptimizationCallback { // struct OptimizationCallback {
// double operator()(const size_t &iterations, const Eigen::VectorXd &x, // double operator()(const size_t &iterations, const Eigen::VectorXd &x,
// const double &fval, Eigen::VectorXd &gradient) const { // const double &fval, Eigen::VectorXd &gradient) const {
@ -117,7 +117,7 @@ double ReducedModelOptimizer::computeError(
const SimulationResults &reducedPatternResults, const SimulationResults &reducedPatternResults,
const Eigen::MatrixX3d &optimalReducedPatternDisplacements) { const Eigen::MatrixX3d &optimalReducedPatternDisplacements) {
double error = 0; double error = 0;
for (const auto reducedFullViPair : g_reducedToFullViMap) { for (const auto reducedFullViPair : global.g_reducedToFullViMap) {
VertexIndex reducedModelVi = reducedFullViPair.first; VertexIndex reducedModelVi = reducedFullViPair.first;
// const auto pos = // const auto pos =
// g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP(); // g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP();
@ -144,7 +144,9 @@ double ReducedModelOptimizer::computeError(
void updateMesh(long n, const double *x) { void updateMesh(long n, const double *x) {
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh = std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]->pMesh; global
.g_reducedPatternSimulationJob[global.g_simulationScenarioIndices[0]]
->pMesh;
// const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0]; // const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0];
// std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << " // std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << "
// " // "
@ -158,11 +160,12 @@ void updateMesh(long n, const double *x) {
// e.properties.E = g_initialParameters * x[ei]; // e.properties.E = g_initialParameters * x[ei];
// e.properties.E = g_initialParameters(0) * x[0]; // e.properties.E = g_initialParameters(0) * x[0];
// e.properties.G = g_initialParameters(1) * x[1]; // e.properties.G = g_initialParameters(1) * x[1];
e.setDimensions(RectangularBeamDimensions( e.setDimensions(
g_initialParameters(0) * x[0], RectangularBeamDimensions(global.g_initialParameters(0) * x[0],
g_initialParameters(0) * x[0] / (g_initialParameters(1) * x[1]))); global.g_initialParameters(0) * x[0] /
(global.g_initialParameters(1) * x[1])));
e.setMaterial(ElementMaterial(e.material.poissonsRatio, e.setMaterial(ElementMaterial(e.material.poissonsRatio,
g_initialParameters(2) * x[2])); global.g_initialParameters(2) * x[2]));
// e.properties.A = g_initialParameters(0) * x[0]; // e.properties.A = g_initialParameters(0) * x[0];
// e.properties.J = g_initialParameters(1) * x[1]; // e.properties.J = g_initialParameters(1) * x[1];
// e.properties.I2 = g_initialParameters(2) * x[2]; // e.properties.I2 = g_initialParameters(2) * x[2];
@ -186,12 +189,12 @@ void updateMesh(long n, const double *x) {
// e.secondBendingConstFactor // e.secondBendingConstFactor
// << std::endl; // << std::endl;
if (g_optimizeInnerHexagonSize) { if (global.g_optimizeInnerHexagonSize) {
assert(pReducedPatternSimulationMesh->EN() == 12); assert(pReducedPatternSimulationMesh->EN() == 12);
for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN(); for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN();
vi += 2) { vi += 2) {
pReducedPatternSimulationMesh->vert[vi].P() = pReducedPatternSimulationMesh->vert[vi].P() =
g_innerHexagonVectors[vi / 2] * x[n - 1]; global.g_innerHexagonVectors[vi / 2] * x[n - 1];
} }
pReducedPatternSimulationMesh->reset(); pReducedPatternSimulationMesh->reset();
@ -219,7 +222,8 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
// std::cout << "x[" + std::to_string(parameterIndex) + "]=" // std::cout << "x[" + std::to_string(parameterIndex) + "]="
// << x[parameterIndex] << std::endl; // << x[parameterIndex] << std::endl;
// } // }
const Element &e = g_reducedPatternSimulationJob[0]->pMesh->elements[0]; const Element &e =
global.g_reducedPatternSimulationJob[0]->pMesh->elements[0];
// std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
// << e.firstBendingConstFactor << " " << // << e.firstBendingConstFactor << " " <<
// e.secondBendingConstFactor // e.secondBendingConstFactor
@ -232,11 +236,12 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
// run simulations // run simulations
double error = 0; double error = 0;
FormFinder simulator;
FormFinder::Settings simulationSettings; FormFinder::Settings simulationSettings;
// simulationSettings.shouldDraw = true; // simulationSettings.shouldDraw = true;
for (const int simulationScenarioIndex : g_simulationScenarioIndices) { for (const int simulationScenarioIndex : global.g_simulationScenarioIndices) {
SimulationResults reducedModelResults = simulator.executeSimulation( SimulationResults reducedModelResults = simulator.executeSimulation(
g_reducedPatternSimulationJob[simulationScenarioIndex], global.g_reducedPatternSimulationJob[simulationScenarioIndex],
simulationSettings); simulationSettings);
std::string filename; std::string filename;
if (!reducedModelResults.converged /*&& if (!reducedModelResults.converged /*&&
@ -267,18 +272,21 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
// simulator.executeSimulation( // simulator.executeSimulation(
// g_reducedPatternSimulationJob[simulationScenarioIndex], // g_reducedPatternSimulationJob[simulationScenarioIndex],
// simulationSettings); // simulationSettings);
numOfSimulationCrashes++; global.numOfSimulationCrashes++;
} else { } else {
error += computeError( error += computeError(
reducedModelResults, reducedModelResults,
g_optimalReducedModelDisplacements[simulationScenarioIndex]); global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);
filename = "/home/iason/Coding/Projects/Approximating shapes with flat " filename = "/home/iason/Coding/Projects/Approximating shapes with flat "
"patterns/RodModelOptimizationForPatterns/build/" "patterns/RodModelOptimizationForPatterns/build/"
"ProblematicSimulationJobs/conv_dimensions.txt"; "ProblematicSimulationJobs/conv_dimensions.txt";
} }
std::ofstream out(filename, std::ios_base::app); std::ofstream out(filename, std::ios_base::app);
auto pMesh = auto pMesh =
g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]->pMesh; global
.g_reducedPatternSimulationJob[global
.g_simulationScenarioIndices[0]]
->pMesh;
for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
out << "x[" + std::to_string(parameterIndex) + "]=" << x[parameterIndex] out << "x[" + std::to_string(parameterIndex) + "]=" << x[parameterIndex]
@ -296,16 +304,17 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
"\n\n"; "\n\n";
out.close(); out.close();
} }
if (error < minY) { if (error < global.minY) {
minY = error; global.minY = error;
minX.assign(x, x + n); global.minX.assign(x, x + n);
} }
if(++numberOfFunctionCalls%50==0){ if (++global.numberOfFunctionCalls % 50 == 0) {
std::cout << "Number of function calls:"<<numberOfFunctionCalls << std::endl; std::cout << "Number of function calls:" << global.numberOfFunctionCalls
<< std::endl;
} }
// compute error and return it // compute error and return it
gObjectiveValueHistory.push_back(error); global.gObjectiveValueHistory.push_back(error);
// auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), // auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
// gObjectiveValueHistory.size()); // gObjectiveValueHistory.size());
// std::vector<double> colors(gObjectiveValueHistory.size(), 2); // std::vector<double> colors(gObjectiveValueHistory.size(), 2);
@ -349,18 +358,18 @@ void ReducedModelOptimizer::computeMaps(
// std::endl; // std::endl;
// Save excluded edges // Save excluded edges
g_reducedPatternExludedEdges.clear(); global.g_reducedPatternExludedEdges.clear();
const size_t fanSize = 6; const size_t fanSize = 6;
const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN(); const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN();
for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
for (const size_t ei : reducedModelExcludedEges) { for (const size_t ei : reducedModelExcludedEges) {
g_reducedPatternExludedEdges.insert( global.g_reducedPatternExludedEdges.insert(
fanIndex * reducedBaseTriangleNumberOfEdges + ei); fanIndex * reducedBaseTriangleNumberOfEdges + ei);
} }
} }
// Construct reduced->full and full->reduced interface vi map // Construct reduced->full and full->reduced interface vi map
reducedToFullInterfaceViMap.clear(); global.reducedToFullInterfaceViMap.clear();
vcg::tri::Allocator<FlatPattern>::PointerUpdater<FlatPattern::VertexPointer> vcg::tri::Allocator<FlatPattern>::PointerUpdater<FlatPattern::VertexPointer>
pu_reducedModel; pu_reducedModel;
reducedPattern.deleteDanglingVertices(pu_reducedModel); reducedPattern.deleteDanglingVertices(pu_reducedModel);
@ -370,13 +379,14 @@ void ReducedModelOptimizer::computeMaps(
reducedPattern.VN() - 1 /*- reducedModelBaseTriangleInterfaceVi*/; reducedPattern.VN() - 1 /*- reducedModelBaseTriangleInterfaceVi*/;
reducedPattern.createFan(); reducedPattern.createFan();
for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex + global.reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset *
reducedModelBaseTriangleInterfaceVi] = fanIndex +
reducedModelBaseTriangleInterfaceVi] =
fullModelBaseTriangleInterfaceVi + fullModelBaseTriangleInterfaceVi +
fanIndex * fullPatternInterfaceVertexOffset; fanIndex * fullPatternInterfaceVertexOffset;
} }
m_fullToReducedInterfaceViMap.clear(); m_fullToReducedInterfaceViMap.clear();
constructInverseMap(reducedToFullInterfaceViMap, constructInverseMap(global.reducedToFullInterfaceViMap,
m_fullToReducedInterfaceViMap); m_fullToReducedInterfaceViMap);
// fullPattern.setLabel("FullPattern"); // fullPattern.setLabel("FullPattern");
@ -391,14 +401,14 @@ void ReducedModelOptimizer::computeMaps(
m_fullPatternOppositeInterfaceViMap[vi0] = vi1; m_fullPatternOppositeInterfaceViMap[vi0] = vi1;
} }
g_reducedToFullViMap = reducedToFullInterfaceViMap; global.g_reducedToFullViMap = global.reducedToFullInterfaceViMap;
const bool debugMapping = false; const bool debugMapping = false;
if (debugMapping) { if (debugMapping) {
reducedPattern.registerForDrawing(); reducedPattern.registerForDrawing();
std::vector<glm::vec3> colors_reducedPatternExcludedEdges( std::vector<glm::vec3> colors_reducedPatternExcludedEdges(
reducedPattern.EN(), glm::vec3(0, 0, 0)); reducedPattern.EN(), glm::vec3(0, 0, 0));
for (const size_t ei : g_reducedPatternExludedEdges) { for (const size_t ei : global.g_reducedPatternExludedEdges) {
colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0); colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0);
} }
const std::string label = reducedPattern.getLabel(); const std::string label = reducedPattern.getLabel();
@ -427,11 +437,12 @@ void ReducedModelOptimizer::computeMaps(
std::vector<glm::vec3> nodeColorsReducedToFull_full(fullPattern.VN(), std::vector<glm::vec3> nodeColorsReducedToFull_full(fullPattern.VN(),
glm::vec3(0, 0, 0)); glm::vec3(0, 0, 0));
for (size_t vi = 0; vi < reducedPattern.VN(); vi++) { for (size_t vi = 0; vi < reducedPattern.VN(); vi++) {
if (reducedToFullInterfaceViMap.contains(vi)) { if (global.reducedToFullInterfaceViMap.contains(vi)) {
auto color = polyscope::getNextUniqueColor(); auto color = polyscope::getNextUniqueColor();
nodeColorsReducedToFull_reduced[vi] = color; nodeColorsReducedToFull_reduced[vi] = color;
nodeColorsReducedToFull_full[reducedToFullInterfaceViMap[vi]] = color; nodeColorsReducedToFull_full[global.reducedToFullInterfaceViMap[vi]] =
color;
} }
} }
polyscope::getCurveNetwork(reducedPattern.getLabel()) polyscope::getCurveNetwork(reducedPattern.getLabel())
@ -482,8 +493,8 @@ void ReducedModelOptimizer::initializePatterns(
FlatPattern copyReducedPattern; FlatPattern copyReducedPattern;
copyFullPattern.copy(fullPattern); copyFullPattern.copy(fullPattern);
copyReducedPattern.copy(reducedPattern); copyReducedPattern.copy(reducedPattern);
g_optimizeInnerHexagonSize = copyReducedPattern.EN() == 2; global.g_optimizeInnerHexagonSize = copyReducedPattern.EN() == 2;
if (g_optimizeInnerHexagonSize) { if (global.g_optimizeInnerHexagonSize) {
const double h = copyReducedPattern.getBaseTriangleHeight(); const double h = copyReducedPattern.getBaseTriangleHeight();
double baseTriangle_bottomEdgeSize = 2 * h / tan(vcg::math::ToRad(60.0)); double baseTriangle_bottomEdgeSize = 2 * h / tan(vcg::math::ToRad(60.0));
VectorType baseTriangle_leftBottomNode(-baseTriangle_bottomEdgeSize / 2, -h, VectorType baseTriangle_leftBottomNode(-baseTriangle_bottomEdgeSize / 2, -h,
@ -497,13 +508,13 @@ void ReducedModelOptimizer::initializePatterns(
vcg::RotationMatrix(rotationAxis, vcg::RotationMatrix(rotationAxis,
vcg::math::ToRad(rotationCounter * 60.0)) * vcg::math::ToRad(rotationCounter * 60.0)) *
baseTriangle_leftBottomNode; baseTriangle_leftBottomNode;
g_innerHexagonVectors[rotationCounter] = rotatedVector; global.g_innerHexagonVectors[rotationCounter] = rotatedVector;
} }
const double innerHexagonInitialPos_x = const double innerHexagonInitialPos_x =
copyReducedPattern.vert[0].cP()[0] / g_innerHexagonVectors[0][0]; copyReducedPattern.vert[0].cP()[0] / global.g_innerHexagonVectors[0][0];
const double innerHexagonInitialPos_y = const double innerHexagonInitialPos_y =
copyReducedPattern.vert[0].cP()[1] / g_innerHexagonVectors[0][1]; copyReducedPattern.vert[0].cP()[1] / global.g_innerHexagonVectors[0][1];
g_innerHexagonInitialPos = innerHexagonInitialPos_x; global.g_innerHexagonInitialPos = innerHexagonInitialPos_x;
} }
computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges); computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges);
createSimulationMeshes(copyFullPattern, copyReducedPattern); createSimulationMeshes(copyFullPattern, copyReducedPattern);
@ -513,9 +524,9 @@ void ReducedModelOptimizer::initializePatterns(
void ReducedModelOptimizer::initializeOptimizationParameters( void ReducedModelOptimizer::initializeOptimizationParameters(
const std::shared_ptr<SimulationMesh> &mesh) { const std::shared_ptr<SimulationMesh> &mesh) {
const int numberOfOptimizationParameters = 3; const int numberOfOptimizationParameters = 3;
g_initialParameters.resize(g_optimizeInnerHexagonSize global.g_initialParameters.resize(global.g_optimizeInnerHexagonSize
? numberOfOptimizationParameters + 1 ? numberOfOptimizationParameters + 1
: numberOfOptimizationParameters); : numberOfOptimizationParameters);
// Save save the beam stiffnesses // Save save the beam stiffnesses
// for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) { // for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) {
// Element &e = pReducedModelElementalMesh->elements[ei]; // Element &e = pReducedModelElementalMesh->elements[ei];
@ -529,9 +540,9 @@ void ReducedModelOptimizer::initializeOptimizationParameters(
const double initialB = std::sqrt(mesh->elements[0].A); const double initialB = std::sqrt(mesh->elements[0].A);
const double initialRatio = 1; const double initialRatio = 1;
; ;
g_initialParameters(0) = initialB; global.g_initialParameters(0) = initialB;
g_initialParameters(1) = initialRatio; global.g_initialParameters(1) = initialRatio;
g_initialParameters(2) = mesh->elements[0].material.youngsModulus; global.g_initialParameters(2) = mesh->elements[0].material.youngsModulus;
// g_initialParameters = // g_initialParameters =
// m_pReducedPatternSimulationMesh->elements[0].properties.E; // m_pReducedPatternSimulationMesh->elements[0].properties.E;
// for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) { // for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) {
@ -605,9 +616,9 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(
const Settings &settings, const Settings &settings,
double (*pObjectiveFunction)(long, const double *)) { double (*pObjectiveFunction)(long, const double *)) {
gObjectiveValueHistory.clear(); global.gObjectiveValueHistory.clear();
const size_t n = g_initialParameters.rows(); const size_t n = global.g_initialParameters.rows();
assert(n == 3); assert(n == 3);
// g_optimizeInnerHexagonSize ? 5: 4; // g_optimizeInnerHexagonSize ? 5: 4;
// const size_t npt = (n + 1) * (n + 2) / 2; // const size_t npt = (n + 1) * (n + 2) / 2;
@ -619,8 +630,8 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(
const size_t initialGuess = 1; const size_t initialGuess = 1;
std::vector<double> x(n, initialGuess); std::vector<double> x(n, initialGuess);
if (g_optimizeInnerHexagonSize) { if (global.g_optimizeInnerHexagonSize) {
x[n - 1] = g_innerHexagonInitialPos; x[n - 1] = global.g_innerHexagonInitialPos;
} }
/*if (!initialGuess.empty()) { /*if (!initialGuess.empty()) {
x = g_optimizationInitialGuess; x = g_optimizationInitialGuess;
@ -655,7 +666,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(
xMax(i) = settings.xRanges[i].max; xMax(i) = settings.xRanges[i].max;
} }
numberOfFunctionCalls = 0; global.numberOfFunctionCalls = 0;
double (*objF)(double, double, double) = &objective; double (*objF)(double, double, double) = &objective;
auto start = std::chrono::system_clock::now(); auto start = std::chrono::system_clock::now();
dlib::function_evaluation result = dlib::find_min_global( dlib::function_evaluation result = dlib::find_min_global(
@ -663,11 +674,11 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(
std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy); std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy);
auto end = std::chrono::system_clock::now(); auto end = std::chrono::system_clock::now();
auto elapsed = std::chrono::duration_cast<std::chrono::seconds>(end - start); auto elapsed = std::chrono::duration_cast<std::chrono::seconds>(end - start);
Results results{numOfSimulationCrashes, minX, minY}; Results results{global.numOfSimulationCrashes, global.minX, global.minY};
std::cout << "Finished optimizing." << endl; std::cout << "Finished optimizing." << endl;
// std::cout << "Solution x:" << endl; // std::cout << "Solution x:" << endl;
// std::cout << result.x << endl; // std::cout << result.x << endl;
std::cout << "Objective value:" << minY << endl; std::cout << "Objective value:" << global.minY << endl;
// std::cout << result.y << endl; // std::cout << result.y << endl;
// std::cout << minY << endl; // std::cout << minY << endl;
// std::cout << "Time(sec):" << elapsed.count() << std::endl; // std::cout << "Time(sec):" << elapsed.count() << std::endl;
@ -874,87 +885,89 @@ ReducedModelOptimizer::createScenarios(
return scenarios; return scenarios;
} }
void ReducedModelOptimizer::runBeamOptimization() { // void ReducedModelOptimizer::runBeamOptimization() {
// load beams // // load beams
VCGEdgeMesh fullBeam; // VCGEdgeMesh fullBeam;
fullBeam.loadPly("/home/iason/Models/simple_beam_model_10elem_1m.ply"); // fullBeam.loadPly("/home/iason/Models/simple_beam_model_10elem_1m.ply");
VCGEdgeMesh reducedBeam; // VCGEdgeMesh reducedBeam;
reducedBeam.loadPly("/home/iason/Models/simple_beam_model_4elem_1m.ply"); // reducedBeam.loadPly("/home/iason/Models/simple_beam_model_4elem_1m.ply");
fullBeam.registerForDrawing(); // fullBeam.registerForDrawing();
reducedBeam.registerForDrawing(); // reducedBeam.registerForDrawing();
// polyscope::show(); // // polyscope::show();
// maps // // maps
std::unordered_map<size_t, size_t> displacementReducedToFullMap; // std::unordered_map<size_t, size_t> displacementReducedToFullMap;
displacementReducedToFullMap[reducedBeam.VN() / 2] = fullBeam.VN() / 2; // displacementReducedToFullMap[reducedBeam.VN() / 2] = fullBeam.VN() / 2;
g_reducedToFullViMap = displacementReducedToFullMap; // g_reducedToFullViMap = displacementReducedToFullMap;
std::unordered_map<size_t, size_t> jobFullToReducedMap; // std::unordered_map<size_t, size_t> jobFullToReducedMap;
jobFullToReducedMap[0] = 0; // jobFullToReducedMap[0] = 0;
jobFullToReducedMap[fullBeam.VN() - 1] = reducedBeam.VN() - 1; // jobFullToReducedMap[fullBeam.VN() - 1] = reducedBeam.VN() - 1;
// full model simuilation job // // full model simuilation job
auto pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullBeam); // auto pFullPatternSimulationMesh =
pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02}); // std::make_shared<SimulationMesh>(fullBeam);
pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); // pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.02,
std::unordered_map<VertexIndex, std::unordered_set<int>> fixedVertices; // 0.02}); pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
fixedVertices[0] = ::unordered_set<int>({0, 1, 2, 3, 4, 5}); // std::unordered_map<VertexIndex, std::unordered_set<int>> fixedVertices;
std::unordered_map<VertexIndex, Vector6d> forces; // fixedVertices[0] = ::unordered_set<int>({0, 1, 2, 3, 4, 5});
forces[fullBeam.VN() - 1] = Vector6d({0, 0, 10, 0, 0, 0}); // std::unordered_map<VertexIndex, Vector6d> forces;
const std::string fullBeamSimulationJobLabel = "Pull_Z"; // forces[fullBeam.VN() - 1] = Vector6d({0, 0, 10, 0, 0, 0});
std::shared_ptr<SimulationJob> pFullModelSimulationJob = // const std::string fullBeamSimulationJobLabel = "Pull_Z";
make_shared<SimulationJob>(SimulationJob(pFullPatternSimulationMesh, // std::shared_ptr<SimulationJob> pFullModelSimulationJob =
fullBeamSimulationJobLabel, // make_shared<SimulationJob>(SimulationJob(pFullPatternSimulationMesh,
fixedVertices, forces)); // fullBeamSimulationJobLabel,
FormFinder formFinder; // fixedVertices, forces));
auto fullModelResults = formFinder.executeSimulation(pFullModelSimulationJob); // auto fullModelResults =
// formFinder.executeSimulation(pFullModelSimulationJob);
// Optimal reduced model displacements // // Optimal reduced model displacements
const size_t numberOfSimulationScenarios = 1; // const size_t numberOfSimulationScenarios = 1;
g_optimalReducedModelDisplacements.resize(numberOfSimulationScenarios); // g_optimalReducedModelDisplacements.resize(numberOfSimulationScenarios);
g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1].resize( // g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1].resize(
reducedBeam.VN(), 3); // reducedBeam.VN(), 3);
computeDesiredReducedModelDisplacements( // computeDesiredReducedModelDisplacements(
fullModelResults, displacementReducedToFullMap, // fullModelResults, displacementReducedToFullMap,
g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]); // g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]);
// reduced model simuilation job // // reduced model simuilation job
auto reducedSimulationMesh = std::make_shared<SimulationMesh>(reducedBeam); // auto reducedSimulationMesh = std::make_shared<SimulationMesh>(reducedBeam);
reducedSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02}); // reducedSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02});
reducedSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); // reducedSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
g_reducedPatternSimulationJob.resize(numberOfSimulationScenarios); // g_reducedPatternSimulationJob.resize(numberOfSimulationScenarios);
SimulationJob reducedSimJob; // SimulationJob reducedSimJob;
computeReducedModelSimulationJob(*pFullModelSimulationJob, // computeReducedModelSimulationJob(*pFullModelSimulationJob,
jobFullToReducedMap, reducedSimJob); // jobFullToReducedMap, reducedSimJob);
reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] = // reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] =
reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] * 0.1; // reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] * 0.1;
g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1] = // g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1] =
make_shared<SimulationJob>( // make_shared<SimulationJob>(
reducedSimulationMesh, fullBeamSimulationJobLabel, // reducedSimulationMesh, fullBeamSimulationJobLabel,
reducedSimJob.constrainedVertices, reducedSimJob.nodalExternalForces); // reducedSimJob.constrainedVertices,
initializeOptimizationParameters(reducedSimulationMesh); // reducedSimJob.nodalExternalForces);
// initializeOptimizationParameters(reducedSimulationMesh);
// const std::string simulationJobsPath = "SimulationJobs"; // // const std::string simulationJobsPath = "SimulationJobs";
// std::filesystem::create_directory(simulationJobsPath); // // std::filesystem::create_directory(simulationJobsPath);
// g_reducedPatternSimulationJob[0].save(simulationJobsPath); // // g_reducedPatternSimulationJob[0].save(simulationJobsPath);
// g_reducedPatternSimulationJob[0].load( // // g_reducedPatternSimulationJob[0].load(
// std::filesystem::path(simulationJobsPath) // // std::filesystem::path(simulationJobsPath)
// .append(g_reducedPatternSimulationJob[0].mesh->getLabel() + // // .append(g_reducedPatternSimulationJob[0].mesh->getLabel() +
// "_simScenario.json")); // // "_simScenario.json"));
runOptimization({}, &objective); // runOptimization({}, &objective);
fullModelResults.registerForDrawing(); // fullModelResults.registerForDrawing();
SimulationResults reducedModelResults = simulator.executeSimulation( // SimulationResults reducedModelResults = simulator.executeSimulation(
g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1]); // g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1]);
double error = computeError( // double error = computeError(
reducedModelResults, // reducedModelResults,
g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]); // g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]);
reducedModelResults.registerForDrawing(); // reducedModelResults.registerForDrawing();
std::cout << "Error between beams:" << error << endl; // std::cout << "Error between beams:" << error << endl;
// registerWorldAxes(); // // registerWorldAxes();
polyscope::show(); // polyscope::show();
fullModelResults.unregister(); // fullModelResults.unregister();
reducedModelResults.unregister(); // reducedModelResults.unregister();
} //}
void ReducedModelOptimizer::visualizeResults( void ReducedModelOptimizer::visualizeResults(
const std::vector<std::shared_ptr<SimulationJob>> const std::vector<std::shared_ptr<SimulationJob>>
@ -971,12 +984,12 @@ void ReducedModelOptimizer::visualizeResults(
fullModelResults.registerForDrawing(); fullModelResults.registerForDrawing();
fullModelResults.saveDeformedModel(); fullModelResults.saveDeformedModel();
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob = const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
g_reducedPatternSimulationJob[simulationScenarioIndex]; global.g_reducedPatternSimulationJob[simulationScenarioIndex];
SimulationResults reducedModelResults = SimulationResults reducedModelResults =
simulator.executeSimulation(pReducedPatternSimulationJob); simulator.executeSimulation(pReducedPatternSimulationJob);
double error = computeError( double error = computeError(
reducedModelResults, reducedModelResults,
g_optimalReducedModelDisplacements[simulationScenarioIndex]); global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);
std::cout << "Error of simulation scenario " std::cout << "Error of simulation scenario "
<< simulationScenarioStrings[simulationScenarioIndex] << " is " << simulationScenarioStrings[simulationScenarioIndex] << " is "
<< error << std::endl; << error << std::endl;
@ -1001,9 +1014,9 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
const Settings &xRanges, const Settings &xRanges,
const std::vector<SimulationScenario> &simulationScenarios) { const std::vector<SimulationScenario> &simulationScenarios) {
g_simulationScenarioIndices = simulationScenarios; global.g_simulationScenarioIndices = simulationScenarios;
if (g_simulationScenarioIndices.empty()) { if (global.g_simulationScenarioIndices.empty()) {
g_simulationScenarioIndices = { global.g_simulationScenarioIndices = {
SimulationScenario::Axial, SimulationScenario::Shear, SimulationScenario::Axial, SimulationScenario::Shear,
SimulationScenario::Bending, SimulationScenario::Dome, SimulationScenario::Bending, SimulationScenario::Dome,
SimulationScenario::Saddle}; SimulationScenario::Saddle};
@ -1011,32 +1024,32 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
std::vector<std::shared_ptr<SimulationJob>> simulationJobs = std::vector<std::shared_ptr<SimulationJob>> simulationJobs =
createScenarios(m_pFullPatternSimulationMesh); createScenarios(m_pFullPatternSimulationMesh);
g_optimalReducedModelDisplacements.resize(6); global.g_optimalReducedModelDisplacements.resize(6);
g_reducedPatternSimulationJob.resize(6); global.g_reducedPatternSimulationJob.resize(6);
g_firstRoundIterationIndex = 0; global.g_firstRoundIterationIndex = 0;
minY = std::numeric_limits<double>::max(); global.minY = std::numeric_limits<double>::max();
numOfSimulationCrashes = 0; global.numOfSimulationCrashes = 0;
// polyscope::removeAllStructures(); // polyscope::removeAllStructures();
FormFinder::Settings settings; FormFinder::Settings settings;
// settings.shouldDraw = true; // settings.shouldDraw = true;
for (int simulationScenarioIndex : g_simulationScenarioIndices) { for (int simulationScenarioIndex : global.g_simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob = const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
simulationJobs[simulationScenarioIndex]; simulationJobs[simulationScenarioIndex];
SimulationResults fullModelResults = SimulationResults fullModelResults =
simulator.executeSimulation(pFullPatternSimulationJob, settings); simulator.executeSimulation(pFullPatternSimulationJob, settings);
g_optimalReducedModelDisplacements[simulationScenarioIndex].resize( global.g_optimalReducedModelDisplacements[simulationScenarioIndex].resize(
m_pReducedPatternSimulationMesh->VN(), 3); m_pReducedPatternSimulationMesh->VN(), 3);
computeDesiredReducedModelDisplacements( computeDesiredReducedModelDisplacements(
fullModelResults, g_reducedToFullViMap, fullModelResults, global.g_reducedToFullViMap,
g_optimalReducedModelDisplacements[simulationScenarioIndex]); global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);
SimulationJob reducedPatternSimulationJob; SimulationJob reducedPatternSimulationJob;
reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
computeReducedModelSimulationJob(*pFullPatternSimulationJob, computeReducedModelSimulationJob(*pFullPatternSimulationJob,
m_fullToReducedInterfaceViMap, m_fullToReducedInterfaceViMap,
reducedPatternSimulationJob); reducedPatternSimulationJob);
g_reducedPatternSimulationJob[simulationScenarioIndex] = global.g_reducedPatternSimulationJob[simulationScenarioIndex] =
std::make_shared<SimulationJob>(reducedPatternSimulationJob); std::make_shared<SimulationJob>(reducedPatternSimulationJob);
} }
Results optResults = runOptimization(xRanges, &objective); Results optResults = runOptimization(xRanges, &objective);
updateMesh(optResults.x.size(), optResults.x.data()); updateMesh(optResults.x.size(), optResults.x.data());

1
src/reducedmodeloptimizer.hpp
View File

@ -106,6 +106,7 @@ private:
computeError(const SimulationResults &reducedPatternResults, computeError(const SimulationResults &reducedPatternResults,
const Eigen::MatrixX3d &optimalReducedPatternDisplacements); const Eigen::MatrixX3d &optimalReducedPatternDisplacements);
static double objective(long n, const double *x); static double objective(long n, const double *x);
FormFinder simulator;
}; };
#endif // REDUCEDMODELOPTIMIZER_HPP #endif // REDUCEDMODELOPTIMIZER_HPP