diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/break_edge.c | 78 | ||||
| -rw-r--r-- | src/fracture.c | 43 | ||||
| -rw-r--r-- | src/fracture.h | 1 | ||||
| -rw-r--r-- | src/net.c | 6 |
4 files changed, 79 insertions, 49 deletions
diff --git a/src/break_edge.c b/src/break_edge.c index 54eaf34..01570d8 100644 --- a/src/break_edge.c +++ b/src/break_edge.c @@ -1,46 +1,50 @@ #include "fracture.h" -bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) { +bool break_edge(net_t *instance, uint_t edge, cholmod_common *c) { instance->fuses[edge] = true; - unsigned int v1 = instance->graph->ev_break[2 * edge]; - unsigned int v2 = instance->graph->ev_break[2 * edge + 1]; + if (instance->factor != NULL) { + uint_t w1 = instance->graph->ev_break[2 * edge]; + uint_t w2 = instance->graph->ev_break[2 * edge + 1]; + factor_update(instance->factor, w1, w2, c); + } - if (instance->factor != NULL) factor_update(instance->factor, v1, v2, c); + uint_t v1, v2, s1, s2, dv1, dv2, ds1, ds2; - if (instance->graph->boundary != TORUS_BOUND) { - unsigned int w1 = instance->graph->ev[2 * edge]; - unsigned int w2 = instance->graph->ev[2 * edge + 1]; + v1 = instance->graph->ev[2 * edge]; + v2 = instance->graph->ev[2 * edge + 1]; + dv1 = instance->graph->dev[2 * edge]; + dv2 = instance->graph->dev[2 * edge + 1]; - unsigned int tw1 = w1 > w2 ? w1 : w2; - unsigned int tw2 = w1 > w2 ? w2 : w1; + s1 = v1 > v2 ? v1 : v2; + s2 = v1 > v2 ? v2 : v1; + ds1 = dv1 > dv2 ? dv1 : dv2; + ds2 = dv1 > dv2 ? dv2 : dv1; + { int_t *lap_p = (int_t *)instance->adjacency->p; int_t *lap_i = (int_t *)instance->adjacency->i; double *lap_x = (double *)instance->adjacency->x; - for (int i = 0; i < lap_p[tw1 + 1] - lap_p[tw1]; i++) { - if (lap_i[lap_p[tw1] + i] == tw2) - lap_x[lap_p[tw1] + i] = 0; + for (int i = 0; i < lap_p[s1 + 1] - lap_p[s1]; i++) { + if (lap_i[lap_p[s1] + i] == s2) + lap_x[lap_p[s1] + i] = 0; } - for (int i = 0; i < lap_p[tw2 + 1] - lap_p[tw2]; i++) { - if (lap_i[lap_p[tw2] + i] == tw1) - lap_x[lap_p[tw2] + i] = 0; + for (int i = 0; i < lap_p[s2 + 1] - lap_p[s2]; i++) { + if (lap_i[lap_p[s2] + i] == s1) + lap_x[lap_p[s2] + i] = 0; } + } - int old_num_components = instance->num_components; + int_t old_num_components = instance->num_components; - instance->num_components = update_components( - instance->adjacency, instance->marks, old_num_components, (int)tw1, - (int)tw2, 0); - } + instance->num_components = update_components( + instance->adjacency, instance->marks, old_num_components, s1, s2, 0); - if (instance->graph->boundary == TORUS_BOUND) { - unsigned int dw1 = instance->graph->dev[2 * edge]; - unsigned int dw2 = instance->graph->dev[2 * edge + 1]; - if (instance->dual_marks[dw1] == instance->dual_marks[dw2]) { + if (instance->graph->boundary == TORUS_BOUND) { + if (instance->dual_marks[dv1] == instance->dual_marks[dv2]) { int **cycles = (int **)malloc(4*instance->graph->ne * sizeof(int *)); unsigned int num_cycles = find_cycles(instance->graph->ne, instance->fuses, instance->graph->dev, instance->graph->dvei, instance->graph->dve, cycles); @@ -76,26 +80,24 @@ bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) { } free(cycles); } + } - unsigned int tw1 = dw1 > dw2 ? dw1 : dw2; - unsigned int tw2 = dw1 > dw2 ? dw2 : dw1; - - int_t *lap_p = (int_t *)instance->adjacency->p; - int_t *lap_i = (int_t *)instance->adjacency->i; - double *lap_x = (double *)instance->adjacency->x; + { + int_t *lap_p = (int_t *)instance->dual_adjacency->p; + int_t *lap_i = (int_t *)instance->dual_adjacency->i; + double *lap_x = (double *)instance->dual_adjacency->x; - for (int i = 0; i < lap_p[tw1 + 1] - lap_p[tw1]; i++) { - if (lap_i[lap_p[tw1] + i] == tw2) - lap_x[lap_p[tw1] + i] = 1; + for (int i = 0; i < lap_p[ds1 + 1] - lap_p[ds1]; i++) { + if (lap_i[lap_p[ds1] + i] == ds2) + lap_x[lap_p[ds1] + i] = 1; } - for (int i = 0; i < lap_p[tw2 + 1] - lap_p[tw2]; i++) { - if (lap_i[lap_p[tw2] + i] == tw1) - lap_x[lap_p[tw2] + i] = 1; + for (int i = 0; i < lap_p[ds2 + 1] - lap_p[ds2]; i++) { + if (lap_i[lap_p[ds2] + i] == ds1) + lap_x[lap_p[ds2] + i] = 1; } - - set_connected(instance->adjacency, instance->dual_marks, dw1, instance->dual_marks[dw2], -1, 0); } + set_connected(instance->dual_adjacency, instance->dual_marks, dv1, instance->dual_marks[dv2], -1, 0); return true; } diff --git a/src/fracture.c b/src/fracture.c index f38568f..718fc05 100644 --- a/src/fracture.c +++ b/src/fracture.c @@ -10,7 +10,7 @@ int main(int argc, char *argv[]) { uint32_t N; uint_t L; double beta, inf, cutoff, crack_len; - bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network, + bool crack_growth_crit, save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network, save_crit_stress, save_stress_field, save_voltage_field, save_toughness, save_conductivity, save_damage, save_damage_field; bound_t boundary; @@ -32,6 +32,7 @@ int main(int argc, char *argv[]) { cutoff = 1e-9; boundary = FREE_BOUND; lattice = VORONOI_LATTICE; + crack_growth_crit = false; save_data = false; save_cluster_dist = false; use_voltage_boundaries = false; @@ -54,7 +55,7 @@ int main(int argc, char *argv[]) { // get commandline options - while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:")) != -1) { + while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:g")) != -1) { switch (opt) { case 'n': N = atoi(optarg); @@ -111,6 +112,9 @@ int main(int argc, char *argv[]) { exit(EXIT_FAILURE); } break; + case 'g': + crack_growth_crit = true; + break; case 'd': save_damage = true; break; @@ -294,16 +298,31 @@ int main(int argc, char *argv[]) { } } - if (save_crit_stress) crit_stress[i] = data->extern_field[max_pos]; - - if (save_conductivity) conductivity[i] = data->conductivity[max_pos]; - - if (save_damage) damage[max_pos]++; - uint_t av_size = 0; double cur_val = 0; + uint_t notch_edge; + + if (crack_growth_crit) { + for (uint_t i = 0; i < net->graph->ne; i++) { + if (net->fuses[i]) { + notch_edge = i; + break; + } + } + } + for (uint_t j = 0; j < max_pos; j++) { - break_edge(net, data->break_list[j], &c); + uint_t next_broken = data->break_list[j]; + + bool attached_to_notch = net->dual_marks[next_broken] == net->dual_marks[notch_edge]; + bool grew_forward = net->graph->ex[2 * next_broken] > crack_len; + + if (crack_growth_crit && attached_to_notch && grew_forward) { + max_pos = j; + break; + } + + break_edge(net, next_broken, &c); double val = data->extern_field[j]; if (save_cluster_dist) { @@ -319,6 +338,12 @@ int main(int argc, char *argv[]) { } } + if (save_crit_stress) crit_stress[i] = data->extern_field[max_pos]; + + if (save_conductivity) conductivity[i] = data->conductivity[max_pos]; + + if (save_damage) damage[max_pos]++; + if (save_stress_field || save_voltage_field) { double *tmp_voltages = net_voltages(net, &c); if (save_voltage_field) { diff --git a/src/fracture.h b/src/fracture.h index 5e1b3e4..f0a22e8 100644 --- a/src/fracture.h +++ b/src/fracture.h @@ -81,6 +81,7 @@ typedef struct { bool voltage_bound; unsigned int num_components; cholmod_sparse *adjacency; + cholmod_sparse *dual_adjacency; bool debug_stop; } net_t; @@ -57,8 +57,8 @@ net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, b net->voltage_bound = vb; net->boundary_cond = bound_set(g, vb, notch_len, c); - if (g->boundary != TORUS_BOUND) net->adjacency = gen_adjacency(net, false, false, 0, c); - else net->adjacency = gen_adjacency(net, true, false, 0, c); + net->adjacency = gen_adjacency(net, false, false, 0, c); + net->dual_adjacency = gen_adjacency(net, true, false, 0, c); net->marks = (uint_t *)malloc((net->graph->break_dim) * sizeof(uint_t)); assert(net->marks != NULL); @@ -112,6 +112,7 @@ net_t *net_copy(const net_t *net, cholmod_common *c) { memcpy(net_copy->dual_marks, net->dual_marks, dual_marks_size); net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c); + net_copy->dual_adjacency = CHOL_F(copy_sparse)(net->dual_adjacency, c); net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c); net_copy->factor = CHOL_F(copy_factor)(net->factor, c); @@ -123,6 +124,7 @@ void net_free(net_t *net, cholmod_common *c) { free(net->thres); CHOL_F(free_dense)(&(net->boundary_cond), c); CHOL_F(free_sparse)(&(net->adjacency), c); + CHOL_F(free_sparse)(&(net->dual_adjacency), c); CHOL_F(free_factor)(&(net->factor), c); free(net->marks); free(net->dual_marks); |