# Description:
#   Generate a mesh triangulation of the unit disc.
#
# Arguments:
#   N       Number of nodes in the mesh.
#   theta   Angle of the slice.
#
# 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: 07-10-2019


import numpy as np
import scipy.spatial as spsa


def GetSlice(N,theta):
    # Controlling the input.
    if N < 3:
        raise RuntimeError("GetSlice","Error. N >= 3 reguired.")
    elif (theta <= 0) or (theta >= 2*np.pi):
        raise RuntimeError("GetSlice","Error. 0 < theta < pi reguired for input.")

    # Approximating design.
    M = np.int(np.floor(np.sqrt(2*N/theta)))
    r = np.linspace(0,1,M+1)
    al = np.floor((theta*M)*r)

    # Fine-tuning to get the right amount of DOF.
    al[0] = 1
    i = M
    while sum(al) > N:
        if al[i] > 0:
            al[i] -= 1
        i -= 1
        if i < 1:
            i = M
    while sum(al) < N:
        al[-1] += 1
    alpha = np.zeros(len(al),dtype=np.int)
    for i in range(0,len(al)):
        alpha[i] = np.int(al[i])

    # Special case, small example.
    if (M == 1) or ((M == 2) and (alpha[1] < 3)):
        i = np.arange(0,N-1)
        t = (theta/(N-2))*i
        p = np.zeros((len(i)+1,2))
        for j in range(1,len(i)+1):
            p[j,:] = [np.cos(t[j-1]),np.sin(t[j-1])]
        E = np.array([np.arange(2,N),np.arange(3,N+1)])
        tri = np.zeros((N-2,3),dtype=np.int)
        edge = np.zeros((N,2),dtype=np.int)
        for j in range(0,N-2):
            tri[j,0] = 1
            tri[j,1] = E[0,j]
            tri[j,2] = E[1,j]
            edge[j,0] = E[0,j]
            edge[j,1] = E[1,j]
        edge[N-2,0] = N
        edge[N-2,1] = 1
        edge[N-1,0] = 1
        edge[N-1,1] = 2

    # General case, large example.
    else:
        p = np.zeros((N,2))
        edge = []
        k = 1
        for i in range(1,M+1):
            t = 0
            for j in range(0,alpha[i]):
                p[k,:] = [np.cos(t)*r[i],np.sin(t)*r[i]]
                if alpha[i]-1 != 0:
                    t += theta/(alpha[i]-1)
                k += 1
            edge += [[k,k-alpha[i]]]
            edge += [[k-alpha[i]-alpha[i-1]+1,k-alpha[i]+1]]
        k = N-alpha[-1]+1
        for i in range(k,N):
            edge += [[i,i+1]]
        edge = np.array(edge)
        mesh = spsa.Delaunay(p)
        tri = mesh.simplices

    return p,tri,edge