diff --git a/CMakeLists.txt b/CMakeLists.txt index 0f04512..3bc7eec 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -54,7 +54,6 @@ download_project(PROJ TBB option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF) add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR}) link_directories(${TBB_BINARY_DIR}) -message(${TBB_BINARY_DIR}) ###Eigen 3 NOTE: Eigen is required on the system the code is ran @@ -77,7 +76,7 @@ target_include_directories(${PROJECT_NAME} ) if(${USE_POLYSCOPE}) - target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot polyscope glad ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread) + target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot polyscope glad ThreedBeamFEA tbb pthread) else() target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread) endif() diff --git a/csvfile.hpp b/csvfile.hpp index bf17814..7c9dffd 100755 --- a/csvfile.hpp +++ b/csvfile.hpp @@ -32,14 +32,10 @@ public: fs_.clear(std::cout.rdstate()); fs_.basic_ios::rdbuf(std::cout.rdbuf()); } else { - if (!std::filesystem::exists( - std::filesystem::path("../OptimizationResults") - .append("statistics.csv"))) { - std::ofstream outfile(std::filesystem::path("../OptimizationResults") - .append("statistics.csv") - .string()); - outfile.close(); - } + if (!std::filesystem::exists(std::filesystem::path(filename))) { + std::ofstream outfile(filename); + outfile.close(); + } overwrite ? fs_.open(filename, std::ios::trunc) : fs_.open(filename, std::ios::app); } diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index abfbad3..03af7fe 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -1030,11 +1030,6 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( } pMesh->previousTotalResidualForcesNorm = pMesh->totalResidualForcesNorm; pMesh->totalResidualForcesNorm = totalResidualForcesNorm; - if (mSettings.beVerbose) { - if (minTotalResidualForcesNorm > pMesh->totalResidualForcesNorm) { - minTotalResidualForcesNorm = pMesh->totalResidualForcesNorm; - } - } pMesh->averageResidualForcesNorm = totalResidualForcesNorm / pMesh->VN(); // pMesh->averageResidualForcesNorm = sumOfResidualForces.norm() / pMesh->VN(); @@ -2277,6 +2272,17 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr pMesh->totalResidualForcesNorm) { + minTotalResidualForcesNorm = pMesh->totalResidualForcesNorm; + static int lastSavedStep = mCurrentSimulationStep; + if (mSettings.saveIntermediateBestStates.has_value() + && mSettings.saveIntermediateBestStates.value() + && mCurrentSimulationStep - lastSavedStep > 20000) { + lastSavedStep = mCurrentSimulationStep; + computeResults(pJob).save("./IntermediateResults_" + pJob->getLabel()); + } + } + //normalized sum of displacements // double sumOfDisplacementsNorm = 0; // for (size_t vi = 0; vi < pMesh->VN(); vi++) { @@ -2693,6 +2699,16 @@ void DRMSimulationModel::Settings::save(const std::filesystem::path &folderPath) json[jsonLabels.linearGuessForceScaleFactor] = linearGuessForceScaleFactor.value(); } + // if (intermediateResultsSaveStep.has_value()) { + // json[GET_VARIABLE_NAME(intermediateResultsSaveStep)] = intermediateResultsSaveStep.value(); + // } + + if (saveIntermediateBestStates.has_value()) { + json[GET_VARIABLE_NAME(saveIntermediateBestStates)] = saveIntermediateBestStates.value() + ? "true" + : "false"; + } + const std::string jsonFilename = "drmSettings.json"; std::ofstream jsonFile(std::filesystem::path(folderPath).append(jsonFilename)); jsonFile << json; @@ -2756,5 +2772,16 @@ bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePat linearGuessForceScaleFactor = json.at(jsonLabels.linearGuessForceScaleFactor); } + // if (json.contains(GET_VARIABLE_NAME(intermediateResultsSaveStep))) { + // intermediateResultsSaveStep = json.at(GET_VARIABLE_NAME(intermediateResultsSaveStep)); + // } + + if (json.contains(GET_VARIABLE_NAME(saveIntermediateBestStates))) { + saveIntermediateBestStates = json.at(GET_VARIABLE_NAME(saveIntermediateBestStates)) + == "true" + ? true + : false; + } + return true; } diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index 905c465..31d31b0 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -45,6 +45,8 @@ public: std::optional totalTranslationalKineticEnergyThreshold; std::optional averageResidualForcesCriterionThreshold; std::optional linearGuessForceScaleFactor; + // std::optional intermediateResultsSaveStep; + std::optional saveIntermediateBestStates; Settings() {} void save(const std::filesystem::path &folderPath) const; bool load(const std::filesystem::path &filePath); diff --git a/simulation_structs.hpp b/simulation_structs.hpp index 9d3b1df..344779d 100755 --- a/simulation_structs.hpp +++ b/simulation_structs.hpp @@ -195,6 +195,7 @@ public: bool load(const std::string &jsonFilename, const bool &shouldLoadMesh = true) { + std::cout << jsonFilename << std::endl; label = "empty_job"; constrainedVertices.clear(); nodalExternalForces.clear(); @@ -256,8 +257,8 @@ public: if (json.contains(jsonLabels.nodalForces)) { auto f = json[jsonLabels.nodalForces].get>>(); - for (const auto &forces : f) { - nodalExternalForces[forces.first] = Vector6d(forces.second); + for (const auto &force : f) { + nodalExternalForces[force.first] = Vector6d(force.second); } if (beVerbose) { std::cout << "Loaded forces. Number of forces found:" << nodalExternalForces.size() @@ -345,6 +346,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m } std::ofstream jsonFile(jsonFilename); jsonFile << json; + jsonFile.close(); // std::cout << "Saved simulation job as:" << jsonFilename << std::endl; return returnValue; @@ -518,9 +520,12 @@ struct SimulationResults const std::string filename(getLabel() + "_displacements.eigenBin"); Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements); - Eigen::write_binary(std::filesystem::path(outputFolderPath).append(filename).string(), m); + const std::filesystem::path resultsFolderPath( + std::filesystem::path(outputFolder).append("Results")); + std::filesystem::create_directories(resultsFolderPath); + Eigen::write_binary(std::filesystem::path(resultsFolderPath).append(filename).string(), m); - saveDeformedModel(outputFolderPath.string()); + saveDeformedModel(resultsFolderPath.string()); } // The comparison of the results happens comparing the 6-dof nodal diff --git a/topologyenumerator.cpp b/topologyenumerator.cpp index 19ced43..4bad741 100755 --- a/topologyenumerator.cpp +++ b/topologyenumerator.cpp @@ -1,11 +1,12 @@ #include "topologyenumerator.hpp" #include +#include #include #include #include #include -const bool debugIsOn{true}; +const bool debugIsOn{false}; const bool savePlyFiles{true}; // size_t binomialCoefficient(size_t n, size_t m) { @@ -123,7 +124,7 @@ void TopologyEnumerator::computeValidPatterns(const std::vector &reduced // statistics.numberOfValidEdges = validEdges.size(); // Find pairs of intersecting edges - std::unordered_map> intersectingEdges + const std::unordered_map> intersectingEdges = patternAllValidEdges.getIntersectingEdges(statistics.numberOfIntersectingEdgePairs); if (debugIsOn) { auto intersectingEdgesPath = std::filesystem::path(resultsPath) @@ -170,19 +171,19 @@ void TopologyEnumerator::computeValidPatterns(const std::vector &reduced // std::filesystem::path(resultsPath).append("patterns.patt")); // } if (numberOfDesiredEdges == -1) { - for (size_t numberOfEdges = 2; numberOfEdges < validEdges.size(); numberOfEdges++) { + for (size_t numberOfEdges = 2; numberOfEdges <= validEdges.size(); numberOfEdges++) { std::cout << "Computing " + setupString << " with " << numberOfEdges << " edges." << std::endl; auto perEdgeResultPath = std::filesystem::path(resultsPath) .append(std::to_string(numberOfEdges)); if (std::filesystem::exists(perEdgeResultPath)) { - if (debugIsOn) { - std::filesystem::remove_all(perEdgeResultPath); + // if (debugIsOn) { + std::filesystem::remove_all(perEdgeResultPath); - } else { - continue; - } + // } else { + // continue; + // } } std::filesystem::create_directory(perEdgeResultPath); computeValidPatterns(numberOfNodesPerSlot, @@ -434,7 +435,7 @@ void TopologyEnumerator::computeValidPatterns( // std::string previousPatternBinaryRepresentation(validEdges.size(),'0'); size_t patternIndex = 0; bool validPatternsExist = false; - const bool exportTilledPattern = false; + const bool exportTilledPattern = debugIsOn; const bool saveCompressedFormat = false; do { patternIndex++; @@ -479,6 +480,7 @@ void TopologyEnumerator::computeValidPatterns( const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold = patternGeometry.hasAngleSmallerThanThreshold(numberOfNodesPerSlot, 15); if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) { + statistics.numberOfPatternsViolatingAngleThreshold++; if (debugIsOn /*|| savePlyFiles*/) { if (savePlyFiles) { exportPattern(std::filesystem::path(resultsPath) @@ -494,6 +496,7 @@ void TopologyEnumerator::computeValidPatterns( const bool tiledPatternHasNodeWithValenceGreaterThanDesired = patternGeometry.hasValenceGreaterThan(numberOfNodesPerSlot, 6); if (tiledPatternHasNodeWithValenceGreaterThanDesired) { + statistics.numberOfPatternsViolatingValenceThreshold++; if (debugIsOn) { if (savePlyFiles) { auto highValencePath = std::filesystem::path(resultsPath) @@ -510,8 +513,82 @@ void TopologyEnumerator::computeValidPatterns( numberOfNodesPerSlot); // marks the nodes with valence>=1 // Create the tiled geometry of the pattern const bool hasFloatingComponents = !patternGeometry.isFullyConnectedWhenTiled(); - FlatPatternTopology topology(numberOfNodesPerSlot, patternEdges); - const bool hasArticulationPoints = topology.containsArticulationPoints(); + + PatternGeometry fanPatternGeometry = PatternGeometry::createFan(patternGeometry); + const int interfaceNodeVi = 3; + std::vector connectedEdges; + vcg::edge::VEStarVE(&fanPatternGeometry.vert[interfaceNodeVi], connectedEdges); + if (!connectedEdges.empty()) { + for (int i = 1; i < 6; i++) { + vcg::tri::Allocator::AddEdge(fanPatternGeometry, + interfaceNodeVi + + (i - 1) * patternGeometry.VN(), + interfaceNodeVi + + i * patternGeometry.VN()); + } + } + vcg::tri::Clean::MergeCloseVertex(fanPatternGeometry, 0.0000005); + vcg::tri::Allocator::CompactEveryVector(fanPatternGeometry); + vcg::tri::UpdateTopology::VertexEdge(fanPatternGeometry); + vcg::tri::UpdateTopology::EdgeEdge(fanPatternGeometry); + // for (PatternGeometry::VertexType &v : tilledPatternGeometry.vert) { + // std::vector connectedEdges; + // vcg::edge::VEStarVE(&v, connectedEdges); + // if (connectedEdges.size() == 1) { + // vcg::tri::Allocator::DeleteVertex(tilledPatternGeometry, v); + // vcg::tri::Allocator::DeleteEdge(tilledPatternGeometry, + // *connectedEdges[0]); + // } + // } + // // vcg::tri::Allocator::CompactEveryVector(tilledPatternGeometry); + // fanPatternGeometry.updateEigenEdgeAndVertices(); + + BoostGraph fanPatternGraph(fanPatternGeometry.VN()); + // std::cout << "Edges:"; + for (const PatternGeometry::EdgeType &e : fanPatternGeometry.edge) { + if (e.IsD() || e.cV(0)->IsD() || e.cV(1)->IsD()) { + continue; + } + const int vi0 = fanPatternGeometry.getIndex(e.cV(0)); + const int vi1 = fanPatternGeometry.getIndex(e.cV(1)); + boost::add_edge(vi0, vi1, fanPatternGraph); + // std::cout << vi0 << "," << vi1 << " "; + } + // std::cout << std::endl; + + std::vector articulationPoints; + boost::articulation_points(fanPatternGraph, std::back_inserter(articulationPoints)); + const bool hasArticulationPoints = !articulationPoints.empty(); + // if (hasArticulationPoints /*&& !patternContainsIntersectingEdges + // && !tiledPatternHasDanglingEdges && !hasFloatingComponents + // && !tiledPatternHasNodeWithValenceGreaterThanDesired + // && !tiledPatternHasEdgesWithAngleSmallerThanThreshold*/) { + // for (PatternGeometry::VertexType &v : patternGeometry.vert) { + // v.C() = vcg::Color4b::Yellow; + // } + // // std::cout << "AP:"; + // for (const int articulationPointVi : articulationPoints) { + // if (articulationPointVi >= patternGeometry.VN()) { + // continue; + // } + // // std::cout << articulationPointVi << " "; + // patternGeometry.vert[articulationPointVi].C() = vcg::Color4b::Red; + // } + // PatternGeometry tilledPatternGeometry = PatternGeometry::createTile(patternGeometry); + // // std::cout << std::endl; + // std::vector fanVertexColors(tilledPatternGeometry.VN(), glm::vec3(0, 0, 1)); + // for (const PatternGeometry::VertexType &v : tilledPatternGeometry.vert) { + // const auto vColor = glm::vec3(v.cC()[0] / 255, v.cC()[1] / 255, v.cC()[2] / 255); + // const auto vi = tilledPatternGeometry.getIndex(v); + // fanVertexColors[vi] = vColor; + // } + // // tilledPatternGeometry.updateEigenEdgeAndVertices(); + // // tilledPatternGeometry.registerForDrawing() + // // ->addNodeColorQuantity("ap_tilled", fanVertexColors) + // // ->setEnabled(true); + // // polyscope::show(); + // // tilledPatternGeometry.unregister(); + // } // duplicated here // Check dangling edges with vcg method // const bool vcg_tiledPatternHasDangling = @@ -618,7 +695,7 @@ void TopologyEnumerator::computeValidPatterns( // if(patternName=='2055'){ // PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( // patternGeometry); // the marked nodes of hasDanglingEdges are - // tiledPatternGeometry.registerForDrawing(); + // tiledPatternGeometry.registerForDrawing(std::array{0, 0, 1}); // polyscope::show(); // tiledPatternGeometry.unregister(); // } @@ -656,6 +733,7 @@ void TopologyEnumerator::computeValidPatterns( } } } + std::vector TopologyEnumerator::getCoincideEdges( const std::vector &edgeNodeIndices) const { diff --git a/topologyenumerator.hpp b/topologyenumerator.hpp index 7871ce6..0de1444 100755 --- a/topologyenumerator.hpp +++ b/topologyenumerator.hpp @@ -62,6 +62,8 @@ class TopologyEnumerator { size_t numberOfPatternsWithMoreThanASingleCC{0}; size_t numberOfPatternsWithADanglingEdgeOrNode{0}; size_t numberOfPatternsWithArticulationPoints{0}; + size_t numberOfPatternsViolatingAngleThreshold{0}; + size_t numberOfPatternsViolatingValenceThreshold{0}; size_t numberOfValidPatterns{0}; nlohmann::json convertToJson() const { @@ -79,6 +81,8 @@ class TopologyEnumerator { json["numPatternsWithDangling"] = numberOfPatternsWithADanglingEdgeOrNode; json["numPatternsWithArticulationPoints"] = numberOfPatternsWithArticulationPoints; json["numValidPatterns"] = numberOfValidPatterns; + json["violatingAngle"] = numberOfPatternsViolatingAngleThreshold; + json["violatingValence"] = numberOfPatternsViolatingValenceThreshold; return json; } @@ -141,6 +145,8 @@ class TopologyEnumerator { numberOfPatternsWithADanglingEdgeOrNode = 0; numberOfPatternsWithArticulationPoints = 0; numberOfValidPatterns = 0; + numberOfPatternsViolatingAngleThreshold = 0; + numberOfPatternsViolatingValenceThreshold = 0; } }; diff --git a/trianglepatterngeometry.cpp b/trianglepatterngeometry.cpp index 921a768..b8f4b22 100755 --- a/trianglepatterngeometry.cpp +++ b/trianglepatterngeometry.cpp @@ -116,6 +116,7 @@ PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern) { // tile.vert[vi].C() = vcg::Color4b::Blue; // } + tile.setLabel("tilled_" + pattern.getLabel()); tile.updateEigenEdgeAndVertices(); return tile; } @@ -674,23 +675,22 @@ bool PatternGeometry::hasIntersectingEdges( bool containsIntersectingEdges = false; for (size_t validEdgeIndex = 0; validEdgeIndex < patternBinaryRepresentation.size(); validEdgeIndex++) { - if (patternBinaryRepresentation[validEdgeIndex] == '1' && - intersectingEdges.count(validEdgeIndex) != 0) { - for (auto edgeIndexIt = - intersectingEdges.find(validEdgeIndex)->second.begin(); - edgeIndexIt != intersectingEdges.find(validEdgeIndex)->second.end(); - edgeIndexIt++) { - if (patternBinaryRepresentation[*edgeIndexIt] == '1') { - containsIntersectingEdges = true; - // statistics.numberOfIntersectingEdgesOverAllPatterns++; - break; // should be uncommented in order to improve - // performance - } + if (patternBinaryRepresentation[validEdgeIndex] == '1' + && intersectingEdges.contains(validEdgeIndex)) { + for (auto edgeIndexIt = intersectingEdges.find(validEdgeIndex)->second.begin(); + edgeIndexIt != intersectingEdges.find(validEdgeIndex)->second.end(); + edgeIndexIt++) { + if (patternBinaryRepresentation[*edgeIndexIt] == '1') { + containsIntersectingEdges = true; + // statistics.numberOfIntersectingEdgesOverAllPatterns++; + break; // should be uncommented in order to improve + // performance + } + } + if (containsIntersectingEdges) { + break; // should be uncommented in order to improve performance + } } - if (containsIntersectingEdges) { - break; // should be uncommented in order to improve performance - } - } } return containsIntersectingEdges; diff --git a/trianglepattterntopology.cpp b/trianglepattterntopology.cpp index 9686396..b9a863a 100755 --- a/trianglepattterntopology.cpp +++ b/trianglepattterntopology.cpp @@ -5,64 +5,71 @@ FlatPatternTopology::FlatPatternTopology() {} -FlatPatternTopology::FlatPatternTopology( - const std::vector &numberOfNodesPerSlot, - const std::vector &edges) - : numberOfNodesPerSlot(numberOfNodesPerSlot) { - pattern = BoostGraph(std::accumulate(numberOfNodesPerSlot.begin(), - numberOfNodesPerSlot.end(), 0)); - for (const vcg::Point2i &e : edges) { - boost::add_edge(e[0], e[1], pattern); - } +FlatPatternTopology::FlatPatternTopology(const std::vector &numberOfNodesPerSlot, + const std::vector &edges) + : numberOfNodesPerSlot(numberOfNodesPerSlot) +{ + pattern = BoostGraph( + std::accumulate(numberOfNodesPerSlot.begin(), numberOfNodesPerSlot.end(), 0)); + isAdjacentTo.resize(boost::num_vertices(pattern), + std::vector(boost::num_vertices(pattern), false)); + for (const vcg::Point2i &e : edges) { + boost::add_edge(e[0], e[1], pattern); + isAdjacentTo[e[0]][e[1]] = true; + isAdjacentTo[e[1]][e[0]] = true; + } - constructNodeToSlotMap(); - constructSlotToNodeMap(); - constructCorresponginNodeMap(); + constructNodeToSlotMap(); + constructSlotToNodeMap(); + constructCorresponginNodeMap(); } -bool FlatPatternTopology::containsArticulationPoints() const { - assert(numberOfNodesPerSlot.size() == 7 && numberOfNodesPerSlot[4] < 2 && - !nodeToSlot.empty() && !correspondingNode.empty()); - BoostGraph copyOfPattern(pattern); - // std::cout << std::endl; - std::vector componentsBefore(boost::num_vertices(copyOfPattern)); - size_t num_components = - boost::connected_components(copyOfPattern, &componentsBefore[0]); - // std::cout << "Number of cc before:" << num_components << std::endl; - // printGraph(copyOfPattern); +bool FlatPatternTopology::containsArticulationPoints(std::vector &articulationPointsVi) const +{ + assert(numberOfNodesPerSlot.size() == 7 && numberOfNodesPerSlot[4] < 2 && !nodeToSlot.empty() + && !correspondingNode.empty()); + BoostGraph copyOfPattern(pattern); + // std::cout << std::endl; + // std::vector componentsBefore(boost::num_vertices(copyOfPattern)); + // size_t num_components = boost::connected_components(copyOfPattern, &componentsBefore[0]); + // std::cout << "Number of cc before:" << num_components << std::endl; + // printGraph(copyOfPattern); - copyOfPattern = constructRotationallySymmetricPattern( - copyOfPattern, slotToNode, nodeToSlot, correspondingNode); - // // Remove edges connected to the bottom edge node - // assert(slotToNode.find(4) != slotToNode.end()); - // std::unordered_set bottomEdgeNodeSet = - // (slotToNode.find(4))->second; - // size_t bottomEdgeNodeIndex = *bottomEdgeNodeSet.begin(); - // boost::clear_vertex(bottomEdgeNodeIndex, copyOfPattern); + copyOfPattern = constructRotationallySymmetricPattern(copyOfPattern, + slotToNode, + nodeToSlot, + correspondingNode); + // // Remove edges connected to the bottom edge node + // assert(slotToNode.find(4) != slotToNode.end()); + // std::unordered_set bottomEdgeNodeSet = + // (slotToNode.find(4))->second; + // size_t bottomEdgeNodeIndex = *bottomEdgeNodeSet.begin(); + // boost::clear_vertex(bottomEdgeNodeIndex, copyOfPattern); - std::vector componentsAfter(boost::num_vertices(copyOfPattern)); - num_components = - boost::connected_components(copyOfPattern, &componentsAfter[0]); - // std::cout << "Number of cc after:" << num_components << std::endl; - // printGraph(copyOfPattern); + // std::vector componentsAfter(boost::num_vertices(copyOfPattern)); + // num_components = boost::connected_components(copyOfPattern, &componentsAfter[0]); + // std::cout << "Number of cc after:" << num_components << std::endl; + // printGraph(copyOfPattern); - // Compute articulation points on the edited graph - std::vector articulationPoints; - boost::articulation_points(copyOfPattern, - std::back_inserter(articulationPoints)); - // std::cout << "Found " << articulationPoints.size() - // << " articulation points.\n"; - // size_t numberOfNonValidArticulationPoints = 0; - for (std::size_t i = 0; i < articulationPoints.size(); ++i) { + // Compute articulation points on the edited graph + std::vector articulationPoints; + boost::articulation_points(copyOfPattern, std::back_inserter(articulationPoints)); + for (const auto apVi : articulationPoints) { + articulationPointsVi.push_back(static_cast(apVi)); + } + // std::cout << "Found " << articulationPoints.size() + // << " articulation points.\n"; + // size_t numberOfNonValidArticulationPoints = 0; + // for (std::size_t i = 0; i < articulationPoints.size(); ++i) { // std::cout << articulationPoints[i] << std::endl; // if (boost::out_degree(articulationPoints[i], copyOfPattern) < 3) { // numberOfNonValidArticulationPoints++; // } - } - // if (numberOfNonValidArticulationPoints == articulationPoints.size()) { - // return false; - // } - return !articulationPoints.empty(); + // } + // if (numberOfNonValidArticulationPoints == articulationPoints.size()) { + // return false; + // } + return !articulationPoints.empty(); } void FlatPatternTopology::constructNodeToSlotMap( @@ -139,6 +146,29 @@ void FlatPatternTopology::constructCorresponginNodeMap() { } } +bool FlatPatternTopology::pathExists(int src, int dest) const +{ + const int N = boost::num_vertices(pattern); + std::vector visited(N, false); + visited[src] = true; + std::stack next; + next.push(src); + + while (!next.empty()) { + int cv = next.top(); + next.pop(); + + for (int nv = 0; nv < N; ++nv) { + if (!visited[nv] && isAdjacentTo[cv][nv] == 1) { + visited[nv] = true; + next.push(nv); + } + } + } + // dest was reached from src? + return visited[dest]; +} + /* * In this function I create an "extended" pattern of the one in the base * triangle by: @@ -152,147 +182,179 @@ BoostGraph FlatPatternTopology::constructRotationallySymmetricPattern( const BoostGraph &pattern, const std::unordered_map> &slotToNodes, const std::unordered_map &nodeToSlot, - const std::unordered_map &correspondingNode) { - BoostGraph rotationallySymmetricPattern(pattern); + const std::unordered_map &correspondingNode) const +{ + BoostGraph rotationallySymmetricPattern(pattern); - boost::graph_traits::out_edge_iterator ei, ei_end; - // Copy edges that lay on the left edge to the right edge - const auto slot3NodesPairIt = slotToNodes.find(3); - if (slot3NodesPairIt != slotToNodes.end()) { - for (const size_t &nodeIndex : slot3NodesPairIt->second) { - for (boost::tie(ei, ei_end) = boost::out_edges(nodeIndex, pattern); - ei != ei_end; ++ei) { - auto vt = boost::target(*ei, pattern); - const auto vtNodeSlotPairIt = nodeToSlot.find(vt); - assert(vtNodeSlotPairIt != nodeToSlot.end()); - const size_t vtSlot = vtNodeSlotPairIt->second; - if (vtSlot == 3 || vtSlot == 1 || vtSlot == 0) { - // Connect the corresponding nodes on the opposite edge - auto correspondingNodeIndexIt = correspondingNode.find(nodeIndex); - assert(correspondingNodeIndexIt != correspondingNode.end()); - auto correspondingVtIt = correspondingNode.find(vt); - assert(correspondingVtIt != correspondingNode.end() || vtSlot == 0); - const size_t &correspondingNodeIndex = - correspondingNodeIndexIt->second; - size_t correspondingVt = 0; - if (correspondingVtIt != correspondingNode.end()) { - correspondingVt = correspondingVtIt->second; - } - boost::add_edge(correspondingNodeIndex, correspondingVt, - rotationallySymmetricPattern); + // for (const std::pair &correspondingPair : correspondingNode) { + // const auto v0 = boost::vertex(correspondingPair.first, pattern); + // const auto v1 = boost::vertex(correspondingPair.second, pattern); + // if (boost::degree(v0, pattern) != 0 || boost::degree(v1, pattern) != 0) { + // boost::add_edge(v0, v1, rotationallySymmetricPattern); + // } + // } + + boost::graph_traits::out_edge_iterator ei, ei_end; + // Copy edges that lay on the left edge to the right edge + const auto slot3NodesPairIt = slotToNodes.find(3); + if (slot3NodesPairIt != slotToNodes.end()) { + for (const size_t &nodeIndex : slot3NodesPairIt->second) { + for (boost::tie(ei, ei_end) = boost::out_edges(nodeIndex, pattern); ei != ei_end; ++ei) { + auto vt = boost::target(*ei, pattern); + const auto vtNodeSlotPairIt = nodeToSlot.find(vt); + assert(vtNodeSlotPairIt != nodeToSlot.end()); + const size_t vtSlot = vtNodeSlotPairIt->second; + if (vtSlot == 3 || vtSlot == 1 || vtSlot == 0) { + // Connect the corresponding nodes on the opposite edge + auto correspondingNodeIndexIt = correspondingNode.find(nodeIndex); + assert(correspondingNodeIndexIt != correspondingNode.end()); + auto correspondingVtIt = correspondingNode.find(vt); + assert(correspondingVtIt != correspondingNode.end() || vtSlot == 0); + const size_t &correspondingNodeIndex = correspondingNodeIndexIt->second; + size_t correspondingVt = 0; + if (correspondingVtIt != correspondingNode.end()) { + correspondingVt = correspondingVtIt->second; + } + boost::add_edge(correspondingNodeIndex, + correspondingVt, + rotationallySymmetricPattern); + } + } } - } } - } - // Copy edges that lay on the right edge to the left edge - const auto slot5NodesPairIt = slotToNodes.find(5); - if (slot5NodesPairIt != slotToNodes.end()) { - for (const size_t &nodeIndex : slot5NodesPairIt->second) { - for (boost::tie(ei, ei_end) = boost::out_edges(nodeIndex, pattern); - ei != ei_end; ++ei) { - auto vt = boost::target(*ei, pattern); - const auto vtNodeSlotPairIt = nodeToSlot.find(vt); - assert(vtNodeSlotPairIt != nodeToSlot.end()); - const size_t vtSlot = vtNodeSlotPairIt->second; - if (vtSlot == 5 || vtSlot == 2 || vtSlot == 0) { - // Connect the corresponding nodes on the opposite edge - auto correspondingNodeIndexIt = correspondingNode.find(nodeIndex); - assert(correspondingNodeIndexIt != correspondingNode.end()); - auto correspondingVtIt = correspondingNode.find(vt); - assert(correspondingVtIt != correspondingNode.end() || vtSlot == 0); - const size_t &correspondingNodeIndex = - correspondingNodeIndexIt->second; - size_t correspondingVt = 0; - if (correspondingVtIt != correspondingNode.end()) { - correspondingVt = correspondingVtIt->second; - } - boost::add_edge(correspondingNodeIndex, correspondingVt, - rotationallySymmetricPattern); + // // Copy edges that lay on the right edge to the left edge + const auto slot5NodesPairIt = slotToNodes.find(5); + if (slot5NodesPairIt != slotToNodes.end()) { + for (const size_t &nodeIndex : slot5NodesPairIt->second) { + for (boost::tie(ei, ei_end) = boost::out_edges(nodeIndex, pattern); ei != ei_end; ++ei) { + auto vt = boost::target(*ei, pattern); + const auto vtNodeSlotPairIt = nodeToSlot.find(vt); + assert(vtNodeSlotPairIt != nodeToSlot.end()); + const size_t vtSlot = vtNodeSlotPairIt->second; + if (vtSlot == 5 || vtSlot == 2 || vtSlot == 0) { + // Connect the corresponding nodes on the opposite edge + auto correspondingNodeIndexIt = correspondingNode.find(nodeIndex); + assert(correspondingNodeIndexIt != correspondingNode.end()); + auto correspondingVtIt = correspondingNode.find(vt); + assert(correspondingVtIt != correspondingNode.end() || vtSlot == 0); + const size_t &correspondingNodeIndex = correspondingNodeIndexIt->second; + size_t correspondingVt = 0; + if (correspondingVtIt != correspondingNode.end()) { + correspondingVt = correspondingVtIt->second; + } + boost::add_edge(correspondingNodeIndex, + correspondingVt, + rotationallySymmetricPattern); + } + } } - } } - } - // NOTE: The problem with that approach is that I connect !all! "external" - // nodes with each other, which might not be entirely true. If the number of - // cc of the tilled configuration is not 1 this might not label patterns as - // having an articulation node although they might have one. Create set of - // nodes connected with "external" edges - std::unordered_set externallyConnectedNodes; - // Mark the star node as external - const auto slot0NodesPairIt = slotToNodes.find(0); - if (slot0NodesPairIt != slotToNodes.end()) { - externallyConnectedNodes.insert(slot0NodesPairIt->second.begin(), - slot0NodesPairIt->second.end()); - } - // Mark all bottom nodes as external since they are allways connected to the - // south-neighbouring pattern - const auto slot4NodesPairIt = slotToNodes.find(4); - if (slot4NodesPairIt != slotToNodes.end()) { - externallyConnectedNodes.insert(slot4NodesPairIt->second.begin(), - slot4NodesPairIt->second.end()); - } - // Add all slot3 nodes that have a connection to the "inside" - if (slot3NodesPairIt != slotToNodes.end()) { - for (const size_t &nodeIndex : slot3NodesPairIt->second) { - for (boost::tie(ei, ei_end) = boost::out_edges(nodeIndex, pattern); - ei != ei_end; ++ei) { - auto vt = boost::target(*ei, pattern); - const auto vtNodeSlotPairIt = nodeToSlot.find(vt); - assert(vtNodeSlotPairIt != nodeToSlot.end()); - const size_t vtSlot = vtNodeSlotPairIt->second; - if (vtSlot != 3) { - auto correspondingNodePairIt = correspondingNode.find(nodeIndex); - assert(correspondingNodePairIt != correspondingNode.end()); - externallyConnectedNodes.insert(correspondingNodePairIt->second); + // NOTE: The problem with that approach is that I connect !all! "external" + // nodes with each other, which might not be entirely true. If the number of + // cc of the tilled configuration is not 1 this might not label patterns as + // having an articulation node although they might have one. Create set of + // nodes connected with "external" edges + std::unordered_set externallyConnectedNodes; + // Mark the star node as external + const auto slot0NodesPairIt = slotToNodes.find(0); + if (slot0NodesPairIt != slotToNodes.end()) { + externallyConnectedNodes.insert(slot0NodesPairIt->second.begin(), + slot0NodesPairIt->second.end()); + } + // Mark all bottom nodes as external since they are allways connected to the + // south-neighbouring pattern + const auto slot4NodesPairIt = slotToNodes.find(4); + if (slot4NodesPairIt != slotToNodes.end()) { + externallyConnectedNodes.insert(slot4NodesPairIt->second.begin(), + slot4NodesPairIt->second.end()); + } + // Add all slot3 nodes that have a connection to the "inside" + if (slot3NodesPairIt != slotToNodes.end()) { + externallyConnectedNodes.insert(slot3NodesPairIt->second.begin(), + slot3NodesPairIt->second.end()); + for (const size_t &nodeIndex : slot3NodesPairIt->second) { + auto correspondingNodePairIt = correspondingNode.find(nodeIndex); + // for (boost::tie(ei, ei_end) = boost::out_edges(nodeIndex, pattern); ei != ei_end; ++ei) { + // auto vt = boost::target(*ei, pattern); + // const auto vtNodeSlotPairIt = nodeToSlot.find(vt); + // assert(vtNodeSlotPairIt != nodeToSlot.end()); + // const size_t vtSlot = vtNodeSlotPairIt->second; + // if (vtSlot != 3) { + // assert(correspondingNodePairIt != correspondingNode.end()); + // externallyConnectedNodes.insert(correspondingNodePairIt->second); + // // boost::add_edge(correspondingNodePairIt->second, + // // vt, + // // rotationallySymmetricPattern); + // // boost::add_edge(nodeIndex, vt, rotationallySymmetricPattern); + // } + // } + boost::add_edge(correspondingNodePairIt->second, + nodeIndex, + rotationallySymmetricPattern); } - } } - } - // Add all slot5 nodes that have a connection to the "inside" - if (slot5NodesPairIt != slotToNodes.end()) { - for (const size_t &nodeIndex : slot5NodesPairIt->second) { - for (boost::tie(ei, ei_end) = boost::out_edges(nodeIndex, pattern); - ei != ei_end; ++ei) { - auto vt = boost::target(*ei, pattern); - const auto vtNodeSlotPairIt = nodeToSlot.find(vt); - assert(vtNodeSlotPairIt != nodeToSlot.end()); - const size_t vtSlot = vtNodeSlotPairIt->second; - if (vtSlot != 5) { - auto correspondingNodePairIt = correspondingNode.find(nodeIndex); - assert(correspondingNodePairIt != correspondingNode.end()); - externallyConnectedNodes.insert(correspondingNodePairIt->second); + // Add all slot5 nodes that have a connection to the "inside" + if (slot5NodesPairIt != slotToNodes.end()) { + for (const size_t &nodeIndex : slot5NodesPairIt->second) { + auto correspondingNodePairIt = correspondingNode.find(nodeIndex); + // for (boost::tie(ei, ei_end) = boost::out_edges(nodeIndex, pattern); ei != ei_end; ++ei) { + // auto vt = boost::target(*ei, pattern); + // const auto vtNodeSlotPairIt = nodeToSlot.find(vt); + // assert(vtNodeSlotPairIt != nodeToSlot.end()); + // const size_t vtSlot = vtNodeSlotPairIt->second; + // if (vtSlot != 5) { + // assert(correspondingNodePairIt != correspondingNode.end()); + // externallyConnectedNodes.insert(correspondingNodePairIt->second); + // boost::add_edge(correspondingNodePairIt->second, + // vt, + // rotationallySymmetricPattern); + // } + // } + boost::add_edge(correspondingNodePairIt->second, + nodeIndex, + rotationallySymmetricPattern); } - } + externallyConnectedNodes.insert(slot5NodesPairIt->second.begin(), + slot5NodesPairIt->second.end()); } - } - // connecting all is wrong. Maybe I should check whether the external nodes - // are connected via a single node? If so this node is an articulation point. - // I could test this by checking if it filters out only the falsely labeled - // pattern 2367 - const size_t &n = externallyConnectedNodes.size(); - const size_t numberOfExternalEdges = n * (n - 1) / 2; - // Connect all external nodes with each other - for (size_t edgeIndex = 0; edgeIndex < numberOfExternalEdges; edgeIndex++) { - const int sei0 = - n - 2 - - std::floor(std::sqrt(-8 * edgeIndex + 4 * n * (n - 1) - 7) / 2.0 - 0.5); - const int sei1 = edgeIndex + sei0 + 1 - n * (n - 1) / 2 + - (n - sei0) * ((n - sei0) - 1) / 2; - const size_t ni0 = *std::next(externallyConnectedNodes.begin(), sei0); - const size_t ni1 = *std::next(externallyConnectedNodes.begin(), sei1); - boost::add_edge(ni0, ni1, rotationallySymmetricPattern); - } + // connecting all is wrong. Maybe I should check whether the external nodes + // are connected via a single node? If so this node is an articulation point. + // I could test this by checking if it filters out only the falsely labeled + // pattern 2367 + const size_t &n = externallyConnectedNodes.size(); + const size_t numberOfExternalEdges = n * (n - 1) / 2; + // Connect all external nodes with each other + for (size_t edgeIndex = 0; edgeIndex < numberOfExternalEdges; edgeIndex++) { + const int sei0 = n - 2 + - std::floor(std::sqrt(-8 * edgeIndex + 4 * n * (n - 1) - 7) / 2.0 - 0.5); + const int sei1 = edgeIndex + sei0 + 1 - n * (n - 1) / 2 + (n - sei0) * ((n - sei0) - 1) / 2; + const size_t ni0 = *std::next(externallyConnectedNodes.begin(), sei0); + const size_t ni1 = *std::next(externallyConnectedNodes.begin(), sei1); + if (correspondingNode.contains(ni0) || correspondingNode.contains(ni1)) { + if (correspondingNode.contains(ni0) && correspondingNode.contains(ni1) + && pathExists(correspondingNode.at(ni0), correspondingNode.at(ni1))) { + boost::add_edge(ni0, ni1, rotationallySymmetricPattern); + } else if (!correspondingNode.contains(ni0) && correspondingNode.contains(ni1) + && pathExists(ni0, correspondingNode.at(ni1))) { + boost::add_edge(ni0, ni1, rotationallySymmetricPattern); + } else if (correspondingNode.contains(ni0) && !correspondingNode.contains(ni1) + && pathExists(correspondingNode.at(ni0), ni1)) { + boost::add_edge(ni0, ni1, rotationallySymmetricPattern); + } + } + } - return rotationallySymmetricPattern; + return rotationallySymmetricPattern; } -void FlatPatternTopology::printGraph(const BoostGraph &g) const { - boost::graph_traits::edge_iterator ei, ei_end; - for (boost::tie(ei, ei_end) = boost::edges(g); ei != ei_end; ++ei) - std::cout << (char)(boost::source(*ei, g) + 'A') << " -- " - << (char)(boost::target(*ei, g) + 'A'); - std::cout << std::endl; +void FlatPatternTopology::printGraph(const BoostGraph &g) +{ + boost::graph_traits::edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = boost::edges(g); ei != ei_end; ++ei) + std::cout << (char) (boost::source(*ei, g) + 'A') << " -- " + << (char) (boost::target(*ei, g) + 'A'); + std::cout << std::endl; } diff --git a/trianglepattterntopology.hpp b/trianglepattterntopology.hpp index 81844ea..a404739 100755 --- a/trianglepattterntopology.hpp +++ b/trianglepattterntopology.hpp @@ -13,7 +13,7 @@ using vertex_t = boost::graph_traits::vertex_descriptor; class FlatPatternTopology { public: - bool containsArticulationPoints() const; + bool containsArticulationPoints(std::vector &articulationPointsVi) const; FlatPatternTopology(const std::vector &numberOfNodesPerSlot, const std::vector &edges); static void @@ -22,28 +22,30 @@ public: static void constructSlotToNodeMap( const std::unordered_map &nodeToSlot, std::unordered_map> &slotToNode); + static void printGraph(const BoostGraph &g); FlatPatternTopology(); - -private: - BoostGraph pattern; + BoostGraph constructRotationallySymmetricPattern( + const BoostGraph &pattern, + const std::unordered_map> &slotToNodes, + const std::unordered_map &nodeToSlot, + const std::unordered_map &correspondingNode) const; std::vector numberOfNodesPerSlot; std::unordered_map nodeToSlot; std::unordered_map> slotToNode; std::unordered_map correspondingNode; + + private: + std::vector> isAdjacentTo; + BoostGraph pattern; void constructCorresponginNodeMap(); /* * Creates a pattern which is a copy of the input pattern but with added edges * that result * */ - void printGraph(const BoostGraph &g) const; - static BoostGraph constructRotationallySymmetricPattern( - const BoostGraph &pattern, - const std::unordered_map> &slotToNodes, - const std::unordered_map &nodeToSlot, - const std::unordered_map &correspondingNode); void constructNodeToSlotMap(); void constructSlotToNodeMap(); + bool pathExists(int src, int dest) const; }; #endif // FLATPATTTERNTOPOLOGY_HPP diff --git a/utilities.hpp b/utilities.hpp index 1d2dd6e..aac3f0a 100755 --- a/utilities.hpp +++ b/utilities.hpp @@ -9,6 +9,8 @@ #include #include +#define GET_VARIABLE_NAME(Variable) (#Variable) + struct Vector6d : public std::array { Vector6d() { for (size_t i = 0; i < 6; i++) {