summaryrefslogtreecommitdiff
path: root/src/current_scaling.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/current_scaling.c')
-rw-r--r--src/current_scaling.c300
1 files changed, 0 insertions, 300 deletions
diff --git a/src/current_scaling.c b/src/current_scaling.c
deleted file mode 100644
index ef36b0e..0000000
--- a/src/current_scaling.c
+++ /dev/null
@@ -1,300 +0,0 @@
-
-#include "fracture.h"
-
-int main(int argc, char *argv[]) {
- int opt;
-
- // defining variables to be (potentially) set by command line flags
- int iter = 1;
- int num = 100;
- int width = 16;
- double crack_len = 8;
- double beta = .3;
- double inf = 1e-8;
- double cutoff = 1e-8;
- bool beta_shift = false;
- bool supplied_bound = false;
- bool ash_beta = false;
- char *bound_file;
- bool voltage_bound = false;
- bool use_first = false;
- bool save_stress = false;
- bool save_bound = false;
- bool save_damage = false;
- bool save_strength = false;
-
- while ((opt = getopt(argc, argv, "n:w:b:l:i:Bf:aVFsSed")) != -1) {
- switch (opt) {
- case 'n':
- num = atoi(optarg);
- break;
- case 'w':
- width = atoi(optarg);
- break;
- case 'b':
- beta = atof(optarg);
- break;
- case 'l':
- crack_len = atof(optarg);
- break;
- case 'i':
- iter = atoi(optarg);
- break;
- case 'B':
- beta_shift = true;
- break;
- case 'a':
- ash_beta = true;
- break;
- case 'V':
- voltage_bound = true;
- break;
- case 'F':
- use_first = true;
- break;
- case 's':
- save_stress = true;
- break;
- case 'd':
- save_damage = true;
- break;
- case 'e':
- save_bound = true;
- break;
- case 'S':
- save_strength = true;
- break;
- case 'f':
- supplied_bound = true;
- bound_file = optarg;
- break;
- default: /* '?' */
- exit(EXIT_FAILURE);
- }
- }
-
- // 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 (voltage_bound) {
- (&c)->supernodal = CHOLMOD_SUPERNODAL;
- } else {
- (&c)->supernodal = CHOLMOD_SIMPLICIAL;
- }
-
- graph_t *network = ini_square_network(width, false, true, &c);
- net_t *perm_instance =
- create_instance(network, inf, voltage_bound, false, &c);
- gen_crack(perm_instance, crack_len, &c);
- finish_instance(perm_instance, &c);
-
- if (voltage_bound) {
- (&c)->supernodal = CHOLMOD_SIMPLICIAL;
- net_t *tmp_instance = create_instance(network, inf, false, false, &c);
- gen_crack(tmp_instance, crack_len, &c);
- finish_instance(tmp_instance, &c);
- double *voltage = get_voltage(tmp_instance, &c);
-
- for (int i = 0; i < network->num_bounds; i++) {
- for (int j = 0; j < network->bound_inds[i + 1] - network->bound_inds[i];
- j++) {
- ((double *)perm_instance->boundary_cond
- ->x)[network->bound_verts[network->bound_inds[i] + j]] =
- voltage[network->bound_verts[network->bound_inds[i] + j]];
- }
- }
- (&c)->supernodal = CHOLMOD_SUPERNODAL;
- }
-
- if (supplied_bound) {
- FILE *bound_f = fopen(bound_file, "r");
- for (int i = 0; i < network->nv; i++) {
- double tmp;
- fscanf(bound_f, "%lg ", &tmp);
- ((double *)perm_instance->boundary_cond->x)[i] = tmp;
- }
-
- ((double *)perm_instance->boundary_cond->x)[network->nv] = 0;
- ((double *)perm_instance->boundary_cond->x)[network->nv + 1] = 0;
- fclose(bound_f);
- }
-
- printf("\n");
- for (int DUMB2 = 0; DUMB2 < iter; DUMB2++) {
-
- double *strength;
- if (save_strength) {
- strength = (double *)malloc(num * sizeof(double));
- }
-
- double *damage;
- if (save_damage) {
- damage = (double *)calloc(network->ne, sizeof(double));
- }
- double *avg_current = (double *)calloc(network->ne, sizeof(double));
- unsigned int *num_current_skipped =
- (unsigned int *)calloc(network->ne, sizeof(unsigned int));
- double *avg_voltage = (double *)calloc(network->nv, sizeof(double));
- unsigned int *num_voltage_skipped =
- (unsigned int *)calloc(network->nv, sizeof(unsigned int));
-
- for (int DUMB = 0; DUMB < num; DUMB++) {
- printf("\033[F\033[JCURRENT_SCALING: ITERATION %0*d: %0*d / %d\n",
- (int)log10(iter) + 1, DUMB2 + 1, (int)log10(num) + 1, DUMB + 1,
- num);
-
- data_t *breaking_data = NULL;
- while (breaking_data == NULL) {
- double *fuse_thres = gen_fuse_thres(
- network->ne, network->ex, beta, beta_scaling_flat);
- net_t *instance = copy_instance(perm_instance, &c);
- breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
- free_instance(instance, &c);
- free(fuse_thres);
- }
-
- unsigned int min_pos = 0;
- double min_val = DBL_MAX;
-
- for (unsigned int j = 0; j < breaking_data->num_broken; j++) {
- double val = fabs(breaking_data->extern_field[j]);
- if (val < min_val) {
- min_pos = j;
- min_val = val;
- }
- }
-
- if (save_strength) {
- strength[DUMB] = fabs(breaking_data->extern_field[min_pos]);
- }
-
- net_t *tmp_instance = copy_instance(perm_instance, &c);
-
- unsigned int until = min_pos;
- if (use_first)
- until = 1;
- for (unsigned int i = 0; i < until; i++) {
- break_edge(tmp_instance, breaking_data->break_list[i], &c);
- if (save_damage) {
- damage[breaking_data->break_list[i]] += 1. / num;
- }
- }
-
- double *voltage = get_voltage(tmp_instance, &c);
- double *current = get_current(tmp_instance, &c);
-
- for (unsigned int i = 0; i < network->ne; i++) {
- avg_current[i] += current[i];
- if (current[i] == 0)
- num_current_skipped[i]++;
- }
-
- for (unsigned int i = 0; i < network->nv; i++) {
- if (tmp_instance->marks[i] == tmp_instance->marks[network->nv]) {
- avg_voltage[i] += voltage[i];
- } else {
- num_voltage_skipped[i]++;
- }
- }
-
- free(current);
- free(voltage);
- free_instance(tmp_instance, &c);
- free(breaking_data->break_list);
- free(breaking_data->extern_field);
- free(breaking_data);
- }
-
- for (int i = 0; i < network->ne; i++) {
- if (num_current_skipped[i] < num) {
- avg_current[i] /= num - num_current_skipped[i];
- }
- }
-
- for (int i = 0; i < network->nv; i++) {
- if (num_voltage_skipped[i] < num) {
- avg_voltage[i] /= num - num_voltage_skipped[i];
- }
- }
-
- double *avg_field;
- if (voltage_bound)
- avg_field = avg_voltage;
- else
- avg_field = avg_current;
-
- update_boundary(perm_instance, avg_field);
-
- if (save_stress) {
- char *c_filename = (char *)malloc(100 * sizeof(char));
- snprintf(c_filename, 100, "current_%d_%g_%d_%g.txt", width, crack_len,
- iter, beta);
- FILE *outfile = fopen(c_filename, "w");
- for (int i = 0; i < network->ne; i++) {
- fprintf(outfile, "%g ", avg_current[i]);
- }
- fclose(outfile);
- free(c_filename);
- }
-
- if (save_damage) {
- char *c_filename = (char *)malloc(100 * sizeof(char));
- snprintf(c_filename, 100, "damage_%d_%g_%d_%g.txt", width, crack_len,
- iter, beta);
- FILE *outfile = fopen(c_filename, "w");
- for (int i = 0; i < network->ne; i++) {
- fprintf(outfile, "%g ", damage[i]);
- }
- fclose(outfile);
- free(c_filename);
- }
-
- if (save_strength) {
- char *s_filename = (char *)malloc(100 * sizeof(char));
- snprintf(s_filename, 100, "strength_%d_%g_%d_%g.txt", width, crack_len, iter, beta);
- FILE *f = fopen(s_filename, "a");
- for (int i = 0; i < num; i++) {
- fprintf(f, "%g ", strength[i]);
- }
- fclose(f);
- free(s_filename);
- }
-
- if (save_bound) {
- char *b_filename = (char *)malloc(100 * sizeof(char));
- snprintf(b_filename, 100, "bounds_%d_%g_%d_%g.txt", width, crack_len,
- iter, beta);
- FILE *outfile = fopen(b_filename, "w");
- for (int i = 0; i < network->nv; i++) {
- fprintf(outfile, "%g ", ((double *)perm_instance->boundary_cond->x)[i]);
- }
- fclose(outfile);
- free(b_filename);
- }
-
- free(avg_current);
- free(avg_voltage);
- if (save_damage) free(damage);
- free(num_current_skipped);
- free(num_voltage_skipped);
- if (save_strength) {
- free(strength);
- }
-
- printf(
- "\033[F\033[JCURRENT_SCALING: ITERATION %0*d COMPLETE, BETA = %.2f\n\n",
- (int)log10(iter) + 1, DUMB2 + 1, beta);
- }
-
- free_instance(perm_instance, &c);
- free_net(network, &c);
-
- CHOL_F(finish)(&c);
-
- return 0;
-}