Dec-04-2024, 04:06 PM
Hi everybody, we are french students in mechanics and we have a project in Python, we have to code finite element method. For now we are restricted to a simple dimension and studying the displacement of a rod divided in 3 parts, the first one has young module of k1, second of k2 and third of k1 again. We wrote this code and our proffessor told us to I quote :
"You can leave the functions you already have where they are. So you need one more function which takes x, connectivity, k_func and F as parameters and returns the displacement. Inside this function you can call other functions you already have. The idea is that from now on the calculation method will remain unchanged, but we will have to evaluate solutions for different meshes and different heterogeneous properties. So to make this easier we wrap the resolution in a function."
And we have absolutely no idea of what to do. Does someone can help our little lost souls please ?
"You can leave the functions you already have where they are. So you need one more function which takes x, connectivity, k_func and F as parameters and returns the displacement. Inside this function you can call other functions you already have. The idea is that from now on the calculation method will remain unchanged, but we will have to evaluate solutions for different meshes and different heterogeneous properties. So to make this easier we wrap the resolution in a function."
And we have absolutely no idea of what to do. Does someone can help our little lost souls please ?
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse.linalg import spsolve
from scipy.sparse import coo_array
# P1 shape functions
N_1 = lambda x: (1 - x) / 2
N_2 = lambda x: (x + 1) / 2
# 2-points quadrature positions
x1 = 1/np.sqrt(3)
x2 = -1/np.sqrt(3)
# Derivatives of P1 shape functions
dN_1 = lambda x: -1. / 2
dN_2 = lambda x: 1. / 2
def extend_dx(dx):
"""Extend dimension of element sizes array to work with local matrices"""
if dx.ndim == 1:
dx = dx[:, np.newaxis, np.newaxis]
return dx
def local_stiffness_matrix(dx,k):
"""
Make array of local stiffness matrices ∫ k(x) Ni'Nj' dx with 2-points Gauss rule
"""
#on doit mettre la dedasn le E dans chaque coef voir page 4 recto
dN = np.array([
[dN_1(x1)*dN_1(x1)*k[0] + dN_1(x2)*dN_1(x2)*k[1], dN_1(x1)*dN_2(x1)*k[0]+ dN_1(x2)*dN_2(x2)*k[1]],
[dN_1(x1)*dN_2(x1)*k[0] + dN_1(x2)*dN_2(x2)*k[1], dN_2(x1)*dN_2(x1)*k[0] + dN_2(x2)*dN_2(x2)*k[1]],
])
return 2 / dx * dN.T
def assemble(x, connectivity, local_matrix_func,k):
"""
Compute local matrices and assemble into a sparse global matrix
Parameters
----------
- x: position of nodes
- connectivity: element definition (node i, node j) pair
- local_matrix_func: function to compute local stiffness matrices
"""
N = x.size
dx = extend_dx(x[1:] - x[:-1])
# Local matrices
element_matrices = local_matrix_func(dx,k)
i = np.zeros_like(element_matrices, int)
j = np.zeros_like(element_matrices, int)
i[:] = connectivity[:, :, np.newaxis]
j[:] = connectivity[:, np.newaxis, :]
# Boundary condition u(0) = 0
element_matrices[(i == 0) | (j == 0)] = 0
element_matrices[(i == 0) & (j == 0)] = 1
# Make sparse matrix
M = coo_array((element_matrices.flatten(),
(i.flatten(), j.flatten()))).tocsc()
return M
"""""
def quad_positions(x, connectivity, N):
x = np.linspace(0, L, N)
connectivity = np.zeros((N - 1, 2), int)
connectivity[:, 0] = np.arange(x - 1)
connectivity[:, 1] = np.arange(x - 1) + 1
xq_tab = np.zeros(2, N-1)
for i in range(N-1):
xq_tab[0,i] = x[i] * N_1(x1) + x[i + 1] * N_2(x1)#calcul du point de quadrature
xq_tab[1,i] = x[i] * N_1(x2) + x[i + 1] * N_2(x2)
return xq_tab
x_moy = (x[:-1] + x[1:]) / 2 #on regarde à la moitié de l'elt
k = k_func(x_q1,x_q2)
k = np.zeros(N - 1)#initialisation de k plein de zéros
pareil mais pour le tableau
tableau avec les points de quadrature qui se remplie au dir et a mesure de la boucle
k[i] = np.mean(k_func(k_tab))#calcul de k en faisant la moyenne des k1 et k2 Ã chaque point de quadrature
k = k_func(x_q1,x_q2) # on regarde la valeur de k à cet endroit
"""
def make_matrices(N):
"""
Make mass and stiffness matrices for a domain of (-pi, pi) and unit
coefficients.
"""
# Domain
x = np.linspace(0, L, N)
# Simple 1D connectivity [(0, 1), (1, 2), ...]
connectivity = np.zeros((N - 1, 2), int)
connectivity[:, 0] = np.arange(N - 1)
connectivity[:, 1] = np.arange(N - 1) + 1
# Make matrices
xq_tab = quad_positions(x, connectivity, N_elements)
k_f = k_func(xq_tab)
K = assemble(x, connectivity, local_stiffness_matrix,k_f)
return x, K, connectivity # On retourne maintenant connectivity
def quad_positions(x, connectivity, N):
"""
on range nos points de quadrature dans un tableau avec N colonne et 2 lignes pour xq1 et xq2
"""
xq_tab = np.zeros((2, N - 1))
for i in range(N - 1):
xq_tab[0, i] = x[connectivity[i, 0]] * N_1(x1) + x[connectivity[i, 1]] * N_2(x1)
xq_tab[1, i] = x[connectivity[i, 0]] * N_1(x2) + x[connectivity[i, 1]] * N_2(x2)
return xq_tab
# on def les k pour les 3 zones de la barre avec k1=k3
def k_func(x):
k = np.zeros_like(x)
k[x < 1/3] = k1
k[(x >= 1/3) & (x <= 2/3)] = k2
k[x > 2/3] = k1
return k
"""
def E_func(x):
return 1+np.cos(2*np.pi*x)**2
"""
N_elements = 9
L = 1
k1 = 10
k2 = 1
F = np.zeros(N_elements)
F0=1
F[-1] = F0 # assemblage de F (zeros sauf dernière composante)
x, K, connectivity = make_matrices(N_elements) # je sors connectivity pour pouvoir l'utiliser après
"""""
#plot la variation de k le long de la barre
plt.figure(figsize=(8, 4))
x_points = (x[:-1] + x[1:]) / 2
plt.plot(x_points, k_func(x_points), label="variation de l'hétéro de k")
plt.xlabel('x')
plt.ylabel('k(x)')
plt.title("avec hétéro")
plt.legend()
plt.grid(True)
"""
# Solution Ku = F
u = spsolve(K, F)
x_ana=np.linspace(0,L,N_elements)
def u_ana(x):#on définit la solution analytique
u1=F0/k1*x
u2=F0/k2*x+F0*L/(3*k1*k2)*(k2-k1)
u3=u1+2*F0*L/(3*k2)+F0*L/(3*k1*k2)*(k2-k1)-2*F0*L/(3*k1)
uana=np.zeros_like(x)
uana[x < 1/3] = u1[x < 1/3]
uana[(x >= 1/3) & (x <= 2/3)] = u2[(x >= 1/3) & (x <= 2/3)]
uana[x > 2/3] = u3[x > 2/3]
return uana
plt.figure(figsize=(8, 5))
plt.plot(x_ana, u_ana(x_ana), label='Solution analytique u(x)', marker='')
plt.plot(x, u, label='Solution u(x)', marker='')
plt.xlabel('x')
plt.ylabel('u(x)')
plt.title('Solution Ku = F hetero')
plt.legend()
plt.grid(True)
plt.show()
#test pour l'erreur pas sûr du tout
import scipy.integrate as integration #ref :https://docs.scipy.org/doc/scipy/tutorial/integrate.html
import scipy.interpolate as interpolation #ref: https://docs.scipy.org/doc/scipy/reference/interpolate.html
e = u_ana(x_ana)- u
# Interpolation polynomiale de e
p = 1
#help(np.polyfit) #degré du polynôme qui va approximer e
coef = np.polyfit(x, e, p)
poly = np.poly1d(coef)
e_deri = np.polyder(poly)
erreur = e_deri(x)**2 + e**2 #ça c'est l'expression dans l'intégrale de la formule de l'erreur énergétique
integrale = np.trapz(erreur, x) #grace à scipy on peut intégrer avec la méthode des trapèzes comme ça
e_norme = np.sqrt(integrale)#norme énergétique de e
print(f"Norme énergétique : {e_norme}")
plt.figure(figsize=(8, 5))
plt.plot(x, e, label="Erreur e(x)")
plt.plot(x, e_deri(x), label="Dérivée de l'erreur e'(x)")
plt.xlabel("x")
plt.ylabel("Valeurs")
plt.title("Erreur et dérivée de l'erreur")
plt.legend()
plt.grid(True)
plt.show()
# Ã faire
#pour chaque élémetn de calculker la vraie position des points de quadrature FAIT
# implementer une mesure de cette erreur : comment l'erreur évolue avec le nombre d'elt
# vecteur force
#on regarde pas au niveau du noeud mais dans l'élément
)