diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index 660e29b..bec098a 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -41,7 +41,7 @@ public: // double residualForcesMovingAverageNormThreshold{1e-8}; double Dtini{0.1}; double xi{0.9969}; - std::optional shouldUseTranslationalKineticEnergyThreshold; + std::optional translationalKineticEnergyThreshold; int gradualForcedDisplacementSteps{50}; // int desiredGradualExternalLoadsSteps{1}; double gamma{0.8}; diff --git a/reducedmodelevaluator.cpp b/reducedmodelevaluator.cpp index 89ee988..97ac713 100644 --- a/reducedmodelevaluator.cpp +++ b/reducedmodelevaluator.cpp @@ -10,9 +10,85 @@ ReducedModelEvaluator::ReducedModelEvaluator() { } +ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( + ReducedModelOptimization::Results &optimizationResult) +{ + const std::filesystem::path tileInto_triMesh_filePath + = "/home/iason/Coding/build/PatternTillingReducedModel/Meshes/" + "instantMeshes_plane_34.ply"; + const std::filesystem::path scenariosDirectoryPath + = "/home/iason/Coding/Projects/Approximating shapes with flat " + "patterns/ReducedModelEvaluator/Scenarios"; + const std::filesystem::path reducedPatternFilePath + = "/home/iason/Coding/Projects/Approximating shapes with flat " + "patterns/ReducedModelOptimization/TestSet/ReducedPatterns/single_reduced.ply"; + const std::filesystem::path fullPatternTessellatedResultsDirectoryPath + = "/home/iason/Coding/build/ReducedModelOptimization/" + "IntermediateResults/TessellatedResults"; -std::vector ReducedModelEvaluator::evaluateReducedModel( - ReducedModelOptimization::Results &optimizationResults) + return evaluateReducedModel(optimizationResult, + tileInto_triMesh_filePath, + scenariosDirectoryPath, + reducedPatternFilePath, + fullPatternTessellatedResultsDirectoryPath); +} + +void ReducedModelEvaluator::printResults(const Results &evaluationResults, + const std::string &resultsLabel, + const std::filesystem::path &resultsOutputDirPath) +{ + //report distance + // csvFile csv_results("", flse); + if (!resultsOutputDirPath.empty()) { + std::filesystem::create_directories(resultsOutputDirPath); + } +#ifdef POLYSCOPE_DEFINED + csvFile csv_results({}, false); +#else + std::filesystem::path csvOutputFilePath = resultsOutputDirPath.empty() + ? "" + : std::filesystem::path(resultsOutputDirPath) + .append("distances_" + resultsLabel + ".csv") + .string(); + csvFile csv_results(csvOutputFilePath, true); +#endif + csv_results << resultsLabel; + csv_results << endrow; +#ifdef POLYSCOPE_DEFINED + csv_results /*<< "Job Label"*/ + << "drm2Reduced" + << "norm_drm2Reduced"; +#else + csv_results << "Job Label" + << "drm2Reduced" + << "norm_drm2Reduced"; + +#endif + csv_results << endrow; + // double sumOfNormalizedFull2Reduced = 0; + for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size(); + jobIndex++) { + const std::string &jobLabel = ReducedModelEvaluator::scenariosTestSetLabels[jobIndex]; + const double &distance_fullDrmToReduced = evaluationResults.distances_drm2reduced[jobIndex]; + const double &distance_normalizedFullDrmToReduced + = evaluationResults.distances_normalizedDrm2reduced[jobIndex]; +#ifdef POLYSCOPE_DEFINED + csv_results /*<< jobLabel*/ << distance_fullDrmToReduced + << distance_normalizedFullDrmToReduced; +#else + csv_results << jobLabel << distance_fullDrmToReduced << distance_normalizedFullDrmToReduced; +#endif + csv_results << endrow; + // sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced; + } +} + +ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( + ReducedModelOptimization::Results &optimizationResult, + const std::filesystem::path &tileInto_triMesh_filePath, + const std::filesystem::path &scenariosDirectoryPath, + const std::filesystem::path &reducedPatternFilePath, + const std::filesystem::path &fullPatternTessellatedResultsDirectoryPath) { // std::shared_ptr pTileIntoSurface = std::make_shared(); // std::filesystem::path tileIntoSurfaceFilePath{ @@ -22,71 +98,31 @@ std::vector ReducedModelEvaluator::evaluateReducedModel( // pTileIntoSurface->load(tileIntoSurfaceFilePath); //Set required file paths - const std::filesystem::path tileInto_triMesh_filename - = "/home/iason/Coding/build/PatternTillingReducedModel/Meshes/" - "instantMeshes_plane_34.ply"; - const std::filesystem::path reducedPatternFilePath - = "/home/iason/Coding/Projects/Approximating shapes with flat " - "patterns/ReducedModelOptimization/TestSet/ReducedPatterns/single_reduced.ply"; - const std::filesystem::path intermediateResultsDirectoryPath - = "/home/iason/Coding/build/ReducedModelOptimization/IntermediateResults"; // const std::filesystem::path drmSettingsFilePath // = "/home/iason/Coding/Projects/Approximating shapes with flat " // "patterns/ReducedModelOptimization/DefaultSettings/DRMSettings/" // "defaultDRMSimulationSettings.json"; DRMSimulationModel::Settings drmSimulationSettings; + drmSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3; // drmSimulationSettings.load(drmSettingsFilePath); - drmSimulationSettings.linearGuessForceScaleFactor = 1; - drmSimulationSettings.debugModeStep = 1000; + drmSimulationSettings.linearGuessForceScaleFactor = 0.5; + drmSimulationSettings.debugModeStep = 10000; drmSimulationSettings.beVerbose = true; + drmSimulationSettings.maxDRMIterations = 1e6; + // drmSimulationSettings.translationalKineticEnergyThreshold = 1e-15; + // drmSimulationSettings.viscousDampingFactor = 1e-2; constexpr bool shouldRerunFullPatternSimulation = false; - const std::vector scenariosTestSetLabels{"22Hex_randomBending0", - "22Hex_randomBending1", - "22Hex_randomBending2", - "22Hex_randomBending3", - "22Hex_randomBending4", - "22Hex_randomBending5", - "22Hex_randomBending6", - "22Hex_randomBending7", - "22Hex_randomBending8", - "22Hex_randomBending9", - "22Hex_randomBending10", - "22Hex_randomBending11", - "22Hex_randomBending12", - "22Hex_randomBending13", - "22Hex_randomBending14", - "22Hex_randomBending15", - "22Hex_randomBending16", - "22Hex_randomBending17", - "22Hex_randomBending18", - "22Hex_randomBending19", - "22Hex_randomBending20", - "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"}; //Load surface std::shared_ptr pTileIntoSurface = [&]() { VCGTriMesh tileInto_triMesh; - const bool surfaceLoadSuccessfull = tileInto_triMesh.load(tileInto_triMesh_filename); + const bool surfaceLoadSuccessfull = tileInto_triMesh.load(tileInto_triMesh_filePath); assert(surfaceLoadSuccessfull); return PolygonalRemeshing::remeshWithPolygons(tileInto_triMesh); }(); - const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResults.baseTriangle.cP(0), - (optimizationResults.baseTriangle.cP(1) - + optimizationResults.baseTriangle.cP( + const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResult.baseTriangle.cP(0), + (optimizationResult.baseTriangle.cP(1) + + optimizationResult.baseTriangle.cP( 2)) / 2); pTileIntoSurface->moveToCenter(); @@ -98,7 +134,7 @@ std::vector ReducedModelEvaluator::evaluateReducedModel( //Tile full pattern into surface std::vector fullPatterns(1); - fullPatterns[0].copy(optimizationResults.baseTriangleFullPattern); + fullPatterns[0].copy(optimizationResult.baseTriangleFullPattern); //// Base triangle pattern might contain dangling vertices.Remove those fullPatterns[0].interfaceNodeIndex = 3; fullPatterns[0].deleteDanglingVertices(); @@ -118,12 +154,13 @@ std::vector ReducedModelEvaluator::evaluateReducedModel( PatternGeometry reducedPattern; reducedPattern.load(reducedPatternFilePath); reducedPattern.deleteDanglingVertices(); + reducedPattern.interfaceNodeIndex = 1; assert(reducedPattern.interfaceNodeIndex == 1); std::vector reducedPatterns(1); reducedPatterns[0].copy(reducedPattern); const auto reducedPatternBaseTriangle = reducedPattern.computeBaseTriangle(); ReducedModelOptimization::Results::applyOptimizationResults_innerHexagon( - optimizationResults, reducedPatternBaseTriangle, reducedPatterns[0]); + optimizationResult, reducedPatternBaseTriangle, reducedPatterns[0]); std::vector> perPatternIndexTiledReducedPatternEdgeIndices; std::vector tileIntoEdgeToTiledReducedVi; @@ -160,12 +197,12 @@ std::vector ReducedModelEvaluator::evaluateReducedModel( 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}); - if (optimizationResults.fullPatternYoungsModulus == 0) { + if (optimizationResult.fullPatternYoungsModulus == 0) { std::cerr << "Full pattern's young modulus not found." << std::endl; std::terminate(); } pTiledFullPattern_simulationMesh->setBeamMaterial(0.3, - optimizationResults.fullPatternYoungsModulus); + optimizationResult.fullPatternYoungsModulus); pTiledFullPattern_simulationMesh->reset(); ////Tessellated reduced pattern simulation mesh @@ -174,20 +211,24 @@ std::vector ReducedModelEvaluator::evaluateReducedModel( const std::vector &tiledPatternElementIndicesForReducedPattern = perPatternIndexTiledReducedPatternEdgeIndices[0]; ReducedModelOptimization::Results::applyOptimizationResults_elements( - optimizationResults, pTiledReducedPattern_simulationMesh); + optimizationResult, pTiledReducedPattern_simulationMesh); pTiledReducedPattern_simulationMesh->reset(); - std::vector distances_drm2reduced(scenariosTestSetLabels.size(), 0); - std::vector distances_normalizedDrm2reduced(scenariosTestSetLabels.size(), 0); + Results evaluationResults; + evaluationResults.distances_drm2reduced; + evaluationResults.distances_normalizedDrm2reduced; for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) { const std::string &jobLabel = scenariosTestSetLabels[jobIndex]; - const std::filesystem::path scenariosDirectoryPath - = "/home/iason/Coding/build/PatternTillingReducedModel/Scenarios/"; const std::filesystem::path tiledReducedPatternJobFilePath = std::filesystem::path(scenariosDirectoryPath) + .append(pTileIntoSurface->getLabel()) .append(jobLabel) - .append("SimulationJobs") - .append("Reduced") + .append("ReducedJob") .append(SimulationJob::jsonDefaultFileName); + if (!std::filesystem::exists(tiledReducedPatternJobFilePath)) { + std::cerr << "Scenario " << jobLabel + << " not found in:" << tiledReducedPatternJobFilePath << std::endl; + continue; + } //set jobs std::shared_ptr pJob_tiledReducedPattern; pJob_tiledReducedPattern = std::make_shared(SimulationJob()); @@ -201,11 +242,10 @@ std::vector ReducedModelEvaluator::evaluateReducedModel( // pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern->getLabel()); // polyscope::show(); //Save reduced job - const std::filesystem::path tesellatedResultsFolderPath - = std::filesystem::path(intermediateResultsDirectoryPath).append("TessellatedResults"); - const std::filesystem::path surfaceFolderPath = std::filesystem::path( - tesellatedResultsFolderPath) - .append(pTileIntoSurface->getLabel()); + + const std::filesystem::path surfaceFolderPath + = std::filesystem::path(fullPatternTessellatedResultsDirectoryPath) + .append(pTileIntoSurface->getLabel()); const std::string scenarioLabel = pJob_tiledFullPattern->getLabel(); const std::filesystem::path scenarioDirectoryPath = std::filesystem::path(surfaceFolderPath) .append(scenarioLabel); @@ -214,21 +254,8 @@ std::vector ReducedModelEvaluator::evaluateReducedModel( std::filesystem::create_directories(reducedJobDirectoryPath); pJob_tiledReducedPattern->save(reducedJobDirectoryPath); //Run scenario - ////Full - - const std::string patternLabel = [&]() { - const std::string patternLabel = optimizationResults.baseTriangleFullPattern.getLabel(); - const int numberOfOccurences = std::count_if(patternLabel.begin(), - patternLabel.end(), - [](char c) { return c == '#'; }); - if (numberOfOccurences == 0) { - return std::to_string(optimizationResults.baseTriangleFullPattern.EN()) + "#" - + optimizationResults.baseTriangleFullPattern.getLabel(); - } else if (numberOfOccurences == 1) { - return optimizationResults.baseTriangleFullPattern.getLabel(); - } - }(); + const std::string patternLabel = optimizationResult.baseTriangleFullPattern.getLabel(); const auto fullResultsFolderPath = std::filesystem::path(scenarioDirectoryPath).append(patternLabel).append("Results"); if (shouldRerunFullPatternSimulation && std::filesystem::exists(fullResultsFolderPath)) { @@ -255,10 +282,13 @@ std::vector ReducedModelEvaluator::evaluateReducedModel( simulationResults_tiledFullPattern_drm.setLabelPrefix("DRM"); } std::filesystem::create_directories(fullResultsFolderPath); - simulationResults_tiledFullPattern_drm.save( - std::filesystem::path(scenarioDirectoryPath).append(patternLabel)); + const std::filesystem::path drmResultsOutputPath + = std::filesystem::path(scenarioDirectoryPath).append(patternLabel); + simulationResults_tiledFullPattern_drm.save(drmResultsOutputPath); if (!simulationResults_tiledFullPattern_drm.converged) { std::cerr << "Full pattern simulation failed." << std::endl; + std::cerr << "Saved failed simulation to:" << drmResultsOutputPath << std::endl; + continue; } LinearSimulationModel linearSimulationModel; @@ -281,37 +311,10 @@ std::vector ReducedModelEvaluator::evaluateReducedModel( } const double distance_normalizedFullDrmToReduced = distance_fullDrmToReduced / distance_fullSumOfAllVerts; - distances_drm2reduced[jobIndex] = distance_fullDrmToReduced; - distances_normalizedDrm2reduced[jobIndex] = distance_normalizedFullDrmToReduced; - } - //#ifndef POLYSCOPE_DEFINED - // return distances_drm2reduced; - //#else - - //report distance - // csvFile csv_results("", flse); - csvFile csv_results({}, false); - // csvFile csv_results(std::filesystem::path(dirPath_thisOptimization) - // .append("results.csv") - // .string(), - // false); - csv_results /*<< "Job Label"*/ - << "drm2Reduced" - << "norm_drm2Reduced"; - csv_results << endrow; - double sumOfNormalizedFull2Reduced = 0; - for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) { - const std::string &jobLabel = scenariosTestSetLabels[jobIndex]; - const double &distance_fullDrmToReduced = distances_drm2reduced[jobIndex]; - const double &distance_normalizedFullDrmToReduced = distances_normalizedDrm2reduced[jobIndex]; - csv_results /*<< jobLabel*/ << distance_fullDrmToReduced - << distance_normalizedFullDrmToReduced; - csv_results << endrow; - sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced; + evaluationResults.distances_drm2reduced[jobIndex] = distance_fullDrmToReduced; + evaluationResults.distances_normalizedDrm2reduced[jobIndex] + = distance_normalizedFullDrmToReduced; } - std::cout << "Average normalized error per scenario:" - << sumOfNormalizedFull2Reduced / scenariosTestSetLabels.size() << std::endl; - return distances_normalizedDrm2reduced; - //#endif + return evaluationResults; } diff --git a/reducedmodelevaluator.hpp b/reducedmodelevaluator.hpp index dd06918..4e5102d 100644 --- a/reducedmodelevaluator.hpp +++ b/reducedmodelevaluator.hpp @@ -6,9 +6,63 @@ class ReducedModelEvaluator { public: + inline static constexpr int NumberOfEvaluationScenarios{36}; + struct Results + { + std::array distances_drm2reduced; + std::array distances_normalizedDrm2reduced; + std::array evaluationScenarioLabels; + }; ReducedModelEvaluator(); - static std::vector evaluateReducedModel( - ReducedModelOptimization::Results &optimizationResults); + static Results evaluateReducedModel( + ReducedModelOptimization::Results &optimizationResult, + const std::filesystem::path &tileInto_triMesh_filePath, + const std::filesystem::path &scenariosDirectoryPath, + const std::filesystem::path &reducedPatternFilePath, + const std::filesystem::path &fullPatternTessellatedResultsDirectoryPath); + static Results evaluateReducedModel(ReducedModelOptimization::Results &optimizationResult); + static void printResults( + const ReducedModelEvaluator::Results &evaluationResults, + const std::string &resultsLabel, + const std::filesystem::path &resultsOutputDirPath = std::filesystem::path()); + + inline static std::array + scenariosTestSetLabels{"22Hex_randomBending0", + "22Hex_randomBending1", + "22Hex_randomBending2", + "22Hex_randomBending3", + "22Hex_randomBending4", + "22Hex_randomBending5", + "22Hex_randomBending6", + "22Hex_randomBending7", + "22Hex_randomBending8", + "22Hex_randomBending9", + "22Hex_randomBending10", + "22Hex_randomBending11", + "22Hex_randomBending12", + "22Hex_randomBending13", + "22Hex_randomBending14", + "22Hex_randomBending15", + "22Hex_randomBending16", + "22Hex_randomBending17", + "22Hex_randomBending18", + "22Hex_randomBending19", + "22Hex_randomBending20", + "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"}; }; #endif // REDUCEDMODELEVALUATOR_HPP diff --git a/trianglepatterngeometry.cpp b/trianglepatterngeometry.cpp index b189371..f8aa2e3 100755 --- a/trianglepatterngeometry.cpp +++ b/trianglepatterngeometry.cpp @@ -154,6 +154,7 @@ bool PatternGeometry::load(const std::filesystem::__cxx11::path &meshFilePath) if (!VCGEdgeMesh::load(meshFilePath)) { return false; } + addNormals(); vcg::tri::UpdateTopology::VertexEdge(*this); baseTriangleHeight = computeBaseTriangleHeight(); baseTriangle = computeBaseTriangle(); diff --git a/trianglepatterngeometry.hpp b/trianglepatterngeometry.hpp index 150912f..3309446 100755 --- a/trianglepatterngeometry.hpp +++ b/trianglepatterngeometry.hpp @@ -40,6 +40,7 @@ private: const int &vi); public: + inline static VectorType DefaultNormal{0.0, 0.0, 1.0}; PatternGeometry(); /*The following function should be a copy constructor with * a const argument but this is impossible due to the @@ -50,11 +51,11 @@ private: void add(const std::vector &vertices); void add(const std::vector &edges); void add(const std::vector &vertices, const std::vector &edges); - void add(const std::vector &numberOfNodesPerSlot, - const std::vector &edges); - static std::vector - constructVertexVector(const std::vector &numberOfNodesPerSlot, - const size_t &fanSize, const double &triangleEdgeSize); + void add(const std::vector &numberOfNodesPerSlot, const std::vector &edges); + static std::vector constructVertexVector( + const std::vector &numberOfNodesPerSlot, + const size_t &fanSize, + const double &triangleEdgeSize); bool hasDanglingEdges(const std::vector &numberOfNodesPerSlot); bool hasValenceGreaterThan(const std::vector &numberOfNodesPerSlot, const size_t &valenceThreshold); @@ -63,12 +64,13 @@ private: static PatternGeometry createTile(PatternGeometry &pattern); double getTriangleEdgeSize() const; bool hasUntiledDanglingEdges(); - std::unordered_map> - getIntersectingEdges(size_t &numberOfIntersectingEdgePairs) const; + std::unordered_map> getIntersectingEdges( + size_t &numberOfIntersectingEdgePairs) const; - static size_t binomialCoefficient(size_t n, size_t m) { - assert(n >= m); - return tgamma(n + 1) / (tgamma(m + 1) * tgamma(n - m + 1)); + static size_t binomialCoefficient(size_t n, size_t m) + { + assert(n >= m); + return tgamma(n + 1) / (tgamma(m + 1) * tgamma(n - m + 1)); } bool isFullyConnectedWhenFanned();