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