summaryrefslogtreecommitdiff
path: root/src/graph_create.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/graph_create.c')
-rw-r--r--src/graph_create.c717
1 files changed, 167 insertions, 550 deletions
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);
- }
-}