# 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: October 2019 import numpy as np import scipy.spatial as spsa def getDisc(N): # Defining auxiliary variables. if N < 4: print("Error. N > 3 reguired for input.") return M = np.floor(np.sqrt(N/np.pi)) # Number of circles outward, (not counting origin). M = np.int(M) r = np.linspace(0,1,M+1) # Radius of the different circles. alpha_temp = np.floor((2*np.pi*M)*r) # Number of DOF in each circle. 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 # Generating the 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][0] = np.cos(t)*r[i] p[k][1] = np.sin(t)*r[i] t += 2*np.pi/alpha[i] k += 1 # Generating the elements. mesh = spsa.Delaunay(p) tri = mesh.simplices # Generating the boundary elements. 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,0] = E[i] edge[i,1] = E[i]+1 edge[-1,-1] = N-alpha[-1]+1 return p,tri,edge def freeBoundary(mesh): # Auxiliary function for generating boundary nodes. ans = [] for ind, neigh in zip(mesh.simplices, mesh.neighbors): for j in range(3): if neigh[j] == -1: ans += [[ind[j-1],ind[j-2]]] return np.array(ans)