more refactoring

1 files changed, 130 insertions, 0 deletions

diff --git a/src/net.c b/src/net.c
new file mode 100644
index 0000000..9f6965b
--- /dev/null
+++ b/src/net.c
@@ -0,0 +1,130 @@
+
+#include "fracture.h"
+
+double *get_thres(uint_t ne, double beta) {
+	double *thres = (double *)malloc(ne * sizeof(double));
+	assert(thres != NULL);
+
+	gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN);
+	gsl_rng_set(r, rand_seed());
+
+	for (uint_t i = 0; i < ne; i++) {
+		thres[i] = rand_dist_pow(r, beta);
+	}
+
+	gsl_rng_free(r);
+
+	return thres;
+}
+
+void net_notch(net_t *net, double notch_len, cholmod_common *c) {
+	for (uint_t i = 0; i < net->graph->ne; i++) {
+		uint_t v1, v2;
+		double v1x, v1y, v2x, v2y, dy;
+		bool crosses_center, not_wrapping, correct_length;
+
+		v1 = net->graph->ev[2 * i];
+		v2 = net->graph->ev[2 * i + 1];
+
+		v1x = net->graph->vx[2 * v1];
+		v1y = net->graph->vx[2 * v1 + 1];
+		v2x = net->graph->vx[2 * v2];
+		v2y = net->graph->vx[2 * v2 + 1];
+
+		dy = v1y - v2y;
+
+		crosses_center = (v1y >= 0.5 && v2y <= 0.5) || (v1y <= 0.5 && v2y >= 0.5);
+		not_wrapping = fabs(dy) < 0.5;
+		//correct_length = v1x + dx / dy * (v1y - 0.5) <= notch_len;
+		correct_length = v1x < notch_len && v2x < notch_len;
+
+		if (crosses_center && not_wrapping && correct_length) {
+			break_edge(net, i, c);
+		}
+	}
+}
+
+net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c) {
+	net_t *net = (net_t *)calloc(1, sizeof(net_t));
+	assert(net != NULL);
+
+	net->graph = g;
+	net->fuses = (bool *)calloc(g->ne, sizeof(bool));
+	assert(net->fuses != NULL);
+	net->thres = get_thres(g->ne, beta);
+	net->inf = inf;
+
+	net->voltage_bound = vb;
+	net->boundary_cond = bound_set(g, vb, notch_len, c);
+
+	if (g->boundary != TORUS_BOUND) net->adjacency = gen_adjacency(net, false, false, 0, c);
+	else net->adjacency = gen_adjacency(net, true, false, 0, c);
+
+	net->marks = (uint_t *)malloc((net->graph->break_dim) * sizeof(uint_t));
+	assert(net->marks != NULL);
+
+	net->dual_marks = (uint_t *)malloc((net->graph->dnv) * sizeof(uint_t));
+	assert(net->dual_marks != NULL);
+
+	for (uint_t i = 0; i < (net->graph->break_dim); i++) {
+		net->marks[i] = 1;
+	}
+	for (uint_t i = 0; i < (net->graph->dnv); i++) {
+		net->dual_marks[i] = i+1;
+	}
+	net->num_components = 1;
+
+	net_notch(net, notch_len, c);
+
+	{
+		cholmod_sparse *laplacian = gen_laplacian(net, c, true);
+		net->factor = CHOL_F(analyze)(laplacian, c);
+		CHOL_F(factorize)(laplacian, net->factor, c);
+		CHOL_F(free_sparse)(&laplacian, c);
+	}
+
+	return net;
+}
+
+net_t *net_copy(const net_t *net, cholmod_common *c) {
+	net_t *net_copy = (net_t *)calloc(1, sizeof(net_t));
+	assert(net_copy != NULL);
+	memcpy(net_copy, net, sizeof(net_t));
+
+	size_t fuses_size = (net->graph)->ne * sizeof(bool);
+	net_copy->fuses = (bool *)malloc(fuses_size);
+	assert(net_copy->fuses != NULL);
+	memcpy(net_copy->fuses, net->fuses, fuses_size);
+
+	size_t thres_size = (net->graph)->ne * sizeof(double);
+	net_copy->thres = (double *)malloc(thres_size);
+	assert(net_copy->thres != NULL);
+	memcpy(net_copy->thres, net->thres, thres_size);
+
+	size_t marks_size = (net->graph->break_dim) * sizeof(uint_t);
+	net_copy->marks = (uint_t *)malloc(marks_size);
+	assert(net_copy->marks != NULL);
+	memcpy(net_copy->marks, net->marks, marks_size);
+
+	size_t dual_marks_size = (net->graph->dnv) * sizeof(uint_t);
+	net_copy->dual_marks = (uint_t *)malloc(dual_marks_size);
+	assert(net_copy->dual_marks != NULL);
+	memcpy(net_copy->dual_marks, net->dual_marks, dual_marks_size);
+
+	net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c);
+	net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c);
+	net_copy->factor = CHOL_F(copy_factor)(net->factor, c);
+
+	return net_copy;
+}
+
+void net_free(net_t *net, cholmod_common *c) {
+	free(net->fuses);
+	free(net->thres);
+	CHOL_F(free_dense)(&(net->boundary_cond), c);
+	CHOL_F(free_sparse)(&(net->adjacency), c);
+	CHOL_F(free_factor)(&(net->factor), c);
+	free(net->marks);
+	free(net->dual_marks);
+	free(net);
+}