started repo again without all the data files gunking the works

1 files changed, 300 insertions, 0 deletions

diff --git a/src/current_scaling.c b/src/current_scaling.c
new file mode 100644
index 0000000..6837f6b
--- /dev/null
+++ b/src/current_scaling.c
@@ -0,0 +1,300 @@
+
+#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;
+	}
+
+	fnet *network = ini_square_network(width, false, true, &c);
+	finst *perm_instance =
+			create_instance(network, inf, voltage_bound, false, &c);
+	gen_crack(perm_instance, crack_len, 1, &c);
+	finish_instance(perm_instance, &c);
+
+	if (voltage_bound) {
+		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
+		finst *tmp_instance = create_instance(network, inf, false, false, &c);
+		gen_crack(tmp_instance, crack_len, 1, &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->num_verts; i++) {
+			double tmp;
+			fscanf(bound_f, "%lg ", &tmp);
+			((double *)perm_instance->boundary_cond->x)[i] = tmp;
+		}
+
+		((double *)perm_instance->boundary_cond->x)[network->num_verts] = 0;
+		((double *)perm_instance->boundary_cond->x)[network->num_verts + 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->num_edges, sizeof(double));
+		}
+		double *avg_current = (double *)calloc(network->num_edges, sizeof(double));
+		unsigned int *num_current_skipped =
+				(unsigned int *)calloc(network->num_edges, sizeof(unsigned int));
+		double *avg_voltage = (double *)calloc(network->num_verts, sizeof(double));
+		unsigned int *num_voltage_skipped =
+				(unsigned int *)calloc(network->num_verts, 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);
+
+			break_data *breaking_data = NULL;
+			while (breaking_data == NULL) {
+				double *fuse_thres = gen_fuse_thres(
+						network->num_edges, network->edge_coords, beta, beta_scaling_flat);
+				finst *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]);
+			}
+
+			finst *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->num_edges; i++) {
+				avg_current[i] += current[i];
+				if (current[i] == 0)
+					num_current_skipped[i]++;
+			}
+
+			for (unsigned int i = 0; i < network->num_verts; i++) {
+				if (tmp_instance->marks[i] == tmp_instance->marks[network->num_verts]) {
+					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->num_edges; i++) {
+			if (num_current_skipped[i] < num) {
+				avg_current[i] /= num - num_current_skipped[i];
+			}
+		}
+
+		for (int i = 0; i < network->num_verts; 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->num_edges; 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->num_edges; 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->num_verts; 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_fnet(network, &c);
+
+	CHOL_F(finish)(&c);
+
+	return 0;
+}