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

Splitted the update of the reduced pattern modular by splitting the geometry from the material parameters.

This commit is contained in:
iasonmanolas 2021-11-08 10:47:35 +02:00
parent b58ae0e13b
commit b763ca92e7
3 changed files with 369 additions and 168 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
src/main.cpp
View File

@ -38,7 +38,7 @@ int main(int argc, char *argv[]) {
reducedPattern.scale(0.03, interfaceNodeIndex); reducedPattern.scale(0.03, interfaceNodeIndex);
// Set the optization settings // Set the optization settings
ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000}; ReducedPatternOptimization::xRange beamE{"E", 0.1, 100};
ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000}; ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000};
ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000}; ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000};
ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000}; ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000};
@ -46,27 +46,38 @@ int main(int argc, char *argv[]) {
ReducedPatternOptimization::xRange innerHexagonSize{"HexSize", 0.05, 0.95}; ReducedPatternOptimization::xRange innerHexagonSize{"HexSize", 0.05, 0.95};
ReducedPatternOptimization::xRange innerHexagonAngle{"HexAngle", -30.0, 30.0}; ReducedPatternOptimization::xRange innerHexagonAngle{"HexAngle", -30.0, 30.0};
ReducedPatternOptimization::Settings settings_optimization; ReducedPatternOptimization::Settings settings_optimization;
settings_optimization.xRanges = {beamE,beamA,beamJ,beamI2,beamI3, // settings_optimization.xRanges
innerHexagonSize, innerHexagonAngle}; // = {beamE, beamA, beamJ, beamI2, beamI3, innerHexagonSize, innerHexagonAngle};
settings_optimization.xRanges = {beamE, beamA, beamI2, innerHexagonSize, innerHexagonAngle};
const bool input_numberOfFunctionCallsDefined = argc >= 4; const bool input_numberOfFunctionCallsDefined = argc >= 4;
settings_optimization.numberOfFunctionCalls = settings_optimization.numberOfFunctionCalls =
input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100; input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100;
settings_optimization.normalizationStrategy settings_optimization.normalizationStrategy
= ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon; = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
settings_optimization.normalizationParameter = 0.0003; settings_optimization.translationNormalizationParameter = 1e-3;
settings_optimization.solverAccuracy = 0.001; settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0);
// settings_optimization.translationNormalizationParameter = 1e-15;
// settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(1e-15);
// settings_optimization.solverAccuracy = 1e-3;
settings_optimization.solverAccuracy = 1e-3;
settings_optimization.objectiveWeights.translational = std::atof(argv[4]); settings_optimization.objectiveWeights.translational = std::atof(argv[4]);
settings_optimization.objectiveWeights.rotational = 2 - std::atof(argv[4]); settings_optimization.objectiveWeights.rotational = 2 - std::atof(argv[4]);
// Optimize pair std::string xConcatNames;
for (const auto &x : settings_optimization.xRanges) {
xConcatNames.append(x.label + "_");
}
xConcatNames.pop_back();
// Optimize pairthere
const std::string pairName = fullPattern.getLabel() + "@" + reducedPattern.getLabel(); const std::string pairName = fullPattern.getLabel() + "@" + reducedPattern.getLabel();
const std::string optimizationName = pairName + "(" const std::string optimizationName = pairName + "("
+ std::to_string(settings_optimization.numberOfFunctionCalls) + std::to_string(settings_optimization.numberOfFunctionCalls)
+ "_" + "_"
+ to_string_with_precision( + to_string_with_precision(
settings_optimization.objectiveWeights.translational) settings_optimization.objectiveWeights.translational)
+ ")"; + ")" + "_" + xConcatNames;
const bool input_resultDirectoryDefined = argc >= 6; const bool input_resultDirectoryDefined = argc >= 6;
const std::string optimizationResultsDirectory = input_resultDirectoryDefined const std::string optimizationResultsDirectory = input_resultDirectoryDefined
? argv[5] ? argv[5]
@ -91,13 +102,14 @@ int main(int argc, char *argv[]) {
} }
ReducedPatternOptimization::Results optimizationResults; ReducedPatternOptimization::Results optimizationResults;
bool optimizationAlreadyComputed = optimizationResultFolderExists; constexpr bool shouldReoptimize = true;
// bool optimizationAlreadyComputed = false; bool optimizationAlreadyComputed = false;
// if (optimizationResultFolderExists) { if (!shouldReoptimize && optimizationResultFolderExists) {
// const bool resultsWereSuccessfullyLoaded = optimizationResults.load(resultsOutputDir); const bool resultsWereSuccessfullyLoaded = optimizationResults.load(resultsOutputDir);
// if (resultsWereSuccessfullyLoaded && optimizationResults.settings == settings_optimization) { if (resultsWereSuccessfullyLoaded && optimizationResults.settings == settings_optimization) {
// } optimizationAlreadyComputed = true;
// } }
}
if (!optimizationAlreadyComputed) { if (!optimizationAlreadyComputed) {
auto start = std::chrono::system_clock::now(); auto start = std::chrono::system_clock::now();
@ -135,6 +147,11 @@ int main(int argc, char *argv[]) {
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
// optimizationResults.saveMeshFiles(); // optimizationResults.saveMeshFiles();
std::cout << "E:"
<< optimizationResults.reducedPatternSimulationJobs[0]
->pMesh->elements[0]
.material.youngsModulus
<< std::endl;
optimizationResults.draw(); optimizationResults.draw();
#endif #endif

471
src/reducedmodeloptimizer.cpp
View File

@ -52,6 +52,8 @@ double ReducedModelOptimizer::computeDisplacementError(
const double rawError = computeRawTranslationalError(fullPatternDisplacements, const double rawError = computeRawTranslationalError(fullPatternDisplacements,
reducedPatternDisplacements, reducedPatternDisplacements,
reducedToFullInterfaceViMap); reducedToFullInterfaceViMap);
// std::cout << "raw trans error:" << rawError << std::endl;
// std::cout << "raw trans error:" << normalizationFactor << std::endl;
return rawError / normalizationFactor; return rawError / normalizationFactor;
} }
@ -120,6 +122,9 @@ double ReducedModelOptimizer::computeError(
simulationResults_reducedPattern.displacements, simulationResults_reducedPattern.displacements,
reducedToFullInterfaceViMap, reducedToFullInterfaceViMap,
normalizationFactor_translationalDisplacement); normalizationFactor_translationalDisplacement);
// std::cout << "normalization factor:" << normalizationFactor_rotationalDisplacement << std::endl;
// std::cout << "trans error:" << translationalError << std::endl;
const double rotationalError const double rotationalError
= computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, = computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion,
simulationResults_reducedPattern.rotationalDisplacementQuaternion, simulationResults_reducedPattern.rotationalDisplacementQuaternion,
@ -129,18 +134,19 @@ double ReducedModelOptimizer::computeError(
+ global.optimizationSettings.objectiveWeights.rotational * rotationalError; + global.optimizationSettings.objectiveWeights.rotational * rotationalError;
} }
double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3, //double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3,
double innerHexagonSize, // double innerHexagonSize,
double innerHexagonRotationAngle) { // double innerHexagonRotationAngle) {
std::vector<double> x{E,A,J,I2,I3, innerHexagonSize, innerHexagonRotationAngle}; // std::vector<double> x{E,A,J,I2,I3, innerHexagonSize, innerHexagonRotationAngle};
return ReducedModelOptimizer::objective(x.size(), x.data()); // return ReducedModelOptimizer::objective(x.size(), x.data());
} //}
double ReducedModelOptimizer::objective(long n, const double *x) { double ReducedModelOptimizer::objective(const dlib::matrix<double, 0, 1> &x)
// std::cout.precision(17); {
// for (int i = 0; i < n; i++) { // std::cout.precision(17);
// std::cout << x[i] << " "; // for (int i = 0; i < x.size(); i++) {
// } // std::cout << x(i) << " ";
// }
// std::cout << std::endl; // std::cout << std::endl;
// std::cout << x[n - 2] << " " << x[n - 1] << std::endl; // std::cout << x[n - 2] << " " << x[n - 1] << std::endl;
@ -150,7 +156,10 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
// << e.firstBendingConstFactor << " " << // << e.firstBendingConstFactor << " " <<
// e.secondBendingConstFactor // e.secondBendingConstFactor
// << std::endl; // << std::endl;
updateMesh(n, x); const int n = x.size();
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh
= global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh;
function_updateReducedPattern(x, pReducedPatternSimulationMesh);
// global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(); // global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing();
// global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(); // global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing();
// polyscope::show(); // polyscope::show();
@ -168,6 +177,9 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
for (const int simulationScenarioIndex : global.simulationScenarioIndices) { for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &reducedJob const std::shared_ptr<SimulationJob> &reducedJob
= global.reducedPatternSimulationJobs[simulationScenarioIndex]; = global.reducedPatternSimulationJobs[simulationScenarioIndex];
//#ifdef POLYSCOPE_DEFINED
// std::cout << reducedJob->getLabel() << ":" << std::endl;
//#endif
SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob); SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob);
// std::string filename; // std::string filename;
if (!reducedModelResults.converged) { if (!reducedModelResults.converged) {
@ -177,7 +189,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
std::cout << "Failed simulation" << std::endl; std::cout << "Failed simulation" << std::endl;
#endif #endif
} else { } else {
const bool usePotentialEnergy = false; constexpr bool usePotentialEnergy = false;
double simulationScenarioError; double simulationScenarioError;
if (usePotentialEnergy) { if (usePotentialEnergy) {
simulationScenarioError = std::abs( simulationScenarioError = std::abs(
@ -191,17 +203,17 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
global.translationalDisplacementNormalizationValues[simulationScenarioIndex], global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
} }
#ifdef POLYSCOPE_DEFINED //#ifdef POLYSCOPE_DEFINED
// std::cout << reducedJob->getLabel() << " sim error:" << simulationScenarioError // reducedJob->pMesh->registerForDrawing(Colors::reducedInitial);
// << std::endl; // reducedModelResults.registerForDrawing(Colors::reducedDeformed);
// reducedJob->pMesh->registerForDrawing(Colors::reducedInitial); // global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed);
// reducedModelResults.registerForDrawing(Colors::reducedDeformed); // global.fullPatternResults[simulationScenarioIndex].registerForDrawing(
// global.fullPatternResults[simulationScenarioIndex].registerForDrawing( // Colors::fullDeformed);
// Colors::fullDeformed); // polyscope::show();
// polyscope::show(); // reducedModelResults.unregister();
// reducedModelResults.unregister(); // global.pFullPatternSimulationMesh->unregister();
// global.fullPatternResults[simulationScenarioIndex].unregister(); // global.fullPatternResults[simulationScenarioIndex].unregister();
#endif //#endif
// if (global.optimizationSettings.normalizationStrategy != // if (global.optimizationSettings.normalizationStrategy !=
// NormalizationStrategy::Epsilon && // NormalizationStrategy::Epsilon &&
// simulationScenarioError > 1) { // simulationScenarioError > 1) {
@ -231,35 +243,37 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
// global.reducedToFullInterfaceViMap, true); // global.reducedToFullInterfaceViMap, true);
// polyscope::removeAllStructures(); // polyscope::removeAllStructures();
// #endif // POLYSCOPE_DEFINED // #endif // POLYSCOPE_DEFINED
totalError += simulationScenarioError; totalError += simulationScenarioError * simulationScenarioError;
} }
} }
// std::cout << error << std::endl; // std::cout << error << std::endl;
if (totalError < global.minY) { if (totalError < global.minY) {
global.minY = totalError; global.minY = totalError;
global.minX.assign(x, x + n); global.minX.assign(x.begin(), x.begin() + n);
} }
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
++global.numberOfFunctionCalls; ++global.numberOfFunctionCalls;
global.objectiveValueHistory.push_back(
std::sqrt(totalError / global.simulationScenarioIndices.size()));
if (global.optimizationSettings.numberOfFunctionCalls >= 100 if (global.optimizationSettings.numberOfFunctionCalls >= 100
&& global.numberOfFunctionCalls % (global.optimizationSettings.numberOfFunctionCalls / 100) && global.numberOfFunctionCalls % (global.optimizationSettings.numberOfFunctionCalls / 100)
== 0) { == 0) {
std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl; std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl;
auto xPlot = matplot::linspace(0,
global.objectiveValueHistory.size(),
global.objectiveValueHistory.size());
// std::vector<double> colors(global.gObjectiveValueHistory.size(), 2);
// if (global.g_firstRoundIterationIndex != 0) {
// for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(),
// [](double &c) { c = 0.7; });
// }
global.gPlotHandle = matplot::scatter(xPlot, global.objectiveValueHistory);
SimulationResultsReporter::createPlot("Number of Steps",
"Objective value",
global.objectiveValueHistory);
} }
#endif #endif
// compute error and return it // compute error and return it
// global.objectiveValueHistory.push_back(totalError);
// auto xPlot = matplot::linspace(0, global.objectiveValueHistory.size(),
// global.objectiveValueHistory.size());
// std::vector<double> colors(global.gObjectiveValueHistory.size(), 2);
// if (global.g_firstRoundIterationIndex != 0) {
// for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(),
// [](double &c) { c = 0.7; });
// }
// global.gPlotHandle = matplot::scatter(xPlot, global.objectiveValueHistory);
// SimulationResultsReporter::createPlot("Number of Steps", "Objective value",
// global.objectiveValueHistory);
return totalError; return totalError;
} }
@ -271,18 +285,19 @@ void ReducedModelOptimizer::createSimulationMeshes(
std::cerr << "Error: A rectangular cross section is expected." << std::endl; std::cerr << "Error: A rectangular cross section is expected." << std::endl;
terminate(); terminate();
} }
const double ym = 1 * 1e9;
// Full pattern // Full pattern
pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel); pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
pFullPatternSimulationMesh->setBeamCrossSection( pFullPatternSimulationMesh->setBeamCrossSection(
CrossSectionType{0.002, 0.002}); CrossSectionType{0.002, 0.002});
pFullPatternSimulationMesh->setBeamMaterial(0.3, 2.3465 * 1e9); pFullPatternSimulationMesh->setBeamMaterial(0.3, ym);
// Reduced pattern // Reduced pattern
pReducedPatternSimulationMesh = pReducedPatternSimulationMesh =
std::make_shared<SimulationMesh>(reducedModel); std::make_shared<SimulationMesh>(reducedModel);
pReducedPatternSimulationMesh->setBeamCrossSection( pReducedPatternSimulationMesh->setBeamCrossSection(
CrossSectionType{0.002, 0.002}); CrossSectionType{0.002, 0.002});
pReducedPatternSimulationMesh->setBeamMaterial(0.3, 2.3465 * 1e9); pReducedPatternSimulationMesh->setBeamMaterial(0.3, ym);
} }
void ReducedModelOptimizer::createSimulationMeshes( void ReducedModelOptimizer::createSimulationMeshes(
@ -480,79 +495,141 @@ void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern,
initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters); initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters);
} }
void updateMesh(long n, const double *x) {
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
->pMesh;
const double E=global.initialParameters(0)*x[0];
const double A=global.initialParameters(1) * x[1];
const double beamWidth=std::sqrt(A);
const double beamHeight=beamWidth;
const double J=global.initialParameters(2) * x[2];
const double I2 = global.initialParameters(3) * x[3];
const double I3=global.initialParameters(4) * x[4];
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.setDimensions(
RectangularBeamDimensions(beamWidth,
beamHeight));
e.setMaterial(ElementMaterial(e.material.poissonsRatio,
E));
e.J = J;
e.I2 = I2;
e.I3 = I3;
}
assert(pReducedPatternSimulationMesh->EN() == 12);
assert(n >= 2);
// CoordType center_barycentric(1, 0, 0);
// CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5);
// CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric)
// * x[n - 2]);
CoordType movableVertex_barycentric(1 - x[n - 2], x[n - 2] / 2, x[n - 2] / 2);
CoordType baseTriangleMovableVertexPosition = global.baseTriangle.cP(0)
* movableVertex_barycentric[0]
+ global.baseTriangle.cP(1)
* movableVertex_barycentric[1]
+ global.baseTriangle.cP(2)
* movableVertex_barycentric[2];
baseTriangleMovableVertexPosition
= vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, vcg::math::ToRad(x[n - 1]))
* baseTriangleMovableVertexPosition;
for (int rotationCounter = 0;
rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) {
pReducedPatternSimulationMesh->vert[2 * rotationCounter].P()
= vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
vcg::math::ToRad(60.0 * rotationCounter))
* baseTriangleMovableVertexPosition;
}
pReducedPatternSimulationMesh->reset();
//#ifdef POLYSCOPE_DEFINED
// pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
// pReducedPatternSimulationMesh->registerForDrawing();
// std::cout << "Angle:" + std::to_string(x[n - 1]) + " size:" + std::to_string(x[n - 2])
// << std::endl;
// std::cout << "Verts:" << pReducedPatternSimulationMesh->VN() << std::endl;
// polyscope::show();
//#endif
}
void ReducedModelOptimizer::initializeOptimizationParameters( void ReducedModelOptimizer::initializeOptimizationParameters(
const std::shared_ptr<SimulationMesh> &mesh,const int& optimizationParamters) { const std::shared_ptr<SimulationMesh> &mesh, const int &optimizationParamters)
{
global.numberOfOptimizationParameters = optimizationParamters; global.numberOfOptimizationParameters = optimizationParamters;
global.initialParameters.resize(global.numberOfOptimizationParameters); global.initialParameters.resize(global.numberOfOptimizationParameters);
global.initialParameters(0) = mesh->elements[0].material.youngsModulus; switch (optimizationParamters) {
global.initialParameters(1) = mesh->elements[0].A; case 2:
global.initialParameters(2) = mesh->elements[0].J; function_updateReducedPattern_material =
global.initialParameters(3) = mesh->elements[0].I2; [](const dlib::matrix<double, 0, 1> &x,
global.initialParameters(4) = mesh->elements[0].I3; std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {};
break;
case 3:
global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
function_updateReducedPattern_material =
[&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
const double E = global.initialParameters(0) * x(0);
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
}
};
break;
case 4:
global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
global.initialParameters(1) = mesh->elements[0].A;
function_updateReducedPattern_material =
[&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
const double E = global.initialParameters(0) * x(0);
const double A = global.initialParameters(1) * x(1);
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight));
e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
}
};
break;
case 5:
global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
global.initialParameters(1) = mesh->elements[0].A;
global.initialParameters(2) = mesh->elements[0].inertia.I2;
function_updateReducedPattern_material =
[&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
const double E = global.initialParameters(0) * x(0);
const double A = global.initialParameters(1) * x(1);
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
const double I = global.initialParameters(2) * x(2);
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight));
e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
e.inertia.I2 = I;
e.inertia.I3 = I;
}
};
break;
case 7:
global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
global.initialParameters(1) = mesh->elements[0].A;
global.initialParameters(2) = mesh->elements[0].inertia.J;
global.initialParameters(3) = mesh->elements[0].inertia.I2;
global.initialParameters(4) = mesh->elements[0].inertia.I3;
function_updateReducedPattern_material = [&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh>
&pReducedPatternSimulationMesh) {
const double E = global.initialParameters(0) * x(0);
const double A = global.initialParameters(1) * x(1);
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
const double J = global.initialParameters(2) * x(2);
const double I2 = global.initialParameters(3) * x(3);
const double I3 = global.initialParameters(4) * x(4);
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight));
e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
e.inertia.J = J;
e.inertia.I2 = I2;
e.inertia.I3 = I3;
}
//#ifdef POLYSCOPE_DEFINED
// pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
// pReducedPatternSimulationMesh->registerForDrawing();
// std::cout << "Angle:" + std::to_string(x[n - 1]) + " size:" + std::to_string(x[n - 2])
// << std::endl;
// std::cout << "Verts:" << pReducedPatternSimulationMesh->VN() << std::endl;
// polyscope::show();
//#endif
};
break;
default:
std::cerr << "Wrong number of parameters:" << optimizationParamters << std::endl;
std::terminate();
break;
}
function_updateReducedPattern_geometry = [&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh>
&pReducedPatternSimulationMesh) {
const int n = x.size();
assert(n >= 2);
// CoordType center_barycentric(1, 0, 0);
// CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5);
// CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric)
// * x[n - 2]);
CoordType movableVertex_barycentric(1 - x(n - 2), x(n - 2) / 2, x(n - 2) / 2);
CoordType baseTriangleMovableVertexPosition = global.baseTriangle.cP(0)
* movableVertex_barycentric[0]
+ global.baseTriangle.cP(1)
* movableVertex_barycentric[1]
+ global.baseTriangle.cP(2)
* movableVertex_barycentric[2];
baseTriangleMovableVertexPosition
= vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
vcg::math::ToRad(x(n - 1)))
* baseTriangleMovableVertexPosition;
for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize;
rotationCounter++) {
pReducedPatternSimulationMesh->vert[2 * rotationCounter].P()
= vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
vcg::math::ToRad(60.0 * rotationCounter))
* baseTriangleMovableVertexPosition;
}
};
} }
void ReducedModelOptimizer::computeReducedModelSimulationJob( void ReducedModelOptimizer::computeReducedModelSimulationJob(
@ -715,7 +792,7 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza
results.optimalXNameValuePairs.resize(settings.xRanges.size()); results.optimalXNameValuePairs.resize(settings.xRanges.size());
std::vector<double> optimalX(settings.xRanges.size()); std::vector<double> optimalX(settings.xRanges.size());
for (int xVariableIndex = 0; xVariableIndex < settings.xRanges.size(); xVariableIndex++) { for (int xVariableIndex = 0; xVariableIndex < settings.xRanges.size(); xVariableIndex++) {
if (xVariableIndex < 5) { if (xVariableIndex < settings.xRanges.size() - 2) {
results.optimalXNameValuePairs[xVariableIndex] results.optimalXNameValuePairs[xVariableIndex]
= std::make_pair(settings.xRanges[xVariableIndex].label, = std::make_pair(settings.xRanges[xVariableIndex].label,
global.minX[xVariableIndex] global.minX[xVariableIndex]
@ -739,7 +816,10 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza
#endif #endif
// Compute obj value per simulation scenario and the raw objective value // Compute obj value per simulation scenario and the raw objective value
updateMesh(optimalX.size(), optimalX.data()); // updateMeshFunction(optimalX);
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh
= global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh;
function_updateReducedPattern(optimizationResult_dlib.x, pReducedPatternSimulationMesh);
// results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios); // results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios);
//TODO:use push_back it will make the code more readable //TODO:use push_back it will make the code more readable
LinearSimulationModel simulator; LinearSimulationModel simulator;
@ -754,6 +834,10 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza
global.simulationScenarioIndices.size()); global.simulationScenarioIndices.size());
results.objectiveValue.perSimulationScenario_total.resize( results.objectiveValue.perSimulationScenario_total.resize(
global.simulationScenarioIndices.size()); global.simulationScenarioIndices.size());
#ifdef POLYSCOPE_DEFINED
global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed);
#endif
results.perScenario_fullPatternPotentialEnergy.resize(global.simulationScenarioIndices.size()); results.perScenario_fullPatternPotentialEnergy.resize(global.simulationScenarioIndices.size());
for (int i = 0; i < global.simulationScenarioIndices.size(); i++) { for (int i = 0; i < global.simulationScenarioIndices.size(); i++) {
const int simulationScenarioIndex = global.simulationScenarioIndices[i]; const int simulationScenarioIndex = global.simulationScenarioIndices[i];
@ -828,17 +912,14 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza
std::cout << "Total Error value:" << results.objectiveValue.perSimulationScenario_total[i] std::cout << "Total Error value:" << results.objectiveValue.perSimulationScenario_total[i]
<< std::endl; << std::endl;
std::cout << std::endl; std::cout << std::endl;
#endif
}
const bool printDebugInfo = false; // reducedModelResults.registerForDrawing(Colors::reducedDeformed);
if (printDebugInfo) { // global.fullPatternResults[simulationScenarioIndex].registerForDrawing(Colors::fullDeformed);
std::cout << "Finished optimizing." << endl; // polyscope::show();
std::cout << "Total optimal objective value:" << results.objectiveValue.total << std::endl; // reducedModelResults.unregister();
assert(global.minY == optimizationResult_dlib.y); // global.fullPatternResults[simulationScenarioIndex].unregister();
if (global.minY != optimizationResult_dlib.y) {
std::cerr << "ERROR in objective value" << std::endl; #endif
}
} }
for (int simulationScenarioIndex : global.simulationScenarioIndices) { for (int simulationScenarioIndex : global.simulationScenarioIndices) {
@ -852,6 +933,7 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza
// global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->registerForDrawing(); // global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->registerForDrawing();
// global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel(temp); // global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel(temp);
} }
// results.draw();
} }
std::vector<std::pair<BaseSimulationScenario, double>> std::vector<std::pair<BaseSimulationScenario, double>>
@ -890,10 +972,15 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces(
fullPatternSimulationScenarioMaxMagnitudes fullPatternSimulationScenarioMaxMagnitudes
= static_cast<std::vector<std::pair<BaseSimulationScenario, double>>>( = static_cast<std::vector<std::pair<BaseSimulationScenario, double>>>(
json.at("maxMagn")); json.at("maxMagn"));
} else { const bool shouldRecompute = fullPatternSimulationScenarioMaxMagnitudes.size()
!= desiredBaseSimulationScenarioIndices.size();
if (!shouldRecompute) {
return fullPatternSimulationScenarioMaxMagnitudes;
}
}
#endif #endif
fullPatternSimulationScenarioMaxMagnitudes = computeFullPatternMaxSimulationForces( fullPatternSimulationScenarioMaxMagnitudes = computeFullPatternMaxSimulationForces(
desiredBaseSimulationScenarioIndices); desiredBaseSimulationScenarioIndices);
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
nlohmann::json json; nlohmann::json json;
@ -902,34 +989,110 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces(
std::filesystem::create_directories(forceMagnitudesDirectoryPath); std::filesystem::create_directories(forceMagnitudesDirectoryPath);
std::ofstream jsonFile(patternMaxForceMagnitudesFilePath.string()); std::ofstream jsonFile(patternMaxForceMagnitudesFilePath.string());
jsonFile << json; jsonFile << json;
}
#endif
assert(fullPatternSimulationScenarioMaxMagnitudes.size()
== desiredBaseSimulationScenarioIndices.size());
return fullPatternSimulationScenarioMaxMagnitudes; #endif
assert(fullPatternSimulationScenarioMaxMagnitudes.size()
== desiredBaseSimulationScenarioIndices.size());
return fullPatternSimulationScenarioMaxMagnitudes;
} }
void ReducedModelOptimizer::runOptimization(const Settings &settings, void ReducedModelOptimizer::runOptimization(const Settings &settings,
ReducedPatternOptimization::Results &results) ReducedPatternOptimization::Results &results)
{ {
global.objectiveValueHistory.clear(); global.objectiveValueHistory.clear();
dlib::matrix<double, 0, 1> xMin(global.numberOfOptimizationParameters); double (*objF)(const dlib::matrix<double, 0, 1> &) = &objective;
dlib::matrix<double, 0, 1> xMax(global.numberOfOptimizationParameters);
for (int i = 0; i < global.numberOfOptimizationParameters; i++) { //Geometry optimization of the reduced pattern
xMin(i) = settings.xRanges[i].min; constexpr int numberOfGeometryOptimizationParameters = 2;
xMax(i) = settings.xRanges[i].max; dlib::matrix<double, 0, 1> xGeometryMin(numberOfGeometryOptimizationParameters);
dlib::matrix<double, 0, 1> xGeometryMax(numberOfGeometryOptimizationParameters);
const int geometryParametersOffset
= global.numberOfOptimizationParameters
- numberOfGeometryOptimizationParameters; //I assume the geometry parameters come at the end
int paramIndex = 0;
for (int i = geometryParametersOffset; i < global.numberOfOptimizationParameters;
i++, paramIndex++) {
xGeometryMin(paramIndex) = settings.xRanges[i].min;
xGeometryMax(paramIndex) = settings.xRanges[i].max;
} }
dlib::function_evaluation result_dlib; //Set reduced pattern update functions to be used during optimization
double (*objF)(double, double, double, double, double,double,double) = &objective; function_updateReducedPattern =
result_dlib = dlib::find_min_global(objF, [&](const dlib::matrix<double, 0, 1> &x,
xMin, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
xMax, function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
dlib::max_function_calls(settings.numberOfFunctionCalls), pReducedPatternSimulationMesh->reset();
std::chrono::hours(24 * 365 * 290), };
settings.solverAccuracy);
getResults(result_dlib, settings, results); dlib::function_evaluation result_dlib_geometry
= dlib::find_min_global(objF,
xGeometryMin,
xGeometryMax,
dlib::max_function_calls(settings.numberOfFunctionCalls),
std::chrono::hours(24 * 365 * 290),
settings.solverAccuracy);
//Material optimization of the reduced pattern
const int numberOfMaterialOptimizationParameters = global.numberOfOptimizationParameters
- numberOfGeometryOptimizationParameters;
std::cout << "opt size:" << result_dlib_geometry.x(0) << std::endl;
std::cout << "opt size:" << global.minX[0] << std::endl;
dlib::function_evaluation result_dlib_material;
if (numberOfMaterialOptimizationParameters != 0) {
dlib::matrix<double, 0, 1> xMaterialMin(numberOfMaterialOptimizationParameters);
dlib::matrix<double, 0, 1> xMaterialMax(numberOfMaterialOptimizationParameters);
for (int i = 0; i < numberOfMaterialOptimizationParameters; i++) {
xMaterialMin(i) = settings.xRanges[i].min;
xMaterialMax(i) = settings.xRanges[i].max;
}
function_updateReducedPattern =
[&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
pReducedPatternSimulationMesh->reset();
};
result_dlib_material = dlib::find_min_global(objF,
xMaterialMin,
xMaterialMax,
dlib::max_function_calls(
settings.numberOfFunctionCalls),
std::chrono::hours(24 * 365 * 290),
settings.solverAccuracy);
}
// constexpr bool useBOBYQA = false;
// if (useBOBYQA) {
// const size_t npt = 2 * global.numberOfOptimizationParameters;
// // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
// const double rhobeg = 0.1;
// // const double rhobeg = 10;
// const double rhoend = rhobeg * 1e-6;
// // const size_t maxFun = 10 * (x.size() ^ 2);
// const size_t bobyqa_maxFunctionCalls = 200;
// dlib::find_min_bobyqa(objF,
// result_dlib.x,
// npt,
// xMaterialMin,
// xMaterialMax,
// rhobeg,
// rhoend,
// bobyqa_maxFunctionCalls);
// }
dlib::function_evaluation optimalResult;
optimalResult.x.set_size(global.numberOfOptimizationParameters);
std::copy(result_dlib_material.x.begin(), result_dlib_material.x.end(), optimalResult.x.begin());
std::copy(result_dlib_geometry.x.begin(),
result_dlib_geometry.x.end(),
optimalResult.x.begin() + geometryParametersOffset);
function_updateReducedPattern =
[&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
pReducedPatternSimulationMesh->reset();
};
optimalResult.y = objective(optimalResult.x);
getResults(optimalResult, settings, results);
} }
void ReducedModelOptimizer::constructAxialSimulationScenario( void ReducedModelOptimizer::constructAxialSimulationScenario(
@ -1186,7 +1349,6 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce(
{ {
double forceMagnitude = 1; double forceMagnitude = 1;
double minimumError; double minimumError;
double translationalOptimizationEpsilon;
bool wasSuccessful = false; bool wasSuccessful = false;
global.constructScenarioFunction = constructBaseScenarioFunctions[scenario]; global.constructScenarioFunction = constructBaseScenarioFunctions[scenario];
const double baseTriangleHeight = vcg::Distance(global.baseTriangle.cP(0), const double baseTriangleHeight = vcg::Distance(global.baseTriangle.cP(0),
@ -1194,6 +1356,7 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce(
+ global.baseTriangle.cP(2)) + global.baseTriangle.cP(2))
/ 2); / 2);
std::function<double(const double &)> objectiveFunction; std::function<double(const double &)> objectiveFunction;
double translationalOptimizationEpsilon{baseTriangleHeight * 0.001};
double objectiveEpsilon = translationalOptimizationEpsilon; double objectiveEpsilon = translationalOptimizationEpsilon;
objectiveFunction = &fullPatternMaxSimulationForceTranslationalObjective; objectiveFunction = &fullPatternMaxSimulationForceTranslationalObjective;
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first; global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
@ -1359,12 +1522,13 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
for (int simulationScenarioIndex : global.simulationScenarioIndices) { for (int simulationScenarioIndex : global.simulationScenarioIndices) {
if (global.optimizationSettings.normalizationStrategy == if (global.optimizationSettings.normalizationStrategy ==
Settings::NormalizationStrategy::Epsilon) { Settings::NormalizationStrategy::Epsilon) {
const double epsilon_translationalDisplacement = global.optimizationSettings const double epsilon_translationalDisplacement
.normalizationParameter; = global.optimizationSettings.translationNormalizationParameter;
global.translationalDisplacementNormalizationValues[simulationScenarioIndex] global.translationalDisplacementNormalizationValues[simulationScenarioIndex]
= std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex], = std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex],
epsilon_translationalDisplacement); epsilon_translationalDisplacement);
const double epsilon_rotationalDisplacement = vcg::math::ToRad(3.0); const double epsilon_rotationalDisplacement = global.optimizationSettings
.rotationNormalizationParameter;
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
= std::max(fullPatternAngularDistance[simulationScenarioIndex], = std::max(fullPatternAngularDistance[simulationScenarioIndex],
epsilon_rotationalDisplacement); epsilon_rotationalDisplacement);
@ -1415,10 +1579,11 @@ void ReducedModelOptimizer::optimize(
results.baseTriangle = global.baseTriangle; results.baseTriangle = global.baseTriangle;
DRMSimulationModel::Settings simulationSettings; DRMSimulationModel::Settings simulationSettings;
simulationSettings.maxDRMIterations = 200000; // simulationSettings.maxDRMIterations = 200000;
simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8; // simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8;
simulationSettings.viscousDampingFactor = 5e-3; // simulationSettings.viscousDampingFactor = 5e-3;
simulationSettings.useKineticDamping = true; // simulationSettings.useKineticDamping = true;
// simulationSettings.averageResidualForcesCriterionThreshold = 1e-5; // simulationSettings.averageResidualForcesCriterionThreshold = 1e-5;
// simulationSettings.viscousDampingFactor = 1e-3; // simulationSettings.viscousDampingFactor = 1e-3;
// simulationSettings.beVerbose = true; // simulationSettings.beVerbose = true;
@ -1426,7 +1591,7 @@ void ReducedModelOptimizer::optimize(
// simulationSettings.isDebugMode = true; // simulationSettings.isDebugMode = true;
// simulationSettings.debugModeStep = 100000; // simulationSettings.debugModeStep = 100000;
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
const bool drawFullPatternSimulationResults = false; constexpr bool drawFullPatternSimulationResults = false;
if (drawFullPatternSimulationResults) { if (drawFullPatternSimulationResults) {
global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
ReducedPatternOptimization::Colors::fullInitial); ReducedPatternOptimization::Colors::fullInitial);

21
src/reducedmodeloptimizer.hpp
View File

@ -177,6 +177,18 @@ public:
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>> const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
&oppositeInterfaceViPairs, &oppositeInterfaceViPairs,
SimulationJob &job); SimulationJob &job);
static std::function<void(const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &m)>
function_updateReducedPattern;
static std::function<void(const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &m)>
function_updateReducedPattern_material;
static std::function<void(const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &m)>
function_updateReducedPattern_material_E;
static std::function<void(const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &m)>
function_updateReducedPattern_geometry;
private: private:
static void computeDesiredReducedModelDisplacements( static void computeDesiredReducedModelDisplacements(
@ -211,6 +223,13 @@ private:
getFullPatternMaxSimulationForces( getFullPatternMaxSimulationForces(
const std::vector<ReducedPatternOptimization::BaseSimulationScenario> const std::vector<ReducedPatternOptimization::BaseSimulationScenario>
&desiredBaseSimulationScenarioIndices); &desiredBaseSimulationScenarioIndices);
static double objective(const dlib::matrix<double, 0, 1> &x);
}; };
void updateMesh(long n, const double *x); inline std::function<void(const dlib::matrix<double, 0, 1> &x, std::shared_ptr<SimulationMesh> &m)>
ReducedModelOptimizer::function_updateReducedPattern;
inline std::function<void(const dlib::matrix<double, 0, 1> &x, std::shared_ptr<SimulationMesh> &m)>
ReducedModelOptimizer::function_updateReducedPattern_material;
inline std::function<void(const dlib::matrix<double, 0, 1> &x, std::shared_ptr<SimulationMesh> &m)>
ReducedModelOptimizer::function_updateReducedPattern_geometry;
#endif // REDUCEDMODELOPTIMIZER_HPP #endif // REDUCEDMODELOPTIMIZER_HPP