From 7a0124155cf51410fc82934d112b9133f6caab18 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Tue, 19 Jul 2022 14:59:05 +0300 Subject: [PATCH] Uses chronos simulation model. Refactoring --- reducedmodelevaluator.cpp | 77 ++++++++++++++++--------------- reducedmodeloptimizer.cpp | 25 ++++++---- reducedmodeloptimizer.hpp | 17 +++---- reducedmodeloptimizer_structs.hpp | 45 +++++++++--------- 4 files changed, 90 insertions(+), 74 deletions(-) diff --git a/reducedmodelevaluator.cpp b/reducedmodelevaluator.cpp index 726fafc..51a83ef 100644 --- a/reducedmodelevaluator.cpp +++ b/reducedmodelevaluator.cpp @@ -331,19 +331,19 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( pTiledFullPattern->setLabel("Tiled_full_patterns"); // pTiledFullPattern->registerForDrawing(); //Tile reduced pattern into surface - // PatternGeometry reducedPattern; ReducedModel reducedModel; - // reducedModel.registerForDrawing(); - // polyscope::show(); reducedModel.deleteDanglingVertices(); - // reducedPattern.interfaceNodeIndex = 1; - // assert(reducedPattern.interfaceNodeIndex == 1); + std::unordered_map + optimalXVariables_set(optimizationResult.optimalXNameValuePairs.begin(), + optimizationResult.optimalXNameValuePairs.end()); + reducedModel.updateBaseTriangleGeometry_R(optimalXVariables_set.at("R")); + reducedModel.updateBaseTriangleGeometry_theta(optimalXVariables_set.at("Theta")); + // const auto reducedPatternBaseTriangle = reducedModel.computeBaseTriangle(); + // ReducedModelOptimization::Results::applyOptimizationResults_reducedModel_nonFanned( + // optimizationResult, reducedPatternBaseTriangle, reducedModel); + std::vector reducedPatterns(1); reducedPatterns[0].copy(reducedModel); - const auto reducedPatternBaseTriangle = reducedModel.computeBaseTriangle(); - ReducedModelOptimization::Results::applyOptimizationResults_reducedModel_nonFanned( - optimizationResult, reducedPatternBaseTriangle, reducedPatterns[0]); - std::vector> perPatternIndexTiledReducedPatternEdgeIndices; std::vector tileIntoEdgeToTiledReducedVi; std::shared_ptr pTiledReducedPattern @@ -354,8 +354,9 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( tileIntoEdgeToTiledReducedVi, perPatternIndexTiledReducedPatternEdgeIndices); pTiledReducedPattern->setLabel("Tiled_reduced_patterns"); - // pTiledReducedPattern->registerForDrawing(); - // polyscope::show(); +#ifdef POLYSCOPE_DEFINED + pTiledReducedPattern->registerForDrawing(); +#endif std::unordered_map fullToReducedViMap; //of only the common vertices @@ -379,9 +380,10 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( std::shared_ptr pTiledFullPattern_simulationMesh; pTiledFullPattern_simulationMesh = std::make_shared(*pTiledFullPattern); //NOTE: Those should be derived from the optimization results instead of hardcoding them - constexpr double beamWidth = 0.001; - constexpr double beamHeight = 0.001; - pTiledFullPattern_simulationMesh->setBeamCrossSection(CrossSectionType{beamWidth, beamHeight}); + // const double beamWidth = optimizationResult.settings.patternBeamDimensions.getWidth(); + // const double beamHeight = optimizationResult.settings.patternBeamDimensions.getHeight(); + pTiledFullPattern_simulationMesh->setBeamCrossSection( + optimizationResult.settings.patternBeamDimensions); if (optimizationResult.fullPatternYoungsModulus == 0) { std::cerr << "Full pattern's young modulus not found." << std::endl; std::terminate(); @@ -419,7 +421,7 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( #endif constexpr bool shouldRerunFullPatternSimulation = false; enum PatternSimulationModelTag { DRM, Chronos }; - const PatternSimulationModelTag simulationModelTag{DRM}; + const PatternSimulationModelTag simulationModelTag{Chronos}; const std::string simulationModelLabel = [&]() { switch (simulationModelTag) { case DRM: @@ -431,7 +433,7 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( // for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) { std::for_each( //#ifndef POLYSCOPE_DEFINED - std::execution::par_unseq, + // std::execution::par_unseq, //#endif scenariosTestSetLabels.begin(), scenariosTestSetLabels.end(), @@ -460,9 +462,9 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( pJob_tiledFullPattern = std::make_shared(SimulationJob()); pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh; pJob_tiledReducedPattern->remap(reducedToFullViMap, *pJob_tiledFullPattern); - // pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern->getLabel()); - // pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern->getLabel()); - // polyscope::show(); + // pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern->getLabel()); + // pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern->getLabel()); + // polyscope::show(); const std::filesystem::path surfaceFolderPath = std::filesystem::path(fullPatternTessellatedResultsDirectoryPath) .append(simulationModelLabel + "_" + pTileIntoSurface->getLabel()); @@ -501,12 +503,12 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( scenarioDirectoryPath) .append(patternLabel) .append("SimulationJob"); - SimulationResults simulationResults_tiledPattern; + SimulationResults simulationResults_tilledPattern; if (std::filesystem::exists(tilledPatternResultsFolderPath)) { //Load full pattern results assert(std::filesystem::exists(fullPatternJobFolderPath)); - simulationResults_tiledPattern.load(tilledPatternResultsFolderPath, - fullPatternJobFolderPath); + simulationResults_tilledPattern.load(tilledPatternResultsFolderPath, + fullPatternJobFolderPath); //#ifdef POLYSCOPE_DEFINED // std::array resultsColor({28.0, 99.0, 227.0}); // simulationResults_tiledFullPattern_drm.registerForDrawing(resultsColor); @@ -532,7 +534,7 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( // simulationResults_tiledFullPattern_drm.unregister(); //#endif - simulationResults_tiledPattern.converged = true; + simulationResults_tilledPattern.converged = true; } else { std::cout << "Tilled pattern simulation results not found in:" << tilledPatternResultsFolderPath << std::endl; @@ -541,20 +543,20 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( switch (simulationModelTag) { case DRM: { DRMSimulationModel drmSimulationModel; - simulationResults_tiledPattern = drmSimulationModel - .executeSimulation(pJob_tiledFullPattern, - drmSimulationSettings); + simulationResults_tilledPattern = drmSimulationModel + .executeSimulation(pJob_tiledFullPattern, + drmSimulationSettings); break; } case Chronos: { ChronosEulerSimulationModel chronosSimulationModel; - simulationResults_tiledPattern = chronosSimulationModel.executeSimulation( + simulationResults_tilledPattern = chronosSimulationModel.executeSimulation( pJob_tiledFullPattern); break; } } } - if (!simulationResults_tiledPattern.converged) { + if (!simulationResults_tilledPattern.converged) { std::cerr << "Full pattern simulation failed." << std::endl; std::cerr << "Not saving results" << std::endl; // continue; @@ -563,22 +565,25 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( std::filesystem::create_directories(tilledPatternResultsFolderPath); const std::filesystem::path drmResultsOutputPath = std::filesystem::path(scenarioDirectoryPath).append(patternLabel); - simulationResults_tiledPattern.save(drmResultsOutputPath); + simulationResults_tilledPattern.save(drmResultsOutputPath); LinearSimulationModel linearSimulationModel; - SimulationResults simulationResults_tiledReducedPattern + SimulationResults simulationResults_tiledReducedModel = linearSimulationModel.executeSimulation(pJob_tiledReducedPattern); - // simulationResults_tiledReducedPattern.registerForDrawing(); - // simulationResults_tiledFullPattern_drm.registerForDrawing(); - // polyscope::show(); +#ifdef POLYSCOPE_DEFINED + simulationResults_tiledReducedModel.registerForDrawing(); + simulationResults_tilledPattern.registerForDrawing(); + polyscope::show(); +#endif //compute the full2reduced distance - const double distance_patternToReduced = simulationResults_tiledPattern.computeDistance( - simulationResults_tiledReducedPattern, fullToReducedViMap); + const double distance_patternToReduced + = simulationResults_tilledPattern + .computeDistance(simulationResults_tiledReducedModel, fullToReducedViMap); double distance_patternSumOfAllVerts = 0; for (std::pair fullToReducedPair : fullToReducedViMap) { - distance_patternSumOfAllVerts += simulationResults_tiledPattern + distance_patternSumOfAllVerts += simulationResults_tilledPattern .displacements[fullToReducedPair.first] .getTranslation() .norm(); diff --git a/reducedmodeloptimizer.cpp b/reducedmodeloptimizer.cpp index 9d63351..19402d2 100644 --- a/reducedmodeloptimizer.cpp +++ b/reducedmodeloptimizer.cpp @@ -362,8 +362,9 @@ double ReducedModelOptimizer::objective(const std::vector &x, } void ReducedModelOptimizer::createSimulationMeshes( - PatternGeometry &fullModel, + PatternGeometry &pattern, PatternGeometry &reducedModel, + const RectangularBeamDimensions &beamDimensions, std::shared_ptr &pFullPatternSimulationMesh, std::shared_ptr &pReducedPatternSimulationMesh) { @@ -373,23 +374,25 @@ void ReducedModelOptimizer::createSimulationMeshes( } // Full pattern - pFullPatternSimulationMesh = std::make_shared(fullModel); - pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{beamWidth, beamHeight}); + pFullPatternSimulationMesh = std::make_shared(pattern); + pFullPatternSimulationMesh->setBeamCrossSection(beamDimensions); pFullPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus); pFullPatternSimulationMesh->reset(); // Reduced pattern pReducedPatternSimulationMesh = std::make_shared(reducedModel); - pReducedPatternSimulationMesh->setBeamCrossSection(CrossSectionType{beamWidth, beamHeight}); + pReducedPatternSimulationMesh->setBeamCrossSection(beamDimensions); pReducedPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus); pReducedPatternSimulationMesh->reset(); } void ReducedModelOptimizer::createSimulationMeshes(PatternGeometry &fullModel, - PatternGeometry &reducedModel) + PatternGeometry &reducedModel, + const RectangularBeamDimensions &beamDimensions) { ReducedModelOptimizer::createSimulationMeshes(fullModel, reducedModel, + beamDimensions, m_pFullPatternSimulationMesh, m_pReducedModelSimulationMesh); } @@ -558,9 +561,11 @@ ReducedModelOptimizer::ReducedModelOptimizer() void ReducedModelOptimizer::initializePatterns( PatternGeometry &fullPattern, ReducedModel &reducedModel, - const std::array &optimizationParameters) + const ReducedModelOptimization::Settings &optimizationSettings) { assert(fullPattern.VN() == reducedModel.VN() && fullPattern.EN() >= reducedModel.EN()); + const std::array &optimizationParameters + = optimizationSettings.variablesRanges; #ifdef POLYSCOPE_DEFINED polyscope::removeAllStructures(); @@ -580,7 +585,9 @@ void ReducedModelOptimizer::initializePatterns( optimizationState.fullPatternOppositeInterfaceViPairs = m_fullPatternOppositeInterfaceViPairs; g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs = m_fullPatternOppositeInterfaceViPairs; - createSimulationMeshes(copyFullPattern, copyReducedModel); + createSimulationMeshes(copyFullPattern, + copyReducedModel, + optimizationSettings.patternBeamDimensions); initializeUpdateReducedPatternFunctions(); initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters); } @@ -596,7 +603,7 @@ void ReducedModelOptimizer::optimize(ConstPatternGeometry &fullPattern, auto start = std::chrono::system_clock::now(); optimizationResults.label = fullPattern.getLabel(); optimizationResults.baseTriangleFullPattern.copy(fullPattern); - initializePatterns(fullPattern, reducedModel, optimizationSettings.variablesRanges); + initializePatterns(fullPattern, reducedModel, optimizationSettings); optimize(optimizationSettings, optimizationResults); optimizationResults.settings = optimizationSettings; //NOTE: currently I set the max force base magn in optimize thus this must be called after the optimize function @@ -668,7 +675,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() const double beamHeight = beamWidth; for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; - e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight)); + e.setDimensions(CrossSectionType(beamWidth, beamHeight)); } }; diff --git a/reducedmodeloptimizer.hpp b/reducedmodeloptimizer.hpp index a307eb4..7de94c9 100644 --- a/reducedmodeloptimizer.hpp +++ b/reducedmodeloptimizer.hpp @@ -90,8 +90,8 @@ private: std::vector scenarioIsSymmetrical; int fullPatternNumberOfEdges; constexpr static double youngsModulus{1 * 1e9}; - constexpr static double beamWidth{0.001}; - constexpr static double beamHeight{0.001}; + constexpr static double defaultBeamWidth{0.002}; + constexpr static double defaultBeamHeight{0.002}; std::string fullPatternLabel; // ReducedModelOptimization::Settings optimizationSettings; @@ -147,8 +147,9 @@ public: const std::unordered_map &fullPatternOppositeInterfaceViMap); static void createSimulationMeshes( - PatternGeometry &fullModel, + PatternGeometry &pattern, PatternGeometry &reducedModel, + const RectangularBeamDimensions &beamDimensions, std::shared_ptr &pFullPatternSimulationMesh, std::shared_ptr &pReducedPatternSimulationMesh); void computeMaps(const std::unordered_map> &slotToNode, @@ -288,10 +289,8 @@ private: ReducedModelOptimization::S})); void initializePatterns(PatternGeometry &fullPattern, - ReducedModel &reducedPattern, - const std::array - &optimizationParameters); + ReducedModel &reducedModel, + const ReducedModelOptimization::Settings &optimizationSettings); static void computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, const std::unordered_map &displacementsReducedToFullMap, @@ -299,7 +298,9 @@ private: void runOptimization(const ReducedModelOptimization::Settings &settings, ReducedModelOptimization::Results &results); void computeMaps(PatternGeometry &fullModel, ReducedModel &reducedModel); - void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); + void createSimulationMeshes(PatternGeometry &fullModel, + PatternGeometry &reducedModel, + const RectangularBeamDimensions &beamDimensions); void initializeOptimizationParameters( const std::shared_ptr &mesh, const std::array baseScenarioMaxMagnitudes{ - 0.590241 / 6, 0.588372 / 6, 0.368304 / 2, 0.05, 1.18 / 4, 0}; //final b,h= 0.001 + // std::array baseScenarioMaxMagnitudes{ + // 0.590241 / 6, 0.588372 / 6, 0.368304 / 2, 0.05, 1.18 / 4, 0}; //final b,h= 0.001 - // std::array baseScenarioMaxMagnitudes{0.590241 / 3, - // 0.588372 / 3, - // 0.368304, - // 0.1, - // 1.18 / 2, - // 0}; //final b,h= 0.002 + std::array baseScenarioMaxMagnitudes{ + 0.590241 / 3, 0.588372 / 3, 0.368304, 0.1, 1.18 / 2, 0}; //final b,h= 0.002 // std::array baseScenarioMaxMagnitudes{0, 0, 0, 0.1, 0}; // std::array @@ -233,6 +229,7 @@ struct Settings // double angularDistanceEpsilon{vcg::math::ToRad(2.0)}; double angularDistanceEpsilon{vcg::math::ToRad(0.0)}; double targetBaseTriangleSize{0.03}; + RectangularBeamDimensions patternBeamDimensions{0.002, 0.002}; std::filesystem::path intermediateResultsDirectoryPath; struct ObjectiveWeights { @@ -310,6 +307,9 @@ struct Settings json[GET_VARIABLE_NAME(angularDistanceEpsilon)] = vcg::math::ToDeg(angularDistanceEpsilon); json[GET_VARIABLE_NAME(targetBaseTriangleSize)] = targetBaseTriangleSize; json[GET_VARIABLE_NAME(baseScenarioMaxMagnitudes)] = baseScenarioMaxMagnitudes; + nlohmann::json json_dimensions; + patternBeamDimensions.to_json(json_dimensions, patternBeamDimensions); + json.update(json_dimensions); #ifdef USE_ENSMALLEN #ifdef USE_PSO @@ -412,6 +412,8 @@ struct Settings pso.numberOfParticles = static_cast(json[GET_VARIABLE_NAME(pso.numberOfParticles)]); } + patternBeamDimensions.from_json(json, patternBeamDimensions); + // perBaseScenarioObjectiveWeights = json.at(JsonKeys::ObjectiveWeights); // objectiveWeights.translational = json.at(JsonKeys::ObjectiveWeights); // objectiveWeights.rotational = 2 - objectiveWeights.translational; @@ -555,6 +557,7 @@ struct Settings std::vector> reducedPatternSimulationJobs; double fullPatternYoungsModulus{0}; + //Full pattern PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern vcg::Triangle3 baseTriangle; // std::string notes; @@ -808,9 +811,9 @@ struct Settings JsonKeys::FullPatternLabel); if (!baseTriangleFullPattern.load( std::filesystem::path(loadFromPath).append(fullPatternLabel + ".ply").string())) { - !baseTriangleFullPattern.load(std::filesystem::path(loadFromPath) - .append(loadFromPath.stem().string() + ".ply") - .string()); + baseTriangleFullPattern.load(std::filesystem::path(loadFromPath) + .append(loadFromPath.stem().string() + ".ply") + .string()); } std::vector baseTriangleVertices = json_optimizationResults.at( @@ -872,7 +875,7 @@ struct Settings static void applyOptimizationResults_reducedModel_nonFanned( const ReducedModelOptimization::Results &reducedPattern_optimizationResults, const vcg::Triangle3 &patternBaseTriangle, - MeshType &reducedPattern) + MeshType &reducedModel) { std::unordered_map optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(), @@ -884,13 +887,13 @@ struct Settings applyOptimizationResults_reducedModel_nonFanned(optimalXVariables.at("HexSize"), optimalXVariables.at("HexAngle"), patternBaseTriangle, - reducedPattern); + reducedModel); return; } applyOptimizationResults_reducedModel_nonFanned(optimalXVariables.at("R"), optimalXVariables.at("Theta"), patternBaseTriangle, - reducedPattern); + reducedModel); } template @@ -898,7 +901,7 @@ struct Settings const double &hexSize, const double &hexAngle, const vcg::Triangle3 &patternBaseTriangle, - MeshType &reducedPattern) + MeshType &reducedModel) { //Set optimal geometrical params of the reduced pattern // CoordType center_barycentric(1, 0, 0); @@ -907,13 +910,13 @@ struct Settings // (center_barycentric * (1 - hexSize) + interfaceEdgeMiddle_barycentric)); CoordType movableVertex_barycentric(1 - hexSize, hexSize / 2, hexSize / 2); - reducedPattern.vert[0].P() = patternBaseTriangle.cP(0) * movableVertex_barycentric[0] - + patternBaseTriangle.cP(1) * movableVertex_barycentric[1] - + patternBaseTriangle.cP(2) * movableVertex_barycentric[2]; + reducedModel.vert[0].P() = patternBaseTriangle.cP(0) * movableVertex_barycentric[0] + + patternBaseTriangle.cP(1) * movableVertex_barycentric[1] + + patternBaseTriangle.cP(2) * movableVertex_barycentric[2]; if (hexAngle != 0) { - reducedPattern.vert[0].P() = vcg::RotationMatrix(CoordType{0, 0, 1}, - vcg::math::ToRad(hexAngle)) - * reducedPattern.vert[0].cP(); + reducedModel.vert[0].P() = vcg::RotationMatrix(CoordType{0, 0, 1}, + vcg::math::ToRad(hexAngle)) + * reducedModel.vert[0].cP(); } // for (int rotationCounter = 0; // rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) {