summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/CMakeLists.txt7
-rw-r--r--src/anal_process.c135
-rw-r--r--src/big_anal_process.c158
-rw-r--r--src/cen_anal_process.c157
-rw-r--r--src/corr_test.c27
-rw-r--r--src/fracture.c1068
-rw-r--r--src/fracture.cpp59
-rw-r--r--src/fracture.h123
-rw-r--r--src/long_anal_process.c158
-rw-r--r--src/measurements.cpp248
-rw-r--r--src/measurements.hpp54
m---------src/randutils0
12 files changed, 368 insertions, 1826 deletions
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 <gsl/gsl_blas.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_sf_erf.h>
-#include <gsl/gsl_sf_laguerre.h>
-#include <gsl/gsl_vector.h>
-
-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 <gsl/gsl_blas.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_sf_erf.h>
-#include <gsl/gsl_sf_laguerre.h>
-#include <gsl/gsl_vector.h>
-#include <sys/stat.h>
-
-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 <gsl/gsl_blas.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_sf_erf.h>
-#include <gsl/gsl_sf_laguerre.h>
-#include <gsl/gsl_vector.h>
-
-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 <fftw3.h>
-
-#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 <random>
+#include <iostream>
+
+#include <cholmod.h>
+
+#include "randutils/randutils.hpp"
+
+#include <graph.hpp>
+#include <network.hpp>
+#include <hooks.hpp>
+#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 <assert.h>
-#include <cholmod.h>
-#include <float.h>
-#include <getopt.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_randist.h>
-#include <gsl/gsl_rng.h>
-#include <gsl/gsl_sf_erf.h>
-#include <gsl/gsl_sf_exp.h>
-#include <gsl/gsl_sf_log.h>
-#include <inttypes.h>
-#include <math.h>
-#include <stdbool.h>
-#include <stdint.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <sys/types.h>
-#include <unistd.h>
-
-#include <jst/graph.h>
-#include <jst/rand.h>
-
-// 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 <gsl/gsl_blas.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_sf_erf.h>
-#include <gsl/gsl_sf_laguerre.h>
-#include <gsl/gsl_vector.h>
-#include <sys/stat.h>
-
-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<boost::undirected_tag,unsigned long>) {
+ return true;
+}
+
+void update_distribution_data(std::string id, const std::vector<unsigned int>& 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<uint64_t> 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 <<std::fixed<< N_old + N << "\n";
+ for (unsigned int i = 0; i < data.size(); i++) {
+ file_out <<std::fixed<< data_old[i] + data[i] << " ";
+ }
+
+ file_out.close();
+}
+
+void update_field_data(std::string id, double tot, const std::vector<double>& 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<uint64_t> data_old(data.size(), 0);
+ std::vector<uint64_t> 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 <<std::fixed<< N_old + 1 << " " << (uint64_t)(tot_old + tot) << " " << (uint64_t)(tot_2_old + pow(tot, 2))<< "\n";
+ for (unsigned int i = 0; i < data.size(); i++) {
+ file_out << (uint64_t)(data_old[i] + data[i]) << " ";
+ }
+ file_out <<std::fixed<< "\n";
+ for (unsigned int i = 0; i < data.size(); i++) {
+ file_out << (uint64_t)(data_old_2[i] + pow(data[i], 2)) << " ";
+ }
+
+ file_out.close();
+}
+
+std::vector<fftw_complex> data_transform(unsigned int L, const std::vector<double>& 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<fftw_complex> 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<double> correlation(unsigned int L, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& 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<double> 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<double, std::vector<double>>& 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<unsigned int> component(boost::num_vertices(G));
+ unsigned int num = boost::connected_components(G, &component[0]);
+
+ std::vector<unsigned int> 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<double> 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<fftw_complex> t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out);
+ std::vector<double> 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 <vector>
+#include <algorithm>
+#include <cmath>
+#include <list>
+#include <fstream>
+#include <string>
+#include <cinttypes>
+#include <sstream>
+
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/connected_components.hpp>
+
+#include <fftw3.h>
+
+#include <network.hpp>
+#include <hooks.hpp>
+
+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 <boost::listS, boost::vecS, boost::undirectedS> G;
+
+ public:
+ unsigned int N;
+ unsigned int L;
+ unsigned int as;
+ double lv;
+ double beta;
+
+
+ std::vector<unsigned int> ad;
+ std::vector<unsigned int> cd;
+ std::vector<unsigned int> bin_counts;
+ std::list<std::list<unsigned int>> 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<double, std::vector<double>>& cur, unsigned int i) override;
+ void post_fracture(network &n) override;
+};
+
diff --git a/src/randutils b/src/randutils
new file mode 160000
+Subproject 8486a610a954a8248c12485fb4cfc390a5f5f85