summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/bound_set.c62
-rw-r--r--src/break_edge.c121
-rw-r--r--src/factor_update.c30
-rw-r--r--src/fracture.c16
-rw-r--r--src/fracture.h58
-rw-r--r--src/frame_create.c167
-rw-r--r--src/gen_laplacian.c253
-rw-r--r--src/get_dual_clusters.c2
-rw-r--r--src/graph_create.c717
-rw-r--r--src/graph_free.c13
-rw-r--r--src/graph_genfunc.c20
-rw-r--r--src/net.c74
-rw-r--r--src/net_currents.c18
-rw-r--r--src/net_fracture.c10
-rw-r--r--src/net_voltages.c24
15 files changed, 645 insertions, 940 deletions
diff --git a/src/bound_set.c b/src/bound_set.c
index f513aa1..65f1e6f 100644
--- a/src/bound_set.c
+++ b/src/bound_set.c
@@ -26,35 +26,73 @@ void bound_set_embedded(double *bound, const graph_t *g, double notch_len) {
y3 = (2. * i + 1.) / L - 0.5; x3 = 0.5 - 1.;
y4 = (2. * i + 1.) / L - 0.5; x4 = 0.5 - 0.;
- bound[g->bound_verts[g->bound_inds[0] + i]] = u_y(x1, y1);
- bound[g->bound_verts[g->bound_inds[1] + i]] = u_y(x2, y2);
- bound[g->bound_verts[g->bound_inds[2] + i]] = u_y(x3, y3);
- bound[g->bound_verts[g->bound_inds[3] + i]] = u_y(x4, y4);
+ bound[g->b[g->bi[0] + i]] = u_y(x1, y1);
+ bound[g->b[g->bi[1] + i]] = u_y(x2, y2);
+ bound[g->b[g->bi[2] + i]] = u_y(x3, y3);
+ bound[g->b[g->bi[3] + i]] = u_y(x4, y4);
}
}
+bool is_in(uint_t len, uint_t *list, uint_t element) {
+ for (uint_t i = 0; i < len; i++) {
+ if (list[i] == element) {
+ return true;
+ }
+ }
+ return false;
+}
+
cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_common *c) {
- cholmod_dense *boundary = CHOL_F(zeros)(g->break_dim, 1, CHOLMOD_REAL, c);
+
+ uint_t dim = g->nv;
+
+ if (vb && g->boundary != TORUS_BOUND) {
+ dim -= g->bi[g->nb];
+ } else if (!vb) {
+ dim += 2;
+ }
+
+ cholmod_dense *boundary = CHOL_F(zeros)(dim, 1, CHOLMOD_REAL, c);
double *bound = (double *)boundary->x;
switch (g->boundary) {
case TORUS_BOUND:
- for (uint_t i = 0; i < g->bound_inds[1]; i++) {
- bound[g->bound_verts[i]] = 1;
- bound[g->nv + i] = -1;
+ for (uint_t i = 0; i < g->bi[1]; i++) {
+ uint_t be = g->b[i];
+ uint_t v1 = g->ev[2 * be];
+ uint_t v2 = g->ev[2 * be + 1];
+ double v1y = g->vx[2 * v1 + 1];
+ double v2y = g->vx[2 * v2 + 1];
+
+ uint_t ind = v1y < v2y ? 0 : 1;
+
+ bound[g->ev[2 * be + ind]] += 1;
+ bound[g->ev[2 * be + !ind]] -= 1;
}
- bound[g->bound_verts[0]] = 1;
break;
+ /*
case EMBEDDED_BOUND:
bound_set_embedded(bound, g, notch_len);
break;
+ */
default:
if (vb) {
- for (uint_t i = 0; i < g->bound_inds[1]; i++) {
- bound[g->bound_verts[i]] = 1;
+ for (uint_t i = 0; i < dim; i++) {
+ uint_t v = g->nbi[i];
+ for (uint_t j = 0; j < g->vei[v+1] - g->vei[v]; j++) {
+ uint_t e = g->ve[g->vei[v] + j];
+ uint_t v0 = g->ev[2 * e];
+ uint_t v1 = g->ev[2 * e + 1];
+
+ if (g->bq[v0] || g->bq[v1]) {
+ uint_t vv = v0 == v ? v1 : v0;
+ if (is_in(g->bi[1], g->b, vv)) {
+ bound[i]++;
+ }
+ }
+ }
}
} else {
- bound[0] = 1;
bound[g->nv] = 1;
bound[g->nv + 1] = -1;
}
diff --git a/src/break_edge.c b/src/break_edge.c
index 4e8bdcf..2f112c2 100644
--- a/src/break_edge.c
+++ b/src/break_edge.c
@@ -1,104 +1,51 @@
#include "fracture.h"
-bool break_edge(net_t *instance, uint_t edge, cholmod_common *c) {
- instance->fuses[edge] = true;
- instance->num_broken++;
- if (instance->factor != NULL) {
- uint_t w1 = instance->graph->ev_break[2 * edge];
- uint_t w2 = instance->graph->ev_break[2 * edge + 1];
- factor_update(instance->factor, w1, w2, c);
- }
-
- uint_t v1, v2, s1, s2, dv1, dv2, ds1, ds2;
+bool break_edge(net_t *net, uint_t edge, cholmod_common *c) {
+ net->fuses[edge] = true;
+ net->num_broken++;
- v1 = instance->graph->ev[2 * edge];
- v2 = instance->graph->ev[2 * edge + 1];
- dv1 = instance->graph->dev[2 * edge];
- dv2 = instance->graph->dev[2 * edge + 1];
+ double *x = (double *)net->boundary_cond->x;
+ const graph_t *g = net->graph;
- s1 = v1 > v2 ? v1 : v2;
- s2 = v1 > v2 ? v2 : v1;
- ds1 = dv1 > dv2 ? dv1 : dv2;
- ds2 = dv1 > dv2 ? dv2 : dv1;
+ uint_t v1 = net->graph->ev[2 * edge];
+ uint_t v2 = net->graph->ev[2 * edge + 1];
- {
- int_t *lap_p = (int_t *)instance->adjacency->p;
- int_t *lap_i = (int_t *)instance->adjacency->i;
- double *lap_x = (double *)instance->adjacency->x;
+ uint_t s1 = v1 < v2 ? v1 : v2;
+ uint_t s2 = v1 < v2 ? v2 : v1;
- for (int i = 0; i < lap_p[s1 + 1] - lap_p[s1]; i++) {
- if (lap_i[lap_p[s1] + i] == s2)
- lap_x[lap_p[s1] + i] = 0;
+ if (net->graph->boundary == TORUS_BOUND) {
+ if (net->graph->bq[edge]) {
+ double v1y = g->vx[2 * v1 + 1];
+ double v2y = g->vx[2 * v2 + 1];
+ uint_t ind = v1y < v2y ? 0 : 1;
+ x[g->ev[2 * edge + ind]] -= 1;
+ x[g->ev[2 * edge + !ind]] += 1;
}
- for (int i = 0; i < lap_p[s2 + 1] - lap_p[s2]; i++) {
- if (lap_i[lap_p[s2] + i] == s1)
- lap_x[lap_p[s2] + i] = 0;
+ if (net->factor != NULL) {
+ factor_update(net->factor, s1, s2, c);
}
- }
-
- int_t old_num_components = instance->num_components;
-
- instance->num_components = update_components(
- instance->adjacency, instance->marks, old_num_components, s1, s2, 0);
-
-
- if (instance->graph->boundary == TORUS_BOUND) {
- if (instance->dual_marks[dv1] == instance->dual_marks[dv2]) {
- int **cycles = (int **)malloc(4*instance->graph->ne * sizeof(int *));
- unsigned int num_cycles = find_cycles(instance->graph->ne, instance->fuses, instance->graph->dev, instance->graph->dvei, instance->graph->dve, cycles);
-
- for (unsigned int i = 0; i < num_cycles; i++) {
- int x_num_crossings = 0;
- int y_num_crossings = 0;
- int ee; unsigned int j = 0;
- while ((ee = cycles[2*i][j]) >= 0) {
- int side = cycles[2*i+1][j];
- j++;
- unsigned int v1, v2;
- double v1x, v1y, v2x, v2y;
- v1 = instance->graph->dev[2 * ee + !side];
- v2 = instance->graph->dev[2 * ee + side];
- v1x = instance->graph->dvx[2 * v1];
- v1y = instance->graph->dvx[2 * v1 + 1];
- v2x = instance->graph->dvx[2 * v2];
- v2y = instance->graph->dvx[2 * v2 + 1];
- double dx = v1x - v2x;
- double dy = v1y - v2y;
- if (((v1x > 0.5 && v2x < 0.5) || (v1x < 0.5 && v2x > 0.5)) && fabs(dx) < 0.5) {
- x_num_crossings += dx > 0 ? 1 : -1;
- }
- if (((v1y > 0.5 && v2y < 0.5) || (v1y < 0.5 && v2y > 0.5)) && fabs(dy) < 0.5) {
- y_num_crossings += dy > 0 ? 1 : -1;
- }
- }
- if ((abs(y_num_crossings) == 0 && abs(x_num_crossings) > 0) || (abs(y_num_crossings) > 0 && abs(x_num_crossings) > 0 && num_cycles > 1)) {
- instance->num_components = 2;
- }
- free(cycles[2*i]);
- free(cycles[2*i+1]);
+ } else if (net->voltage_bound) {
+ if (g->bq[v1] || g->bq[v2]) {
+ uint_t vv = g->bq[v1] ? v2 : v1;
+ uint_t uu = v1 == vv ? v2 : v1;
+ if (is_in(g->bi[1], g->b, uu)) {
+ x[g->bni[vv]] -= 1;
+ }
+ if (net->factor != NULL) {
+ factor_update2(net->factor, g->bni[vv], c);
+ }
+ } else {
+ if (net->factor != NULL) {
+ factor_update(net->factor, g->bni[s1], g->bni[s2], c);
}
- free(cycles);
- }
- }
-
- {
- int_t *lap_p = (int_t *)instance->dual_adjacency->p;
- int_t *lap_i = (int_t *)instance->dual_adjacency->i;
- double *lap_x = (double *)instance->dual_adjacency->x;
-
- for (int i = 0; i < lap_p[ds1 + 1] - lap_p[ds1]; i++) {
- if (lap_i[lap_p[ds1] + i] == ds2)
- lap_x[lap_p[ds1] + i] = 1;
}
- for (int i = 0; i < lap_p[ds2 + 1] - lap_p[ds2]; i++) {
- if (lap_i[lap_p[ds2] + i] == ds1)
- lap_x[lap_p[ds2] + i] = 1;
+ } else {
+ if (net->factor != NULL) {
+ factor_update(net->factor, s1, s2, c);
}
}
- set_connected(instance->dual_adjacency, instance->dual_marks, dv1, instance->dual_marks[dv2], -1, 0);
-
return true;
}
diff --git a/src/factor_update.c b/src/factor_update.c
index f27971b..a51705a 100644
--- a/src/factor_update.c
+++ b/src/factor_update.c
@@ -36,3 +36,33 @@ void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common
CHOL_F(free_sparse)(&perm_update_mat, c);
CHOL_F(free_sparse)(&update_mat, c);
}
+
+void factor_update2(cholmod_factor *factor, uint_t v, cholmod_common *c) {
+ uint_t n = factor->n;
+
+ cholmod_sparse *update_mat =
+ CHOL_F(allocate_sparse)(n, n, 1, true, true, 0, CHOLMOD_REAL, c);
+
+ int_t *pp = (int_t *)update_mat->p;
+ int_t *ii = (int_t *)update_mat->i;
+ double *xx = (double *)update_mat->x;
+
+ for (uint_t i = 0; i <= v; i++) {
+ pp[i] = 0;
+ }
+
+ for (uint_t i = v + 1; i <= n; i++) {
+ pp[i] = 1;
+ }
+
+ ii[0] = v;
+ xx[0] = 1;
+
+ cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
+ update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
+
+ CHOL_F(updown)(false, perm_update_mat, factor, c);
+
+ CHOL_F(free_sparse)(&perm_update_mat, c);
+ CHOL_F(free_sparse)(&update_mat, c);
+}
diff --git a/src/fracture.c b/src/fracture.c
index 1ed594f..6ad2f26 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -263,13 +263,6 @@ int main(int argc, char *argv[]) {
graph_t *g = graph_create(lattice, boundary, L, use_dual, &c);
net_t *net = net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c);
- if (net->marks[g->nv] != net->marks[g->nv+1]) {
- if (save_crit_stress) crit_stress[i] = 0;
- if (save_conductivity) conductivity[i] = 0;
- if (save_damage) damage[net->num_broken]++;
- if (save_energy) energy[i] = 0;
- if (save_threshold) thresholds[i] = 0;
- } else {
net_t *tmp_net = net_copy(net, &c);
data_t *data = net_fracture(tmp_net, &c, cutoff);
net_free(tmp_net, &c);
@@ -325,14 +318,6 @@ int main(int argc, char *argv[]) {
}
}
- FILE *testout = fopen("test.txt", "w");
- double *tmp_voltage = net_voltages(net, &c);
- for (uint_t j = 0; j < g->nv; j++) {
- fprintf(testout, "%g ", tmp_voltage[j]);
- }
- fclose(testout);
- free(tmp_voltage);
-
if (save_damage) {
uint_t would_break = 0;
double *tmp_voltage = net_voltages(net, &c);
@@ -382,7 +367,6 @@ int main(int argc, char *argv[]) {
}
data_free(data);
- }
if (save_network) {
FILE *net_out = fopen("network.txt", "w");
for (uint_t j = 0; j < g->nv; j++) {
diff --git a/src/fracture.h b/src/fracture.h
index 8f883bc..0a3a687 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -48,26 +48,35 @@ typedef struct {
uint_t ne;
uint_t nv;
uint_t dnv;
- uint_t nv_break;
- uint_t num_bounds;
+ uint_t *ev;
+ uint_t *dev;
+ double *vx;
+ double *dvx;
+} frame_t;
+
+typedef struct {
+ uint_t L;
+ bound_t boundary;
+ uint_t ne;
+ uint_t nv;
+ uint_t dnv;
+ uint_t nb;
uint_t *ev;
- uint_t *ev_break;
+ uint_t *dev;
+ double *vx;
+ double *dvx;
uint_t *vei;
uint_t *ve;
- uint_t *bound_inds;
- uint_t *bound_verts;
- double *vx;
- double *ex;
- uint_t *dev;
uint_t *dvei;
uint_t *dve;
- double *dvx;
- uint_t num_spanning_edges;
+ uint_t *bi;
+ uint_t *b;
+ uint_t *nbi;
+ uint_t *bni;
+ bool *bq;
uint_t *spanning_edges;
- uint_t L;
- uint_t break_dim;
+ uint_t num_spanning_edges;
cholmod_sparse *voltcurmat;
- bound_t boundary;
} graph_t;
typedef struct {
@@ -77,14 +86,11 @@ typedef struct {
double inf;
cholmod_dense *boundary_cond;
cholmod_factor *factor;
- uint_t *marks;
- uint_t *dual_marks;
bool voltage_bound;
- uint_t num_components;
- cholmod_sparse *adjacency;
- cholmod_sparse *dual_adjacency;
- bool debug_stop;
uint_t num_broken;
+ uint_t dim;
+ uint_t nep;
+ uint_t *evp;
} net_t;
typedef struct {
@@ -101,11 +107,9 @@ int update_components(const cholmod_sparse *laplacian, uint_t *marks,
uint_t *find_components(const cholmod_sparse *laplacian, uint_t skip);
-cholmod_sparse *gen_adjacency(const net_t *instance, bool dual, bool breakv,
- uint_t pad, cholmod_common *c);
+cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, bool symmetric, cholmod_common *c);
-cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c,
- bool symmetric);
+cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c);
int edge_to_verts(uint_t width, bool periodic, uint_t edge,
bool index);
@@ -160,8 +164,8 @@ uint_t *get_clusters(net_t *instance, cholmod_common *c);
uint_t *get_cluster_dist(net_t *instance, cholmod_common *c);
-double *genfunc_uniform(uint_t L, bound_t boundary, gsl_rng *r, uint_t *num);
-double *genfunc_hyperuniform(uint_t L, bound_t boundary, gsl_rng *r, uint_t *num);
+double *genfunc_uniform(uint_t num, gsl_rng *r);
+double *genfunc_hyperuniform(uint_t num, gsl_rng *r);
void randfunc_flat(gsl_rng *r, double *x, double *y);
void randfunc_gaus(gsl_rng *r, double *x, double *y);
@@ -188,3 +192,7 @@ unsigned long int rand_seed();
long double rand_dist_pow(const gsl_rng *r, double beta);
double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate);
+
+frame_t *frame_create(lattice_t lattice, uint_t L, bool dual);
+void frame_free(frame_t *frame);
+bool is_in(uint_t len, uint_t *list, uint_t element);
diff --git a/src/frame_create.c b/src/frame_create.c
new file mode 100644
index 0000000..89ce69d
--- /dev/null
+++ b/src/frame_create.c
@@ -0,0 +1,167 @@
+
+#include "fracture.h"
+
+uint_t *get_voro_dual_edges(uint_t num_edges,
+ uint_t num_verts, uint_t *edges,
+ uint_t *triangles) {
+ uint_t *dual_edges =
+ (uint_t *)malloc(2 * num_edges * sizeof(uint_t));
+ uint_t place = 0;
+ for (uint_t i = 0; i < num_edges; i++) {
+ uint_t v1, v2;
+ v1 = edges[2 * i];
+ v2 = edges[2 * i + 1];
+ if (v1 < num_verts && v2 < num_verts) {
+ bool found_match = false;
+ for (uint_t j = 0; j < 3; j++) {
+ for (uint_t k = 0; k < 3; k++) {
+ uint_t t11, t12, t21, t22;
+ t11 = triangles[3 * v1 + j];
+ t12 = triangles[3 * v1 + ((j + 1) % 3)];
+ t21 = triangles[3 * v2 + k];
+ t22 = triangles[3 * v2 + ((k + 1) % 3)];
+ if ((t11 == t21 && t12 == t22) || (t11 == t22 && t12 == t21)) {
+ dual_edges[2 * place] = t11 < t12 ? t11 : t12;
+ dual_edges[2 * place + 1] = t11 < t12 ? t12 : t11;
+ place++;
+ found_match = true;
+ break;
+ }
+ }
+ if (found_match)
+ break;
+ }
+ }
+ }
+
+ return dual_edges;
+}
+
+frame_t *frame_create_voronoi(uint_t L, bool dual, bool hyperuniform) {
+ double *dvx = NULL;
+ uint_t dnv = 2 * pow(L / 2, 2);
+
+ {
+ gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+ gsl_rng_set(r, rand_seed());
+
+ if (hyperuniform) {
+ dvx = genfunc_hyperuniform(dnv, r);
+ } else {
+ dvx = genfunc_uniform(dnv, r);
+ }
+
+ gsl_rng_free(r);
+ }
+
+ // retrieve a periodic voronoi tesselation of the lattice
+ intptr_t *vout = run_voronoi(dnv, dvx, true, 0, 1, 0, 1);
+
+ uint_t nv = ((uint_t *)vout[0])[0];
+ uint_t ne = ((uint_t *)vout[0])[1];
+ double *vx = (double *)vout[2];
+ uint_t *ev = (uint_t *)vout[3];
+ uint_t *voro_tris = (uint_t *)vout[5];
+
+ free((void *)vout[0]);
+ free((void *)vout[1]);
+ free((void *)vout[4]);
+ free(vout);
+
+ // get dual edges of the fully periodic graph
+ uint_t *dev = get_voro_dual_edges(ne, nv, ev, voro_tris);
+
+ frame_t *frame = (frame_t *)malloc(sizeof(frame_t));
+ frame->ne = ne;
+
+ // when use_dual is specificed, the edge and vertex sets are swapped with the
+ // dual edge and dual vertex sets. once formally relabelled, everything
+ // works the same way
+ if (dual) {
+ frame->nv = dnv;
+ frame->dnv = nv;
+ frame->ev = dev;
+ frame->dev = ev;
+ frame->vx = dvx;
+ frame->dvx = vx;
+ } else {
+ frame->nv = nv;
+ frame->dnv = dnv;
+ frame->ev = ev;
+ frame->dev = dev;
+ frame->vx = vx;
+ frame->dvx = dvx;
+ }
+
+ return frame;
+}
+
+frame_t *frame_create_square(uint_t L, bool dual) {
+ uint_t nv = 2 * pow(L / 2, 2);
+ uint_t ne = pow(L, 2);
+
+ uint_t *ev = (uint_t *)malloc(2 * ne * sizeof(uint_t));
+ uint_t *dev = (uint_t *)malloc(2 * ne * sizeof(uint_t));
+ double *vx = (double *)malloc(2 * nv * sizeof(double));
+ double *dvx = (double *)malloc(2 * nv * sizeof(double));
+
+ for (uint_t i = 0; i < ne; i++) {
+ uint_t x = i / L;
+ uint_t y = i % L;
+
+ ev[2 * i] = (L * x) / 2 + ((y + x % 2) / 2) % (L / 2);
+ ev[2 * i + 1] = ((L * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2)) % nv;
+
+ dev[2 * i] = (L * x) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2);
+ dev[2 * i + 1] = ((L * (x + 1)) / 2 + ((y + x % 2) / 2) % (L / 2)) % nv;
+ }
+
+ double dx = 1. / L;
+
+ for (uint_t i = 0; i < nv; i++) {
+ vx[2 * i] = dx * ((1 + i / (L / 2)) % 2 + 2 * (i % (L / 2)));
+ vx[2 * i + 1] = dx * (i / (L / 2));
+
+ dvx[2 * i] = dx * ((i / (L / 2)) % 2 + 2 * (i % (L / 2)));
+ dvx[2 * i + 1] = dx * (i / (L / 2));
+ }
+
+ frame_t *frame = (frame_t *)malloc(sizeof(frame_t));
+ frame->ne = ne;
+ frame->nv = nv;
+ frame->dnv = nv;
+
+ if (dual) {
+ frame->ev = dev;
+ frame->dev = ev;
+ frame->vx = dvx;
+ frame->dvx = vx;
+ } else {
+ frame->ev = ev;
+ frame->dev = dev;
+ frame->vx = vx;
+ frame->dvx = dvx;
+ }
+
+ return frame;
+}
+
+frame_t *frame_create(lattice_t lattice, uint_t L, bool dual) {
+ switch (lattice) {
+ case (SQUARE_LATTICE):
+ return frame_create_square(L, dual);
+ case (VORONOI_LATTICE):
+ return frame_create_voronoi(L, dual, false);
+ case (VORONOI_HYPERUNIFORM_LATTICE):
+ return frame_create_voronoi(L, dual, true);
+ }
+}
+
+void frame_free(frame_t *frame) {
+ free(frame->ev);
+ free(frame->dev);
+ free(frame->vx);
+ free(frame->dvx);
+ free(frame);
+}
+
diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c
index 411734c..1cc93a4 100644
--- a/src/gen_laplacian.c
+++ b/src/gen_laplacian.c
@@ -1,118 +1,79 @@
#include "fracture.h"
-cholmod_sparse *gen_adjacency(const net_t *instance, bool dual, bool breakv,
- uint_t pad, cholmod_common *c) {
- uint_t ne = instance->graph->ne;
- if (!breakv && instance->graph->boundary != TORUS_BOUND && !dual) {
- ne += instance->graph->bound_inds[2];
+cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, bool symmetric, cholmod_common *c) {
+ const graph_t *g = net->graph;
+ uint_t nv;
+ uint_t ne;
+ uint_t nre;
+ uint_t *ev;
+
+ if (use_gp) {
+ nv = net->dim;
+ ne = net->nep;
+ nre = (int_t)net->nep - (int_t)net->num_broken;
+ ev = net->evp;
+ } else if (dual) {
+ nv = g->dnv;
+ ne = g->ne;
+ nre = net->num_broken;
+ ev = g->dev;
+ } else {
+ nv = g->nv;
+ ne = g->ne;
+ nre = (int_t)g->ne - (int_t)net->num_broken;
+ ev = g->ev;
}
- uint_t nnz = 2 * ne;
+ uint_t nnz = nre;
- uint_t nv;
- if (dual)
- nv = instance->graph->dnv;
- else {
- if (breakv) nv = instance->graph->nv_break;
- else nv = instance->graph->nv;
- if (instance->graph->boundary != TORUS_BOUND && !breakv) {
- nv += 2;
- }
- }
+ cholmod_triplet *t = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);
- cholmod_triplet *t =
- CHOL_F(allocate_triplet)(nv + pad, nv + pad, nnz, 0, CHOLMOD_REAL, c);
int_t *ri = (int_t *)t->i;
int_t *ci = (int_t *)t->j;
double *ai = (double *)t->x;
t->nnz = nnz;
- uint_t *etv;
- if (dual)
- etv = instance->graph->dev;
- else {
- if (breakv) etv = instance->graph->ev_break;
- else etv = instance->graph->ev;
- }
-
- uint_t count = 0;
-
- for (uint_t i = 0; i < instance->graph->ne; i++) {
- if ((instance->fuses[i] && dual) || (!instance->fuses[i] && !dual)) {
- uint_t v1 = etv[2 * i];
- uint_t v2 = etv[2 * i + 1];
+ uint_t a = 0;
- ri[2 * count] = v1;
- ci[2 * count] = v2;
- ai[2 * count] = 1;
+ for (uint_t i = 0; i < ne; i++) {
+ if ((net->fuses[i] && dual) || (!net->fuses[i] && !dual)) {
+ uint_t v1 = ev[2 * i];
+ uint_t v2 = ev[2 * i + 1];
- ri[2 * count + 1] = v2;
- ci[2 * count + 1] = v1;
- ai[2 * count + 1] = 1;
+ uint_t s1 = v1 < v2 ? v1 : v2;
+ uint_t s2 = v1 < v2 ? v2 : v1;
- count++;
- } else {
- uint_t v1 = etv[2 * i];
- uint_t v2 = etv[2 * i + 1];
+ ri[a] = s2;
+ ci[a] = s1;
+ ai[a] = 1;
- ri[2 * count] = v1;
- ci[2 * count] = v2;
- ai[2 * count] = 0;
-
- ri[2 * count + 1] = v2;
- ci[2 * count + 1] = v1;
- ai[2 * count + 1] = 0;
-
- count++;
+ a++;
}
}
- if (!breakv && instance->graph->boundary != TORUS_BOUND && !dual) {
- for (uint_t i = 0; i < 2; i++) {
- for (uint_t j = 0; j < instance->graph->bound_inds[i+1] - instance->graph->bound_inds[i]; j++) {
- ri[2*count] = instance->graph->nv + i;
- ci[2*count] = instance->graph->bound_verts[instance->graph->bound_inds[i] + j];
- ai[2*count] = 1;
- ri[2*count+1] = instance->graph->bound_verts[instance->graph->bound_inds[i] + j];
- ci[2*count+1] = instance->graph->nv + i;
- ai[2*count+1] = 1;
- count++;
- }
- }
- }
-
-
cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c);
CHOL_F(free_triplet)(&t, c);
+ if (!symmetric) {
+ cholmod_sparse *tmp_s = CHOL_F(copy)(s, 0, 1, c);
+ CHOL_F(free_sparse)(&s, c);
+ s = tmp_s;
+ }
+
return s;
}
-cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c,
- bool symmetric) {
- const graph_t *network = instance->graph;
- uint_t num_verts = network->nv_break;
- uint_t num_bounds = network->num_bounds;
- double inf = instance->inf;
- bool voltage_bound = instance->voltage_bound;
- uint_t *bound_inds = network->bound_inds;
- uint_t *bound_verts = network->bound_verts;
-
- uint_t num_gedges;
- if (voltage_bound) {
- num_gedges = network->bound_inds[num_bounds];
- } else {
- num_gedges = network->bound_inds[num_bounds] + (num_bounds - 2) * 2;
- }
- if (network->boundary == TORUS_BOUND) num_gedges = bound_inds[1];
- uint_t num_gverts = network->break_dim;
+cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c) {
+ const graph_t *g = net->graph;
+ uint_t nv = net->dim;
+ uint_t ne = net->nep;
+ uint_t *ev = net->evp;
- uint_t nnz = num_gverts + 2 * num_gedges;
+ uint_t nnz = nv;
- cholmod_triplet *temp_m =
- CHOL_F(allocate_triplet)(num_gverts, num_gverts, nnz, 0, CHOLMOD_REAL, c);
+ cholmod_triplet *temp_m = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);
int_t *rowind = (int_t *)temp_m->i;
int_t *colind = (int_t *)temp_m->j;
@@ -120,126 +81,56 @@ cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c,
temp_m->nnz = nnz;
- for (uint_t i = 0; i < num_gverts; i++) {
+ for (uint_t i = 0; i < nv; i++) {
rowind[i] = i;
colind[i] = i;
acoo[i] = 0;
}
- cholmod_sparse *adjacency = gen_adjacency(instance, false, true, (int)num_gverts - (int)num_verts, c);
- int_t *ia = (int_t *)adjacency->p;
- double *aa = (double *)adjacency->x;
+ cholmod_sparse *adjacency = gen_adjacency(net, false, true, true, c);
- for (uint_t i = 0; i < num_verts; i++) {
- for (uint_t j = 0; j < ia[i + 1] - ia[i]; j++) {
- acoo[i] += aa[ia[i] + j];
- }
- }
+ for (uint_t i = 0; i < ne; i++) {
+ if (!net->fuses[i]){
+ uint_t v1 = ev[2 * i];
+ uint_t v2 = ev[2 * i + 1];
- if (network->boundary == TORUS_BOUND) {
- for (uint_t i = 0; i < network->bound_inds[1]; i++) {
- uint_t vv = network->bound_verts[i];
- rowind[num_gverts + 2*i] = vv;
- colind[num_gverts + 2*i] = network->nv + i;
- acoo[num_gverts + 2*i] = -1;
- rowind[num_gverts + 2*i+1] = network->nv + i;
- colind[num_gverts + 2*i+1] = vv;
- acoo[num_gverts + 2*i+1] = -1;
- acoo[vv]++;
- acoo[network->nv + i]++;
+ acoo[v1]++;
+ acoo[v2]++;
}
- acoo[bound_verts[0]]++;
}
- if (network->boundary != TORUS_BOUND) {
- if (voltage_bound) {
- for (uint_t i = 0; i < num_bounds; i++) {
- for (uint_t j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) {
- acoo[bound_verts[bound_inds[i] + j]]++;
-
- rowind[nnz - 1 - 2 * (bound_inds[i] + j)] = bound_verts[bound_inds[i] + j];
- colind[nnz - 1 - 2 * (bound_inds[i] + j)] = num_gverts - 1;
- acoo[nnz - 1 - 2 * (bound_inds[i] + j)] = -1;
- rowind[nnz - 1 - 2 * (bound_inds[i] + j) - 1] = num_gverts - 1;
- colind[nnz - 1 - 2 * (bound_inds[i] + j) - 1] = bound_verts[bound_inds[i] + j];
- acoo[nnz - 1 - 2 * (bound_inds[i] + j) - 1] = -1;
- acoo[num_gverts - 1]++;
- acoo[bound_verts[bound_inds[i] + j]]++;
- }
- }
- } else {
- for (uint_t i = 0; i < num_bounds; i++) {
- if (network->boundary != EMBEDDED_BOUND) acoo[0]++;
-
- rowind[num_verts + i] = num_verts + i;
- colind[num_verts + i] = num_verts + i;
- if (i < 2)
- acoo[num_verts + i] = (network->bound_inds[i + 1] -
- network->bound_inds[i] + network->num_bounds - 2) *
- inf;
- else
- acoo[num_verts + i] =
- (network->bound_inds[i + 1] - network->bound_inds[i] + 2) * inf;
- }
+ if (net->voltage_bound && g->boundary != TORUS_BOUND) {
+ for (uint_t i = 0; i < net->dim; i++) {
+ uint_t v = g->nbi[i];
+ for (uint_t j = 0; j < g->vei[v+1] - g->vei[v]; j++) {
+ uint_t e = g->ve[g->vei[v] + j];
+ uint_t v0 = g->ev[2 * e];
+ uint_t v1 = g->ev[2 * e + 1];
- uint_t start = num_gverts;
-
- for (uint_t i = 0; i < 2; i++) {
- for (uint_t j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) {
- rowind[start + bound_inds[i] + j] = bound_verts[bound_inds[i] + j];
- colind[start + bound_inds[i] + j] = num_verts + i;
- acoo[start + bound_inds[i] + j] = -inf;
- acoo[bound_verts[bound_inds[i] + j]] += inf;
- colind[start + bound_inds[num_bounds] + bound_inds[i] + j] =
- bound_verts[bound_inds[i] + j];
- rowind[start + bound_inds[num_bounds] + bound_inds[i] + j] =
- num_verts + i;
- acoo[start + bound_inds[num_bounds] + bound_inds[i] + j] = -inf;
+ if (g->bq[v0] || g->bq[v1]) {
+ acoo[i]++;
}
}
-
- for (uint_t i = 0; i < num_bounds - 2; i++) {
- rowind[nnz - 2 * i - 1] = num_verts;
- colind[nnz - 2 * i - 1] = num_verts + 2 + i;
- acoo[nnz - 2 * i - 1] = -inf;
- rowind[nnz - 2 * i - 2] = num_verts + 1;
- colind[nnz - 2 * i - 2] = num_verts + 2 + i;
- acoo[nnz - 2 * i - 2] = -inf;
- colind[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = num_verts;
- rowind[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = num_verts + 2 + i;
- acoo[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = -inf;
- colind[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = num_verts + 1;
- rowind[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = num_verts + 2 + i;
- acoo[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = -inf;
- }
}
+ } else {
+ acoo[0]++;
}
-
- for (uint_t i = 0; i < num_gverts; i++) {
- if (instance->marks[i] != instance->marks[network->nv])
- acoo[i]++;
+ for (uint_t i = 0; i < nv; i++) {
+ if (acoo[i] == 0) acoo[i]++;
}
- assert(CHOL_F(check_triplet)(temp_m, c));
+ //assert(CHOL_F(check_triplet)(temp_m, c));
cholmod_sparse *t_out = CHOL_F(triplet_to_sparse)(temp_m, nnz, c);
- assert(CHOL_F(check_sparse)(t_out, c));
+ //assert(CHOL_F(check_sparse)(t_out, c));
double alpha[2] = {1, 0};
double beta[2] = {-1, 0};
- cholmod_sparse *t_laplacian =
- CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c);
-
- cholmod_sparse *laplacian;
- if (symmetric)
- laplacian = CHOL_F(copy)(t_laplacian, 1, 1, c);
- else
- laplacian = CHOL_F(copy)(t_laplacian, 0, 1, c);
+ cholmod_sparse *laplacian = CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c);
CHOL_F(free_sparse)(&t_out, c);
CHOL_F(free_sparse)(&adjacency, c);
- CHOL_F(free_sparse)(&t_laplacian, c);
CHOL_F(free_triplet)(&temp_m, c);
return laplacian;
diff --git a/src/get_dual_clusters.c b/src/get_dual_clusters.c
index 3a51c38..78bf185 100644
--- a/src/get_dual_clusters.c
+++ b/src/get_dual_clusters.c
@@ -2,7 +2,7 @@
#include "fracture.h"
unsigned int *get_clusters(net_t *instance, cholmod_common *c) {
- cholmod_sparse *s_dual = gen_adjacency(instance, true, false, 0, c);
+ cholmod_sparse *s_dual = gen_adjacency(instance, true, false, true, c);
unsigned int *dual_marks = find_components(s_dual, 0);
CHOL_F(free_sparse)(&s_dual, c);
diff --git a/src/graph_create.c b/src/graph_create.c
index c1f556c..09a3aed 100644
--- a/src/graph_create.c
+++ b/src/graph_create.c
@@ -97,586 +97,203 @@ uint_t *get_verts_to_edges(uint_t num_verts, uint_t num_edges,
return output;
}
-graph_t *ini_square_network(uint_t width, bound_t boundary, bool side_bounds,
- cholmod_common *c) {
- graph_t *network = (graph_t *)calloc(1, sizeof(graph_t));
-
- network->boundary = boundary;
- bool periodic = (boundary == CYLINDER_BOUND) || (boundary == TORUS_BOUND) ? true : false;
-
- network->ne = pow(width, 2);
- if (boundary == CYLINDER_BOUND) {
- assert(width % 2 == 0);
- assert(!side_bounds);
- network->nv = (width / 2) * (width + 1);
- network->nv_break = (width / 2) * (width + 1);
- network->dnv = (width / 2 + 1) * (width / 2) + pow(width / 2, 2);
- network->num_bounds = 2;
- network->bound_inds = (uint_t *)malloc((network->num_bounds + 1) *
- sizeof(uint_t));
- network->bound_inds[0] = 0;
- network->bound_inds[1] = width / 2;
- network->bound_inds[2] = width;
- network->bound_verts = (uint_t *)calloc(width, sizeof(uint_t));
- network->break_dim = network->nv + network->num_bounds;
- } else if (boundary == FREE_BOUND || boundary == EMBEDDED_BOUND) {
- network->nv = 2 * ((width + 1) / 2) * (width / 2 + 1);
- network->nv_break = 2 * ((width + 1) / 2) * (width / 2 + 1);
- network->dnv = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2);
- if (side_bounds)
- network->num_bounds = 4;
- else
- network->num_bounds = 2;
- network->bound_inds = (uint_t *)malloc((network->num_bounds + 1) *
- sizeof(uint_t));
- for (uint_t i = 0; i < network->num_bounds + 1; i++) {
- network->bound_inds[i] = i * ((width + 1) / 2);
- }
- network->bound_verts = (uint_t *)malloc(
- network->num_bounds * ((width + 1) / 2) * sizeof(uint_t));
- if (side_bounds) {
- for (uint_t i = 0; i < (width + 1) / 2; i++) {
- network->bound_verts[2 * ((width + 1) / 2) + i] =
- (width + 1) / 2 + i * (width + 1);
- if (width % 2) {
- network->bound_verts[3 * ((width + 1) / 2) + i] =
- (width + 1) / 2 + i * (width + 1) - 1;
- } else {
- network->bound_verts[3 * ((width + 1) / 2) + i] =
- (width + 1) / 2 + i * (width + 1) + width / 2;
- }
- }
- }
- network->break_dim = network->nv + 2;
- } else if (boundary == TORUS_BOUND) {
- network->nv = (width / 2) * (width + 1) - (width / 2);
- network->nv_break = (width / 2) * (width + 1);
- network->dnv = (width / 2 + 1) * (width / 2) + pow(width / 2, 2) - (width / 2);
- network->num_bounds = 1;
- network->bound_inds = (uint_t *)malloc((network->num_bounds + 1) *
- sizeof(uint_t));
- network->bound_inds[0] = 0;
- network->bound_inds[1] = width / 2;
- network->bound_verts = (uint_t *)calloc(width / 2, sizeof(uint_t));
- network->break_dim = network->nv_break + 1;
- }
- if (boundary != TORUS_BOUND) {
- for (uint_t i = 0; i < (width + 1) / 2; i++) {
- network->bound_verts[i] = i;
- network->bound_verts[(width + 1) / 2 + i] = network->nv - (width + 1) / 2 + i;
- }
- } else {
- for (uint_t i = 0; i < width / 2; i++) {
- network->bound_verts[i] = i;
- }
- }
- network->ev_break =
- (uint_t *)calloc(2 * network->ne, sizeof(uint_t));
- network->ev =
- (uint_t *)calloc(2 * network->ne, sizeof(uint_t));
- for (uint_t i = 0; i < network->ne; i++) {
- network->ev_break[2 * i] = edge_to_verts(width, periodic, i, 1);
- network->ev_break[2 * i + 1] = edge_to_verts(width, periodic, i, 0);
- network->ev[2 * i] = network->ev_break[2 * i] % network->nv;
- network->ev[2 * i + 1] = network->ev_break[2 * i + 1] % network->nv;
- }
- network->vei =
- (uint_t *)calloc(network->nv + 1, sizeof(uint_t));
- network->vei[0] = 0;
- uint_t pos1 = 0;
- for (uint_t i = 0; i < network->nv; i++) {
- bool in_bound = false;
- for (uint_t j = 0; j < network->num_bounds; j++) {
- for (uint_t k = 0;
- k < network->bound_inds[j + 1] - network->bound_inds[j]; k++) {
- if (i == network->bound_verts[network->bound_inds[j] + k]) {
- in_bound = true;
- break;
- }
+uint_t get_cut_edges(uint_t ne, const uint_t *ev, const double *vx, bool both, uint_t *ce) {
+ uint_t nce = 0;
+
+ for (uint_t i = 0; i < ne; i++) {
+ uint_t v1 = ev[2 * i];
+ uint_t v2 = ev[2 * i + 1];
+
+ double v1y = vx[2 * v1 + 1];
+ double v2y = vx[2 * v2 + 1];
+
+ if (fabs(v1y - v2y) > 0.5) {
+ ce[nce] = i;
+ nce++;
+ } else if (both) {
+ double v1x = vx[2 * v1];
+ double v2x = vx[2 * v2];
+
+ if (fabs(v1x - v2x) > 0.5) {
+ ce[nce] = i;
+ nce++;
}
}
- if (in_bound)
- pos1 += 2;
- else
- pos1 += 4;
-
- network->vei[i + 1] = pos1;
- }
- network->ve = (uint_t *)calloc(
- network->vei[network->nv], sizeof(uint_t));
- uint_t *vert_counts =
- (uint_t *)calloc(network->nv, sizeof(uint_t));
- for (uint_t i = 0; i < network->ne; i++) {
- uint_t v0 = network->ev[2 * i];
- uint_t v1 = network->ev[2 * i + 1];
- network->ve[network->vei[v0] + vert_counts[v0]] =
- i;
- network->ve[network->vei[v1] + vert_counts[v1]] =
- i;
- vert_counts[v0]++;
- vert_counts[v1]++;
- }
- free(vert_counts);
-
- network->vx =
- (double *)malloc(2 * network->nv * sizeof(double));
- for (uint_t i = 0; i < network->nv; i++) {
- if (!periodic) {
- network->vx[2 * i + 1] = ((double)((2 * i + 1) / (width + 1)))/width;
- network->vx[2 * i] = ((double)((2 * i + 1) % (width + 1)))/width;
- } else {
- network->vx[2 * i + 1] = ((double)((2 * i + 1) / (width)))/width;
- network->vx[2 * i] = ((double)((((2 * i + 1) / (width)+1) % 2) + ((2 * i) % (width))))/width;
- }
}
- network->L = width;
- network->ex = get_edge_coords(
- network->ne, network->vx, network->ev);
-
- network->dev =
- (uint_t *)malloc(2 * network->ne * sizeof(uint_t));
- for (uint_t i = 0; i < network->ne; i++) {
- network->dev[2 * i] =
- dual_edge_to_verts(width, periodic, i, 0) % network->dnv;
- network->dev[2 * i + 1] =
- dual_edge_to_verts(width, periodic, i, 1) % network->dnv;
- }
- network->dvx =
- (double *)malloc(2 * network->dnv * sizeof(double));
- for (uint_t i = 0; i < network->dnv; i++) {
- network->dvx[2 * i + 1] =
- dual_vert_to_coord(width, periodic, i, 0) / width;
- network->dvx[2 * i] =
- dual_vert_to_coord(width, periodic, i, 1) / width;
- }
+ return nce;
+}
- network->voltcurmat = gen_voltcurmat(network->ne,
- network->break_dim,
- network->ev_break, c);
+graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c) {
+ graph_t *g = (graph_t *)calloc(1, sizeof(graph_t));
+ frame_t *f = frame_create(lattice, L, dual);
- network->dvei =
- get_verts_to_edges_ind(network->dnv, network->ne,
- network->dev);
- network->dve = get_verts_to_edges(
- network->dnv, network->ne, network->dev,
- network->dvei);
- network->spanning_edges = get_spanning_edges(network->ne, network->ev, network->vx, 0.51, &(network->num_spanning_edges));
+ g->L = L;
+ g->boundary = bound;
+ g->ne = f->ne;
+
+ if (bound == TORUS_BOUND) {
+ g->nv = f->nv;
+ g->dnv = f->dnv;
+ g->nb = 1;
+
+ g->ev = f->ev;
+ f->ev = NULL;
+ g->dev = f->dev;
+ f->dev = NULL;
+ g->vx = f->vx;
+ f->vx = NULL;
+ g->dvx = f->dvx;
+ f->dvx = NULL;
+
+ // the boundary for the torus consists of a list of edges required to cut
+ // the torus into a cylinder
+ g->bi = (uint_t *)calloc(2, sizeof(uint_t));
+ g->b = (uint_t *)malloc(g->ne * sizeof(uint_t));
+ g->bi[1] = get_cut_edges(g->ne, g->ev, g->vx, false, g->b);
+ g->bq = (bool *)calloc(g->ne, sizeof(bool));
+ for (uint_t i = 0; i < g->bi[1]; i++) {
+ g->bq[g->b[i]] = true;
+ }
+ } else {
+ uint_t *ce = (uint_t *)malloc(g->ne * sizeof(uint_t));
+ uint_t nce = 0;
- return network;
-}
+ if (bound == CYLINDER_BOUND) {
+ g->nb = 2;
+ nce = get_cut_edges(f->ne, f->ev, f->vx, false, ce);
+ } else {
+ g->nb = 4;
+ nce = get_cut_edges(f->ne, f->ev, f->vx, true, ce);
+ }
-uint_t *get_voro_dual_edges(uint_t num_edges,
- uint_t num_verts, uint_t *edges,
- uint_t *triangles) {
- uint_t *dual_edges =
- (uint_t *)malloc(2 * num_edges * sizeof(uint_t));
- uint_t place = 0;
- for (uint_t i = 0; i < num_edges; i++) {
- uint_t v1, v2;
- v1 = edges[2 * i];
- v2 = edges[2 * i + 1];
- if (v1 < num_verts && v2 < num_verts) {
- bool found_match = false;
- for (uint_t j = 0; j < 3; j++) {
- for (uint_t k = 0; k < 3; k++) {
- uint_t t11, t12, t21, t22;
- t11 = triangles[3 * v1 + j];
- t12 = triangles[3 * v1 + ((j + 1) % 3)];
- t21 = triangles[3 * v2 + k];
- t22 = triangles[3 * v2 + ((k + 1) % 3)];
- if ((t11 == t21 && t12 == t22) || (t11 == t22 && t12 == t21)) {
- dual_edges[2 * place] = t11 < t12 ? t11 : t12;
- dual_edges[2 * place + 1] = t11 < t12 ? t12 : t11;
- place++;
- found_match = true;
- break;
- }
- }
- if (found_match)
- break;
+ g->nv = f->nv;
+ g->dnv = f->dnv;
+ g->vx = (double *)malloc(2 * (f->nv + nce) * sizeof(double));
+ g->dvx = (double *)malloc(2 * (f->dnv + nce) * sizeof(double));
+ g->ev = f->ev;
+ g->dev = f->dev;
+ f->ev = NULL;
+ f->dev = NULL;
+ memcpy(g->vx, f->vx, 2 * f->nv * sizeof(double));
+ memcpy(g->dvx, f->dvx, 2 * f->dnv * sizeof(double));
+
+ uint_t nbv = 0;
+ uint_t *bv = (uint_t *)calloc(f->nv, sizeof(uint_t));
+ uint_t *dbv = (uint_t *)calloc(f->dnv, sizeof(uint_t));
+ uint_t nside = 0;
+ bool *side = (bool *)calloc(f->nv, sizeof(bool));
+
+ for (uint_t i = 0; i < nce; i++) {
+ uint_t cv1 = g->ev[2 * ce[i]];
+ uint_t cv2 = g->ev[2 * ce[i] + 1];
+ uint_t dv1 = g->dev[2 * ce[i]];
+ uint_t dv2 = g->dev[2 * ce[i] + 1];
+
+ uint_t cin = 1;
+
+ if (bound == FREE_BOUND && (f->vx[2 * cv2] - f->vx[2 * cv1]) > 0.5) {
+ cin = 0;
}
- }
- }
- return dual_edges;
-}
+ uint_t vin = f->vx[2 * cv1 + cin] < f->vx[2 * cv2 + cin] ? 0 : 1;
+ uint_t dvin = f->dvx[2 * dv1 + cin] < f->dvx[2 * dv2 + cin] ? 0 : 1;
-graph_t *ini_voro_graph(uint_t L, bound_t boundary, bool use_dual,
- double *(*genfunc)(uint_t, bound_t, gsl_rng *, uint_t *),
- cholmod_common *c) {
- graph_t *g = (graph_t *)calloc(1, sizeof(graph_t));
+ if (bv[g->ev[2 * ce[i] + vin]] == 0) {
+ nbv++;
+ if (cin == 0) {
+ nside++;
+ side[g->ev[2 * ce[i] + vin]] = true;
+ }
- // generate the dual lattice
- double *lattice;
- uint_t num;
- {
- 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);
- lattice = genfunc(L, boundary, r, &num);
- gsl_rng_free(r);
- }
+ bv[g->ev[2 * ce[i] + vin]] = g->nv;
- // retrieve a periodic voronoi tesselation of the lattice
- bool run_periodic;
- if (boundary == EMBEDDED_BOUND) run_periodic = false;
- else run_periodic = true;
- intptr_t *vout = run_voronoi(num, lattice, run_periodic, 0, 1, 0, 1);
-
- uint_t tmp_num_verts = ((uint_t *)vout[0])[0];
- uint_t tmp_num_edges = ((uint_t *)vout[0])[1];
- double *tmp_vert_coords = (double *)vout[2];
- uint_t *tmp_edges = (uint_t *)vout[3];
- uint_t *tmp_tris = (uint_t *)vout[5];
-
- free((void *)vout[0]);
- free((void *)vout[1]);
- free((void *)vout[4]);
- free(vout);
-
- // get dual edges of the fully periodic graph
- uint_t *tmp_dual_edges =
- get_voro_dual_edges(tmp_num_edges, tmp_num_verts, tmp_edges, tmp_tris);
-
- // when use_dual is specificed, the edge and vertex sets are swapped with the
- // dual edge and dual vertex sets. once formally relabelled, everything
- // works the same way
- if (use_dual) {
- uint_t *tmp_tmp_dual_edges = tmp_edges;
- double *tmp_lattice = tmp_vert_coords;
- uint_t tmp_num = tmp_num_verts;
-
- tmp_edges = tmp_dual_edges;
- tmp_dual_edges = tmp_tmp_dual_edges;
-
- tmp_vert_coords = lattice;
- lattice = tmp_lattice;
-
- tmp_num_verts = num;
- num = tmp_num;
- }
+ g->vx[2 * g->nv + cin] = 1 + f->vx[2 * g->ev[2 * ce[i] + vin] + cin];
+ g->vx[2 * g->nv + !cin] = f->vx[2 * g->ev[2 * ce[i] + vin] + !cin];
+ g->ev[2 * ce[i] + vin] = g->nv;
- // prune the edges of the lattice and assign boundary vertices based on the
- // desired boundary conditions
- uint_t num_bounds;
- uint_t num_verts;
- double *vert_coords;
- uint_t *bound_inds;
- uint_t *bound_verts;
- uint_t num_edges;
- uint_t *edges;
- uint_t *dual_edges;
- switch (boundary) {
- case FREE_BOUND: {
- num_bounds = 4;
- bound_inds = (uint_t *)malloc((1 + num_bounds) * sizeof(uint_t));
- bound_inds[0] = 0;
- vert_coords = tmp_vert_coords;
- num_verts = tmp_num_verts;
- num_edges = 0;
- edges = (uint_t *)malloc(2 * tmp_num_edges * sizeof(uint_t));
- dual_edges = (uint_t *)malloc(2 * tmp_num_edges * sizeof(uint_t));
- uint_t num_t, num_b, num_l, num_r;
- bool *bound_top, *bound_b, *bound_l, *bound_r;
- num_t = 0; num_b = 0; num_l = 0; num_r = 0;
- bound_top = (bool *)calloc(num_verts, sizeof(bool));
- bound_b = (bool *)calloc(num_verts, sizeof(bool));
- bound_l = (bool *)calloc(num_verts, sizeof(bool));
- bound_r = (bool *)calloc(num_verts, sizeof(bool));
- for (uint_t i = 0; i < tmp_num_edges; i++) {
- uint_t v1, v2;
- double v1x, v1y, v2x, v2y, dx, dy;
- v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1];
- v1x = vert_coords[2 * v1]; v1y = vert_coords[2 * v1 + 1];
- v2x = vert_coords[2 * v2]; v2y = vert_coords[2 * v2 + 1];
- dx = v1x - v2x; dy = v1y - v2y;
- if (fabs(dy) > 0.5) {
- if (dy > 0) {
- if (!bound_top[v1] && !bound_l[v1] && !bound_r[v1]) {
- bound_top[v1] = true; num_t++;
- }
- if (!bound_b[v2] && !bound_l[v2] && !bound_r[v2]) {
- bound_b[v2] = true; num_b++;
- }
- } else {
- if (!bound_top[v2] && !bound_l[v2] && !bound_r[v2]) {
- bound_top[v2] = true; num_t++;
- }
- if (!bound_b[v1] && !bound_l[v1] && !bound_r[v1]) {
- bound_b[v1] = true; num_b++;
- }
- }
- } else if (fabs(dx) > 0.5) {
- if (dx > 0) {
- if (!bound_r[v1] && !bound_top[v1] && !bound_b[v1]) {
- bound_r[v1] = true; num_r++;
- }
- if (!bound_l[v2] && !bound_top[v2] && !bound_b[v2]) {
- bound_l[v2] = true; num_l++;
- }
- } else {
- if (!bound_r[v2] && !bound_top[v2] && !bound_b[v2]) {
- bound_r[v2] = true; num_r++;
- }
- if (!bound_l[v1] && !bound_top[v1] && !bound_b[v1]) {
- bound_l[v1] = true; num_l++;
- }
- }
- } else {
- edges[2 * num_edges] = v1 < v2 ? v1 : v2;
- edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1;
- uint_t d1 = tmp_dual_edges[2 * i];
- uint_t d2 = tmp_dual_edges[2 * i + 1];
- dual_edges[2 * num_edges] = d1 < d2 ? d1 : d2;
- dual_edges[2 * num_edges + 1] = d1 < d2 ? d2 : d1;
- num_edges++;
- }
- }
- bound_verts = malloc((num_t + num_b + num_l + num_r) * sizeof(uint_t));
- bound_inds[1] = num_t; bound_inds[2] = num_t + num_b;
- bound_inds[3] = num_l + num_t + num_b;
- bound_inds[4] = num_t + num_b + num_r + num_l;
- uint_t pos_t, pos_b, pos_l, pos_r;
- pos_t = 0; pos_b = 0; pos_l = 0; pos_r = 0;
- for (uint_t i = 0; i < num_verts; i++) {
- if (bound_top[i]) {
- bound_verts[pos_t] = i; pos_t++;
- } else if (bound_b[i]) {
- bound_verts[num_t + pos_b] = i; pos_b++;
- } else if (bound_l[i]) {
- bound_verts[num_t + num_b + pos_l] = i; pos_l++;
- } else if (bound_r[i]) {
- bound_verts[num_t + num_b + num_l + pos_r] = i; pos_r++;
- }
+ g->nv++;
+ } else {
+ g->ev[2 * ce[i] + vin] = bv[g->ev[2 * ce[i] + vin]];
}
- free(bound_l); free(bound_r); free(bound_top); free(bound_b);
- free(tmp_edges);
- free(tmp_dual_edges);
- num_bounds = 2;
- g->ev_break = edges;
- g->ev = edges;
- g->nv_break = num_verts;
- g->nv = num_verts;
- g->break_dim = num_verts + num_bounds;
- break;
- }
- case CYLINDER_BOUND: {
- num_bounds = 2;
- bound_inds = (uint_t *)malloc((1 + num_bounds) * sizeof(uint_t));
- bound_inds[0] = 0;
- vert_coords = tmp_vert_coords;
- num_verts = tmp_num_verts;
- num_edges = 0;
- edges = (uint_t *)malloc(2 * tmp_num_edges * sizeof(uint_t));
- dual_edges = (uint_t *)malloc(2 * tmp_num_edges * sizeof(uint_t));
- uint_t num_t, num_b;
- bool *bound_top, *bound_b;
- num_t = 0; num_b = 0;
- bound_top = (bool *)calloc(num_verts, sizeof(bool));
- bound_b = (bool *)calloc(num_verts, sizeof(bool));
- for (uint_t i = 0; i < tmp_num_edges; i++) {
- uint_t v1, v2;
- double v1y, v2y, dy;
- v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1];
- v1y = vert_coords[2 * v1 + 1]; v2y = vert_coords[2 * v2 + 1];
- dy = v1y - v2y;
- if (fabs(dy) > 0.5) {
- if (dy > 0) {
- if (!bound_top[v1]) {
- bound_top[v1] = true; num_t++;
- }
- if (!bound_b[v2]) {
- bound_b[v2] = true; num_b++;
- }
- } else {
- if (!bound_top[v2]) {
- bound_top[v2] = true; num_t++;
- }
- if (!bound_b[v1]) {
- bound_b[v1] = true; num_b++;
- }
- }
+ if (dbv[g->dev[2 * ce[i] + dvin]] == 0) {
+ dbv[g->dev[2 * ce[i] + dvin]] = g->dnv;
+
+ if (f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin] < 0.5) {
+ g->dvx[2 * g->dnv + cin] = 1 + f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin];
} else {
- edges[2 * num_edges] = v1 < v2 ? v1 : v2;
- edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1;
- uint_t d1 = tmp_dual_edges[2 * i];
- uint_t d2 = tmp_dual_edges[2 * i + 1];
- dual_edges[2 * num_edges] = d1 < d2 ? d1 : d2;
- dual_edges[2 * num_edges + 1] = d1 < d2 ? d2 : d1;
- num_edges++;
- }
- }
- bound_verts = malloc((num_t + num_b) * sizeof(uint_t));
- bound_inds[1] = num_t; bound_inds[2] = num_t + num_b;
- uint_t pos_t, pos_b;
- pos_t = 0; pos_b = 0;
- for (uint_t i = 0; i < num_verts; i++) {
- if (bound_top[i]) {
- bound_verts[pos_t] = i; pos_t++;
- } else if (bound_b[i]) {
- bound_verts[num_t + pos_b] = i; pos_b++;
+ g->dvx[2 * g->dnv + cin] = f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin];
}
+ g->dvx[2 * g->dnv + !cin] = f->dvx[2 * g->dev[2 * ce[i] + dvin] + !cin];
+ g->dev[2 * ce[i] + dvin] = g->dnv;
+
+ g->dnv++;
+ } else {
+ g->dev[2 * ce[i] + dvin] = dbv[g->dev[2 * ce[i] + dvin]];
}
- free(bound_top); free(bound_b);
- free(tmp_edges);
- free(tmp_dual_edges);
- g->ev_break = edges;
- g->ev = edges;
- g->nv_break = num_verts;
- g->nv = num_verts;
- g->break_dim = num_verts + num_bounds;
- break;
- }
- case TORUS_BOUND: {
- num_bounds = 1;
- bound_inds = (uint_t *)malloc((1 + num_bounds) * sizeof(uint_t));
- bound_inds[0] = 0;
- num_edges = tmp_num_edges;
- edges = (uint_t *)malloc(2* num_edges*sizeof(uint_t));
- for (uint_t i = 0; i < num_edges; i++) {
- edges[2*i] = tmp_edges[2*i];
- edges[2*i+1] = tmp_edges[2*i+1];
- }
- dual_edges = tmp_dual_edges;
- bool *bound_top = (bool *)calloc(tmp_num_verts, sizeof(bool));
- int *edge_change = (int *)calloc(num_edges, sizeof(int));
- uint_t num_t = 0;
- for (uint_t i = 0; i < num_edges; i++) {
- uint_t v1, v2;
- double v1y, v2y, dy;
- v1 = edges[2 * i]; v2 = edges[2 * i + 1];
- v1y = tmp_vert_coords[2 * v1 + 1];
- v2y = tmp_vert_coords[2 * v2 + 1];
- dy = v1y - v2y;
- if (fabs(dy) > 0.5) {
- if (dy > 0) {
- edge_change[i] = 1;
- if (!bound_top[v1]) {
- bound_top[v1] = true; num_t++;
- }
- } else {
- edge_change[i] = 2;
- if (!bound_top[v2]) {
- bound_top[v2] = true; num_t++;
- }
- }
- }
- }
- num_verts = tmp_num_verts + num_t;
- vert_coords = (double *)malloc(2 * num_verts * sizeof(double));
- bound_verts = malloc(num_t * sizeof(uint_t));
- bound_inds[1] = num_t;
- uint_t pos_t = 0;
- for (uint_t i = 0; i < tmp_num_verts; i++) {
- vert_coords[2*i] = tmp_vert_coords[2*i];
- vert_coords[2*i+1] = tmp_vert_coords[2*i+1];
- if (bound_top[i]) {
- bound_verts[pos_t] = i;
- vert_coords[2*(tmp_num_verts + pos_t)] = tmp_vert_coords[2*i];
- vert_coords[2*(tmp_num_verts + pos_t)+1] = tmp_vert_coords[2*i+1];
- pos_t++;
- }
- }
- for (uint_t i = 0; i < num_edges; i++) {
- if (edge_change[i]) {
- for (uint_t j = 0; j < num_t; j++) {
- if (edges[2*i+(edge_change[i]-1)] == bound_verts[j]) {
- edges[2*i+(edge_change[i]-1)] = tmp_num_verts + j;
- break;
- }
- }
- }
- }
- free(tmp_vert_coords);
- free(bound_top);
- free(edge_change);
- g->nv_break = num_verts;
- g->nv = tmp_num_verts;
- g->ev_break = edges;
- g->ev = tmp_edges;
- g->break_dim = num_verts + 1;
- break;
- }
- case EMBEDDED_BOUND: {
- num_bounds = 4;
- bound_inds = (uint_t *)malloc(5 * sizeof(uint_t));
- bound_verts = (uint_t *)malloc(2 * L * sizeof(uint_t));
- for (uint_t i = 0; i < 5; i++) bound_inds[i] = i * L / 2;
- for (uint_t i = 0; i < 2 * L; i++) bound_verts[i] = i;
- uint_t num_away = 0;
- for (uint_t i = 0; i < tmp_num_edges; i++) {
- if (tmp_dual_edges[2*i] > num || tmp_dual_edges[2*i+1] > num) num_away++;
- }
- num_edges = (int)tmp_num_edges - (int)num_away;
- num_verts = tmp_num_verts;
- edges = tmp_edges;
- dual_edges = tmp_dual_edges;
- vert_coords = tmp_vert_coords;
- g->nv_break = num_verts;
- g->nv = num_verts;
- g->ev_break = edges;
- g->ev = edges;
- g->break_dim = num_verts + 2;
}
- }
- g->boundary = boundary;
- g->num_bounds = num_bounds;
- g->bound_inds = bound_inds;
- g->bound_verts = bound_verts;
- g->ne = num_edges;
- g->dev = dual_edges;
- g->vx = vert_coords;
+ g->bi = (uint_t *)calloc(1 + g->nb, sizeof(uint_t));
+ g->b = (uint_t *)malloc(2 * nbv * sizeof(uint_t));
- g->vei = get_verts_to_edges_ind(g->nv, g->ne, g->ev);
- g->ve =
- get_verts_to_edges(g->nv, g->ne,
- g->ev, g->vei);
+ g->bi[1] = ((int_t)nbv - (int_t)nside);
+ g->bi[g->nb] = 2 * nbv;
- g->L = L;
+ if (bound == FREE_BOUND) {
+ g->bi[2] = 2 * ((int_t)nbv - (int_t)nside);
+ g->bi[3] = 2 * (int_t)nbv - (int_t)nside;
+ }
- g->ex = get_edge_coords(
- g->ne, g->vx, g->ev);
+ uint_t n_ins0 = 0;
+ uint_t n_ins1 = 0;
+
+ g->nbi = (uint_t *)malloc(((int_t)g->nv - (int_t)g->bi[g->nb]) * sizeof(uint_t));
+ g->bni = (uint_t *)malloc((g->nv + 1) * sizeof(uint_t));
+ g->bq = (bool *)calloc(g->nv, sizeof(bool));
+ uint_t n_tmp = 0;
+
+ for (uint_t i = 0; i < f->nv; i++) {
+ g->bni[i] = n_tmp;
+ if (bv[i] != 0) {
+ g->bq[i] = true;
+ g->bq[bv[i]] = true;
+ if (side[i]) {
+ g->b[g->bi[2] + n_ins1] = i;
+ g->b[g->bi[3] + n_ins1] = bv[i];
+ n_ins1++;
+ } else {
+ g->b[g->bi[0] + n_ins0] = i;
+ g->b[g->bi[1] + n_ins0] = bv[i];
+ n_ins0++;
+ }
+ } else {
+ g->nbi[n_tmp] = i;
+ n_tmp++;
+ }
+ }
- free(tmp_tris);
+ for (uint_t i = 0; i < g->nv - f->nv + 1; i++) {
+ g->bni[f->nv + i] = n_tmp;
+ }
- g->dvx = lattice;
- g->dnv = num;
+ free(bv);
+ free(dbv);
+ free(side);
+ }
- g->voltcurmat = gen_voltcurmat(g->ne,
- g->break_dim,
- g->ev_break, c);
+ g->vei = get_verts_to_edges_ind(g->nv, g->ne, g->ev);
+ g->ve = get_verts_to_edges(g->nv, g->ne, g->ev, g->vei);
+ g->dvei = get_verts_to_edges_ind(g->dnv, g->ne, g->dev);
+ g->dve = get_verts_to_edges(g->dnv, g->ne, g->dev, g->dvei);
- g->dvei =
- get_verts_to_edges_ind(g->dnv, g->ne,
- g->dev);
- g->dve = get_verts_to_edges(
- g->dnv, g->ne, g->dev,
- g->dvei);
+ g->voltcurmat = gen_voltcurmat(g->ne, g->nv, g->ev, c);
+ uint_t num_spanning_edges;
+ g->spanning_edges = get_spanning_edges(g->ne, g->ev, g->vx, 0.5, &num_spanning_edges);
+ g->num_spanning_edges = num_spanning_edges;
- g->spanning_edges = get_spanning_edges(g->ne, g->ev_break, g->vx, 0.5, &(g->num_spanning_edges));
+ frame_free(f);
return g;
}
-graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c) {
- bool side_bounds;
- switch (lattice) {
- case (VORONOI_LATTICE):
- return ini_voro_graph(L, bound, dual, genfunc_uniform, c);
- case (SQUARE_LATTICE):
- if (bound == EMBEDDED_BOUND) side_bounds = true;
- else side_bounds = false;
- return ini_square_network(L, bound, side_bounds, c);
- case (VORONOI_HYPERUNIFORM_LATTICE):
- return ini_voro_graph(L, bound, dual, genfunc_hyperuniform, c);
- default:
- printf("lattice type unsupported\n");
- exit(EXIT_FAILURE);
- }
-}
diff --git a/src/graph_free.c b/src/graph_free.c
index e9c55d7..ee30c99 100644
--- a/src/graph_free.c
+++ b/src/graph_free.c
@@ -3,20 +3,17 @@
void graph_free(graph_t *g, cholmod_common *c) {
free(g->ev);
- if (g->ev_break != g->ev) {
- free(g->ev_break);
- }
free(g->vei);
free(g->ve);
- free(g->bound_inds);
- free(g->bound_verts);
+ free(g->bi);
+ free(g->b);
free(g->vx);
- free(g->ex);
- free(g->dev);
free(g->dvx);
+ free(g->dev);
free(g->dvei);
free(g->dve);
- free(g->spanning_edges);
+ free(g->nbi);
+ free(g->bq);
CHOL_F(free_sparse)(&(g->voltcurmat), c);
free(g);
}
diff --git a/src/graph_genfunc.c b/src/graph_genfunc.c
index 8dcd0df..d396206 100644
--- a/src/graph_genfunc.c
+++ b/src/graph_genfunc.c
@@ -1,26 +1,18 @@
#include "fracture.h"
-double *genfunc_uniform(unsigned int L, bound_t boundary, gsl_rng *r, unsigned int *num) {
- *num = 2 * pow(L / 2, 2);
+double *genfunc_uniform(uint_t num, gsl_rng *r) {
+ double *lattice = (double *)malloc(2 * num * sizeof(double));
- double *lattice = (double *)malloc(2 * (*num) * sizeof(double));
- for (unsigned int i = 0; i < (*num); i++) {
- lattice[2*i] = gsl_ran_flat(r, 0, 1);
- lattice[2*i+1] = gsl_ran_flat(r, 0, 1);
+ for (uint_t i = 0; i < 2 * num; i++) {
+ lattice[i] = gsl_rng_uniform(r);
}
return lattice;
}
-double g(double rho, double dist) {
- return 1 - gsl_sf_exp(-M_PI * rho * dist);
-}
-
-double *genfunc_hyperuniform(unsigned int L, bound_t boundary, gsl_rng *r, unsigned int *num) {
- *num = 2 * pow(L / 2, 2);
-
- double *lattice = spheres(*num, 0.8, 0.5, 0.9, 100, 100000);
+double *genfunc_hyperuniform(uint_t num, gsl_rng *r) {
+ double *lattice = spheres(num, 0.8, 0.5, 0.9, 100, 100000);
return lattice;
}
diff --git a/src/net.c b/src/net.c
index b61c9ea..a1047e2 100644
--- a/src/net.c
+++ b/src/net.c
@@ -50,35 +50,53 @@ net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, b
net->graph = g;
net->num_broken = 0;
+
net->fuses = (bool *)calloc(g->ne, sizeof(bool));
assert(net->fuses != NULL);
+
net->thres = get_thres(g->ne, beta);
net->inf = inf;
+ net->dim = g->nv;
+
+ if (g->boundary == TORUS_BOUND) {
+ net->nep = g->ne;
+ net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
+ memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
+ } else {
+ if (vb) {
+ net->dim -= g->bi[g->nb];
+ net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
+ net->nep = 0;
+ for (uint_t i = 0; i < g->ne; i++) {
+ if (!(g->bq[g->ev[2*i]] || g->bq[g->ev[2*i+1]])) {
+ net->evp[2*net->nep] = g->bni[g->ev[2*i]];
+ net->evp[2*net->nep + 1] = g->bni[g->ev[2*i + 1]];
+ net->nep++;
+ }
+ }
+ } else {
+ net->dim += 2;
+ net->evp = (uint_t *)malloc(2 * (g->ne + g->bi[2]) * sizeof(uint_t));
+ memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
+ net->nep = g->ne + g->bi[2];
+
+ for (uint_t i = 0; i < 2; i++) {
+ for (uint_t j = 0; j < g->bi[i+1] - g->bi[i]; j++){
+ net->evp[2 * (g->ne + g->bi[i] + j)] = g->b[g->bi[i] + j];
+ net->evp[2 * (g->ne + g->bi[i] + j) + 1] = g->nv + i;
+ }
+ }
+ }
+ }
+
net->voltage_bound = vb;
net->boundary_cond = bound_set(g, vb, notch_len, c);
- net->adjacency = gen_adjacency(net, false, false, 0, c);
- net->dual_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);
+ cholmod_sparse *laplacian = gen_laplacian(net, c);
net->factor = CHOL_F(analyze)(laplacian, c);
CHOL_F(factorize)(laplacian, net->factor, c);
CHOL_F(free_sparse)(&laplacian, c);
@@ -102,18 +120,11 @@ net_t *net_copy(const net_t *net, cholmod_common *c) {
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);
+ size_t evp_size = 2 * net->nep * sizeof(uint_t);
+ net_copy->evp = (uint_t *)malloc(thres_size);
+ assert(net_copy->evp != NULL);
+ memcpy(net_copy->evp, net->evp, evp_size);
- net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c);
- net_copy->dual_adjacency = CHOL_F(copy_sparse)(net->dual_adjacency, c);
net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c);
net_copy->factor = CHOL_F(copy_factor)(net->factor, c);
@@ -124,10 +135,7 @@ 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_sparse)(&(net->dual_adjacency), c);
CHOL_F(free_factor)(&(net->factor), c);
- free(net->marks);
- free(net->dual_marks);
+ free(net->evp);
free(net);
}
diff --git a/src/net_currents.c b/src/net_currents.c
index 431818f..9bb2874 100644
--- a/src/net_currents.c
+++ b/src/net_currents.c
@@ -3,7 +3,7 @@
double *net_currents(const net_t *net, const double *voltages, cholmod_common *c) {
uint_t ne = net->graph->ne;
- uint_t dim = net->graph->break_dim;
+ uint_t dim = net->graph->nv;
cholmod_sparse *voltcurmat = net->graph->voltcurmat;
cholmod_dense *x = CHOL_F(allocate_dense)(dim, 1, dim, CHOLMOD_REAL, c);
@@ -27,6 +27,22 @@ double *net_currents(const net_t *net, const double *voltages, cholmod_common *c
}
}
+ if (net->graph->boundary == TORUS_BOUND) {
+ for (uint_t i = 0; i < net->graph->bi[1]; i++) {
+ uint_t e = net->graph->b[i];
+ uint_t v1 = net->graph->ev[2 * e];
+ uint_t v2 = net->graph->ev[2 * e + 1];
+ double v1y = net->graph->vx[2 * v1 + 1];
+ double v2y = net->graph->vx[2 * v2 + 1];
+
+ if (v1y > v2y) {
+ currents[e] += 1;
+ } else {
+ currents[e] -= 1;
+ }
+ }
+ }
+
x->x = tmp_x;
CHOL_F(free_dense)(&x, c);
CHOL_F(free_dense)(&y, c);
diff --git a/src/net_fracture.c b/src/net_fracture.c
index 71e8c4b..e7f18fc 100644
--- a/src/net_fracture.c
+++ b/src/net_fracture.c
@@ -38,7 +38,7 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
double conductivity = net_conductivity(net, voltages);
- if (conductivity < 1e-12 && net->graph->boundary == TORUS_BOUND) {
+ if (conductivity < cutoff) {
free(voltages);
free(currents);
break;
@@ -60,14 +60,6 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
free(currents);
break_edge(net, last_broke, c);
-
- if (net->num_components > 1 && net->graph->boundary == TORUS_BOUND) {
- break;
- }
-
- if (net->marks[net->graph->nv] != net->marks[net->graph->nv + 1] && net->graph->boundary != TORUS_BOUND) {
- break;
- }
}
return data;
diff --git a/src/net_voltages.c b/src/net_voltages.c
index dedf5b2..c3537a5 100644
--- a/src/net_voltages.c
+++ b/src/net_voltages.c
@@ -12,11 +12,29 @@ double *net_voltages(const net_t *net, cholmod_common *c) {
exit(EXIT_FAILURE);
}
- double *voltages = (double *)x->x;
+ double *t_voltages = (double *)x->x;
x->x = NULL;
-
CHOL_F(free_dense)(&x, c);
- return voltages;
+ graph_t *g = net->graph;
+
+ if (g->boundary == TORUS_BOUND) {
+ return t_voltages;
+ } else if (net->voltage_bound) {
+ double *voltages = (double *)malloc(g->nv * sizeof(double));
+ for (uint_t i = 0; i < g->nv - g->bi[g->nb]; i++) {
+ voltages[g->nbi[i]] = t_voltages[i];
+ }
+ for (uint_t i = 0; i < 2; i++) {
+ for (uint_t j = 0; j < g->bi[i + 1] - g->bi[i]; j++) {
+ voltages[g->b[g->bi[i] + j]] = 1 - i;
+ }
+ }
+
+ free(t_voltages);
+ return voltages;
+ } else {
+ return t_voltages;
+ }
}