diff --git a/CMakeLists.txt b/CMakeLists.txt index 85d0601..ee6f18c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -37,7 +37,6 @@ download_project(PROJ MYSOURCES endif() - #Polyscope download_project(PROJ POLYSCOPE GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git @@ -89,4 +88,5 @@ target_include_directories(${PROJECT_NAME} ) link_directories(${MYSOURCES_SOURCE_DIR}/boost_graph/libs) -target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib) +target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen matplot dlib::dlib) +#target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib) diff --git a/src/main.cpp b/src/main.cpp index 07b512b..6cf66d0 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -2,9 +2,6 @@ #include "csvfile.hpp" #include "edgemesh.hpp" #include "flatpattern.hpp" -#include "polyscope/curve_network.h" -#include "polyscope/point_cloud.h" -#include "polyscope/polyscope.h" #include "reducedmodeloptimizer.hpp" #include "simulationhistoryplotter.hpp" #include "trianglepattterntopology.hpp" @@ -16,6 +13,11 @@ #include #include +#ifdef POLYSCOPE_DEFINED +#include "polyscope/curve_network.h" +#include "polyscope/point_cloud.h" +#include "polyscope/polyscope.h" +#endif int main(int argc, char *argv[]) { if (argc < 3) { std::cerr << "Specify at least the two pattern filepaths to be " @@ -46,10 +48,11 @@ int main(int argc, char *argv[]) { ReducedModelOptimizer::Settings settings_optimization; settings_optimization.xRanges = {beamWidth, beamDimensionsRatio, beamE, innerHexagonSize}; - const bool input_numberOfFunctionCallsDefined = argc >= 5; + const bool input_numberOfFunctionCallsDefined = argc >= 4; settings_optimization.numberOfFunctionCalls = - input_numberOfFunctionCallsDefined ? std::atoi(argv[4]) : 100; - settings_optimization.normalizeObjectiveValue = true; + input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100; + settings_optimization.shouldNormalizeObjectiveValue = true; + settings_optimization.solutionAccuracy = 0.1; // Optimize pair const std::string pairName = @@ -62,9 +65,9 @@ int main(int argc, char *argv[]) { optimizer.optimize(settings_optimization); // Export results - const bool input_resultDirectoryDefined = argc >= 6; + const bool input_resultDirectoryDefined = argc >= 5; std::string optimizationResultsDirectory = - input_resultDirectoryDefined ? argv[5] : "OptimizationResults"; + input_resultDirectoryDefined ? argv[4] : "OptimizationResults"; if (optimizationResults.numberOfSimulationCrashes != 0) { const auto crashedJobsDirPath = std::filesystem::path(optimizationResultsDirectory) diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index d8623e6..34d904c 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -9,7 +9,7 @@ struct GlobalOptimizationVariables { std::vector g_optimalReducedModelDisplacements; std::vector> fullPatternDisplacements; - std::vector fullPatternDisplacementNormSum; + std::vector reducedPatternMaximumDisplacementNormSum; std::vector> fullPatternSimulationJobs; std::vector> reducedPatternSimulationJobs; std::unordered_map @@ -35,26 +35,49 @@ struct GlobalOptimizationVariables { ReducedModelOptimizer::Settings optimizationSettings; } global; +std::vector reducedPatternMaximumDisplacementSimulationJobs; + double ReducedModelOptimizer::computeError( const std::vector &reducedPatternDisplacements, const std::vector &fullPatternDisplacements, - const double &interfaceDisplacementsNormSum, + const std::unordered_map + &reducedToFullInterfaceViMap, + const double &normalizationFactor) { + const double rawError = + computeRawError(reducedPatternDisplacements, fullPatternDisplacements, + reducedToFullInterfaceViMap); + if (global.optimizationSettings.shouldNormalizeObjectiveValue) { + assert(rawError <= normalizationFactor); + if (rawError > normalizationFactor) { + std::cout << "Normalized error is bigger than one." << std::endl; + } + return rawError / + normalizationFactor; // std::max(interfaceDisplacementsNormSum, + // 0.00003); + } + + return rawError; +} + +double ReducedModelOptimizer::computeRawError( + const std::vector &reducedPatternDisplacements, + const std::vector &fullPatternDisplacements, const std::unordered_map &reducedToFullInterfaceViMap) { double error = 0; for (const auto reducedFullViPair : reducedToFullInterfaceViMap) { - VertexIndex reducedModelVi = reducedFullViPair.first; // const auto pos = // g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP(); // std::cout << "Interface vi " << reducedModelVi << " is at position " // << pos[0] << " " << pos[1] << " " << pos[2] << std::endl; Eigen::Vector3d reducedVertexDisplacement( - reducedPatternDisplacements[reducedModelVi][0], - reducedPatternDisplacements[reducedModelVi][1], - reducedPatternDisplacements[reducedModelVi][2]); + reducedPatternDisplacements[reducedFullViPair.first][0], + reducedPatternDisplacements[reducedFullViPair.first][1], + reducedPatternDisplacements[reducedFullViPair.first][2]); if (!std::isfinite(reducedVertexDisplacement[0]) || !std::isfinite(reducedVertexDisplacement[1]) || !std::isfinite(reducedVertexDisplacement[2])) { + std::cout << "Displacements are not finite" << std::endl; std::terminate(); } Eigen::Vector3d fullVertexDisplacement( @@ -67,9 +90,6 @@ double ReducedModelOptimizer::computeError( error += errorVector.norm(); } - if (global.optimizationSettings.normalizeObjectiveValue) { - return error / std::max(interfaceDisplacementsNormSum, 0.00003); - } return error; } @@ -132,8 +152,8 @@ void updateMesh(long n, const double *x) { pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); } -double ReducedModelOptimizer::objective(double b, double h, double E) { - std::vector x{b, h, E}; +double ReducedModelOptimizer::objective(double b, double r, double E) { + std::vector x{b, r, E}; return ReducedModelOptimizer::objective(x.size(), x.data()); } @@ -146,12 +166,9 @@ double ReducedModelOptimizer::objective(double b, double h, double E, double ReducedModelOptimizer::objective(long n, const double *x) { // std::cout.precision(17); - // for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { - // std::cout << "x[" + std::to_string(parameterIndex) + "]=" - // << x[parameterIndex] << std::endl; - // } - const Element &e = global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; - // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " + // const Element &e = + // global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout << + // e.axialConstFactor << " " << e.torsionConstFactor << " " // << e.firstBendingConstFactor << " " << // e.secondBendingConstFactor // << std::endl; @@ -167,51 +184,79 @@ double ReducedModelOptimizer::objective(long n, const double *x) { FormFinder::Settings simulationSettings; // simulationSettings.shouldDraw = true; for (const int simulationScenarioIndex : global.simulationScenarioIndices) { + // std::cout << "Executing " + // << simulationScenarioStrings[simulationScenarioIndex] + // << std::endl; SimulationResults reducedModelResults = simulator.executeSimulation( global.reducedPatternSimulationJobs[simulationScenarioIndex], simulationSettings); - std::string filename; + // std::string filename; if (!reducedModelResults.converged) { - error += std::numeric_limits::max(); global.numOfSimulationCrashes++; -if(beVerbose){ - std::cout << "Failed simulation" << std::endl; - } }else { - error += computeError( - reducedModelResults.displacements, - global.fullPatternDisplacements[simulationScenarioIndex], - global.fullPatternDisplacementNormSum[simulationScenarioIndex], - global.reducedToFullInterfaceViMap); - } - std::ofstream out(filename, std::ios_base::app); - auto pMesh = - global - .reducedPatternSimulationJobs[global.simulationScenarioIndices[0]] - ->pMesh; + bool beVerbose = false; // TODO: use an extern or data member for that + if (beVerbose) { + std::cout << "Failed simulation" << std::endl; + } + } else { + const double thisSimulationScenarioError = + computeError(reducedModelResults.displacements, + global.fullPatternDisplacements[simulationScenarioIndex], + global.reducedToFullInterfaceViMap, + global.reducedPatternMaximumDisplacementNormSum + [simulationScenarioIndex]); + if (thisSimulationScenarioError > 1) { + std::cout << "Simulation scenario " + << simulationScenarioStrings[simulationScenarioIndex] + << " results in an error bigger than one." << std::endl; + for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { + std::cout << "x[" + std::to_string(parameterIndex) + "]=" + << x[parameterIndex] << std::endl; + } - for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { - out << "x[" + std::to_string(parameterIndex) + "]=" << x[parameterIndex] - << std::endl; - } - out << pMesh->elements[0].dimensions.toString() + "\n" + - pMesh->elements[0].material.toString() + " \nA=" - << pMesh->elements[0].A << " \nratio=" - << pMesh->elements[0].dimensions.b / pMesh->elements[0].dimensions.h - << " \naxialRig:" << pMesh->elements[0].rigidity.axial - << " \ntorsionalRig:" << pMesh->elements[0].rigidity.torsional - << " \nfirstBendingRig:" << pMesh->elements[0].rigidity.firstBending - << " \nsecondBendingRig:" << pMesh->elements[0].rigidity.secondBending - << " \nscenario:" + simulationScenarioStrings[simulationScenarioIndex] + - "\n\n"; - out.close(); #ifdef POLYSCOPE_DEFINED - ReducedModelOptimizer::visualizeResults( - global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs, - std::vector{ - static_cast(simulationScenarioIndex)}, - global.reducedToFullInterfaceViMap); + ReducedModelOptimizer::visualizeResults( + global.fullPatternSimulationJobs[simulationScenarioIndex], + global.reducedPatternSimulationJobs[simulationScenarioIndex], + global.reducedToFullInterfaceViMap, false); + + ReducedModelOptimizer::visualizeResults( + global.fullPatternSimulationJobs[simulationScenarioIndex], + std::make_shared( + reducedPatternMaximumDisplacementSimulationJobs + [simulationScenarioIndex]), + global.reducedToFullInterfaceViMap, true); + polyscope::removeAllStructures(); #endif // POLYSCOPE_DEFINED + } + error += thisSimulationScenarioError; + } + // std::ofstream out(filename, std::ios_base::app); + // auto pMesh = + // global + // .reducedPatternSimulationJobs[global.simulationScenarioIndices[0]] + // ->pMesh; + // for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) + // { + // out << "x[" + std::to_string(parameterIndex) + "]=" << + // x[parameterIndex] + // << std::endl; + // } + // out << pMesh->elements[0].dimensions.toString() + "\n" + + // pMesh->elements[0].material.toString() + " \nA=" + // << pMesh->elements[0].A << " \nratio=" + // << pMesh->elements[0].dimensions.b / + // pMesh->elements[0].dimensions.h + // << " \naxialRig:" << pMesh->elements[0].rigidity.axial + // << " \ntorsionalRig:" << pMesh->elements[0].rigidity.torsional + // << " \nfirstBendingRig:" << + // pMesh->elements[0].rigidity.firstBending + // << " \nsecondBendingRig:" << + // pMesh->elements[0].rigidity.secondBending + // << " \nscenario:" + + // simulationScenarioStrings[simulationScenarioIndex] + + // "\n\n"; + // out.close(); } // std::cout << error << std::endl; if (error < global.minY) { @@ -533,69 +578,50 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob( #if POLYSCOPE_DEFINED void ReducedModelOptimizer::visualizeResults( - const std::vector> - &fullPatternSimulationJobs, - const std::vector> - &reducedPatternSimulationJobs, - const std::vector &simulationScenarios, + const std::shared_ptr &pFullPatternSimulationJob, + const std::shared_ptr &pReducedPatternSimulationJob, const std::unordered_map - &reducedToFullInterfaceViMap) { + &reducedToFullInterfaceViMap, + const bool draw = true) { FormFinder simulator; std::shared_ptr pFullPatternSimulationMesh = - fullPatternSimulationJobs[0]->pMesh; + pFullPatternSimulationJob->pMesh; pFullPatternSimulationMesh->registerForDrawing(); - pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() + - "_undeformed.ply"); - reducedPatternSimulationJobs[0]->pMesh->savePly( - reducedPatternSimulationJobs[0]->pMesh->getLabel() + "_undeformed.ply"); - double totalError = 0; - for (const int simulationScenarioIndex : simulationScenarios) { - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationScenarioIndex]; - pFullPatternSimulationJob->registerForDrawing( - pFullPatternSimulationMesh->getLabel()); - SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob); - fullModelResults.registerForDrawing(); - // fullModelResults.saveDeformedModel(); - const std::shared_ptr &pReducedPatternSimulationJob = - reducedPatternSimulationJobs[simulationScenarioIndex]; - SimulationResults reducedModelResults = - simulator.executeSimulation(pReducedPatternSimulationJob); - double interfaceDisplacementNormSum = 0; - for (const auto &interfaceViPair : reducedToFullInterfaceViMap) { - const int fullPatternInterfaceIndex = interfaceViPair.second; - Eigen::Vector3d fullPatternDisplacementVector( - fullModelResults.displacements[fullPatternInterfaceIndex][0], - fullModelResults.displacements[fullPatternInterfaceIndex][1], - fullModelResults.displacements[fullPatternInterfaceIndex][2]); - interfaceDisplacementNormSum += fullPatternDisplacementVector.norm(); - } - reducedModelResults.saveDeformedModel(); - fullModelResults.saveDeformedModel(); - double error = computeError( - reducedModelResults.displacements, fullModelResults.displacements, - interfaceDisplacementNormSum, reducedToFullInterfaceViMap); - std::cout << "Error of simulation scenario " - << simulationScenarioStrings[simulationScenarioIndex] << " is " - << error << std::endl; - totalError += error; - reducedModelResults.registerForDrawing(); - // firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing(); - // registerWorldAxes(); + // pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() + // + + // "_undeformed.ply"); + // reducedPatternSimulationJobs[0]->pMesh->savePly( + // reducedPatternSimulationJobs[0]->pMesh->getLabel() + + // "_undeformed.ply"); + pFullPatternSimulationJob->registerForDrawing( + pFullPatternSimulationMesh->getLabel()); + SimulationResults fullModelResults = + simulator.executeSimulation(pFullPatternSimulationJob); + fullModelResults.registerForDrawing(); + // fullModelResults.saveDeformedModel(); + SimulationResults reducedModelResults = + simulator.executeSimulation(pReducedPatternSimulationJob); + // reducedModelResults.saveDeformedModel(); + // fullModelResults.saveDeformedModel(); + double error = computeRawError( + reducedModelResults.displacements, fullModelResults.displacements, + reducedToFullInterfaceViMap/*, + global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/); + std::cout << "Error is " << error << std::endl; + reducedModelResults.registerForDrawing(); + if (draw) { + polyscope::show(); const std::string screenshotFilename = "/home/iason/Coding/Projects/Approximating shapes with flat " "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" "Images/" + - pFullPatternSimulationMesh->getLabel() + "_" + - simulationScenarioStrings[simulationScenarioIndex]; - polyscope::show(); + pFullPatternSimulationMesh->getLabel() + "_" + "noScenarioName"; + // simulationScenarioStrings[simulationScenarioIndex]; polyscope::screenshot(screenshotFilename, false); fullModelResults.unregister(); reducedModelResults.unregister(); - // firstOptimizationRoundResults[simulationScenarioIndex].unregister(); + pFullPatternSimulationMesh->unregister(); } - std::cout << "Total error:" << totalError << std::endl; } #endif // POLYSCOPE_DEFINED @@ -694,7 +720,19 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) { results.numberOfSimulationCrashes = global.numOfSimulationCrashes; results.x = global.minX; results.objectiveValue = global.minY; + if (settings.shouldNormalizeObjectiveValue) { + global.optimizationSettings.shouldNormalizeObjectiveValue = false; + if (global.optimizeInnerHexagonSize) { + results.objectiveValue = objective(4, results.x.data()); + + } else { + results.objectiveValue = + objective(results.x[0], results.x[1], results.x[2]); + } + global.optimizationSettings.shouldNormalizeObjectiveValue = true; + } results.time = elapsed.count() / 1000.0; + const bool printDebugInfo = false; if (printDebugInfo) { std::cout << "Finished optimizing." << endl; // std::cout << "Solution x:" << endl; @@ -1006,7 +1044,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( global.g_optimalReducedModelDisplacements.resize(6); global.reducedPatternSimulationJobs.resize(6); global.fullPatternDisplacements.resize(6); - global.fullPatternDisplacementNormSum.resize(6); + global.reducedPatternMaximumDisplacementNormSum.resize(6); global.g_firstRoundIterationIndex = 0; global.minY = std::numeric_limits::max(); global.numOfSimulationCrashes = 0; @@ -1014,34 +1052,18 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( global.optimizationSettings = optimizationSettings; global.fullPatternSimulationJobs = createScenarios(m_pFullPatternSimulationMesh); + reducedPatternMaximumDisplacementSimulationJobs.resize(6); // polyscope::removeAllStructures(); - FormFinder::Settings settings; + FormFinder::Settings simulationSettings; // settings.shouldDraw = true; for (int simulationScenarioIndex : global.simulationScenarioIndices) { const std::shared_ptr &pFullPatternSimulationJob = global.fullPatternSimulationJobs[simulationScenarioIndex]; - SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob, settings); + SimulationResults fullModelResults = simulator.executeSimulation( + pFullPatternSimulationJob, simulationSettings); global.fullPatternDisplacements[simulationScenarioIndex] = fullModelResults.displacements; - double interfaceDisplacementNormSum = 0; - for (const auto &interfaceViPair : global.reducedToFullInterfaceViMap) { - const int fullPatternInterfaceIndex = interfaceViPair.second; - Eigen::Vector3d fullPatternDisplacementVector( - fullModelResults.displacements[fullPatternInterfaceIndex][0], - fullModelResults.displacements[fullPatternInterfaceIndex][1], - fullModelResults.displacements[fullPatternInterfaceIndex][2]); - interfaceDisplacementNormSum += fullPatternDisplacementVector.norm(); - } - global.fullPatternDisplacementNormSum[simulationScenarioIndex] = - interfaceDisplacementNormSum; - // global.g_optimalReducedModelDisplacements[simulationScenarioIndex].resize( - // m_pReducedPatternSimulationMesh->VN(), 3); - // computeDesiredReducedModelDisplacements( - // fullModelResults, global.reducedToFullInterfaceViMap, - // global.g_optimalReducedModelDisplacements[simulationScenarioIndex]); - SimulationJob reducedPatternSimulationJob; reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; computeReducedModelSimulationJob(*pFullPatternSimulationJob, @@ -1049,6 +1071,54 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( reducedPatternSimulationJob); global.reducedPatternSimulationJobs[simulationScenarioIndex] = std::make_shared(reducedPatternSimulationJob); + + if (optimizationSettings.shouldNormalizeObjectiveValue) { + // Compute the quantities for normalizing the obj value + const double minB = optimizationSettings.xRanges[0].min; + const double maxRatio = optimizationSettings.xRanges[1].max; + const double minE = optimizationSettings.xRanges[2].min; + const double minHS = 0.3; + std::vector mostFlexibleOptimizationParameters{minB, maxRatio, + minE, minHS}; + if (global.optimizeInnerHexagonSize) { + updateMesh(4, mostFlexibleOptimizationParameters.data()); + } else { + updateMesh(3, mostFlexibleOptimizationParameters.data()); + } + reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex] = + global.reducedPatternSimulationJobs[simulationScenarioIndex] + ->getCopy(); + reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex] + .pMesh->setLabel("reduced_maxDisplacement"); + SimulationResults reducedModelResults = simulator.executeSimulation( + global.reducedPatternSimulationJobs[simulationScenarioIndex], + simulationSettings); + const double errorOfMaxDisplacedReduced = computeRawError( + reducedModelResults.displacements, fullModelResults.displacements, + global.reducedToFullInterfaceViMap); + const double errorOfNonDisplacedReduced = computeRawError( + std::vector(reducedModelResults.displacements.size(), + Vector6d(0)), + fullModelResults.displacements, global.reducedToFullInterfaceViMap); + double displacementMultiplier = 2; + // if (global.optimizeInnerHexagonSize) { + // displacementMultiplier = 2; + // } else { + // displacementMultiplier = 2; + // } + global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex] = + displacementMultiplier * + std::max(errorOfMaxDisplacedReduced, errorOfNonDisplacedReduced); + // if (errorOfMaxDisplacedReduced > errorOfNonDisplacedReduced) { + // std::cout << "Max disp results in a bigger error for scenario " + // << simulationScenarioStrings[simulationScenarioIndex] + // << std::endl; + // } else { + // std::cout << "Zero disp results in a bigger error for scenario " + // << simulationScenarioStrings[simulationScenarioIndex] + // << std::endl; + // } + } } Results optResults = runOptimization(optimizationSettings); for (int simulationScenarioIndex : global.simulationScenarioIndices) { @@ -1061,7 +1131,10 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( // optResults.draw(); // visualizeResults(simulationJobs, global.simulationScenarioIndices); - // visualizeResults(simulationJobs, global.reducedPatternSimulationJobs, - // global.simulationScenarioIndices,global.reducedToFullInterfaceViMap); +#ifdef POLYSCOPE_DEFINED +// visualizeResults( +// global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs, +// global.simulationScenarioIndices, global.reducedToFullInterfaceViMap); +#endif // POLYSCOPE_DEFINED return optResults; } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 66913c4..f3f1f81 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -6,6 +6,7 @@ #include "edgemesh.hpp" #include "elementalmesh.hpp" #include "matplot/matplot.h" +#include "optimizationstructs.hpp" #include using FullPatternVertexIndex = VertexIndex; @@ -32,6 +33,7 @@ public: NumberOfSimulationScenarios }; + struct Results; struct xRange { std::string label; double min; @@ -41,13 +43,12 @@ public: "]"; } }; - struct Results; struct Settings { std::vector xRanges; int numberOfFunctionCalls{100}; double solutionAccuracy{1e-2}; - bool normalizeObjectiveValue{true}; + bool shouldNormalizeObjectiveValue{true}; std::string toString() const { std::string settingsString; @@ -58,9 +59,10 @@ public: } settingsString += xRangesString; } - settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls) + - " Accuracy=" + std::to_string(solutionAccuracy) + - " Norm=" + (normalizeObjectiveValue ? "yes" : "no"); + settingsString += + "FuncCalls=" + std::to_string(numberOfFunctionCalls) + + " Accuracy=" + std::to_string(solutionAccuracy) + + " Norm=" + (shouldNormalizeObjectiveValue ? "yes" : "no"); return settingsString; } @@ -88,9 +90,189 @@ public: os << solutionAccuracy; } }; + struct Results { + double time{-1}; + int numberOfSimulationCrashes{0}; + std::vector x; + double objectiveValue; + std::vector> fullPatternSimulationJobs; + std::vector> reducedPatternSimulationJobs; + + void save(const string &saveToPath) const { + assert(std::filesystem::is_directory(saveToPath)); + + const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); + for (int simulationJobIndex = 0; + simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { + const std::shared_ptr &pFullPatternSimulationJob = + fullPatternSimulationJobs[simulationJobIndex]; + std::filesystem::path simulationJobFolderPath( + std::filesystem::path(saveToPath) + .append(pFullPatternSimulationJob->label)); + std::filesystem::create_directory(simulationJobFolderPath); + const auto fullPatternDirectoryPath = + std::filesystem::path(simulationJobFolderPath).append("Full"); + std::filesystem::create_directory(fullPatternDirectoryPath); + pFullPatternSimulationJob->save(fullPatternDirectoryPath.string()); + const std::shared_ptr &pReducedPatternSimulationJob = + reducedPatternSimulationJobs[simulationJobIndex]; + const auto reducedPatternDirectoryPath = + std::filesystem::path(simulationJobFolderPath).append("Reduced"); + if (!std::filesystem::exists(reducedPatternDirectoryPath)) { + std::filesystem::create_directory(reducedPatternDirectoryPath); + } + pReducedPatternSimulationJob->save( + reducedPatternDirectoryPath.string()); + } + } + void load(const string &loadFromPath) { + assert(std::filesystem::is_directory(loadFromPath)); + + for (const auto &directoryEntry : + filesystem::directory_iterator(loadFromPath)) { + const auto simulationScenarioPath = directoryEntry.path(); + if (!std::filesystem::is_directory(simulationScenarioPath)) { + continue; + } + // Load reduced pattern files + for (const auto &fileEntry : filesystem::directory_iterator( + std::filesystem::path(simulationScenarioPath) + .append("Full"))) { + const auto filepath = fileEntry.path(); + if (filepath.extension() == ".json") { + SimulationJob job; + job.load(filepath.string()); + fullPatternSimulationJobs.push_back( + std::make_shared(job)); + } + } + + // Load full pattern files + for (const auto &fileEntry : filesystem::directory_iterator( + std::filesystem::path(simulationScenarioPath) + .append("Reduced"))) { + const auto filepath = fileEntry.path(); + if (filepath.extension() == ".json") { + SimulationJob job; + job.load(filepath.string()); + reducedPatternSimulationJobs.push_back( + std::make_shared(job)); + } + } + } + } +#if POLYSCOPE_DEFINED + void draw() const { + initPolyscope(); + FormFinder simulator; + assert(fullPatternSimulationJobs.size() == + reducedPatternSimulationJobs.size()); + fullPatternSimulationJobs[0]->pMesh->registerForDrawing(); + reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(); + + const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); + for (int simulationJobIndex = 0; + simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { + // Drawing of full pattern results + const std::shared_ptr &pFullPatternSimulationJob = + fullPatternSimulationJobs[simulationJobIndex]; + pFullPatternSimulationJob->registerForDrawing( + fullPatternSimulationJobs[0]->pMesh->getLabel()); + SimulationResults fullModelResults = + simulator.executeSimulation(pFullPatternSimulationJob); + fullModelResults.registerForDrawing(); + // Drawing of reduced pattern results + const std::shared_ptr &pReducedPatternSimulationJob = + reducedPatternSimulationJobs[simulationJobIndex]; + SimulationResults reducedModelResults = + simulator.executeSimulation(pReducedPatternSimulationJob); + reducedModelResults.registerForDrawing(); + polyscope::show(); + // Save a screensh + // const std::string screenshotFilename = + // "/home/iason/Coding/Projects/Approximating shapes with flat " + // "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" + // + m_pFullPatternSimulationMesh->getLabel() + "_" + + // simulationScenarioStrings[simulationScenarioIndex]; + // polyscope::screenshot(screenshotFilename, false); + fullModelResults.unregister(); + reducedModelResults.unregister(); + // double error = computeError( + // reducedModelResults, + // global.g_optimalReducedModelDisplacements[simulationScenarioIndex]); + // std::cout << "Error of simulation scenario " + // << simulationScenarioStrings[simulationScenarioIndex] << + // " is " + // << error << std::endl; + } + } +#endif // POLYSCOPE_DEFINED + void saveMeshFiles() const { + const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); + assert(numberOfSimulationJobs != 0 && + fullPatternSimulationJobs.size() == + reducedPatternSimulationJobs.size()); + + fullPatternSimulationJobs[0]->pMesh->savePly( + "FullPattern_undeformed.ply"); + reducedPatternSimulationJobs[0]->pMesh->savePly( + "ReducedPattern_undeformed.ply"); + FormFinder simulator; + for (int simulationJobIndex = 0; + simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { + // Drawing of full pattern results + const std::shared_ptr &pFullPatternSimulationJob = + fullPatternSimulationJobs[simulationJobIndex]; + SimulationResults fullModelResults = + simulator.executeSimulation(pFullPatternSimulationJob); + fullModelResults.saveDeformedModel(); + + // Drawing of reduced pattern results + const std::shared_ptr &pReducedPatternSimulationJob = + reducedPatternSimulationJobs[simulationJobIndex]; + SimulationResults reducedModelResults = + simulator.executeSimulation(pReducedPatternSimulationJob); + reducedModelResults.saveDeformedModel(); + } + } + + void + writeHeaderTo(const ReducedModelOptimizer::Settings &settings_optimization, + csvFile &os) { + os << "Obj value"; + for (const ReducedModelOptimizer::xRange &range : + settings_optimization.xRanges) { + os << range.label; + } + os << "Time(s)"; + os << "#Crashes"; + } + + void + writeResultsTo(const ReducedModelOptimizer::Settings &settings_optimization, + csvFile &os) const { + os << objectiveValue; + for (const double &optimalX : x) { + os << optimalX; + } + + for (int unusedXVarCounter = 0; + unusedXVarCounter < settings_optimization.xRanges.size() - x.size(); + unusedXVarCounter++) { + os << "-"; + } + + os << time; + if (numberOfSimulationCrashes == 0) { + os << "-"; + } else { + os << numberOfSimulationCrashes; + } + } + }; inline static const std::string simulationScenarioStrings[] = { - "Axial", "Shear", "Bending", "Double", "Saddle"}; + "Axial", "Shear", "Bending", "Dome", "Saddle"}; Results optimize(const Settings &xRanges, const std::vector &simulationScenarios = std::vector()); @@ -117,7 +299,7 @@ public: std::vector &x); static double objective(double x0, double x1, double x2, double x3); - static double objective(double b, double h, double E); + static double objective(double b, double r, double E); static std::vector> createScenarios(const std::shared_ptr &pMesh, @@ -137,22 +319,28 @@ public: &fullToReducedInterfaceViMap, std::unordered_map &fullPatternOppositeInterfaceViMap); - static void - visualizeResults(const std::vector> - &fullPatternSimulationJobs, - const std::vector> - &reducedPatternSimulationJobs, - const std::vector &simulationScenarios, - const std::unordered_map - &reducedToFullInterfaceViMap); + static void visualizeResults( + const std::shared_ptr &pFullPatternSimulationJob, + const std::shared_ptr &pReducedPatternSimulationJob, + const std::unordered_map + &reducedToFullInterfaceViMap, + const bool dontShow); - static double computeError(const std::vector &reducedPatternResults, - const std::vector &fullPatternResults, - const double &interfaceDisplacementNormSum, - const std::unordered_map - &reducedToFullInterfaceViMap); + static double + computeRawError(const std::vector &reducedPatternDisplacements, + const std::vector &fullPatternDisplacements, + const std::unordered_map + &reducedToFullInterfaceViMap); + + static double + computeError(const std::vector &reducedPatternDisplacements, + const std::vector &fullPatternDisplacements, + const std::unordered_map + &reducedToFullInterfaceViMap, + const double &normalizationFactor); private: void @@ -173,188 +361,8 @@ private: static void initializeOptimizationParameters(const std::shared_ptr &mesh); - static double - computeError(const SimulationResults &reducedPatternResults, - const Eigen::MatrixX3d &optimalReducedPatternDisplacements); static double objective(long n, const double *x); FormFinder simulator; }; -struct ReducedModelOptimizer::Results { - double time{-1}; - int numberOfSimulationCrashes{0}; - std::vector x; - double objectiveValue; - std::vector> fullPatternSimulationJobs; - std::vector> reducedPatternSimulationJobs; - - void save(const string &saveToPath) const { - assert(std::filesystem::is_directory(saveToPath)); - - const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); - for (int simulationJobIndex = 0; - simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationJobIndex]; - std::filesystem::path simulationJobFolderPath( - std::filesystem::path(saveToPath) - .append(pFullPatternSimulationJob->label)); - std::filesystem::create_directory(simulationJobFolderPath); - const auto fullPatternDirectoryPath = - std::filesystem::path(simulationJobFolderPath).append("Full"); - std::filesystem::create_directory(fullPatternDirectoryPath); - pFullPatternSimulationJob->save(fullPatternDirectoryPath.string()); - const std::shared_ptr &pReducedPatternSimulationJob = - reducedPatternSimulationJobs[simulationJobIndex]; - const auto reducedPatternDirectoryPath = - std::filesystem::path(simulationJobFolderPath).append("Reduced"); - if (!std::filesystem::exists(reducedPatternDirectoryPath)) { - std::filesystem::create_directory(reducedPatternDirectoryPath); - } - pReducedPatternSimulationJob->save(reducedPatternDirectoryPath.string()); - } - } - void load(const string &loadFromPath) { - assert(std::filesystem::is_directory(loadFromPath)); - - for (const auto &directoryEntry : - filesystem::directory_iterator(loadFromPath)) { - const auto simulationScenarioPath = directoryEntry.path(); - if (!std::filesystem::is_directory(simulationScenarioPath)) { - continue; - } - // Load reduced pattern files - for (const auto &fileEntry : filesystem::directory_iterator( - std::filesystem::path(simulationScenarioPath).append("Full"))) { - const auto filepath = fileEntry.path(); - if (filepath.extension() == ".json") { - SimulationJob job; - job.load(filepath.string()); - fullPatternSimulationJobs.push_back( - std::make_shared(job)); - } - } - - // Load full pattern files - for (const auto &fileEntry : filesystem::directory_iterator( - std::filesystem::path(simulationScenarioPath) - .append("Reduced"))) { - const auto filepath = fileEntry.path(); - if (filepath.extension() == ".json") { - SimulationJob job; - job.load(filepath.string()); - reducedPatternSimulationJobs.push_back( - std::make_shared(job)); - } - } - } - } -#if POLYSCOPE_DEFINED - void draw() const { - initPolyscope(); - FormFinder simulator; - assert(fullPatternSimulationJobs.size() == - reducedPatternSimulationJobs.size()); - fullPatternSimulationJobs[0]->pMesh->registerForDrawing(); - reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(); - - const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); - for (int simulationJobIndex = 0; - simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { - // Drawing of full pattern results - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationJobIndex]; - pFullPatternSimulationJob->registerForDrawing( - fullPatternSimulationJobs[0]->pMesh->getLabel()); - SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob); - fullModelResults.registerForDrawing(); - // Drawing of reduced pattern results - const std::shared_ptr &pReducedPatternSimulationJob = - reducedPatternSimulationJobs[simulationJobIndex]; - SimulationResults reducedModelResults = - simulator.executeSimulation(pReducedPatternSimulationJob); - reducedModelResults.registerForDrawing(); - polyscope::show(); - // Save a screensh - // const std::string screenshotFilename = - // "/home/iason/Coding/Projects/Approximating shapes with flat " - // "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" - // + m_pFullPatternSimulationMesh->getLabel() + "_" + - // simulationScenarioStrings[simulationScenarioIndex]; - // polyscope::screenshot(screenshotFilename, false); - fullModelResults.unregister(); - reducedModelResults.unregister(); - // double error = computeError( - // reducedModelResults, - // global.g_optimalReducedModelDisplacements[simulationScenarioIndex]); - // std::cout << "Error of simulation scenario " - // << simulationScenarioStrings[simulationScenarioIndex] << " - // is " - // << error << std::endl; - } - } -#endif // POLYSCOPE_DEFINED - void saveMeshFiles() const { - const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); - assert(numberOfSimulationJobs != 0 && - fullPatternSimulationJobs.size() == - reducedPatternSimulationJobs.size()); - - fullPatternSimulationJobs[0]->pMesh->savePly("FullPattern_undeformed.ply"); - reducedPatternSimulationJobs[0]->pMesh->savePly( - "ReducedPattern_undeformed.ply"); - FormFinder simulator; - for (int simulationJobIndex = 0; - simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { - // Drawing of full pattern results - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationJobIndex]; - SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob); - fullModelResults.saveDeformedModel(); - - // Drawing of reduced pattern results - const std::shared_ptr &pReducedPatternSimulationJob = - reducedPatternSimulationJobs[simulationJobIndex]; - SimulationResults reducedModelResults = - simulator.executeSimulation(pReducedPatternSimulationJob); - reducedModelResults.saveDeformedModel(); - } - } - - void - writeHeaderTo(const ReducedModelOptimizer::Settings &settings_optimization, - csvFile &os) { - os << "Obj value"; - for (const ReducedModelOptimizer::xRange &range : - settings_optimization.xRanges) { - os << range.label; - } - os << "Time(s)"; - os << "#Crashes"; - } - - void - writeResultsTo(const ReducedModelOptimizer::Settings &settings_optimization, - csvFile &os) const { - os << objectiveValue; - for (const double &optimalX : x) { - os << optimalX; - } - - for (int unusedXVarCounter = 0; - unusedXVarCounter < settings_optimization.xRanges.size() - x.size(); - unusedXVarCounter++) { - os << "-"; - } - - os << time; - if (numberOfSimulationCrashes == 0) { - os << "-"; - } else { - os << numberOfSimulationCrashes; - } - } -}; #endif // REDUCEDMODELOPTIMIZER_HPP