diff options
46 files changed, 5826 insertions, 0 deletions
diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a8276ab --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +bin/* +obj/* +data/* diff --git a/README.md b/README.md new file mode 100644 index 0000000..2c1ec50 --- /dev/null +++ b/README.md @@ -0,0 +1,4 @@ +fracture +======== + +tools for simulating fracture in fuse networks
\ No newline at end of file diff --git a/makefile b/makefile new file mode 100644 index 0000000..eb9fc4d --- /dev/null +++ b/makefile @@ -0,0 +1,31 @@ + +CC = clang +CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -march=native #-fopenmp +LDFLAGS = -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lmetis -lgsl -lprofiler -ltcmalloc + +OBJ = break_data voronoi_bound_ini bin_values correlations beta_scales randfuncs instance get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry fracture_network get_current cracking_ini update_factor update_boundary get_file update_beta gen_voltcurmat ini_network free_network fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory get_conductivity +BIN = fracture current_scaling course_grain_square corr_test homo_voronoi_fracture homo_square_fracture compare_voronoi_fracture + +all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} + +opt: + @echo build options: + @echo "CC = ${CC}" + @echo "CFLAGS = ${CFLAGS}" + @echo "LDFLAGS = ${LDFLAGS}" + +obj/%.o: src/%.c + @echo CC -c -o $@ + @${CC} -c -o $@ $< ${CFLAGS} + +bin/%: obj/%.o ${OBJ:%=obj/%.o} + @echo CC -o $@ + @${CC} -o $@ $< ${OBJ:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS} + +clean: + @echo cleaning: + @echo rm -f ${OBJ:%=obj/%} ${BIN:%=obj/%.o} ${BIN:%=bin/%} + @rm -f ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} + +.PHONY: all clean + diff --git a/makefile_ept b/makefile_ept new file mode 100644 index 0000000..826fde8 --- /dev/null +++ b/makefile_ept @@ -0,0 +1,31 @@ + +CC = clang +CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -I/usr/local/cuda/include -I/home/jaron/.local/include -I/usr/lib/gcc/x86_64-linux-gnu/5/include -march=native -DGPU_BLAS -fopenmp=libiomp5 +LDFLAGS = -L/usr/local/cuda/lib64 -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lmetis -lgsl -lprofiler -ltcmalloc + +OBJ = break_data voronoi_bound_ini bin_values correlations beta_scales randfuncs instance get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry fracture_network get_current cracking_ini update_factor update_boundary get_file update_beta gen_voltcurmat ini_network free_network fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory get_conductivity +BIN = fracture current_scaling course_grain_square corr_test homo_voronoi_fracture + +all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} + +opt: + @echo build options: + @echo "CC = ${CC}" + @echo "CFLAGS = ${CFLAGS}" + @echo "LDFLAGS = ${LDFLAGS}" + +obj/%.o: src/%.c + @echo CC -c -o $@ + @${CC} -c -o $@ $< ${CFLAGS} + +bin/%: obj/%.o ${OBJ:%=obj/%.o} + @echo CC -o $@ + @${CC} -o $@ $< ${OBJ:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS} + +clean: + @echo cleaning: + @echo rm -f ${OBJ:%=obj/%} ${BIN:%=obj/%.o} ${BIN:%=bin/%} + @rm -f ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} + +.PHONY: all clean + diff --git a/makefile_hal b/makefile_hal new file mode 100644 index 0000000..73cb7a8 --- /dev/null +++ b/makefile_hal @@ -0,0 +1,31 @@ + +CC = clang +CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -fopenmp +LDFLAGS = -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lmetis -lgsl -lprofiler -ltcmalloc + +OBJ = break_data voronoi_bound_ini bin_values correlations beta_scales randfuncs instance get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry fracture_network get_current cracking_ini update_factor update_boundary get_file update_beta gen_voltcurmat ini_network free_network fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory get_conductivity +BIN = fracture current_scaling course_grain_square corr_test homo_voronoi_fracture homo_square_fracture compare_voronoi_fracture + +all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} + +opt: + @echo build options: + @echo "CC = ${CC}" + @echo "CFLAGS = ${CFLAGS}" + @echo "LDFLAGS = ${LDFLAGS}" + +obj/%.o: src/%.c + @echo CC -c -o $@ + @${CC} -c -o $@ $< ${CFLAGS} + +bin/%: obj/%.o ${OBJ:%=obj/%.o} + @echo CC -o $@ + @${CC} -o $@ $< ${OBJ:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS} + +clean: + @echo cleaning: + @echo rm -f ${OBJ:%=obj/%} ${BIN:%=obj/%.o} ${BIN:%=bin/%} + @rm -f ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} + +.PHONY: all clean + diff --git a/src/beta_scales.c b/src/beta_scales.c new file mode 100644 index 0000000..8cbc5d1 --- /dev/null +++ b/src/beta_scales.c @@ -0,0 +1,21 @@ + +#include "fracture.h" + +double beta_scaling_flat(double beta, double x, double y) { return beta; } + +double beta_mag(double beta) { + double aa = -7.52579; + double bb = 9.63706; + double cc = 0.692515; + double dd = -2.47638; + + return gsl_sf_exp(aa + bb * gsl_sf_erf(cc * (gsl_sf_log(beta) - dd))); +} + +double beta_scaling_gaus(double beta, double x, double y) { + double sigma = 0.25; + double nu_f = 1.56; + double betap = beta * gsl_sf_exp((gsl_pow_2(x - 0.5) + gsl_pow_2(y - 0.5)) / + (2 * sigma * 2 * nu_f)); + return betap; +} diff --git a/src/bin_values.c b/src/bin_values.c new file mode 100644 index 0000000..2cb9cf6 --- /dev/null +++ b/src/bin_values.c @@ -0,0 +1,23 @@ + +#include "fracture.h" + +double *bin_values(fnet *network, unsigned int width, double *values) { + double *binned = calloc(pow(width, 2), sizeof(double)); + unsigned int *num_binned = calloc(pow(width, 2), sizeof(unsigned int)); + for (unsigned int i = 0; i < network->num_edges; i++) { + if (values[i] != 0) { + unsigned int x = ((int)(network->edge_coords[2 * i] * width)) % width; + unsigned int y = ((int)(network->edge_coords[2 * i + 1] * width)) % width; + binned[width * x + y] += fabs(values[i]); + num_binned[width * x + y]++; + } + } + + for (unsigned int i = 0; i < pow(width, 2); i++) { + if (num_binned[i] != 0) { + binned[i] /= num_binned[i]; + } + } + + return binned; +} diff --git a/src/break_data.c b/src/break_data.c new file mode 100644 index 0000000..0ce6730 --- /dev/null +++ b/src/break_data.c @@ -0,0 +1,34 @@ + +#include "fracture.h" + +break_data *alloc_break_data(unsigned int num_edges) { + break_data *data = malloc(1 * sizeof(break_data)); assert(data != NULL); + + data->num_broken = 0; + + data->break_list = (unsigned int *)malloc(num_edges * sizeof(unsigned int)); + assert(data->break_list != NULL); + + data->extern_field = (double *)malloc(num_edges * sizeof(double)); + assert(data->extern_field != NULL); + + data->conductivity = (double *)malloc(num_edges * sizeof(double)); + assert(data->conductivity != NULL); + + return data; +} + +void free_break_data(break_data *data) { + free(data->break_list); + free(data->extern_field); + free(data->conductivity); + free(data); +} + +void update_break_data(break_data *data, unsigned int last_broke, double strength, double conductivity) { + data->break_list[data->num_broken] = last_broke; + data->extern_field[data->num_broken] = strength; + data->conductivity[data->num_broken] = conductivity; + data->num_broken++; +} + diff --git a/src/break_edge.c b/src/break_edge.c new file mode 100644 index 0000000..0bd0277 --- /dev/null +++ b/src/break_edge.c @@ -0,0 +1,89 @@ + +#include "fracture.h" + +bool break_edge(finst *instance, unsigned int edge, cholmod_common *c) { + instance->fuses[edge] = true; + instance->num_remaining_edges--; + + unsigned int v1 = instance->network->edges_to_verts_break[2 * edge]; + unsigned int v2 = instance->network->edges_to_verts_break[2 * edge + 1]; + + update_factor(instance->factor, v1, v2, c); + + + unsigned int w1 = instance->network->edges_to_verts[2 * edge]; + unsigned int w2 = instance->network->edges_to_verts[2 * edge + 1]; + unsigned int dw1 = instance->network->dual_edges_to_verts[2 * edge]; + unsigned int dw2 = instance->network->dual_edges_to_verts[2 * edge + 1]; + + if (instance->network->boundary != TORUS_BOUND) { + + unsigned int tw1 = w1 > w2 ? w1 : w2; + unsigned int tw2 = w1 > w2 ? w2 : w1; + + CHOL_INT *lap_p = (CHOL_INT *)instance->adjacency->p; + CHOL_INT *lap_i = (CHOL_INT *)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[tw2 + 1] - lap_p[tw2]; i++) { + if (lap_i[lap_p[tw2] + i] == tw1) + lap_x[lap_p[tw2] + i] = 0; + } + + int old_num_components = instance->num_components; + + instance->num_components = update_components( + instance->adjacency, instance->marks, old_num_components, (int)tw1, + (int)tw2, 0); + } + + if (instance->network->boundary == TORUS_BOUND) { + if (instance->dual_marks[dw1] == instance->dual_marks[dw2]) { + int **cycles = (int **)malloc(2*instance->network->num_edges * sizeof(int *)); + unsigned int num_cycles = find_cycles(instance->network->num_edges, instance->fuses, instance->network->dual_edges_to_verts, instance->network->dual_verts_to_edges_ind, instance->network->dual_verts_to_edges, cycles); + + for (unsigned int i = 0; i < num_cycles; i++) { + int x_num_crossings = 0; + int y_num_crossings = 0; + int ee; unsigned int j = 0; + while ((ee = cycles[2*i][j]) >= 0) { + int side = cycles[2*i+1][j]; + j++; + unsigned int v1, v2; + double v1x, v1y, v2x, v2y; + v1 = instance->network->dual_edges_to_verts[2 * ee + !side]; + v2 = instance->network->dual_edges_to_verts[2 * ee + side]; + v1x = instance->network->dual_vert_coords[2 * v1]; + v1y = instance->network->dual_vert_coords[2 * v1 + 1]; + v2x = instance->network->dual_vert_coords[2 * v2]; + v2y = instance->network->dual_vert_coords[2 * v2 + 1]; + double dx = v1x - v2x; + double dy = v1y - v2y; + if (((v1x > 0.5 && v2x < 0.5) || (v1x < 0.5 && v2x > 0.5)) && fabs(dx) < 0.5) { + x_num_crossings += dx > 0 ? 1 : -1; + } + if (((v1y > 0.5 && v2y < 0.5) || (v1y < 0.5 && v2y > 0.5)) && fabs(dy) < 0.5) { + y_num_crossings += dy > 0 ? 1 : -1; + } + } + if ((abs(y_num_crossings) == 0 && abs(x_num_crossings) > 0) || (abs(y_num_crossings) > 0 && abs(x_num_crossings) > 0 && num_cycles > 1)) { + instance->num_components = 2; + } + free(cycles[2*i]); + free(cycles[2*i+1]); + } + free(cycles); + } + + unsigned int *tmp_marks = get_clusters(instance, c); + free(instance->dual_marks); + instance->dual_marks = tmp_marks; + } + + + return true; +} diff --git a/src/compare_voronoi_fracture.c b/src/compare_voronoi_fracture.c new file mode 100644 index 0000000..5578987 --- /dev/null +++ b/src/compare_voronoi_fracture.c @@ -0,0 +1,124 @@ + + +#include "fracture.h" + +int main(int argc, char *argv[]) { + int opt; + + // defining variables to be (potentially) set by command line flags + unsigned int N, L, filename_len; + double beta, inf, cutoff; + boundary_type boundary; + + filename_len = 100; + + N = 100; + L = 16; + beta = .3; + inf = 1e10; + cutoff = 1e-10; + boundary = FREE_BOUND; + + int boundary_int; + char boundc2 = 'f'; + + while ((opt = getopt(argc, argv, "n:L:b:B:dcoNsCrt")) != -1) { + switch (opt) { + case 'n': + N = atoi(optarg); + break; + case 'L': + L = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + case 'B': + boundary_int = atoi(optarg); + switch (boundary_int) { + case 0: + boundary = FREE_BOUND; + boundc2 = 'f'; + break; + case 1: + boundary = CYLINDER_BOUND; + boundc2 = 'c'; + break; + default: + printf("boundary specifier must be 0 (FREE_BOUND) or 1 (CYLINDER_BOUND).\n"); + } + break; + default: /* '?' */ + exit(EXIT_FAILURE); + } + } + + char *break_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(break_filename, filename_len, "breaks_v_vc_%c_%u_%g.txt", boundc2, L, beta); + FILE *break_out = fopen(break_filename, "a"); + free(break_filename); + + + // start cholmod + cholmod_common c; + CHOL_F(start)(&c); + + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + + + fnet *network = ini_voronoi_network(L, boundary, genfunc_hyperuniform, &c); + finst *perm_voltage_instance = create_instance(network, inf, true, true, &c); + finst *perm_current_instance = create_instance(network, inf, false, true, &c); + double *fuse_thres = gen_fuse_thres(network->num_edges, network->edge_coords, beta, beta_scaling_flat); + finst *voltage_instance = copy_instance(perm_voltage_instance, &c); + finst *current_instance = copy_instance(perm_current_instance, &c); + break_data *breaking_data_voltage = fracture_network(voltage_instance, fuse_thres, &c, cutoff); + break_data *breaking_data_current = fracture_network(current_instance, fuse_thres, &c, cutoff); + free_instance(voltage_instance, &c); + free_instance(current_instance, &c); + free_instance(perm_voltage_instance, &c); + free_instance(perm_current_instance, &c); + free(fuse_thres); + + FILE *net_out = fopen("network.txt", "w"); + for (unsigned int j = 0; j < network->num_verts; j++) { + fprintf(net_out, "%f %f ", network->vert_coords[2 * j], + network->vert_coords[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < network->num_edges; j++) { + fprintf(net_out, "%u %u ", network->edges_to_verts[2 * j], + network->edges_to_verts[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < network->num_dual_verts; j++) { + fprintf(net_out, "%f %f ", network->dual_vert_coords[2 * j], + network->dual_vert_coords[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < network->num_edges; j++) { + fprintf(net_out, "%u %u ", network->dual_edges_to_verts[2 * j], + network->dual_edges_to_verts[2 * j + 1]); + } + + free_fnet(network, &c); + + for (unsigned int j = 0; j < breaking_data_voltage->num_broken; j++) { + fprintf(break_out, "%u %g %g ", breaking_data_voltage->break_list[j], + breaking_data_voltage->extern_field[j], breaking_data_voltage->conductivity[j]); + } + fprintf(break_out, "\n"); + for (unsigned int j = 0; j < breaking_data_current->num_broken; j++) { + fprintf(break_out, "%u %g %g ", breaking_data_current->break_list[j], + breaking_data_current->extern_field[j], breaking_data_current->conductivity[j]); + } + fprintf(break_out, "\n"); + + free_break_data(breaking_data_voltage); + free_break_data(breaking_data_current); + fclose(break_out); + + CHOL_F(finish)(&c); + + return 0; +} diff --git a/src/corr_test.c b/src/corr_test.c new file mode 100644 index 0000000..2c93b8e --- /dev/null +++ b/src/corr_test.c @@ -0,0 +1,30 @@ + +#include "fracture.h" + +int main() { + cholmod_common c; + CHOL_F(start)(&c); + + unsigned int width = 64; + unsigned int n = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2); + + fnet *network = ini_square_network(width, true, false, &c); + finst *instance = create_instance(network, 1e-14, true, true, &c); + double *fuse_thres = gen_fuse_thres(network->num_edges, network->edge_coords, + 0.001, beta_scaling_flat); + fracture_network(instance, fuse_thres, &c, 1e-10); + + double *corr = get_corr(instance, NULL, &c); + + for (int i = 0; i < 2 * width; i++) { + printf("%f ", corr[i]); + } + printf("\n"); + + free_instance(instance, &c); + free_fnet(network, &c); + + CHOL_F(finish)(&c); + + return 0; +} diff --git a/src/correlations.c b/src/correlations.c new file mode 100644 index 0000000..bd05daf --- /dev/null +++ b/src/correlations.c @@ -0,0 +1,345 @@ + +#include "fracture.h" + +struct fibn { + unsigned int value; + unsigned int priority; + unsigned int rank; + bool marked; + struct fibn *first_child; + struct fibn *parent; + struct fibn *prev; + struct fibn *next; +}; + +typedef struct { + struct fibn *min; + unsigned int rank; + unsigned int trees; + struct fibn **vertex_to_node; +} fibh; + +unsigned int fib_findmin(fibh *heap) { return heap->min->value; } + +void ll_setparent(struct fibn *ll, struct fibn *parent) { + struct fibn *curnode = ll->next; + ll->parent = parent; + while (curnode != ll) { + curnode->parent = parent; + curnode = curnode->next; + } +} + +struct fibn *ll_merge(struct fibn *ll1, struct fibn *ll2) { + if (ll1 == NULL) + return ll2; + if (ll2 == NULL) + return ll1; + + // link the beginning of list one to the end of list two and vice versa + struct fibn *ll1_beg = ll1; + struct fibn *ll1_end = ll1->prev; + struct fibn *ll2_beg = ll2; + struct fibn *ll2_end = ll2->prev; + + ll1_beg->prev = ll2_end; + ll1_end->next = ll2_beg; + ll2_beg->prev = ll1_end; + ll2_end->next = ll1_beg; + + return ll1; +} + +struct fibn *ll_delroot(struct fibn *ll) { + if (ll == NULL) + return NULL; + + if (ll->next == ll) { + return NULL; + } + + struct fibn *ll_beg = ll->next; + struct fibn *ll_end = ll->prev; + + ll_beg->prev = ll_end; + ll_end->next = ll_beg; + + ll->next = ll; + ll->prev = ll; + + return ll_beg; +} + +void fib_insert(fibh *heap, unsigned int value, unsigned int priority) { + struct fibn *newnode = calloc(1, sizeof(struct fibn)); + newnode->value = value; + newnode->priority = priority; + newnode->next = newnode; + newnode->prev = newnode; + + heap->min = ll_merge(heap->min, newnode); + + if (priority < heap->min->priority) { + heap->min = newnode; + } + + heap->vertex_to_node[value] = newnode; + + heap->trees++; +} + +void fib_deletemin(fibh *heap) { + unsigned int min_rank = heap->min->rank; + struct fibn *min_children = heap->min->first_child; + + struct fibn *trees = ll_delroot(heap->min); + heap->vertex_to_node[heap->min->value] = NULL; + free(heap->min); + + if (min_children != NULL) + ll_setparent(min_children, NULL); + trees = ll_merge(trees, min_children); + + heap->trees += min_rank - 1; + + if (trees == NULL) { + // min had no children and was only tree, return empty heap + heap->min = NULL; + heap->rank = 0; + } else { + // min had children or there were other trees + unsigned int min_val = UINT_MAX; + struct fibn *min_ptr = NULL; + unsigned int max_rank = 0; + struct fibn *curnode = trees; + for (unsigned int i = 0; i < heap->trees; i++) { + if (curnode->priority < min_val) { + min_ptr = curnode; + min_val = curnode->priority; + } + if (curnode->rank > max_rank) + max_rank = curnode->rank; + curnode = curnode->next; + } + + if (min_ptr == NULL) + min_ptr = trees; + heap->min = min_ptr; + heap->rank = max_rank; + + struct fibn **rankslots = calloc(max_rank + 1, sizeof(struct fibn *)); + curnode = heap->min; + while (curnode != rankslots[curnode->rank]) { + if (rankslots[curnode->rank] == NULL) { + rankslots[curnode->rank] = curnode; + curnode = curnode->next; + } else { + struct fibn *oldnode = rankslots[curnode->rank]; + rankslots[curnode->rank] = NULL; + struct fibn *smaller = + curnode->priority < oldnode->priority || curnode == heap->min + ? curnode + : oldnode; + struct fibn *larger = smaller == curnode ? oldnode : curnode; + ll_delroot(larger); + ll_setparent(larger, smaller); + struct fibn *smaller_children = smaller->first_child; + smaller->first_child = ll_merge(smaller_children, larger); + heap->trees--; + smaller->rank++; + if (smaller->rank > heap->rank) { + heap->rank = smaller->rank; + rankslots = + realloc(rankslots, (heap->rank + 1) * sizeof(struct fibn *)); + rankslots[heap->rank] = NULL; + } + curnode = smaller; + } + } + free(rankslots); + } +} + +void fib_decreasekey(fibh *heap, unsigned int value, + unsigned int new_priority) { + struct fibn *node = heap->vertex_to_node[value]; + if (node != NULL) { + node->priority = new_priority; + if (node->parent != NULL) { + if (node->priority < node->parent->priority) { + + struct fibn *curnode = node; + curnode->marked = true; + while (curnode->parent != NULL && curnode->marked) { + struct fibn *oldparent = curnode->parent; + oldparent->rank--; + oldparent->first_child = ll_delroot(curnode); + ll_setparent(curnode, NULL); + ll_merge(heap->min, curnode); + heap->trees++; + if (curnode->marked) { + curnode->marked = false; + } + + if (oldparent->marked) { + curnode = oldparent; + } else { + oldparent->marked = true; + break; + } + } + } + } + + if (new_priority < heap->min->priority) { + heap->min = node; + } + } +} + +unsigned int *dijkstra(fnet *network, unsigned int source) { + unsigned int nv = network->num_dual_verts; + unsigned int *vei = network->dual_verts_to_edges_ind; + unsigned int *ve = network->dual_verts_to_edges; + unsigned int *ev = network->dual_edges_to_verts; + + unsigned int *dist = (unsigned int *)calloc(nv, sizeof(unsigned int)); + fibh *Q = (fibh *)calloc(1, sizeof(fibh)); + Q->vertex_to_node = (struct fibn **)calloc(nv, sizeof(struct fibn *)); + + for (unsigned int i = 0; i < nv; i++) { + if (i != source) { + dist[i] = UINT_MAX; + } + + fib_insert(Q, i, dist[i]); + } + + while (Q->min != NULL) { + unsigned int u = fib_findmin(Q); + fib_deletemin(Q); + + for (unsigned int i = 0; i < vei[u + 1] - vei[u]; i++) { + unsigned int e = ve[vei[u] + i]; + unsigned int v = ev[2 * e] == u ? ev[2 * e + 1] : ev[2 * e]; + unsigned int alt = dist[u] + 1; + if (alt < dist[v]) { + dist[v] = alt; + fib_decreasekey(Q, v, alt); + } + } + } + + free(Q->vertex_to_node); + free(Q); + + return dist; +} + +unsigned int **get_dists(fnet *network) { + unsigned int nv = network->num_dual_verts; + unsigned int **dists = (unsigned int **)malloc(nv * sizeof(unsigned int *)); + +#pragma omp parallel for + for (unsigned int i = 0; i < nv; i++) { + dists[i] = dijkstra(network, i); + } + + return dists; +} + +double *get_corr(finst *instance, unsigned int **dists, cholmod_common *c) { + unsigned int nv = instance->network->num_dual_verts; + unsigned int ne = instance->network->num_edges; + unsigned int *ev = instance->network->dual_edges_to_verts; + bool nulldists = false; + if (dists == NULL) { + dists = get_dists(instance->network); + nulldists = true; + } + double *corr = calloc(nv, sizeof(double)); + unsigned int *marks = get_clusters(instance, c); + unsigned int *numat = calloc(nv, sizeof(unsigned int)); + + for (unsigned int i = 0; i < ne; i++) { + unsigned int v1 = ev[2 * i]; + unsigned int v2 = ev[2 * i + 1]; + for (unsigned int j = 0; j < ne; j++) { + unsigned int v3 = ev[2 * j]; + unsigned int v4 = ev[2 * j + 1]; + unsigned int dist1 = dists[v1][v3]; + unsigned int dist2 = dists[v1][v4]; + unsigned int dist3 = dists[v2][v3]; + unsigned int dist4 = dists[v2][v4]; + unsigned int dist = (dist1 + dist2 + dist3 + dist4) / 4; + corr[dist] += instance->fuses[i] && instance->fuses[j]; + numat[dist]++; + } + } + + for (unsigned int i = 0; i < nv; i++) { + if (numat[i] > 0) { + corr[i] /= numat[i]; + } + } + + if (nulldists) { + for (int i = 0; i < nv; i++) { + free(dists[i]); + } + free(dists); + } + + free(marks); + + return corr; +} + +/* +double *get_space_corr(finst *instance, cholmod_common *c) { + unsigned int nv = instance->network->num_dual_verts; + double *vc = instance->network->dual_vert_coords; + unsigned int ne = instance->network->num_edges; + unsigned int *ev = instance->network->dual_edges_to_verts; + double *corr = calloc(nv, sizeof(unsigned int)); + unsigned int *numat = calloc(nv, sizeof(unsigned int)); + + for (unsigned int i = 0; i < ne; i++) { + unsigned int v1 = ev[2 * i]; + unsigned int v2 = ev[2 * i + 1]; + double v1x = vc[2 * v1]; + double v1y = vc[2 * v1 + 1]; + double v2x = vc[2 * v2]; + double v2y = vc[2 * v2 + 1]; + double e1x = (v1x + v2x) / 2; + double e1y = (v1y + v2y) / 2; + for (unsigned int j = 0; j < ne; j++) { + unsigned int v3 = ev[2 * j]; + unsigned int v4 = ev[2 * j + 1]; + double v1x = vc[2 * v1]; + double v1y = vc[2 * v1 + 1]; + double v2x = vc[2 * v2]; + double v2y = vc[2 * v2 + 1]; + double e1x = (v1x + v2x) / 2; + double e1y = (v1y + v2y) / 2; + corr[dist] += instance->fuses[i] && instance->fuses[j]; + numat[dist]++; + } + } + + for (unsigned int i = 0; i < nv; i++) { + if (numat[i] > 0) { + corr[i] /= numat[i]; + } + } + + for (int i = 0; i < nv; i++) { + free(dists[i]); + } + + free(marks); + free(dists); + + return corr; +} +*/ diff --git a/src/course_grain_square.c b/src/course_grain_square.c new file mode 100644 index 0000000..9035cf7 --- /dev/null +++ b/src/course_grain_square.c @@ -0,0 +1,168 @@ + +#include "fracture.h" + +int main(int argc, char *argv[]) { + int opt; + + // defining variables to be (potentially) set by command line flags + int num = 100; + int width = 16; + double beta = .3; + bool save_clusters = false; + bool voltage_bound = false; + bool output_break_data = false; + + while ((opt = getopt(argc, argv, "n:w:b:cVo")) != -1) { + switch (opt) { + case 'n': + num = atoi(optarg); + break; + case 'w': + width = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + case 'c': + save_clusters = true; + break; + case 'V': + voltage_bound = true; + break; + case 'o': + output_break_data = true; + break; + default: /* '?' */ + exit(EXIT_FAILURE); + } + } + + FILE *break_out; + if (output_break_data) { + char *break_filename = (char *)malloc(100 * sizeof(char)); + snprintf(break_filename, 100, "breaks_%d_%.3f_%d.txt", width, beta, num); + break_out = fopen(break_filename, "w"); + free(break_filename); + } + + bool periodic = true; + double inf = 1; + double cutoff = 1e-10; + + // start cholmod + cholmod_common c; + CHOL_F(start)(&c); + + /* if we use voltage boundary conditions, the laplacian matrix is positive + * definite and we can use a supernodal LL decomposition. otherwise we need + * to use the simplicial LDL decomposition + */ + if (voltage_bound) { + (&c)->supernodal = CHOLMOD_SUPERNODAL; + } else { + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + } + + fnet *network = ini_square_network(width, periodic, false, &c); + fnet *network_p = ini_square_network(width / 2, periodic, false, &c); + finst *perm_instance = create_instance(network, inf, voltage_bound, true, &c); + unsigned int c_dist_size = network->num_dual_verts; + unsigned int c_p_dist_size = network_p->num_dual_verts; + + // define arrays for saving cluster and avalanche distributions + unsigned int *cluster_size_dist = + (unsigned int *)calloc(c_dist_size, sizeof(unsigned int)); + unsigned int *cluster_p_size_dist = + (unsigned int *)calloc(c_p_dist_size, sizeof(unsigned int)); + + printf("\n"); + for (int DUMB = 0; DUMB < num; DUMB++) { + printf("\033[F\033[JCOURSEGRAIN_SQUARE: %0*d / %d\n", (int)log10(num) + 1, + DUMB + 1, num); + + break_data *breaking_data = NULL; + finst *instance = NULL; + while (breaking_data == NULL) { + double *fuse_thres = gen_fuse_thres( + network->num_edges, network->edge_coords, beta, beta_scaling_flat); + instance = copy_instance(perm_instance, &c); + breaking_data = fracture_network(instance, fuse_thres, &c, cutoff); + free_instance(instance, &c); + free(fuse_thres); + } + + unsigned int min_pos = 0; + double min_val = DBL_MAX; + + for (unsigned int j = 0; j < breaking_data->num_broken; j++) { + double val = fabs(breaking_data->extern_field[j]); + if (val < min_val) { + min_pos = j; + min_val = val; + } + if (val > 10 * min_val) + break; + } + + finst *tmp_instance = copy_instance(perm_instance, &c); + + for (unsigned int i = 0; i < breaking_data->num_broken; i++) { + break_edge(tmp_instance, breaking_data->break_list[i], &c); + } + + unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c); + for (unsigned int i = 0; i < network->num_dual_verts; i++) { + cluster_size_dist[i] += tmp_cluster_dist[i]; + } + free(tmp_cluster_dist); + + finst *instance_p = coursegrain_square(tmp_instance, network_p, &c); + + unsigned int *tmp_cluster_p_dist = get_cluster_dist(instance_p, &c); + for (unsigned int i = 0; i < network_p->num_dual_verts; i++) { + cluster_p_size_dist[i] += tmp_cluster_p_dist[i]; + } + free(tmp_cluster_p_dist); + + free_instance(tmp_instance, &c); + free_instance(instance_p, &c); + if (output_break_data) { + for (unsigned int i = 0; i < breaking_data->num_broken; i++) { + fprintf(break_out, "%u %f ", breaking_data->break_list[i], + breaking_data->extern_field[i]); + } + fprintf(break_out, "\n"); + } + free(breaking_data->break_list); + free(breaking_data->extern_field); + free(breaking_data); + } + + printf("\033[F\033[JCURRENT_SCALING: COMPLETE"); + + if (save_clusters) { + FILE *cluster_out = get_file("cluster", width, 0, beta, 1, 1, num, false); + for (int i = 0; i < c_dist_size; i++) { + fprintf(cluster_out, "%u ", cluster_size_dist[i]); + } + fprintf(cluster_out, "\n"); + for (int i = 0; i < c_p_dist_size; i++) { + fprintf(cluster_out, "%u ", cluster_p_size_dist[i]); + } + fclose(cluster_out); + } + + if (output_break_data) { + fclose(break_out); + } + + free(cluster_size_dist); + free(cluster_p_size_dist); + free_instance(perm_instance, &c); + free_fnet(network, &c); + free_fnet(network_p, &c); + + CHOL_F(finish)(&c); + + return 0; +} diff --git a/src/coursegrain.c b/src/coursegrain.c new file mode 100644 index 0000000..bd84e1c --- /dev/null +++ b/src/coursegrain.c @@ -0,0 +1,40 @@ + +#include "fracture.h" + +finst *coursegrain_square(finst *instance, fnet *network_p, cholmod_common *c) { + unsigned int width = sqrt(instance->network->num_edges); + assert(width % 4 == 0); + + finst *instance_p = create_instance(network_p, instance->inf, + instance->voltage_bound, true, c); + + unsigned int width_p = width / 2; + bool *fuses = instance->fuses; + + for (unsigned int i = 0; i < network_p->num_edges; i++) { + int xp = i / width_p; + int yp = i % width_p; + unsigned int x1 = 2 * xp; + unsigned int y1 = (2 * yp - 1) % width; + unsigned int x2 = 2 * xp; + unsigned int y2 = 2 * yp; + unsigned int x3 = 2 * xp + 1; + unsigned int y3 = (2 * yp - 1) % width; + unsigned int x4 = 2 * xp + 1; + unsigned int y4 = 2 * yp; + bool f1 = fuses[width * x1 + y1]; + bool f2 = fuses[width * x2 + y2]; + bool f3 = fuses[width * x3 + y3]; + bool f4 = fuses[width * x4 + y4]; + + if ((f1 && f2) || (f3 && f4)) { + // instance_p->fuses[i] = true; + // instance_p->num_remaining_edges--; + break_edge(instance_p, i, c); + } + } + + // fin_instance(instance_p, c); + + return instance_p; +} diff --git a/src/cracking_ini.c b/src/cracking_ini.c new file mode 100644 index 0000000..9ab4a9d --- /dev/null +++ b/src/cracking_ini.c @@ -0,0 +1,53 @@ + +#include "fracture.h" + +double *gen_fuse_thres(unsigned int num_edges, double *edge_coords, double beta, + double (*beta_scale)(double, double, double)) { + unsigned int size = num_edges; + assert(beta > 0); + + double *fuse_thres = (double *)malloc(size * sizeof(double)); + assert(fuse_thres != NULL); + + gsl_set_error_handler_off(); + + gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); + FILE *rf = fopen("/dev/urandom", "r"); + unsigned long int seed; + fread(&seed, sizeof(unsigned long int), 1, rf); + fclose(rf); + gsl_rng_set(r, seed); + + for (unsigned int i = 0; i < size; i++) { + double x = edge_coords[2 * i]; + double y = edge_coords[2 * i + 1]; + while ((fuse_thres[i] = gsl_sf_exp(gsl_sf_log(gsl_ran_flat(r, 0, 1)) / + beta_scale(beta, x, y))) == 0.0) + ; + } + + gsl_rng_free(r); + + return fuse_thres; +} + +bool gen_crack(finst *instance, double crack_len, double crack_width, + cholmod_common *c) { + assert(instance != NULL); + bool *fuses = instance->fuses; + assert(fuses != NULL); + const fnet *network = instance->network; + unsigned int num_edges = network->num_edges; + double *edge_coords = network->edge_coords; + + for (unsigned int j = 0; j < num_edges; j++) { + if (edge_coords[2 * j + 1] < crack_len && + fabs(edge_coords[2 * j] - network->L / 2) < crack_width) { + instance->fuses[j] = true; + instance->num_remaining_edges--; + } else + fuses[j] = false; + } + + return true; +} diff --git a/src/current_scaling.c b/src/current_scaling.c new file mode 100644 index 0000000..6837f6b --- /dev/null +++ b/src/current_scaling.c @@ -0,0 +1,300 @@ + +#include "fracture.h" + +int main(int argc, char *argv[]) { + int opt; + + // defining variables to be (potentially) set by command line flags + int iter = 1; + int num = 100; + int width = 16; + double crack_len = 8; + double beta = .3; + double inf = 1e-8; + double cutoff = 1e-8; + bool beta_shift = false; + bool supplied_bound = false; + bool ash_beta = false; + char *bound_file; + bool voltage_bound = false; + bool use_first = false; + bool save_stress = false; + bool save_bound = false; + bool save_damage = false; + bool save_strength = false; + + while ((opt = getopt(argc, argv, "n:w:b:l:i:Bf:aVFsSed")) != -1) { + switch (opt) { + case 'n': + num = atoi(optarg); + break; + case 'w': + width = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + case 'l': + crack_len = atof(optarg); + break; + case 'i': + iter = atoi(optarg); + break; + case 'B': + beta_shift = true; + break; + case 'a': + ash_beta = true; + break; + case 'V': + voltage_bound = true; + break; + case 'F': + use_first = true; + break; + case 's': + save_stress = true; + break; + case 'd': + save_damage = true; + break; + case 'e': + save_bound = true; + break; + case 'S': + save_strength = true; + break; + case 'f': + supplied_bound = true; + bound_file = optarg; + break; + default: /* '?' */ + exit(EXIT_FAILURE); + } + } + + // start cholmod + cholmod_common c; + CHOL_F(start)(&c); + + /* if we use voltage boundary conditions, the laplacian matrix is positive + * definite and we can use a supernodal LL decomposition. otherwise we need + * to use the simplicial LDL decomposition + */ + if (voltage_bound) { + (&c)->supernodal = CHOLMOD_SUPERNODAL; + } else { + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + } + + fnet *network = ini_square_network(width, false, true, &c); + finst *perm_instance = + create_instance(network, inf, voltage_bound, false, &c); + gen_crack(perm_instance, crack_len, 1, &c); + finish_instance(perm_instance, &c); + + if (voltage_bound) { + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + finst *tmp_instance = create_instance(network, inf, false, false, &c); + gen_crack(tmp_instance, crack_len, 1, &c); + finish_instance(tmp_instance, &c); + double *voltage = get_voltage(tmp_instance, &c); + + for (int i = 0; i < network->num_bounds; i++) { + for (int j = 0; j < network->bound_inds[i + 1] - network->bound_inds[i]; + j++) { + ((double *)perm_instance->boundary_cond + ->x)[network->bound_verts[network->bound_inds[i] + j]] = + voltage[network->bound_verts[network->bound_inds[i] + j]]; + } + } + (&c)->supernodal = CHOLMOD_SUPERNODAL; + } + + if (supplied_bound) { + FILE *bound_f = fopen(bound_file, "r"); + for (int i = 0; i < network->num_verts; i++) { + double tmp; + fscanf(bound_f, "%lg ", &tmp); + ((double *)perm_instance->boundary_cond->x)[i] = tmp; + } + + ((double *)perm_instance->boundary_cond->x)[network->num_verts] = 0; + ((double *)perm_instance->boundary_cond->x)[network->num_verts + 1] = 0; + fclose(bound_f); + } + + printf("\n"); + for (int DUMB2 = 0; DUMB2 < iter; DUMB2++) { + + double *strength; + if (save_strength) { + strength = (double *)malloc(num * sizeof(double)); + } + + double *damage; + if (save_damage) { + damage = (double *)calloc(network->num_edges, sizeof(double)); + } + double *avg_current = (double *)calloc(network->num_edges, sizeof(double)); + unsigned int *num_current_skipped = + (unsigned int *)calloc(network->num_edges, sizeof(unsigned int)); + double *avg_voltage = (double *)calloc(network->num_verts, sizeof(double)); + unsigned int *num_voltage_skipped = + (unsigned int *)calloc(network->num_verts, sizeof(unsigned int)); + + for (int DUMB = 0; DUMB < num; DUMB++) { + printf("\033[F\033[JCURRENT_SCALING: ITERATION %0*d: %0*d / %d\n", + (int)log10(iter) + 1, DUMB2 + 1, (int)log10(num) + 1, DUMB + 1, + num); + + break_data *breaking_data = NULL; + while (breaking_data == NULL) { + double *fuse_thres = gen_fuse_thres( + network->num_edges, network->edge_coords, beta, beta_scaling_flat); + finst *instance = copy_instance(perm_instance, &c); + breaking_data = fracture_network(instance, fuse_thres, &c, cutoff); + free_instance(instance, &c); + free(fuse_thres); + } + + unsigned int min_pos = 0; + double min_val = DBL_MAX; + + for (unsigned int j = 0; j < breaking_data->num_broken; j++) { + double val = fabs(breaking_data->extern_field[j]); + if (val < min_val) { + min_pos = j; + min_val = val; + } + } + + if (save_strength) { + strength[DUMB] = fabs(breaking_data->extern_field[min_pos]); + } + + finst *tmp_instance = copy_instance(perm_instance, &c); + + unsigned int until = min_pos; + if (use_first) + until = 1; + for (unsigned int i = 0; i < until; i++) { + break_edge(tmp_instance, breaking_data->break_list[i], &c); + if (save_damage) { + damage[breaking_data->break_list[i]] += 1. / num; + } + } + + double *voltage = get_voltage(tmp_instance, &c); + double *current = get_current(tmp_instance, &c); + + for (unsigned int i = 0; i < network->num_edges; i++) { + avg_current[i] += current[i]; + if (current[i] == 0) + num_current_skipped[i]++; + } + + for (unsigned int i = 0; i < network->num_verts; i++) { + if (tmp_instance->marks[i] == tmp_instance->marks[network->num_verts]) { + avg_voltage[i] += voltage[i]; + } else { + num_voltage_skipped[i]++; + } + } + + free(current); + free(voltage); + free_instance(tmp_instance, &c); + free(breaking_data->break_list); + free(breaking_data->extern_field); + free(breaking_data); + } + + for (int i = 0; i < network->num_edges; i++) { + if (num_current_skipped[i] < num) { + avg_current[i] /= num - num_current_skipped[i]; + } + } + + for (int i = 0; i < network->num_verts; i++) { + if (num_voltage_skipped[i] < num) { + avg_voltage[i] /= num - num_voltage_skipped[i]; + } + } + + double *avg_field; + if (voltage_bound) + avg_field = avg_voltage; + else + avg_field = avg_current; + + update_boundary(perm_instance, avg_field); + + if (save_stress) { + char *c_filename = (char *)malloc(100 * sizeof(char)); + snprintf(c_filename, 100, "current_%d_%g_%d_%g.txt", width, crack_len, + iter, beta); + FILE *outfile = fopen(c_filename, "w"); + for (int i = 0; i < network->num_edges; i++) { + fprintf(outfile, "%g ", avg_current[i]); + } + fclose(outfile); + free(c_filename); + } + + if (save_damage) { + char *c_filename = (char *)malloc(100 * sizeof(char)); + snprintf(c_filename, 100, "damage_%d_%g_%d_%g.txt", width, crack_len, + iter, beta); + FILE *outfile = fopen(c_filename, "w"); + for (int i = 0; i < network->num_edges; i++) { + fprintf(outfile, "%g ", damage[i]); + } + fclose(outfile); + free(c_filename); + } + + if (save_strength) { + char *s_filename = (char *)malloc(100 * sizeof(char)); + snprintf(s_filename, 100, "strength_%d_%g_%d_%g.txt", width, crack_len, iter, beta); + FILE *f = fopen(s_filename, "a"); + for (int i = 0; i < num; i++) { + fprintf(f, "%g ", strength[i]); + } + fclose(f); + free(s_filename); + } + + if (save_bound) { + char *b_filename = (char *)malloc(100 * sizeof(char)); + snprintf(b_filename, 100, "bounds_%d_%g_%d_%g.txt", width, crack_len, + iter, beta); + FILE *outfile = fopen(b_filename, "w"); + for (int i = 0; i < network->num_verts; i++) { + fprintf(outfile, "%g ", ((double *)perm_instance->boundary_cond->x)[i]); + } + fclose(outfile); + free(b_filename); + } + + free(avg_current); + free(avg_voltage); + if (save_damage) free(damage); + free(num_current_skipped); + free(num_voltage_skipped); + if (save_strength) { + free(strength); + } + + printf( + "\033[F\033[JCURRENT_SCALING: ITERATION %0*d COMPLETE, BETA = %.2f\n\n", + (int)log10(iter) + 1, DUMB2 + 1, beta); + } + + free_instance(perm_instance, &c); + free_fnet(network, &c); + + CHOL_F(finish)(&c); + + return 0; +} diff --git a/src/fortune/defs.h b/src/fortune/defs.h new file mode 100644 index 0000000..b186581 --- /dev/null +++ b/src/fortune/defs.h @@ -0,0 +1,134 @@ +#ifndef NULL +#define NULL 0 +#endif +#define DELETED -2 +#include <stdlib.h> +#include <limits.h> +#include <stdbool.h> + +extern int triangulate, sorted, plot, debug; + +struct Freenode { + struct Freenode *nextfree; +}; +struct Freelist { + struct Freenode *head; + int nodesize; +}; +char *getfree(); +char *myalloc(); + +extern double xmin, xmax, ymin, ymax, deltax, deltay; + +struct Point { + double x, y; +}; + +/* structure used both for sites and for vertices */ +struct Site { + struct Point coord; + int sitenbr; + int refcnt; +}; + +extern struct Site *sites; +extern int nsites; +extern int siteidx; +extern int sqrt_nsites; +extern int nvertices; +extern struct Freelist sfl; +extern struct Site *bottomsite; + +char **alloclist; +int allocnum; + +struct Edge { + double a, b, c; + struct Site *ep[2]; + struct Site *reg[2]; + int edgenbr; +}; +#define le 0 +#define re 1 +extern int nedges; +extern struct Freelist efl; + +int has_endpoint(), right_of(); +struct Site *intersect(); +double dist(); +struct Point PQ_min(); +struct Halfedge *PQextractmin(); +struct Edge *bisect(); + +struct Halfedge { + struct Halfedge *ELleft, *ELright; + struct Edge *ELedge; + int ELrefcnt; + char ELpm; + struct Site *vertex; + double ystar; + struct Halfedge *PQnext; +}; + +extern struct Freelist hfl; +extern struct Halfedge *ELleftend, *ELrightend; +extern int ELhashsize; +extern struct Halfedge **ELhash; +struct Halfedge *HEcreate(), *ELleft(), *ELright(), *ELleftbnd(); +struct Site *leftreg(), *rightreg(); + +extern int PQhashsize; +extern struct Halfedge *PQhash; +extern struct Halfedge *PQfind(); +extern int PQcount; +extern int PQmin; +int PQempty(); + +extern double *site_list; +extern double *vert_list; +extern double *line_list; +extern unsigned int *edge_list; +extern unsigned int *dual_list; + +extern int vert_count; +extern int edge_count; +extern int line_count; +extern int dual_count; +extern int site_count; + +void makefree(struct Freenode *curr, struct Freelist *fl); +void voronoi(int triangulate, struct Site *(*nextsite)()); +char *getfree(struct Freelist *fl); +char *myalloc(unsigned n); +void freeinit(struct Freelist *fl, int size); +void ELinitialize(); +struct Halfedge *HEcreate(struct Edge *e, int pm); +void ELinsert(struct Halfedge *lb, struct Halfedge *new); +struct Halfedge *ELgethash(int b); +struct Halfedge *ELleftbnd(struct Point *p); +void ELdelete(struct Halfedge *he); +struct Halfedge *ELright(struct Halfedge *he); +struct Halfedge *ELleft(struct Halfedge *he); +struct Site *leftreg(struct Halfedge *he); +struct Site *rightreg(struct Halfedge *he); +void geominit(); +struct Edge *bisect(struct Site *s1, struct Site *s2); +struct Site *intersect(struct Halfedge *el1, struct Halfedge *el2); +int right_of(struct Halfedge *el, struct Point *p); +void endpoint(struct Edge *e, int lr, struct Site *s); +double dist(struct Site *s, struct Site *t); +void makevertex(struct Site *v); +void deref(struct Site *v); +void ref(struct Site *v); +void PQinsert(struct Halfedge *he, struct Site *v, double offset); +void PQdelete(struct Halfedge *he); +int PQbucket(struct Halfedge *he); +int PQempty(); +struct Point PQ_min(); +struct Halfedge *PQextractmin(); +void PQinitialize(); +void out_bisector(struct Edge *e); +void out_ep(struct Edge *e); +void out_vertex(struct Site *v); +void out_site(struct Site *s); +void out_triple(struct Site *s1, struct Site *s2, struct Site *s3); diff --git a/src/fortune/edgelist.c b/src/fortune/edgelist.c new file mode 100644 index 0000000..dfd85a7 --- /dev/null +++ b/src/fortune/edgelist.c @@ -0,0 +1,136 @@ +# +#include "defs.h" +int nedges; +struct Freelist efl; +struct Freelist hfl; +struct Halfedge *ELleftend, *ELrightend; +int ELhashsize; +struct Halfedge **ELhash; + +int ntry, totalsearch; + +void ELinitialize() { + int i; + freeinit(&hfl, sizeof **ELhash); + ELhashsize = 2 * sqrt_nsites; + ELhash = (struct Halfedge **)myalloc(sizeof *ELhash * ELhashsize); + for (i = 0; i < ELhashsize; i += 1) + ELhash[i] = (struct Halfedge *)NULL; + ELleftend = HEcreate((struct Edge *)NULL, 0); + ELrightend = HEcreate((struct Edge *)NULL, 0); + ELleftend->ELleft = (struct Halfedge *)NULL; + ELleftend->ELright = ELrightend; + ELrightend->ELleft = ELleftend; + ELrightend->ELright = (struct Halfedge *)NULL; + ELhash[0] = ELleftend; + ELhash[ELhashsize - 1] = ELrightend; +} + +struct Halfedge *HEcreate(e, pm) struct Edge *e; +int pm; +{ + struct Halfedge *answer; + answer = (struct Halfedge *)getfree(&hfl); + answer->ELedge = e; + answer->ELpm = pm; + answer->PQnext = (struct Halfedge *)NULL; + answer->vertex = (struct Site *)NULL; + answer->ELrefcnt = 0; + return (answer); +} + +void ELinsert(struct Halfedge *lb, struct Halfedge *new) { + new->ELleft = lb; + new->ELright = lb->ELright; + (lb->ELright)->ELleft = new; + lb->ELright = new; +} + +/* Get entry from hash table, pruning any deleted nodes */ +struct Halfedge *ELgethash(b) int b; +{ + struct Halfedge *he; + + if (b < 0 || b >= ELhashsize) + return ((struct Halfedge *)NULL); + he = ELhash[b]; + if (he == (struct Halfedge *)NULL || he->ELedge != (struct Edge *)DELETED) + return (he); + + /* Hash table points to deleted half edge. Patch as necessary. */ + ELhash[b] = (struct Halfedge *)NULL; + if ((he->ELrefcnt -= 1) == 0) + makefree((struct Freenode *)he, &hfl); + return ((struct Halfedge *)NULL); +} + +struct Halfedge *ELleftbnd(p) struct Point *p; +{ + int i, bucket; + struct Halfedge *he; + + /* Use hash table to get close to desired halfedge */ + bucket = (p->x - xmin) / deltax * ELhashsize; + if (bucket < 0) + bucket = 0; + if (bucket >= ELhashsize) + bucket = ELhashsize - 1; + he = ELgethash(bucket); + if (he == (struct Halfedge *)NULL) { + for (i = 1; 1; i += 1) { + if ((he = ELgethash(bucket - i)) != (struct Halfedge *)NULL) + break; + if ((he = ELgethash(bucket + i)) != (struct Halfedge *)NULL) + break; + }; + totalsearch += i; + }; + ntry += 1; + /* Now search linear list of halfedges for the corect one */ + if (he == ELleftend || (he != ELrightend && right_of(he, p))) { + do { + he = he->ELright; + } while (he != ELrightend && right_of(he, p)); + he = he->ELleft; + } else + do { + he = he->ELleft; + } while (he != ELleftend && !right_of(he, p)); + + /* Update hash table and reference counts */ + if (bucket > 0 && bucket < ELhashsize - 1) { + if (ELhash[bucket] != (struct Halfedge *)NULL) + ELhash[bucket]->ELrefcnt -= 1; + ELhash[bucket] = he; + ELhash[bucket]->ELrefcnt += 1; + }; + return (he); +} + +/* This delete routine can't reclaim node, since pointers from hash + table may be present. */ +void ELdelete(struct Halfedge *he) { + (he->ELleft)->ELright = he->ELright; + (he->ELright)->ELleft = he->ELleft; + he->ELedge = (struct Edge *)DELETED; +} + +struct Halfedge *ELright(he) struct Halfedge *he; +{ return (he->ELright); } + +struct Halfedge *ELleft(he) struct Halfedge *he; +{ return (he->ELleft); } + +struct Site *leftreg(he) struct Halfedge *he; +{ + if (he->ELedge == (struct Edge *)NULL) + return (bottomsite); + return (he->ELpm == le ? he->ELedge->reg[le] : he->ELedge->reg[re]); +} + +struct Site *rightreg(he) struct Halfedge *he; +{ + if (he->ELedge == (struct Edge *)NULL) + return (bottomsite); + return (he->ELpm == le ? he->ELedge->reg[re] : he->ELedge->reg[le]); +} diff --git a/src/fortune/geometry.c b/src/fortune/geometry.c new file mode 100644 index 0000000..131cacf --- /dev/null +++ b/src/fortune/geometry.c @@ -0,0 +1,184 @@ +# +#include "defs.h" +#include <math.h> + +void geominit() { + struct Edge e; + double sn; + + freeinit(&efl, sizeof e); + nvertices = 0; + nedges = 0; + sn = nsites + 4; + sqrt_nsites = sqrt(sn); + deltay = ymax - ymin; + deltax = xmax - xmin; +} + +struct Edge *bisect(s1, s2) struct Site *s1, *s2; +{ + double dx, dy, adx, ady; + struct Edge *newedge; + + newedge = (struct Edge *)getfree(&efl); + + newedge->reg[0] = s1; + newedge->reg[1] = s2; + ref(s1); + ref(s2); + newedge->ep[0] = (struct Site *)NULL; + newedge->ep[1] = (struct Site *)NULL; + + dx = s2->coord.x - s1->coord.x; + dy = s2->coord.y - s1->coord.y; + adx = dx > 0 ? dx : -dx; + ady = dy > 0 ? dy : -dy; + newedge->c = s1->coord.x * dx + s1->coord.y * dy + (dx * dx + dy * dy) * 0.5; + if (adx > ady) { + newedge->a = 1.0; + newedge->b = dy / dx; + newedge->c /= dx; + } else { + newedge->b = 1.0; + newedge->a = dx / dy; + newedge->c /= dy; + }; + + newedge->edgenbr = nedges; + out_bisector(newedge); + nedges += 1; + return (newedge); +} + +struct Site *intersect(el1, el2) struct Halfedge *el1, *el2; +{ + struct Edge *e1, *e2, *e; + struct Halfedge *el; + double d, xint, yint; + int right_of_site; + struct Site *v; + + e1 = el1->ELedge; + e2 = el2->ELedge; + if (e1 == (struct Edge *)NULL || e2 == (struct Edge *)NULL) + return ((struct Site *)NULL); + if (e1->reg[1] == e2->reg[1]) + return ((struct Site *)NULL); + + d = e1->a * e2->b - e1->b * e2->a; + /* printf("intersect: d=%g\n", d); */ + if (-1.0e-20 < d && d < 1.0e-20) { + return ((struct Site *)NULL); + }; + + xint = (e1->c * e2->b - e2->c * e1->b) / d; + yint = (e2->c * e1->a - e1->c * e2->a) / d; + + if ((e1->reg[1]->coord.y < e2->reg[1]->coord.y) || + (e1->reg[1]->coord.y == e2->reg[1]->coord.y && + e1->reg[1]->coord.x < e2->reg[1]->coord.x)) { + el = el1; + e = e1; + } else { + el = el2; + e = e2; + }; + right_of_site = xint >= e->reg[1]->coord.x; + if ((right_of_site && el->ELpm == le) || (!right_of_site && el->ELpm == re)) + return ((struct Site *)NULL); + + v = (struct Site *)getfree(&sfl); + v->refcnt = 0; + v->coord.x = xint; + v->coord.y = yint; + return (v); +} + +/* returns 1 if p is to right of halfedge e */ +int right_of(el, p) struct Halfedge *el; +struct Point *p; +{ + struct Edge *e; + struct Site *topsite; + int right_of_site, above, fast; + double dxp, dyp, dxs, t1, t2, t3, yl; + + e = el->ELedge; + topsite = e->reg[1]; + right_of_site = p->x > topsite->coord.x; + if (right_of_site && el->ELpm == le) + return (1); + if (!right_of_site && el->ELpm == re) + return (0); + + if (e->a == 1.0) { + dyp = p->y - topsite->coord.y; + dxp = p->x - topsite->coord.x; + fast = 0; + if ((!right_of_site && e->b < 0.0) | (right_of_site && e->b >= 0.0)) { + above = dyp >= e->b * dxp; + fast = above; + } else { + above = p->x + p->y * e->b > e->c; + if (e->b < 0.0) + above = !above; + if (!above) + fast = 1; + }; + if (!fast) { + dxs = topsite->coord.x - (e->reg[0])->coord.x; + above = e->b * (dxp * dxp - dyp * dyp) < + dxs * dyp * (1.0 + 2.0 * dxp / dxs + e->b * e->b); + if (e->b < 0.0) + above = !above; + }; + } else /*e->b==1.0 */ + { + yl = e->c - e->a * p->x; + t1 = p->y - yl; + t2 = p->x - topsite->coord.x; + t3 = yl - topsite->coord.y; + above = t1 * t1 > t2 * t2 + t3 * t3; + }; + return (el->ELpm == le ? above : !above); +} + +void endpoint(e, lr, s) struct Edge *e; +int lr; +struct Site *s; +{ + e->ep[lr] = s; + ref(s); + if (e->ep[re - lr] == (struct Site *)NULL) { + } else { + out_ep(e); + deref(e->reg[le]); + deref(e->reg[re]); + makefree((struct Freenode *)e, &efl); + } +} + +double dist(s, t) struct Site *s, *t; +{ + double dx, dy; + dx = s->coord.x - t->coord.x; + dy = s->coord.y - t->coord.y; + return (sqrt(dx * dx + dy * dy)); +} + +void makevertex(v) struct Site *v; +{ + v->sitenbr = nvertices; + nvertices += 1; + out_vertex(v); +} + +void deref(v) struct Site *v; +{ + v->refcnt -= 1; + if (v->refcnt == 0) + makefree((struct Freenode *)v, &sfl); +} + +void ref(v) struct Site *v; +{ v->refcnt += 1; } diff --git a/src/fortune/heap.c b/src/fortune/heap.c new file mode 100644 index 0000000..2ebbd58 --- /dev/null +++ b/src/fortune/heap.c @@ -0,0 +1,94 @@ +# +#include "defs.h" +int PQhashsize; +struct Halfedge *PQhash; +struct Halfedge *PQfind(); +int PQcount; +int PQmin; +int PQempty(); + +void PQinsert(he, v, offset) struct Halfedge *he; +struct Site *v; +double offset; +{ + struct Halfedge *last, *next; + he->vertex = v; + ref(v); + he->ystar = v->coord.y + offset; + last = &PQhash[PQbucket(he)]; + while ((next = last->PQnext) != (struct Halfedge *)NULL && + (he->ystar > next->ystar || + (he->ystar == next->ystar && v->coord.x > next->vertex->coord.x))) { + last = next; + }; + he->PQnext = last->PQnext; + last->PQnext = he; + PQcount += 1; +} + +void PQdelete(he) struct Halfedge *he; +{ + struct Halfedge *last; + + if (he->vertex != (struct Site *)NULL) { + last = &PQhash[PQbucket(he)]; + while (last->PQnext != he) + last = last->PQnext; + last->PQnext = he->PQnext; + PQcount -= 1; + deref(he->vertex); + he->vertex = (struct Site *)NULL; + }; +} + +int PQbucket(he) struct Halfedge *he; +{ + int bucket; + + if (he->ystar < ymin) + bucket = 0; + else if (he->ystar >= ymax) + bucket = PQhashsize - 1; + else + bucket = (he->ystar - ymin) / deltay * PQhashsize; + if (bucket < 0) + bucket = 0; + if (bucket >= PQhashsize) + bucket = PQhashsize - 1; + if (bucket < PQmin) + PQmin = bucket; + return (bucket); +} + +int PQempty() { return (PQcount == 0); } + +struct Point PQ_min() { + struct Point answer; + + while (PQhash[PQmin].PQnext == (struct Halfedge *)NULL) { + PQmin += 1; + }; + answer.x = PQhash[PQmin].PQnext->vertex->coord.x; + answer.y = PQhash[PQmin].PQnext->ystar; + return (answer); +} + +struct Halfedge *PQextractmin() { + struct Halfedge *curr; + curr = PQhash[PQmin].PQnext; + PQhash[PQmin].PQnext = curr->PQnext; + PQcount -= 1; + return (curr); +} + +void PQinitialize() { + int i; +// struct Point *s; + + PQcount = 0; + PQmin = 0; + PQhashsize = 4 * sqrt_nsites; + PQhash = (struct Halfedge *)myalloc(PQhashsize * sizeof *PQhash); + for (i = 0; i < PQhashsize; i += 1) + PQhash[i].PQnext = (struct Halfedge *)NULL; +} diff --git a/src/fortune/main.c b/src/fortune/main.c new file mode 100644 index 0000000..7571945 --- /dev/null +++ b/src/fortune/main.c @@ -0,0 +1,284 @@ +# +#include <stdio.h> +#include <stdint.h> +#include "defs.h" +#include <math.h> +struct Site *nextone(); + +int triangulate, sorted, plot, debug; +double xmin, xmax, ymin, ymax, deltax, deltay; +int vert_count, edge_count, line_count, dual_count, site_count; +double *site_list, *vert_list, *line_list; +unsigned int *edge_list, *dual_list; + +/* sort sites on y, then x, coord */ +int scomp(s1, s2) struct Point *s1, *s2; +{ + if (s1->y < s2->y) + return (-1); + if (s1->y > s2->y) + return (1); + if (s1->x < s2->x) + return (-1); + if (s1->x > s2->x) + return (1); + return (0); +} + +/* read all sites, sort, and compute xmin, xmax, ymin, ymax */ +void readsites(double *lattice) { + + sites = (struct Site *)myalloc(nsites * sizeof *sites); + for (unsigned int j = 0; j < nsites; j++) { + sites[j].coord.x = lattice[2 * j]; + sites[j].coord.y = lattice[2 * j + 1]; + sites[j].sitenbr = j; + sites[j].refcnt = 0; + } + qsort(sites, nsites, sizeof(struct Site), scomp); + xmin = sites[0].coord.x; + xmax = sites[0].coord.x; + for (unsigned int i = 1; i < nsites; i += 1) { + if (sites[i].coord.x < xmin) + xmin = sites[i].coord.x; + if (sites[i].coord.x > xmax) + xmax = sites[i].coord.x; + }; + ymin = sites[0].coord.y; + ymax = sites[nsites - 1].coord.y; +} + +unsigned int delete_duplicates(unsigned int ne, unsigned int *etv) { + unsigned int ndup = 0; + bool *duplicates = (bool *)calloc(ne, sizeof(bool)); + for (unsigned int i = 0; i < ne; i++) { + unsigned int v1 = etv[2 * i]; + unsigned int v2 = etv[2 * i + 1]; + for (unsigned int j = 0; j < i; j++) { + unsigned int w1 = etv[2 * j]; + unsigned int w2 = etv[2 * j + 1]; + if (w1 == v1 && w2 == v2) { + duplicates[j] = true; + ndup++; + break; + } + } + } + unsigned int newsize = (int)ne - (int)ndup; + unsigned int count = 0; + for (unsigned int i = 0; i < ne; i++) { + if (!duplicates[i]) { + etv[2 * count] = etv[2 * i]; + etv[2 * count + 1] = etv[2 * i + 1]; + count++; + } + } + free(duplicates); + return newsize; +} + +intptr_t *run_voronoi(unsigned int num, double *lattice, double xmin, double xmax, double ymin, double ymax) +{ + struct Site *(*next)(); + + sorted = 0; + triangulate = 0; + plot = 0; + debug = 0; + + alloclist = (char **)malloc(9 * num * sizeof(char *)); + allocnum = 0; + + freeinit(&sfl, sizeof *sites); + + unsigned int eff_num = 9 * num; + double *eff_lattice = (double *)malloc(2 * eff_num * sizeof(double)); + + for (unsigned int i = 0; i < num; i++) { + // original sites - our baby boys + eff_lattice[2*i] = lattice[2*i]; + eff_lattice[2*i+1] = lattice[2*i+1]; + + // sites to the right + eff_lattice[2*(num+i)+1] = lattice[2*i+1] - xmin + xmax; + eff_lattice[2*(num+i)] = lattice[2*i]; + + // sites to the left + eff_lattice[2*(2*num+i)+1] = lattice[2*i+1] - xmax + xmin; + eff_lattice[2*(2*num+i)] = lattice[2*i]; + + // sites to the top + eff_lattice[2*(3*num+i)+1] = lattice[2*i+1]; + eff_lattice[2*(3*num+i)] = lattice[2*i] - ymin + ymax; + + // sites to the bottom + eff_lattice[2*(4*num+i)+1] = lattice[2*i+1]; + eff_lattice[2*(4*num+i)] = lattice[2*i] - ymax + ymin; + + // sites to the upper right + eff_lattice[2*(5*num+i)+1] = lattice[2*i+1] - xmin + xmax; + eff_lattice[2*(5*num+i)] = lattice[2*i] - ymin + ymax; + + // sites to the upper left + eff_lattice[2*(6*num+i)+1] = lattice[2*i+1] - xmax + xmin; + eff_lattice[2*(6*num+i)] = lattice[2*i] - ymin + ymax; + + // sites to the lower left + eff_lattice[2*(7*num+i)+1] = lattice[2*i+1] - xmax + xmin; + eff_lattice[2*(7*num+i)] = lattice[2*i] + ymin - ymax; + + // sites to the lower right + eff_lattice[2*(8*num+i)+1] = lattice[2*i+1] + xmax - xmin; + eff_lattice[2*(8*num+i)] = lattice[2*i] + ymin - ymax; + } + + nsites = eff_num; + readsites(eff_lattice); + free(eff_lattice); + next = nextone; + + siteidx = 0; + geominit(); + + site_list = malloc(2 * nsites * sizeof(double)); + vert_list = NULL; + line_list = NULL; + edge_list = NULL; + dual_list = NULL; + site_count = 0; + vert_count = 0; + edge_count = 0; + line_count = 0; + dual_count = 0; + + voronoi(triangulate, next); + + unsigned int real_vert_count = vert_count; + unsigned int real_edge_count = edge_count; + unsigned int real_dual_count = dual_count; + double *real_vert_list = vert_list; + unsigned int *real_edge_list = edge_list; + unsigned int *real_dual_list = dual_list; + + { + real_vert_count = 0; + real_edge_count = 0; + real_dual_count = 0; + real_vert_list = (double *)malloc(2 * vert_count * sizeof(double)); + real_edge_list = (unsigned int *)malloc(2 * edge_count * sizeof(unsigned int)); + real_dual_list = (unsigned int *)malloc(3 * dual_count * sizeof(unsigned int)); + unsigned int *triangle_analyzed = (unsigned int *)malloc(4 * dual_count * sizeof(unsigned int)); + unsigned int q = 0; + int *new_vert_ind = (int *)malloc(vert_count * sizeof(int)); + + for (unsigned int i = 0; i < vert_count; i++) { + unsigned int t1, t2, t3; + t1 = dual_list[3*i]; t2 = dual_list[3*i+1]; t3 = dual_list[3*i+2]; + if (t1 >= num && t2 >= num && t3 >= num) { + new_vert_ind[i] = -1; + } + else if (t1 < num && t2 < num && t3 < num) { + real_vert_list[2*real_vert_count] = vert_list[2*i]; + real_vert_list[2*real_vert_count+1] = vert_list[2*i+1]; + real_dual_list[3*real_vert_count] = t1; + real_dual_list[3*real_vert_count+1] = t2; + real_dual_list[3*real_vert_count+2] = t3; + new_vert_ind[i] = real_vert_count; + real_vert_count++; + } + else { + unsigned int tt1, tt2, tt3; + tt1 = t1 % num; tt2 = t2 % num; tt3 = t3 % num; + unsigned int s1, s2, s3; + s1 = tt1 < tt2 ? (tt1 < tt3 ? tt1 : tt3) : (tt2 < tt3 ? tt2 : tt3); + s3 = tt1 > tt2 ? (tt1 > tt3 ? tt1 : tt3) : (tt2 > tt3 ? tt2 : tt3); + s2 = tt1 < tt2 ? (tt1 > tt3 ? tt1 : (tt2 < tt3 ? tt2 : tt3)) : (tt1 < tt3 ? tt1 : (tt2 < tt3 ? tt3: tt2)); + + bool matched = false; + unsigned int ass_vert; + + for (unsigned int j = 0; j < q; j++) { + if (s1 == triangle_analyzed[4*j] && s2 == triangle_analyzed[4*j+1] && s3 == triangle_analyzed[4*j+2]) { + matched = true; + ass_vert = triangle_analyzed[4*j+3]; + break; + } + } + + if (matched) { + new_vert_ind[i] = ass_vert; + } else { + triangle_analyzed[4*q] = s1; + triangle_analyzed[4*q+1] = s2; + triangle_analyzed[4*q+2] = s3; + triangle_analyzed[4*q+3] = real_vert_count; + q++; + real_vert_list[2*real_vert_count+1] = fmod(2*xmax+vert_list[2*i+1], xmax); + real_vert_list[2*real_vert_count] = fmod(2*ymax+vert_list[2*i], ymax); + real_dual_list[3*real_vert_count] = t1 % num; + real_dual_list[3*real_vert_count+1] = t2 % num; + real_dual_list[3*real_vert_count+2] = t3 % num; + new_vert_ind[i] = real_vert_count; + real_vert_count++; + } + } + } + for (unsigned int i = 0; i < edge_count; i++) { + unsigned int v1, v2; + v1 = edge_list[2 * i]; + v2 = edge_list[2 * i+1]; + if (v1 != UINT_MAX && v2 != UINT_MAX) { + if (new_vert_ind[v1] >= 0 && new_vert_ind[v2] >=0) { + unsigned int nv1 = new_vert_ind[v1]; + unsigned int nv2 = new_vert_ind[v2]; + real_edge_list[2*real_edge_count] = nv1 < nv2 ? nv1 : nv2; + real_edge_list[2*real_edge_count+1] = nv1 < nv2 ? nv2 : nv1; + real_edge_count++; + } + } + } + unsigned int new_edge_count = delete_duplicates(real_edge_count, real_edge_list); + real_edge_count = new_edge_count; + free(triangle_analyzed); + free(new_vert_ind); + free(dual_list); + free(edge_list); + free(vert_list); + } + + + intptr_t *output = (intptr_t *)malloc(6 * sizeof(intptr_t)); + output[0] = (intptr_t)malloc(4 * sizeof(unsigned int)); + ((unsigned int *)output[0])[0] = real_vert_count; + ((unsigned int *)output[0])[1] = real_edge_count; + ((unsigned int *)output[0])[2] = line_count; + ((unsigned int *)output[0])[3] = real_vert_count; + output[1] = (intptr_t)site_list; + output[2] = (intptr_t)real_vert_list; + output[3] = (intptr_t)real_edge_list; + output[4] = (intptr_t)line_list; + output[5] = (intptr_t)real_dual_list; + + free(ELhash); + free(PQhash); + free(sites); + for (unsigned int i = 0; i < allocnum; i++) { + free(alloclist[i]); + } + free(alloclist); + + return output; +} + +/* return a single in-storage site */ +struct Site *nextone() { + for (; siteidx < nsites; siteidx += 1) { + if (siteidx == 0 || sites[siteidx].coord.x != sites[siteidx - 1].coord.x || + sites[siteidx].coord.y != sites[siteidx - 1].coord.y) { + siteidx += 1; + return (&sites[siteidx - 1]); + }; + }; + return ((struct Site *)NULL); +} + diff --git a/src/fortune/memory.c b/src/fortune/memory.c new file mode 100644 index 0000000..807d1b0 --- /dev/null +++ b/src/fortune/memory.c @@ -0,0 +1,44 @@ +# +#include "defs.h" +#include <stdio.h> + +void freeinit(struct Freelist *fl, int size) { + fl->head = (struct Freenode *)NULL; + fl->nodesize = size; +} + +char *getfree(fl) struct Freelist *fl; +{ + int i; + struct Freenode *t; + if (fl->head == (struct Freenode *)NULL) { + t = (struct Freenode *)myalloc(sqrt_nsites * fl->nodesize); + alloclist[allocnum] = t; + allocnum++; + for (i = 0; i < sqrt_nsites; i += 1) + makefree((struct Freenode *)((char *)t + i * fl->nodesize), fl); + }; + t = fl->head; + fl->head = (fl->head)->nextfree; + return ((char *)t); +} + +void makefree(curr, fl) struct Freenode *curr;struct Freelist *fl; +{ + curr->nextfree = fl->head; + fl->head = curr; +} + +int total_alloc; +char *myalloc(n) unsigned n; +{ + char *t; + if ((t = malloc(n)) == (char *)0) { + fprintf(stderr, + "Insufficient memory processing site %d (%d bytes in use)\n", + siteidx, total_alloc); + exit(1); + }; + total_alloc += n; + return (t); +} diff --git a/src/fortune/output.c b/src/fortune/output.c new file mode 100644 index 0000000..d496feb --- /dev/null +++ b/src/fortune/output.c @@ -0,0 +1,46 @@ +# +#include "defs.h" +#include <stdio.h> +double pxmin, pxmax, pymin, pymax, cradius; + +void out_bisector(e) struct Edge *e; +{ + if (line_count % nsites == 0) line_list = realloc(line_list, 3 * (nsites + line_count) * sizeof(double)); + line_list[3*line_count] = e->a; + line_list[3*line_count+1] = e->b; + line_list[3*line_count+2] = e->c; + line_count++; +} + +void out_ep(e) struct Edge *e; +{ + if (edge_count % nsites == 0) edge_list = realloc(edge_list, 2 * (nsites + edge_count) * sizeof(unsigned int)); + edge_list[2 * edge_count] = (e->ep[le] != (struct Site *)NULL ? e->ep[le]->sitenbr : UINT_MAX); + edge_list[2 * edge_count + 1] = (e->ep[re] != (struct Site *)NULL ? e->ep[re]->sitenbr : UINT_MAX); + edge_count++; +} + +void out_vertex(v) struct Site *v; +{ + if (vert_count % nsites == 0) vert_list = realloc(vert_list, 2 * (nsites + vert_count) * sizeof(double)); + vert_list[2 * vert_count] = v->coord.x; + vert_list[2 * vert_count + 1] = v->coord.y; + vert_count++; +} + +void out_site(s) struct Site *s; +{ + site_list[2 * site_count] = s->coord.x; + site_list[2 * site_count + 1] = s->coord.y; + site_count++; +} + +void out_triple(s1, s2, s3) struct Site *s1, *s2, *s3; +{ + if (dual_count % nsites == 0) dual_list = realloc(dual_list, 3 * (nsites + dual_count) * sizeof(unsigned int)); + dual_list[dual_count * 3] = s1->sitenbr; + dual_list[dual_count * 3 + 1] = s2->sitenbr; + dual_list[dual_count * 3 + 2] = s3->sitenbr; + dual_count++; +} + diff --git a/src/fortune/voronoi.c b/src/fortune/voronoi.c new file mode 100644 index 0000000..dc9945f --- /dev/null +++ b/src/fortune/voronoi.c @@ -0,0 +1,103 @@ +# +#include "defs.h" + +struct Site *sites; +int nsites; +int siteidx; +int sqrt_nsites; +int nvertices; +struct Freelist sfl; +struct Site *bottomsite; + +/* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax, + deltax, deltay (can all be estimates). + Performance suffers if they are wrong; better to make nsites, + deltax, and deltay too big than too small. (?) */ + +void voronoi(triangulate, nextsite) int triangulate; +struct Site *(*nextsite)(); +{ + struct Site *newsite, *bot, *top, *temp, *p; + struct Site *v; + struct Point newintstar; + int pm; + struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector; + struct Edge *e; + + PQinitialize(); + bottomsite = (*nextsite)(); + out_site(bottomsite); + ELinitialize(); + + newsite = (*nextsite)(); + while (1) { + if (!PQempty()) + newintstar = PQ_min(); + + if (newsite != (struct Site *)NULL && + (PQempty() || newsite->coord.y < newintstar.y || + (newsite->coord.y == newintstar.y && + newsite->coord.x < newintstar.x))) {/* new site is smallest */ + out_site(newsite); + lbnd = ELleftbnd(&(newsite->coord)); + rbnd = ELright(lbnd); + bot = rightreg(lbnd); + e = bisect(bot, newsite); + bisector = HEcreate(e, le); + ELinsert(lbnd, bisector); + if ((p = intersect(lbnd, bisector)) != (struct Site *)NULL) { + PQdelete(lbnd); + PQinsert(lbnd, p, dist(p, newsite)); + }; + lbnd = bisector; + bisector = HEcreate(e, re); + ELinsert(lbnd, bisector); + if ((p = intersect(bisector, rbnd)) != (struct Site *)NULL) { + PQinsert(bisector, p, dist(p, newsite)); + }; + newsite = (*nextsite)(); + } else if (!PQempty()) + /* intersection is smallest */ + { + lbnd = PQextractmin(); + llbnd = ELleft(lbnd); + rbnd = ELright(lbnd); + rrbnd = ELright(rbnd); + bot = leftreg(lbnd); + top = rightreg(rbnd); + out_triple(bot, top, rightreg(lbnd)); + v = lbnd->vertex; + makevertex(v); + endpoint(lbnd->ELedge, lbnd->ELpm, v); + endpoint(rbnd->ELedge, rbnd->ELpm, v); + ELdelete(lbnd); + PQdelete(rbnd); + ELdelete(rbnd); + pm = le; + if (bot->coord.y > top->coord.y) { + temp = bot; + bot = top; + top = temp; + pm = re; + } + e = bisect(bot, top); + bisector = HEcreate(e, pm); + ELinsert(llbnd, bisector); + endpoint(e, re - pm, v); + deref(v); + if ((p = intersect(llbnd, bisector)) != (struct Site *)NULL) { + PQdelete(llbnd); + PQinsert(llbnd, p, dist(p, bot)); + }; + if ((p = intersect(bisector, rrbnd)) != (struct Site *)NULL) { + PQinsert(bisector, p, dist(p, bot)); + }; + } else + break; + }; + + for (lbnd = ELright(ELleftend); lbnd != ELrightend; lbnd = ELright(lbnd)) { + e = lbnd->ELedge; + out_ep(e); + }; +} diff --git a/src/fracture.c b/src/fracture.c new file mode 100644 index 0000000..6798231 --- /dev/null +++ b/src/fracture.c @@ -0,0 +1,578 @@ + + +#include "fracture.h" + +int main(int argc, char *argv[]) { + int opt; + + // defining variables to be (potentially) set by command line flags + unsigned int num, width, filename_len; + double crack_len, crack_width, beta, inf, cutoff; + boundary_type periodic; + bool include_breaking, supplied_bound, save_clusters, voronoi, + side_bounds, voltage_bound, repeat_voronoi, network_out, save_avg_stress, save_avg_ztress, + save_app_stress, save_corr, save_cond, use_crit, use_first, save_damage, save_homo_damage; + char *bound_filename; + + filename_len = 100; + + num = 100; + width = 16; + crack_len = 16; + beta = .3; + inf = 1e-8; + cutoff = 1e-8; + periodic = FREE_BOUND; + include_breaking = false; + supplied_bound = false; + save_clusters = false; + voronoi = false; + side_bounds = false; + voltage_bound = false; + repeat_voronoi = false; + network_out = false; + save_avg_stress = false; + save_avg_ztress = false; + save_app_stress = false; + save_damage = false; + save_homo_damage = false; + save_corr = false; + save_cond = false; + use_crit = true; + use_first = false; + + int periodic_int; + + while ((opt = getopt(argc, argv, "n:w:b:l:f:ocp:vVsrNCaSFtzdH")) != -1) { + switch (opt) { + case 'n': + num = atoi(optarg); + break; + case 'w': + width = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + case 'l': + crack_len = atof(optarg); + break; + case 'p': + periodic_int = atoi(optarg); + switch (periodic_int) { + case 0: + periodic = FREE_BOUND; + break; + case 1: + periodic = CYLINDER_BOUND; + break; + case 2: + periodic = TORUS_BOUND; + break; + } + break; + case 'C': + save_avg_stress = true; + use_crit = true; + break; + case 'd': + save_damage = true; + break; + case 'H': + save_homo_damage = true; + break; + case 'z': + save_avg_ztress = true; + use_crit = true; + break; + case 'v': + voronoi = true; + break; + case 'F': + use_first = true; + break; + case 'r': + repeat_voronoi = true; + break; + case 'V': + voltage_bound = true; + break; + case 's': + side_bounds = true; + break; + case 'c': + save_clusters = true; + break; + case 'o': + include_breaking = true; + break; + case 'N': + network_out = true; + break; + case 'a': + save_app_stress = true; + break; + case 'S': + save_corr = true; + break; + case 't': + save_cond = true; + inf = 1; + break; + case 'f': + supplied_bound = true; + bound_filename = optarg; + break; + default: /* '?' */ + exit(EXIT_FAILURE); + } + } + + FILE *break_out; + if (include_breaking) { + char *break_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(break_filename, filename_len, "breaks_%u_%g_%g_%u.txt", width, + crack_len, beta, num); + break_out = fopen(break_filename, "w"); + free(break_filename); + } + + // start cholmod + cholmod_common c; + CHOL_F(start)(&c); + + /* if we use voltage boundary conditions, the laplacian matrix is positive + * definite and we can use a supernodal LL decomposition. otherwise we need + * to use the simplicial LDL decomposition + */ + if (voltage_bound) { + (&c)->supernodal = CHOLMOD_SUPERNODAL; + } else { + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + } + + fnet *network; + finst *perm_instance; + unsigned int c_dist_size; + unsigned int a_dist_size; + unsigned int voronoi_num = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2); + + if (voronoi) { + crack_width = 1. / (2 * width); + c_dist_size = 2 * voronoi_num; + a_dist_size = 2 * voronoi_num; + if (repeat_voronoi) { + while ((network = ini_voronoi_network(width, periodic, + genfunc_uniform, &c)) == NULL) + ; + perm_instance = create_instance(network, inf, voltage_bound, false, &c); + gen_crack(perm_instance, crack_len, crack_width, &c); + finish_instance(perm_instance, &c); + } + } else { + network = ini_square_network(width, periodic, side_bounds, &c); + crack_width = 1; + perm_instance = create_instance(network, inf, voltage_bound, false, &c); + gen_crack(perm_instance, crack_len, crack_width, &c); + finish_instance(perm_instance, &c); + c_dist_size = network->num_dual_verts; + a_dist_size = network->num_verts; + } + + // define arrays for saving cluster and avalanche distributions + unsigned int *cluster_size_dist; + unsigned int *avalanche_size_dist; + char *c_filename; + if (save_clusters) { + cluster_size_dist = + (unsigned int *)calloc(c_dist_size, sizeof(unsigned int)); + avalanche_size_dist = + (unsigned int *)calloc(a_dist_size, sizeof(unsigned int)); + + c_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(c_filename, filename_len, "cluster_%d_%g_%g.txt", width, crack_len, + beta); + + FILE *cluster_out = fopen(c_filename, "r"); + + if (cluster_out != NULL) { + for (unsigned int i = 0; i < c_dist_size; i++) { + unsigned int tmp; + fscanf(cluster_out, "%u ", &tmp); + cluster_size_dist[i] = tmp; + } + fscanf(cluster_out, "\n"); + for (unsigned int i = 0; i < a_dist_size; i++) { + unsigned int tmp; + fscanf(cluster_out, "%u ", &tmp); + avalanche_size_dist[i] = tmp; + } + fclose(cluster_out); + } + } + + double *app_stress; + double *app_stress_crit; + if (save_app_stress) { + app_stress = (double *)malloc(num * sizeof(double)); + app_stress_crit = (double *)malloc(num * sizeof(double)); + } + + double *avg_corr; + unsigned int **dists = NULL; + if (save_corr) { + avg_corr = (double *)calloc(2 * voronoi_num, sizeof(double)); + if (!voronoi) + dists = get_dists(network); + } + + double *conductivity; + if (save_cond) { + conductivity = (double *)malloc(num * sizeof(double)); + } + + double *avg_stress; + unsigned int *num_stress; + if (save_avg_stress) { + avg_stress = (double *)calloc(pow(width, 2), sizeof(double)); + num_stress = (unsigned int *)calloc(pow(width, 2), sizeof(unsigned int)); + } + double *avg_ztress; + if (save_avg_ztress) { + avg_ztress = (double *)calloc(pow(width, 2), sizeof(double)); + } + + double *damage; + if (save_damage) { + damage = (double *)calloc(pow(width, 2), sizeof(double)); + } + + double *homo_damage; + if (save_homo_damage) { + homo_damage = (double *)calloc(num, sizeof(double)); + } + + double *square_bound; + unsigned int square_num_verts = 2 * ((width + 1) / 2) * (width / 2 + 1); + if (supplied_bound) { + FILE *bound_file = fopen(bound_filename, "r"); + square_bound = (double *)calloc(square_num_verts, sizeof(double)); + for (unsigned int i = 0; i < square_num_verts; i++) { + double tmp; + fscanf(bound_file, "%lg ", &tmp); + square_bound[i] = tmp; + } + fclose(bound_file); + if (!voronoi) { + double total = 0; + for (unsigned int i = 0; i < square_num_verts; i++) { + ((double *)perm_instance->boundary_cond->x)[i] = square_bound[i]; + total += square_bound[i]; + } + ((double *)perm_instance->boundary_cond->x)[square_num_verts] = 0; + ((double *)perm_instance->boundary_cond->x)[square_num_verts + 1] = 0; + } + } + + printf("\n"); + for (int DUMB = 0; DUMB < num; DUMB++) { + printf("\033[F\033[JFRACTURE: %0*d / %d\n", (int)log10(num) + 1, DUMB + 1, + num); + + break_data *breaking_data = NULL; + while (breaking_data == NULL) { + if (voronoi && !repeat_voronoi) { + while ((network = ini_voronoi_network(width, periodic, + genfunc_uniform, &c)) == NULL) + ; + perm_instance = create_instance(network, inf, voltage_bound, false, &c); + gen_crack(perm_instance, crack_len, crack_width, &c); + finish_instance(perm_instance, &c); + if (supplied_bound) { + voronoi_bound_ini(perm_instance, square_bound, width); + } + } + double *fuse_thres = gen_fuse_thres( + network->num_edges, network->edge_coords, beta, beta_scaling_flat); + finst *instance = copy_instance(perm_instance, &c); + breaking_data = fracture_network(instance, fuse_thres, &c, cutoff); + free_instance(instance, &c); + free(fuse_thres); + } + + unsigned int min_pos = 0; + double min_val = DBL_MAX; + + for (unsigned int j = 0; j < breaking_data->num_broken; j++) { + double val = fabs(breaking_data->extern_field[j]); + + + if (val < min_val) { + min_pos = j; + min_val = val; + } + } + + if (save_app_stress) { + app_stress[DUMB] = + 1 / fabs(breaking_data->extern_field[breaking_data->num_broken - 1]); + app_stress_crit[DUMB] = 1 / fabs(breaking_data->extern_field[min_pos]); + } + + finst *tmp_instance = copy_instance(perm_instance, &c); + + int stop_at = breaking_data->num_broken; + if (use_crit) + stop_at = min_pos; + if (use_first) + stop_at = 0; + + unsigned int av_size = 0; + double cur_val = DBL_MAX; + for (unsigned int i = 0; i < min_pos; i++) { + double val = fabs(breaking_data->extern_field[i]); + if (save_clusters) { + if (val > cur_val) { + av_size++; + } + + if (val < cur_val) { + avalanche_size_dist[av_size]++; + av_size = 0; + cur_val = val; + } + } + } + + for (unsigned int i = 0; i < stop_at; i++) { + break_edge(tmp_instance, breaking_data->break_list[i], &c); + } + + if (save_cond) { + double *tmp_voltage = get_voltage(tmp_instance, &c); + conductivity[DUMB] = fabs(tmp_voltage[tmp_instance->network->num_verts + 1] - tmp_voltage[tmp_instance->network->num_verts]); + free(tmp_voltage); + } + + if (save_homo_damage) { + homo_damage[DUMB] = ((double)min_pos) / tmp_instance->network->num_edges; + } + + if (save_avg_stress || save_avg_ztress) { + double *tmp_stress = get_current(tmp_instance, &c); + if (voronoi) { + double *tmp_stress_2 = + bin_values(tmp_instance->network, width, tmp_stress); + free(tmp_stress); + tmp_stress = tmp_stress_2; + } + for (unsigned int i = 0; i < pow(width, 2); i++) { + double sigma = tmp_stress[i]; + if (sigma != 0 && sigma == sigma && save_avg_stress) { + avg_stress[i] += sigma; + num_stress[i]++; + } + if (sigma == sigma && save_avg_ztress) { + avg_ztress[i] += sigma; + } + } + free(tmp_stress); + } + + if (save_damage) { + double *tmp_damage = (double *)calloc(tmp_instance->network->num_edges, sizeof(double)); + for (unsigned int i = 0; i < stop_at; i++) { + tmp_damage[breaking_data->break_list[i]] += 1; + } + if (voronoi) { + double *tmp_damage_2 = bin_values(tmp_instance->network, width, tmp_damage); + free(tmp_damage); + tmp_damage = tmp_damage_2; + } + for (unsigned int i = 0; i < pow(width, 2); i++) { + damage[i] += tmp_damage[i]; + } + free(tmp_damage); + } + + if (save_clusters) { + unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c); + for (unsigned int i = 0; i < network->num_dual_verts; i++) { + cluster_size_dist[i] += tmp_cluster_dist[i]; + } + free(tmp_cluster_dist); + } + + if (save_corr) { + double *tmp_corr = get_corr(tmp_instance, dists, &c); + for (unsigned int i = 0; i < tmp_instance->network->num_dual_verts; i++) { + avg_corr[i] += tmp_corr[i] / num; + } + free(tmp_corr); + } + + if (network_out) { + FILE *net_out = fopen("network.txt", "w"); + for (unsigned int i = 0; i < network->num_verts; i++) { + fprintf(net_out, "%f %f ", network->vert_coords[2 * i], + network->vert_coords[2 * i + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int i = 0; i < network->num_edges; i++) { + fprintf(net_out, "%u %u ", network->edges_to_verts[2 * i], + network->edges_to_verts[2 * i + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int i = 0; i < network->num_dual_verts; i++) { + fprintf(net_out, "%f %f ", network->dual_vert_coords[2 * i], + network->dual_vert_coords[2 * i + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int i = 0; i < network->num_edges; i++) { + fprintf(net_out, "%u %u ", network->dual_edges_to_verts[2 * i], + network->dual_edges_to_verts[2 * i + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int i = 0; i < network->num_edges; i++) { + fprintf(net_out, "%d ", tmp_instance->fuses[i]); + } + fclose(net_out); + } + + free_instance(tmp_instance, &c); + if (voronoi && !repeat_voronoi) { + free_fnet(network, &c); + free_instance(perm_instance, &c); + } + if (include_breaking) { + for (unsigned int i = 0; i < breaking_data->num_broken; i++) { + fprintf(break_out, "%u %f ", breaking_data->break_list[i], + breaking_data->extern_field[i]); + } + fprintf(break_out, "\n"); + } + free_break_data(breaking_data); + } + + printf("\033[F\033[JFRACTURE: COMPLETE\n"); + + if (save_clusters) { + FILE *cluster_out = fopen(c_filename, "w"); + + for (int i = 0; i < c_dist_size; i++) { + fprintf(cluster_out, "%u ", cluster_size_dist[i]); + } + fprintf(cluster_out, "\n"); + for (int i = 0; i < a_dist_size; i++) { + fprintf(cluster_out, "%u ", avalanche_size_dist[i]); + } + fclose(cluster_out); + + free(cluster_size_dist); + free(avalanche_size_dist); + } + + if (save_corr) { + char *corr_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(corr_filename, filename_len, "corr_%d_%g_%g.txt", width, crack_len, + beta); + FILE *corr_file = fopen(corr_filename, "w"); + for (unsigned int i = 0; i < 2 * voronoi_num; i++) { + fprintf(corr_file, "%g ", avg_corr[i]); + } + fclose(corr_file); + free(corr_filename); + } + + if (save_cond) { + char *cond_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(cond_filename, filename_len, "conductivity_%d_%g_%g.txt", width, crack_len, beta); + FILE *cond_file = fopen(cond_filename, "a"); + for (unsigned int i = 0; i < num; i++) { + fprintf(cond_file, "%g ", conductivity[i]); + } + fclose(cond_file); + free(cond_filename); + free(conductivity); + } + + if (save_homo_damage) { + char *hdam_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(hdam_filename, filename_len, "damagep_%d_%g.txt", width, beta); + FILE *hdam_file = fopen(hdam_filename, "a"); + for (unsigned int i = 0; i < num; i++) { + fprintf(hdam_file, "%g ", homo_damage[i]); + } + fclose(hdam_file); + free(hdam_filename); + free(homo_damage); + } + + if (!voronoi || repeat_voronoi) { + free_instance(perm_instance, &c); + //free_fnet(network, &c); + } + + if (include_breaking) { + fclose(break_out); + } + + if (save_avg_stress) { + char *stress_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(stress_filename, filename_len, "current_%d_%g_%g_%u.txt", width, + crack_len, beta, num); + FILE *stress_file = fopen(stress_filename, "w"); + for (unsigned int i = 0; i < pow(width, 2); i++) { + if (num_stress[i] != 0) + avg_stress[i] /= num_stress[i]; + fprintf(stress_file, "%g ", avg_stress[i]); + } + fclose(stress_file); + free(stress_filename); + free(avg_stress); + } + if (save_avg_ztress) { + char *ztress_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(ztress_filename, filename_len, "zurrent_%d_%g_%g_%u.txt", width, + crack_len, beta, num); + FILE *stress_file = fopen(ztress_filename, "w"); + for (unsigned int i = 0; i < pow(width, 2); i++) { + fprintf(stress_file, "%g ", avg_ztress[i] / num); + } + fclose(stress_file); + free(ztress_filename); + free(avg_ztress); + } + + if (save_app_stress) { + char *a_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(a_filename, filename_len, "astress_%d_%g_%g.txt", width, crack_len, + beta); + FILE *a_file = fopen(a_filename, "a"); + for (int i = 0; i < num; i++) { + fprintf(a_file, "%g %g\n", app_stress[i], app_stress_crit[i]); + } + fclose(a_file); + free(a_filename); + free(app_stress); + free(app_stress_crit); + } + + if (save_damage) { + char *damage_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(damage_filename, filename_len, "damage_%d_%g_%g_%u.txt", width, + crack_len, beta, num); + FILE *damage_file = fopen(damage_filename, "w"); + for (unsigned int i = 0; i < pow(width, 2); i++) { + fprintf(damage_file, "%g ", damage[i] / num); + } + fclose(damage_file); + free(damage_filename); + free(damage); + + } + + CHOL_F(finish)(&c); + + return 0; +} diff --git a/src/fracture.h b/src/fracture.h new file mode 100644 index 0000000..4ee7050 --- /dev/null +++ b/src/fracture.h @@ -0,0 +1,182 @@ + +#pragma once + +#include <assert.h> +#include <cholmod.h> +#include <float.h> +#include <getopt.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_randist.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_sf_erf.h> +#include <gsl/gsl_sf_exp.h> +#include <gsl/gsl_sf_log.h> +#include <inttypes.h> +#include <math.h> +#include <omp.h> +#include <stdbool.h> +#include <stdint.h> +#include <stdio.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <sys/types.h> +#include <unistd.h> + +// these two defs allow me to switch to long int cholmod in a sitch +#define CHOL_INT long int +#define CHOL_F(x) cholmod_l_##x + +typedef enum boundary_type { + FREE_BOUND, + CYLINDER_BOUND, + TORUS_BOUND +} boundary_type; + +typedef struct { + unsigned int num_edges; + unsigned int num_verts; + unsigned int num_verts_break; + unsigned int num_bounds; + unsigned int *edges_to_verts; + unsigned int *edges_to_verts_break; + unsigned int *verts_to_edges_ind; + unsigned int *verts_to_edges; + unsigned int *bound_inds; + unsigned int *bound_verts; + double *vert_coords; + double *edge_coords; + unsigned int *dual_edges_to_verts; + unsigned int *dual_verts_to_edges_ind; + unsigned int *dual_verts_to_edges; + double *dual_vert_coords; + unsigned int num_spanning_edges; + unsigned int *spanning_edges; + double L; + unsigned int num_dual_verts; + unsigned int break_dim; + cholmod_sparse *voltcurmat; + boundary_type boundary; +} fnet; + +typedef struct { + fnet *network; + unsigned int num_remaining_edges; + bool *fuses; + double inf; + cholmod_dense *boundary_cond; + cholmod_factor *factor; + unsigned int *marks; + unsigned int *dual_marks; + bool voltage_bound; + unsigned int num_components; + cholmod_sparse *adjacency; + bool debug_stop; +} finst; + +typedef struct { + unsigned int num_broken; + unsigned int *break_list; + double *conductivity; + double *extern_field; +} break_data; + +intptr_t *run_voronoi(unsigned int num_coords, double *coords, double xmin, double xmax, double ymin, double ymax); + +int update_components(const cholmod_sparse *laplacian, unsigned int *marks, + int old_num_components, int v1, int v2, int exclude); + +unsigned int *find_components(const cholmod_sparse *laplacian, unsigned int skip); + +cholmod_sparse *gen_adjacency(const finst *instance, bool dual, bool breakv, + unsigned int pad, cholmod_common *c); + +cholmod_sparse *gen_laplacian(const finst *instance, cholmod_common *c, + bool symmetric); + +int edge_to_verts(unsigned int width, bool periodic, unsigned int edge, + bool index); + +int dual_edge_to_verts(unsigned int width, bool periodic, unsigned int edge, + bool index); + +double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert, + bool index); + +bool update_factor(cholmod_factor *factor, unsigned int v1, unsigned int v2, + cholmod_common *c); + +break_data *fracture_network(finst *instance, double *fuse_thres, + cholmod_common *c, double cutoff); + +double *get_current(const finst *instance, cholmod_common *c); +double *get_current_v(const finst *instance, double *voltages, cholmod_common *c); +double *get_voltage(const finst *instance, cholmod_common *c); + +double *gen_fuse_thres(unsigned int num_edges, double *edge_coords, double beta, + double (*beta_scale)(double, double, double)); + +bool gen_crack(finst *instance, double crack_len, double crack_width, + cholmod_common *c); + +void update_boundary(finst *instance, const double *avg_field); + +FILE *get_file(const char *prefix, unsigned int width, unsigned int crack, + double beta, unsigned int iter, unsigned int num_iter, + unsigned int num, bool read); + +double update_beta(double beta, unsigned int width, const double *stress, + const double *damage, double bound_total); + +cholmod_sparse *gen_voltcurmat(unsigned int num_edges, unsigned int num_verts, + unsigned int *edges_to_verts, cholmod_common *c); + +finst *copy_instance(const finst *instance, cholmod_common *c); + +fnet *ini_square_network(unsigned int width, boundary_type boundary, bool side_bounds, + cholmod_common *c); + +void free_fnet(fnet *network, cholmod_common *c); +void free_instance(finst *instance, cholmod_common *c); + +finst *create_instance(fnet *network, double inf, bool voltage_bound, + bool startnow, cholmod_common *c); + +fnet *ini_voronoi_network(unsigned int L, boundary_type boundary, + double *(*genfunc)(unsigned int, gsl_rng *, unsigned int *), + cholmod_common *c); + +bool check_instance(const finst *instance, cholmod_common *c); + +bool break_edge(finst *instance, unsigned int edge, cholmod_common *c); + +void finish_instance(finst *instance, cholmod_common *c); + +finst *coursegrain_square(finst *instance, fnet *network_p, cholmod_common *c); + +unsigned int *get_clusters(finst *instance, cholmod_common *c); + +unsigned int *get_cluster_dist(finst *instance, cholmod_common *c); + +double *genfunc_uniform(unsigned int L, gsl_rng *r, unsigned int *num); +double *genfunc_hyperuniform(unsigned int L, gsl_rng *r, unsigned int *num); +void randfunc_flat(gsl_rng *r, double *x, double *y); +void randfunc_gaus(gsl_rng *r, double *x, double *y); +double beta_scaling_flat(double beta, double x, double y); +double beta_scaling_gaus(double beta, double x, double y); +double beta_mag(double beta); + +unsigned int *dijkstra(fnet *network, unsigned int source); +unsigned int **get_dists(fnet *network); +double *get_corr(finst *instance, unsigned int **dists, cholmod_common *c); + +double *bin_values(fnet *network, unsigned int width, double *values); + +void voronoi_bound_ini(finst *instance, double *square_bound, + unsigned int width); + +break_data *alloc_break_data(unsigned int num_edges); +void free_break_data(break_data *data); +void update_break_data(break_data *data, unsigned int last_broke, double strength, double conductivity); + +double get_conductivity(finst *inst, double *current, cholmod_common *c); diff --git a/src/fracture_network.c b/src/fracture_network.c new file mode 100644 index 0000000..148b08c --- /dev/null +++ b/src/fracture_network.c @@ -0,0 +1,61 @@ + +#include "fracture.h" + +int inc_break_fuses(finst *instance, double *thres, double *field, + double cutoff) { + unsigned int size = (instance->network)->num_edges; + + int min_pos = -1; + long double min_val = -1; + + for (unsigned int i = 0; i < size; i++) { + if (!instance->fuses[i] && fabs(field[i]) > cutoff) { + double val = fabs(field[i]) / thres[i]; + if (min_val < val) { + min_val = val; min_pos = i; + } + } + } + + return min_pos; +} + +break_data *fracture_network(finst *instance, double *fuse_thres, + cholmod_common *c, double cutoff) { + unsigned int num_edges = instance->network->num_edges; + unsigned int num_verts = instance->network->num_verts; + + break_data *breaking_data = alloc_break_data(num_edges); + + while (true) { + double *voltages = get_voltage(instance, c); + double *field = get_current_v(instance, voltages, c); + + double conductivity = get_conductivity(instance, voltages, c); + if (conductivity < 1e-12 && instance->voltage_bound) { + free(voltages); + free(field); + break; + } + + int last_broke = inc_break_fuses(instance, fuse_thres, field, cutoff); + if (last_broke > num_edges || last_broke < -1) { + printf("%g \n", conductivity); + getchar(); + } + + update_break_data(breaking_data, last_broke, fabs(conductivity * fuse_thres[last_broke] / field[last_broke]), conductivity); + + free(voltages); + free(field); + + break_edge(instance, last_broke, c); + + if (instance->num_components > 1) { + break; + } + } + + return breaking_data; +} + diff --git a/src/free_network.c b/src/free_network.c new file mode 100644 index 0000000..651d3ae --- /dev/null +++ b/src/free_network.c @@ -0,0 +1,22 @@ + +#include "fracture.h" + +void free_fnet(fnet *network, cholmod_common *c) { + free(network->edges_to_verts); + if (network->edges_to_verts_break != network->edges_to_verts) { + free(network->edges_to_verts_break); + } + free(network->verts_to_edges_ind); + free(network->verts_to_edges); + free(network->bound_inds); + free(network->bound_verts); + free(network->vert_coords); + free(network->edge_coords); + free(network->dual_edges_to_verts); + free(network->dual_vert_coords); + free(network->dual_verts_to_edges_ind); + free(network->dual_verts_to_edges); + free(network->spanning_edges); + CHOL_F(free_sparse)(&(network->voltcurmat), c); + free(network); +} diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c new file mode 100644 index 0000000..6f4263a --- /dev/null +++ b/src/gen_laplacian.c @@ -0,0 +1,231 @@ + +#include "fracture.h" + +cholmod_sparse *gen_adjacency(const finst *instance, bool dual, bool breakv, + unsigned int pad, cholmod_common *c) { + unsigned int ne; + if (dual) + ne = ((int)instance->network->num_edges - + (int)instance->num_remaining_edges); + else { + ne = instance->num_remaining_edges; + if (!breakv && instance->network->boundary != TORUS_BOUND) { + ne += instance->network->bound_inds[2]; + } + } + + unsigned int nnz = 2 * ne; + + unsigned int nv; + if (dual) + nv = instance->network->num_dual_verts; + else { + if (breakv) nv = instance->network->num_verts_break; + else nv = instance->network->num_verts; + if (instance->network->boundary != TORUS_BOUND && !breakv) { + nv += 2; + } + } + + cholmod_triplet *t = + CHOL_F(allocate_triplet)(nv + pad, nv + pad, nnz, 0, CHOLMOD_REAL, c); + CHOL_INT *ri = (CHOL_INT *)t->i; + CHOL_INT *ci = (CHOL_INT *)t->j; + double *ai = (double *)t->x; + + t->nnz = nnz; + + unsigned int *etv; + if (dual) + etv = instance->network->dual_edges_to_verts; + else { + if (breakv) etv = instance->network->edges_to_verts_break; + else etv = instance->network->edges_to_verts; + } + + unsigned int count = 0; + + for (unsigned int i = 0; i < instance->network->num_edges; i++) { + if ((instance->fuses[i] && dual) || (!instance->fuses[i] && !dual)) { + unsigned int v1 = etv[2 * i]; + unsigned int v2 = etv[2 * i + 1]; + + ri[2 * count] = v1; + ci[2 * count] = v2; + ai[2 * count] = 1; + + ri[2 * count + 1] = v2; + ci[2 * count + 1] = v1; + ai[2 * count + 1] = 1; + + count++; + } + } + + if (!breakv && instance->network->boundary != TORUS_BOUND && !dual) { + for (unsigned int i = 0; i < 2; i++) { + for (unsigned int j = 0; j < instance->network->bound_inds[i+1] - instance->network->bound_inds[i]; j++) { + ri[2*count] = instance->network->num_verts + i; + ci[2*count] = instance->network->bound_verts[instance->network->bound_inds[i] + j]; + ai[2*count] = 1; + ri[2*count+1] = instance->network->bound_verts[instance->network->bound_inds[i] + j]; + ci[2*count+1] = instance->network->num_verts + i; + ai[2*count+1] = 1; + count++; + } + } + } + + + cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c); + CHOL_F(free_triplet)(&t, c); + + return s; +} + +cholmod_sparse *gen_laplacian(const finst *instance, cholmod_common *c, + bool symmetric) { + const fnet *network = instance->network; + unsigned int num_verts = network->num_verts_break; + double *vert_coords = network->vert_coords; + unsigned int num_bounds = network->num_bounds; + double inf = instance->inf; + bool voltage_bound = instance->voltage_bound; + unsigned int *bound_inds = network->bound_inds; + unsigned int *bound_verts = network->bound_verts; + + unsigned int num_gedges; + if (voltage_bound) { + num_gedges = 0; + } else { + num_gedges = network->bound_inds[num_bounds] + (num_bounds - 2) * 2; + } + if (network->boundary == TORUS_BOUND) num_gedges = bound_inds[1]; + unsigned int num_gverts = network->break_dim; + + unsigned int nnz = num_gverts + 2 * num_gedges; + + cholmod_triplet *temp_m = + CHOL_F(allocate_triplet)(num_gverts, num_gverts, nnz, 0, CHOLMOD_REAL, c); + + CHOL_INT *rowind = (CHOL_INT *)temp_m->i; + CHOL_INT *colind = (CHOL_INT *)temp_m->j; + double *acoo = (double *)temp_m->x; + + temp_m->nnz = nnz; + + for (unsigned int i = 0; i < num_gverts; i++) { + rowind[i] = i; + colind[i] = i; + acoo[i] = 0; + } + + cholmod_sparse *adjacency = gen_adjacency(instance, false, true, (int)num_gverts - (int)num_verts, c); + CHOL_INT *ia = (CHOL_INT *)adjacency->p; + CHOL_INT *ja = (CHOL_INT *)adjacency->i; + double *aa = (double *)adjacency->x; + + for (unsigned int i = 0; i < num_verts; i++) { + for (unsigned int j = 0; j < ia[i + 1] - ia[i]; j++) { + acoo[i] += aa[ia[i] + j]; + } + } + + if (network->boundary == TORUS_BOUND) { + for (unsigned int i = 0; i < network->bound_inds[1]; i++) { + unsigned int vv = network->bound_verts[i]; + rowind[num_gverts + 2*i] = vv; + colind[num_gverts + 2*i] = network->num_verts + i; + acoo[num_gverts + 2*i] = -1; + rowind[num_gverts + 2*i+1] = network->num_verts + i; + colind[num_gverts + 2*i+1] = vv; + acoo[num_gverts + 2*i+1] = -1; + acoo[vv]++; + acoo[network->num_verts + i]++; + } + acoo[bound_verts[0]]++; + } + + if (network->boundary != TORUS_BOUND) { + if (voltage_bound) { + for (unsigned int i = 0; i < 2; i++) { + for (unsigned int j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) { + acoo[bound_verts[bound_inds[i] + j]]++; + } + } + } else { + acoo[0]++; + for (unsigned int i = 0; i < num_bounds; i++) { + rowind[num_verts + i] = num_verts + i; + colind[num_verts + i] = num_verts + i; + if (i < 2) + acoo[num_verts + i] = (network->bound_inds[i + 1] - + network->bound_inds[i] + network->num_bounds - 2) * + inf; + else + acoo[num_verts + i] = + (network->bound_inds[i + 1] - network->bound_inds[i] + 2) * inf; + } + + unsigned int start = num_gverts; + + for (unsigned int i = 0; i < num_bounds; i++) { + for (unsigned int j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) { + rowind[start + bound_inds[i] + j] = bound_verts[bound_inds[i] + j]; + colind[start + bound_inds[i] + j] = num_verts + i; + acoo[start + bound_inds[i] + j] = -inf; + acoo[bound_verts[bound_inds[i] + j]] += inf; + colind[start + bound_inds[num_bounds] + bound_inds[i] + j] = + bound_verts[bound_inds[i] + j]; + rowind[start + bound_inds[num_bounds] + bound_inds[i] + j] = + num_verts + i; + acoo[start + bound_inds[num_bounds] + bound_inds[i] + j] = -inf; + } + } + + for (unsigned int i = 0; i < num_bounds - 2; i++) { + rowind[nnz - 2 * i - 1] = num_verts; + colind[nnz - 2 * i - 1] = num_verts + 2 + i; + acoo[nnz - 2 * i - 1] = -inf; + rowind[nnz - 2 * i - 2] = num_verts + 1; + colind[nnz - 2 * i - 2] = num_verts + 2 + i; + acoo[nnz - 2 * i - 2] = -inf; + colind[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = num_verts; + rowind[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = num_verts + 2 + i; + acoo[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = -inf; + colind[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = num_verts + 1; + rowind[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = num_verts + 2 + i; + acoo[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = -inf; + } + } + } + + + for (unsigned int i = 0; i < num_gverts; i++) { + if (acoo[i] == 0) + acoo[i] = 1; + } + + assert(CHOL_F(check_triplet)(temp_m, c)); + + cholmod_sparse *t_out = CHOL_F(triplet_to_sparse)(temp_m, nnz, c); + assert(CHOL_F(check_sparse)(t_out, c)); + + double alpha[2] = {1, 0}; + double beta[2] = {-1, 0}; + cholmod_sparse *t_laplacian = + CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c); + + cholmod_sparse *laplacian; + if (symmetric) + laplacian = CHOL_F(copy)(t_laplacian, 1, 1, c); + else + laplacian = CHOL_F(copy)(t_laplacian, 0, 1, c); + + CHOL_F(free_sparse)(&t_out, c); + CHOL_F(free_sparse)(&adjacency, c); + CHOL_F(free_sparse)(&t_laplacian, c); + CHOL_F(free_triplet)(&temp_m, c); + + return laplacian; +} diff --git a/src/gen_voltcurmat.c b/src/gen_voltcurmat.c new file mode 100644 index 0000000..ab16150 --- /dev/null +++ b/src/gen_voltcurmat.c @@ -0,0 +1,36 @@ + +#include "fracture.h" + +cholmod_sparse *gen_voltcurmat(unsigned int num_edges, unsigned int num_verts, + unsigned int *edges_to_verts, + cholmod_common *c) { + + cholmod_triplet *t_m = CHOL_F(allocate_triplet)( + num_edges, num_verts, num_edges * 2, 0, CHOLMOD_REAL, c); + assert(t_m != NULL); + + CHOL_INT *rowind = (CHOL_INT *)t_m->i; + CHOL_INT *colind = (CHOL_INT *)t_m->j; + double *acoo = (double *)t_m->x; + + for (unsigned int i = 0; i < num_edges; i++) { + unsigned int v1 = edges_to_verts[2 * i]; + unsigned int v2 = edges_to_verts[2 * i + 1]; + rowind[2 * i] = i; + rowind[2 * i + 1] = i; + colind[2 * i] = v1; + colind[2 * i + 1] = v2; + acoo[2 * i] = 1; + acoo[2 * i + 1] = -1; + } + + t_m->nnz = num_edges * 2; + + assert(CHOL_F(check_triplet)(t_m, c)); + + cholmod_sparse *m = CHOL_F(triplet_to_sparse)(t_m, num_edges * 2, c); + + CHOL_F(free_triplet)(&t_m, c); + + return m; +} diff --git a/src/geometry.c b/src/geometry.c new file mode 100644 index 0000000..83ef0bf --- /dev/null +++ b/src/geometry.c @@ -0,0 +1,55 @@ + +#include "fracture.h" + +int edge_to_verts(unsigned int width, bool periodic, unsigned int edge, + bool index) { + assert(edge < pow(width, 2)); + + int x = edge / width + 1; + int y = edge % width + 1; + + if (periodic) { + return (((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) - 1) % + width + + (x - index) * width) / + 2; + } else { + return ((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) + + (x - index) * (width + 1) - 1) / + 2; + } +} + +int dual_edge_to_verts(unsigned int width, bool periodic, unsigned int edge, + bool index) { + assert(edge < pow(width, 2)); + + int x = edge / width + 1; + int y = edge % width + 1; + + if (periodic) { + return (((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) - 1) % + width + + (x - index) * width) / + 2; + } else { + return ((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) + + (x - index) * (width + 1) - 1) / + 2; + } +} + +double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert, + bool index) { + if (periodic) { + if (index) + return vert % (width / 2) + (vert / (width / 2)) % 2; + else + return vert / (width / 2); + } else { + if (index) + return (2 * vert) % (width + 1); + else + return (2 * vert) / (width + 1); + } +} diff --git a/src/get_conductivity.c b/src/get_conductivity.c new file mode 100644 index 0000000..7ef08a8 --- /dev/null +++ b/src/get_conductivity.c @@ -0,0 +1,24 @@ + +#include "fracture.h" + +double get_conductivity(finst *inst, double *voltage, cholmod_common *c) { + if (inst->voltage_bound) { + double tot_cur = 0; + for (unsigned int i = 0; i < inst->network->num_spanning_edges; i++) { + unsigned int e = inst->network->spanning_edges[i]; + if (!inst->fuses[e]) { + unsigned int v1 = inst->network->edges_to_verts[2*e]; + unsigned int v2 = inst->network->edges_to_verts[2*e+1]; + double v1y = inst->network->vert_coords[2 * v1 + 1]; + double v2y = inst->network->vert_coords[2 * v2 + 1]; + unsigned int s1 = v1y < v2y ? v1 : v2; + unsigned int s2 = v1y < v2y ? v2 : v1; + tot_cur += voltage[s1] - voltage[s2]; + } + } + + return fabs(tot_cur); + } else { + return 1 / fabs(voltage[inst->network->num_verts] - voltage[inst->network->num_verts + 1]); + } +} diff --git a/src/get_current.c b/src/get_current.c new file mode 100644 index 0000000..5360054 --- /dev/null +++ b/src/get_current.c @@ -0,0 +1,99 @@ + +#include "fracture.h" + +double *get_voltage(const finst *instance, cholmod_common *c) { + cholmod_dense *b = instance->boundary_cond; + cholmod_factor *factor = instance->factor; + + cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); + + if (((double *)x->x)[0] != ((double *)x->x)[0]) { + printf("ERROR: GET_VOLTAGE FAILED\n\n"); + CHOL_F(free_dense)(&x, c); + return NULL; + } + + double *field = (double *)x->x; + x->x = NULL; + + CHOL_F(free_dense)(&x, c); + + + return field; +} + +double *get_current(const finst *instance, cholmod_common *c) { + unsigned int num_edges = instance->network->num_edges; + unsigned int num_verts = instance->network->num_verts_break; + unsigned int num_gverts = instance->network->break_dim; + unsigned int num_bounds = instance->network->num_bounds; + cholmod_sparse *voltcurmat = instance->network->voltcurmat; + + double *voltages = get_voltage(instance, c); + if (voltages == NULL) { + return NULL; + } + cholmod_dense *x = CHOL_F(allocate_dense)( + num_gverts, 1, num_gverts, CHOLMOD_REAL, c); + double *tmp_x = x->x; + x->x = voltages; + + cholmod_dense *y = + CHOL_F(allocate_dense)(num_edges, 1, num_edges, CHOLMOD_REAL, c); + + double alpha[2] = {1, 0}; + double beta[2] = {0, 0}; + CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); + + double *field = (double *)malloc(num_edges * sizeof(double)); + + for (int i = 0; i < num_edges; i++) { + if (instance->fuses[i]) + field[i] = 0; + else + field[i] = ((double *)y->x)[i]; + } + + x->x = tmp_x; + free(voltages); + CHOL_F(free_dense)(&x, c); + CHOL_F(free_dense)(&y, c); + + return field; +} + + +double *get_current_v(const finst *instance, double *voltages, cholmod_common *c) { + unsigned int num_edges = instance->network->num_edges; + unsigned int num_verts = instance->network->num_verts_break; + unsigned int num_gverts = instance->network->break_dim; + unsigned int num_bounds = instance->network->num_bounds; + cholmod_sparse *voltcurmat = instance->network->voltcurmat; + + cholmod_dense *x = CHOL_F(allocate_dense)( + num_gverts, 1, num_gverts, CHOLMOD_REAL, c); + double *tmp_x = x->x; + x->x = voltages; + + cholmod_dense *y = + CHOL_F(allocate_dense)(num_edges, 1, num_edges, CHOLMOD_REAL, c); + + double alpha[2] = {1, 0}; + double beta[2] = {0, 0}; + CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); + + double *field = (double *)malloc(num_edges * sizeof(double)); + + for (int i = 0; i < num_edges; i++) { + if (instance->fuses[i]) + field[i] = 0; + else + field[i] = ((double *)y->x)[i]; + } + + x->x = tmp_x; + CHOL_F(free_dense)(&x, c); + CHOL_F(free_dense)(&y, c); + + return field; +} diff --git a/src/get_dual_clusters.c b/src/get_dual_clusters.c new file mode 100644 index 0000000..488793f --- /dev/null +++ b/src/get_dual_clusters.c @@ -0,0 +1,36 @@ + +#include "fracture.h" + +unsigned int *get_clusters(finst *instance, cholmod_common *c) { + cholmod_sparse *s_dual = gen_adjacency(instance, true, false, 0, c); + + unsigned int *dual_marks = find_components(s_dual, 0); + CHOL_F(free_sparse)(&s_dual, c); + + return dual_marks; +} + +unsigned int *get_cluster_dist(finst *instance, cholmod_common *c) { + unsigned int *clusters = get_clusters(instance, c); + unsigned int *cluster_dist = (unsigned int *)calloc( + instance->network->num_dual_verts, sizeof(unsigned int)); + + unsigned int cur_mark = 0; + while (true) { + cur_mark++; + unsigned int num_in_cluster = 0; + for (unsigned int i = 0; i < instance->network->num_dual_verts; i++) { + if (clusters[i] == cur_mark) + num_in_cluster++; + } + + if (num_in_cluster == 0) + break; + + cluster_dist[num_in_cluster - 1]++; + } + + free(clusters); + + return cluster_dist; +} diff --git a/src/get_file.c b/src/get_file.c new file mode 100644 index 0000000..9141631 --- /dev/null +++ b/src/get_file.c @@ -0,0 +1,38 @@ + +#include "fracture.h" + +FILE *get_file(const char *prefix, unsigned int width, unsigned int crack, + double beta, unsigned int iter, unsigned int num_iter, + unsigned int num, bool read) { + int prefix_len = strlen(prefix); + int width_len = 1 + (int)log10(width); + int crack_len; + if (crack != 0) { + crack_len = 1 + (int)log10(crack); + } else { + crack_len = 1; + } + int beta_len; + if (beta > 1) { + beta_len = 1 + (int)log10(beta) + 3; + } else { + beta_len = 4; + } + int iter_len = 1 + (int)log10(num_iter); + int num_len = 1 + (int)log10(num); + int num_underscores = 5; + + int len = prefix_len + width_len + crack_len + beta_len + iter_len + num_len + + num_underscores + 4; + + char filename[len + 1]; + snprintf(filename, sizeof(filename), "%s_%u_%u_%.2f_%0*u_%u.txt", prefix, + width, crack, beta, iter_len, iter + 1, num); + + char *mode = "w"; + if (read) { + mode = "r"; + } + FILE *file = fopen(filename, mode); + return file; +} diff --git a/src/graph_components.c b/src/graph_components.c new file mode 100644 index 0000000..036f2fe --- /dev/null +++ b/src/graph_components.c @@ -0,0 +1,185 @@ + +#include "fracture.h" + +bool set_connected(const cholmod_sparse *laplacian, unsigned int *marks, + int vertex, int label, int stop_at, int exclude) { + unsigned int num_verts = (CHOL_INT)laplacian->nrow; + + CHOL_INT *ia = (CHOL_INT *)laplacian->p; + CHOL_INT *ja = (CHOL_INT *)laplacian->i; + double *a = (double *)laplacian->x; + + bool stopped = false; + + int *queue = (int *)malloc((num_verts) * sizeof(int)); + assert(queue != NULL); + + int queue_pos = 0; + + int *component_list = (int *)malloc((num_verts) * sizeof(int)); + assert(component_list != NULL); + + int list_pos = 0; + queue[0] = vertex; + component_list[0] = vertex; + + unsigned int old_label = marks[vertex]; + marks[vertex] = label; + + while (queue_pos >= 0) { + int next_vertex = queue[queue_pos]; + queue_pos--; + for (int j = 0; j < ia[next_vertex + 1] - ia[next_vertex]; j++) { + if (ja[ia[next_vertex] + j] == stop_at && a[ia[next_vertex] + j] != 0) { + // if we run into the other vertex, the components haven't changed + stopped = true; queue_pos = -1; + break; + } + if (marks[ja[ia[next_vertex] + j]] != label && + a[ia[next_vertex] + j] != 0 && + ja[ia[next_vertex] + j] < num_verts - exclude) { + marks[ja[ia[next_vertex] + j]] = label; + queue_pos++; + queue[queue_pos] = ja[ia[next_vertex] + j]; + list_pos++; + component_list[list_pos] = ja[ia[next_vertex] + j]; + } + } + } + + if (stopped) { + for (unsigned int i = 0; i <= list_pos; i++) { + marks[component_list[i]] = old_label; + } + } + + free(queue); + free(component_list); + + if (stopped) + return false; + else + return true; +} + +int update_components(const cholmod_sparse *laplacian, unsigned int *marks, + int old_num_components, int v1, int v2, int exclude) { + assert(laplacian != NULL); + assert(marks != NULL); + + if (set_connected(laplacian, marks, v1, old_num_components + 1, v2, + exclude)) { + return old_num_components + 1; + } else + return old_num_components; +} + +unsigned int *find_components(const cholmod_sparse *laplacian, unsigned int skip) { + size_t num_verts = laplacian->nrow; + unsigned int *marks = (unsigned int *)calloc(num_verts, sizeof(unsigned int)); + + int num_components = 0; + + for (int i = 0; i < num_verts; i++) { + if (marks[i] == 0) { + num_components++; + set_connected(laplacian, marks, i, num_components, -1, skip); + } + } + + return marks; +} + +unsigned int find_cycles(unsigned int num_edges, const bool *fuses, const unsigned int *ev, const unsigned + int *vei, const unsigned int *ve, int **cycles) { + unsigned int num_cycles = 0; + unsigned int *elist = (unsigned int *)malloc(num_edges * sizeof(unsigned int)); + unsigned int *side = (unsigned int *)calloc(num_edges, sizeof(unsigned int)); + unsigned int *num = (unsigned int *)malloc(num_edges * sizeof(unsigned int)); + for (unsigned int i = 0; i < num_edges; i++) { + if (fuses[i]) { + bool in_cycle = false; + for (unsigned int j = 0; j < num_cycles; j++) { + unsigned int k = 0; + int e; + while ((e = cycles[j][k]) >= 0) { + if (e == i) { + in_cycle = true; + break; + } + k++; + } + if (in_cycle) { + break; + } + } + if (!in_cycle) { + unsigned int pos = 0; + unsigned int cur_e = i; + bool *included = (bool *)calloc(num_edges, sizeof(bool)); + included[cur_e] = true; + elist[0] = cur_e; + side[0] = 1; + num[0] = 0; + while (true) { + unsigned int cv = ev[2 * cur_e + side[pos]]; + unsigned int new_e = cur_e; + unsigned int j; + for (j = num[pos]; j < vei[cv+1] - vei[cv]; j++) { + new_e = ve[vei[cv] + j]; + if (new_e != cur_e && fuses[new_e]) break; + } + if (new_e != cur_e && fuses[new_e]) { + if (new_e == elist[0] && cv == ev[2 * elist[0] + !side[0]]) { + pos++; + cycles[2*num_cycles] = (int *)malloc((pos + 1) * sizeof(int)); + cycles[2*num_cycles+1] = (int *)malloc((pos + 1) * sizeof(int)); + for (unsigned int i = 0; i < pos; i++) { + cycles[2*num_cycles][i] = elist[i]; + cycles[2*num_cycles+1][i] = side[i]; + } + cycles[2*num_cycles][pos] = -1; + cycles[2*num_cycles+1][pos] = -1; + num_cycles++; + pos--; + num[pos]++; + included[pos] = false; + cur_e = elist[pos]; + } else if (included[new_e]) { + if (pos == 0) break; + included[cur_e] = false; + pos--; + num[pos]++; + cur_e = elist[pos]; + } else { + num[pos] = j; + pos++; + elist[pos] = new_e; + included[new_e] = true; + unsigned int nv1 = ev[2 * new_e]; + unsigned int nv2 = ev[2 * new_e + 1]; + unsigned int v = nv1 == cv ? nv2 : nv1; + side[pos] = v == nv1 ? 0 : 1; + num[pos] = 0; + cur_e = new_e; + } + } else { + if (pos == 0) break; + included[cur_e] = false; + pos--; + num[pos]++; + cur_e = elist[pos]; + } + } + free(included); + } + } + } + + free(elist); + free(side); + free(num); + + return num_cycles; +} + diff --git a/src/homo_square_fracture.c b/src/homo_square_fracture.c new file mode 100644 index 0000000..a0dd305 --- /dev/null +++ b/src/homo_square_fracture.c @@ -0,0 +1,418 @@ + + +#include "fracture.h" + +int main(int argc, char *argv[]) { + int opt; + + // defining variables to be (potentially) set by command line flags + unsigned int N, L, filename_len; + double beta, inf, cutoff; + boundary_type boundary; + bool include_breaking, save_cluster_dist, use_voltage_boundaries, save_network, + save_crit_stress, save_toughness, save_corr, save_conductivity, + save_damage; + + filename_len = 100; + + N = 100; + L = 16; + beta = .3; + inf = 1e-8; + cutoff = 1e-10; + boundary = FREE_BOUND; + include_breaking = false; + save_cluster_dist = false; + use_voltage_boundaries = false; + save_network = false; + save_crit_stress = false; + save_damage = false; + save_corr = false; + save_conductivity = false; + save_toughness = false; + + int boundary_int; + char boundc2 = 'f'; + + while ((opt = getopt(argc, argv, "n:L:b:B:dVcoNsCrt")) != -1) { + switch (opt) { + case 'n': + N = atoi(optarg); + break; + case 'L': + L = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + case 'B': + boundary_int = atoi(optarg); + switch (boundary_int) { + case 0: + boundary = FREE_BOUND; + boundc2 = 'f'; + break; + case 1: + boundary = CYLINDER_BOUND; + boundc2 = 'c'; + break; + case 2: + boundary = TORUS_BOUND; + use_voltage_boundaries = true; + boundc2 = 't'; + break; + default: + printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n"); + } + break; + case 'd': + save_damage = true; + break; + case 'V': + use_voltage_boundaries = true; + break; + case 'c': + save_cluster_dist = true; + break; + case 'o': + include_breaking = true; + break; + case 'N': + save_network = true; + break; + case 's': + save_crit_stress = true; + break; + case 'C': + save_corr = true; + break; + case 'r': + save_conductivity = true; + inf = 1; + break; + case 't': + save_toughness = true; + break; + default: /* '?' */ + exit(EXIT_FAILURE); + } + } + + char boundc; + if (use_voltage_boundaries) boundc = 'v'; + else boundc = 'c'; + + FILE *break_out; + if (include_breaking) { + char *break_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(break_filename, filename_len, "breaks_s_%c_%c_%u_%g.txt", boundc, boundc2, L, beta); + break_out = fopen(break_filename, "a"); + free(break_filename); + } + + // start cholmod + cholmod_common c; + CHOL_F(start)(&c); + + /* if we use voltage boundary conditions, the laplacian matrix is positive + * definite and we can use a supernodal LL decomposition. otherwise we need + * to use the simplicial LDL decomposition + */ + if (use_voltage_boundaries) { + //(&c)->supernodal = CHOLMOD_SUPERNODAL; + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + } else { + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + } + + + fnet *network = ini_square_network(L, boundary, false, &c); + finst *perm_instance = create_instance(network, inf, use_voltage_boundaries, true, &c); + unsigned int c_dist_size = network->num_dual_verts; + unsigned int a_dist_size = network->num_verts; + + // define arrays for saving cluster and avalanche distributions + unsigned int *cluster_size_dist; + unsigned int *avalanche_size_dist; + char *c_filename; + if (save_cluster_dist) { + cluster_size_dist = + (unsigned int *)calloc(c_dist_size, sizeof(unsigned int)); + avalanche_size_dist = + (unsigned int *)calloc(a_dist_size, sizeof(unsigned int)); + + c_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(c_filename, filename_len, "cstr_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); + + FILE *cluster_out = fopen(c_filename, "r"); + + if (cluster_out != NULL) { + for (unsigned int i = 0; i < c_dist_size; i++) { + unsigned int tmp; + fscanf(cluster_out, "%u ", &tmp); + cluster_size_dist[i] = tmp; + } + fscanf(cluster_out, "\n"); + for (unsigned int i = 0; i < a_dist_size; i++) { + unsigned int tmp; + fscanf(cluster_out, "%u ", &tmp); + avalanche_size_dist[i] = tmp; + } + fclose(cluster_out); + } + } + + double *crit_stress; + if (save_crit_stress) { + crit_stress = (double *)malloc(N * sizeof(double)); + } + + double *avg_corr; + unsigned int **dists = NULL; + if (save_corr) { + avg_corr = (double *)calloc(c_dist_size, sizeof(double)); + } + + double *conductivity; + if (save_conductivity) { + conductivity = (double *)malloc(N * sizeof(double)); + } + + double *damage; + if (save_damage) { + damage = (double *)malloc(N * sizeof(double)); + } + + double *toughness; + if (save_toughness) { + toughness = (double *)malloc(N * sizeof(double)); + } + + + printf("\n"); + for (unsigned int i = 0; i < N; i++) { + printf("\033[F\033[JFRACTURE: %0*d / %d\n", (int)log10(N) + 1, i + 1, N); + + double *fuse_thres = gen_fuse_thres(network->num_edges, network->edge_coords, beta, beta_scaling_flat); + finst *instance = copy_instance(perm_instance, &c); + break_data *breaking_data = fracture_network(instance, fuse_thres, &c, cutoff); + free_instance(instance, &c); + free(fuse_thres); + + unsigned int max_pos = 0; + double max_val = 0; + + for (unsigned int j = 0; j < breaking_data->num_broken; j++) { + double val = breaking_data->extern_field[j]; + + if (val > max_val) { + max_pos = j; + max_val = val; + } + } + + if (save_crit_stress) { + crit_stress[i] = + breaking_data->extern_field[max_pos]; + } + + finst *tmp_instance = copy_instance(perm_instance, &c); + + unsigned int av_size = 0; + double cur_val = DBL_MAX; + for (unsigned int j = 0; j < max_pos; j++) { + break_edge(tmp_instance, breaking_data->break_list[j], &c); + + double val = fabs(breaking_data->extern_field[j]); + if (save_cluster_dist) { + if (val < cur_val) { + av_size++; + } + + if (val > cur_val) { + avalanche_size_dist[av_size]++; + av_size = 0; + cur_val = val; + } + } + } + + if (save_conductivity) { + if (!use_voltage_boundaries) { + double *tmp_voltage = get_voltage(tmp_instance, &c); + conductivity[i] = 1/fabs(tmp_voltage[tmp_instance->network->num_verts + 1] - tmp_voltage[tmp_instance->network->num_verts]); + free(tmp_voltage); + } else { + conductivity[i] = breaking_data->conductivity[max_pos]; + } + } + + if (save_toughness) { + double tmp_toughness = 0; + if (max_pos > 0) { + double sigma1 = breaking_data->extern_field[0]; + double epsilon1 = sigma1 / breaking_data->conductivity[0]; + for (unsigned int j = 0; j < max_pos - 1; j++) { + double sigma2 = breaking_data->extern_field[j+1]; + double epsilon2 = sigma2 / breaking_data->conductivity[j+1]; + if (epsilon2 > epsilon1) { + tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2; + sigma1 = sigma2; epsilon1 = epsilon2; + } + } + } + toughness[i] = tmp_toughness; + } + + if (save_damage) { + damage[i] = ((double)max_pos) / tmp_instance->network->num_edges; + } + + if (save_cluster_dist) { + unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c); + for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) { + cluster_size_dist[j] += tmp_cluster_dist[j]; + } + free(tmp_cluster_dist); + } + + if (save_corr) { + double *tmp_corr = get_corr(tmp_instance, dists, &c); + for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) { + avg_corr[i] += tmp_corr[j] / N; + } + free(tmp_corr); + } + + if (save_network) { + FILE *net_out = fopen("network.txt", "w"); + for (unsigned int j = 0; j < network->num_verts; j++) { + fprintf(net_out, "%f %f ", network->vert_coords[2 * j], + tmp_instance->network->vert_coords[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) { + fprintf(net_out, "%u %u ", tmp_instance->network->edges_to_verts[2 * j], + tmp_instance->network->edges_to_verts[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) { + fprintf(net_out, "%f %f ", tmp_instance->network->dual_vert_coords[2 * j], + tmp_instance->network->dual_vert_coords[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) { + fprintf(net_out, "%u %u ", tmp_instance->network->dual_edges_to_verts[2 * j], + tmp_instance->network->dual_edges_to_verts[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) { + fprintf(net_out, "%d ", tmp_instance->fuses[j]); + } + fclose(net_out); + } + + free_instance(tmp_instance, &c); + + if (include_breaking) { + for (unsigned int j = 0; j < breaking_data->num_broken; j++) { + fprintf(break_out, "%u %g %g ", breaking_data->break_list[j], + breaking_data->extern_field[j], breaking_data->conductivity[j]); + } + fprintf(break_out, "\n"); + } + + free_break_data(breaking_data); + } + + printf("\033[F\033[JFRACTURE: COMPLETE\n"); + + free_instance(perm_instance, &c); + free_fnet(network, &c); + + if (save_cluster_dist) { + FILE *cluster_out = fopen(c_filename, "w"); + + for (int i = 0; i < c_dist_size; i++) { + fprintf(cluster_out, "%u ", cluster_size_dist[i]); + } + fprintf(cluster_out, "\n"); + for (int i = 0; i < a_dist_size; i++) { + fprintf(cluster_out, "%u ", avalanche_size_dist[i]); + } + fclose(cluster_out); + + free(c_filename); + free(cluster_size_dist); + free(avalanche_size_dist); + } + + if (save_corr) { + char *corr_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(corr_filename, filename_len, "corr_s_%c_%c_%d_%g.txt", boundc, boundc2, L, + beta); + FILE *corr_file = fopen(corr_filename, "w"); + for (unsigned int i = 0; i < c_dist_size; i++) { + fprintf(corr_file, "%g ", avg_corr[i]); + } + fclose(corr_file); + free(corr_filename); + free(avg_corr); + } + + if (save_conductivity) { + char *cond_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(cond_filename, filename_len, "cond_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); + FILE *cond_file = fopen(cond_filename, "a"); + for (unsigned int i = 0; i < N; i++) { + fprintf(cond_file, "%g ", conductivity[i]); + } + fclose(cond_file); + free(cond_filename); + free(conductivity); + } + + if (save_toughness) { + char *tough_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(tough_filename, filename_len, "tuff_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); + FILE *tough_file = fopen(tough_filename, "a"); + for (unsigned int i = 0; i < N; i++) { + fprintf(tough_file, "%g ", toughness[i]); + } + fclose(tough_file); + free(tough_filename); + free(toughness); + } + + if (save_damage) { + char *hdam_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(hdam_filename, filename_len, "damg_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); + FILE *hdam_file = fopen(hdam_filename, "a"); + for (unsigned int i = 0; i < N; i++) { + fprintf(hdam_file, "%g ", damage[i]); + } + fclose(hdam_file); + free(hdam_filename); + free(damage); + } + + if (include_breaking) { + fclose(break_out); + } + + if (save_crit_stress) { + char *a_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(a_filename, filename_len, "strs_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); + FILE *a_file = fopen(a_filename, "a"); + for (int i = 0; i < N; i++) { + fprintf(a_file, "%g ", crit_stress[i]); + } + fclose(a_file); + free(a_filename); + free(crit_stress); + } + + CHOL_F(finish)(&c); + + return 0; +} diff --git a/src/homo_voronoi_fracture.c b/src/homo_voronoi_fracture.c new file mode 100644 index 0000000..1d9d0e6 --- /dev/null +++ b/src/homo_voronoi_fracture.c @@ -0,0 +1,420 @@ + + +#include "fracture.h" + +int main(int argc, char *argv[]) { + int opt; + + // defining variables to be (potentially) set by command line flags + unsigned int N, L, filename_len; + double beta, inf, cutoff; + boundary_type boundary; + bool include_breaking, save_cluster_dist, use_voltage_boundaries, save_network, + save_crit_stress, save_toughness, save_corr, save_conductivity, + save_damage; + + filename_len = 100; + + N = 100; + L = 16; + beta = .3; + inf = 1e10; + cutoff = 1e-9; + boundary = FREE_BOUND; + include_breaking = false; + save_cluster_dist = false; + use_voltage_boundaries = false; + save_network = false; + save_crit_stress = false; + save_damage = false; + save_corr = false; + save_conductivity = false; + save_toughness = false; + + int boundary_int; + char boundc2 = 'f'; + + while ((opt = getopt(argc, argv, "n:L:b:B:dVcoNsCrt")) != -1) { + switch (opt) { + case 'n': + N = atoi(optarg); + break; + case 'L': + L = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + case 'B': + boundary_int = atoi(optarg); + switch (boundary_int) { + case 0: + boundary = FREE_BOUND; + boundc2 = 'f'; + break; + case 1: + boundary = CYLINDER_BOUND; + boundc2 = 'c'; + break; + case 2: + boundary = TORUS_BOUND; + use_voltage_boundaries = true; + boundc2 = 't'; + break; + default: + printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n"); + } + break; + case 'd': + save_damage = true; + break; + case 'V': + use_voltage_boundaries = true; + break; + case 'c': + save_cluster_dist = true; + break; + case 'o': + include_breaking = true; + break; + case 'N': + save_network = true; + break; + case 's': + save_crit_stress = true; + break; + case 'C': + save_corr = true; + break; + case 'r': + save_conductivity = true; + //inf = 1; + break; + case 't': + save_toughness = true; + break; + default: /* '?' */ + exit(EXIT_FAILURE); + } + } + + char boundc; + if (use_voltage_boundaries) boundc = 'v'; + else boundc = 'c'; + + FILE *break_out; + if (include_breaking) { + char *break_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(break_filename, filename_len, "breaks_v_%c_%c_%u_%g.txt", boundc, boundc2, L, beta); + break_out = fopen(break_filename, "a"); + free(break_filename); + } + + unsigned int voronoi_max_verts = 4 * pow(L, 2); + unsigned int c_dist_size = voronoi_max_verts; + unsigned int a_dist_size = voronoi_max_verts; + + // define arrays for saving cluster and avalanche distributions + unsigned int *cluster_size_dist; + unsigned int *avalanche_size_dist; + char *c_filename; + if (save_cluster_dist) { + cluster_size_dist = + (unsigned int *)calloc(c_dist_size, sizeof(unsigned int)); + avalanche_size_dist = + (unsigned int *)calloc(a_dist_size, sizeof(unsigned int)); + + c_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(c_filename, filename_len, "cstr_v_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); + + FILE *cluster_out = fopen(c_filename, "r"); + + if (cluster_out != NULL) { + for (unsigned int i = 0; i < c_dist_size; i++) { + unsigned int tmp; + fscanf(cluster_out, "%u ", &tmp); + cluster_size_dist[i] = tmp; + } + fscanf(cluster_out, "\n"); + for (unsigned int i = 0; i < a_dist_size; i++) { + unsigned int tmp; + fscanf(cluster_out, "%u ", &tmp); + avalanche_size_dist[i] = tmp; + } + fclose(cluster_out); + } + } + + double *crit_stress; + if (save_crit_stress) { + crit_stress = (double *)malloc(N * sizeof(double)); + } + + double *avg_corr; + unsigned int **dists = NULL; + if (save_corr) { + avg_corr = (double *)calloc(voronoi_max_verts, sizeof(double)); + } + + double *conductivity; + if (save_conductivity) { + conductivity = (double *)malloc(N * sizeof(double)); + } + + double *damage; + if (save_damage) { + damage = (double *)malloc(N * sizeof(double)); + } + + double *toughness; + if (save_toughness) { + toughness = (double *)malloc(N * sizeof(double)); + } + + + // start cholmod + cholmod_common c; + CHOL_F(start)(&c); + + /* if we use voltage boundary conditions, the laplacian matrix is positive + * definite and we can use a supernodal LL decomposition. otherwise we need + * to use the simplicial LDL decomposition + */ + if (use_voltage_boundaries) { + //(&c)->supernodal = CHOLMOD_SUPERNODAL; + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + } else { + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + } + + + printf("\n"); + for (unsigned int i = 0; i < N; i++) { + printf("\033[F\033[JFRACTURE: %0*d / %d\n", (int)log10(N) + 1, i + 1, N); + + fnet *network = ini_voronoi_network(L, boundary, genfunc_hyperuniform, &c); + finst *perm_instance = create_instance(network, inf, use_voltage_boundaries, true, &c); + double *fuse_thres = gen_fuse_thres(network->num_edges, network->edge_coords, beta, beta_scaling_flat); + finst *instance = copy_instance(perm_instance, &c); + break_data *breaking_data = fracture_network(instance, fuse_thres, &c, cutoff); + bool debug_stop = instance->debug_stop; + free_instance(instance, &c); + free(fuse_thres); + + unsigned int max_pos = 0; + double max_val = 0; + + for (unsigned int j = 0; j < breaking_data->num_broken; j++) { + double val = breaking_data->extern_field[j]; + + if (val > max_val) { + max_pos = j; + max_val = val; + } + } + + if (save_crit_stress) { + crit_stress[i] = + breaking_data->extern_field[max_pos]; + } + + finst *tmp_instance = copy_instance(perm_instance, &c); + + unsigned int av_size = 0; + double cur_val = 0; + for (unsigned int j = 0; j < max_pos; j++) { + break_edge(tmp_instance, breaking_data->break_list[j], &c); + + double val = breaking_data->extern_field[j]; + if (save_cluster_dist) { + if (val < cur_val) { + av_size++; + } + + if (val > cur_val) { + avalanche_size_dist[av_size]++; + av_size = 0; + cur_val = val; + } + } + } + + if (save_conductivity) { + if (!use_voltage_boundaries) { + double *tmp_voltage = get_voltage(tmp_instance, &c); + conductivity[i] = 1/fabs(tmp_voltage[tmp_instance->network->num_verts + 1] - tmp_voltage[tmp_instance->network->num_verts]); + free(tmp_voltage); + } else { + conductivity[i] = breaking_data->conductivity[max_pos]; + } + } + + if (save_toughness) { + double tmp_toughness = 0; + if (max_pos > 0) { + double sigma1 = breaking_data->extern_field[0]; + double epsilon1 = sigma1 / breaking_data->conductivity[0]; + for (unsigned int j = 0; j < max_pos - 1; j++) { + double sigma2 = breaking_data->extern_field[j+1]; + double epsilon2 = sigma2 / breaking_data->conductivity[j+1]; + if (epsilon2 > epsilon1) { + tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2; + sigma1 = sigma2; epsilon1 = epsilon2; + } + } + } + toughness[i] = tmp_toughness; + } + + if (save_damage) { + damage[i] = ((double)max_pos) / tmp_instance->network->num_edges; + } + + if (save_cluster_dist) { + unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c); + for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) { + cluster_size_dist[j] += tmp_cluster_dist[j]; + } + free(tmp_cluster_dist); + } + + if (save_corr) { + double *tmp_corr = get_corr(tmp_instance, dists, &c); + for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) { + avg_corr[i] += tmp_corr[j] / N; + } + free(tmp_corr); + } + + if (save_network) { + FILE *net_out = fopen("network.txt", "w"); + for (unsigned int j = 0; j < network->num_verts; j++) { + fprintf(net_out, "%f %f ", network->vert_coords[2 * j], + tmp_instance->network->vert_coords[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) { + fprintf(net_out, "%u %u ", tmp_instance->network->edges_to_verts[2 * j], + tmp_instance->network->edges_to_verts[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) { + fprintf(net_out, "%f %f ", tmp_instance->network->dual_vert_coords[2 * j], + tmp_instance->network->dual_vert_coords[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) { + fprintf(net_out, "%u %u ", tmp_instance->network->dual_edges_to_verts[2 * j], + tmp_instance->network->dual_edges_to_verts[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) { + fprintf(net_out, "%d ", tmp_instance->fuses[j]); + } + fclose(net_out); + } + + free_instance(tmp_instance, &c); + free_instance(perm_instance, &c); + free_fnet(network, &c); + + if (include_breaking) { + for (unsigned int j = 0; j < breaking_data->num_broken; j++) { + fprintf(break_out, "%u %g %g ", breaking_data->break_list[j], + breaking_data->extern_field[j], breaking_data->conductivity[j]); + } + fprintf(break_out, "\n"); + } + + free_break_data(breaking_data); + if (debug_stop) break; + } + + printf("\033[F\033[JFRACTURE: COMPLETE\n"); + + if (save_cluster_dist) { + FILE *cluster_out = fopen(c_filename, "w"); + + for (int i = 0; i < c_dist_size; i++) { + fprintf(cluster_out, "%u ", cluster_size_dist[i]); + } + fprintf(cluster_out, "\n"); + for (int i = 0; i < a_dist_size; i++) { + fprintf(cluster_out, "%u ", avalanche_size_dist[i]); + } + fclose(cluster_out); + + free(c_filename); + free(cluster_size_dist); + free(avalanche_size_dist); + } + + if (save_corr) { + char *corr_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(corr_filename, filename_len, "corr_v_%c_%c_%d_%g.txt", boundc, boundc2, L, + beta); + FILE *corr_file = fopen(corr_filename, "w"); + for (unsigned int i = 0; i < voronoi_max_verts; i++) { + fprintf(corr_file, "%g ", avg_corr[i]); + } + fclose(corr_file); + free(corr_filename); + free(avg_corr); + } + + if (save_conductivity) { + char *cond_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(cond_filename, filename_len, "cond_v_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); + FILE *cond_file = fopen(cond_filename, "a"); + for (unsigned int i = 0; i < N; i++) { + fprintf(cond_file, "%g ", conductivity[i]); + } + fclose(cond_file); + free(cond_filename); + free(conductivity); + } + + if (save_toughness) { + char *tough_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(tough_filename, filename_len, "tuff_v_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); + FILE *tough_file = fopen(tough_filename, "a"); + for (unsigned int i = 0; i < N; i++) { + fprintf(tough_file, "%g ", toughness[i]); + } + fclose(tough_file); + free(tough_filename); + free(toughness); + } + + if (save_damage) { + char *hdam_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(hdam_filename, filename_len, "damg_v_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); + FILE *hdam_file = fopen(hdam_filename, "a"); + for (unsigned int i = 0; i < N; i++) { + fprintf(hdam_file, "%g ", damage[i]); + } + fclose(hdam_file); + free(hdam_filename); + free(damage); + } + + if (include_breaking) { + fclose(break_out); + } + + if (save_crit_stress) { + char *a_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(a_filename, filename_len, "strs_v_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); + FILE *a_file = fopen(a_filename, "a"); + for (int i = 0; i < N; i++) { + fprintf(a_file, "%g ", crit_stress[i]); + } + fclose(a_file); + free(a_filename); + free(crit_stress); + } + + CHOL_F(finish)(&c); + + return 0; +} diff --git a/src/ini_network.c b/src/ini_network.c new file mode 100644 index 0000000..0434d98 --- /dev/null +++ b/src/ini_network.c @@ -0,0 +1,620 @@ + +#include "fracture.h" + +unsigned int *get_spanning_edges(unsigned int num_edges, unsigned int *edges_to_verts, double *vert_coords, double cut, unsigned int *n) { + unsigned int *spanning_edges = (unsigned int *)malloc(num_edges * sizeof(unsigned int)); + (*n) = 0; + for (unsigned int i = 0; i < num_edges; i++) { + unsigned int v1, v2; + v1 = edges_to_verts[2 * i]; + v2 = edges_to_verts[2 * i + 1]; + double v1y, v2y; + v1y = vert_coords[2 * v1 + 1]; + v2y = vert_coords[2 * v2 + 1]; + if ((fabs(v1y - v2y) < 0.25) && ((v1y < cut && v2y > cut) || (v1y > cut && v2y < cut))) { + spanning_edges[*n] = i; + (*n)++; + } + } + return spanning_edges; +} + +double *get_edge_coords(unsigned int num_edges, double *vert_coords, + unsigned int *edges_to_verts) { + double *output = (double *)malloc(2 * num_edges * sizeof(double)); + +#pragma omp parallel for + for (unsigned int i = 0; i < num_edges; i++) { + unsigned int v1, v2; + double v1x, v1y, v2x, v2y, dx, dy; + v1 = edges_to_verts[2 * i]; + v2 = edges_to_verts[2 * i + 1]; + v1x = vert_coords[2 * v1]; + v1y = vert_coords[2 * v1 + 1]; + v2x = vert_coords[2 * v2]; + v2y = vert_coords[2 * v2 + 1]; + dx = v1x - v2x; + dy = v1y - v2y; + if (fabs(dx) > 0.5) { + if (dx > 0) { + v2x = v2x + 1; + } else { + v1x = v1x + 1; + } + } + if (fabs(dy) > 0.5) { + if (dy > 0) { + v2y = v2y + 1; + } else { + v1y = v1y + 1; + } + } + output[2 * i] = (v1x + v2x) / 2; + output[2 * i + 1] = (v1y + v2y) / 2; + } + + return output; +} + +unsigned int *get_verts_to_edges_ind(unsigned int num_verts, + unsigned int num_edges, + const unsigned int *edges_to_verts) { + unsigned int *output = + (unsigned int *)calloc(num_verts + 1, sizeof(unsigned int)); + assert(output != NULL); + + for (unsigned int i = 0; i < 2 * num_edges; i++) { + if (edges_to_verts[i] < num_verts) { + output[edges_to_verts[i] + 1]++; + } + } + + for (unsigned int i = 0; i < num_verts; i++) { + output[i + 1] += output[i]; + } + + return output; +} + +unsigned int *get_verts_to_edges(unsigned int num_verts, unsigned int num_edges, + const unsigned int *edges_to_verts, + const unsigned int *verts_to_edges_ind) { + unsigned int *output = (unsigned int *)calloc(verts_to_edges_ind[num_verts], + sizeof(unsigned int)); + unsigned int *counts = + (unsigned int *)calloc(num_verts, sizeof(unsigned int)); + for (int i = 0; i < 2 * num_edges; i++) { + if (edges_to_verts[i] < num_verts) { + output[verts_to_edges_ind[edges_to_verts[i]] + + counts[edges_to_verts[i]]] = i / 2; + counts[edges_to_verts[i]]++; + } + } + + free(counts); + + return output; +} + +fnet *ini_square_network(unsigned int width, boundary_type boundary, bool side_bounds, + cholmod_common *c) { + fnet *network = (fnet *)calloc(1, sizeof(fnet)); + + network->boundary = boundary; + bool periodic = (boundary == CYLINDER_BOUND) || (boundary == TORUS_BOUND) ? true : false; + + network->num_edges = pow(width, 2); + if (boundary == CYLINDER_BOUND) { + assert(width % 2 == 0); + assert(!side_bounds); + network->num_verts = (width / 2) * (width + 1); + network->num_verts_break = (width / 2) * (width + 1); + network->num_dual_verts = (width / 2 + 1) * (width / 2) + pow(width / 2, 2); + network->num_bounds = 2; + network->bound_inds = (unsigned int *)malloc((network->num_bounds + 1) * + sizeof(unsigned int)); + network->bound_inds[0] = 0; + network->bound_inds[1] = width / 2; + network->bound_inds[2] = width; + network->bound_verts = (unsigned int *)calloc(width, sizeof(unsigned int)); + network->break_dim = network->num_verts + network->num_bounds; + } else if (boundary == FREE_BOUND) { + network->num_verts = 2 * ((width + 1) / 2) * (width / 2 + 1); + network->num_verts_break = 2 * ((width + 1) / 2) * (width / 2 + 1); + network->num_dual_verts = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2); + if (side_bounds) + network->num_bounds = 4; + else + network->num_bounds = 2; + network->bound_inds = (unsigned int *)malloc((network->num_bounds + 1) * + sizeof(unsigned int)); + for (unsigned int i = 0; i < network->num_bounds + 1; i++) { + network->bound_inds[i] = i * ((width + 1) / 2); + } + network->bound_verts = (unsigned int *)calloc( + network->num_bounds * ((width + 1) / 2), sizeof(unsigned int)); + if (side_bounds) { + for (unsigned int i = 0; i < (width + 1) / 2; i++) { + network->bound_verts[2 * ((width + 1) / 2) + i] = + (width + 1) / 2 + i * (width + 1); + if (width % 2) { + network->bound_verts[3 * ((width + 1) / 2) + i] = + (width + 1) / 2 + i * (width + 1) - 1; + } else { + network->bound_verts[3 * ((width + 1) / 2) + i] = + (width + 1) / 2 + i * (width + 1) + width / 2; + } + } + } + network->break_dim = network->num_verts + network->num_bounds; + } else if (boundary == TORUS_BOUND) { + network->num_verts = (width / 2) * (width + 1) - (width / 2); + network->num_verts_break = (width / 2) * (width + 1); + network->num_dual_verts = (width / 2 + 1) * (width / 2) + pow(width / 2, 2) - (width / 2); + network->num_bounds = 1; + network->bound_inds = (unsigned int *)malloc((network->num_bounds + 1) * + sizeof(unsigned int)); + network->bound_inds[0] = 0; + network->bound_inds[1] = width / 2; + network->bound_verts = (unsigned int *)calloc(width / 2, sizeof(unsigned int)); + network->break_dim = network->num_verts_break; + } + if (boundary != TORUS_BOUND) { + for (unsigned int i = 0; i < (width + 1) / 2; i++) { + network->bound_verts[i] = i; + network->bound_verts[(width + 1) / 2 + i] = network->num_verts - 1 - i; + } + } else { + for (unsigned int i = 0; i < width / 2; i++) { + network->bound_verts[i] = i; + } + } + network->edges_to_verts_break = + (unsigned int *)calloc(2 * network->num_edges, sizeof(unsigned int)); + network->edges_to_verts = + (unsigned int *)calloc(2 * network->num_edges, sizeof(unsigned int)); + for (unsigned int i = 0; i < network->num_edges; i++) { + network->edges_to_verts_break[2 * i] = edge_to_verts(width, periodic, i, 1); + network->edges_to_verts_break[2 * i + 1] = edge_to_verts(width, periodic, i, 0); + network->edges_to_verts[2 * i] = network->edges_to_verts_break[2 * i] % network->num_verts; + network->edges_to_verts[2 * i + 1] = network->edges_to_verts_break[2 * i + 1] % network->num_verts; + } + network->verts_to_edges_ind = + (unsigned int *)calloc(network->num_verts + 1, sizeof(unsigned int)); + network->verts_to_edges_ind[0] = 0; + unsigned int pos1 = 0; + for (unsigned int i = 0; i < network->num_verts; i++) { + bool in_bound = false; + for (unsigned int j = 0; j < network->num_bounds; j++) { + for (unsigned int k = 0; + k < network->bound_inds[j + 1] - network->bound_inds[j]; k++) { + if (i == network->bound_verts[network->bound_inds[j] + k]) { + in_bound = true; + break; + } + } + } + if (in_bound) + pos1 += 2; + else + pos1 += 4; + + network->verts_to_edges_ind[i + 1] = pos1; + } + network->verts_to_edges = (unsigned int *)calloc( + network->verts_to_edges_ind[network->num_verts], sizeof(unsigned int)); + unsigned int *vert_counts = + (unsigned int *)calloc(network->num_verts, sizeof(unsigned int)); + for (unsigned int i = 0; i < network->num_edges; i++) { + unsigned int v0 = network->edges_to_verts[2 * i]; + unsigned int v1 = network->edges_to_verts[2 * i + 1]; + network->verts_to_edges[network->verts_to_edges_ind[v0] + vert_counts[v0]] = + i; + network->verts_to_edges[network->verts_to_edges_ind[v1] + vert_counts[v1]] = + i; + vert_counts[v0]++; + vert_counts[v1]++; + } + free(vert_counts); + + network->vert_coords = + (double *)malloc(2 * network->num_verts * sizeof(double)); + for (unsigned int i = 0; i < network->num_verts; i++) { + if (!periodic) { + network->vert_coords[2 * i] = ((double)((2 * i + 1) / (width + 1)))/width; + network->vert_coords[2 * i + 1] = ((double)((2 * i + 1) % (width + 1)))/width; + } else { + network->vert_coords[2 * i] = ((double)((2 * i + 1) / (width)))/width; + network->vert_coords[2 * i + 1] = ((double)((((2 * i + 1) / (width)+1) % 2) + ((2 * i) % (width))))/width; + } + } + + network->L = width; + network->edge_coords = get_edge_coords( + network->num_edges, network->vert_coords, network->edges_to_verts); + + network->dual_edges_to_verts = + (unsigned int *)malloc(2 * network->num_edges * sizeof(unsigned int)); + for (unsigned int i = 0; i < network->num_edges; i++) { + network->dual_edges_to_verts[2 * i] = + dual_edge_to_verts(width, periodic, i, 0) % network->num_dual_verts; + network->dual_edges_to_verts[2 * i + 1] = + dual_edge_to_verts(width, periodic, i, 1) % network->num_dual_verts; + } + network->dual_vert_coords = + (double *)malloc(2 * network->num_dual_verts * sizeof(double)); + for (unsigned int i = 0; i < network->num_dual_verts; i++) { + network->dual_vert_coords[2 * i] = + 2*dual_vert_to_coord(width, periodic, i, 0); + network->dual_vert_coords[2 * i + 1] = + 2*dual_vert_to_coord(width, periodic, i, 1); + } + + network->voltcurmat = gen_voltcurmat(network->num_edges, + network->break_dim, + network->edges_to_verts_break, c); + + network->dual_verts_to_edges_ind = + get_verts_to_edges_ind(network->num_dual_verts, network->num_edges, + network->dual_edges_to_verts); + network->dual_verts_to_edges = get_verts_to_edges( + network->num_dual_verts, network->num_edges, network->dual_edges_to_verts, + network->dual_verts_to_edges_ind); + network->spanning_edges = get_spanning_edges(network->num_edges, network->edges_to_verts_break, network->vert_coords, 0.51, &(network->num_spanning_edges)); + + return network; +} + +unsigned int *get_voro_dual_edges(unsigned int num_edges, + unsigned int num_verts, unsigned int *edges, + unsigned int *triangles) { + unsigned int *dual_edges = + (unsigned int *)malloc(2 * num_edges * sizeof(unsigned int)); +#pragma omp parallel for + for (unsigned int i = 0; i < num_edges; i++) { + unsigned int v1, v2; + v1 = edges[2 * i]; + v2 = edges[2 * i + 1]; + if (v1 < num_verts && v2 < num_verts) { + bool found_match = false; + for (unsigned int j = 0; j < 3; j++) { + for (unsigned int k = 0; k < 3; k++) { + unsigned int t11, t12, t21, t22; + t11 = triangles[3 * v1 + j]; + t12 = triangles[3 * v1 + ((j + 1) % 3)]; + t21 = triangles[3 * v2 + k]; + t22 = triangles[3 * v2 + ((k + 1) % 3)]; + if ((t11 == t21 && t12 == t22) || (t11 == t22 && t12 == t21)) { + dual_edges[2 * i] = t11 < t12 ? t11 : t12; + dual_edges[2 * i + 1] = t11 < t12 ? t12 : t11; + found_match = true; + break; + } + } + if (found_match) + break; + } + } + } + + return dual_edges; +} + +fnet *ini_voronoi_network(unsigned int L, boundary_type boundary, + double *(*genfunc)(unsigned int, gsl_rng *, unsigned int *), + cholmod_common *c) { + fnet *network = (fnet *)calloc(1, sizeof(fnet)); + + // generate the dual lattice + double *lattice; + unsigned int num; + { + gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); + FILE *rf = fopen("/dev/urandom", "r"); + unsigned long int seed; + fread(&seed, sizeof(unsigned long int), 1, rf); + fclose(rf); + gsl_rng_set(r, seed); + lattice = genfunc(L, r, &num); + gsl_rng_free(r); + } + + // retrieve a periodic voronoi tesselation of the lattice + intptr_t *vout = run_voronoi(num, lattice, 0, 1, 0, 1); + + unsigned int tmp_num_verts = ((unsigned int *)vout[0])[0]; + unsigned int tmp_num_edges = ((unsigned int *)vout[0])[1]; + double *tmp_vert_coords = (double *)vout[2]; + unsigned int *tmp_edges = (unsigned int *)vout[3]; + unsigned int *tmp_tris = (unsigned int *)vout[5]; + + free((void *)vout[0]); + free((void *)vout[1]); + free((void *)vout[4]); + free(vout); + + // get dual edges of the fully periodic graph + unsigned int *tmp_dual_edges = + get_voro_dual_edges(tmp_num_edges, tmp_num_verts, tmp_edges, tmp_tris); + + // prune the edges of the lattice and assign boundary vertices based on the + // desired boundary conditions + unsigned int num_bounds; + unsigned num_verts; + double *vert_coords; + unsigned int *bound_inds; + unsigned int *bound_verts; + unsigned int num_edges; + unsigned int *edges; + unsigned int *dual_edges; + switch (boundary) { + case FREE_BOUND: { + num_bounds = 4; + bound_inds = (unsigned int *)malloc((1 + num_bounds) * sizeof(unsigned int)); + bound_inds[0] = 0; + vert_coords = tmp_vert_coords; + num_verts = tmp_num_verts; + num_edges = 0; + edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int)); + dual_edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int)); + unsigned int num_t, num_b, num_l, num_r; + bool *bound_t, *bound_b, *bound_l, *bound_r; + num_t = 0; num_b = 0; num_l = 0; num_r = 0; + bound_t = (bool *)calloc(num_verts, sizeof(bool)); + bound_b = (bool *)calloc(num_verts, sizeof(bool)); + bound_l = (bool *)calloc(num_verts, sizeof(bool)); + bound_r = (bool *)calloc(num_verts, sizeof(bool)); + for (unsigned int i = 0; i < tmp_num_edges; i++) { + unsigned int v1, v2; + double v1x, v1y, v2x, v2y, dx, dy; + v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1]; + v1x = vert_coords[2 * v1]; v1y = vert_coords[2 * v1 + 1]; + v2x = vert_coords[2 * v2]; v2y = vert_coords[2 * v2 + 1]; + dx = v1x - v2x; dy = v1y - v2y; + if (fabs(dy) > 0.5) { + if (dy > 0) { + if (!bound_t[v1] && !bound_l[v1] && !bound_r[v1]) { + bound_t[v1] = true; num_t++; + } + if (!bound_b[v2] && !bound_l[v2] && !bound_r[v2]) { + bound_b[v2] = true; num_b++; + } + } else { + if (!bound_t[v2] && !bound_l[v2] && !bound_r[v2]) { + bound_t[v2] = true; num_t++; + } + if (!bound_b[v1] && !bound_l[v1] && !bound_r[v1]) { + bound_b[v1] = true; num_b++; + } + } + } else if (fabs(dx) > 0.5) { + if (dx > 0) { + if (!bound_r[v1] && !bound_t[v1] && !bound_b[v1]) { + bound_r[v1] = true; num_r++; + } + if (!bound_l[v2] && !bound_t[v2] && !bound_b[v2]) { + bound_l[v2] = true; num_l++; + } + } else { + if (!bound_r[v2] && !bound_t[v2] && !bound_b[v2]) { + bound_r[v2] = true; num_r++; + } + if (!bound_l[v1] && !bound_t[v1] && !bound_b[v1]) { + bound_l[v1] = true; num_l++; + } + } + } else { + edges[2 * num_edges] = v1 < v2 ? v1 : v2; + edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1; + unsigned int d1 = tmp_dual_edges[2 * i]; + unsigned int d2 = tmp_dual_edges[2 * i + 1]; + dual_edges[2 * num_edges] = d1 < d2 ? d1 : d2; + dual_edges[2 * num_edges + 1] = d1 < d2 ? d2 : d1; + num_edges++; + } + } + bound_verts = malloc((num_t + num_b + num_l + num_r) * sizeof(unsigned int)); + bound_inds[1] = num_t; bound_inds[2] = num_t + num_b; + bound_inds[3] = num_l + num_t + num_b; + bound_inds[4] = num_t + num_b + num_r + num_l; + unsigned int pos_t, pos_b, pos_l, pos_r; + pos_t = 0; pos_b = 0; pos_l = 0; pos_r = 0; + for (unsigned int i = 0; i < num_verts; i++) { + if (bound_t[i]) { + bound_verts[pos_t] = i; pos_t++; + } else if (bound_b[i]) { + bound_verts[num_t + pos_b] = i; pos_b++; + } else if (bound_l[i]) { + bound_verts[num_t + num_b + pos_l] = i; pos_l++; + } else if (bound_r[i]) { + bound_verts[num_t + num_b + num_l + pos_r] = i; pos_r++; + } + } + free(bound_l); free(bound_r); free(bound_t); free(bound_b); + free(tmp_edges); + free(tmp_dual_edges); + network->edges_to_verts_break = edges; + network->edges_to_verts = edges; + network->num_verts_break = num_verts; + network->num_verts = num_verts; + network->break_dim = num_verts + num_bounds; + break; + } + case CYLINDER_BOUND: { + num_bounds = 2; + bound_inds = (unsigned int *)malloc((1 + num_bounds) * sizeof(unsigned int)); + bound_inds[0] = 0; + vert_coords = tmp_vert_coords; + num_verts = tmp_num_verts; + num_edges = 0; + edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int)); + dual_edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int)); + unsigned int num_t, num_b; + bool *bound_t, *bound_b; + num_t = 0; num_b = 0; + bound_t = (bool *)calloc(num_verts, sizeof(bool)); + bound_b = (bool *)calloc(num_verts, sizeof(bool)); + for (unsigned int i = 0; i < tmp_num_edges; i++) { + unsigned int v1, v2; + double v1x, v1y, v2x, v2y, dx, dy; + v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1]; + v1y = vert_coords[2 * v1 + 1]; v2y = vert_coords[2 * v2 + 1]; + dy = v1y - v2y; + if (fabs(dy) > 0.5) { + if (dy > 0) { + if (!bound_t[v1]) { + bound_t[v1] = true; num_t++; + } + if (!bound_b[v2]) { + bound_b[v2] = true; num_b++; + } + } else { + if (!bound_t[v2]) { + bound_t[v2] = true; num_t++; + } + if (!bound_b[v1]) { + bound_b[v1] = true; num_b++; + } + } + } else { + edges[2 * num_edges] = v1 < v2 ? v1 : v2; + edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1; + unsigned int d1 = tmp_dual_edges[2 * i]; + unsigned int d2 = tmp_dual_edges[2 * i + 1]; + dual_edges[2 * num_edges] = d1 < d2 ? d1 : d2; + dual_edges[2 * num_edges + 1] = d1 < d2 ? d2 : d1; + num_edges++; + } + } + bound_verts = malloc((num_t + num_b) * sizeof(unsigned int)); + bound_inds[1] = num_t; bound_inds[2] = num_t + num_b; + unsigned int pos_t, pos_b; + pos_t = 0; pos_b = 0; + for (unsigned int i = 0; i < num_verts; i++) { + if (bound_t[i]) { + bound_verts[pos_t] = i; pos_t++; + } else if (bound_b[i]) { + bound_verts[num_t + pos_b] = i; pos_b++; + } + } + free(bound_t); free(bound_b); + free(tmp_edges); + free(tmp_dual_edges); + network->edges_to_verts_break = edges; + network->edges_to_verts = edges; + network->num_verts_break = num_verts; + network->num_verts = num_verts; + network->break_dim = num_verts + num_bounds; + break; + } + case TORUS_BOUND: { + num_bounds = 1; + bound_inds = (unsigned int *)malloc((1 + num_bounds) * sizeof(unsigned int)); + bound_inds[0] = 0; + num_edges = tmp_num_edges; + edges = (unsigned int *)malloc(2* num_edges*sizeof(unsigned int)); + for (unsigned int i = 0; i < num_edges; i++) { + edges[2*i] = tmp_edges[2*i]; + edges[2*i+1] = tmp_edges[2*i+1]; + } + dual_edges = tmp_dual_edges; + bool *bound_t = (bool *)calloc(tmp_num_verts, sizeof(bool)); + int *edge_change = (int *)calloc(num_edges, sizeof(int)); + unsigned int num_t = 0; + for (unsigned int i = 0; i < num_edges; i++) { + unsigned int v1, v2; + double v1x, v1y, v2x, v2y, dx, dy; + v1 = edges[2 * i]; v2 = edges[2 * i + 1]; + v1x = tmp_vert_coords[2 * v1]; v1y = tmp_vert_coords[2 * v1 + 1]; + v2x = tmp_vert_coords[2 * v2]; v2y = tmp_vert_coords[2 * v2 + 1]; + dy = v1y - v2y; + if (fabs(dy) > 0.5) { + if (dy > 0) { + edge_change[i] = 1; + if (!bound_t[v1]) { + bound_t[v1] = true; num_t++; + } + } else { + edge_change[i] = 2; + if (!bound_t[v2]) { + bound_t[v2] = true; num_t++; + } + } + } + } + num_verts = tmp_num_verts + num_t; + vert_coords = (double *)malloc(2 * num_verts * sizeof(double)); + bound_verts = malloc(num_t * sizeof(unsigned int)); + bound_inds[1] = num_t; + unsigned int pos_t = 0; + for (unsigned int i = 0; i < tmp_num_verts; i++) { + vert_coords[2*i] = tmp_vert_coords[2*i]; + vert_coords[2*i+1] = tmp_vert_coords[2*i+1]; + if (bound_t[i]) { + bound_verts[pos_t] = i; + vert_coords[2*(tmp_num_verts + pos_t)] = tmp_vert_coords[2*i]; + vert_coords[2*(tmp_num_verts + pos_t)+1] = tmp_vert_coords[2*i+1]; + pos_t++; + } + } + for (unsigned int i = 0; i < num_edges; i++) { + if (edge_change[i]) { + for (unsigned int j = 0; j < num_t; j++) { + if (edges[2*i+(edge_change[i]-1)] == bound_verts[j]) { + edges[2*i+(edge_change[i]-1)] = tmp_num_verts + j; + break; + } + } + } + } + free(tmp_vert_coords); + free(bound_t); + free(edge_change); + network->num_verts_break = num_verts; + network->num_verts = tmp_num_verts; + network->edges_to_verts_break = edges; + network->edges_to_verts = tmp_edges; + network->break_dim = num_verts; + break; + } + } + + network->boundary = boundary; + network->num_bounds = num_bounds; + network->bound_inds = bound_inds; + network->bound_verts = bound_verts; + network->num_edges = num_edges; + network->dual_edges_to_verts = dual_edges; + network->vert_coords = vert_coords; + + network->verts_to_edges_ind = get_verts_to_edges_ind( + network->num_verts, network->num_edges, network->edges_to_verts); + network->verts_to_edges = + get_verts_to_edges(network->num_verts, network->num_edges, + network->edges_to_verts, network->verts_to_edges_ind); + + network->L = 1; + + network->edge_coords = get_edge_coords( + network->num_edges, network->vert_coords, network->edges_to_verts); + + free(tmp_tris); + + network->dual_vert_coords = lattice; + network->num_dual_verts = num; + + network->voltcurmat = gen_voltcurmat(network->num_edges, + network->break_dim, + network->edges_to_verts_break, c); + + network->dual_verts_to_edges_ind = + get_verts_to_edges_ind(network->num_dual_verts, network->num_edges, + network->dual_edges_to_verts); + network->dual_verts_to_edges = get_verts_to_edges( + network->num_dual_verts, network->num_edges, network->dual_edges_to_verts, + network->dual_verts_to_edges_ind); + + network->spanning_edges = get_spanning_edges(network->num_edges, network->edges_to_verts_break, network->vert_coords, 0.5, &(network->num_spanning_edges)); + + return network; +} diff --git a/src/instance.c b/src/instance.c new file mode 100644 index 0000000..5b3d8fb --- /dev/null +++ b/src/instance.c @@ -0,0 +1,118 @@ + +#include "fracture.h" + +finst *create_instance(fnet *network, double inf, bool voltage_bound, + bool startnow, cholmod_common *c) { + finst *instance = (finst *)calloc(1, sizeof(finst)); + assert(instance != NULL); + + instance->network = network; + instance->num_remaining_edges = network->num_edges; + instance->fuses = (bool *)calloc(network->num_edges, sizeof(bool)); + assert(instance->fuses != NULL); + instance->inf = inf; + instance->voltage_bound = voltage_bound; + instance->boundary_cond = CHOL_F(zeros)( + network->break_dim, 1, CHOLMOD_REAL, c); + if (network->boundary != TORUS_BOUND) { + if (voltage_bound) { + for (unsigned int i = 0; i < network->bound_inds[1]; i++) { + ((double *)instance->boundary_cond->x)[network->bound_verts[i]] =1; + } + } else { + ((double *)instance->boundary_cond->x)[0] = 1; + ((double *)instance->boundary_cond->x)[network->num_verts] = 1; + ((double *)instance->boundary_cond->x)[network->num_verts + 1] = -1; + } + } else { + for (unsigned int i = 0; i < network->bound_inds[1]; i++) { + ((double *)instance->boundary_cond->x)[network->bound_verts[i]] = 1; + ((double *)instance->boundary_cond->x)[network->num_verts + i] = -1; + } + ((double *)instance->boundary_cond->x)[network->bound_verts[0]] = 1; + } + + instance->adjacency = gen_adjacency(instance, false, false, 0, c); + + if (startnow) { + cholmod_sparse *laplacian = gen_laplacian(instance, c, true); + instance->factor = CHOL_F(analyze)(laplacian, c); + CHOL_F(factorize)(laplacian, instance->factor, c); + CHOL_F(free_sparse)(&laplacian, c); + } + + instance->marks = (unsigned int *)malloc( + (instance->network->break_dim) * + sizeof(unsigned int)); + instance->dual_marks = (unsigned int *)malloc( + (instance->network->num_dual_verts) * + sizeof(unsigned int)); + assert(instance->marks != NULL); + + for (unsigned int i = 0; + i < (instance->network->break_dim); + i++) { + instance->marks[i] = 1; + } + for (unsigned int i = 0; + i < (instance->network->num_dual_verts); + i++) { + instance->dual_marks[i] = i+1; + } + instance->num_components = 1; + + return instance; +} + +void finish_instance(finst *instance, cholmod_common *c) { + cholmod_sparse *laplacian = gen_laplacian(instance, c, true); + instance->factor = CHOL_F(analyze)(laplacian, c); + CHOL_F(factorize)(laplacian, instance->factor, c); + CHOL_F(free_sparse)(&laplacian, c); +} + +finst *copy_instance(const finst *instance, cholmod_common *c) { + finst *instance_copy = (finst *)calloc(1, sizeof(finst)); + memcpy(instance_copy, instance, sizeof(finst)); + + size_t fuses_size = (instance->network)->num_edges * sizeof(bool); + instance_copy->fuses = (bool *)malloc(fuses_size); + memcpy(instance_copy->fuses, instance->fuses, fuses_size); + + size_t marks_size = + (instance->network->break_dim) * + sizeof(unsigned int); + instance_copy->marks = (unsigned int *)malloc(marks_size); + memcpy(instance_copy->marks, instance->marks, marks_size); + + size_t dual_marks_size = + (instance->network->num_dual_verts) * + sizeof(unsigned int); + instance_copy->dual_marks = (unsigned int *)malloc(dual_marks_size); + memcpy(instance_copy->dual_marks, instance->dual_marks, dual_marks_size); + + instance_copy->adjacency = CHOL_F(copy_sparse)(instance->adjacency, c); + instance_copy->boundary_cond = CHOL_F(copy_dense)(instance->boundary_cond, c); + instance_copy->factor = CHOL_F(copy_factor)(instance->factor, c); + + return instance_copy; +} + +void free_instance(finst *instance, cholmod_common *c) { + free(instance->fuses); + CHOL_F(free_dense)(&(instance->boundary_cond), c); + CHOL_F(free_sparse)(&(instance->adjacency), c); + CHOL_F(free_factor)(&(instance->factor), c); + free(instance->marks); + free(instance->dual_marks); + free(instance); +} + +bool check_instance(const finst *instance, cholmod_common *c) { + assert(instance != NULL); + assert(instance->fuses != NULL); + assert(CHOL_F(check_dense)(instance->boundary_cond, c)); + assert(CHOL_F(check_factor)(instance->factor, c)); + + return true; +} diff --git a/src/randfuncs.c b/src/randfuncs.c new file mode 100644 index 0000000..260c047 --- /dev/null +++ b/src/randfuncs.c @@ -0,0 +1,71 @@ + +#include "fracture.h" + +double *genfunc_uniform(unsigned int L, gsl_rng *r, unsigned int *num) { + *num = pow(L / 2 + 1, 2) + pow((L + 1) / 2, 2); + + double *lattice = (double *)malloc(2 * (*num) * sizeof(double)); + for (unsigned int i = 0; i < (*num); i++) { + lattice[2*i] = gsl_ran_flat(r, 0, 1); + lattice[2*i+1] = gsl_ran_flat(r, 0, 1); + } + + return lattice; +} + +double *genfunc_hyperuniform(unsigned int L, gsl_rng *r, unsigned int *num) { + *num = pow(L / 2 + 1, 2) + pow((L + 1) / 2, 2); + + // necessary to prevent crashing when underflow occurs + gsl_set_error_handler_off(); + + double *lattice = (double *)malloc(2 * (*num) * sizeof(double)); + double rho = *num; + for (unsigned int i = 0; i < (*num); i++) { + bool reject = true; + double x, y; + while(reject) { + x = gsl_ran_flat(r, 0, 1); + y = gsl_ran_flat(r, 0, 1); + reject = false; + for (unsigned int j = 0; j < i; j++) { + double *ds = (double *)malloc(5 * sizeof(double)); + ds[0] = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1],2); + ds[1] = pow(x-lattice[2*j] + 1,2)+pow(y-lattice[2*j+1],2); + ds[2] = pow(x-lattice[2*j] - 1,2)+pow(y-lattice[2*j+1],2); + ds[3] = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1] + 1,2); + ds[4] = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1] - 1,2); + unsigned int min_pos = 0; double min_val = 100; + for (unsigned int k = 0; k < 5; k++) { + if (min_val > ds[k]) { + min_pos = k; min_val = ds[k]; + } + } + if (1-gsl_sf_exp(-M_PI * rho * min_val) < gsl_ran_flat(r, 0, 1)) { + reject = true; + break; + } + } + } + lattice[2 * i] = x; + lattice[2 * i + 1] = y; + } + + return lattice; +} + +void randfunc_flat(gsl_rng *r, double *x, double *y) { + *x = gsl_ran_flat(r, 0, 1); + *y = gsl_ran_flat(r, 0, 1); +} + +void randfunc_gaus(gsl_rng *r, double *x, double *y) { + *x = 100; + *y = 100; + double sigma = 0.25; + while (fabs(*x) > 0.5 || fabs(*y) > 0.5) { + gsl_ran_bivariate_gaussian(r, sigma, sigma, 0, x, y); + } + *x += 0.5; + *y += 0.5; +} diff --git a/src/toy_setup.c b/src/toy_setup.c new file mode 100644 index 0000000..228acbc --- /dev/null +++ b/src/toy_setup.c @@ -0,0 +1,37 @@ + +#include "fracture.h" + +bool break_fuses(unsigned int width, double *fuse_thres, double *field, + bool *fuses) { + assert(fuse_thres != NULL); + assert(field != NULL); + assert(fuses != NULL); + unsigned int size = pow(width, 2); + + for (unsigned int i = 0; i < size; i++) { + if (fuses[i]) { + if (1 / fuse_thres[i] > field[i]) { + fuses[i] = 0; + } + } + } + + return true; +} + +bool gen_toy_field(unsigned int width, double strength, double *field) { + assert(field != NULL); + unsigned int size = pow(width, 2); + + for (unsigned intint i = 0; i < size; i++) { + int x = i / width + 1; + int y = i % width + 1; + if (y < (width + 1) / 2.) + field[i] = fabs((width) / 2. - x) / strength; + else + field[i] = + sqrt(pow((width) / 2. - x, 2) + pow((width) / 2. - y, 2)) / strength; + } + + return true; +} diff --git a/src/update_beta.c b/src/update_beta.c new file mode 100644 index 0000000..4c1bf65 --- /dev/null +++ b/src/update_beta.c @@ -0,0 +1,34 @@ + +#include "fracture.h" + +double f(double damage) { + assert(damage <= 1 && damage >= 0); + return sqrt(1 - sqrt(1 - damage)); + // return 0.5 - 0.68182 * (0.5 - damage); +} + +double update_beta(double beta, unsigned int width, const double *stress, + const double *damage, double bound_total) { + + double total = 0; + unsigned int num_totaled = 0; + + for (unsigned int i = 0; i < pow(width, 2); i++) { + unsigned int stress_index = + width / 4 + (width / 4 + (i / width) / 2) * width + (i % width) / 2; + double outer_damage = f(pow(fabs(stress[i]), beta)); + double inner_stress = fabs(2 * stress[stress_index] / bound_total); + + if (outer_damage > 0 && inner_stress > 0 && inner_stress != 1) { + total += log(outer_damage) / log(inner_stress); + num_totaled++; + } + } + + assert(num_totaled > 0); + assert(total == total); + + double new_beta = total / num_totaled; + + return new_beta; +} diff --git a/src/update_boundary.c b/src/update_boundary.c new file mode 100644 index 0000000..3d2a086 --- /dev/null +++ b/src/update_boundary.c @@ -0,0 +1,86 @@ + +#include "fracture.h" + +void update_boundary(finst *instance, const double *avg_field) { + int size = instance->network->num_edges; + int width = sqrt(size); + int num_verts = instance->network->num_verts; + double *boundary = (double *)instance->boundary_cond->x; + + boundary[num_verts] = 0; + boundary[num_verts + 1] = 0; + + if (instance->voltage_bound) { + for (int i = 0; i < width / 4; i++) { + int t_ind = (width + 1) * (width / 8) + width / 4 - 1 + i; + double t_val = avg_field[t_ind]; + boundary[2 * i] = t_val; + boundary[2 * i + 1] = t_val; + + int b_ind = num_verts - (width + 1) * (width / 8) - width / 4 - i; + double b_val = avg_field[b_ind]; + boundary[num_verts - 1 - 2 * i] = b_val; + boundary[num_verts - 2 - 2 * i] = b_val; + + int l_ind = (width + 1) * (width / 8) + width / 2 + width / 4 - 1 + + i * (width + 1); + double l_val = avg_field[l_ind]; + boundary[width / 2 + 2 * i * (width + 1)] = l_val; + boundary[width / 2 + (2 * i + 1) * (width + 1)] = l_val; + + int r_ind = (width + 1) * (width / 8) + width - 1 + i * (width + 1); + double r_val = avg_field[r_ind]; + boundary[width + 2 * i * (width + 1)] = r_val; + boundary[width + (2 * i + 1) * (width + 1)] = r_val; + } + } else { + const double *stress = avg_field; + + for (int i = 0; i < width / 4; i++) { + boundary[2 * i] = stress[(width / 4 - 1) * width + width / 4 + 2 * i] + + stress[(width / 4 - 1) * width + width / 4 + 1 + 2 * i]; + boundary[2 * i + 1] = + stress[(width / 4 - 1) * width + width / 4 + 2 * i] + + stress[(width / 4 - 1) * width + width / 4 + 1 + 2 * i]; + + boundary[num_verts - 2 * i - 1] = + -stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 2 * i)] - + stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 1 + 2 * i)]; + boundary[num_verts - 2 * i - 2] = + -stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 2 * i)] - + stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 1 + 2 * i)]; + + boundary[(width + 1) / 2 + 2 * i * (width + 1)] = + stress[((width / 4) + 2 * i) * width + width / 4 - 1] - + stress[((width / 4) + 2 * i + 1) * width + width / 4 - 1]; + boundary[(width + 1) / 2 + (2 * i + 1) * (width + 1)] = + stress[((width / 4) + 2 * i) * width + width / 4 - 1] - + stress[((width / 4) + 2 * i + 1) * width + width / 4 - 1]; + + boundary[(width + 1) / 2 + 2 * i * (width + 1) + width / 2] = + stress[((width / 4) + 2 * i) * width + width / 4 + width / 2] - + stress[((width / 4) + 2 * i + 1) * width + width / 4 + width / 2]; + boundary[(width + 1) / 2 + (2 * i + 1) * (width + 1) + width / 2] = + stress[((width / 4) + 2 * i) * width + width / 4 + width / 2] - + stress[((width / 4) + 2 * i + 1) * width + width / 4 + width / 2]; + } + + double s_total = 0; + double total = 0; + for (int i = 0; i < num_verts; i++) { + if (boundary[i] != boundary[i]) { + boundary[i] = 0; + } else { + total += fabs(boundary[i]); + s_total += boundary[i]; + } + } + + // assert(fabs(s_total) < inf); + boundary[num_verts] -= s_total; + + for (int i = 0; i < num_verts + 1; i++) { + boundary[i] /= total; + } + } +} diff --git a/src/update_factor.c b/src/update_factor.c new file mode 100644 index 0000000..6b8e05e --- /dev/null +++ b/src/update_factor.c @@ -0,0 +1,44 @@ + +#include "fracture.h" + +bool update_factor(cholmod_factor *factor, unsigned int v1, unsigned int v2, + cholmod_common *c) { + int n = factor->n; + assert(v1 < n); + assert(v2 < n); + + cholmod_sparse *update_mat = + CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); + + unsigned int v3, v4; + v3 = v1 < v2 ? v1 : v2; + v4 = v1 > v2 ? v1 : v2; + + for (int i = 0; i < n; i++) { + if (i <= v3) + ((CHOL_INT *)update_mat->p)[i] = 0; + else + ((CHOL_INT *)update_mat->p)[i] = 2; + } + ((CHOL_INT *)update_mat->p)[n] = 2; + ((CHOL_INT *)update_mat->i)[0] = v3; + ((CHOL_INT *)update_mat->i)[1] = v4; + ((double *)update_mat->x)[0] = 1; + ((double *)update_mat->x)[1] = -1; + + // assert(CHOL_F(check_sparse)(update_mat, c)); + + cholmod_sparse *perm_update_mat = CHOL_F(submatrix)( + update_mat, factor->Perm, factor->n, NULL, -1, true, true, c); + + // assert(CHOL_F(check_sparse)(perm_update_mat, c)); + + CHOL_F(updown)(false, perm_update_mat, factor, c); + + // assert(CHOL_F(check_factor)(factor, c)); + + CHOL_F(free_sparse)(&perm_update_mat, c); + CHOL_F(free_sparse)(&update_mat, c); + + return true; +} diff --git a/src/voronoi_bound_ini.c b/src/voronoi_bound_ini.c new file mode 100644 index 0000000..38f41cc --- /dev/null +++ b/src/voronoi_bound_ini.c @@ -0,0 +1,79 @@ + +#include "fracture.h" + +void voronoi_bound_ini(finst *instance, double *square_bound, + unsigned int width) { + unsigned int square_num_verts = 2 * ((width + 1) / 2) * (width / 2 + 1); + double total1 = 0; + double total2 = 0; + double total3 = 0; + double total4 = 0; + for (unsigned int i = 0; i < instance->network->num_bounds; i++) { + for (unsigned int j = 0; j < instance->network->bound_inds[i + 1] - + instance->network->bound_inds[i]; + j++) { + double x = instance->network + ->vert_coords[2 * + instance->network->bound_verts + [instance->network->bound_inds[i] + j]]; + double y = + instance->network + ->vert_coords[2 * + instance->network->bound_verts + [instance->network->bound_inds[i] + j] + + 1]; + + unsigned int vw = (width + 1) / 2; + unsigned int xp = (unsigned int)(x * vw); + unsigned int yp = (unsigned int)(y * vw); + + if (i == 0) { + ((double *)instance->boundary_cond + ->x)[instance->network + ->bound_verts[instance->network->bound_inds[i] + j]] = + square_bound[xp]; + total1 += square_bound[xp]; + } + if (i == 1) { + ((double *)instance->boundary_cond + ->x)[instance->network + ->bound_verts[instance->network->bound_inds[i] + j]] = + square_bound[square_num_verts - vw - 1 + xp]; + total2 += square_bound[square_num_verts - vw - 1 + xp]; + } + if (i == 2) { + ((double *)instance->boundary_cond + ->x)[instance->network + ->bound_verts[instance->network->bound_inds[i] + j]] = + square_bound[(width + 1) / 2 + yp * (width + 1)]; + total3 += square_bound[(width + 1) / 2 + yp * (width + 1)]; + } + if (i == 3) { + ((double *)instance->boundary_cond + ->x)[instance->network + ->bound_verts[instance->network->bound_inds[i] + j]] = + square_bound[(width + 1) / 2 + yp * (width + 1) + width / 2]; + total4 += square_bound[(width + 1) / 2 + yp * (width + 1) + width / 2]; + } + } + } + + ((double *)instance->boundary_cond->x)[instance->network->num_verts] = + - total1 / 2 - total2 / 2; + ((double *)instance->boundary_cond->x)[instance->network->num_verts + 1] = + - total1 / 2 - total2 / 2; + ((double *)instance->boundary_cond->x)[instance->network->num_verts + 2] = + -total3; + ((double *)instance->boundary_cond->x)[instance->network->num_verts + 3] = + -total4; + /* + ((double *)instance->boundary_cond->x)[instance->network->num_verts] = + 0; + ((double *)instance->boundary_cond->x)[instance->network->num_verts + 1] = + 0; + ((double *)instance->boundary_cond->x)[instance->network->num_verts + 2] = + 0; + ((double *)instance->boundary_cond->x)[instance->network->num_verts + 3] = + 0; + */ +} |