From 858282c859dc4f5b1d47f16c6f02c6fcc52aa3fb Mon Sep 17 00:00:00 2001 From: Iason Date: Mon, 1 Feb 2021 16:10:24 +0200 Subject: [PATCH] Parallelization on the set of valid patterns. --- src/main.cpp | 85 ++++---- src/reducedmodeloptimizer.cpp | 351 ++++++++++++++++++---------------- src/reducedmodeloptimizer.hpp | 1 + 3 files changed, 232 insertions(+), 205 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index b11ee25..b8be993 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -11,10 +11,10 @@ #include #include #include +#include #include #include #include -#include int main(int argc, char *argv[]) { @@ -43,11 +43,35 @@ int main(int argc, char *argv[]) { std::vector reducedModels{&singleBarReducedModel, &CWReducedModel, &CCWReducedModel}; - ReducedModelOptimizer optimizer(numberOfNodesPerSlot); + // Test set of full patterns + std::string fullPatternsTestSetDirectory = "../TestSet"; + // "/home/iason/Documents/PhD/Research/Approximating shapes with flat " + // "patterns/Pattern_enumerator/Results/1v_0v_2e_1e_1c_6fan/3/Valid"; + std::vector> patternPairs; + for (const auto &entry : + filesystem::directory_iterator(fullPatternsTestSetDirectory)) { + const auto filepath = + // std::filesystem::path(fullPatternsTestSetDirectory).append("305.ply"); + entry.path(); + const std::string filepathString = filepath.string(); + const std::string tiledSuffix = "_tiled.ply"; + if (filepathString.compare(filepathString.size() - tiledSuffix.size(), + tiledSuffix.size(), tiledSuffix) == 0) { + continue; + } + + FlatPattern *pFullPattern = new FlatPattern(filepathString); + pFullPattern->setLabel(filepath.stem().string()); + pFullPattern->scale(0.03); + FlatPattern *pReducedPattern = new FlatPattern(); + pReducedPattern->copy(*reducedModels[0]); + patternPairs.push_back(std::make_pair(pFullPattern, pReducedPattern)); + } + // for (double rangeOffset = 0.15; rangeOffset <= 0.95; rangeOffset += 0.05) // { ReducedModelOptimizer::Settings settings; - for (settings.maxSimulations = 2200; settings.maxSimulations < 5000; + for (settings.maxSimulations = 2600; settings.maxSimulations < 5000; settings.maxSimulations += 200) { ReducedModelOptimizer::xRange beamWidth{"B", 0.5, 1.5}; ReducedModelOptimizer::xRange beamDimensionsRatio{"bOverh", 0.7, 1.3}; @@ -68,40 +92,17 @@ int main(int argc, char *argv[]) { double totalError = 0; int totalNumberOfSimulationCrashes = 0; - std::vector errors; - std::string fullPatternsTestSetDirectory = "../TestSet"; - // "/home/iason/Documents/PhD/Research/Approximating shapes with flat " - // "patterns/Pattern_enumerator/Results/1v_0v_2e_1e_1c_6fan/3/Valid"; std::vector> - resultsPerPattern; + resultsPerPattern(patternPairs.size()); auto start = std::chrono::high_resolution_clock::now(); - int patternsOptimized = 0; - //extract paths - std::vector fullPatternSetFilenames; - for (const auto& entry : - filesystem::directory_iterator(fullPatternsTestSetDirectory)) { - const auto filepath = - // std::filesystem::path(fullPatternsTestSetDirectory).append("305.ply"); - entry.path(); - fullPatternSetFilenames.push_back(filepath); - } - -//#pragma omp parallel for //schedule(static) num_threads(8) - for (int fullPatternIndex = 0; fullPatternIndex < fullPatternSetFilenames.size();fullPatternIndex++) { - const auto& filepath = fullPatternSetFilenames[fullPatternIndex]; - const std::string& filepathString = fullPatternSetFilenames[fullPatternIndex].string(); - //const auto filepathString = filepath.string(); +#pragma omp parallel for // schedule(static) num_threads(8) + for (int patternPairIndex = 0; patternPairIndex < patternPairs.size(); + patternPairIndex++) { + FlatPattern *pPattern = patternPairs[patternPairIndex].first; + // const auto filepathString = filepath.string(); // Use only the base triangle version - const std::string tiledSuffix = "_tiled.ply"; - if (filepathString.compare(filepathString.size() - tiledSuffix.size(), - tiledSuffix.size(), tiledSuffix) == 0) { - continue; - } // std::cout << "Full pattern:" << filepathString << std::endl; - FlatPattern pattern(filepathString); - pattern.setLabel(filepath.stem().string()); - pattern.scale(0.03); // for (int reducedPatternIndex = 0; // reducedPatternIndex < reducedModels.size(); // reducedPatternIndex++) { @@ -117,7 +118,10 @@ int main(int argc, char *argv[]) { // } // FlatPattern cp; // cp.copy(*reducedModels[0]); - optimizer.initializePatterns(pattern, *reducedModels[0], + const std::vector numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; + ReducedModelOptimizer optimizer(numberOfNodesPerSlot); + optimizer.initializePatterns(*pPattern, + *patternPairs[patternPairIndex].second, optimizationExcludedEi); // optimizer.optimize({ReducedModelOptimizer::Axial}); ReducedModelOptimizer::Results optimizationResults = @@ -139,11 +143,16 @@ int main(int argc, char *argv[]) { // thisOptimizationStatistics << endrow; totalError += optimizationResults.objectiveValue; - resultsPerPattern.push_back( - std::make_pair(filepath.stem().string(), optimizationResults)); + resultsPerPattern[patternPairIndex] = + std::make_pair(pPattern->getLabel(), optimizationResults); totalNumberOfSimulationCrashes += optimizationResults.numberOfSimulationCrashes; - std::cout << "Have optimized " <<++patternsOptimized<<"/"<< static_cast(std::distance(std::filesystem::directory_iterator(fullPatternsTestSetDirectory),std::filesystem::directory_iterator()))<<" patterns." << std::endl; + // std::cout << "Have optimized " << ++patternsOptimized << "/" + // << static_cast( + // std::distance(std::filesystem::directory_iterator( + // fullPatternsTestSetDirectory), + // std::filesystem::directory_iterator())) + // << " patterns." << std::endl; // } } auto end = std::chrono::high_resolution_clock::now(); @@ -177,5 +186,9 @@ int main(int argc, char *argv[]) { statistics << endrow; } + for(auto patternPair:patternPairs){ + delete patternPair.first; + delete patternPair.second; + } return 0; } diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 86346f8..e55c567 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -7,35 +7,35 @@ #include #include -const bool gShouldDraw = true; - -FormFinder simulator; -std::vector g_optimalReducedModelDisplacements; -std::vector g_fullPatternSimulationJob; -std::vector> g_reducedPatternSimulationJob; -std::unordered_map - reducedToFullInterfaceViMap; -std::unordered_map - g_reducedToFullViMap; -matplot::line_handle gPlotHandle; -std::vector gObjectiveValueHistory; -Eigen::Vector2d g_initialX; -std::unordered_set g_reducedPatternExludedEdges; -// double g_initialParameters; -Eigen::VectorXd g_initialParameters; -std::vector - g_simulationScenarioIndices; -std::vector g_innerHexagonVectors{6, VectorType(0, 0, 0)}; -double g_innerHexagonInitialPos = 0; -bool g_optimizeInnerHexagonSize{false}; -std::vector firstOptimizationRoundResults; -int g_firstRoundIterationIndex{0}; -double minY{std::numeric_limits::max()}; -std::vector minX; -std::vector> failedSimulationsXRatio; -int numOfSimulationCrashes{false}; -int numberOfFunctionCalls{ 0 }; - +struct GlobalOptimizationVariables { + const bool gShouldDraw = true; + std::vector g_optimalReducedModelDisplacements; + std::vector g_fullPatternSimulationJob; + std::vector> g_reducedPatternSimulationJob; + std::unordered_map + reducedToFullInterfaceViMap; + std::unordered_map + g_reducedToFullViMap; + matplot::line_handle gPlotHandle; + std::vector gObjectiveValueHistory; + Eigen::Vector2d g_initialX; + std::unordered_set g_reducedPatternExludedEdges; + // double g_initialParameters; + Eigen::VectorXd g_initialParameters; + std::vector + g_simulationScenarioIndices; + std::vector g_innerHexagonVectors{6, VectorType(0, 0, 0)}; + double g_innerHexagonInitialPos = 0; + bool g_optimizeInnerHexagonSize{false}; + std::vector firstOptimizationRoundResults; + int g_firstRoundIterationIndex{0}; + double minY{std::numeric_limits::max()}; + std::vector minX; + std::vector> failedSimulationsXRatio; + int numOfSimulationCrashes{false}; + int numberOfFunctionCalls{0}; +} global; +#pragma omp threadprivate(global) // struct OptimizationCallback { // double operator()(const size_t &iterations, const Eigen::VectorXd &x, // const double &fval, Eigen::VectorXd &gradient) const { @@ -117,7 +117,7 @@ double ReducedModelOptimizer::computeError( const SimulationResults &reducedPatternResults, const Eigen::MatrixX3d &optimalReducedPatternDisplacements) { double error = 0; - for (const auto reducedFullViPair : g_reducedToFullViMap) { + for (const auto reducedFullViPair : global.g_reducedToFullViMap) { VertexIndex reducedModelVi = reducedFullViPair.first; // const auto pos = // g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP(); @@ -144,7 +144,9 @@ double ReducedModelOptimizer::computeError( void updateMesh(long n, const double *x) { std::shared_ptr &pReducedPatternSimulationMesh = - g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]->pMesh; + global + .g_reducedPatternSimulationJob[global.g_simulationScenarioIndices[0]] + ->pMesh; // const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0]; // std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << " // " @@ -158,11 +160,12 @@ void updateMesh(long n, const double *x) { // e.properties.E = g_initialParameters * x[ei]; // e.properties.E = g_initialParameters(0) * x[0]; // e.properties.G = g_initialParameters(1) * x[1]; - e.setDimensions(RectangularBeamDimensions( - g_initialParameters(0) * x[0], - g_initialParameters(0) * x[0] / (g_initialParameters(1) * x[1]))); + e.setDimensions( + RectangularBeamDimensions(global.g_initialParameters(0) * x[0], + global.g_initialParameters(0) * x[0] / + (global.g_initialParameters(1) * x[1]))); e.setMaterial(ElementMaterial(e.material.poissonsRatio, - g_initialParameters(2) * x[2])); + global.g_initialParameters(2) * x[2])); // e.properties.A = g_initialParameters(0) * x[0]; // e.properties.J = g_initialParameters(1) * x[1]; // e.properties.I2 = g_initialParameters(2) * x[2]; @@ -186,12 +189,12 @@ void updateMesh(long n, const double *x) { // e.secondBendingConstFactor // << std::endl; - if (g_optimizeInnerHexagonSize) { + if (global.g_optimizeInnerHexagonSize) { assert(pReducedPatternSimulationMesh->EN() == 12); for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN(); vi += 2) { pReducedPatternSimulationMesh->vert[vi].P() = - g_innerHexagonVectors[vi / 2] * x[n - 1]; + global.g_innerHexagonVectors[vi / 2] * x[n - 1]; } pReducedPatternSimulationMesh->reset(); @@ -219,7 +222,8 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // std::cout << "x[" + std::to_string(parameterIndex) + "]=" // << x[parameterIndex] << std::endl; // } - const Element &e = g_reducedPatternSimulationJob[0]->pMesh->elements[0]; + const Element &e = + global.g_reducedPatternSimulationJob[0]->pMesh->elements[0]; // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " // << e.firstBendingConstFactor << " " << // e.secondBendingConstFactor @@ -232,11 +236,12 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // run simulations double error = 0; + FormFinder simulator; FormFinder::Settings simulationSettings; // simulationSettings.shouldDraw = true; - for (const int simulationScenarioIndex : g_simulationScenarioIndices) { + for (const int simulationScenarioIndex : global.g_simulationScenarioIndices) { SimulationResults reducedModelResults = simulator.executeSimulation( - g_reducedPatternSimulationJob[simulationScenarioIndex], + global.g_reducedPatternSimulationJob[simulationScenarioIndex], simulationSettings); std::string filename; if (!reducedModelResults.converged /*&& @@ -267,18 +272,21 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // simulator.executeSimulation( // g_reducedPatternSimulationJob[simulationScenarioIndex], // simulationSettings); - numOfSimulationCrashes++; + global.numOfSimulationCrashes++; } else { error += computeError( reducedModelResults, - g_optimalReducedModelDisplacements[simulationScenarioIndex]); + global.g_optimalReducedModelDisplacements[simulationScenarioIndex]); filename = "/home/iason/Coding/Projects/Approximating shapes with flat " "patterns/RodModelOptimizationForPatterns/build/" "ProblematicSimulationJobs/conv_dimensions.txt"; } std::ofstream out(filename, std::ios_base::app); auto pMesh = - g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]->pMesh; + global + .g_reducedPatternSimulationJob[global + .g_simulationScenarioIndices[0]] + ->pMesh; for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { out << "x[" + std::to_string(parameterIndex) + "]=" << x[parameterIndex] @@ -296,16 +304,17 @@ double ReducedModelOptimizer::objective(long n, const double *x) { "\n\n"; out.close(); } - if (error < minY) { - minY = error; - minX.assign(x, x + n); + if (error < global.minY) { + global.minY = error; + global.minX.assign(x, x + n); } - if(++numberOfFunctionCalls%50==0){ - std::cout << "Number of function calls:"< colors(gObjectiveValueHistory.size(), 2); @@ -349,18 +358,18 @@ void ReducedModelOptimizer::computeMaps( // std::endl; // Save excluded edges - g_reducedPatternExludedEdges.clear(); + global.g_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) { - g_reducedPatternExludedEdges.insert( + global.g_reducedPatternExludedEdges.insert( fanIndex * reducedBaseTriangleNumberOfEdges + ei); } } // Construct reduced->full and full->reduced interface vi map - reducedToFullInterfaceViMap.clear(); + global.reducedToFullInterfaceViMap.clear(); vcg::tri::Allocator::PointerUpdater pu_reducedModel; reducedPattern.deleteDanglingVertices(pu_reducedModel); @@ -370,13 +379,14 @@ void ReducedModelOptimizer::computeMaps( reducedPattern.VN() - 1 /*- reducedModelBaseTriangleInterfaceVi*/; reducedPattern.createFan(); for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { - reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex + - reducedModelBaseTriangleInterfaceVi] = + global.reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * + fanIndex + + reducedModelBaseTriangleInterfaceVi] = fullModelBaseTriangleInterfaceVi + fanIndex * fullPatternInterfaceVertexOffset; } m_fullToReducedInterfaceViMap.clear(); - constructInverseMap(reducedToFullInterfaceViMap, + constructInverseMap(global.reducedToFullInterfaceViMap, m_fullToReducedInterfaceViMap); // fullPattern.setLabel("FullPattern"); @@ -391,14 +401,14 @@ void ReducedModelOptimizer::computeMaps( m_fullPatternOppositeInterfaceViMap[vi0] = vi1; } - g_reducedToFullViMap = reducedToFullInterfaceViMap; + 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 : g_reducedPatternExludedEdges) { + for (const size_t ei : global.g_reducedPatternExludedEdges) { colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0); } const std::string label = reducedPattern.getLabel(); @@ -427,11 +437,12 @@ void ReducedModelOptimizer::computeMaps( std::vector nodeColorsReducedToFull_full(fullPattern.VN(), glm::vec3(0, 0, 0)); for (size_t vi = 0; vi < reducedPattern.VN(); vi++) { - if (reducedToFullInterfaceViMap.contains(vi)) { + if (global.reducedToFullInterfaceViMap.contains(vi)) { auto color = polyscope::getNextUniqueColor(); nodeColorsReducedToFull_reduced[vi] = color; - nodeColorsReducedToFull_full[reducedToFullInterfaceViMap[vi]] = color; + nodeColorsReducedToFull_full[global.reducedToFullInterfaceViMap[vi]] = + color; } } polyscope::getCurveNetwork(reducedPattern.getLabel()) @@ -482,8 +493,8 @@ void ReducedModelOptimizer::initializePatterns( FlatPattern copyReducedPattern; copyFullPattern.copy(fullPattern); copyReducedPattern.copy(reducedPattern); - g_optimizeInnerHexagonSize = copyReducedPattern.EN() == 2; - if (g_optimizeInnerHexagonSize) { + global.g_optimizeInnerHexagonSize = copyReducedPattern.EN() == 2; + if (global.g_optimizeInnerHexagonSize) { const double h = copyReducedPattern.getBaseTriangleHeight(); double baseTriangle_bottomEdgeSize = 2 * h / tan(vcg::math::ToRad(60.0)); VectorType baseTriangle_leftBottomNode(-baseTriangle_bottomEdgeSize / 2, -h, @@ -497,13 +508,13 @@ void ReducedModelOptimizer::initializePatterns( vcg::RotationMatrix(rotationAxis, vcg::math::ToRad(rotationCounter * 60.0)) * baseTriangle_leftBottomNode; - g_innerHexagonVectors[rotationCounter] = rotatedVector; + global.g_innerHexagonVectors[rotationCounter] = rotatedVector; } const double innerHexagonInitialPos_x = - copyReducedPattern.vert[0].cP()[0] / g_innerHexagonVectors[0][0]; + copyReducedPattern.vert[0].cP()[0] / global.g_innerHexagonVectors[0][0]; const double innerHexagonInitialPos_y = - copyReducedPattern.vert[0].cP()[1] / g_innerHexagonVectors[0][1]; - g_innerHexagonInitialPos = innerHexagonInitialPos_x; + copyReducedPattern.vert[0].cP()[1] / global.g_innerHexagonVectors[0][1]; + global.g_innerHexagonInitialPos = innerHexagonInitialPos_x; } computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges); createSimulationMeshes(copyFullPattern, copyReducedPattern); @@ -513,9 +524,9 @@ void ReducedModelOptimizer::initializePatterns( void ReducedModelOptimizer::initializeOptimizationParameters( const std::shared_ptr &mesh) { const int numberOfOptimizationParameters = 3; - g_initialParameters.resize(g_optimizeInnerHexagonSize - ? numberOfOptimizationParameters + 1 - : numberOfOptimizationParameters); + global.g_initialParameters.resize(global.g_optimizeInnerHexagonSize + ? numberOfOptimizationParameters + 1 + : numberOfOptimizationParameters); // Save save the beam stiffnesses // for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) { // Element &e = pReducedModelElementalMesh->elements[ei]; @@ -529,9 +540,9 @@ void ReducedModelOptimizer::initializeOptimizationParameters( const double initialB = std::sqrt(mesh->elements[0].A); const double initialRatio = 1; ; - g_initialParameters(0) = initialB; - g_initialParameters(1) = initialRatio; - g_initialParameters(2) = mesh->elements[0].material.youngsModulus; + global.g_initialParameters(0) = initialB; + global.g_initialParameters(1) = initialRatio; + global.g_initialParameters(2) = mesh->elements[0].material.youngsModulus; // g_initialParameters = // m_pReducedPatternSimulationMesh->elements[0].properties.E; // for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) { @@ -605,9 +616,9 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization( const Settings &settings, double (*pObjectiveFunction)(long, const double *)) { - gObjectiveValueHistory.clear(); + global.gObjectiveValueHistory.clear(); - const size_t n = g_initialParameters.rows(); + const size_t n = global.g_initialParameters.rows(); assert(n == 3); // g_optimizeInnerHexagonSize ? 5: 4; // const size_t npt = (n + 1) * (n + 2) / 2; @@ -619,8 +630,8 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization( const size_t initialGuess = 1; std::vector x(n, initialGuess); - if (g_optimizeInnerHexagonSize) { - x[n - 1] = g_innerHexagonInitialPos; + if (global.g_optimizeInnerHexagonSize) { + x[n - 1] = global.g_innerHexagonInitialPos; } /*if (!initialGuess.empty()) { x = g_optimizationInitialGuess; @@ -655,7 +666,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization( xMax(i) = settings.xRanges[i].max; } - numberOfFunctionCalls = 0; + global.numberOfFunctionCalls = 0; double (*objF)(double, double, double) = &objective; auto start = std::chrono::system_clock::now(); dlib::function_evaluation result = dlib::find_min_global( @@ -663,11 +674,11 @@ 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{numOfSimulationCrashes, minX, minY}; + Results results{global.numOfSimulationCrashes, global.minX, global.minY}; std::cout << "Finished optimizing." << endl; // std::cout << "Solution x:" << endl; // std::cout << result.x << endl; - std::cout << "Objective value:" << minY << endl; + std::cout << "Objective value:" << global.minY << endl; // std::cout << result.y << endl; // std::cout << minY << endl; // std::cout << "Time(sec):" << elapsed.count() << std::endl; @@ -874,87 +885,89 @@ ReducedModelOptimizer::createScenarios( return scenarios; } -void ReducedModelOptimizer::runBeamOptimization() { - // load beams - VCGEdgeMesh fullBeam; - fullBeam.loadPly("/home/iason/Models/simple_beam_model_10elem_1m.ply"); - VCGEdgeMesh reducedBeam; - reducedBeam.loadPly("/home/iason/Models/simple_beam_model_4elem_1m.ply"); - fullBeam.registerForDrawing(); - reducedBeam.registerForDrawing(); - // polyscope::show(); - // maps - std::unordered_map displacementReducedToFullMap; - displacementReducedToFullMap[reducedBeam.VN() / 2] = fullBeam.VN() / 2; - g_reducedToFullViMap = displacementReducedToFullMap; - std::unordered_map jobFullToReducedMap; - jobFullToReducedMap[0] = 0; - jobFullToReducedMap[fullBeam.VN() - 1] = reducedBeam.VN() - 1; +// void ReducedModelOptimizer::runBeamOptimization() { +// // load beams +// VCGEdgeMesh fullBeam; +// fullBeam.loadPly("/home/iason/Models/simple_beam_model_10elem_1m.ply"); +// VCGEdgeMesh reducedBeam; +// reducedBeam.loadPly("/home/iason/Models/simple_beam_model_4elem_1m.ply"); +// fullBeam.registerForDrawing(); +// reducedBeam.registerForDrawing(); +// // polyscope::show(); +// // maps +// std::unordered_map displacementReducedToFullMap; +// displacementReducedToFullMap[reducedBeam.VN() / 2] = fullBeam.VN() / 2; +// g_reducedToFullViMap = displacementReducedToFullMap; +// std::unordered_map jobFullToReducedMap; +// jobFullToReducedMap[0] = 0; +// jobFullToReducedMap[fullBeam.VN() - 1] = reducedBeam.VN() - 1; - // full model simuilation job - auto pFullPatternSimulationMesh = std::make_shared(fullBeam); - pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02}); - pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); - std::unordered_map> fixedVertices; - fixedVertices[0] = ::unordered_set({0, 1, 2, 3, 4, 5}); - std::unordered_map forces; - forces[fullBeam.VN() - 1] = Vector6d({0, 0, 10, 0, 0, 0}); - const std::string fullBeamSimulationJobLabel = "Pull_Z"; - std::shared_ptr pFullModelSimulationJob = - make_shared(SimulationJob(pFullPatternSimulationMesh, - fullBeamSimulationJobLabel, - fixedVertices, forces)); - FormFinder formFinder; - auto fullModelResults = formFinder.executeSimulation(pFullModelSimulationJob); +// // full model simuilation job +// auto pFullPatternSimulationMesh = +// std::make_shared(fullBeam); +// pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, +// 0.02}); pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); +// std::unordered_map> fixedVertices; +// fixedVertices[0] = ::unordered_set({0, 1, 2, 3, 4, 5}); +// std::unordered_map forces; +// forces[fullBeam.VN() - 1] = Vector6d({0, 0, 10, 0, 0, 0}); +// const std::string fullBeamSimulationJobLabel = "Pull_Z"; +// std::shared_ptr pFullModelSimulationJob = +// make_shared(SimulationJob(pFullPatternSimulationMesh, +// fullBeamSimulationJobLabel, +// fixedVertices, forces)); +// auto fullModelResults = +// formFinder.executeSimulation(pFullModelSimulationJob); - // Optimal reduced model displacements - const size_t numberOfSimulationScenarios = 1; - g_optimalReducedModelDisplacements.resize(numberOfSimulationScenarios); - g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1].resize( - reducedBeam.VN(), 3); - computeDesiredReducedModelDisplacements( - fullModelResults, displacementReducedToFullMap, - g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]); +// // Optimal reduced model displacements +// const size_t numberOfSimulationScenarios = 1; +// g_optimalReducedModelDisplacements.resize(numberOfSimulationScenarios); +// g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1].resize( +// reducedBeam.VN(), 3); +// computeDesiredReducedModelDisplacements( +// fullModelResults, displacementReducedToFullMap, +// g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]); - // reduced model simuilation job - auto reducedSimulationMesh = std::make_shared(reducedBeam); - reducedSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02}); - reducedSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); - g_reducedPatternSimulationJob.resize(numberOfSimulationScenarios); - SimulationJob reducedSimJob; - computeReducedModelSimulationJob(*pFullModelSimulationJob, - jobFullToReducedMap, reducedSimJob); - reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] = - reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] * 0.1; - g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1] = - make_shared( - reducedSimulationMesh, fullBeamSimulationJobLabel, - reducedSimJob.constrainedVertices, reducedSimJob.nodalExternalForces); - initializeOptimizationParameters(reducedSimulationMesh); +// // reduced model simuilation job +// auto reducedSimulationMesh = std::make_shared(reducedBeam); +// reducedSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02}); +// reducedSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); +// g_reducedPatternSimulationJob.resize(numberOfSimulationScenarios); +// SimulationJob reducedSimJob; +// computeReducedModelSimulationJob(*pFullModelSimulationJob, +// jobFullToReducedMap, reducedSimJob); +// reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] = +// reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] * 0.1; +// g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1] = +// make_shared( +// reducedSimulationMesh, fullBeamSimulationJobLabel, +// reducedSimJob.constrainedVertices, +// reducedSimJob.nodalExternalForces); +// initializeOptimizationParameters(reducedSimulationMesh); - // const std::string simulationJobsPath = "SimulationJobs"; - // std::filesystem::create_directory(simulationJobsPath); - // g_reducedPatternSimulationJob[0].save(simulationJobsPath); - // g_reducedPatternSimulationJob[0].load( - // std::filesystem::path(simulationJobsPath) - // .append(g_reducedPatternSimulationJob[0].mesh->getLabel() + - // "_simScenario.json")); +// // const std::string simulationJobsPath = "SimulationJobs"; +// // std::filesystem::create_directory(simulationJobsPath); +// // g_reducedPatternSimulationJob[0].save(simulationJobsPath); +// // g_reducedPatternSimulationJob[0].load( +// // std::filesystem::path(simulationJobsPath) +// // .append(g_reducedPatternSimulationJob[0].mesh->getLabel() + +// // "_simScenario.json")); - runOptimization({}, &objective); +// runOptimization({}, &objective); - fullModelResults.registerForDrawing(); - SimulationResults reducedModelResults = simulator.executeSimulation( - g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1]); - double error = computeError( - reducedModelResults, - g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]); - reducedModelResults.registerForDrawing(); - std::cout << "Error between beams:" << error << endl; - // registerWorldAxes(); - polyscope::show(); - fullModelResults.unregister(); - reducedModelResults.unregister(); -} +// fullModelResults.registerForDrawing(); +// SimulationResults reducedModelResults = simulator.executeSimulation( +// g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1]); +// double error = computeError( +// reducedModelResults, +// g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]); +// reducedModelResults.registerForDrawing(); +// std::cout << "Error between beams:" << error << endl; +// // registerWorldAxes(); +// polyscope::show(); +// fullModelResults.unregister(); +// reducedModelResults.unregister(); +//} void ReducedModelOptimizer::visualizeResults( const std::vector> @@ -971,12 +984,12 @@ void ReducedModelOptimizer::visualizeResults( fullModelResults.registerForDrawing(); fullModelResults.saveDeformedModel(); const std::shared_ptr &pReducedPatternSimulationJob = - g_reducedPatternSimulationJob[simulationScenarioIndex]; + global.g_reducedPatternSimulationJob[simulationScenarioIndex]; SimulationResults reducedModelResults = simulator.executeSimulation(pReducedPatternSimulationJob); double error = computeError( reducedModelResults, - g_optimalReducedModelDisplacements[simulationScenarioIndex]); + global.g_optimalReducedModelDisplacements[simulationScenarioIndex]); std::cout << "Error of simulation scenario " << simulationScenarioStrings[simulationScenarioIndex] << " is " << error << std::endl; @@ -1001,9 +1014,9 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( const Settings &xRanges, const std::vector &simulationScenarios) { - g_simulationScenarioIndices = simulationScenarios; - if (g_simulationScenarioIndices.empty()) { - g_simulationScenarioIndices = { + global.g_simulationScenarioIndices = simulationScenarios; + if (global.g_simulationScenarioIndices.empty()) { + global.g_simulationScenarioIndices = { SimulationScenario::Axial, SimulationScenario::Shear, SimulationScenario::Bending, SimulationScenario::Dome, SimulationScenario::Saddle}; @@ -1011,32 +1024,32 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( std::vector> simulationJobs = createScenarios(m_pFullPatternSimulationMesh); - g_optimalReducedModelDisplacements.resize(6); - g_reducedPatternSimulationJob.resize(6); - g_firstRoundIterationIndex = 0; - minY = std::numeric_limits::max(); - numOfSimulationCrashes = 0; + global.g_optimalReducedModelDisplacements.resize(6); + global.g_reducedPatternSimulationJob.resize(6); + global.g_firstRoundIterationIndex = 0; + global.minY = std::numeric_limits::max(); + global.numOfSimulationCrashes = 0; // polyscope::removeAllStructures(); FormFinder::Settings settings; // settings.shouldDraw = true; - for (int simulationScenarioIndex : g_simulationScenarioIndices) { + for (int simulationScenarioIndex : global.g_simulationScenarioIndices) { const std::shared_ptr &pFullPatternSimulationJob = simulationJobs[simulationScenarioIndex]; SimulationResults fullModelResults = simulator.executeSimulation(pFullPatternSimulationJob, settings); - g_optimalReducedModelDisplacements[simulationScenarioIndex].resize( + global.g_optimalReducedModelDisplacements[simulationScenarioIndex].resize( m_pReducedPatternSimulationMesh->VN(), 3); computeDesiredReducedModelDisplacements( - fullModelResults, g_reducedToFullViMap, - g_optimalReducedModelDisplacements[simulationScenarioIndex]); + fullModelResults, global.g_reducedToFullViMap, + global.g_optimalReducedModelDisplacements[simulationScenarioIndex]); SimulationJob reducedPatternSimulationJob; reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; computeReducedModelSimulationJob(*pFullPatternSimulationJob, m_fullToReducedInterfaceViMap, reducedPatternSimulationJob); - g_reducedPatternSimulationJob[simulationScenarioIndex] = - std::make_shared(reducedPatternSimulationJob); + global.g_reducedPatternSimulationJob[simulationScenarioIndex] = + std::make_shared(reducedPatternSimulationJob); } Results optResults = runOptimization(xRanges, &objective); updateMesh(optResults.x.size(), optResults.x.data()); diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 2c6dd26..81a53e0 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -106,6 +106,7 @@ private: computeError(const SimulationResults &reducedPatternResults, const Eigen::MatrixX3d &optimalReducedPatternDisplacements); static double objective(long n, const double *x); + FormFinder simulator; }; #endif // REDUCEDMODELOPTIMIZER_HPP