From 07906baa42470bad14d2c40f57967691f6118969 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 1 Nov 2018 12:33:37 -0400 Subject: revamped and simplied fracture code with c++ --- src/CMakeLists.txt | 7 + src/anal_process.c | 135 ------ src/big_anal_process.c | 158 ------- src/cen_anal_process.c | 157 ------- src/corr_test.c | 27 -- src/fracture.c | 1068 ----------------------------------------------- src/fracture.cpp | 59 +++ src/fracture.h | 123 ------ src/long_anal_process.c | 158 ------- src/measurements.cpp | 248 +++++++++++ src/measurements.hpp | 54 +++ src/randutils | 1 + 12 files changed, 369 insertions(+), 1826 deletions(-) create mode 100644 src/CMakeLists.txt delete mode 100644 src/anal_process.c delete mode 100644 src/big_anal_process.c delete mode 100644 src/cen_anal_process.c delete mode 100644 src/corr_test.c delete mode 100644 src/fracture.c create mode 100644 src/fracture.cpp delete mode 100644 src/fracture.h delete mode 100644 src/long_anal_process.c create mode 100644 src/measurements.cpp create mode 100644 src/measurements.hpp create mode 160000 src/randutils (limited to 'src') diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..57a817b --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,7 @@ + +add_library(measurements measurements.cpp) +target_link_libraries(measurements frac) + +add_executable(fracture fracture.cpp) +target_link_libraries(fracture frac measurements fftw3 cholmod) + diff --git a/src/anal_process.c b/src/anal_process.c deleted file mode 100644 index de27571..0000000 --- a/src/anal_process.c +++ /dev/null @@ -1,135 +0,0 @@ - -#include "fracture.h" -#include -#include -#include -#include -#include - -void mom(uint_t len, uint_t n, uint32_t *data, double result[2]) { - uint_t total = 0; - double mom = 0; - double mom_err = 0; - - for (uint_t i = 0; i < len; i++) { - uint32_t datai = data[i]; - double in = pow(i, n); - - total += datai; - mom += in * datai; - mom_err += gsl_pow_2(in) * datai; - } - - double momf = mom / total; - double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total); - - result[0] = momf; - result[1] = momf_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)); - double *vals_c1 = (double *)malloc(nc * sizeof(double)); - double *errors_c1 = (double *)malloc(nc * sizeof(double)); - double *vals_c2 = (double *)malloc(nc * sizeof(double)); - double *errors_c2 = (double *)malloc(nc * sizeof(double)); - double *vals_c3 = (double *)malloc(nc * sizeof(double)); - double *errors_c3 = (double *)malloc(nc * sizeof(double)); - double *vals_c4 = (double *)malloc(nc * sizeof(double)); - double *errors_c4 = (double *)malloc(nc * sizeof(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); - - uint_t dist_len = 4 * pow(Ls_c[i], 2); - uint32_t *dist = malloc(dist_len * sizeof(uint32_t)); - FILE *dist_file = fopen(fn, "rb"); - fread(dist, sizeof(uint32_t), dist_len, dist_file); - fclose(dist_file); - { - double mom1[2]; - mom(dist_len, 1, dist, mom1); - vals_c1[i] = mom1[0]; - errors_c1[i] = mom1[1]; - - double mom2[2]; - mom(dist_len, 2, dist, mom2); - vals_c2[i] = mom2[0]; - errors_c2[i] = mom2[1]; - - double mom3[2]; - mom(dist_len, 3, dist, mom3); - vals_c3[i] = mom3[0]; - errors_c3[i] = mom3[1]; - - double mom4[2]; - mom(dist_len, 4, dist, 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 %g %g %g %g %g %g %g %g\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); - - return 0; -} diff --git a/src/big_anal_process.c b/src/big_anal_process.c deleted file mode 100644 index 8c1d8ba..0000000 --- a/src/big_anal_process.c +++ /dev/null @@ -1,158 +0,0 @@ - -#include "fracture.h" -#include -#include -#include -#include -#include -#include - -void get_mean(uint_t len, double *data, double result[2]) { - double total = 0; - - for (uint_t i = 0; i < len; i++) { - total += data[i]; - } - - double mean = total / len; - double total_sq_diff = 0; - - for (uint_t i = 0; i < len; i++) { - total_sq_diff += pow(mean - data[i], 2); - } - - double mean_err = sqrt(total_sq_diff) / len; - - result[0] = mean; - result[1] = mean_err; -} - -void get_mom(uint_t len, uint_t n, double *data, double mean[2], - double result[2]) { - double total_n_diff = 0; - double total_n2_diff = 0; - - for (uint_t i = 0; i < len; i++) { - total_n_diff += pow(fabs(mean[0] - data[i]), n); - total_n2_diff += pow(fabs(mean[0] - data[i]), 2 * n); - } - - double mom = total_n_diff / len; - 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)); - double *vals_c1 = (double *)malloc(nc * sizeof(double)); - double *errors_c1 = (double *)malloc(nc * sizeof(double)); - double *vals_c2 = (double *)malloc(nc * sizeof(double)); - double *errors_c2 = (double *)malloc(nc * sizeof(double)); - double *vals_c3 = (double *)malloc(nc * sizeof(double)); - double *errors_c3 = (double *)malloc(nc * sizeof(double)); - double *vals_c4 = (double *)malloc(nc * sizeof(double)); - double *errors_c4 = (double *)malloc(nc * sizeof(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(double); - - double *dist = malloc(dist_len * sizeof(double)); - FILE *dist_file = fopen(fn, "rb"); - fread(dist, sizeof(double), dist_len, dist_file); - fclose(dist_file); - { - double mom1[2]; - get_mean(dist_len, dist, mom1); - vals_c1[i] = mom1[0]; - errors_c1[i] = mom1[1]; - - double mom2[2]; - get_mom(dist_len, 2, dist, mom1, mom2); - vals_c2[i] = mom2[0]; - errors_c2[i] = mom2[1]; - - double mom3[2]; - get_mom(dist_len, 3, dist, mom1, mom3); - vals_c3[i] = mom3[0]; - errors_c3[i] = mom3[1]; - - 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 %g %g %g %g %g %g %g %g\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/cen_anal_process.c b/src/cen_anal_process.c deleted file mode 100644 index ee2b029..0000000 --- a/src/cen_anal_process.c +++ /dev/null @@ -1,157 +0,0 @@ - -#include "fracture.h" -#include -#include -#include -#include -#include - -void get_mean(uint_t len, uint32_t *data, double result[2]) { - uint_t total = 0; - double mean = 0; - double mean_err = 0; - - for (uint_t i = 0; i < len; i++) { - uint32_t datai = data[i]; - - total += datai; - mean += i * datai; - mean_err += gsl_pow_2(i) * datai; - } - - double meanf = mean / total; - double meanf_err = meanf * sqrt(mean_err / gsl_pow_2(mean) + 1 / total); - - result[0] = meanf; - result[1] = meanf_err; -} - -void get_mom(uint_t len, uint_t n, uint32_t *data, double mean[2], - double result[2]) { - uint_t total = 0; - double mom = 0; - double mom_err = 0; - - for (uint_t i = 0; i < len; i++) { - uint32_t datai = data[i]; - double in = pow(fabs(((double)i) - mean[0]), n); - - total += datai; - mom += in * datai; - mom_err += gsl_pow_2(in) * - datai; // + gsl_pow_2(n * mean[1] / (((double)i) - mean[0]))); - } - - double momf = mom / total; - double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total); - - result[0] = momf; - result[1] = momf_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)); - double *vals_c1 = (double *)malloc(nc * sizeof(double)); - double *errors_c1 = (double *)malloc(nc * sizeof(double)); - double *vals_c2 = (double *)malloc(nc * sizeof(double)); - double *errors_c2 = (double *)malloc(nc * sizeof(double)); - double *vals_c3 = (double *)malloc(nc * sizeof(double)); - double *errors_c3 = (double *)malloc(nc * sizeof(double)); - double *vals_c4 = (double *)malloc(nc * sizeof(double)); - double *errors_c4 = (double *)malloc(nc * sizeof(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); - - uint_t dist_len = 4 * pow(Ls_c[i], 2); - uint32_t *dist = malloc(dist_len * sizeof(uint32_t)); - FILE *dist_file = fopen(fn, "rb"); - fread(dist, sizeof(uint32_t), dist_len, dist_file); - fclose(dist_file); - { - double mom1[2]; - get_mean(dist_len, dist, mom1); - vals_c1[i] = mom1[0]; - errors_c1[i] = mom1[1]; - - double mom2[2]; - get_mom(dist_len, 2, dist, mom1, mom2); - vals_c2[i] = mom2[0]; - errors_c2[i] = mom2[1]; - - double mom3[2]; - get_mom(dist_len, 3, dist, mom1, mom3); - vals_c3[i] = mom3[0]; - errors_c3[i] = mom3[1]; - - 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 %g %g %g %g %g %g %g %g\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); - - return 0; -} diff --git a/src/corr_test.c b/src/corr_test.c deleted file mode 100644 index cb00212..0000000 --- a/src/corr_test.c +++ /dev/null @@ -1,27 +0,0 @@ - -#include "fracture.h" - -int main() { - cholmod_common c; - CHOL_F(start)(&c); - - unsigned int width = 64; - - graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 128, false); - net_t *instance = net_create(network, 1e-14, 3, 0, true, &c); - net_fracture(instance, &c, 1e-10); - - double *corr = get_corr(instance, NULL, &c); - - for (int i = 0; i < 2 * width; i++) { - printf("%g ", corr[i]); - } - printf("\n"); - - net_free(instance, &c); - graph_free(network); - - CHOL_F(finish)(&c); - - return 0; -} diff --git a/src/fracture.c b/src/fracture.c deleted file mode 100644 index 5b44238..0000000 --- a/src/fracture.c +++ /dev/null @@ -1,1068 +0,0 @@ - -#include - -#include "fracture.h" - -int main(int argc, char *argv[]) { - int opt; - - fftw_set_timelimit(1); - - // 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 save_data, save_cluster_dist, use_voltage_boundaries, use_dual, - save_network, save_crit_stress, save_energy, save_conductivity, - save_damage, save_threshold, save_current_load, save_correlations; - bound_t boundary; - lattice_t lattice; - - // assume filenames will be less than 100 characters - - filename_len = 100; - - // set default values - - N = 100; - L = 16; - crack_len = 0.; - beta = .3; - inf = 1e10; - cutoff = 1e-9; - boundary = FREE_BOUND; - lattice = VORONOI_LATTICE; - save_data = false; - save_cluster_dist = false; - use_voltage_boundaries = false; - use_dual = false; - save_network = false; - save_crit_stress = false; - save_damage = false; - save_conductivity = false; - save_energy = false; - save_threshold = false; - save_current_load = false; - save_correlations = false; - - uint8_t bound_i; - char boundc2 = 'f'; - uint8_t lattice_i; - char lattice_c = 'v'; - char dual_c = 'o'; - - // get commandline options - - while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TEz")) != -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 'q': - lattice_i = atoi(optarg); - switch (lattice_i) { - case 0: - lattice = VORONOI_LATTICE; - lattice_c = 'v'; - break; - case 1: - lattice = DIAGONAL_LATTICE; - lattice_c = 'd'; - break; - case 2: - lattice = VORONOI_HYPERUNIFORM_LATTICE; - lattice_c = 'h'; - break; - case 3: - lattice = TRIANGULAR_LATTICE; - lattice_c = 't'; - break; - case 4: - lattice = SQUARE_LATTICE; - lattice_c = 's'; - break; - default: - printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 " - "(DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n"); - exit(EXIT_FAILURE); - } - break; - case 'd': - save_damage = true; - break; - case 'V': - use_voltage_boundaries = true; - break; - case 'D': - use_dual = true; - dual_c = 'd'; - break; - case 'c': - save_cluster_dist = true; - break; - case 'o': - save_data = true; - break; - case 'N': - save_network = true; - break; - case 's': - save_crit_stress = true; - break; - case 'r': - save_conductivity = true; - break; - case 'E': - save_energy = true; - break; - case 'T': - save_threshold = true; - break; - case 'C': - save_current_load = true; - break; - case 'z': - save_correlations = true; - break; - default: /* '?' */ - exit(EXIT_FAILURE); - } - } - - char boundc; - if (use_voltage_boundaries) { - boundc = 'v'; - } else { - boundc = 'c'; - } - - double *dd_correlations; // damage-damage - double *dc_correlations; // damage-crack - double *db_correlations; // damage-backbone - double *ds_correlations; // damage-stress - double *dA_correlations; // damage-avalanche - double *cc_correlations; // crack-crack - double *cb_correlations; // crack-backbone - double *cs_correlations; // crack-stress - double *cA_correlations; // crack-avalanche - double *bb_correlations; // backbone-backbone - double *bs_correlations; // backbone-stress - double *bA_correlations; // backbone-avalanche - double *ss_correlations; // stress-stress - double *AA_correlations; // avalanche-avalanche - double *DD_correlations; // after-crack damage-damage - double *DC_correlations; // after-crack damage-crack - double *DB_correlations; // after-crack damage-backbone - double *CC_correlations; // after-crack crack-crack - double *CB_correlations; // after-crack crack-backbone - double *BB_correlations; // after-crack backbone-backbone - double *fftw_forward_in; - fftw_complex *fftw_forward_out; - fftw_complex *fftw_reverse_in; - double *fftw_reverse_out; - fftw_plan forward_plan; - fftw_plan reverse_plan; - uint64_t N_averaged = 0; - double mean_D = 0; - double mean_A = 0; - double mean_B = 0; - double mean_C = 0; - double mean_d = 0; - double mean_c = 0; - double mean_b = 0; - double mean_s = 0; - char *correlations_filename; - if (save_correlations) { - assert(lattice == DIAGONAL_LATTICE); - dd_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - dc_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - db_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - ds_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - dA_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - cc_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - cb_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - cs_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - cA_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - bb_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - bs_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - bA_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - ss_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - AA_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - DD_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - DC_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - DB_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - CC_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - CB_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - BB_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - fftw_forward_in = (double *)fftw_malloc(pow(L, 2) * sizeof(double)); - fftw_forward_out = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex)); - fftw_reverse_in = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex)); - fftw_reverse_out = (double *)fftw_malloc(pow(L, 2) * sizeof(double)); - forward_plan = fftw_plan_dft_r2c_2d(L, L, fftw_forward_in, fftw_forward_out, 0); - reverse_plan = fftw_plan_dft_c2r_2d(L, L, fftw_reverse_in, fftw_reverse_out, 0); - - correlations_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(correlations_filename, filename_len, "corr_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - - FILE *correlations_out = fopen(correlations_filename, "rb"); - - if (correlations_out != NULL) { - fread(&N_averaged, sizeof(uint64_t), 1, correlations_out); - fread(&mean_d, sizeof(double), 1, correlations_out); - fread(&mean_c, sizeof(double), 1, correlations_out); - fread(&mean_b, sizeof(double), 1, correlations_out); - fread(&mean_s, sizeof(double), 1, correlations_out); - fread(&mean_A, sizeof(double), 1, correlations_out); - fread(&mean_D, sizeof(double), 1, correlations_out); - fread(&mean_C, sizeof(double), 1, correlations_out); - fread(&mean_B, sizeof(double), 1, correlations_out); - fread(dd_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(dc_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(db_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(ds_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(dA_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(cc_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(cb_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(cs_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(cA_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(bb_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(bs_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(bA_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(ss_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(AA_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(DD_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(DC_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(DB_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(CC_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(CB_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(BB_correlations, sizeof(double), pow(L, 2), correlations_out); - fclose(correlations_out); - } - } - - FILE *data_out; - - if (save_data) { - char *data_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(data_filename, filename_len, "data_%c_%c_%c_%c_%u_%g_%g.txt", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - data_out = fopen(data_filename, "a"); - free(data_filename); - } - - uint_t max_verts, max_edges; - - // these are very liberal estimates - max_verts = 4 * pow(L, 2); - max_edges = 4 * pow(L, 2); - - if (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 *)calloc(max_verts, sizeof(uint32_t)); - avalanche_size_dist = (uint32_t *)calloc(max_edges, 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_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - snprintf(a_filename, filename_len, "avln_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - - 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), max_verts, cluster_out); - fclose(cluster_out); - } - if (avalanche_out != NULL) { - fread(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out); - fclose(avalanche_out); - } - } - - long double *crit_stress; - if (save_crit_stress) { - crit_stress = (long double *)malloc(N * sizeof(long 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 *)calloc(max_edges, sizeof(uint32_t)); - - d_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(d_filename, filename_len, "damg_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - - FILE *damage_out = fopen(d_filename, "rb"); - - if (damage_out != NULL) { - fread(damage, sizeof(uint32_t), max_edges, damage_out); - fclose(damage_out); - } - } - - long double *energy; - if (save_energy) { - energy = (long double *)malloc(N * sizeof(long double)); - } - - long double *thresholds; - if (save_threshold) { - thresholds = (long double *)malloc(N * sizeof(long double)); - } - - long double *loads; - if (save_current_load) { - loads = (long double *)malloc(N * sizeof(long 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); - - graph_t *g = graph_create(lattice, boundary, L, use_dual); - net_t *net = - net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c); - net_t *tmp_net = net_copy(net, &c); - data_t *data = net_fracture(tmp_net, &c, cutoff); - - uint_t max_pos = 0; - long double max_val = 0; - - double cond0; - { - double *tmp_voltages = net_voltages(net, &c); - cond0 = net_conductivity(net, tmp_voltages); - free(tmp_voltages); - } - - for (uint_t j = 0; j < data->num_broken; j++) { - long double val = data->extern_field[j]; - - if (val > max_val) { - max_pos = j; - max_val = val; - } - } - - uint_t av_size = 0; - long double cur_val = 0; - - for (uint_t j = 0; j < max_pos; j++) { - uint_t next_broken = data->break_list[j]; - - break_edge(net, next_broken, &c); - - long double val = 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_correlations) { - uint32_t damage1 = 0; - for (uint32_t j = 0; j < g->ne; j++) { - if (tmp_net->fuses[j]) { - fftw_forward_in[j] = 1.0; - damage1++; - } else { - fftw_forward_in[j] = 0.0; - } - } - - fftw_execute(forward_plan); - fftw_complex *D_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex)); - memcpy(D_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - dll_t *cycle = find_cycles(g, tmp_net->fuses); - components_t *comp = get_clusters(tmp_net); - - uint32_t in_crack1 = 0; - for (uint32_t j = 0; j < g->ne; j++) { - if (tmp_net->fuses[j] && comp->labels[g->dev[2 * cycle->e]] == comp->labels[g->dev[2 * j]]) { - fftw_forward_in[j] = 1.0; - in_crack1++; - } else { - fftw_forward_in[j] = 0.0; - } - } - - fftw_execute(forward_plan); - fftw_complex *C_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex)); - memcpy(C_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - for (uint32_t j = 0; j < g->ne; j++) { - fftw_forward_in[j] = 0.0; - } - - uint32_t in_backbone1 = 0; - dll_t *tmp_cycle = cycle; - while (tmp_cycle != NULL) { - fftw_forward_in[tmp_cycle->e] = 1.0; - in_backbone1++; - tmp_cycle = tmp_cycle->right; - } - - fftw_execute(forward_plan); - fftw_complex *B_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex)); - memcpy(B_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - uint32_t in_avalanche = 0; - for (uint32_t j = 0; j < g->ne; j++) { - if (tmp_net->fuses[j] && !(net->fuses[j])) { - fftw_forward_in[j] = 1.0; - } else { - fftw_forward_in[j] = 0.0; - } - } - - fftw_execute(forward_plan); - fftw_complex *A_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex)); - memcpy(A_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - uint32_t damage2 = 0; - for (uint32_t j = 0; j < g->ne; j++) { - if (net->fuses[j]) { - fftw_forward_in[j] = 1.0; - damage2++; - } else { - fftw_forward_in[j] = 0.0; - } - } - - fftw_execute(forward_plan); - fftw_complex *d_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex)); - memcpy(d_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - uint32_t in_crack2 = 0; - for (uint32_t j = 0; j < g->ne; j++) { - if (net->fuses[j] && comp->labels[2 * g->dev[cycle->e]] == comp->labels[g->dev[2 * j]]) { - fftw_forward_in[j] = 1.0; - in_crack2++; - } else { - fftw_forward_in[j] = 0.0; - } - } - - fftw_execute(forward_plan); - fftw_complex *c_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex)); - memcpy(c_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - for (uint32_t j = 0; j < g->ne; j++) { - fftw_forward_in[j] = 0.0; - } - - uint32_t in_backbone2 = 0; - tmp_cycle = cycle; - while (tmp_cycle != NULL) { - if (net->fuses[tmp_cycle->e]) { - fftw_forward_in[tmp_cycle->e] = 1.0; - in_backbone2++; - } - tmp_cycle = tmp_cycle->right; - } - - fftw_execute(forward_plan); - fftw_complex *b_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex)); - memcpy(b_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - graph_components_free(comp); - while (cycle != NULL) { - dll_t *old = cycle; - cycle = cycle->right; - free(old); - } - - double *tmp_voltage = net_voltages(net, &c); - double *tmp_current = net_currents(net, tmp_voltage, &c); - free(tmp_voltage); - - double t_stress = 0; - for (uint32_t j = 0; j < g->ne; j++) { - fftw_forward_in[j] = tmp_current[j]; - t_stress += tmp_current[j]; - } - - free(tmp_current); - - fftw_execute(forward_plan); - fftw_complex *s_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex)); - memcpy(s_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - mean_D = (double)damage1 / (1.0 + N_averaged) + (double)N_averaged * mean_D / (N_averaged + 1.0); - mean_C = (double)in_crack1 / (1.0 + N_averaged) + (double)N_averaged * mean_C / (N_averaged + 1.0); - mean_B = (double)in_backbone1 / (1.0 + N_averaged) + (double)N_averaged * mean_B / (N_averaged + 1.0); - mean_A = (double)in_avalanche / (1.0 + N_averaged) + (double)N_averaged * mean_A / (N_averaged + 1.0); - mean_d = (double)damage2 / (1.0 + N_averaged) + (double)N_averaged * mean_d / (N_averaged + 1.0); - mean_c = (double)in_crack2 / (1.0 + N_averaged) + (double)N_averaged * mean_c / (N_averaged + 1.0); - mean_b = (double)in_backbone2 / (1.0 + N_averaged) + (double)N_averaged * mean_b / (N_averaged + 1.0); - mean_s = (double)t_stress / (1.0 + N_averaged) + (double)N_averaged * mean_s / (N_averaged + 1.0); - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = D_transform[j][0] * D_transform[j][0] + D_transform[j][1] * D_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - DD_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DD_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = D_transform[j][0] * C_transform[j][0] + D_transform[j][1] * C_transform[j][1]; - fftw_reverse_in[j][1] = D_transform[j][0] * C_transform[j][1] - D_transform[j][1] * C_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - DC_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DC_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = D_transform[j][0] * B_transform[j][0] + D_transform[j][1] * B_transform[j][1]; - fftw_reverse_in[j][1] = D_transform[j][0] * B_transform[j][1] - D_transform[j][1] * B_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - DB_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DB_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = C_transform[j][0] * C_transform[j][0] + C_transform[j][1] * C_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - CC_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * CC_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = C_transform[j][0] * B_transform[j][0] + C_transform[j][1] * B_transform[j][1]; - fftw_reverse_in[j][1] = C_transform[j][0] * B_transform[j][1] - C_transform[j][1] * B_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - CB_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * CB_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = B_transform[j][0] * B_transform[j][0] + B_transform[j][1] * B_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - BB_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * BB_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = A_transform[j][0] * A_transform[j][0] + A_transform[j][1] * A_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - AA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * AA_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = d_transform[j][0] * d_transform[j][0] + d_transform[j][1] * d_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - dd_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dd_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = d_transform[j][0] * c_transform[j][0] + d_transform[j][1] * c_transform[j][1]; - fftw_reverse_in[j][1] = d_transform[j][0] * c_transform[j][1] - d_transform[j][1] * c_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - dc_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dc_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = d_transform[j][0] * s_transform[j][0] + d_transform[j][1] * s_transform[j][1]; - fftw_reverse_in[j][1] = d_transform[j][0] * s_transform[j][1] - d_transform[j][1] * s_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - ds_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * ds_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = d_transform[j][0] * b_transform[j][0] + d_transform[j][1] * b_transform[j][1]; - fftw_reverse_in[j][1] = d_transform[j][0] * b_transform[j][1] - d_transform[j][1] * b_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - db_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * db_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = d_transform[j][0] * A_transform[j][0] + d_transform[j][1] * A_transform[j][1]; - fftw_reverse_in[j][1] = d_transform[j][0] * A_transform[j][1] - d_transform[j][1] * A_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - dA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dA_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = c_transform[j][0] * c_transform[j][0] + c_transform[j][1] * c_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - cc_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cc_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = c_transform[j][0] * s_transform[j][0] + c_transform[j][1] * s_transform[j][1]; - fftw_reverse_in[j][1] = c_transform[j][0] * s_transform[j][1] - c_transform[j][1] * s_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - cs_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cs_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = c_transform[j][0] * b_transform[j][0] + c_transform[j][1] * b_transform[j][1]; - fftw_reverse_in[j][1] = c_transform[j][0] * b_transform[j][1] - c_transform[j][1] * b_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - cb_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cb_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = c_transform[j][0] * A_transform[j][0] + c_transform[j][1] * A_transform[j][1]; - fftw_reverse_in[j][1] = c_transform[j][0] * A_transform[j][1] - c_transform[j][1] * A_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - cA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cA_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = b_transform[j][0] * b_transform[j][0] + b_transform[j][1] * b_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - bb_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * bb_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = b_transform[j][0] * s_transform[j][0] + s_transform[j][1] * s_transform[j][1]; - fftw_reverse_in[j][1] = b_transform[j][0] * s_transform[j][1] - s_transform[j][1] * s_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - bs_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * bs_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = A_transform[j][0] * b_transform[j][0] + A_transform[j][1] * b_transform[j][1]; - fftw_reverse_in[j][1] = A_transform[j][0] * b_transform[j][1] - A_transform[j][1] * b_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - bA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * bA_correlations[j] / (N_averaged + 1.0); - } - - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = s_transform[j][0] * s_transform[j][0] + s_transform[j][1] * s_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - ss_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * ss_correlations[j] / (N_averaged + 1.0); - } - - free(D_transform); - free(C_transform); - free(A_transform); - free(B_transform); - free(b_transform); - free(d_transform); - free(c_transform); - free(s_transform); - N_averaged++; - } - - net_free(tmp_net, &c); - - if (save_crit_stress) { - crit_stress[i] = data->extern_field[max_pos]; - } - - if (save_conductivity) { - if (max_pos > 0) { - conductivity[i] = data->conductivity[max_pos - 1]; - } else { - conductivity[i] = cond0; - } - } - - if (save_damage) { - uint_t would_break = 0; - double *tmp_voltage = net_voltages(net, &c); - double *tmp_current = net_currents(net, tmp_voltage, &c); - free(tmp_voltage); - for (uint_t j = 0; j < g->ne; j++) { - bool broken = net->fuses[j]; - bool under_thres = - net->thres[j] < net->thres[data->break_list[max_pos]]; - bool zero_field = fabs(tmp_current[j]) < cutoff; - if (!broken && under_thres && zero_field) { - break_edge(net, j, &c); - } - } - damage[net->num_broken]++; - free(tmp_current); - } - - if (save_energy) { - long double tmp_energy = 0; - if (max_pos > 0) { - long double sigma1 = data->extern_field[0]; - double cond1 = cond0; - for (uint_t j = 0; j < max_pos - 1; j++) { - 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; - sigma1 = sigma2; - cond1 = cond2; - } - } - } - energy[i] = tmp_energy; - } - - if (save_threshold) { - thresholds[i] = net->thres[data->break_list[max_pos]]; - } - - if (save_current_load) { - loads[i] = - data->extern_field[max_pos] / net->thres[data->break_list[max_pos]]; - } - - if (save_data) { - for (uint_t j = 0; j < data->num_broken; j++) { - fprintf(data_out, "%u %Lg %g ", data->break_list[j], - data->extern_field[j], data->conductivity[j]); - } - fprintf(data_out, "\n"); - } - - data_free(data); - if (save_network) { - FILE *net_out = fopen("network.txt", "w"); - for (uint_t j = 0; j < g->nv; j++) { - fprintf(net_out, "%f %f ", g->vx[2 * j], g->vx[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (uint_t j = 0; j < g->ne; j++) { - fprintf(net_out, "%u %u ", g->ev[2 * j], g->ev[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (uint_t j = 0; j < g->dnv; j++) { - fprintf(net_out, "%f %f ", g->dvx[2 * j], g->dvx[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (uint_t j = 0; j < g->ne; j++) { - fprintf(net_out, "%u %u ", g->dev[2 * j], g->dev[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (uint_t j = 0; j < g->ne; j++) { - fprintf(net_out, "%d ", net->fuses[j]); - } - fclose(net_out); - } - - if (save_cluster_dist) { - uint_t *tmp_cluster_dist = get_cluster_dist(net); - for (uint_t j = 0; j < g->dnv; j++) { - cluster_size_dist[j] += tmp_cluster_dist[j]; - } - free(tmp_cluster_dist); - } - - net_free(net, &c); - graph_free(g); - } - - 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), max_verts, cluster_out); - fwrite(avalanche_size_dist, sizeof(uint32_t), max_edges, 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_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - FILE *cond_file = fopen(cond_filename, "ab"); - fwrite(conductivity, sizeof(double), N, cond_file); - fclose(cond_file); - free(cond_filename); - free(conductivity); - } - - 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(long double), N, tough_file); - fclose(tough_file); - free(tough_filename); - free(energy); - } - - if (save_threshold) { - 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(long double), N, thres_file); - fclose(thres_file); - free(thres_filename); - free(thresholds); - } - - if (save_current_load) { - char *load_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(load_filename, filename_len, "load_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - FILE *load_file = fopen(load_filename, "ab"); - fwrite(loads, sizeof(long double), N, load_file); - fclose(load_file); - free(load_filename); - free(loads); - } - - if (save_damage) { - FILE *hdam_file = fopen(d_filename, "wb"); - fwrite(damage, sizeof(uint32_t), max_edges, hdam_file); - fclose(hdam_file); - free(d_filename); - free(damage); - } - - if (save_data) { - fclose(data_out); - } - - if (save_crit_stress) { - 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(long double), N, str_file); - fclose(str_file); - free(str_filename); - free(crit_stress); - } - - if (save_correlations) { - FILE *correlations_out = fopen(correlations_filename, "wb"); - fwrite(&N_averaged, sizeof(uint64_t), 1, correlations_out); - fwrite(&mean_d, sizeof(double), 1, correlations_out); - fwrite(&mean_c, sizeof(double), 1, correlations_out); - fwrite(&mean_b, sizeof(double), 1, correlations_out); - fwrite(&mean_s, sizeof(double), 1, correlations_out); - fwrite(&mean_A, sizeof(double), 1, correlations_out); - fwrite(&mean_D, sizeof(double), 1, correlations_out); - fwrite(&mean_C, sizeof(double), 1, correlations_out); - fwrite(&mean_B, sizeof(double), 1, correlations_out); - fwrite(dd_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(dc_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(db_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(ds_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(dA_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(cc_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(cb_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(cs_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(cA_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(bb_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(bs_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(bA_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(ss_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(AA_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(DD_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(DC_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(DB_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(CC_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(CB_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(BB_correlations, sizeof(double), pow(L, 2), correlations_out); - fclose(correlations_out); - - free(dd_correlations); - free(dc_correlations); - free(db_correlations); - free(ds_correlations); - free(dA_correlations); - free(cc_correlations); - free(cb_correlations); - free(cs_correlations); - free(cA_correlations); - free(bb_correlations); - free(bs_correlations); - free(bA_correlations); - free(ss_correlations); - free(AA_correlations); - free(DD_correlations); - free(DC_correlations); - free(DB_correlations); - free(CC_correlations); - free(CB_correlations); - free(BB_correlations); - fftw_free(fftw_forward_in); - fftw_free(fftw_forward_out); - fftw_free(fftw_reverse_in); - fftw_free(fftw_reverse_out); - fftw_destroy_plan(forward_plan); - fftw_destroy_plan(reverse_plan); - free(correlations_filename); - } - - fftw_cleanup(); - - CHOL_F(finish)(&c); - - return 0; -} diff --git a/src/fracture.cpp b/src/fracture.cpp new file mode 100644 index 0000000..77af253 --- /dev/null +++ b/src/fracture.cpp @@ -0,0 +1,59 @@ + +#include +#include + +#include + +#include "randutils/randutils.hpp" + +#include +#include +#include +#include "measurements.hpp" + +int main(int argc, char* argv[]) { + int opt; + + unsigned int N = 1; + unsigned int L = 16; + double beta = 0.5; + + while ((opt = getopt(argc, argv, "N:L:b:")) != -1) { + switch (opt) { + case 'N': + N = (unsigned int)atof(optarg); + break; + case 'L': + L = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + default: + exit(1); + } + } + + cholmod_common c; + CHOL_F(start)(&c); + c.supernodal = CHOLMOD_SUPERNODAL; + + + graph G(L); + network base_network(G, &c); + ma meas(N, L, beta); + + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + for (unsigned int trial = 0; trial < N; trial++) { + network tmp_network(base_network); + tmp_network.set_thresholds(beta, rng); + tmp_network.fracture(meas); + } + + CHOL_F(finish)(&c); + + return 0; +} + diff --git a/src/fracture.h b/src/fracture.h deleted file mode 100644 index 5eb0a1d..0000000 --- a/src/fracture.h +++ /dev/null @@ -1,123 +0,0 @@ - -#pragma once - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -// these defs allow me to switch to long int cholmod in a sitch -#define int_t int -#define uint_t unsigned int -#define CINT_MAX INT_MAX -#define CHOL_F(x) cholmod_##x - -#define GSL_RAND_GEN gsl_rng_mt19937 - -typedef struct { - const graph_t *graph; - bool *fuses; - long double *thres; - double inf; - cholmod_dense *boundary_cond; - cholmod_factor *factor; - bool voltage_bound; - uint_t num_broken; - uint_t dim; - uint_t nep; - uint_t *evp; - cholmod_sparse *voltcurmat; -} net_t; - -typedef struct { - uint_t num_broken; - uint_t *break_list; - double *conductivity; - 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); - -cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, - bool symmetric, cholmod_common *c); - -cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c); - -int edge_to_verts(uint_t width, bool periodic, uint_t edge, bool index); - -int dual_edge_to_verts(uint_t width, bool periodic, uint_t edge, bool index); - -double dual_vert_to_coord(uint_t width, bool periodic, uint_t vert, bool index); - -void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, - cholmod_common *c); -void factor_update2(cholmod_factor *factor, uint_t v1, cholmod_common *c); - -void net_notch(net_t *net, double notch_len, cholmod_common *c); -data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff); -double *net_voltages(const net_t *net, cholmod_common *c); -double *net_currents(const net_t *net, const double *voltages, - cholmod_common *c); -double net_conductivity(const net_t *net, const double *voltages); - -void update_boundary(net_t *instance, const double *avg_field); - -FILE *get_file(const char *prefix, uint_t width, uint_t crack, double beta, - uint_t iter, uint_t num_iter, uint_t num, bool read); - -double update_beta(double beta, uint_t width, const double *stress, - const double *damage, double bound_total); - -cholmod_sparse *gen_voltcurmat(uint_t num_edges, uint_t num_verts, - uint_t *edges_to_verts, cholmod_common *c); - -net_t *net_copy(const net_t *net, cholmod_common *c); - -void net_free(net_t *instance, cholmod_common *c); - -net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, - bool vb, cholmod_common *c); - -bool break_edge(net_t *instance, uint_t edge, cholmod_common *c); - -components_t *get_clusters(net_t *instance); - -uint_t *get_cluster_dist(net_t *instance); - -void randfunc_flat(gsl_rng *r, double *x, double *y); -void randfunc_gaus(gsl_rng *r, double *x, double *y); - -double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c); - -double *bin_values(graph_t *network, uint_t width, double *values); - -cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, - cholmod_common *c); - -data_t *data_create(uint_t num_edges); -void data_free(data_t *data); -void data_update(data_t *data, uint_t last_broke, long double strength, - double conductivity); - -long double rand_dist_pow(const gsl_rng *r, double beta); - -bool is_in(uint_t len, uint_t *list, uint_t element); diff --git a/src/long_anal_process.c b/src/long_anal_process.c deleted file mode 100644 index d4ec4e0..0000000 --- a/src/long_anal_process.c +++ /dev/null @@ -1,158 +0,0 @@ - -#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/measurements.cpp b/src/measurements.cpp new file mode 100644 index 0000000..7d5d539 --- /dev/null +++ b/src/measurements.cpp @@ -0,0 +1,248 @@ + +#include "measurements.hpp" + +bool trivial(boost::detail::edge_desc_impl) { + return true; +} + +void update_distribution_data(std::string id, const std::vector& data, unsigned int N, unsigned int L, double beta) { + std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat"; + std::ifstream file(filename); + + uint64_t N_old = 0; + std::vector data_old(data.size(), 0); + + if (file.is_open()) { + file >> N_old; + for (unsigned int i = 0; i < data.size(); i++) { + uint64_t num; + file >> num; + data_old[i] = num; + } + + file.close(); + } + + std::ofstream file_out(filename); + + file_out <& data, unsigned int L, double beta) { + std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat"; + std::ifstream file(filename); + + uint64_t N_old = 0; + uint64_t tot_old = 0; + uint64_t tot_2_old = 0; + std::vector data_old(data.size(), 0); + std::vector data_old_2(data.size(), 0); + + if (file.is_open()) { + file >> N_old; + file >> tot_old; + file >> tot_2_old; + for (unsigned int i = 0; i < data.size(); i++) { + uint64_t num; + file >> num; + data_old[i] = num; + } + for (unsigned int i = 0; i < data.size(); i++) { + uint64_t num; + file >> num; + data_old_2[i] = num; + } + + file.close(); + } + + std::ofstream file_out(filename); + + file_out < data_transform(unsigned int L, const std::vector& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) { + for (unsigned int i = 0; i < pow(L, 2); i++) { + fftw_forward_in[i] = data[i]; + } + + fftw_execute(forward_plan); + + std::vector output(pow(L, 2)); + + for (unsigned int i = 0; i < pow(L, 2); i++) { + output[i][0] = fftw_forward_out[i][0]; + output[i][1] = fftw_forward_out[i][1]; + } + + return output; +} + +std::vector correlation(unsigned int L, const std::vector& tx1, const std::vector& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) { + for (unsigned int i = 0; i < pow(L, 2); i++) { + fftw_reverse_in[i][0] = tx1[i][0] * tx2[i][0] + tx1[i][1] * tx2[i][1]; + fftw_reverse_in[i][1] = tx1[i][0] * tx2[i][1] - tx1[i][1] * tx2[i][0]; + } + + fftw_execute(reverse_plan); + + std::vector output(pow(L / 2 + 1, 2)); + + for (unsigned int i = 0; i < pow(L / 2 + 1, 2); i++) { + output[i] = fftw_reverse_out[L * (i / (L / 2 + 1)) + i % (L / 2 + 1)] / pow(L, 2); + } + + return output; +} + +ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2)), bin_counts(log2(L) + 1, 0), N(N), beta(beta) { + ad.resize(pow(L, 2), 0); + cd.resize(pow(L, 2), 0); + + // FFTW setup for correlation functions + fftw_set_timelimit(1); + + fftw_forward_in = (double *)fftw_malloc(pow(L, 2) * sizeof(double)); + fftw_forward_out = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex)); + fftw_reverse_in = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex)); + fftw_reverse_out = (double *)fftw_malloc(pow(L, 2) * sizeof(double)); + + forward_plan = fftw_plan_dft_r2c_2d(L, L, fftw_forward_in, fftw_forward_out, 0); + reverse_plan = fftw_plan_dft_c2r_2d(L, L, fftw_reverse_in, fftw_reverse_out, 0); +} + +ma::~ma() { + // clean up FFTW objects + fftw_free(fftw_forward_in); + fftw_free(fftw_forward_out); + fftw_free(fftw_reverse_in); + fftw_free(fftw_reverse_out); + fftw_destroy_plan(forward_plan); + fftw_destroy_plan(reverse_plan); + fftw_cleanup(); + + update_distribution_data("na", ad, N, L, beta); + update_distribution_data("nc", cd, N, L, beta); + update_distribution_data("bc", bin_counts, N, L, beta); + +} + +void ma::pre_fracture(const network &) { + lv = 0; + as = 0; + avalanches.clear(); + boost::remove_edge_if(trivial, G); +} + +void ma::bond_broken(const network& net, const std::pair>& cur, unsigned int i) { + if (cur.first / fabs(cur.second[i]) * net.thresholds[i] > lv) { + ad[as]++; + as = 0; + lv = cur.first / fabs(cur.second[i]) * net.thresholds[i]; + avalanches.push_back({}); + } else { + as++; + avalanches.back().push_back(i); + } + + boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], G); +} + +void ma::post_fracture(network &n) { + std::vector component(boost::num_vertices(G)); + unsigned int num = boost::connected_components(G, &component[0]); + + std::vector comp_sizes(num, 0); + + for (unsigned int c : component) { + comp_sizes[c]++; + } + + unsigned int max_i = 0; + unsigned int max_size = 0; + + for (unsigned int i = 0; i < num; i++) { + if (comp_sizes[i] > max_size) { + max_i = i; + max_size = comp_sizes[i]; + } + } + + for (unsigned int be = 0; be < log2(L); be++) { + unsigned int bin = pow(2, be); + + for (unsigned int i = 0; i < pow(L / bin, 2); i++) { + bool in_bin = false; + for (unsigned int j = 0; j < pow(bin, 2); j++) { + unsigned int edge = L * (bin * ((i * bin) / L) + j / bin) + (i * bin) % L + j % bin; + if (!n.fuses[edge] && max_i == component[n.G.dual_edges[edge][0]]) { + in_bin = true; + break; + } + } + + if (in_bin) { + bin_counts[be]++; + } + } + } + + bin_counts[log2(L)]++; + + std::vector crack_damage(pow(L, 2), 0.0); + + double damage_tot = 0; + for (unsigned int i = 0; i < pow(L, 2); i++) { + if (!n.fuses[i] && max_i == component[n.G.dual_edges[i][0]]) { + damage_tot++; + crack_damage[i] = 1.0; + } + } + + std::vector t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out); + std::vector Ccc = correlation(L, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out); + + update_field_data("Ccc", damage_tot, Ccc, L, beta); + + for (auto e = avalanches.back().rbegin(); e != avalanches.back().rend(); e++) { + boost::remove_edge(n.G.dual_edges[*e][0], n.G.dual_edges[*e][1], G); + } + + num = boost::connected_components(G, &component[0]); + + for (unsigned int i = 0; i < num; i++) { + double size = 0; + std::fill(crack_damage.begin(), crack_damage.end(), 0.0); + + for (unsigned int j = 0; j < pow(L, 2); j++) { + if (component[n.G.edges[j][0]] == i && n.fuses[j]) { + size++; + crack_damage[j] = 1.0; + } + } + + if (size > 0) { + cd[size - 1]++; + t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out); + Ccc = correlation(L, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out); + + update_field_data("Cclcl", size, Ccc, L, beta); + } + } + +} + diff --git a/src/measurements.hpp b/src/measurements.hpp new file mode 100644 index 0000000..817481b --- /dev/null +++ b/src/measurements.hpp @@ -0,0 +1,54 @@ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include + +class ma : public hooks { + // need: + // - measurements for correlation functions (nice looking?) + // - interface for reading and writing from files (nice looking?) + // - interface for turning on and off specific measurements + // + private: + double *fftw_forward_in; + fftw_complex *fftw_forward_out; + fftw_complex *fftw_reverse_in; + double *fftw_reverse_out; + fftw_plan forward_plan; + fftw_plan reverse_plan; + boost::adjacency_list G; + + public: + unsigned int N; + unsigned int L; + unsigned int as; + double lv; + double beta; + + + std::vector ad; + std::vector cd; + std::vector bin_counts; + std::list> avalanches;; + + ma(unsigned int N, unsigned int L, double beta); + ~ma(); + + void pre_fracture(const network &) override; + void bond_broken(const network& net, const std::pair>& cur, unsigned int i) override; + void post_fracture(network &n) override; +}; + diff --git a/src/randutils b/src/randutils new file mode 160000 index 0000000..8486a61 --- /dev/null +++ b/src/randutils @@ -0,0 +1 @@ +Subproject commit 8486a610a954a8248c12485fb4cfc390a5f5f854 -- cgit v1.2.3-54-g00ecf