# Description: # Generate a mesh triangulation on the L-shape domain. # # Arguments: # N Number of nodes on each local patch. # # 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: Abdullah Abdulhaque # Last edit: 28-10-2019 import numpy as np import scipy.spatial as spsa def GetLShape(N): # Controlling the input. if N < 2: raise RuntimeError("GetLShape","Error. N >= 2 reguired.") # Generating nodal points. M = 2*N-1 c = np.linspace(-1,1,M) p =[] for i in c: for j in c: if (i <= 0) or ((i > 0) and (j <= 0)): p += [[i,j]] p = np.array(p) # Generating elements. mesh = spsa.Delaunay(p) tri_temp = mesh.simplices tri = [] for i in tri_temp: pts = p[i,:] if not np.all(np.mean(p[i,:],axis=0) > 0): tri.append(i) tri = np.array(tri) # Generating edges. edge = [] for i in range(1,M): edge += [[i,i+1]] for i in range(1,N): edge += [[i*M,(i+1)*M]] for i in range(1,N): edge += [[N*M-i+1,N*M-i]] edge += [[(N-1)*M+N,N*M+N]] for i in range(2,N): edge += [[N*M+(i-1)*N,N*M+i*N]] for i in range(1,N): edge += [[M*N+N*(N-1)-i+1,M*N+N*(N-1)-i]] for i in range(2,N): edge += [[M*N+N*(N-i)+1,M*N+N*(N-i-1)+1]] for i in range(1,N+2): edge += [[M*(N-i+1)+1,M*(N-i)+1]] edge = np.array(edge) edge -= 1 return p,tri,edge