1 files changed, 69 insertions, 0 deletions
diff --git a/src/net_create.c b/src/net_create.c
new file mode 100644
index 0000000..fedcb3a
--- /dev/null
+++ b/src/net_create.c
@@ -0,0 +1,69 @@
+
+#include "fracture.h"
+
+double *get_thres(uint_t ne, double beta) {
+	assert(beta > 0);
+
+	double *thres = (double *)malloc(ne * sizeof(double));
+	assert(thres != NULL);
+
+	gsl_set_error_handler_off();
+
+	gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+	{
+		FILE *rf = fopen("/dev/urandom", "r");
+		unsigned long int seed;
+		fread(&seed, sizeof(unsigned long int), 1, rf);
+		fclose(rf);
+		gsl_rng_set(r, seed);
+	}
+
+	for (uint_t i = 0; i < ne; i++) {
+		while ((thres[i] = exp(log(gsl_ran_flat(r, 0, 1)) / beta)) == 0.0);
+	}
+
+	gsl_rng_free(r);
+
+	return thres;
+}
+
+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));
+	net->dual_marks = (uint_t *)malloc((net->graph->dnv) * sizeof(uint_t));
+	assert(net->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;
+}
+