diff --git a/CMakeLists.txt b/CMakeLists.txt index 1209bb6..525f8db 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -24,8 +24,8 @@ download_project(PROJ MATPLOTPLUSPLUS ) add_subdirectory(${MATPLOTPLUSPLUS_SOURCE_DIR}) -if (EXISTS ${MYSOURCES_SOURCE_DIR}) set(MYSOURCES_SOURCE_DIR "/home/iason/Coding/Libraries/MySources") +if (EXISTS ${MYSOURCES_SOURCE_DIR}) else() ##MySources download_project(PROJ MYSOURCES diff --git a/src/main.cpp b/src/main.cpp index 7b225b1..1eac1ab 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -30,13 +30,13 @@ int main(int argc, char *argv[]) { std::vector CWreducedModelEdges{vcg::Point2i(1, 5), vcg::Point2i(3, 1)}; FlatPattern CWReducedModel(numberOfNodesPerSlot, CWreducedModelEdges); - CWReducedModel.setLabel("CW_reduced"); + CWReducedModel.setLabel("CCW_reduced"); CWReducedModel.scale(0.03); std::vector CCWreducedModelEdges{vcg::Point2i(1, 5), vcg::Point2i(3, 5)}; FlatPattern CCWReducedModel(numberOfNodesPerSlot, CCWreducedModelEdges); - CCWReducedModel.setLabel("CCW_reduced"); + CCWReducedModel.setLabel("CW_reduced"); CCWReducedModel.scale(0.03); std::vector reducedModels{&singleBarReducedModel, @@ -48,7 +48,7 @@ int main(int argc, char *argv[]) { ReducedModelOptimizer::xRange beamE{"E", 0.1, 1.9}; ReducedModelOptimizer::xRange innerHexagonSize{"HexagonSize", 0.1, 0.9}; // Test set of full patterns - std::string fullPatternsTestSetDirectory = "../../TestSet"; + std::string fullPatternsTestSetDirectory = "../TestSet"; if (!std::filesystem::exists( std::filesystem::path(fullPatternsTestSetDirectory))) { std::cerr << "Full pattern directory does not exist: " @@ -77,6 +77,7 @@ int main(int argc, char *argv[]) { pFullPattern->scale(0.03); FlatPattern *pReducedPattern = new FlatPattern(); pReducedPattern->copy(*reducedModels[reducedPatternIndex]); + // pReducedPattern->copy(*reducedModels[2]); patternPairs.push_back(std::make_pair(pFullPattern, pReducedPattern)); } } @@ -90,7 +91,7 @@ int main(int argc, char *argv[]) { // for (settings_optimization.numberOfFunctionCalls = 100; // settings_optimization.numberOfFunctionCalls < 5000; // settings_optimization.numberOfFunctionCalls += 100) { - settings_optimization.numberOfFunctionCalls = 1000; + settings_optimization.numberOfFunctionCalls = 10; const std::string optimizationSettingsString = settings_optimization.toString(); std::filesystem::path thisOptimizationDirectory( @@ -129,47 +130,60 @@ int main(int argc, char *argv[]) { auto runtime_ms = std::chrono::duration_cast(end - start); - // for (int patternPairIndex = 0; patternPairIndex < patternPairs.size(); - // patternPairIndex++) { - // std::filesystem::path fullPatternPath( - // std::filesystem::path(thisOptimizationDirectory) - // .append(patternPairs[patternPairIndex].first->getLabel())); - // std::filesystem::create_directories(std::filesystem::path(fullPatternPath)); + for (int patternPairIndex = 0; patternPairIndex < patternPairs.size(); + patternPairIndex++) { + std::filesystem::path saveToPath( + std::filesystem::path(thisOptimizationDirectory) + .append(patternPairs[patternPairIndex].first->getLabel() + "@" + + patternPairs[patternPairIndex].second->getLabel())); + std::filesystem::create_directories(std::filesystem::path(saveToPath)); - // std::filesystem::path reducedPatternPath(fullPatternPath.append( - // patternPairs[patternPairIndex].second->getLabel())); - // std::filesystem::create_directories( - // std::filesystem::path(reducedPatternPath)); - - // optimizationResults_testSet[patternPairIndex].save(reducedPatternPath); - // } + optimizationResults_testSet[patternPairIndex].save(saveToPath); + // optimizationResults_testSet[patternPairIndex].draw(); + } csvFile statistics(std::filesystem::path(thisOptimizationDirectory) .append("statistics.csv") .string(), false); - for (const auto &range : settings_optimization.xRanges) { - statistics << range.min << range.max; - } - statistics << settings_optimization.numberOfFunctionCalls; - if (totalNumberOfSimulationCrashes == 0) { - statistics << "No crashes"; - } else { - statistics << totalNumberOfSimulationCrashes; + // Write header to csv + statistics << "FullPattern@ReducedPattern" + << "Obj value"; + for (const ReducedModelOptimizer::xRange &range : + settings_optimization.xRanges) { + statistics << range.label; } + statistics << "Time(s)"; + statistics << "#Crashes"; + statistics << endrow; - statistics << totalError; + // Write data for (int patternPairIndex = 0; patternPairIndex < patternPairs.size(); patternPairIndex++) { - statistics << patternPairs[patternPairIndex].first->getLabel(); + statistics << patternPairs[patternPairIndex].first->getLabel() + "@" + + patternPairs[patternPairIndex].second->getLabel(); statistics << optimizationResults_testSet[patternPairIndex].objectiveValue; for (const double &optimalX : optimizationResults_testSet[patternPairIndex].x) { statistics << optimalX; } - } - statistics << runtime_ms.count() / 1000.0; + for (int unusedXVarCounter = 0; + unusedXVarCounter < + settings_optimization.xRanges.size() - + optimizationResults_testSet[patternPairIndex].x.size(); + unusedXVarCounter++) { + statistics << "-"; + } + + statistics << optimizationResults_testSet[patternPairIndex].time; + if (totalNumberOfSimulationCrashes == 0) { + statistics << "No crashes"; + } else { + statistics << totalNumberOfSimulationCrashes; + } + statistics << endrow; + } statistics << endrow; // } diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index b91bb1c..b81fb03 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -10,7 +10,7 @@ struct GlobalOptimizationVariables { std::vector g_optimalReducedModelDisplacements; std::vector g_fullPatternSimulationJob; - std::vector> g_reducedPatternSimulationJob; + std::vector> reducedPatternSimulationJobs; std::unordered_map reducedToFullInterfaceViMap; std::unordered_map @@ -18,7 +18,7 @@ struct GlobalOptimizationVariables { matplot::line_handle gPlotHandle; std::vector gObjectiveValueHistory; Eigen::Vector2d g_initialX; - std::unordered_set g_reducedPatternExludedEdges; + std::unordered_set reducedPatternExludedEdges; Eigen::VectorXd g_initialParameters; std::vector simulationScenarioIndices; @@ -127,7 +127,7 @@ double ReducedModelOptimizer::computeError( const Eigen::MatrixX3d &optimalReducedPatternDisplacements) { double error = 0; auto &global = tls[omp_get_thread_num()]; - for (const auto reducedFullViPair : global.g_reducedToFullViMap) { + for (const auto reducedFullViPair : global.reducedToFullInterfaceViMap) { VertexIndex reducedModelVi = reducedFullViPair.first; // const auto pos = // g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP(); @@ -152,10 +152,43 @@ double ReducedModelOptimizer::computeError( return error; } +double ReducedModelOptimizer::computeError( + 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]); + if (!std::isfinite(reducedVertexDisplacement[0]) || + !std::isfinite(reducedVertexDisplacement[1]) || + !std::isfinite(reducedVertexDisplacement[2])) { + std::terminate(); + } + Eigen::Vector3d fullVertexDisplacement( + fullPatternDisplacements[reducedFullViPair.second][0], + fullPatternDisplacements[reducedFullViPair.second][1], + fullPatternDisplacements[reducedFullViPair.second][2]); + Eigen::Vector3d errorVector = + fullVertexDisplacement - reducedVertexDisplacement; + // error += errorVector.squaredNorm(); + error += errorVector.norm(); + } + return error; +} + void updateMesh(long n, const double *x) { auto &global = tls[omp_get_thread_num()]; std::shared_ptr &pReducedPatternSimulationMesh = - global.g_reducedPatternSimulationJob[global.simulationScenarioIndices[0]] + global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]] ->pMesh; // const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0]; // std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << " @@ -233,8 +266,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // std::cout << "x[" + std::to_string(parameterIndex) + "]=" // << x[parameterIndex] << std::endl; // } - const Element &e = - global.g_reducedPatternSimulationJob[0]->pMesh->elements[0]; + const Element &e = global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " // << e.firstBendingConstFactor << " " << // e.secondBendingConstFactor @@ -252,7 +284,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // simulationSettings.shouldDraw = true; for (const int simulationScenarioIndex : global.simulationScenarioIndices) { SimulationResults reducedModelResults = simulator.executeSimulation( - global.g_reducedPatternSimulationJob[simulationScenarioIndex], + global.reducedPatternSimulationJobs[simulationScenarioIndex], simulationSettings); std::string filename; if (!reducedModelResults.converged /*&& @@ -295,7 +327,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) { std::ofstream out(filename, std::ios_base::app); auto pMesh = global - .g_reducedPatternSimulationJob[global.simulationScenarioIndices[0]] + .reducedPatternSimulationJobs[global.simulationScenarioIndices[0]] ->pMesh; for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { @@ -314,6 +346,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) { "\n\n"; out.close(); } + // std::cout << error << std::endl; if (error < global.minY) { global.minY = error; global.minX.assign(x, x + n); @@ -341,9 +374,47 @@ double ReducedModelOptimizer::objective(long n, const double *x) { return error; } +void ReducedModelOptimizer::createSimulationMeshes( + FlatPattern &fullModel, FlatPattern &reducedModel, + std::shared_ptr &pFullPatternSimulationMesh, + std::shared_ptr &pReducedPatternSimulationMesh) { + if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) { + std::cerr << "Error: A rectangular cross section is expected." << std::endl; + terminate(); + } + // Full pattern + pFullPatternSimulationMesh = std::make_shared(fullModel); + pFullPatternSimulationMesh->setBeamCrossSection( + CrossSectionType{0.002, 0.002}); + pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); + + // Reduced pattern + pReducedPatternSimulationMesh = + std::make_shared(reducedModel); + pReducedPatternSimulationMesh->setBeamCrossSection( + CrossSectionType{0.002, 0.002}); + pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); +} + +void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel, + FlatPattern &reducedModel) { + ReducedModelOptimizer::createSimulationMeshes( + fullModel, reducedModel, m_pFullPatternSimulationMesh, + m_pReducedPatternSimulationMesh); +} + void ReducedModelOptimizer::computeMaps( + const std::unordered_set &reducedModelExcludedEdges, + const std::unordered_map> &slotToNode, FlatPattern &fullPattern, FlatPattern &reducedPattern, - const std::unordered_set &reducedModelExcludedEges) { + std::unordered_map + &reducedToFullInterfaceViMap, + std::unordered_map + &fullToReducedInterfaceViMap, + std::unordered_map + &fullPatternOppositeInterfaceViMap) { + + auto &global = tls[omp_get_thread_num()]; // Compute the offset between the interface nodes const size_t interfaceSlotIndex = 4; // bottom edge assert(slotToNode.find(interfaceSlotIndex) != slotToNode.end() && @@ -368,20 +439,20 @@ void ReducedModelOptimizer::computeMaps( // std::cout << "Dups in fan pair:" << duplicateVerticesPerFanPair << // std::endl; - // Save excluded edges - auto &global = tls[omp_get_thread_num()]; - global.g_reducedPatternExludedEdges.clear(); + // Save excluded edges TODO:this changes the global object. Should this be a + // function parameter? + global.reducedPatternExludedEdges.clear(); const size_t fanSize = 6; const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN(); for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { - for (const size_t ei : reducedModelExcludedEges) { - global.g_reducedPatternExludedEdges.insert( + for (const size_t ei : reducedModelExcludedEdges) { + global.reducedPatternExludedEdges.insert( fanIndex * reducedBaseTriangleNumberOfEdges + ei); } } // Construct reduced->full and full->reduced interface vi map - global.reducedToFullInterfaceViMap.clear(); + reducedToFullInterfaceViMap.clear(); vcg::tri::Allocator::PointerUpdater pu_reducedModel; reducedPattern.deleteDanglingVertices(pu_reducedModel); @@ -391,36 +462,30 @@ void ReducedModelOptimizer::computeMaps( reducedPattern.VN() - 1 /*- reducedModelBaseTriangleInterfaceVi*/; reducedPattern.createFan(); for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { - global.reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * - fanIndex + - reducedModelBaseTriangleInterfaceVi] = + reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex + + reducedModelBaseTriangleInterfaceVi] = fullModelBaseTriangleInterfaceVi + fanIndex * fullPatternInterfaceVertexOffset; } - m_fullToReducedInterfaceViMap.clear(); - constructInverseMap(global.reducedToFullInterfaceViMap, - m_fullToReducedInterfaceViMap); + fullToReducedInterfaceViMap.clear(); + constructInverseMap(reducedToFullInterfaceViMap, fullToReducedInterfaceViMap); - // fullPattern.setLabel("FullPattern"); - // reducedPattern.setLabel("ReducedPattern"); // Create opposite vertex map - m_fullPatternOppositeInterfaceViMap.clear(); + fullPatternOppositeInterfaceViMap.clear(); for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) { const size_t vi0 = fullModelBaseTriangleInterfaceVi + fanIndex * fullPatternInterfaceVertexOffset; const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset; assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN()); - m_fullPatternOppositeInterfaceViMap[vi0] = vi1; + fullPatternOppositeInterfaceViMap[vi0] = vi1; } - global.g_reducedToFullViMap = global.reducedToFullInterfaceViMap; - const bool debugMapping = false; if (debugMapping) { reducedPattern.registerForDrawing(); std::vector colors_reducedPatternExcludedEdges( reducedPattern.EN(), glm::vec3(0, 0, 0)); - for (const size_t ei : global.g_reducedPatternExludedEdges) { + for (const size_t ei : global.reducedPatternExludedEdges) { colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0); } const std::string label = reducedPattern.getLabel(); @@ -433,7 +498,7 @@ void ReducedModelOptimizer::computeMaps( std::vector nodeColorsOpposite(fullPattern.VN(), glm::vec3(0, 0, 0)); for (const std::pair oppositeVerts : - m_fullPatternOppositeInterfaceViMap) { + fullPatternOppositeInterfaceViMap) { auto color = polyscope::getNextUniqueColor(); nodeColorsOpposite[oppositeVerts.first] = color; nodeColorsOpposite[oppositeVerts.second] = color; @@ -469,21 +534,14 @@ void ReducedModelOptimizer::computeMaps( } } -void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel, - FlatPattern &reducedModel) { - if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) { - std::cerr << "Error: A rectangular cross section is expected." << std::endl; - terminate(); - } - m_pReducedPatternSimulationMesh = - std::make_shared(reducedModel); - m_pReducedPatternSimulationMesh->setBeamCrossSection( - CrossSectionType{0.002, 0.002}); - m_pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); - m_pFullPatternSimulationMesh = std::make_shared(fullModel); - m_pFullPatternSimulationMesh->setBeamCrossSection( - CrossSectionType{0.002, 0.002}); - m_pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); +void ReducedModelOptimizer::computeMaps( + FlatPattern &fullPattern, FlatPattern &reducedPattern, + const std::unordered_set &reducedModelExcludedEdges) { + auto &global = tls[omp_get_thread_num()]; + ReducedModelOptimizer::computeMaps( + reducedModelExcludedEdges, slotToNode, fullPattern, reducedPattern, + global.reducedToFullInterfaceViMap, m_fullToReducedInterfaceViMap, + m_fullPatternOppositeInterfaceViMap); } ReducedModelOptimizer::ReducedModelOptimizer( @@ -636,8 +694,6 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization( global.gObjectiveValueHistory.clear(); - const size_t n = global.g_initialParameters.rows(); - assert(n == 3); // g_optimizeInnerHexagonSize ? 5: 4; // const size_t npt = (n + 1) * (n + 2) / 2; // // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; @@ -701,8 +757,13 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization( std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy); } auto end = std::chrono::system_clock::now(); - auto elapsed = std::chrono::duration_cast(end - start); - Results results{global.numOfSimulationCrashes, global.minX, global.minY}; + auto elapsed = + std::chrono::duration_cast(end - start); + Results results; + results.numberOfSimulationCrashes = global.numOfSimulationCrashes; + results.x = global.minX; + results.objectiveValue = global.minY; + results.time = elapsed.count() / 1000.0; std::cout << "Finished optimizing." << endl; // std::cout << "Solution x:" << endl; // std::cout << result.x << endl; @@ -1000,28 +1061,36 @@ ReducedModelOptimizer::createScenarios( void ReducedModelOptimizer::visualizeResults( const std::vector> &fullPatternSimulationJobs, - const std::vector &simulationScenarios) { - m_pFullPatternSimulationMesh->registerForDrawing(); + const std::vector> + &reducedPatternSimulationJobs, + const std::vector &simulationScenarios, + const std::unordered_map + &reducedToFullInterfaceViMap) { + FormFinder simulator; + std::shared_ptr pFullPatternSimulationMesh = + fullPatternSimulationJobs[0]->pMesh; + pFullPatternSimulationMesh->registerForDrawing(); + double totalError = 0; for (const int simulationScenarioIndex : simulationScenarios) { const std::shared_ptr &pFullPatternSimulationJob = fullPatternSimulationJobs[simulationScenarioIndex]; pFullPatternSimulationJob->registerForDrawing( - m_pFullPatternSimulationMesh->getLabel()); + pFullPatternSimulationMesh->getLabel()); SimulationResults fullModelResults = simulator.executeSimulation(pFullPatternSimulationJob); fullModelResults.registerForDrawing(); fullModelResults.saveDeformedModel(); - auto &global = tls[omp_get_thread_num()]; const std::shared_ptr &pReducedPatternSimulationJob = - global.g_reducedPatternSimulationJob[simulationScenarioIndex]; + reducedPatternSimulationJobs[simulationScenarioIndex]; SimulationResults reducedModelResults = simulator.executeSimulation(pReducedPatternSimulationJob); - double error = computeError( - reducedModelResults, - global.g_optimalReducedModelDisplacements[simulationScenarioIndex]); + double error = computeError(reducedModelResults.displacements, + fullModelResults.displacements, + reducedToFullInterfaceViMap); std::cout << "Error of simulation scenario " << simulationScenarioStrings[simulationScenarioIndex] << " is " << error << std::endl; + totalError += error; reducedModelResults.registerForDrawing(); // firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing(); // reducedModelResults.saveDeformedModel(); @@ -1029,8 +1098,8 @@ void ReducedModelOptimizer::visualizeResults( const std::string screenshotFilename = "/home/iason/Coding/Projects/Approximating shapes with flat " "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" - "Images" + - m_pFullPatternSimulationMesh->getLabel() + "_" + + "Images/" + + pFullPatternSimulationMesh->getLabel() + "_" + simulationScenarioStrings[simulationScenarioIndex]; polyscope::show(); polyscope::screenshot(screenshotFilename, false); @@ -1038,6 +1107,7 @@ void ReducedModelOptimizer::visualizeResults( reducedModelResults.unregister(); // firstOptimizationRoundResults[simulationScenarioIndex].unregister(); } + std::cout << "Total error:" << totalError << std::endl; } ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( @@ -1056,7 +1126,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( std::vector> simulationJobs = createScenarios(m_pFullPatternSimulationMesh); global.g_optimalReducedModelDisplacements.resize(6); - global.g_reducedPatternSimulationJob.resize(6); + global.reducedPatternSimulationJobs.resize(6); global.g_firstRoundIterationIndex = 0; global.minY = std::numeric_limits::max(); global.numOfSimulationCrashes = 0; @@ -1072,18 +1142,28 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( global.g_optimalReducedModelDisplacements[simulationScenarioIndex].resize( m_pReducedPatternSimulationMesh->VN(), 3); computeDesiredReducedModelDisplacements( - fullModelResults, global.g_reducedToFullViMap, + fullModelResults, global.reducedToFullInterfaceViMap, global.g_optimalReducedModelDisplacements[simulationScenarioIndex]); SimulationJob reducedPatternSimulationJob; reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; computeReducedModelSimulationJob(*pFullPatternSimulationJob, m_fullToReducedInterfaceViMap, reducedPatternSimulationJob); - global.g_reducedPatternSimulationJob[simulationScenarioIndex] = + global.reducedPatternSimulationJobs[simulationScenarioIndex] = std::make_shared(reducedPatternSimulationJob); } Results optResults = runOptimization(xRanges, &objective); - updateMesh(optResults.x.size(), optResults.x.data()); - visualizeResults(simulationJobs, global.simulationScenarioIndices); + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + optResults.fullPatternSimulationJobs.push_back( + simulationJobs[simulationScenarioIndex]); + optResults.reducedPatternSimulationJobs.push_back( + global.reducedPatternSimulationJobs[simulationScenarioIndex]); + } + // updateMesh(optResults.x.size(), optResults.x.data()); + // optResults.draw(); + + // visualizeResults(simulationJobs, global.simulationScenarioIndices); + // visualizeResults(simulationJobs, global.g_reducedPatternSimulationJob, + // global.g_simulationScenarioIndices); return optResults; } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 535f2c6..f2ceaef 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -31,11 +31,7 @@ public: Saddle, NumberOfSimulationScenarios }; - struct Results { - int numberOfSimulationCrashes{0}; - std::vector x; - double objectiveValue; - }; + struct xRange { std::string label; double min; @@ -45,6 +41,7 @@ public: "]"; } }; + struct Results; struct Settings { std::vector xRanges; @@ -120,6 +117,40 @@ public: static double objective(double x0, double x1, double x2, double x3); static double objective(double b, double h, double E); + static std::vector> + createScenarios(const std::shared_ptr &pMesh, + const std::unordered_map + &fullPatternOppositeInterfaceViMap); + static void createSimulationMeshes( + FlatPattern &fullModel, FlatPattern &reducedModel, + std::shared_ptr &pFullPatternSimulationMesh, + std::shared_ptr &pReducedPatternSimulationMesh); + static void computeMaps( + const std::unordered_set &reducedModelExcludedEdges, + const std::unordered_map> &slotToNode, + FlatPattern &fullPattern, FlatPattern &reducedPattern, + std::unordered_map + &reducedToFullInterfaceViMap, + std::unordered_map + &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 double computeError(const std::vector &reducedPatternResults, + const std::vector &fullPatternResults, + const std::unordered_map + &reducedToFullInterfaceViMap); + private: void visualizeResults(const std::vector> @@ -147,5 +178,116 @@ private: 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 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; + } + } + 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("FullPattern"); + std::filesystem::create_directory(fullPatternDirectoryPath); + pFullPatternSimulationJob->save(fullPatternDirectoryPath.string()); + const std::shared_ptr &pReducedPatternSimulationJob = + reducedPatternSimulationJobs[simulationJobIndex]; + const auto reducedPatternDirectoryPath = + std::filesystem::path(simulationJobFolderPath) + .append("ReducedPattern"); + 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(); + // Load reduced pattern files + for (const auto &fileEntry : filesystem::directory_iterator( + std::filesystem::path(simulationScenarioPath) + .append("FullPattern"))) { + 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("ReducedPattern"))) { + const auto filepath = fileEntry.path(); + if (filepath.extension() == ".json") { + SimulationJob job; + job.load(filepath.string()); + reducedPatternSimulationJobs.push_back( + std::make_shared(job)); + } + } + } + } +}; #endif // REDUCEDMODELOPTIMIZER_HPP diff --git a/src/reducedmodeloptimizer.hpp.autosave b/src/reducedmodeloptimizer.hpp.autosave deleted file mode 100644 index 242b8af..0000000 --- a/src/reducedmodeloptimizer.hpp.autosave +++ /dev/null @@ -1,151 +0,0 @@ -#ifndef REDUCEDMODELOPTIMIZER_HPP -#define REDUCEDMODELOPTIMIZER_HPP - -#include "beamformfinder.hpp" -#include "csvfile.hpp" -#include "edgemesh.hpp" -#include "elementalmesh.hpp" -#include "matplot/matplot.h" -#include - -using FullPatternVertexIndex = VertexIndex; -using ReducedPatternVertexIndex = VertexIndex; - -class ReducedModelOptimizer { - std::shared_ptr m_pReducedPatternSimulationMesh; - std::shared_ptr m_pFullPatternSimulationMesh; - std::unordered_map - m_fullToReducedInterfaceViMap; - std::unordered_map - m_fullPatternOppositeInterfaceViMap; - std::unordered_map nodeToSlot; - std::unordered_map> slotToNode; - std::vector initialGuess; - -public: - enum SimulationScenario { - Axial, - Shear, - Bending, - Dome, - Saddle, - NumberOfSimulationScenarios - }; - struct Results { - int numberOfSimulationCrashes{0}; - std::vector x; - double objectiveValue; - }; - struct xRange { - std::string label; - double min; - double max; - std::string toString() const { - return label + "=[" + std::to_string(min) + "," + std::to_string(max) + - "]"; - } - }; - - struct Settings { - std::vector xRanges; - int numberOfFunctionCalls{100}; - double solutionAccuracy{1e-5}; - - std::string toString() const { - std::string settingsString; - if (!xRanges.empty()) { - std::string xRangesString; - for (const xRange &range : xRanges) { - xRangesString += range.toString() + " "; - } - settingsString += xRangesString; - } - settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls) + - " Accuracy=" + std::to_string(solutionAccuracy); - - return settingsString; - } - - void writeTo(csvFile &csv) const { - // Create settings csv header - if (!xRanges.empty()) { - for (const xRange &range : xRanges) { - csv << range.label + " min"; - csv << range.label + " max"; - } - } - csv << "Function Calls"; - csv << "Solution Accuracy"; - csv << endrow; - - if (!xRanges.empty()) { - for (const xRange &range : xRanges) { - csv << range.min; - csv << range.max; - } - } - csv << numberOfFunctionCalls; - csv << solutionAccuracy; - csv << endrow; - } - }; - - inline static const std::string simulationScenarioStrings[] = { - "Axial", "Shear", "Bending", "Double", "Saddle"}; - Results optimize(const Settings &xRanges, - const std::vector &simulationScenarios = - std::vector()); - double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; - - ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot); - static void computeReducedModelSimulationJob( - const SimulationJob &simulationJobOfFullModel, - const std::unordered_map &simulationJobFullToReducedMap, - SimulationJob &simulationJobOfReducedModel); - - SimulationJob - getReducedSimulationJob(const SimulationJob &fullModelSimulationJob); - - void initializePatterns( - FlatPattern &fullPattern, FlatPattern &reducedPatterm, - const std::unordered_set &reducedModelExcludedEges); - - void setInitialGuess(std::vector v); - - static void runBeamOptimization(); - - static void runSimulation(const std::string &filename, - std::vector &x); - - static double objective(double x0, double x1, double x2, double x3); - static double objective(double b, double h, double E); - -private: - void - visualizeResults(const std::vector> - &fullPatternSimulationJobs, - const std::vector &simulationScenarios); - static void computeDesiredReducedModelDisplacements( - const SimulationResults &fullModelResults, - const std::unordered_map &displacementsReducedToFullMap, - Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); - static Results runOptimization(const Settings &settings, - double (*pObjectiveFunction)(long, - const double *)); - std::vector> - createScenarios(const std::shared_ptr &pMesh); - void computeMaps(FlatPattern &fullModel, FlatPattern &reducedPattern, - const std::unordered_set &reducedModelExcludedEges); - void createSimulationMeshes(FlatPattern &fullModel, - FlatPattern &reducedModel); - 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; -}; - -#endif // REDUCEDMODELOPTIMIZER_HPP