# Description:
#   Generate a mesh triangulation of the unit disc.
#
# Arguments:
#   N       Number of nodes in the mesh.
#
# Returns:
#   p		Nodal points, (x,y)-coordinates for point i given in row i.
#   tri   	Elements. Index to the three corners of element i given in row i.
#   edge  	Edge lines. Index list to the two corners of edge line i given in row i.
#
#   Author: Kjetil A. Johannessen, Abdullah Abdulhaque
#   Last edit: 09-10-2019


import numpy as np
import scipy.spatial as spsa


def GetDisc(N):
	# Controlling the input.
	if N < 4:
		print("Error. N >= 4 reguired for input.")
		return

	# Defining auxiliary variables.
	M,r,alpha,theta = CircleData(N)

	# Generating the nodal points.
	p = NodalPoints(M,N,alpha,theta,r)

	# Generating the elements.
	mesh = spsa.Delaunay(p)
	tri = mesh.simplices

	# Generating the boundary elements.
	edge = FreeBoundary(N,alpha)

	return p,tri,edge

def NodalPoints(M,N,alpha,theta,r):
	# Auxiliary function for generating nodal points.
	p = np.zeros((N,2))
	k = 1
	for i in range(1,M+1):
		t = theta[i]
		for j in range(0,alpha[i]):
			p[k,:] = [np.cos(t)*r[i],np.sin(t)*r[i]]
			t += 2*np.pi/alpha[i]
			k += 1

	return p


def FreeBoundary(N,alpha):
	# Auxiliary function for generating boundary nodes.
	E = np.arange(N-alpha[-1]+1,N+1)
	edge = np.zeros((len(E),2),dtype=np.int)
	for i in range(0,len(E)):
		edge[i,:] = [E[i],E[i]+1]
	edge[-1,-1] = N-alpha[-1]+1
	edge -= 1

	return edge


def CircleData(N):
	# Number of outward circles,excluding the origin.
	M = np.int(np.floor(np.sqrt(N/np.pi)))

	# Radius of the different circles.
	r = np.linspace(0,1,M+1)

	# Number of DOF in each circle.
	alpha_temp = np.floor((2*np.pi*M)*r)
	alpha = np.zeros(len(alpha_temp),dtype=np.int)
	for i in range(0,len(alpha_temp)):
		alpha[i] = np.int(alpha_temp[i])

	# Fine-tuning to get the right amount of DOF.
	alpha[0] = 1
	i = 1
	while sum(alpha) > N:
		if alpha[i] > 0:
			alpha[i] -= 1
		i += 1
		if sum(alpha[1:M]) == 0:
			i = M
		elif i > M:
			i = 1
	while sum(alpha) < N:
		alpha[-1] += 1

	# Creating the starting angle.
	theta = np.pi/alpha
	theta[0:len(alpha):2] = 0

	return M,r,alpha,theta