summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--makefile_hal4
-rw-r--r--src/break_edge.c2
-rw-r--r--src/cracked_voronoi_fracture.c427
-rw-r--r--src/cracking_ini.c24
-rw-r--r--src/fracture.h5
-rw-r--r--src/fracture_network.c6
-rw-r--r--src/gen_laplacian.c7
-rw-r--r--src/ini_network.c4
-rw-r--r--src/instance.c19
-rw-r--r--src/voronoi_bound_ini.c98
10 files changed, 506 insertions, 90 deletions
diff --git a/makefile_hal b/makefile_hal
index 0b2f6ff..4559e77 100644
--- a/makefile_hal
+++ b/makefile_hal
@@ -1,10 +1,10 @@
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
+LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -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 = current_scaling course_grain_square corr_test homo_voronoi_fracture homo_square_fracture compare_voronoi_fracture
+BIN = current_scaling course_grain_square corr_test homo_voronoi_fracture homo_square_fracture compare_voronoi_fracture cracked_voronoi_fracture
all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
diff --git a/src/break_edge.c b/src/break_edge.c
index 0ba1ad7..72ccf3e 100644
--- a/src/break_edge.c
+++ b/src/break_edge.c
@@ -8,7 +8,7 @@ bool break_edge(finst *instance, unsigned int edge, cholmod_common *c) {
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);
+ if (instance->factor != NULL) update_factor(instance->factor, v1, v2, c);
if (instance->network->boundary != TORUS_BOUND) {
unsigned int w1 = instance->network->edges_to_verts[2 * edge];
diff --git a/src/cracked_voronoi_fracture.c b/src/cracked_voronoi_fracture.c
new file mode 100644
index 0000000..532c3c3
--- /dev/null
+++ b/src/cracked_voronoi_fracture.c
@@ -0,0 +1,427 @@
+
+
+#include "fracture.h"
+
+int main(int argc, char *argv[]) {
+ int opt;
+
+ // defining variables to be (potentially) set by command line flags
+ uint8_t filename_len;
+ uint32_t N;
+ uint_t L;
+ double beta, inf, cutoff, crack_len;
+ bool include_breaking, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
+ save_crit_stress, save_stress_field, save_voltage_field, save_toughness, save_conductivity,
+ save_damage, save_damage_field;
+ bound_t boundary;
+
+ filename_len = 100;
+
+ N = 100;
+ L = 16;
+ crack_len = 0.;
+ beta = .3;
+ inf = 1e10;
+ cutoff = 1e-9;
+ boundary = FREE_BOUND;
+ include_breaking = false;
+ save_cluster_dist = false;
+ use_voltage_boundaries = false;
+ use_dual = false;
+ save_network = false;
+ save_crit_stress = false;
+ save_stress_field = false;
+ save_voltage_field = false;
+ save_damage = false;
+ save_damage_field = false;
+ save_conductivity = false;
+ save_toughness = false;
+
+ uint8_t bound_i;
+ char boundc2 = 'f';
+
+ while ((opt = getopt(argc, argv, "n:L:b:B:dVcoNsCrtDSvel:")) != -1) {
+ switch (opt) {
+ case 'n':
+ N = atoi(optarg);
+ break;
+ case 'L':
+ L = atoi(optarg);
+ break;
+ case 'b':
+ beta = atof(optarg);
+ break;
+ case 'l':
+ crack_len = atof(optarg);
+ break;
+ case 'B':
+ bound_i = atoi(optarg);
+ switch (bound_i) {
+ 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;
+ case 3:
+ boundary = EMBEDDED_BOUND;
+ boundc2 = 'e';
+ use_dual = true;
+ use_voltage_boundaries = true;
+ break;
+ default:
+ printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n");
+ exit(EXIT_FAILURE);
+ }
+ break;
+ case 'd':
+ save_damage = true;
+ break;
+ case 'e':
+ save_damage_field = true;
+ break;
+ case 'V':
+ use_voltage_boundaries = true;
+ break;
+ case 'D':
+ use_dual = 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 'S':
+ save_stress_field = true;
+ break;
+ case 'v':
+ save_voltage_field = 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);
+ }
+
+ uint_t voronoi_max_verts, c_dist_size, a_dist_size;
+
+ voronoi_max_verts = 4 * pow(L, 2);
+ c_dist_size = voronoi_max_verts;
+ a_dist_size = voronoi_max_verts;
+
+ if (voronoi_max_verts > CINT_MAX) {
+ exit(EXIT_FAILURE);
+ }
+
+ // define arrays for saving cluster and avalanche distributions
+ uint32_t *cluster_size_dist;
+ uint32_t *avalanche_size_dist;
+ char *c_filename;
+ char *a_filename;
+ if (save_cluster_dist) {
+ cluster_size_dist =
+ (uint32_t *)malloc(c_dist_size * sizeof(uint32_t));
+ avalanche_size_dist =
+ (uint32_t *)malloc(a_dist_size * sizeof(uint32_t));
+
+ c_filename = (char *)malloc(filename_len * sizeof(char));
+ a_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(c_filename, filename_len, "cstr_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta);
+ snprintf(a_filename, filename_len, "avln_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta);
+
+ FILE *cluster_out = fopen(c_filename, "rb");
+ FILE *avalanche_out = fopen(a_filename, "rb");
+
+ if (cluster_out != NULL) {
+ fread(cluster_size_dist, sizeof(uint32_t), c_dist_size, cluster_out);
+ fclose(cluster_out);
+ }
+ if (avalanche_out != NULL) {
+ fread(avalanche_size_dist, sizeof(uint32_t), a_dist_size, avalanche_out);
+ fclose(avalanche_out);
+ }
+ }
+
+ double *crit_stress;
+ if (save_crit_stress) {
+ crit_stress = (double *)malloc(N * sizeof(double));
+ }
+
+ double *conductivity;
+ if (save_conductivity) {
+ conductivity = (double *)malloc(N * sizeof(double));
+ }
+
+ // define arrays for saving damage distributions
+ uint32_t *damage;
+ char *d_filename;
+ if (save_damage) {
+ damage =
+ (uint32_t *)malloc(a_dist_size * sizeof(uint32_t));
+
+ d_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(d_filename, filename_len, "damg_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta);
+
+ FILE *damage_out = fopen(d_filename, "rb");
+
+ if (damage_out != NULL) {
+ fread(damage, sizeof(uint32_t), a_dist_size, damage_out);
+ fclose(damage_out);
+ }
+ }
+
+ 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 (uint32_t i = 0; i < N; i++) {
+ printf("\033[F\033[JFRACTURE: %0*d / %d\n", (uint8_t)log10(N) + 1, i + 1, N);
+
+ fnet *network = ini_voronoi_network(L, boundary, use_dual, genfunc_hyperuniform, &c);
+ finst *perm_instance = create_instance(network, inf, use_voltage_boundaries, false, &c);
+ if (boundary == EMBEDDED_BOUND) {
+ voronoi_bound_ini(perm_instance, L, crack_len);
+ }
+ gen_voro_crack(perm_instance, crack_len, &c);
+ finish_instance(perm_instance, &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);
+ free_instance(instance, &c);
+ free(fuse_thres);
+
+ uint_t max_pos = 0;
+ double max_val = 0;
+
+ for (uint_t 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);
+
+ uint_t av_size = 0;
+ double cur_val = 0;
+ for (uint_t 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 (uint_t 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[max_pos]++;
+ }
+
+ if (save_cluster_dist) {
+ uint_t *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c);
+ for (uint_t j = 0; j < tmp_instance->network->num_dual_verts; j++) {
+ cluster_size_dist[j] += tmp_cluster_dist[j];
+ }
+ free(tmp_cluster_dist);
+ }
+
+ if (save_network) {
+ FILE *net_out = fopen("network.txt", "w");
+ for (uint_t 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 (uint_t 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 (uint_t 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 (uint_t 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 (uint_t 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 (uint_t 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");
+
+ if (save_cluster_dist) {
+ FILE *cluster_out = fopen(c_filename, "wb");
+ FILE *avalanche_out = fopen(a_filename, "wb");
+
+ fwrite(cluster_size_dist, sizeof(uint32_t), c_dist_size, cluster_out);
+ fwrite(avalanche_size_dist, sizeof(uint32_t), a_dist_size, avalanche_out);
+
+ fclose(cluster_out);
+ fclose(avalanche_out);
+
+ free(c_filename);
+ free(a_filename);
+ free(cluster_size_dist);
+ free(avalanche_size_dist);
+ }
+
+ if (save_conductivity) {
+ char *cond_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(cond_filename, filename_len, "cond_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta);
+ FILE *cond_file = fopen(cond_filename, "ab");
+ fwrite(conductivity, sizeof(double), N, cond_file);
+ 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.dat", boundc, boundc2, L, beta);
+ FILE *tough_file = fopen(tough_filename, "ab");
+ fwrite(toughness, sizeof(double), N, tough_file);
+ fclose(tough_file);
+ free(tough_filename);
+ free(toughness);
+ }
+
+ if (save_damage) {
+ FILE *hdam_file = fopen(d_filename, "wb");
+ fwrite(damage, sizeof(uint32_t), a_dist_size, hdam_file);
+ fclose(hdam_file);
+ free(d_filename);
+ free(damage);
+ }
+
+ if (include_breaking) {
+ fclose(break_out);
+ }
+
+ if (save_crit_stress) {
+ char *str_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(str_filename, filename_len, "strs_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta);
+ FILE *str_file = fopen(str_filename, "ab");
+ fwrite(crit_stress, sizeof(double), N, str_file);
+ fclose(str_file);
+ free(str_filename);
+ free(crit_stress);
+ }
+
+ CHOL_F(finish)(&c);
+
+ return 0;
+}
diff --git a/src/cracking_ini.c b/src/cracking_ini.c
index 9ab4a9d..988af5b 100644
--- a/src/cracking_ini.c
+++ b/src/cracking_ini.c
@@ -31,6 +31,30 @@ double *gen_fuse_thres(unsigned int num_edges, double *edge_coords, double beta,
return fuse_thres;
}
+void gen_voro_crack(finst *instance, double crack_len, cholmod_common *c) {
+ for (uint_t i = 0; i < instance->network->num_edges; i++) {
+ uint_t v1, v2;
+ double v1x, v1y, v2x, v2y, dx, dy;
+
+ v1 = instance->network->edges_to_verts[2 * i];
+ v2 = instance->network->edges_to_verts[2 * i + 1];
+
+ v1x = instance->network->vert_coords[2 * v1];
+ v1y = instance->network->vert_coords[2 * v1 + 1];
+ v2x = instance->network->vert_coords[2 * v2];
+ v2y = instance->network->vert_coords[2 * v2 + 1];
+
+ dx = v1x - v2x;
+ dy = v1y - v2y;
+
+ if (((v1y > 0.5 && v2y < 0.5) || (v1y < 0.5 && v2y > 0.5)) && fabs(dy) < 0.5) {
+ if (v1x + dx / dy * (v1y - 0.5) <= crack_len) {
+ break_edge(instance, i, c);
+ }
+ }
+ }
+}
+
bool gen_crack(finst *instance, double crack_len, double crack_width,
cholmod_common *c) {
assert(instance != NULL);
diff --git a/src/fracture.h b/src/fracture.h
index 91b2bfc..93e2ed1 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -175,11 +175,12 @@ 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);
+void voronoi_bound_ini(finst *instance, uint_t L, double crack_len);
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);
+
+void gen_voro_crack(finst *instance, double crack_len, cholmod_common *c);
diff --git a/src/fracture_network.c b/src/fracture_network.c
index fe520fd..3f06104 100644
--- a/src/fracture_network.c
+++ b/src/fracture_network.c
@@ -53,7 +53,11 @@ break_data *fracture_network(finst *instance, double *fuse_thres,
break_edge(instance, last_broke, c);
- if (instance->num_components > 1) {
+ if (instance->num_components > 1 && instance->network->boundary == TORUS_BOUND) {
+ break;
+ }
+
+ if (instance->marks[num_verts] != instance->marks[num_verts + 1] && instance->network->boundary != TORUS_BOUND) {
break;
}
}
diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c
index d8f0c6d..97a2c9d 100644
--- a/src/gen_laplacian.c
+++ b/src/gen_laplacian.c
@@ -159,14 +159,15 @@ cholmod_sparse *gen_laplacian(const finst *instance, cholmod_common *c,
}
if (network->boundary != TORUS_BOUND) {
+ if (network->boundary != EMBEDDED_BOUND) acoo[0]++;
+
if (voltage_bound) {
- for (unsigned int i = 0; i < 2; i++) {
+ for (unsigned int i = 0; i < num_bounds; 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;
@@ -181,7 +182,7 @@ cholmod_sparse *gen_laplacian(const finst *instance, cholmod_common *c,
unsigned int start = num_gverts;
- for (unsigned int i = 0; i < num_bounds; i++) {
+ for (unsigned int i = 0; i < 2; 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;
diff --git a/src/ini_network.c b/src/ini_network.c
index 263b1b2..d80f43b 100644
--- a/src/ini_network.c
+++ b/src/ini_network.c
@@ -604,7 +604,7 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
break;
}
case EMBEDDED_BOUND: {
- num_bounds = 2;
+ num_bounds = 4;
bound_inds = (unsigned int *)malloc(5 * sizeof(unsigned int));
bound_verts = (unsigned int *)malloc(2 * L * sizeof(unsigned int));
for (unsigned int i = 0; i < 5; i++) bound_inds[i] = i * L / 2;
@@ -622,7 +622,7 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
network->num_verts = num_verts;
network->edges_to_verts_break = edges;
network->edges_to_verts = edges;
- network->break_dim = num_verts + num_bounds;
+ network->break_dim = num_verts + 2;
}
}
diff --git a/src/instance.c b/src/instance.c
index 098ed92..cc33591 100644
--- a/src/instance.c
+++ b/src/instance.c
@@ -14,7 +14,18 @@ finst *create_instance(fnet *network, double inf, bool voltage_bound,
instance->voltage_bound = voltage_bound;
instance->boundary_cond = CHOL_F(zeros)(
network->break_dim, 1, CHOLMOD_REAL, c);
- if (network->boundary != TORUS_BOUND) {
+ if (network->boundary == TORUS_BOUND) {
+ 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;
+ } else if (network->boundary == EMBEDDED_BOUND) {
+ // do nothing
+ for (unsigned int i = 0; i < network->bound_inds[1]; i++) {
+ ((double *)instance->boundary_cond->x)[network->bound_verts[i]] =1;
+ }
+ } else {
if (voltage_bound) {
for (unsigned int i = 0; i < network->bound_inds[1]; i++) {
((double *)instance->boundary_cond->x)[network->bound_verts[i]] =1;
@@ -24,12 +35,6 @@ finst *create_instance(fnet *network, double inf, bool voltage_bound,
((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;
}
if (network->boundary != TORUS_BOUND) instance->adjacency = gen_adjacency(instance, false, false, 0, c);
diff --git a/src/voronoi_bound_ini.c b/src/voronoi_bound_ini.c
index 38f41cc..d201093 100644
--- a/src/voronoi_bound_ini.c
+++ b/src/voronoi_bound_ini.c
@@ -1,79 +1,33 @@
#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];
+double th_p(double x, double y, double th) {
+ if (x >= 0 && y >= 0) return th;
+ if (x < 0 && y >= 0) return M_PI - th;
+ if (x < 0 && y < 0) return th - M_PI;
+ if (x >= 0 && y < 0) return -th;
+}
- unsigned int vw = (width + 1) / 2;
- unsigned int xp = (unsigned int)(x * vw);
- unsigned int yp = (unsigned int)(y * vw);
+double u_y(double x, double y) {
+ double r = sqrt(pow(x, 2) + pow(y, 2));
+ double th = th_p(x, y, atan(fabs(y / x)));
- 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];
- }
- }
- }
+ return sqrt(r) * sin(th / 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;
- */
+void voronoi_bound_ini(finst *instance, uint_t L, double crack_len) {
+ double *bound = (double *)instance->boundary_cond->x;
+ for (uint_t i = 0; i < L / 2; i++) {
+ double x1, y1, x2, y2, x3, y3, x4, y4;
+ x1 = (2. * i + 1.) / L - crack_len; y1 = 0.5 - 1.;
+ x2 = (2. * i + 1.) / L - crack_len; y2 = 0.5 - 0.;
+ y3 = (2. * i + 1.) / L - 0.5; x3 = 0.5 - 1.;
+ y4 = (2. * i + 1.) / L - 0.5; x4 = 0.5 - 0.;
+
+ bound[i] = u_y(x1, y1);
+ bound[L / 2 + i] = u_y(x2, y2);
+ bound[L + i] = u_y(x3, y3);
+ bound[3 * L / 2 + i] = u_y(x4, y4);
+ }
}
+