fixed voltage and torus conditions, current and free boundaries and broken right now

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);
+}
+