path: root/src/current_scaling.c

#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, 1, &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, 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->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->edge_coords, 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;
}