From 0cf3090b05aa4c35ba20f3773db4f95d4600b55e Mon Sep 17 00:00:00 2001 From: Jaron Date: Tue, 8 Nov 2016 10:18:15 -0500 Subject: restarted all data generation to correct for errors, made certain data (stress, energy, thresholds) longs --- src/data.c | 4 +- src/fracture.c | 106 +++++++------------------------- src/fracture.h | 4 +- src/long_anal_process.c | 156 ++++++++++++++++++++++++++++++++++++++++++++++++ src/net_fracture.c | 4 +- 5 files changed, 182 insertions(+), 92 deletions(-) create mode 100644 src/long_anal_process.c (limited to 'src') diff --git a/src/data.c b/src/data.c index b54236d..2047c44 100644 --- a/src/data.c +++ b/src/data.c @@ -10,7 +10,7 @@ data_t *data_create(uint_t ne) { data->break_list = (uint_t *)malloc(ne * sizeof(uint_t)); assert(data->break_list != NULL); - data->extern_field = (double *)malloc(ne * sizeof(double)); + data->extern_field = (long double *)malloc(ne * sizeof(long double)); assert(data->extern_field != NULL); data->conductivity = (double *)malloc(ne * sizeof(double)); @@ -26,7 +26,7 @@ void data_free(data_t *data) { free(data); } -void data_update(data_t *data, uint_t last_broke, double strength, double conductivity) { +void data_update(data_t *data, uint_t last_broke, long 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; diff --git a/src/fracture.c b/src/fracture.c index 7ee82fb..59bd424 100644 --- a/src/fracture.c +++ b/src/fracture.c @@ -10,8 +10,8 @@ int main(int argc, char *argv[]) { uint32_t N; uint_t L, refactor_every; double beta, inf, cutoff, crack_len; - bool crack_growth_crit, save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network, - save_crit_stress, save_stress_field, save_voltage_field, save_toughness, save_energy, save_conductivity, + bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network, + save_crit_stress, save_stress_field, save_voltage_field, save_energy, save_conductivity, save_damage, save_damage_field, save_threshold; bound_t boundary; lattice_t lattice; @@ -33,7 +33,6 @@ int main(int argc, char *argv[]) { cutoff = 1e-9; boundary = FREE_BOUND; lattice = VORONOI_LATTICE; - crack_growth_crit = false; save_data = false; save_cluster_dist = false; use_voltage_boundaries = false; @@ -45,7 +44,6 @@ int main(int argc, char *argv[]) { save_damage = false; save_damage_field = false; save_conductivity = false; - save_toughness = false; save_energy = true; save_threshold = false; @@ -59,7 +57,7 @@ int main(int argc, char *argv[]) { // get commandline options - while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:gTER:")) != -1) { + while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDSvel:TER:")) != -1) { switch (opt) { case 'n': N = atoi(optarg); @@ -119,9 +117,6 @@ int main(int argc, char *argv[]) { exit(EXIT_FAILURE); } break; - case 'g': - crack_growth_crit = true; - break; case 'd': save_damage = true; break; @@ -156,9 +151,6 @@ int main(int argc, char *argv[]) { case 'r': save_conductivity = true; break; - case 't': - save_toughness = true; - break; case 'E': save_energy = true; break; @@ -222,9 +214,9 @@ int main(int argc, char *argv[]) { } } - double *crit_stress; + long double *crit_stress; if (save_crit_stress) { - crit_stress = (double *)malloc(N * sizeof(double)); + crit_stress = (long double *)malloc(N * sizeof(long double)); } double *stress_field; @@ -268,19 +260,14 @@ int main(int argc, char *argv[]) { } } - double *toughness; - if (save_toughness) { - toughness = (double *)malloc(N * sizeof(double)); - } - - double *energy; + long double *energy; if (save_energy) { - energy = (double *)malloc(N * sizeof(double)); + energy = (long double *)malloc(N * sizeof(long double)); } - double *thresholds; + long double *thresholds; if (save_threshold) { - thresholds = (double *)malloc(N * sizeof(double)); + thresholds = (long double *)malloc(N * sizeof(long double)); } @@ -311,7 +298,7 @@ int main(int argc, char *argv[]) { net_free(tmp_net, &c); uint_t max_pos = 0; - double max_val = 0; + long double max_val = 0; double cond0; { @@ -321,7 +308,7 @@ int main(int argc, char *argv[]) { } for (uint_t j = 0; j < data->num_broken; j++) { - double val = data->extern_field[j]; + long double val = data->extern_field[j]; if (val > max_val) { max_pos = j; @@ -330,41 +317,15 @@ int main(int argc, char *argv[]) { } uint_t av_size = 0; - double cur_val = 0; - uint_t notch_mark = 0; - - if (crack_growth_crit) { - uint_t n_edge = UINT_MAX; - for (uint_t j = 0; j < g->ne; j++) { - // right now the only fuses broken are from a notch - if (net->fuses[j]) { - n_edge = j; - break; - } - } - - if (n_edge < UINT_MAX) { - uint_t n_vert = g->dev[2 * n_edge]; - notch_mark = net->dual_marks[n_vert]; - } - } + long double cur_val = 0; for (uint_t j = 0; j < max_pos; j++) { uint_t next_broken = data->break_list[j]; - bool attached_to_notch = net->dual_marks[g->dev[2*next_broken]] == notch_mark; - bool grew_forward = g->ex[2 * next_broken] > crack_len; - - if (crack_growth_crit && attached_to_notch && grew_forward) { - printf("whoops\n\n"); - max_pos = j; - break; - } - bool refactor = ((j + 1) % refactor_every) == 0; break_edge(net, next_broken, &c, refactor); - double val = data->extern_field[j]; + long double val = data->extern_field[j]; if (save_cluster_dist) { if (val < cur_val) { av_size++; @@ -422,12 +383,12 @@ int main(int argc, char *argv[]) { } if (save_energy) { - double tmp_energy = 0; + long double tmp_energy = 0; if (max_pos > 0) { - double sigma1 = data->extern_field[0]; + long double sigma1 = data->extern_field[0]; double cond1 = cond0; for (uint_t j = 0; j < max_pos - 1; j++) { - double sigma2 = data->extern_field[j+1]; + long double sigma2 = data->extern_field[j+1]; double cond2 = data->conductivity[j]; if (sigma2 > sigma1) { tmp_energy += 0.5 * gsl_pow_2(sigma1) * (1 - cond2 / cond1) / cond1; @@ -439,23 +400,6 @@ int main(int argc, char *argv[]) { energy[i] = tmp_energy; } - if (save_toughness) { - double tmp_toughness = 0; - if (max_pos > 0) { - double sigma1 = data->extern_field[0]; - double epsilon1 = sigma1 / data->conductivity[0]; - for (uint_t j = 0; j < max_pos - 1; j++) { - double sigma2 = data->extern_field[j+1]; - double epsilon2 = sigma2 / data->conductivity[j+1]; - if (epsilon2 > epsilon1) { - tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2; - sigma1 = sigma2; epsilon1 = epsilon2; - } - } - } - toughness[i] = tmp_toughness; - } - if (save_threshold) { thresholds[i] = net->thres[data->break_list[max_pos]]; } @@ -499,7 +443,7 @@ int main(int argc, char *argv[]) { if (save_data) { for (uint_t j = 0; j < data->num_broken; j++) { - fprintf(data_out, "%u %g %g ", data->break_list[j], + fprintf(data_out, "%u %Lg %g ", data->break_list[j], data->extern_field[j], data->conductivity[j]); } fprintf(data_out, "\n"); @@ -566,21 +510,11 @@ int main(int argc, char *argv[]) { free(conductivity); } - if (save_toughness) { - char *tough_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(tough_filename, filename_len, "tuff_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - FILE *tough_file = fopen(tough_filename, "ab"); - fwrite(toughness, sizeof(double), N, tough_file); - fclose(tough_file); - free(tough_filename); - free(toughness); - } - if (save_energy) { char *tough_filename = (char *)malloc(filename_len * sizeof(char)); snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); FILE *tough_file = fopen(tough_filename, "ab"); - fwrite(energy, sizeof(double), N, tough_file); + fwrite(energy, sizeof(long double), N, tough_file); fclose(tough_file); free(tough_filename); free(energy); @@ -590,7 +524,7 @@ int main(int argc, char *argv[]) { char *thres_filename = (char *)malloc(filename_len * sizeof(char)); snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); FILE *thres_file = fopen(thres_filename, "ab"); - fwrite(thresholds, sizeof(double), N, thres_file); + fwrite(thresholds, sizeof(long double), N, thres_file); fclose(thres_file); free(thres_filename); free(thresholds); @@ -612,7 +546,7 @@ int main(int argc, char *argv[]) { char *str_filename = (char *)malloc(filename_len * sizeof(char)); snprintf(str_filename, filename_len, "strs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); FILE *str_file = fopen(str_filename, "ab"); - fwrite(crit_stress, sizeof(double), N, str_file); + fwrite(crit_stress, sizeof(long double), N, str_file); fclose(str_file); free(str_filename); free(crit_stress); diff --git a/src/fracture.h b/src/fracture.h index f122869..556e1e5 100644 --- a/src/fracture.h +++ b/src/fracture.h @@ -89,7 +89,7 @@ typedef struct { uint_t num_broken; uint_t *break_list; double *conductivity; - double *extern_field; + long double *extern_field; } data_t; intptr_t *run_voronoi(uint_t num_coords, double *coords, bool periodic, double xmin, double xmax, double ymin, double ymax); @@ -173,7 +173,7 @@ cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_co data_t *data_create(uint_t num_edges); void data_free(data_t *data); -void data_update(data_t *data, uint_t last_broke, double strength, double conductivity); +void data_update(data_t *data, uint_t last_broke, long double strength, double conductivity); graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c); diff --git a/src/long_anal_process.c b/src/long_anal_process.c new file mode 100644 index 0000000..ba29152 --- /dev/null +++ b/src/long_anal_process.c @@ -0,0 +1,156 @@ + +#include "fracture.h" +#include +#include +#include +#include +#include +#include + +void get_mean(uint_t len, long double *data, long double result[2]) { + long double total = 0; + + for (uint_t i = 0; i < len; i++) { + total += data[i]; + } + + long double mean = total / len; + long double total_sq_diff = 0; + + for (uint_t i = 0; i < len; i++) { + total_sq_diff += pow(mean - data[i], 2); + } + + long double mean_err = sqrt(total_sq_diff) / len; + + result[0] = mean; + result[1] = mean_err; +} + +void get_mom(uint_t len, uint_t n, long double *data, long double mean[2], long double result[2]) { + long double total_n_diff = 0; + long double total_n2_diff = 0; + + for (uint_t i = 0; i < len; i++) { + total_n_diff += pow(fabsl(mean[0] - data[i]), n); + total_n2_diff += pow(fabsl(mean[0] - data[i]), 2 * n); + } + + long double mom = total_n_diff / len; + long double mom_err = sqrt(total_n2_diff) / len; + + result[0] = mom; + result[1] = mom_err; +} + +int main(int argc, char *argv[]) { + uint_t nc = argc - 1; + uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t)); + double *betas_c = (double *)malloc(nc * sizeof(double)); + long double *vals_c1 = (long double *)malloc(nc * sizeof(long double)); + long double *errors_c1 = (long double *)malloc(nc * sizeof(long double)); + long double *vals_c2 = (long double *)malloc(nc * sizeof(long double)); + long double *errors_c2 = (long double *)malloc(nc * sizeof(long double)); + long double *vals_c3 = (long double *)malloc(nc * sizeof(long double)); + long double *errors_c3 = (long double *)malloc(nc * sizeof(long double)); + long double *vals_c4 = (long double *)malloc(nc * sizeof(long double)); + long double *errors_c4 = (long double *)malloc(nc * sizeof(long double)); + + char *out_filename = (char *)malloc(100 * sizeof(char)); + uint_t out_filename_len = 0; + + for (uint_t i = 0; i < nc; i++) { + char *fn = argv[1 + i]; + char *buff = (char *)malloc(20 * sizeof(char)); + uint_t pos = 0; uint_t c = 0; + while (c < 5) { + if (fn[pos] == '_') { + c++; + } + if (i == 0) { + out_filename[pos] = fn[pos]; + out_filename_len++; + } + pos++; + } + c = 0; + while (fn[pos] != '_') { + buff[c] = fn[pos]; + pos++; + c++; + } + buff[c] = '\0'; + Ls_c[i] = atoi(buff); + c = 0; + pos++; + while (fn[pos] != '_') { + buff[c] = fn[pos]; + pos++; + c++; + } + buff[c] = '\0'; + betas_c[i] = atof(buff); + free(buff); + + struct stat info; + stat(fn, &info); + uint_t num_bytes = info.st_size; + uint_t dist_len = (num_bytes * sizeof(char)) / sizeof(long double); + + long double *dist = malloc(dist_len * sizeof(long double)); + FILE *dist_file = fopen(fn, "rb"); + fread(dist, sizeof(long double), dist_len, dist_file); + fclose(dist_file); + { + long double mom1[2]; + get_mean(dist_len, dist, mom1); + vals_c1[i] = mom1[0]; + errors_c1[i] = mom1[1]; + + long double mom2[2]; + get_mom(dist_len, 2, dist, mom1, mom2); + vals_c2[i] = mom2[0]; + errors_c2[i] = mom2[1]; + + long double mom3[2]; + get_mom(dist_len, 3, dist, mom1, mom3); + vals_c3[i] = mom3[0]; + errors_c3[i] = mom3[1]; + + long double mom4[2]; + get_mom(dist_len, 4, dist, mom1, mom4); + vals_c4[i] = mom4[0]; + errors_c4[i] = mom4[1]; + } + free(dist); + } + + out_filename[out_filename_len-1] = '.'; + out_filename[out_filename_len] = 't'; + out_filename[out_filename_len+1] = 'x'; + out_filename[out_filename_len+2] = 't'; + out_filename[out_filename_len+3] = '\0'; + + FILE *out_file = fopen(out_filename, "w"); + + for (uint_t i = 0; i < nc; i++) { + fprintf(out_file, "%u %g %Lg %Lg %Lg %Lg %Lg %Lg %Lg %Lg\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]); + } + + fclose(out_file); + + + free(Ls_c); + free(betas_c); + free(vals_c1); + free(errors_c1); + free(vals_c2); + free(errors_c2); + free(vals_c3); + free(errors_c3); + free(vals_c4); + free(errors_c4); + + return 0; +} + diff --git a/src/net_fracture.c b/src/net_fracture.c index 3bfcd83..b5ab524 100644 --- a/src/net_fracture.c +++ b/src/net_fracture.c @@ -46,7 +46,7 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff, uint_t refact uint_t last_broke = get_next_broken(net, currents, cutoff); - double sim_current; + long double sim_current; if (net->voltage_bound) { sim_current = conductivity; @@ -54,7 +54,7 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff, uint_t refact sim_current = 1; } - data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] / currents[last_broke]), conductivity); + data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] / ((long double)currents[last_broke])), conductivity); free(voltages); free(currents); -- cgit v1.2.3-70-g09d2