diff --git a/CMakeLists.txt b/CMakeLists.txt index 4bc7a5d..85d0601 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -37,6 +37,17 @@ download_project(PROJ MYSOURCES endif() + +#Polyscope +download_project(PROJ POLYSCOPE + GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git + GIT_TAG master + PREFIX ${CMAKE_CURRENT_SOURCE_DIR}/build/external/ + ${UPDATE_DISCONNECTED_IF_AVAILABLE} +) +add_subdirectory(${POLYSCOPE_SOURCE_DIR}) +add_compile_definitions(POLYSCOPE_DEFINED) + #dlib download_project(PROJ DLIB GIT_REPOSITORY https://github.com/davisking/dlib.git diff --git a/src/main.cpp b/src/main.cpp index 7831655..07b512b 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,8 +1,10 @@ #include "beamformfinder.hpp" #include "csvfile.hpp" #include "edgemesh.hpp" -#include "externvariables.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" @@ -14,30 +16,23 @@ #include #include -bool printDebugInfo; - int main(int argc, char *argv[]) { - if (argc < 5) { - std::cerr << "Specify at least if D(ebug) or R(elease) and the pattern " - "pair to be optimized." + if (argc < 3) { + std::cerr << "Specify at least the two pattern filepaths to be " + "optimized.Exiting.." << std::endl; std::terminate(); } - if (argv[1] == "D") { - printDebugInfo = true; - } else { - printDebugInfo = false; - } // Populate the pattern pair to be optimized ////Full pattern - const std::string filepath_fullPattern = argv[2]; + const std::string filepath_fullPattern = argv[1]; FlatPattern fullPattern(filepath_fullPattern); fullPattern.setLabel( std::filesystem::path(filepath_fullPattern).stem().string()); fullPattern.scale(0.03); ////Reduced pattern - const std::string filepath_reducedPattern = argv[3]; + const std::string filepath_reducedPattern = argv[2]; FlatPattern reducedPattern(filepath_reducedPattern); reducedPattern.setLabel( std::filesystem::path(filepath_reducedPattern).stem().string()); @@ -59,9 +54,7 @@ int main(int argc, char *argv[]) { // Optimize pair const std::string pairName = fullPattern.getLabel() + "@" + reducedPattern.getLabel(); - if (printDebugInfo) { - std::cout << "Optimizing " << pairName << std::endl; - } + const std::vector numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; ReducedModelOptimizer optimizer(numberOfNodesPerSlot); optimizer.initializePatterns(fullPattern, reducedPattern, {}); @@ -72,28 +65,20 @@ int main(int argc, char *argv[]) { const bool input_resultDirectoryDefined = argc >= 6; std::string optimizationResultsDirectory = input_resultDirectoryDefined ? argv[5] : "OptimizationResults"; - //// Get current date for creating the results folder - std::time_t now = time(0); - std::tm *ltm = std::localtime(&now); - std::string currentDate = std::to_string(ltm->tm_mday) + "_" + - std::to_string(1 + ltm->tm_mon) + "_" + - std::to_string(1900 + ltm->tm_year); - std::filesystem::path optimizationResultsDirectoryDatePath( - std::filesystem::path(optimizationResultsDirectory).append(currentDate)); if (optimizationResults.numberOfSimulationCrashes != 0) { const auto crashedJobsDirPath = - std::filesystem::path(optimizationResultsDirectoryDatePath) + std::filesystem::path(optimizationResultsDirectory) .append("CrashedJobs") .append(pairName); std::filesystem::create_directories(crashedJobsDirPath); optimizationResults.save(crashedJobsDirPath.string()); } else { - std::filesystem::path dirPath_thisOptimization( - std::filesystem::path(optimizationResultsDirectoryDatePath) + std::filesystem::path convergedJobsDirPath( + std::filesystem::path(optimizationResultsDirectory) .append("ConvergedJobs") .append(pairName)); - std::filesystem::create_directories(dirPath_thisOptimization); - optimizationResults.save(dirPath_thisOptimization.string()); + std::filesystem::create_directories(convergedJobsDirPath); + optimizationResults.save(convergedJobsDirPath.string()); } csvFile csv_results({}, false); // csvFile csv_results(std::filesystem::path(dirPath_thisOptimization) diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 852696e..d8623e6 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -10,7 +10,7 @@ struct GlobalOptimizationVariables { std::vector g_optimalReducedModelDisplacements; std::vector> fullPatternDisplacements; std::vector fullPatternDisplacementNormSum; - std::vector g_fullPatternSimulationJob; + std::vector> fullPatternSimulationJobs; std::vector> reducedPatternSimulationJobs; std::unordered_map reducedToFullInterfaceViMap; @@ -34,84 +34,6 @@ struct GlobalOptimizationVariables { int numberOfOptimizationParameters{3}; ReducedModelOptimizer::Settings optimizationSettings; } global; -//#pragma omp threadprivate(global) - -// struct OptimizationCallback { -// double operator()(const size_t &iterations, const Eigen::VectorXd &x, -// const double &fval, Eigen::VectorXd &gradient) const { -// // run simulation -// // SimulationResults reducedModelResults = -// // simulator.executeSimulation(reducedModelSimulationJob); -// // reducedModelResults.draw(reducedModelSimulationJob); -// gObjectiveValueHistory.push_back(fval); -// auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), -// gObjectiveValueHistory.size()); -// gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory); -// // const std::string plotImageFilename = "objectivePlot.png"; -// // matplot::save(plotImageFilename); -// // if (numberOfOptimizationRounds % 30 == 0) { -// // std::filesystem::copy_file( -// // std::filesystem::path(plotImageFilename), -// // std::filesystem::path("objectivePlot_copy.png")); -// // } -// // std::stringstream ss; -// // ss << x; -// // reducedModelResults.simulationLabel = ss.str(); -// // SimulationResultsReporter resultsReporter; -// // resultsReporter.reportResults( -// // {reducedModelResults}, -// // std::filesystem::current_path().append("Results")); - -// return true; -// } -//}; - -// struct Objective { -// double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const { -// assert(x.rows() == 4); - -// // drawSimulationJob(simulationJob); -// // Set mesh from x -// std::shared_ptr reducedModel = -// g_reducedPatternSimulationJob.mesh; -// for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) { -// if (g_reducedPatternExludedEdges.contains(ei)) { -// continue; -// } -// Element &e = reducedModel->elements[ei]; -// e.axialConstFactor = g_initialStiffnessFactors(ei, 0) * x(0); -// e.torsionConstFactor = g_initialStiffnessFactors(ei, 1) * x(1); -// e.firstBendingConstFactor = g_initialStiffnessFactors(ei, 2) * x(2); -// e.secondBendingConstFactor = g_initialStiffnessFactors(ei, 3) * x(3); -// } -// // run simulation -// SimulationResults reducedModelResults = -// simulator.executeSimulation(g_reducedPatternSimulationJob); -// // std::stringstream ss; -// // ss << x; -// // reducedModelResults.simulationLabel = ss.str(); -// // SimulationResultsReporter resultsReporter; -// // resultsReporter.reportResults( -// // {reducedModelResults}, -// // std::filesystem::current_path().append("Results")); -// // compute error and return it -// double error = 0; -// for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) { -// VertexIndex reducedModelVi = reducedFullViPair.first; -// Eigen::Vector3d vertexDisplacement( -// reducedModelResults.displacements[reducedModelVi][0], -// reducedModelResults.displacements[reducedModelVi][1], -// reducedModelResults.displacements[reducedModelVi][2]); -// Eigen::Vector3d errorVector = -// Eigen::Vector3d( -// g_optimalReducedModelDisplacements.row(reducedModelVi)) - -// vertexDisplacement; -// error += errorVector.norm(); -// } - -// return error; -// } -//}; double ReducedModelOptimizer::computeError( const std::vector &reducedPatternDisplacements, @@ -204,12 +126,10 @@ void updateMesh(long n, const double *x) { pReducedPatternSimulationMesh->vert[vi].P() = global.g_innerHexagonVectors[vi / 2] * x[n - 1]; } - - pReducedPatternSimulationMesh->reset(); - pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); - // pReducedPatternSimulationMesh->registerForDrawing("Optimized - // hexagon"); polyscope::show(); } + + pReducedPatternSimulationMesh->reset(); + pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); } double ReducedModelOptimizer::objective(double b, double h, double E) { @@ -252,9 +172,12 @@ double ReducedModelOptimizer::objective(long n, const double *x) { simulationSettings); std::string filename; if (!reducedModelResults.converged) { + error += std::numeric_limits::max(); global.numOfSimulationCrashes++; - } else { +if(beVerbose){ + std::cout << "Failed simulation" << std::endl; + } }else { error += computeError( reducedModelResults.displacements, global.fullPatternDisplacements[simulationScenarioIndex], @@ -282,6 +205,13 @@ double ReducedModelOptimizer::objective(long n, const double *x) { << " \nscenario:" + simulationScenarioStrings[simulationScenarioIndex] + "\n\n"; out.close(); +#ifdef POLYSCOPE_DEFINED + ReducedModelOptimizer::visualizeResults( + global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs, + std::vector{ + static_cast(simulationScenarioIndex)}, + global.reducedToFullInterfaceViMap); +#endif // POLYSCOPE_DEFINED } // std::cout << error << std::endl; if (error < global.minY) { @@ -414,6 +344,61 @@ void ReducedModelOptimizer::computeMaps( assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN()); fullPatternOppositeInterfaceViMap[vi0] = vi1; } + + const bool debugMapping = false; + if (debugMapping) { +#if POLYSCOPE_DEFINED + reducedPattern.registerForDrawing(); + std::vector colors_reducedPatternExcludedEdges( + reducedPattern.EN(), glm::vec3(0, 0, 0)); + for (const size_t ei : global.reducedPatternExludedEdges) { + colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0); + } + const std::string label = reducedPattern.getLabel(); + polyscope::getCurveNetwork(label) + ->addEdgeColorQuantity("Excluded edges", + colors_reducedPatternExcludedEdges) + ->setEnabled(true); + polyscope::show(); + + std::vector nodeColorsOpposite(fullPattern.VN(), + glm::vec3(0, 0, 0)); + for (const std::pair oppositeVerts : + fullPatternOppositeInterfaceViMap) { + auto color = polyscope::getNextUniqueColor(); + nodeColorsOpposite[oppositeVerts.first] = color; + nodeColorsOpposite[oppositeVerts.second] = color; + } + fullPattern.registerForDrawing(); + polyscope::getCurveNetwork(fullPattern.getLabel()) + ->addNodeColorQuantity("oppositeMap", nodeColorsOpposite) + ->setEnabled(true); + polyscope::show(); + + std::vector nodeColorsReducedToFull_reduced(reducedPattern.VN(), + glm::vec3(0, 0, 0)); + std::vector nodeColorsReducedToFull_full(fullPattern.VN(), + glm::vec3(0, 0, 0)); + for (size_t vi = 0; vi < reducedPattern.VN(); vi++) { + if (global.reducedToFullInterfaceViMap.contains(vi)) { + + auto color = polyscope::getNextUniqueColor(); + nodeColorsReducedToFull_reduced[vi] = color; + nodeColorsReducedToFull_full[global.reducedToFullInterfaceViMap[vi]] = + color; + } + } + polyscope::getCurveNetwork(reducedPattern.getLabel()) + ->addNodeColorQuantity("reducedToFull_reduced", + nodeColorsReducedToFull_reduced) + ->setEnabled(true); + polyscope::getCurveNetwork(fullPattern.getLabel()) + ->addNodeColorQuantity("reducedToFull_full", + nodeColorsReducedToFull_full) + ->setEnabled(true); + polyscope::show(); +#endif + } } void ReducedModelOptimizer::computeMaps( @@ -438,6 +423,9 @@ void ReducedModelOptimizer::initializePatterns( // reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel()); assert(fullPattern.VN() == reducedPattern.VN() && fullPattern.EN() >= reducedPattern.EN()); +#if POLYSCOPE_DEFINED + polyscope::removeAllStructures(); +#endif // Create copies of the input models FlatPattern copyFullPattern; FlatPattern copyReducedPattern; @@ -543,6 +531,74 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob( // reducedModelNodalForcedNormals; } +#if POLYSCOPE_DEFINED +void ReducedModelOptimizer::visualizeResults( + const std::vector> + &fullPatternSimulationJobs, + const std::vector> + &reducedPatternSimulationJobs, + const std::vector &simulationScenarios, + const std::unordered_map + &reducedToFullInterfaceViMap) { + FormFinder simulator; + std::shared_ptr pFullPatternSimulationMesh = + fullPatternSimulationJobs[0]->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(); + const std::string screenshotFilename = + "/home/iason/Coding/Projects/Approximating shapes with flat " + "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" + "Images/" + + pFullPatternSimulationMesh->getLabel() + "_" + + simulationScenarioStrings[simulationScenarioIndex]; + polyscope::show(); + polyscope::screenshot(screenshotFilename, false); + fullModelResults.unregister(); + reducedModelResults.unregister(); + // firstOptimizationRoundResults[simulationScenarioIndex].unregister(); + } + std::cout << "Total error:" << totalError << std::endl; +} +#endif // POLYSCOPE_DEFINED + void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, const std::unordered_map &displacementsReducedToFullMap, @@ -947,8 +1003,6 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( SimulationScenario::Saddle}; } - std::vector> simulationJobs = - createScenarios(m_pFullPatternSimulationMesh); global.g_optimalReducedModelDisplacements.resize(6); global.reducedPatternSimulationJobs.resize(6); global.fullPatternDisplacements.resize(6); @@ -958,13 +1012,15 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( global.numOfSimulationCrashes = 0; global.numberOfFunctionCalls = 0; global.optimizationSettings = optimizationSettings; + global.fullPatternSimulationJobs = + createScenarios(m_pFullPatternSimulationMesh); // polyscope::removeAllStructures(); FormFinder::Settings settings; // settings.shouldDraw = true; for (int simulationScenarioIndex : global.simulationScenarioIndices) { const std::shared_ptr &pFullPatternSimulationJob = - simulationJobs[simulationScenarioIndex]; + global.fullPatternSimulationJobs[simulationScenarioIndex]; SimulationResults fullModelResults = simulator.executeSimulation(pFullPatternSimulationJob, settings); global.fullPatternDisplacements[simulationScenarioIndex] = @@ -997,7 +1053,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( Results optResults = runOptimization(optimizationSettings); for (int simulationScenarioIndex : global.simulationScenarioIndices) { optResults.fullPatternSimulationJobs.push_back( - simulationJobs[simulationScenarioIndex]); + global.fullPatternSimulationJobs[simulationScenarioIndex]); optResults.reducedPatternSimulationJobs.push_back( global.reducedPatternSimulationJobs[simulationScenarioIndex]); } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index b4dc7c6..66913c4 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -248,7 +248,52 @@ struct ReducedModelOptimizer::Results { } } } +#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 &&