diff options
Diffstat (limited to 'src/frame_create.c')
| -rw-r--r-- | src/frame_create.c | 167 |
1 files changed, 167 insertions, 0 deletions
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); +} + |