commit 81a60c686c60d500f8d27c26ded3cb1eef6dd644 Author: iasonmanolas Date: Mon May 31 13:40:00 2021 +0300 Initial commit diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100755 index 0000000..1949bfa --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,61 @@ +cmake_minimum_required(VERSION 2.8) +project(PatternTillingReducedModel) + +#Add the project cmake scripts to the module path +list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake) + +#Download external dependencies NOTE: If the user has one of these libs it shouldn't be downloaded again. +include(${CMAKE_MODULE_PATH}/DownloadProject.cmake) +if (CMAKE_VERSION VERSION_LESS 3.2) + set(UPDATE_DISCONNECTED_IF_AVAILABLE "") +else() + set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1") +endif() + +set(EXTERNAL_DEPS_DIR "/home/iason/Coding/build/external dependencies") +if(NOT EXISTS ${EXTERNAL_DEPS_DIR}) + set(EXTERNAL_DEPS_DIR ${CMAKE_CURRENT_BINARY_DIR}) +endif() +##Create directory for the external libraries +file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR}) + +add_compile_definitions(ENABLE_OPENMP) + +set(USE_POLYSCOPE TRUE) +if(${USE_POLYSCOPE}) +##Polyscope +download_project(PROJ POLYSCOPE + GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git + GIT_TAG master + PREFIX ${EXTERNAL_DEPS_DIR} + ${UPDATE_DISCONNECTED_IF_AVAILABLE} +) +add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR}) +add_compile_definitions(POLYSCOPE_DEFINED) +endif() + +set(MYSOURCES_SOURCE_DIR "/home/iason/Coding/Libraries/MySources") +if (NOT EXISTS ${MYSOURCES_SOURCE_DIR}) +##MySources +download_project(PROJ MYSOURCES + GIT_REPOSITORY https://gitea-s2i2s.isti.cnr.it/manolas/MySources.git + GIT_TAG master + PREFIX ${EXTERNAL_DEPS_DIR} + ${UPDATE_DISCONNECTED_IF_AVAILABLE} + ) +endif() +add_subdirectory(${MYSOURCES_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}/MySourcesBinDir) + +#Add the project sources +file(GLOB SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/*.hpp ${POLYSCOPE_SOURCE_DIR}/deps/imgui/imgui/misc/cpp/imgui_stdlib.cpp) + + +add_executable(${PROJECT_NAME} ${SOURCES} ) +target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_20) + +target_include_directories(${PROJECT_NAME} + PRIVATE ${MYSOURCES_SOURCE_DIR} + ) +target_link_libraries(${PROJECT_NAME} polyscope MySources ) + + diff --git a/src/hexagonremesher.hpp b/src/hexagonremesher.hpp new file mode 100644 index 0000000..d4a56c5 --- /dev/null +++ b/src/hexagonremesher.hpp @@ -0,0 +1,128 @@ +#ifndef HEXAGONREMESHER_HPP +#define HEXAGONREMESHER_HPP +#include "polymesh.hpp" +#include "vcgtrimesh.hpp" +#include +#include +#include +#include +#include +#include + +namespace PolygonalRemeshing { +//std::shared_ptr +std::shared_ptr remeshWithPolygons(VCGTriMesh &surface) +{ + // const std::string tileIntoFilePath( + // "/home/iason/Coding/build/PatternTillingReducedModel/RelWithDebInfo/plane.ply"); + // const std::string tileIntoFilePath("/home/iason/Coding/build/PatternTillingReducedModel/" + // "RelWithDebInfo/hexagon.ply"); + // const std::string tileIntoFilePath("/home/iason/Coding/build/PatternTillingReducedModel/" + // "RelWithDebInfo/instantMeshes_plane_34.ply"); + // const std::string tileIntoFilePath("/home/iason/Coding/build/PatternTillingReducedModel/" + // "RelWithDebInfo/instantMeshes_plane.ply"); + // assert(std::filesystem::exists(tileIntoFilePath)); + // assert(tileInto_polyMesh.load(tileIntoFilePath)); + + ///Load triangle mesh + // VCGTriMesh tileInto_triMesh; + // assert(std::filesystem::exists(tileIntoFilePath)); + // tileInto_triMesh.load(tileIntoFilePath); + + ///Sample the triangle mesh + // std::vector pointVec; + // double radius = 1; + // int sampleNum = 50; + // // // vcg::tri::PoissonSampling(tileInto_triMesh, pointVec, sampleNum, radius); + // // // vcg::tri::TrivialSampler ps; + // // // vcg::tri::SurfaceSampling::VertexUniform(tileInto_triMesh, ps, sampleNum); + // // // pointVec = ps.SampleVec(); + // vcg::tri::MontecarloSampling(tileInto_triMesh, pointVec, sampleNum); + + // vcg::tri::VoronoiProcessingParameter vpp; + // vcg::tri::VoronoiProcessing::PreprocessForVoronoi(tileInto_triMesh, radius, vpp); + // std::vector seedVec; + // vcg::tri::VoronoiProcessing::SeedToVertexConversion(tileInto_triMesh, + // pointVec, + // seedVec); + + // vcg::tri::EuclideanDistance df; + // vpp.geodesicRelaxFlag = false; + // const int iterNum = 100; + // int actualIter = vcg::tri::VoronoiProcessing::VoronoiRelaxing(tileInto_triMesh, + // seedVec, + // iterNum, + // df, + // vpp); + // std::cout << "Voronoi relaxation iterations performed:" << actualIter << std::endl; + + // bool saveWasSuccessful + // = 0 + // == vcg::tri::io::ExporterOBJ::Save(tileInto_triMesh, + // "./tileInto_triMesh.obj", + // vcg::tri::io::Mask::IOM_ALL); + + // VCGTriMesh delaunay; + // vcg::tri::VoronoiProcessing::ConvertDelaunayTriangulationToMesh(tileInto_triMesh, + // delaunay, + // seedVec); + // delaunay.setLabel("Delaunay"); + // delaunay.save(delaunay.getLabel() + ".ply"); + // delaunay.registerForDrawing(); + + // ////Load surface + VCGPolyMesh tileInto_polyMesh, dual; + // vcg::tri::PolygonSupport::ImportFromTriMesh(tileInto_polyMesh, + // delaunay); + vcg::tri::PolygonSupport::ImportFromTriMesh(tileInto_polyMesh, surface); + vcg::tri::DualMeshing::MakeDual(tileInto_polyMesh, dual, false); + + // for(size_t i=0; i< Edges.size(); ++i) + // { + // std::vector fpVec; + // std::vector eiVec; + // face::EFStarFF(Edges[i].f,Edges[i].z,fpVec,eiVec); + // for(size_t j=0;jFEp(eiVec[j])=&(m.edge[i]); + // } + dual.setLabel(surface.getLabel() + "_polyDual"); + + bool saveWasSuccessful + = 0 + == vcg::tri::io::ExporterOBJ::Save(dual, + (dual.getLabel() + ".obj").c_str(), + vcg::tri::io::Mask::IOM_BITPOLYGONAL); + assert(saveWasSuccessful); + vcg::tri::UpdateNormal::PerFaceNormalized(dual); + + // vcg::PolygonalAlgorithm::SmoothReprojectPCA(dual, 1000, false, 0.5, 0, false, false); + // vcg::tri::io::ExporterOBJ::Save(dual, + // "./dual_optimized.obj", + // vcg::tri::io::Mask::IOM_BITPOLYGONAL); + + // if (vcg::tri::VoronoiProcessing::CheckVoronoiTopology(tileInto_triMesh, seedVec)) { + // VCGTriMesh voronoiMesh, voronoiPoly; + // vcg::tri::VoronoiProcessing::ConvertVoronoiDiagramToMesh(tileInto_triMesh, + // voronoiMesh, + // voronoiPoly, + // seedVec, + // vpp); + // vcg::tri::Allocator::CompactEveryVector(voronoiMesh); + // voronoiMesh.setLabel("Voro"); + // voronoiMesh.save(voronoiMesh.getLabel() + ".ply"); + // // voronoiMesh.registerForDrawing(); + // vcg::tri::Allocator::CompactEveryVector(voronoiPoly); + // voronoiPoly.setLabel("Poly"); + // voronoiPoly.save(voronoiPoly.getLabel() + ".ply"); + // // voronoiPoly.registerForDrawing(); + // // polyscope::show(); + // } + // dual.registerForDrawing(); + // polyscope::show(); + // assert(vcg::tri::IsValidPointer(dual, dual.face[0].cFEp(0))); + return std::make_shared(dual); +} + +} // namespace PolygonalRemeshing + +#endif // HEXAGONREMESHER_HPP diff --git a/src/main.cpp b/src/main.cpp new file mode 100755 index 0000000..51d6364 --- /dev/null +++ b/src/main.cpp @@ -0,0 +1,63 @@ +#include "edgemesh.hpp" +#include "hexagonremesher.hpp" +#include "polyscope/curve_network.h" +#include "polyscope/polyscope.h" +#include "polyscope/surface_mesh.h" +#include "reducedmodeloptimizer_structs.hpp" +#include "reducedpatternsimulator.hpp" +#include "simulationmesh.hpp" +#include "trianglepatterngeometry.hpp" +#include + +int main(int argc, char *argv[]) +{ + // DRMSimulationModel::runUnitTests(); + // const std::string tileInto_triMesh_filename( + // "/home/iason/Coding/build/PatternTillingReducedModel/Meshes/" + // "instantMeshes_strip_45.ply"); + // const std::string tileInto_triMesh_filename( + // "/home/iason/Coding/build/PatternTillingReducedModel/Meshes/" + // "instantMeshes_widerStrip_100.ply"); + const std::string tileInto_triMesh_filename( + "/home/iason/Coding/build/PatternTillingReducedModel/Meshes/" + "instantMeshes_plane_34.ply"); + // const std::string tileInto_triMesh_filename( + // "/home/iason/Coding/build/PatternTillingReducedModel/RelWithDebInfo/" + // "instantMeshes_plane4Hexagon_5.ply"); + VCGTriMesh tileInto_triMesh; + tileInto_triMesh.load(tileInto_triMesh_filename); + + std::shared_ptr pTileInto_polyMesh = PolygonalRemeshing::remeshWithPolygons( + tileInto_triMesh); + + //Load optimization results + // const size_t numberOfOptimizationResults = argc - 1; + // std::vector optimizationResults( + // numberOfOptimizationResults); + // for (size_t commandLineParameterIndex = 0; + // commandLineParameterIndex < numberOfOptimizationResults; + // commandLineParameterIndex++) + // optimizationResults[commandLineParameterIndex].load(argv[commandLineParameterIndex + 1]); + + // for (const auto &dirEntry : std::filesystem::directory_iterator( + std::filesystem::path optimizationResultsPath0( + "/home/iason/Coding/build/ReducedModelOptimization/RelWithDebInfo/OptimizationResults/" + "ConvergedJobs"); + std::filesystem::path optimizationResultsPath1( + "/home/iason/Coding/Projects/Approximating shapes with flat " + "patterns/ReducedModelOptimization/Results/selectionOfPatterns_0.2To1.6/" + "selectionOfPatterns/1.2/ConvergedJobs"); + std::vector optimizationResults + = ReducedPatternSimulator::loadOptimizationResults(optimizationResultsPath1); + //Load reduced + const std::filesystem::path patternFilePath( + "/home/iason/Coding/Projects/Approximating shapes with flat " + "patterns/ReducedModelOptimization/TestSet/ReducedPatterns/single_reduced.ply"); + PatternGeometry reducedPattern(patternFilePath.string()); + reducedPattern.setLabel(std::filesystem::path(patternFilePath).stem().string()); + + PolyscopeInterface::init(); + ReducedPatternSimulator reducedPatternSimulator(optimizationResults); + reducedPatternSimulator.simulate(pTileInto_polyMesh, optimizationResults, reducedPattern); + return 0; +} diff --git a/src/reducedpatternsimulator.cpp b/src/reducedpatternsimulator.cpp new file mode 100644 index 0000000..fa2872c --- /dev/null +++ b/src/reducedpatternsimulator.cpp @@ -0,0 +1,1183 @@ +#include "reducedpatternsimulator.hpp" +#include "imgui.h" +#include "misc/cpp/imgui_stdlib.h" +#include "simulationhistoryplotter.hpp" +#include + +void ReducedPatternSimulator::applyOptimizationResults_elements( + const ReducedPatternOptimization::Results &reducedPattern_optimizationResults, + const shared_ptr &pTiledReducedPattern_simulationMesh, + const std::vector &elementIndices) +{ + assert(CrossSectionType::name == RectangularBeamDimensions::name); + // Set reduced parameters fron the optimization results + std::unordered_map + optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(), + reducedPattern_optimizationResults.optimalXNameValuePairs.end()); + + const std::string ALabel = "A"; + assert(optimalXVariables.contains(ALabel)); + const double A = optimalXVariables.at(ALabel); + const double beamWidth = std::sqrt(A); + const double beamHeight = beamWidth; + CrossSectionType elementDimensions(beamWidth, beamHeight); + + const double poissonsRatio = 0.3; + const std::string ymLabel = "E"; + assert(optimalXVariables.contains(ymLabel)); + const double E = optimalXVariables.at(ymLabel); + const ElementMaterial elementMaterial(poissonsRatio, E); + + const std::string JLabel = "J"; + assert(optimalXVariables.contains(JLabel)); + const double J = optimalXVariables.at(JLabel); + + const std::string I2Label = "I2"; + assert(optimalXVariables.contains(I2Label)); + const double I2 = optimalXVariables.at(I2Label); + + const std::string I3Label = "I3"; + assert(optimalXVariables.contains(I3Label)); + const double I3 = optimalXVariables.at(I3Label); + for (size_t ei : elementIndices) { + Element &e = pTiledReducedPattern_simulationMesh->elements[ei]; + e.setDimensions(elementDimensions); + e.setMaterial(elementMaterial); + e.J = J; + e.I2 = I2; + e.I3 = I3; + } +} + +std::pair ReducedPatternSimulator::runSimulations( + const std::string &scenarioLabel, + const bool &shouldRerunFullPatternSimulation, + const bool &showReducedPatternResultsBeforeRunningDRM) +{ + if (pJob_tiledFullPattern->isEmpty() || pJob_tiledFullPattern->isEmpty()) { + std::cout << "Empty simulation job. Load or create a simulation job to run a simulation." + << std::endl; + return std::make_pair(SimulationResults(), SimulationResults()); + } + // std::cout << "Running " << scenarioLabel << std::endl; + pJob_tiledFullPattern->label = scenarioLabel; + pJob_tiledReducedPattern->label = scenarioLabel; + const std::filesystem::path outputPath( + std::filesystem::path("/home/iason/Coding/build/PatternTillingReducedModel/Scenarios/") + .append(scenarioLabel)); + const std::string simJobsOutputFolderPath( + std::filesystem::path(outputPath).append("SimulationJobs").string()); + + std::filesystem::create_directories(simJobsOutputFolderPath); + // if (!std::filesystem::exists(simJobsOutputFolderPath)) { + saveJobs(simJobsOutputFolderPath); + // } + //Reduced + LinearSimulationModel linearSimulationModel; + SimulationResults simulationResults_reducedPattern = linearSimulationModel.executeSimulation( + pJob_tiledReducedPattern); + if (showReducedPatternResultsBeforeRunningDRM) { + simulationResults_reducedPattern.registerForDrawing(); + polyscope::show(); + simulationResults_reducedPattern.unregister(); + } + // const auto reducedResultsFolderPath = std::filesystem::path(outputPath) + // .append("ReducedResults") + // .append(fullPatternsLabel); + // std::filesystem::create_directories(reducedResultsFolderPath); + // simulationResults_reducedPattern.save(reducedResultsFolderPath); + + //Full + const auto fullResultsFolderPath + = std::filesystem::path(outputPath).append("FullResults").append(fullPatternsLabel); + SimulationResults simulationResults_fullPattern; + if (shouldRerunFullPatternSimulation && std::filesystem::exists(fullResultsFolderPath)) { + std::filesystem::remove_all(fullResultsFolderPath); + } + + const std::filesystem::path fullPatternJobFolderPath + = std::filesystem::path(simJobsOutputFolderPath).append(fullPatternsLabel); + if (std::filesystem::exists(fullResultsFolderPath) && !randomTesselationIsEnabled) { + //Load full pattern results + assert(std::filesystem::exists(fullPatternJobFolderPath)); + simulationResults_fullPattern.load(fullResultsFolderPath, fullPatternJobFolderPath); + simulationResults_fullPattern.converged = true; + } else { + //Full + DRMSimulationModel drmSimulationModel; + DRMSimulationModel::Settings drmSimulationSettings; + drmSimulationSettings.isDebugMode = true; + drmSimulationSettings.shouldDraw = true; + drmSimulationSettings.drawingStep = 100000; + drmSimulationSettings.beVerbose = false; + drmSimulationSettings.shouldCreatePlots = true; + drmSimulationSettings.Dtini = 0.1; + drmSimulationSettings.totalExternalForcesNormPercentageTermination = 0.001; + // drmSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true; + // drmSimulationSettings.totalTranslationalKineticEnergyThreshold = 1e-10; + + // PolyscopeInterface::deinitPolyscope(); + simulationResults_fullPattern = drmSimulationModel.executeSimulation(pJob_tiledFullPattern, + drmSimulationSettings); + if (!simulationResults_fullPattern.converged) { + std::cerr << "Full pattern simulation failed." << std::endl; + } + + if (!randomTesselationIsEnabled && simulationResults_fullPattern.converged) { + std::filesystem::create_directories(fullResultsFolderPath); + simulationResults_fullPattern.save(fullResultsFolderPath); + pJob_tiledFullPattern->save(fullPatternJobFolderPath); + } + } + + return std::make_pair(simulationResults_fullPattern, simulationResults_reducedPattern); +} + +void ReducedPatternSimulator::tileReducedPatterns( + const std::vector &optimizationResults, + const std::vector &perSurfaceFacePatternIndices, + std::shared_ptr &pSimulationMesh_tiledReduced, + std::vector &tileIntoEdgeToTiledReducedVi) +{ + std::vector reducedPatterns(optimizationResults.size()); + for (size_t reducedPatternIndex = 0; reducedPatternIndex < optimizationResults.size(); + reducedPatternIndex++) { + reducedPatterns[reducedPatternIndex].copy(reducedPattern); + ReducedPatternOptimization::Results::applyOptimizationResults_innerHexagon( + optimizationResults[reducedPatternIndex], + reducedPatternBaseTriangle, + reducedPatterns[reducedPatternIndex]); + } + std::vector> perPatternIndexTiledReducedPatternEdgeIndices; + std::shared_ptr pTiledReducedPattern + = PatternGeometry::tilePattern(reducedPatterns, + {0}, + *pTileIntoSurface, + perSurfaceFacePatternIndices, + tileIntoEdgeToTiledReducedVi, + perPatternIndexTiledReducedPatternEdgeIndices); + pTiledReducedPattern->setLabel("Tiled_reduced_patterns"); + + std::size_t shuffleHash = computeHashOrdered(perSurfaceFacePatternIndices); + if (!shuffleToNumOfOccur.contains(shuffleHash)) { + shuffleToNumOfOccur[shuffleHash] = 0; + } else { + shuffleToNumOfOccur[shuffleHash]++; + } + if (shuffleToNumOfOccur[shuffleHash] != 0) { + std::cout << shuffleToNumOfOccur[shuffleHash] << std::endl; + } + + // polyscope::show(); + // pTiledReducedPattern->unregister(); + + pSimulationMesh_tiledReduced.reset(); + pSimulationMesh_tiledReduced = std::make_shared(*pTiledReducedPattern); + + for (size_t reducedPatternIndex = 0; reducedPatternIndex < optimizationResults.size(); + reducedPatternIndex++) { + const std::vector &tiledPatternElementIndicesForReducedPattern + = perPatternIndexTiledReducedPatternEdgeIndices[reducedPatternIndex]; + applyOptimizationResults_elements(optimizationResults[reducedPatternIndex], + pTiledReducedPattern_simulationMesh, + tiledPatternElementIndicesForReducedPattern); + } + pTiledReducedPattern_simulationMesh->reset(); + +} + +void ReducedPatternSimulator::shuffleReducedPatterns( + const std::vector &optimizationResults, + const std::vector &tileIntoEdgeToTiledFullPattern, + const std::vector &perSurfaceFacePatternIndices) +{ + std::vector randomShufflePerSurfaceFacePatternIndices + = perSurfaceFacePatternIndices; + std::random_shuffle(randomShufflePerSurfaceFacePatternIndices.begin(), + randomShufflePerSurfaceFacePatternIndices.end()); + std::vector tileIntoEdgesToTiledReducedVi; + tileReducedPatterns(optimizationResults, + randomShufflePerSurfaceFacePatternIndices, + pTiledReducedPattern_simulationMesh, + tileIntoEdgesToTiledReducedVi); + + constructViMaps(tileIntoEdgeToTiledFullPattern, tileIntoEdgesToTiledReducedVi); +} + +double ReducedPatternSimulator::computeDistance(const std::string &scenarioLabel, + const bool &shouldRerunDRMSimulation, + const bool &shouldDraw) +{ + loadScenario(scenarioLabel); + const auto results = runSimulations(scenarioLabel, shouldRerunDRMSimulation); + if (shouldDraw) { + results.first.registerForDrawing(); + results.second.registerForDrawing(); + polyscope::show(); + results.first.unregister(); + results.second.unregister(); + } + const double distance = computeDistance(results.first, results.second, fullToReducedViMap) + ; + return distance; +} + +std::vector ReducedPatternSimulator::computeDistancesPerScenario( + const std::vector &scenarioLabels, + const bool &shouldRerunDRMSimulation, + const bool &shouldDraw) +{ + std::vector fullReducedDistancesPerScenario(scenarioLabels.size()); + for (size_t scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) { + const std::string &scenarioLabel = scenarioLabels[scenarioIndex]; + // std::cout << "Running:" << scenarioLabel << std::endl; + fullReducedDistancesPerScenario[scenarioIndex] = computeDistance(scenarioLabel, + shouldRerunDRMSimulation, + shouldDraw); + // std::cout << "Distance:" << fullReducedDistancesPerScenario[scenarioIndex] << std::endl; + reset(); + } + return fullReducedDistancesPerScenario; +} + +std::vector ReducedPatternSimulator::computePerSurfaceFacePatternsIndices( + const std::vector &optimizationResults) const +{ + std::vector perSurfaceFacePatternIndices(pTileIntoSurface->FN()); + if (randomTesselationIsEnabled) { + for (size_t tileIntoFi = 0; tileIntoFi < pTileIntoSurface->FN(); tileIntoFi++) { + const size_t randomPatternIndex = (rand() + % static_cast(optimizationResults.size())); + perSurfaceFacePatternIndices[tileIntoFi] = randomPatternIndex; + } + } else { + const size_t faceIndexStep = pTileIntoSurface->FN() / optimizationResults.size(); + for (OptimizationResultsIndex optResIndex = 0; optResIndex < optimizationResults.size(); + optResIndex++) { + std::fill_n(perSurfaceFacePatternIndices.begin() + optResIndex * faceIndexStep, + faceIndexStep, + optResIndex); + } + } + return perSurfaceFacePatternIndices; +} + +std::vector ReducedPatternSimulator::loadOptimizationResults( + const std::filesystem::path &optimizationResultsFolderPath) +{ + std::set optimizationResultPaths; + for (const auto &dirEntry : std::filesystem::directory_iterator(optimizationResultsFolderPath)) { + // std::cout << "Loading:" << dirEntry.path() << std::endl; + optimizationResultPaths.insert(dirEntry.path()); + } + std::vector optimizationResults( + optimizationResultPaths.size()); + size_t optimizationResultIndex = 0; + for (const auto &p : optimizationResultPaths) { + optimizationResults[optimizationResultIndex++].load(p); + } + + return optimizationResults; +} +void ReducedPatternSimulator::tileFullPatterns( + std::vector &fullPatterns, + const std::vector &perSurfaceFacePatternIndices, + std::shared_ptr &pTiledFullPattern_simulationMesh, + std::vector &tileIntoEdgeToTiledFullVi) +{ + ///Create tiled configurations + std::vector> perPatternIndexTiledFullPatternEdgeIndices; + std::shared_ptr pTiledFullPattern + = PatternGeometry::tilePattern(fullPatterns, + {}, + *pTileIntoSurface, + perSurfaceFacePatternIndices, + tileIntoEdgeToTiledFullVi, + perPatternIndexTiledFullPatternEdgeIndices); + pTiledFullPattern->setLabel("Tiled_full_patterns"); + //#ifdef POLYSCOPE_DEFINED + pTiledFullPattern->registerForDrawing(); + // polyscope::show(); + // pTiledFullPattern->unregister(); + //#endif + + pTiledFullPattern_simulationMesh.reset(); + pTiledFullPattern_simulationMesh = std::make_shared(*pTiledFullPattern); + //NOTE: Those should be derived from the optimization results instead of hardcoding them + pTiledFullPattern_simulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002}); + pTiledFullPattern_simulationMesh->setBeamMaterial(0.3, 1 * 1e9); +} + +void ReducedPatternSimulator::constructViMaps( + const std::vector &tileIntoEdgeToTiledFullPattern, + const std::vector &tileIntoEdgeToTiledReducedPattern) +{ + fullToReducedViMap.clear(); + for (int ei = 0; ei < pTileIntoSurface->EN(); ei++) { + fullToReducedViMap[tileIntoEdgeToTiledFullPattern[ei]] + = tileIntoEdgeToTiledReducedPattern[ei]; + } + constructInverseMap(fullToReducedViMap, reducedToFullViMap); + + tilledFullPatternInterfaceVi.clear(); + tilledFullPatternInterfaceVi.resize(fullToReducedViMap.size()); + size_t viIndex = 0; + for (std::pair fullToReducedPair : fullToReducedViMap) { + tilledFullPatternInterfaceVi[viIndex++] = fullToReducedPair.first; + } +} + +void ReducedPatternSimulator::createShufflings( + std::vector &optimizationResults) +{ + std::vector shuffledPerSurfaceFacePatternIndices = computePerSurfaceFacePatternsIndices( + optimizationResults); + std::random_shuffle(shuffledPerSurfaceFacePatternIndices.begin(), + shuffledPerSurfaceFacePatternIndices.end()); + + createTiledSimulationMeshes(pTileIntoSurface, + optimizationResults, + shuffledPerSurfaceFacePatternIndices); +} + +void ReducedPatternSimulator::runWeightEvaluation( + const std::filesystem::path &optimizationResultsFolderPath) +{ + assert(pTileIntoSurface->VN() == 62 && pTileIntoSurface->EN() == 83 + && pTileIntoSurface->FN() == 22); + std::filesystem::path bestPerformingFolder; + double minDist = std::numeric_limits::max(); + std::vector> shufflings = csvFile::parse( + std::filesystem::path("/home/iason/Coding/build/PatternTillingReducedModel/RelWithDebInfo") + .append("shufflings.csv")); + const size_t numberOfShufflings = shufflings.size(); + for (auto &p : std::filesystem::directory_iterator(optimizationResultsFolderPath)) { + if (!std::filesystem::is_directory(std::filesystem::path(p))) { + continue; + } + std::vector optimizationResults + = loadOptimizationResults(std::filesystem::path(p).append("ConvergedJobs")); + double averageDistanceOverAllShufflings = 0; + // computePerSurfaceFacePatternsIndices(optimizationResults); + for (std::vector shuffledPerSurfaceFacePatternIndices : shufflings) { + // for (const auto &el : shuffledPerSurfaceFacePatternIndices) { + // std::cout << el << " "; + // } + // std::cout << std::endl; + createTiledSimulationMeshes(pTileIntoSurface, + optimizationResults, + shuffledPerSurfaceFacePatternIndices); + std::vector perScenarioDistances = computeDistancesPerScenario( + scenariosTestSetLabels); + averageDistanceOverAllShufflings += std::accumulate(perScenarioDistances.begin(), + perScenarioDistances.end(), + 0.0) + / numberOfShufflings; + } + if (averageDistanceOverAllShufflings < minDist) { + bestPerformingFolder = p; + minDist = averageDistanceOverAllShufflings; + } + std::cout << "Evaluated:" << p.path().string() << std::endl; + std::cout << "Dist:" << averageDistanceOverAllShufflings << std::endl; + } + + std::cout << "The best performing set was at:" << bestPerformingFolder << std::endl; +} + +void ReducedPatternSimulator::runShufflingEvaluation( + const std::vector &optimizationResults, + const std::vector &tileIntoEiToTiledFullVi) +{ + assert(pTileIntoSurface->VN() == 62 && pTileIntoSurface->EN() == 83 + && pTileIntoSurface->FN() == 22); + std::vector perSurfaceFacePatternIndices = computePerSurfaceFacePatternsIndices( + optimizationResults); + std::vector nonShuffledDistances = computeDistancesPerScenario(scenariosTestSetLabels); + const size_t numberOfShufflings = 500; + std::vector> shuffledDistances(numberOfShufflings, + std::vector( + scenariosTestSetLabels.size())); + for (size_t shufflingIndex = 0; shufflingIndex < numberOfShufflings; shufflingIndex++) { + shuffleReducedPatterns(optimizationResults, + tileIntoEiToTiledFullVi, + perSurfaceFacePatternIndices); + shuffledDistances[shufflingIndex] = computeDistancesPerScenario(scenariosTestSetLabels); + // for (size_t scenarioIndex = 0; scenarioIndex < scenariosTestSetLabels.size(); + // scenarioIndex++) { + // const bool shuffledPerformsBetter = shuffledDistances[shufflingIndex][scenarioIndex] + // < nonShuffledDistances[scenarioIndex]; + // if (shuffledPerformsBetter) { + // computeDistance(scenariosTestSetLabels[scenarioIndex], false, true); + // } + // } + } + + //Compute results + size_t numberOfShuffledTesselationsPerformingBetter = 0; + size_t totalNumberOfScenariosPerformingBetter = 0; + std::vector numberOfBetterPerformingShufflinsPerScenario(scenariosTestSetLabels.size(), + 0); + for (size_t shufflingIndex = 0; shufflingIndex < numberOfShufflings; shufflingIndex++) { + size_t scenariosPerformingBetter = 0; + for (size_t scenarioIndex = 0; scenarioIndex < scenariosTestSetLabels.size(); + scenarioIndex++) { + const bool shuffledPerformsBetter = shuffledDistances[shufflingIndex][scenarioIndex] + < nonShuffledDistances[scenarioIndex]; + // if (shuffledPerformsBetter) { + // int i = 0; + // i++; + // } + scenariosPerformingBetter += shuffledPerformsBetter; + numberOfBetterPerformingShufflinsPerScenario[scenarioIndex] += shuffledPerformsBetter; + } + // const float performingBetterThreshold = 0.8; + // if (static_cast(scenariosPerformingBetter) / scenariosTestSetLabels.size() + // >= performingBetterThreshold) { + // numberOfShuffledTesselationsPerformingBetter++; + // } + + totalNumberOfScenariosPerformingBetter += scenariosPerformingBetter; + } + matplot::bar(numberOfBetterPerformingShufflinsPerScenario); + + std::vector scenarioTestSetSimplifiedLabels(scenariosTestSetLabels.size()); + std::transform( + scenariosTestSetLabels.begin(), + scenariosTestSetLabels.end(), + scenarioTestSetSimplifiedLabels.begin(), + [](const std::string &label) { + std::string simplifiedLabel(label.begin() + 6, label.end()); + simplifiedLabel.erase(std::remove(simplifiedLabel.begin(), simplifiedLabel.end(), '_'), + simplifiedLabel.end()); + // if (simplifiedLabel.size() > 12) { + // simplifiedLabel.insert(simplifiedLabel.begin() + 12, "\n"); + // } + // simplifiedLabel = "a\na"; + const size_t beginPos = simplifiedLabel.find("pullOppositeVerts"); + if (beginPos != simplifiedLabel.npos) { + simplifiedLabel.replace(beginPos, std::string("pullOppositeVerts").size(), "POV"); + } + + return simplifiedLabel; + }); + matplot::gca()->x_axis().ticklabels(scenarioTestSetSimplifiedLabels); + matplot::gca()->y_axis().limits(std::array{0, numberOfShufflings}); + matplot::show(); + + std::cout << "Number of shufflings:" << numberOfShufflings << std::endl; + std::cout << "Average scenarios performing better/shuffling:" + << static_cast(totalNumberOfScenariosPerformingBetter) / numberOfShufflings + << std::endl; + // std::cout << "Number of shufflings performing better:" + // << numberOfShuffledTesselationsPerformingBetter << std::endl; + // std::cout << "Percentage:" + // << static_cast(numberOfShuffledTesselationsPerformingBetter) + // / numberOfShufflings + // << std::endl; +} + +std::pair +ReducedPatternSimulator::getPickedVertices(const std::pair &selection) const +{ + std::shared_ptr &pTiledFullPattern_simulationMesh = pJob_tiledFullPattern->pMesh; + assert(pTiledFullPattern_simulationMesh != nullptr); + std::shared_ptr &pTiledReducedPattern_simulationMesh = pJob_tiledReducedPattern + ->pMesh; + assert(pTiledReducedPattern_simulationMesh != nullptr); + const int &vi = selection.second; + const std::string &pickedStructureName = selection.first; + assert(!pickedStructureName.empty()); + int fullVi = -1; + int reducedVi = -1; + if (pickedStructureName == pTiledFullPattern_simulationMesh->getLabel()) { + fullVi = vi; + reducedVi = fullToReducedViMap.at(vi); + } else if (pickedStructureName == pTiledReducedPattern_simulationMesh->getLabel()) { + reducedVi = vi; + fullVi = reducedToFullViMap.at(vi); + } + + return std::make_pair(fullVi, reducedVi); +} + +void ReducedPatternSimulator::reset() +{ + minInputForceMagnitude = std::numeric_limits::max(); + + //Keep only initial patterns + for (std::pair> + polyscopeStructureCategory : polyscope::state::structures) { + for (std::pair polyscopeStructure : + polyscopeStructureCategory.second) { + if (polyscopeStructure.first != pTiledReducedPattern_simulationMesh->getLabel() + && polyscopeStructure.first != pTiledFullPattern_simulationMesh->getLabel()) { + polyscope::removeStructure(polyscopeStructure.second); + } else { + dynamic_cast(polyscopeStructure.second) + ->removeAllQuantities(); + } + } + } + //Remove simulation job quantities from the initial meshes + pJob_tiledReducedPattern->unregister(pTiledReducedPattern_simulationMesh->getLabel()); + pJob_tiledReducedPattern->clear(); + pJob_tiledFullPattern->unregister(pTiledFullPattern_simulationMesh->getLabel()); + pJob_tiledFullPattern->clear(); + polyscope::requestRedraw(); +} + +double ReducedPatternSimulator::computeDistance( + const SimulationResults &resultsA, + const SimulationResults &resultsB, + const std::unordered_map &resultsAToResultsBViMap) const +{ + double distance = 0; + for (std::pair resultsAToResultsBViPair : resultsAToResultsBViMap) { + const double vertexToVertexDistance + = (resultsA.displacements[resultsAToResultsBViPair.first].getTranslation() + - resultsB.displacements[resultsAToResultsBViPair.second].getTranslation()) + .norm(); + distance += vertexToVertexDistance; + } + return distance; +} + +void ReducedPatternSimulator::saveJobs(const std::filesystem::path &saveToPath) +{ + std::filesystem::create_directories(saveToPath); + const auto fullJobFolderPath = std::filesystem::path(saveToPath).append(fullPatternsLabel); + std::filesystem::create_directories(fullJobFolderPath); + pJob_tiledFullPattern->save(fullJobFolderPath.string()); + const auto reducedJobFolderPath = std::filesystem::path(saveToPath).append("Reduced"); + std::filesystem::create_directories(reducedJobFolderPath); + pJob_tiledReducedPattern->save(reducedJobFolderPath.string()); +} + +void ReducedPatternSimulator::saveResults(const std::string &outputFolderPath, + SimulationResults &fullResults, + SimulationResults &reducedResults) +{ + std::filesystem::create_directories(outputFolderPath); + const auto fullResultsFolderPath = std::filesystem::path(outputFolderPath).append("Full"); + std::filesystem::create_directories(fullResultsFolderPath); + const auto reducedResultsFolderPath = std::filesystem::path(outputFolderPath).append("Reduced"); + std::filesystem::create_directories(reducedResultsFolderPath); + fullResults.save(fullResultsFolderPath.string()); + reducedResults.save(reducedResultsFolderPath.string()); +} + +void ReducedPatternSimulator::generateRandomSimulationScenario( + const std::array &randomScenarioParameters) +{ + const float percentageOfRigidVertices = gui_randomnessParams[0]; + const float percentageOfFixedVertices = gui_randomnessParams[1]; + const float percentageOfNodalForcesVertices = gui_randomnessParams[2]; + assert(percentageOfFixedVertices <= 1 && percentageOfNodalForcesVertices <= 1 + && percentageOfRigidVertices <= 1); + const float maxAbsForceMagnitude = gui_randomnessParams[3]; + assert(maxAbsForceMagnitude >= 0); + std::random_device dev; + std::mt19937 rng(dev()); + std::uniform_int_distribution + distInt(0, tilledFullPatternInterfaceVi.size() - 1); + std::unordered_set usedVertices; + + size_t numberOfFixedVertices = 0; + while (numberOfFixedVertices / static_cast(tilledFullPatternInterfaceVi.size()) + < percentageOfFixedVertices) { + const size_t fullTilledVi = tilledFullPatternInterfaceVi[distInt(rng)]; + if (!usedVertices.contains(fullTilledVi)) { + pJob_tiledFullPattern->constrainedVertices[fullTilledVi] = {0, 1, 2}; + numberOfFixedVertices++; + usedVertices.insert(fullTilledVi); + } + } + + size_t numberOfRigidVertices = 0; + while (numberOfRigidVertices / static_cast(tilledFullPatternInterfaceVi.size()) + < percentageOfRigidVertices) { + const size_t fullTilledVi = tilledFullPatternInterfaceVi[distInt(rng)]; + if (!usedVertices.contains(fullTilledVi)) { + pJob_tiledFullPattern->constrainedVertices[fullTilledVi] = {0, 1, 2, 3, 4, 5}; + numberOfRigidVertices++; + usedVertices.insert(fullTilledVi); + } + } + + std::uniform_real_distribution<> distFloat(-maxAbsForceMagnitude, maxAbsForceMagnitude); + while (pJob_tiledFullPattern->nodalExternalForces.size() + / static_cast(tilledFullPatternInterfaceVi.size()) + < percentageOfNodalForcesVertices) { + Vector6d externalForce({distFloat(rng) / 20, distFloat(rng) / 20, distFloat(rng), 0, 0, 0}); + minInputForceMagnitude = std::min(minInputForceMagnitude, externalForce.norm()); + + const size_t fullTilledVi = tilledFullPatternInterfaceVi[distInt(rng)]; + if (!usedVertices.contains(fullTilledVi)) { + pJob_tiledFullPattern->nodalExternalForces[fullTilledVi] = externalForce; + usedVertices.insert(fullTilledVi); + } + } + + //Apply the reduced job to the full pattern + pJob_tiledFullPattern->remap(fullToReducedViMap, *pJob_tiledReducedPattern); + + pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh; + pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh; +} + +void ReducedPatternSimulator::reportDistances( + const std::tuple + &simulationResults, + const std::string &csvFileName) +{ + DRMFullSimulationResults drmFullTilledResults; + ReducedSimulationResults reducedTilledResults; + LinearFullSimulationResults linearFullTilledResults; + std::tie(drmFullTilledResults, reducedTilledResults, linearFullTilledResults) = simulationResults; + const double distance_fullDrmToReduced = computeDistance(drmFullTilledResults, + reducedTilledResults, + fullToReducedViMap); + // std::cout << "Distance full to reduced" << distance_fullDrmToReduced<< std::endl; + double distance_fullSumOfAllVerts = 0; + for (std::pair fullToReducedPair : fullToReducedViMap) { + distance_fullSumOfAllVerts + += drmFullTilledResults.displacements[fullToReducedPair.first].getTranslation().norm(); + } + // csvFile csv_results({}, false); + // csv_results << "Full Pattern Names" + // << "Simulation job label" + // << "Distance"; + // csv_results.endrow(); + // csv_results << fullPatternsSurfacelabel << gui_jobLabel << std::to_string(distance); + // csv_results.endrow(); + std::unordered_map fullToFullDummyViMap = fullToReducedViMap; + for (auto &viPair : fullToFullDummyViMap) { + viPair.second = viPair.first; + } + const double distance_fullLinearToDrm = computeDistance(drmFullTilledResults, + linearFullTilledResults, + fullToFullDummyViMap) + /*/ (2 * surfaceBaseTriangleHeight * fullToReducedViMap.size())*/; + // std::cout << "Distance tilled full linear vs drm:" << distance_fullLinearToDrm << std::endl; + // std::cout << "Distance full to reduced norm:" + // << distance_fullTilledToReducedTilled / distance_fullSumOfAllVerts << std::endl; + // std::cout << "Distance tilled full linear vs drm norm:" + // << distance_fullLinearToDrm / distance_fullSumOfAllVerts << std::endl; + + // Write results in csv + csvFile csv_results(csvFileName, false); + // csvFile csv_results(std::filesystem::path(dirPath_thisOptimization) + // .append("results.csv") + // .string(), + // false); + csv_results << "Job Label" + << "drm2Reduced" + << "drm2Linear" + << "norm_drm2Reduced" + << "norm_drm2Linear"; + csv_results << endrow; + csv_results << pJob_tiledFullPattern->getLabel() << distance_fullDrmToReduced + << distance_fullLinearToDrm + << distance_fullDrmToReduced / distance_fullSumOfAllVerts + << distance_fullLinearToDrm / distance_fullSumOfAllVerts; + csv_results << endrow; +} + +std::tuple +ReducedPatternSimulator::runAllSimulations(const std::string &jobLabel, + const bool &shouldRerunDRMFullPatternSimulation, + const bool &showReducedPatternResultsBeforeDRM) +{ + std::cout << "Executing scenario:" << jobLabel << std::endl; + std::pair results + = runSimulations(jobLabel, + shouldRerunDRMFullPatternSimulation, + showReducedPatternResultsBeforeDRM); + if (results.first.converged == false || results.second.converged == false) { + std::cerr << "Simulation did not converge" << std::endl; + return std::tuple(); + } + + LinearSimulationModel linearSimulationModel; + SimulationResults simulationResults_fullTilledLinearModel + = linearSimulationModel.executeSimulation(pJob_tiledFullPattern); + simulationResults_fullTilledLinearModel.setLabelPrefix("LinearSimModel"); + + return std::make_tuple(results.first, results.second, simulationResults_fullTilledLinearModel); +} + +void ReducedPatternSimulator::createGuiMenu() +{ + // std::shared_ptr &pTiledFullPattern_simulationMesh = pJob_tiledFullPattern->pMesh; + // assert(pTiledFullPattern_simulationMesh != nullptr); + // std::shared_ptr &pTiledReducedPattern_simulationMesh = pJob_tiledReducedPattern + // ->pMesh; + assert(pTiledReducedPattern_simulationMesh != nullptr); + PolyscopeInterface::addUserCallback([&]() { + // ImGuiIO &io = ImGui::GetIO(); + if (ImGui::Button("Fix vertex")) { + const std::pair selection = PolyscopeInterface::getSelection(); + const std::string &pickedStructureName = selection.first; + if (pickedStructureName.empty()) { + return; + } + std::pair pickedVertices + = getPickedVertices(selection); + pJob_tiledFullPattern->constrainedVertices[pickedVertices.first] = {0, 1, 2, 3, 4, 5}; + pJob_tiledReducedPattern->constrainedVertices[pickedVertices.second] = {0, 1, 2, 3, 4, 5}; + pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(), + true); + pJob_tiledReducedPattern + ->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true); + polyscope::requestRedraw(); + } + static int gui_dof = 0; + ImGui::InputInt("DoF", &gui_dof); + if (ImGui::Button("Constrain DoF of vertex")) { + const std::pair selection = PolyscopeInterface::getSelection(); + const std::string &pickedStructureName = selection.first; + if (pickedStructureName.empty()) { + return; + } + std::pair pickedVertices + = getPickedVertices(selection); + pJob_tiledFullPattern->constrainedVertices[pickedVertices.first] = {gui_dof}; + pJob_tiledReducedPattern->constrainedVertices[pickedVertices.second] = {gui_dof}; + pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(), + true); + pJob_tiledReducedPattern + ->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true); + polyscope::requestRedraw(); + } + + if (ImGui::Button("Constrain translation of vertex")) { + const std::pair selection = PolyscopeInterface::getSelection(); + const std::string &pickedStructureName = selection.first; + if (pickedStructureName.empty()) { + return; + } + std::pair pickedVertices + = getPickedVertices(selection); + pJob_tiledFullPattern->constrainedVertices[pickedVertices.first] = {0, 1, 2}; + pJob_tiledReducedPattern->constrainedVertices[pickedVertices.second] = {0, 1, 2}; + pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(), + true); + pJob_tiledReducedPattern + ->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true); + polyscope::requestRedraw(); + } + if (ImGui::Button("Remove constrains on vertex")) { + const std::pair selection = PolyscopeInterface::getSelection(); + const std::string &pickedStructureName = selection.first; + if (pickedStructureName.empty()) { + return; + } + std::pair pickedVertices + = getPickedVertices(selection); + pJob_tiledFullPattern->constrainedVertices.erase(pickedVertices.first); + pJob_tiledReducedPattern->constrainedVertices.erase(pickedVertices.second); + pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(), + true); + pJob_tiledReducedPattern + ->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true); + polyscope::requestRedraw(); + } + + if (ImGui::Button("Set all forces/moments")) { + const Vector6d desiredExternalLoad({gui_externalForce[0], + gui_externalForce[1], + gui_externalForce[2], + gui_externalMoment[0], + gui_externalMoment[1], + gui_externalMoment[2]}); + + for (auto &externalLoad : pJob_tiledFullPattern->nodalExternalForces) { + externalLoad.second = desiredExternalLoad; + } + pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(), + true); + for (auto &externalLoad : pJob_tiledReducedPattern->nodalExternalForces) { + externalLoad.second = desiredExternalLoad; + } + pJob_tiledReducedPattern + ->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true); + minInputForceMagnitude = std::min(minInputForceMagnitude, externalForce.norm()); + + polyscope::requestRedraw(); + } + if (ImGui::Button("Add force/moment")) { + const std::pair selection = PolyscopeInterface::getSelection(); + const std::string &pickedStructureName = selection.first; + if (pickedStructureName.empty()) { + return; + } + + std::pair pickedVertices + = getPickedVertices(selection); + Vector6d externalForce({gui_externalForce[0], + gui_externalForce[1], + gui_externalForce[2], + gui_externalMoment[0], + gui_externalMoment[1], + gui_externalMoment[2]}); + pJob_tiledFullPattern->nodalExternalForces[pickedVertices.first] = externalForce; + minInputForceMagnitude = std::min(minInputForceMagnitude, externalForce.norm()); + pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(), + true); + pJob_tiledReducedPattern->nodalExternalForces[pickedVertices.second] = externalForce; + minInputForceMagnitude = std::min(minInputForceMagnitude, externalForce.norm()); + pJob_tiledReducedPattern + ->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true); + polyscope::requestRedraw(); + } + ImGui::InputFloat3("Force", gui_externalForce.data()); + ImGui::InputFloat3("Moment", gui_externalMoment.data()); + + if (ImGui::Button("Remove Force/Moment")) { + const std::pair selection = PolyscopeInterface::getSelection(); + const std::string pickedStructureName = selection.first; + if (pickedStructureName.empty()) { + return; + } + std::pair pickedVertices + = getPickedVertices(selection); + pJob_tiledFullPattern->nodalExternalForces.erase(pickedVertices.first); + pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(), + true); + pJob_tiledReducedPattern->nodalExternalForces.erase(pickedVertices.second); + pJob_tiledReducedPattern + ->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true); + polyscope::requestRedraw(); + } + + static std::string gui_jobLabel; + ImGui::InputText("Job label", &gui_jobLabel); + + if (ImGui::Button("Load job")) { + reset(); + loadScenario(gui_jobLabel); + + polyscope::requestRedraw(); + } + + static bool gui_shouldRerunFullPatternSimulation = false; + ImGui::Checkbox("Re-run full pattern simulation", &gui_shouldRerunFullPatternSimulation); + + if (ImGui::Button("Run Simulations")) { + if (pJob_tiledFullPattern->isEmpty() || pJob_tiledFullPattern->isEmpty()) { + std::cerr + << "Empty simulation job. Load or create a simulation job to run a simulation." + << std::endl; + return; + } + DRMFullSimulationResults drmFullTilledResults; + ReducedSimulationResults reducedTilledResults; + LinearFullSimulationResults linearFullTilledResults; + std::tie(drmFullTilledResults, reducedTilledResults, linearFullTilledResults) + = runAllSimulations(gui_jobLabel, gui_shouldRerunFullPatternSimulation, true); + if (!drmFullTilledResults.converged || !reducedTilledResults.converged + || !linearFullTilledResults.converged) { + std::cerr << "A simulation did not converge." << std::endl; + return; + } + reportDistances({drmFullTilledResults, reducedTilledResults, linearFullTilledResults}); + drmFullTilledResults.registerForDrawing(); + reducedTilledResults.registerForDrawing(); + linearFullTilledResults.registerForDrawing(); + + polyscope::requestRedraw(); + } + + ImGui::InputFloat4("Randomness params", gui_randomnessParams.data()); + if (ImGui::Button("Generate random scenario")) { + gui_jobLabel = "22Hex_random"; + reset(); + + generateRandomSimulationScenario(gui_randomnessParams); + pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(), + true); + pJob_tiledReducedPattern + ->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true); + polyscope::requestRedraw(); + } + + if (ImGui::Button("Run random simulations")) { + const size_t numberOfRandomSimulations = 20; + gui_shouldRerunFullPatternSimulation = true; + for (size_t randomSimulationIndex = 0; + randomSimulationIndex < numberOfRandomSimulations; + randomSimulationIndex++) { + reset(); + gui_jobLabel = "22Hex_random" + std::to_string(randomSimulationIndex); + generateRandomSimulationScenario(gui_randomnessParams); + DRMFullSimulationResults drmFullTilledResults; + ReducedSimulationResults reducedTilledResults; + LinearFullSimulationResults linearFullTilledResults; + std::tie(drmFullTilledResults, reducedTilledResults, linearFullTilledResults) + = runAllSimulations(gui_jobLabel, gui_shouldRerunFullPatternSimulation, false); + reportDistances({drmFullTilledResults, + reducedTilledResults, + linearFullTilledResults}, + "distances.csv"); + // drmFullTilledResults.registerForDrawing(); + // reducedTilledResults.registerForDrawing(); + // linearFullTilledResults.registerForDrawing(); + // polyscope::requestRedraw(); + } + } + + if (ImGui::Button("Reset")) { + reset(); + } + }); +} + +void ReducedPatternSimulator::fillFullPatterns( + std::vector &optimizationResults) +{ + for (size_t fullPatternIndex = 0; fullPatternIndex < optimizationResults.size(); + fullPatternIndex++) { + // const vcg::Triangle3 fullPatternBaseTriangle = fullPattern.computeBaseTriangle(); + PatternGeometry &baseTriangleFullPattern = optimizationResults[fullPatternIndex] + .baseTriangleFullPattern; + fullPatterns[fullPatternIndex].copy(baseTriangleFullPattern); + vcg::tri::Allocator::PointerUpdater pu; + fullPatterns[fullPatternIndex].deleteDanglingVertices(pu); + if (!pu.remap.empty()) { + fullPatterns[fullPatternIndex].interfaceNodeIndex = pu.remap[3]; + } else { + fullPatterns[fullPatternIndex].interfaceNodeIndex = 3; + } + } +} + +void ReducedPatternSimulator::createTiledSimulationMeshes( + const std::shared_ptr &pTileIntoSurface, + std::vector &optimizationResults, + const std::vector &perSurfaceFacePatternIndices) +{ + std::vector tileIntoEdgeToTiledFullPattern; + createTiledSimulationMeshes(pTileIntoSurface, + optimizationResults, + perSurfaceFacePatternIndices, + tileIntoEdgeToTiledFullPattern); +} + +void ReducedPatternSimulator::createTiledSimulationMeshes( + const std::shared_ptr &pTileIntoSurface, + std::vector &optimizationResults, + const std::vector &perSurfaceFacePatternIndices, + std::vector &tileIntoEdgeToTiledFullPattern) +{ + assert(perSurfaceFacePatternIndices.size() == pTileIntoSurface->FN()); + + fillFullPatterns(optimizationResults); + tileFullPatterns(fullPatterns, + perSurfaceFacePatternIndices, + pTiledFullPattern_simulationMesh, + tileIntoEdgeToTiledFullPattern); + + std::vector tileIntoEdgeToTiledReducedPattern; + tileReducedPatterns(optimizationResults, + perSurfaceFacePatternIndices, + pTiledReducedPattern_simulationMesh, + tileIntoEdgeToTiledReducedPattern); + constructViMaps(tileIntoEdgeToTiledFullPattern, tileIntoEdgeToTiledReducedPattern); + computeLabels(optimizationResults, perSurfaceFacePatternIndices); + // std::vector reducedPatterns(optimizationResults.size()); + // for (size_t reducedPatternIndex = 0; reducedPatternIndex < optimizationResults.size(); + // reducedPatternIndex++) { + // reducedPatterns[reducedPatternIndex].copy(reducedPattern); + // ReducedPatternOptimization::Results::applyOptimizationResults_innerHexagon( + // optimizationResults[reducedPatternIndex], + // reducedPatternBaseTriangle, + // reducedPatterns[reducedPatternIndex]); + // } + + // std::vector randomShufflePerSurfaceFacePatternIndices + // = perSurfaceFacePatternIndices; + // std::random_shuffle(randomShufflePerSurfaceFacePatternIndices.begin(), + // randomShufflePerSurfaceFacePatternIndices.end()); + // std::vector tileIntoEdgeToTiledReducedPattern; + // std::vector> perPatternIndexTiledReducedPatternEdgeIndices; + // std::shared_ptr pTiledReducedPattern = PatternGeometry::tilePattern( + // reducedPatterns, + // {0}, + // *pTileIntoSurface, + // // randomShufflePerSurfaceFacePatternIndices, + // perSurfaceFacePatternIndices, + // tileIntoEdgeToTiledReducedPattern, + // perPatternIndexTiledReducedPatternEdgeIndices); + // pTiledReducedPattern->setLabel("Tiled_" + reducedPatterns[0].getLabel()); + //#ifdef POLYSCOPE_DEFINED + // pTiledReducedPattern->registerForDrawing(); + //#endif + + // pTiledReducedPattern_simulationMesh = std::make_shared(*pTiledReducedPattern); + + // for (size_t reducedPatternIndex = 0; reducedPatternIndex < optimizationResults.size(); + // reducedPatternIndex++) { + // const std::vector &tiledPatternElementIndicesForReducedPattern + // = perPatternIndexTiledReducedPatternEdgeIndices[reducedPatternIndex]; + // applyOptimizationResults_elements(optimizationResults[reducedPatternIndex], + // pTiledReducedPattern_simulationMesh, + // tiledPatternElementIndicesForReducedPattern); + // } + // pTiledReducedPattern_simulationMesh->reset(); + + // for (int ei = 0; ei < pTileIntoSurface->EN(); ei++) { + // fullToReducedViMap[tileIntoEdgeToTiledFullPattern[ei]] + // = tileIntoEdgeToTiledReducedPattern[ei]; + // } + // constructInverseMap(fullToReducedViMap, reducedToFullViMap); +} + +void ReducedPatternSimulator::loadScenario(const string &scenarioLabel) +{ + std::filesystem::path simulationJobsInputFolderPath( + "/home/iason/Coding/build/PatternTillingReducedModel/SimulationJobs/"); + simulationJobsInputFolderPath.append(scenarioLabel); + + if (!std::filesystem::exists(simulationJobsInputFolderPath) + || !std::filesystem::is_directory(simulationJobsInputFolderPath)) { + simulationJobsInputFolderPath + = std::filesystem::path("/home/iason/Coding/build/PatternTillingReducedModel/Scenarios") + .append(scenarioLabel) + .append("SimulationJobs"); + } + + if (!std::filesystem::is_directory(simulationJobsInputFolderPath)) { + std::cerr << "The given path is not a directory:" << simulationJobsInputFolderPath + << std::endl; + return; + } + // //Full + // std::filesystem::path fullPatternSimulationJobFilePath; + // for (auto &p : std::filesystem::recursive_directory_iterator( + // std::filesystem::path(inputFolderPath).append("Full"))) { + // if (p.path().extension() == ".json") + // fullPatternSimulationJobFilePath = p.path(); + // } + // assert(pTiledFullPattern_simulationMesh != nullptr); + // pJob_tiledFullPattern->load(fullPatternSimulationJobFilePath.string()); + // // pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh; + // assert(pTiledFullPattern_simulationMesh != nullptr); + // pTiledFullPattern_simulationMesh->registerForDrawing(); + // pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(), + // true); + // for (const auto &externalLoad : pJob_tiledFullPattern->nodalExternalForces) { + // minInputForceMagnitude = std::min(minInputForceMagnitude, + // externalLoad.second.norm()); + // } + //Reduced + std::filesystem::path reducedPatternSimulationJobFilePath; + for (auto &p : std::filesystem::directory_iterator(std::filesystem::path( + std::filesystem::path(simulationJobsInputFolderPath).append("Reduced")))) { + if (p.path().extension() == ".json") + reducedPatternSimulationJobFilePath = p.path(); + } + pJob_tiledReducedPattern->load(reducedPatternSimulationJobFilePath.string(), false); + //NOTE: Is there a way to not apply the optimization results again? I only need the Job and not the mesh + // pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh; + // applyOptimizationResults_elements(optimizationResults, + // pTiledReducedPattern_simulationMesh); + pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh; + assert(pTiledReducedPattern_simulationMesh != nullptr); + // pTiledReducedPattern_simulationMesh->registerForDrawing(); + pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), + true); + pJob_tiledFullPattern->label = scenarioLabel; + pJob_tiledReducedPattern->label = scenarioLabel; + + //Apply the reduced job to the full pattern + pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh; + pJob_tiledReducedPattern->remap(reducedToFullViMap, *pJob_tiledFullPattern); + pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(), true); +} + +void ReducedPatternSimulator::computeLabels( + const std::vector &optimizationResults, + const std::vector &perSurfaceFacePatternIndex) +{ + fullPatternsLabel.clear(); + for (OptimizationResultsIndex optResIndex = 0; optResIndex < optimizationResults.size(); + optResIndex++) { + fullPatternsLabel += optimizationResults[optResIndex].baseTriangleFullPattern.getLabel() + + "_"; + } + fullPatternsLabel += std::to_string(computeHashOrdered(perSurfaceFacePatternIndex)); + // fullPatternsLabel += "_" + std::to_string(optimizationResults[0].settings.objectiveWeights.translational); + fullPatternsSurfacelabel = fullPatternsLabel + "_" + pTileIntoSurface->getLabel(); +} + +ReducedPatternSimulator::ReducedPatternSimulator( + const std::vector &optimizationResults) + : fullPatterns(optimizationResults.size()) +{} + +void ReducedPatternSimulator::simulate( + std::shared_ptr &tileIntoSurface, + std::vector &optimizationResults, + PatternGeometry &reducedPattern) +{ + reducedPatternBaseTriangle = reducedPattern.computeBaseTriangle(); + this->reducedPattern.copy(reducedPattern); + this->reducedPattern.deleteDanglingVertices(); + this->reducedPattern.interfaceNodeIndex = 1; + this->pTileIntoSurface = tileIntoSurface; + + ///Preprocess tile into surface + /// NOTE: I am assuming here that the base triangle is the same for all optimization results + const double optimizedBaseTriangleHeight + = vcg::Distance(optimizationResults[0].baseTriangle.cP(0), + (optimizationResults[0].baseTriangle.cP(1) + + optimizationResults[0].baseTriangle.cP(2)) + / 2); + preprocessSurface(optimizedBaseTriangleHeight, *pTileIntoSurface); + vcg::tri::io::ExporterOBJ::Save(*pTileIntoSurface, + "tileIntoSurface.obj", + vcg::tri::io::Mask::IOM_BITPOLYGONAL); + + std::vector perSurfaceFacePatternIndices = computePerSurfaceFacePatternsIndices( + optimizationResults); + std::vector tileIntoEiToTiledFullPatternVi; + createTiledSimulationMeshes(tileIntoSurface, + optimizationResults, + perSurfaceFacePatternIndices, + tileIntoEiToTiledFullPatternVi); + + ///Create simulation jobs + pJob_tiledFullPattern = std::make_shared(SimulationJob()); + // pJob_tiledFullPattern->label="TiledFull"; + pJob_tiledReducedPattern = std::make_shared(SimulationJob()); + // pJob_tiledReducedPattern->setLabel("TiledReduced"); + // createSimulationJobs(externalForce, *pJob_tiledFullPattern, *pJob_tiledReducedPattern); + pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh; + pTiledFullPattern_simulationMesh->registerForDrawing(); + pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh; + pTiledReducedPattern_simulationMesh->registerForDrawing(); + // pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel()); + // pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel()); + + // runShufflingEvaluation(optimizationResults, tileIntoEiToTiledFullPatternVi); + // runWeightEvaluation("/home/iason/Coding/Projects/Approximating shapes with flat " + // "patterns/ReducedModelOptimization/Results/selectionOfPatterns_0.2To1.6/" + // "selectionOfPatterns"); + // computeDistancesPerScenario(scenariosTestSetLabels, true); + createGuiMenu(); + polyscope::show(); +} + +void ReducedPatternSimulator::preprocessSurface(const double &desiredBaseTriangleHeight, + VCGPolyMesh &surface) +{ + surface.moveToCenter(); + // const double desiredAverageBaseTriangleSize = des; //TODO: obtain from optimization results + vcg::tri::UpdatePosition::Scale(surface, + desiredBaseTriangleHeight + / surface.getAverageFaceRadius()); + surfaceBaseTriangleHeight = desiredBaseTriangleHeight; +} diff --git a/src/reducedpatternsimulator.hpp b/src/reducedpatternsimulator.hpp new file mode 100644 index 0000000..a896f61 --- /dev/null +++ b/src/reducedpatternsimulator.hpp @@ -0,0 +1,171 @@ +#ifndef REDUCEDPATTERNSIMULATOR_HPP +#define REDUCEDPATTERNSIMULATOR_HPP + +#include "reducedmodeloptimizer_structs.hpp" +#include "simulationmesh.hpp" +#include "trianglepatterngeometry.hpp" +#include "utilities.hpp" + +using DRMFullSimulationResults = SimulationResults; +using LinearFullSimulationResults = SimulationResults; +using ReducedSimulationResults = SimulationResults; +using FullPatternVertexIndex = int; +using ReducedPatternVertexIndex = int; +class ReducedPatternSimulator +{ + std::shared_ptr pTiledFullPattern_simulationMesh; + std::shared_ptr pTiledReducedPattern_simulationMesh; + std::shared_ptr pJob_tiledReducedPattern; + std::shared_ptr pJob_tiledFullPattern; + const Vector6d externalForce{0, 0, 0.1, 0, 0, 0}; + double minInputForceMagnitude = std::numeric_limits::max(); + std::array gui_externalForce{0, 0, 0}; + std::array gui_externalMoment{0, 0, 0}; + std::array gui_randomnessParams{0.02, 0.02, 0.3, 0.4}; + std::unordered_map fullToReducedViMap; //of only the common vertices + std::unordered_map reducedToFullViMap; //of only the common vertices + std::vector tilledFullPatternInterfaceVi; + double surfaceBaseTriangleHeight{-1}; + std::string fullPatternsSurfacelabel; + std::string fullPatternsLabel; + + vcg::Triangle3 reducedPatternBaseTriangle; + PatternGeometry reducedPattern; + const bool randomTesselationIsEnabled{false}; + const std::vector scenariosTestSetLabels{// "22Hex_random0", + // "22Hex_random1", + // "22Hex_random2", + // "22Hex_random3", + // "22Hex_random4", + // "22Hex_random5", + // "22Hex_random6", + // "22Hex_random7", + // "22Hex_random8", + // "22Hex_random9", + // "22Hex_random10", + // "22Hex_random11", + // "22Hex_random12", + // "22Hex_random17", + // "22Hex_random14" + "22Hex_bending_0.005N", + "22Hex_bending_0.01N", + "22Hex_bending_0.03N", + "22Hex_bending_0.05N", + "22Hex_pullOppositeVerts_0.05N", + "22Hex_pullOppositeVerts_0.1N", + "22Hex_pullOppositeVerts_0.3N", + "22Hex_shear_2N", + "22Hex_shear_5N", + "22Hex_axial_10N", + "22Hex_axial_20N", + "22Hex_cylinder_0.05N", + "22Hex_cylinder_0.1N", + "22Hex_s_0.05N", + "22Hex_s_0.1N"}; + + std::unordered_map shuffleToNumOfOccur; + +public: + ReducedPatternSimulator( + const std::vector &optimizationResults); + void simulate(std::shared_ptr &tileIntoSurface, + std::vector &optimizationResults, + PatternGeometry &reducedPattern); + /* + * centers the surface + * scales it such that its average base triangle size matches a desired one. This is done in order to match the base triangle on which the reduced pattern was optimized on + * */ + void preprocessSurface(const double &desiredBaseTriangleHeight, VCGPolyMesh &surface); + static void applyOptimizationResults_elements( + const ReducedPatternOptimization::Results &reducedPattern_optimizationResults, + const std::shared_ptr &pTiledReducedPattern_simulationMesh, + const std::vector &elementIndices); + static void applyOptimizationResults_innerHexagon( + const ReducedPatternOptimization::Results &reducedPattern_optimizationResults, + const vcg::Triangle3 &patternBaseTriangle, + PatternGeometry &reducedPattern); + static std::vector loadOptimizationResults( + const std::filesystem::path &optimizationResultsFolderPath); + +private: + std::shared_ptr pTileIntoSurface; + std::vector fullPatterns; + static void createSimulationJobs(const Vector6d &externalForce, + SimulationJob &job_fullPattern, + SimulationJob &job_reducedPattern); + void createGuiMenu(); + std::pair getPickedVertices( + const std::pair &selection) const; + void reset(); + void saveJobs(const filesystem::__cxx11::path &outputFolderPath); + void saveResults(const std::string &outputFolderPath, + SimulationResults &fullResults, + SimulationResults &reducedResults); + void createTiledSimulationMeshes(std::vector &fullPatterns, + std::vector &reducedPatterns); + using OptimizationResultsIndex = size_t; + void createTiledSimulationMeshes( + const std::shared_ptr &pTileIntoSurface, + std::vector &optimizationResults, + const std::vector &perSurfaceFacePatternIndices); + void createTiledSimulationMeshes( + const std::shared_ptr &tileIntoSurface, + std::vector &optimizationResults, + const std::vector &perSurfaceFacePatterns, + std::vector &tileIntoEdgeToTiledFullPattern); + void loadScenario(const std::string &scenarioLabel); + std::pair runSimulations( + const std::string &scenarioLabel, + const bool &shouldRerunFullPatternSimulation = false, + const bool &showReducedPatternResultsBeforeRunningDRM = false); + void runShufflingEvaluation( + const std::vector &optimizationResults, + const std::vector &tileIntoEiToTiledFullVi); + void shuffleReducedPatterns( + std::shared_ptr &pSimulationMesh_tiledShuffledReduced, + std::unordered_map &fullToShuffledReduced); + double computeDistance( + const SimulationResults &resultsA, + const SimulationResults &resultsB, + const std::unordered_map &resultAToResultsBViMap) const; + std::vector computeDistancesPerScenario( + const std::vector &scenarioLabels = std::vector(), + const bool &shouldRerunDRMSimulation = false, + const bool &shouldDraw = false); + double computeDistance(const std::string &scenarioLabel, + const bool &shouldRerunDRMSimulation, + const bool &shouldDraw); + void tileReducedPatterns( + const std::vector &optimizationResults, + const std::vector &perSurfaceFacePatternIndices, + std::shared_ptr &pSimulationMesh_tiledReduced, + std::vector &tileIntoEdgeToTiledReduced); + void shuffleReducedPatterns( + const std::vector &optimizationResults, + const std::vector &tileIntoEdgeToTiledFullPattern, + const std::vector &perSurfaceFacePatternIndices); + void runWeightEvaluation(const std::filesystem::path &optimizationResultsFolderPath); + void constructViMaps(const std::vector &tileIntoEdgeToTiledFullPattern, + const std::vector &tileIntoEdgeToTiledReducedPattern); + std::vector computePerSurfaceFacePatternsIndices( + const std::vector &optimizationResults) const; + void tileFullPatterns(std::vector &fullPatterns, + const std::vector &perSurfaceFacePatternIndices, + std::shared_ptr &pTiledFullPattern_simulationMesh, + std::vector &tileIntoEdgeToTiledFullVi); + void fillFullPatterns(std::vector &optimizationResults); + void computeLabels(const std::vector &optimizationResults, + const std::vector &perSurfaceFacePatternIndex); + void createShufflings(std::vector &optimizationResults); + void generateRandomSimulationScenario(const std::array &randomScenarioParameters); + std::tuple + runAllSimulations(const string &jobLabel, + const bool &shouldRerunDRMFullPatternSimulation = false, + const bool &showReducedPatternResultsBeforeDRM = true); + void reportDistances(const std::tuple &simulationResults, + const string &csvFilePath = {}); +}; + +#endif // REDUCEDPATTERNSIMULATOR_HPP