File contents
import numpy as np
import matplotlib.pyplot as plt
def Q_mesh(Nx,Ny):
# Crear nodos en el intervalo [0,1]
x = np.linspace(0, 1, Nx + 1)
y = np.linspace(0, 1, Ny + 1)
nodes = np.array([[xi,yi] for xi in x for yi in y])
# Crear elementos finitos (cuadriláteros)
elements = []
for i in range(Nx):
for j in range(Ny):
node1 = i * (Ny + 1) + j
node2 = (i + 1) * (Ny + 1) + j
node3 = (i + 1) * (Ny + 1) + (j + 1)
node4 = i * (Ny + 1) + (j + 1)
elements.append([node1, node2, node3, node4])
elements=np.array(elements)
# Graficar elementos finitos
plt.figure(figsize=(4,4))
for element in elements:
x_coords = nodes[element, 0]
y_coords = nodes[element, 1]
plt.plot(np.append(x_coords, x_coords[0]), np.append(y_coords, y_coords[0]), 'b-')
plt.title('Malla Rectangular Uniforme')
plt.show()
return nodes,elements,x,y
def shishkin_mesh(Nx,Ny,epsilon):
N=Nx
l=epsilon*np.log(N)
# Crear nodos
x0=np.linspace(0,l,int(N/4)+1)
x1=np.linspace(l,1-l,int(N/2)+1)
x1=x1[1:]
x=np.append(x0,x1)
x2=np.linspace(1-l,1,int(N/4)+1)
x2=x2[1:]
x=np.append(x,x2)
y0 = np.linspace(0,l,int(N/4)+1)
y1 = np.linspace(l,1-l,int(N/2)+1)
y1=y1[1:]
y=np.append(y0,y1)
y2 = np.linspace(1-l,1,int(N/4)+1)
y2=y2[1:]
y=np.append(y,y2)
nodes = np.array([[xi, yi] for xi in x for yi in y])
nx=len(x)-1
ny=len(y)-1
# Crear elementos finitos (cuadriláteros)
elements = []
for i in range(nx):
for j in range(ny):
node1 = i * (ny + 1) + j
node2 = (i + 1) * (ny + 1) + j
node3 = (i + 1) * (ny + 1) + (j + 1)
node4 = i * (ny + 1) + (j + 1)
elements.append([node1, node2, node3, node4])
elements=np.array(elements)
# Graficar elementos finitos
plt.figure(figsize=(4,4))
for element in elements:
x_coords = nodes[element, 0]
y_coords = nodes[element, 1]
plt.plot(np.append(x_coords, x_coords[0]), np.append(y_coords, y_coords[0]), 'b-')
plt.title('Malla Rectangular de Shishkin')
plt.show()
return nodes,elements,x,y
def boundary(nodos): #funcion para identificar los nodos del borde
borde=np.array([])
j=0
for nodo in nodos:
if nodo[0]==0 or nodo[1]==0 or nodo[0]==1 or nodo[1]==1:
borde=np.append(borde,j)
j=j+1
return borde
Click here to get the file
