From 1e1fdfc2e3892667bccaf317a01defd8832041c7 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 16 Jan 2017 01:31:10 -0500 Subject: fixed voltage and torus conditions, current and free boundaries and broken right now --- src/graph_create.c | 717 +++++++++++++---------------------------------------- 1 file changed, 167 insertions(+), 550 deletions(-) (limited to 'src/graph_create.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); - } -} -- cgit v1.2.3-70-g09d2