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.c299
1 files changed, 0 insertions, 299 deletions
diff --git a/src/graph_create.c b/src/graph_create.c
deleted file mode 100644
index 09a3aed..0000000
--- a/src/graph_create.c
+++ /dev/null
@@ -1,299 +0,0 @@
-
-#include "fracture.h"
-
-uint_t *get_spanning_edges(uint_t num_edges, uint_t *edges_to_verts, double *vert_coords, double cut, uint_t *n) {
- uint_t *spanning_edges = (uint_t *)malloc(num_edges * sizeof(uint_t));
- (*n) = 0;
- for (uint_t i = 0; i < num_edges; i++) {
- uint_t v1, v2;
- v1 = edges_to_verts[2 * i];
- v2 = edges_to_verts[2 * i + 1];
- double v1y, v2y;
- v1y = vert_coords[2 * v1 + 1];
- v2y = vert_coords[2 * v2 + 1];
- if ((fabs(v1y - v2y) < 0.5) && ((v1y <= cut && v2y > cut) || (v1y >= cut && v2y < cut))) {
- spanning_edges[*n] = i;
- (*n)++;
- }
- }
- return spanning_edges;
-}
-
-double *get_edge_coords(uint_t num_edges, double *vert_coords,
- uint_t *edges_to_verts) {
- double *output = (double *)malloc(2 * num_edges * sizeof(double));
-
- for (uint_t i = 0; i < num_edges; i++) {
- uint_t v1, v2;
- double v1x, v1y, v2x, v2y, dx, dy;
- v1 = edges_to_verts[2 * i];
- v2 = edges_to_verts[2 * i + 1];
- output[2 * i] = 0;
- output[2 * i + 1] = 0;
- 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(dx) > 0.5) {
- if (dx > 0) {
- v2x = v2x + 1;
- } else {
- v1x = v1x + 1;
- }
- }
- if (fabs(dy) > 0.5) {
- if (dy > 0) {
- v2y = v2y + 1;
- } else {
- v1y = v1y + 1;
- }
- }
- output[2 * i] = (v1x + v2x) / 2;
- output[2 * i + 1] = (v1y + v2y) / 2;
- }
-
- return output;
-}
-
-uint_t *get_verts_to_edges_ind(uint_t num_verts,
- uint_t num_edges,
- const uint_t *edges_to_verts) {
- uint_t *output =
- (uint_t *)calloc(num_verts + 1, sizeof(uint_t));
- assert(output != NULL);
-
- for (uint_t i = 0; i < 2 * num_edges; i++) {
- if (edges_to_verts[i] < num_verts) {
- output[edges_to_verts[i] + 1]++;
- }
- }
-
- for (uint_t i = 0; i < num_verts; i++) {
- output[i + 1] += output[i];
- }
-
- return output;
-}
-
-uint_t *get_verts_to_edges(uint_t num_verts, uint_t num_edges,
- const uint_t *edges_to_verts,
- const uint_t *verts_to_edges_ind) {
- uint_t *output = (uint_t *)calloc(verts_to_edges_ind[num_verts],
- sizeof(uint_t));
- uint_t *counts =
- (uint_t *)calloc(num_verts, sizeof(uint_t));
- for (int i = 0; i < 2 * num_edges; i++) {
- if (edges_to_verts[i] < num_verts) {
- output[verts_to_edges_ind[edges_to_verts[i]] +
- counts[edges_to_verts[i]]] = i / 2;
- counts[edges_to_verts[i]]++;
- }
- }
-
- free(counts);
-
- return output;
-}
-
-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++;
- }
- }
- }
-
- return nce;
-}
-
-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);
-
- 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;
-
- 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);
- }
-
- 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;
- }
-
- 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;
-
- if (bv[g->ev[2 * ce[i] + vin]] == 0) {
- nbv++;
- if (cin == 0) {
- nside++;
- side[g->ev[2 * ce[i] + vin]] = true;
- }
-
- bv[g->ev[2 * ce[i] + vin]] = g->nv;
-
- 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;
-
- g->nv++;
- } else {
- g->ev[2 * ce[i] + vin] = bv[g->ev[2 * ce[i] + vin]];
- }
- 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 {
- 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]];
- }
- }
-
- g->bi = (uint_t *)calloc(1 + g->nb, sizeof(uint_t));
- g->b = (uint_t *)malloc(2 * nbv * sizeof(uint_t));
-
- g->bi[1] = ((int_t)nbv - (int_t)nside);
- g->bi[g->nb] = 2 * nbv;
-
- if (bound == FREE_BOUND) {
- g->bi[2] = 2 * ((int_t)nbv - (int_t)nside);
- g->bi[3] = 2 * (int_t)nbv - (int_t)nside;
- }
-
- 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++;
- }
- }
-
- for (uint_t i = 0; i < g->nv - f->nv + 1; i++) {
- g->bni[f->nv + i] = n_tmp;
- }
-
- free(bv);
- free(dbv);
- free(side);
- }
-
- 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->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;
-
- frame_free(f);
-
- return g;
-}
-
-