summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/break_edge.c78
-rw-r--r--src/fracture.c43
-rw-r--r--src/fracture.h1
-rw-r--r--src/net.c6
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;
diff --git a/src/net.c b/src/net.c
index 9f6965b..f3aeda9 100644
--- a/src/net.c
+++ b/src/net.c
@@ -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);