From b763ca92e7e7d523955dcfeaaa028d150ea56d3a Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Mon, 8 Nov 2021 10:47:35 +0200 Subject: [PATCH] Splitted the update of the reduced pattern modular by splitting the geometry from the material parameters. --- src/main.cpp | 45 +++- src/reducedmodeloptimizer.cpp | 471 +++++++++++++++++++++++----------- src/reducedmodeloptimizer.hpp | 21 +- 3 files changed, 369 insertions(+), 168 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 5ae434d..d9b6dd2 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -38,7 +38,7 @@ int main(int argc, char *argv[]) { reducedPattern.scale(0.03, interfaceNodeIndex); // 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 beamI2{"I2", 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 innerHexagonAngle{"HexAngle", -30.0, 30.0}; ReducedPatternOptimization::Settings settings_optimization; - settings_optimization.xRanges = {beamE,beamA,beamJ,beamI2,beamI3, - innerHexagonSize, innerHexagonAngle}; + // settings_optimization.xRanges + // = {beamE, beamA, beamJ, beamI2, beamI3, innerHexagonSize, innerHexagonAngle}; + settings_optimization.xRanges = {beamE, beamA, beamI2, innerHexagonSize, innerHexagonAngle}; const bool input_numberOfFunctionCallsDefined = argc >= 4; settings_optimization.numberOfFunctionCalls = input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100; settings_optimization.normalizationStrategy = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon; - settings_optimization.normalizationParameter = 0.0003; - settings_optimization.solverAccuracy = 0.001; + settings_optimization.translationNormalizationParameter = 1e-3; + 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.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 optimizationName = pairName + "(" + std::to_string(settings_optimization.numberOfFunctionCalls) + "_" + to_string_with_precision( settings_optimization.objectiveWeights.translational) - + ")"; + + ")" + "_" + xConcatNames; const bool input_resultDirectoryDefined = argc >= 6; const std::string optimizationResultsDirectory = input_resultDirectoryDefined ? argv[5] @@ -91,13 +102,14 @@ int main(int argc, char *argv[]) { } ReducedPatternOptimization::Results optimizationResults; - bool optimizationAlreadyComputed = optimizationResultFolderExists; - // bool optimizationAlreadyComputed = false; - // if (optimizationResultFolderExists) { - // const bool resultsWereSuccessfullyLoaded = optimizationResults.load(resultsOutputDir); - // if (resultsWereSuccessfullyLoaded && optimizationResults.settings == settings_optimization) { - // } - // } + constexpr bool shouldReoptimize = true; + bool optimizationAlreadyComputed = false; + if (!shouldReoptimize && optimizationResultFolderExists) { + const bool resultsWereSuccessfullyLoaded = optimizationResults.load(resultsOutputDir); + if (resultsWereSuccessfullyLoaded && optimizationResults.settings == settings_optimization) { + optimizationAlreadyComputed = true; + } + } if (!optimizationAlreadyComputed) { auto start = std::chrono::system_clock::now(); @@ -135,6 +147,11 @@ int main(int argc, char *argv[]) { #ifdef POLYSCOPE_DEFINED // optimizationResults.saveMeshFiles(); + std::cout << "E:" + << optimizationResults.reducedPatternSimulationJobs[0] + ->pMesh->elements[0] + .material.youngsModulus + << std::endl; optimizationResults.draw(); #endif diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 8f10819..d08078b 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -52,6 +52,8 @@ double ReducedModelOptimizer::computeDisplacementError( const double rawError = computeRawTranslationalError(fullPatternDisplacements, reducedPatternDisplacements, reducedToFullInterfaceViMap); + // std::cout << "raw trans error:" << rawError << std::endl; + // std::cout << "raw trans error:" << normalizationFactor << std::endl; return rawError / normalizationFactor; } @@ -120,6 +122,9 @@ double ReducedModelOptimizer::computeError( simulationResults_reducedPattern.displacements, reducedToFullInterfaceViMap, normalizationFactor_translationalDisplacement); + + // std::cout << "normalization factor:" << normalizationFactor_rotationalDisplacement << std::endl; + // std::cout << "trans error:" << translationalError << std::endl; const double rotationalError = computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, simulationResults_reducedPattern.rotationalDisplacementQuaternion, @@ -129,18 +134,19 @@ double ReducedModelOptimizer::computeError( + global.optimizationSettings.objectiveWeights.rotational * rotationalError; } -double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3, - double innerHexagonSize, - double innerHexagonRotationAngle) { - std::vector x{E,A,J,I2,I3, innerHexagonSize, innerHexagonRotationAngle}; - return ReducedModelOptimizer::objective(x.size(), x.data()); -} +//double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3, +// double innerHexagonSize, +// double innerHexagonRotationAngle) { +// std::vector x{E,A,J,I2,I3, innerHexagonSize, innerHexagonRotationAngle}; +// return ReducedModelOptimizer::objective(x.size(), x.data()); +//} -double ReducedModelOptimizer::objective(long n, const double *x) { - // std::cout.precision(17); - // for (int i = 0; i < n; i++) { - // std::cout << x[i] << " "; - // } +double ReducedModelOptimizer::objective(const dlib::matrix &x) +{ + // std::cout.precision(17); + // for (int i = 0; i < x.size(); i++) { + // std::cout << x(i) << " "; + // } // std::cout << 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.secondBendingConstFactor // << std::endl; - updateMesh(n, x); + const int n = x.size(); + std::shared_ptr &pReducedPatternSimulationMesh + = global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh; + function_updateReducedPattern(x, pReducedPatternSimulationMesh); // global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(); // global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(); // polyscope::show(); @@ -168,6 +177,9 @@ double ReducedModelOptimizer::objective(long n, const double *x) { for (const int simulationScenarioIndex : global.simulationScenarioIndices) { const std::shared_ptr &reducedJob = global.reducedPatternSimulationJobs[simulationScenarioIndex]; + //#ifdef POLYSCOPE_DEFINED + // std::cout << reducedJob->getLabel() << ":" << std::endl; + //#endif SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob); // std::string filename; if (!reducedModelResults.converged) { @@ -177,7 +189,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) { std::cout << "Failed simulation" << std::endl; #endif } else { - const bool usePotentialEnergy = false; + constexpr bool usePotentialEnergy = false; double simulationScenarioError; if (usePotentialEnergy) { simulationScenarioError = std::abs( @@ -191,17 +203,17 @@ double ReducedModelOptimizer::objective(long n, const double *x) { global.translationalDisplacementNormalizationValues[simulationScenarioIndex], global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); } -#ifdef POLYSCOPE_DEFINED -// std::cout << reducedJob->getLabel() << " sim error:" << simulationScenarioError -// << std::endl; -// reducedJob->pMesh->registerForDrawing(Colors::reducedInitial); -// reducedModelResults.registerForDrawing(Colors::reducedDeformed); -// global.fullPatternResults[simulationScenarioIndex].registerForDrawing( -// Colors::fullDeformed); -// polyscope::show(); -// reducedModelResults.unregister(); -// global.fullPatternResults[simulationScenarioIndex].unregister(); -#endif + //#ifdef POLYSCOPE_DEFINED + // reducedJob->pMesh->registerForDrawing(Colors::reducedInitial); + // reducedModelResults.registerForDrawing(Colors::reducedDeformed); + // global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); + // global.fullPatternResults[simulationScenarioIndex].registerForDrawing( + // Colors::fullDeformed); + // polyscope::show(); + // reducedModelResults.unregister(); + // global.pFullPatternSimulationMesh->unregister(); + // global.fullPatternResults[simulationScenarioIndex].unregister(); + //#endif // if (global.optimizationSettings.normalizationStrategy != // NormalizationStrategy::Epsilon && // simulationScenarioError > 1) { @@ -231,35 +243,37 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // global.reducedToFullInterfaceViMap, true); // polyscope::removeAllStructures(); // #endif // POLYSCOPE_DEFINED - totalError += simulationScenarioError; + totalError += simulationScenarioError * simulationScenarioError; } } // std::cout << error << std::endl; if (totalError < global.minY) { global.minY = totalError; - global.minX.assign(x, x + n); + global.minX.assign(x.begin(), x.begin() + n); } #ifdef POLYSCOPE_DEFINED ++global.numberOfFunctionCalls; + global.objectiveValueHistory.push_back( + std::sqrt(totalError / global.simulationScenarioIndices.size())); if (global.optimizationSettings.numberOfFunctionCalls >= 100 && global.numberOfFunctionCalls % (global.optimizationSettings.numberOfFunctionCalls / 100) == 0) { std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl; + auto xPlot = matplot::linspace(0, + global.objectiveValueHistory.size(), + global.objectiveValueHistory.size()); + // std::vector 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 -// compute error and return it -// global.objectiveValueHistory.push_back(totalError); -// auto xPlot = matplot::linspace(0, global.objectiveValueHistory.size(), -// global.objectiveValueHistory.size()); -// std::vector 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); - + // compute error and return it return totalError; } @@ -271,18 +285,19 @@ void ReducedModelOptimizer::createSimulationMeshes( std::cerr << "Error: A rectangular cross section is expected." << std::endl; terminate(); } + const double ym = 1 * 1e9; // Full pattern pFullPatternSimulationMesh = std::make_shared(fullModel); pFullPatternSimulationMesh->setBeamCrossSection( CrossSectionType{0.002, 0.002}); - pFullPatternSimulationMesh->setBeamMaterial(0.3, 2.3465 * 1e9); + pFullPatternSimulationMesh->setBeamMaterial(0.3, ym); // Reduced pattern pReducedPatternSimulationMesh = std::make_shared(reducedModel); pReducedPatternSimulationMesh->setBeamCrossSection( CrossSectionType{0.002, 0.002}); - pReducedPatternSimulationMesh->setBeamMaterial(0.3, 2.3465 * 1e9); + pReducedPatternSimulationMesh->setBeamMaterial(0.3, ym); } void ReducedModelOptimizer::createSimulationMeshes( @@ -480,79 +495,141 @@ void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern, initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters); } -void updateMesh(long n, const double *x) { - std::shared_ptr &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( - const std::shared_ptr &mesh,const int& optimizationParamters) { + const std::shared_ptr &mesh, const int &optimizationParamters) +{ global.numberOfOptimizationParameters = optimizationParamters; global.initialParameters.resize(global.numberOfOptimizationParameters); - global.initialParameters(0) = mesh->elements[0].material.youngsModulus; - global.initialParameters(1) = mesh->elements[0].A; - global.initialParameters(2) = mesh->elements[0].J; - global.initialParameters(3) = mesh->elements[0].I2; - global.initialParameters(4) = mesh->elements[0].I3; + switch (optimizationParamters) { + case 2: + function_updateReducedPattern_material = + [](const dlib::matrix &x, + std::shared_ptr &pReducedPatternSimulationMesh) {}; + break; + + case 3: + global.initialParameters(0) = mesh->elements[0].material.youngsModulus; + function_updateReducedPattern_material = + [&](const dlib::matrix &x, + std::shared_ptr &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 &x, + std::shared_ptr &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 &x, + std::shared_ptr &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 &x, + std::shared_ptr + &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 &x, + std::shared_ptr + &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( @@ -715,7 +792,7 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza results.optimalXNameValuePairs.resize(settings.xRanges.size()); std::vector optimalX(settings.xRanges.size()); for (int xVariableIndex = 0; xVariableIndex < settings.xRanges.size(); xVariableIndex++) { - if (xVariableIndex < 5) { + if (xVariableIndex < settings.xRanges.size() - 2) { results.optimalXNameValuePairs[xVariableIndex] = std::make_pair(settings.xRanges[xVariableIndex].label, global.minX[xVariableIndex] @@ -739,7 +816,10 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza #endif // Compute obj value per simulation scenario and the raw objective value - updateMesh(optimalX.size(), optimalX.data()); + // updateMeshFunction(optimalX); + std::shared_ptr &pReducedPatternSimulationMesh + = global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh; + function_updateReducedPattern(optimizationResult_dlib.x, pReducedPatternSimulationMesh); // results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios); //TODO:use push_back it will make the code more readable LinearSimulationModel simulator; @@ -754,6 +834,10 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza global.simulationScenarioIndices.size()); results.objectiveValue.perSimulationScenario_total.resize( global.simulationScenarioIndices.size()); +#ifdef POLYSCOPE_DEFINED + global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); +#endif + results.perScenario_fullPatternPotentialEnergy.resize(global.simulationScenarioIndices.size()); for (int i = 0; i < global.simulationScenarioIndices.size(); 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::endl; std::cout << std::endl; -#endif - } - const bool printDebugInfo = false; - if (printDebugInfo) { - std::cout << "Finished optimizing." << endl; - std::cout << "Total optimal objective value:" << results.objectiveValue.total << std::endl; - assert(global.minY == optimizationResult_dlib.y); - if (global.minY != optimizationResult_dlib.y) { - std::cerr << "ERROR in objective value" << std::endl; - } + // reducedModelResults.registerForDrawing(Colors::reducedDeformed); + // global.fullPatternResults[simulationScenarioIndex].registerForDrawing(Colors::fullDeformed); + // polyscope::show(); + // reducedModelResults.unregister(); + // global.fullPatternResults[simulationScenarioIndex].unregister(); + +#endif } 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->setLabel(temp); } + // results.draw(); } std::vector> @@ -890,10 +972,15 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces( fullPatternSimulationScenarioMaxMagnitudes = static_cast>>( json.at("maxMagn")); - } else { + const bool shouldRecompute = fullPatternSimulationScenarioMaxMagnitudes.size() + != desiredBaseSimulationScenarioIndices.size(); + if (!shouldRecompute) { + return fullPatternSimulationScenarioMaxMagnitudes; + } + } #endif - fullPatternSimulationScenarioMaxMagnitudes = computeFullPatternMaxSimulationForces( - desiredBaseSimulationScenarioIndices); + fullPatternSimulationScenarioMaxMagnitudes = computeFullPatternMaxSimulationForces( + desiredBaseSimulationScenarioIndices); #ifdef POLYSCOPE_DEFINED nlohmann::json json; @@ -902,34 +989,110 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces( std::filesystem::create_directories(forceMagnitudesDirectoryPath); std::ofstream jsonFile(patternMaxForceMagnitudesFilePath.string()); jsonFile << json; - } -#endif - assert(fullPatternSimulationScenarioMaxMagnitudes.size() - == desiredBaseSimulationScenarioIndices.size()); - return fullPatternSimulationScenarioMaxMagnitudes; +#endif + assert(fullPatternSimulationScenarioMaxMagnitudes.size() + == desiredBaseSimulationScenarioIndices.size()); + + return fullPatternSimulationScenarioMaxMagnitudes; } void ReducedModelOptimizer::runOptimization(const Settings &settings, ReducedPatternOptimization::Results &results) { global.objectiveValueHistory.clear(); - dlib::matrix xMin(global.numberOfOptimizationParameters); - dlib::matrix xMax(global.numberOfOptimizationParameters); - for (int i = 0; i < global.numberOfOptimizationParameters; i++) { - xMin(i) = settings.xRanges[i].min; - xMax(i) = settings.xRanges[i].max; + double (*objF)(const dlib::matrix &) = &objective; + + //Geometry optimization of the reduced pattern + constexpr int numberOfGeometryOptimizationParameters = 2; + dlib::matrix xGeometryMin(numberOfGeometryOptimizationParameters); + dlib::matrix 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; - double (*objF)(double, double, double, double, double,double,double) = &objective; - result_dlib = dlib::find_min_global(objF, - xMin, - xMax, - dlib::max_function_calls(settings.numberOfFunctionCalls), - std::chrono::hours(24 * 365 * 290), - settings.solverAccuracy); - getResults(result_dlib, settings, results); + //Set reduced pattern update functions to be used during optimization + function_updateReducedPattern = + [&](const dlib::matrix &x, + std::shared_ptr &pReducedPatternSimulationMesh) { + function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); + pReducedPatternSimulationMesh->reset(); + }; + + 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 xMaterialMin(numberOfMaterialOptimizationParameters); + dlib::matrix 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 &x, + std::shared_ptr &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 &x, + std::shared_ptr &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( @@ -1186,7 +1349,6 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( { double forceMagnitude = 1; double minimumError; - double translationalOptimizationEpsilon; bool wasSuccessful = false; global.constructScenarioFunction = constructBaseScenarioFunctions[scenario]; const double baseTriangleHeight = vcg::Distance(global.baseTriangle.cP(0), @@ -1194,6 +1356,7 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( + global.baseTriangle.cP(2)) / 2); std::function objectiveFunction; + double translationalOptimizationEpsilon{baseTriangleHeight * 0.001}; double objectiveEpsilon = translationalOptimizationEpsilon; objectiveFunction = &fullPatternMaxSimulationForceTranslationalObjective; global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first; @@ -1359,12 +1522,13 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() { for (int simulationScenarioIndex : global.simulationScenarioIndices) { if (global.optimizationSettings.normalizationStrategy == Settings::NormalizationStrategy::Epsilon) { - const double epsilon_translationalDisplacement = global.optimizationSettings - .normalizationParameter; + const double epsilon_translationalDisplacement + = global.optimizationSettings.translationNormalizationParameter; global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex], epsilon_translationalDisplacement); - const double epsilon_rotationalDisplacement = vcg::math::ToRad(3.0); + const double epsilon_rotationalDisplacement = global.optimizationSettings + .rotationNormalizationParameter; global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = std::max(fullPatternAngularDistance[simulationScenarioIndex], epsilon_rotationalDisplacement); @@ -1415,10 +1579,11 @@ void ReducedModelOptimizer::optimize( results.baseTriangle = global.baseTriangle; DRMSimulationModel::Settings simulationSettings; - simulationSettings.maxDRMIterations = 200000; - simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8; - simulationSettings.viscousDampingFactor = 5e-3; - simulationSettings.useKineticDamping = true; +// simulationSettings.maxDRMIterations = 200000; +// simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8; +// simulationSettings.viscousDampingFactor = 5e-3; +// simulationSettings.useKineticDamping = true; + // simulationSettings.averageResidualForcesCriterionThreshold = 1e-5; // simulationSettings.viscousDampingFactor = 1e-3; // simulationSettings.beVerbose = true; @@ -1426,7 +1591,7 @@ void ReducedModelOptimizer::optimize( // simulationSettings.isDebugMode = true; // simulationSettings.debugModeStep = 100000; #ifdef POLYSCOPE_DEFINED - const bool drawFullPatternSimulationResults = false; + constexpr bool drawFullPatternSimulationResults = false; if (drawFullPatternSimulationResults) { global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( ReducedPatternOptimization::Colors::fullInitial); diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index a8280e3..b5d08fb 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -177,6 +177,18 @@ public: const std::vector> &oppositeInterfaceViPairs, SimulationJob &job); + static std::function &x, + std::shared_ptr &m)> + function_updateReducedPattern; + static std::function &x, + std::shared_ptr &m)> + function_updateReducedPattern_material; + static std::function &x, + std::shared_ptr &m)> + function_updateReducedPattern_material_E; + static std::function &x, + std::shared_ptr &m)> + function_updateReducedPattern_geometry; private: static void computeDesiredReducedModelDisplacements( @@ -211,6 +223,13 @@ private: getFullPatternMaxSimulationForces( const std::vector &desiredBaseSimulationScenarioIndices); + static double objective(const dlib::matrix &x); }; -void updateMesh(long n, const double *x); +inline std::function &x, std::shared_ptr &m)> + ReducedModelOptimizer::function_updateReducedPattern; +inline std::function &x, std::shared_ptr &m)> + ReducedModelOptimizer::function_updateReducedPattern_material; +inline std::function &x, std::shared_ptr &m)> + ReducedModelOptimizer::function_updateReducedPattern_geometry; + #endif // REDUCEDMODELOPTIMIZER_HPP