summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/data.c4
-rw-r--r--src/fracture.c106
-rw-r--r--src/fracture.h4
-rw-r--r--src/long_anal_process.c156
-rw-r--r--src/net_fracture.c4
5 files changed, 182 insertions, 92 deletions
diff --git a/src/data.c b/src/data.c
index b54236d..2047c44 100644
--- a/src/data.c
+++ b/src/data.c
@@ -10,7 +10,7 @@ data_t *data_create(uint_t ne) {
data->break_list = (uint_t *)malloc(ne * sizeof(uint_t));
assert(data->break_list != NULL);
- data->extern_field = (double *)malloc(ne * sizeof(double));
+ data->extern_field = (long double *)malloc(ne * sizeof(long double));
assert(data->extern_field != NULL);
data->conductivity = (double *)malloc(ne * sizeof(double));
@@ -26,7 +26,7 @@ void data_free(data_t *data) {
free(data);
}
-void data_update(data_t *data, uint_t last_broke, double strength, double conductivity) {
+void data_update(data_t *data, uint_t last_broke, long double strength, double conductivity) {
data->break_list[data->num_broken] = last_broke;
data->extern_field[data->num_broken] = strength;
data->conductivity[data->num_broken] = conductivity;
diff --git a/src/fracture.c b/src/fracture.c
index 7ee82fb..59bd424 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -10,8 +10,8 @@ int main(int argc, char *argv[]) {
uint32_t N;
uint_t L, refactor_every;
double beta, inf, cutoff, crack_len;
- bool crack_growth_crit, save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
- save_crit_stress, save_stress_field, save_voltage_field, save_toughness, save_energy, save_conductivity,
+ bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
+ save_crit_stress, save_stress_field, save_voltage_field, save_energy, save_conductivity,
save_damage, save_damage_field, save_threshold;
bound_t boundary;
lattice_t lattice;
@@ -33,7 +33,6 @@ int main(int argc, char *argv[]) {
cutoff = 1e-9;
boundary = FREE_BOUND;
lattice = VORONOI_LATTICE;
- crack_growth_crit = false;
save_data = false;
save_cluster_dist = false;
use_voltage_boundaries = false;
@@ -45,7 +44,6 @@ int main(int argc, char *argv[]) {
save_damage = false;
save_damage_field = false;
save_conductivity = false;
- save_toughness = false;
save_energy = true;
save_threshold = false;
@@ -59,7 +57,7 @@ int main(int argc, char *argv[]) {
// get commandline options
- while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:gTER:")) != -1) {
+ while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDSvel:TER:")) != -1) {
switch (opt) {
case 'n':
N = atoi(optarg);
@@ -119,9 +117,6 @@ int main(int argc, char *argv[]) {
exit(EXIT_FAILURE);
}
break;
- case 'g':
- crack_growth_crit = true;
- break;
case 'd':
save_damage = true;
break;
@@ -156,9 +151,6 @@ int main(int argc, char *argv[]) {
case 'r':
save_conductivity = true;
break;
- case 't':
- save_toughness = true;
- break;
case 'E':
save_energy = true;
break;
@@ -222,9 +214,9 @@ int main(int argc, char *argv[]) {
}
}
- double *crit_stress;
+ long double *crit_stress;
if (save_crit_stress) {
- crit_stress = (double *)malloc(N * sizeof(double));
+ crit_stress = (long double *)malloc(N * sizeof(long double));
}
double *stress_field;
@@ -268,19 +260,14 @@ int main(int argc, char *argv[]) {
}
}
- double *toughness;
- if (save_toughness) {
- toughness = (double *)malloc(N * sizeof(double));
- }
-
- double *energy;
+ long double *energy;
if (save_energy) {
- energy = (double *)malloc(N * sizeof(double));
+ energy = (long double *)malloc(N * sizeof(long double));
}
- double *thresholds;
+ long double *thresholds;
if (save_threshold) {
- thresholds = (double *)malloc(N * sizeof(double));
+ thresholds = (long double *)malloc(N * sizeof(long double));
}
@@ -311,7 +298,7 @@ int main(int argc, char *argv[]) {
net_free(tmp_net, &c);
uint_t max_pos = 0;
- double max_val = 0;
+ long double max_val = 0;
double cond0;
{
@@ -321,7 +308,7 @@ int main(int argc, char *argv[]) {
}
for (uint_t j = 0; j < data->num_broken; j++) {
- double val = data->extern_field[j];
+ long double val = data->extern_field[j];
if (val > max_val) {
max_pos = j;
@@ -330,41 +317,15 @@ int main(int argc, char *argv[]) {
}
uint_t av_size = 0;
- double cur_val = 0;
- uint_t notch_mark = 0;
-
- if (crack_growth_crit) {
- uint_t n_edge = UINT_MAX;
- for (uint_t j = 0; j < g->ne; j++) {
- // right now the only fuses broken are from a notch
- if (net->fuses[j]) {
- n_edge = j;
- break;
- }
- }
-
- if (n_edge < UINT_MAX) {
- uint_t n_vert = g->dev[2 * n_edge];
- notch_mark = net->dual_marks[n_vert];
- }
- }
+ long double cur_val = 0;
for (uint_t j = 0; j < max_pos; j++) {
uint_t next_broken = data->break_list[j];
- bool attached_to_notch = net->dual_marks[g->dev[2*next_broken]] == notch_mark;
- bool grew_forward = g->ex[2 * next_broken] > crack_len;
-
- if (crack_growth_crit && attached_to_notch && grew_forward) {
- printf("whoops\n\n");
- max_pos = j;
- break;
- }
-
bool refactor = ((j + 1) % refactor_every) == 0;
break_edge(net, next_broken, &c, refactor);
- double val = data->extern_field[j];
+ long double val = data->extern_field[j];
if (save_cluster_dist) {
if (val < cur_val) {
av_size++;
@@ -422,12 +383,12 @@ int main(int argc, char *argv[]) {
}
if (save_energy) {
- double tmp_energy = 0;
+ long double tmp_energy = 0;
if (max_pos > 0) {
- double sigma1 = data->extern_field[0];
+ long double sigma1 = data->extern_field[0];
double cond1 = cond0;
for (uint_t j = 0; j < max_pos - 1; j++) {
- double sigma2 = data->extern_field[j+1];
+ long double sigma2 = data->extern_field[j+1];
double cond2 = data->conductivity[j];
if (sigma2 > sigma1) {
tmp_energy += 0.5 * gsl_pow_2(sigma1) * (1 - cond2 / cond1) / cond1;
@@ -439,23 +400,6 @@ int main(int argc, char *argv[]) {
energy[i] = tmp_energy;
}
- if (save_toughness) {
- double tmp_toughness = 0;
- if (max_pos > 0) {
- double sigma1 = data->extern_field[0];
- double epsilon1 = sigma1 / data->conductivity[0];
- for (uint_t j = 0; j < max_pos - 1; j++) {
- double sigma2 = data->extern_field[j+1];
- double epsilon2 = sigma2 / data->conductivity[j+1];
- if (epsilon2 > epsilon1) {
- tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2;
- sigma1 = sigma2; epsilon1 = epsilon2;
- }
- }
- }
- toughness[i] = tmp_toughness;
- }
-
if (save_threshold) {
thresholds[i] = net->thres[data->break_list[max_pos]];
}
@@ -499,7 +443,7 @@ int main(int argc, char *argv[]) {
if (save_data) {
for (uint_t j = 0; j < data->num_broken; j++) {
- fprintf(data_out, "%u %g %g ", data->break_list[j],
+ fprintf(data_out, "%u %Lg %g ", data->break_list[j],
data->extern_field[j], data->conductivity[j]);
}
fprintf(data_out, "\n");
@@ -566,21 +510,11 @@ int main(int argc, char *argv[]) {
free(conductivity);
}
- if (save_toughness) {
- char *tough_filename = (char *)malloc(filename_len * sizeof(char));
- snprintf(tough_filename, filename_len, "tuff_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
- FILE *tough_file = fopen(tough_filename, "ab");
- fwrite(toughness, sizeof(double), N, tough_file);
- fclose(tough_file);
- free(tough_filename);
- free(toughness);
- }
-
if (save_energy) {
char *tough_filename = (char *)malloc(filename_len * sizeof(char));
snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
FILE *tough_file = fopen(tough_filename, "ab");
- fwrite(energy, sizeof(double), N, tough_file);
+ fwrite(energy, sizeof(long double), N, tough_file);
fclose(tough_file);
free(tough_filename);
free(energy);
@@ -590,7 +524,7 @@ int main(int argc, char *argv[]) {
char *thres_filename = (char *)malloc(filename_len * sizeof(char));
snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
FILE *thres_file = fopen(thres_filename, "ab");
- fwrite(thresholds, sizeof(double), N, thres_file);
+ fwrite(thresholds, sizeof(long double), N, thres_file);
fclose(thres_file);
free(thres_filename);
free(thresholds);
@@ -612,7 +546,7 @@ int main(int argc, char *argv[]) {
char *str_filename = (char *)malloc(filename_len * sizeof(char));
snprintf(str_filename, filename_len, "strs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
FILE *str_file = fopen(str_filename, "ab");
- fwrite(crit_stress, sizeof(double), N, str_file);
+ fwrite(crit_stress, sizeof(long double), N, str_file);
fclose(str_file);
free(str_filename);
free(crit_stress);
diff --git a/src/fracture.h b/src/fracture.h
index f122869..556e1e5 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -89,7 +89,7 @@ typedef struct {
uint_t num_broken;
uint_t *break_list;
double *conductivity;
- double *extern_field;
+ long double *extern_field;
} data_t;
intptr_t *run_voronoi(uint_t num_coords, double *coords, bool periodic, double xmin, double xmax, double ymin, double ymax);
@@ -173,7 +173,7 @@ cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_co
data_t *data_create(uint_t num_edges);
void data_free(data_t *data);
-void data_update(data_t *data, uint_t last_broke, double strength, double conductivity);
+void data_update(data_t *data, uint_t last_broke, long double strength, double conductivity);
graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c);
diff --git a/src/long_anal_process.c b/src/long_anal_process.c
new file mode 100644
index 0000000..ba29152
--- /dev/null
+++ b/src/long_anal_process.c
@@ -0,0 +1,156 @@
+
+#include "fracture.h"
+#include <gsl/gsl_sf_erf.h>
+#include <gsl/gsl_sf_laguerre.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.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/net_fracture.c b/src/net_fracture.c
index 3bfcd83..b5ab524 100644
--- a/src/net_fracture.c
+++ b/src/net_fracture.c
@@ -46,7 +46,7 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff, uint_t refact
uint_t last_broke = get_next_broken(net, currents, cutoff);
- double sim_current;
+ long double sim_current;
if (net->voltage_bound) {
sim_current = conductivity;
@@ -54,7 +54,7 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff, uint_t refact
sim_current = 1;
}
- data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] / currents[last_broke]), conductivity);
+ data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] / ((long double)currents[last_broke])), conductivity);
free(voltages);
free(currents);