summaryrefslogtreecommitdiff
path: root/src/frame_create.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/frame_create.c')
-rw-r--r--src/frame_create.c167
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);
+}
+