From fd14c5e39d962be94a1f68b0d4cacb7a4aa9c3e7 Mon Sep 17 00:00:00 2001 From: pants Date: Fri, 2 Sep 2016 15:24:34 -0400 Subject: embedded systems with crack fully supported --- src/break_edge.c | 2 +- src/cracked_voronoi_fracture.c | 427 +++++++++++++++++++++++++++++++++++++++++ src/cracking_ini.c | 24 +++ src/fracture.h | 5 +- src/fracture_network.c | 6 +- src/gen_laplacian.c | 7 +- src/ini_network.c | 4 +- src/instance.c | 19 +- src/voronoi_bound_ini.c | 98 +++------- 9 files changed, 504 insertions(+), 88 deletions(-) create mode 100644 src/cracked_voronoi_fracture.c (limited to 'src') 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); + } } + -- cgit v1.2.3-70-g09d2