1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
#include "analysis_tools.hpp"
template <class T>
bool is_shorter(const std::list<T> &l1, const std::list<T> &l2) {
return l1.size() < l2.size();
}
bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>) {
return true;
}
std::list<unsigned> find_minimal_crack(const Graph& G, const network& n) {
Graph Gtmp(n.G.vertices.size());
std::list<unsigned> removed_edges;
class add_tree_edges : public boost::default_dfs_visitor {
public:
Graph& G;
std::list<unsigned>& E;
add_tree_edges(Graph& G, std::list<unsigned>& E) : G(G), E(E) {}
void tree_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
boost::add_edge(boost::source(e, g), boost::target(e, g), g[e], G);
}
void back_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
if (!(boost::edge(boost::source(e, g), boost::target(e, g), G).second)) {
E.push_back(g[e].index);
}
}
};
add_tree_edges ate(Gtmp, removed_edges);
boost::depth_first_search(G, visitor(ate));
class find_cycle : public boost::default_dfs_visitor {
public:
std::list<unsigned>& E;
unsigned end;
struct done{};
find_cycle(std::list<unsigned>& E, unsigned end) : E(E), end(end) {}
void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) {
if (v == end) {
throw done{};
}
}
void examine_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
E.push_back(g[e].index);
}
void finish_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
E.erase(std::find(E.begin(), E.end(), g[e].index));
}
};
std::list<std::list<unsigned>> cycles;
for (auto edge : removed_edges) {
std::list<unsigned> cycle = {edge};
find_cycle vis(cycle, n.G.dual_edges[edge].v[1]);
std::vector<boost::default_color_type> new_color_map(boost::num_vertices(Gtmp));
try {
boost::depth_first_visit(Gtmp, n.G.dual_edges[edge].v[0], vis, boost::make_iterator_property_map(new_color_map.begin(), boost::get(boost::vertex_index, Gtmp), new_color_map[0]));
} catch(find_cycle::done const&) {
cycles.push_back(cycle);
}
}
if (cycles.size() > 1) {
std::list<std::valarray<uint8_t>> bool_cycles;
for (auto cycle : cycles) {
std::valarray<uint8_t> bool_cycle(n.G.edges.size());
for (auto v : cycle) {
bool_cycle[v] = 1;
}
bool_cycles.push_back(bool_cycle);
}
// generate all possible cycles by taking xor of the edge sets of the known cycles
for (auto it1 = bool_cycles.begin(); it1 != std::prev(bool_cycles.end()); it1++) {
for (auto it2 = std::next(it1); it2 != bool_cycles.end(); it2++) {
std::valarray<uint8_t> new_bool_cycle = (*it1) ^ (*it2);
std::list<unsigned> new_cycle;
unsigned pos = 0;
for (uint8_t included : new_bool_cycle) {
if (included) {
new_cycle.push_back(pos);
}
pos++;
}
cycles.push_back(new_cycle);
}
}
// find the cycle representing the crack by counting boundary crossings
for (auto cycle : cycles) {
std::array<unsigned, 2> crossing_count{0,0};
for (auto edge : cycle) {
if (n.G.dual_edges[edge].crossings[0]) {
crossing_count[0]++;
}
if (n.G.dual_edges[edge].crossings[1]) {
crossing_count[1]++;
}
}
if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
return cycle;
}
}
} else {
return cycles.front();
}
exit(5);
}
|