In [1]:
import numpy as np
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve
from scipy.spatial import Delaunay
import matplotlib.pyplot as plt
import sympy as sp
In [18]:
class FEMSolver:
def __init__(self, dimension, domain, f_expr, bc, mesh_density=10, orientation="null", T_max=5.0, dt=0.01):
"""
Initialise le solveur FEM.
:param dimension: Dimension du problème (1 pour 1D, 2 pour 2D).
:param domain: Domaine de résolution.
- En 1D : tuple (a, b) où [a, b] est l'intervalle.
- En 2D : tuple (xmin, xmax, ymin, ymax).
:param f_expr: Expression sympy de f(x[, y]) pour le second membre.
:param bc: Conditions aux bords.
- En 1D : dict {x_val: u_val} pour imposer u(a) ou u(b).
- En 2D : dict {(x, y): u_val} pour les points sur le bord.
:param mesh_density: Densité du maillage.
:param orientation: Orientation des diagonales pour le maillage 2D ("left", "right", "null").
:param T_max: Temps maximum pour la simulation.
:param dt: Pas de temps pour la discrétisation temporelle.
"""
self.dimension = dimension
self.domain = domain
self.f_expr = f_expr
self.bc = bc
self.mesh_density = mesh_density
self.orientation = orientation
self.T_max = T_max
self.dt = dt
self.points = None
self.simplices = None
self.solution = None
def generate_mesh(self):
"""
Génère un maillage adapté à la dimension.
"""
if self.dimension == 1:
a, b = self.domain
self.points = np.linspace(a, b, self.mesh_density)
elif self.dimension == 2:
xmin, xmax, ymin, ymax = self.domain
x = np.linspace(xmin, xmax, self.mesh_density)
y = np.linspace(ymin, ymax, self.mesh_density)
X, Y = np.meshgrid(x, y)
self.points = np.c_[X.ravel(), Y.ravel()]
simplices = []
for i in range(self.mesh_density - 1):
for j in range(self.mesh_density - 1):
p0 = i * self.mesh_density + j
p1 = i * self.mesh_density + j + 1
p2 = (i + 1) * self.mesh_density + j
p3 = (i + 1) * self.mesh_density + j + 1
if self.orientation == "left":
simplices.append([p0, p2, p3])
simplices.append([p0, p3, p1])
elif self.orientation == "right":
simplices.append([p0, p2, p1])
simplices.append([p2, p3, p1])
else:
simplices.extend(Delaunay(self.points).simplices)
self.simplices = simplices
return
self.simplices = np.array(simplices)
def assemble_system(self):
"""
Assemble le système linéaire (matrice et vecteur de charges).
"""
if self.dimension == 1:
num_points = len(self.points)
A = lil_matrix((num_points, num_points))
b = np.zeros(num_points)
h = self.points[1] - self.points[0]
f_func = sp.lambdify('x', self.f_expr, 'numpy')
for i in range(num_points - 1):
A[i, i] += 1 / h
A[i, i + 1] -= 1 / h
A[i + 1, i] -= 1 / h
A[i + 1, i + 1] += 1 / h
b[i] += h / 2 * f_func(self.points[i])
b[i + 1] += h / 2 * f_func(self.points[i + 1])
for x, u in self.bc.items():
idx = np.where(np.isclose(self.points, x))[0][0]
A[idx, :] = 0
A[idx, idx] = 1
b[idx] = u
return A, b
elif self.dimension == 2:
num_points = len(self.points)
A = lil_matrix((num_points, num_points))
b = np.zeros(num_points)
f_func = sp.lambdify(('x', 'y'), self.f_expr, 'numpy')
def compute_area(p0, p1, p2):
# Calcul de l'aire du triangle (determinant de la matrice J)
J = np.array([[p1[0] - p0[0], p2[0] - p0[0]],
[p1[1] - p0[1], p2[1] - p0[1]]])
detJ = np.abs(np.linalg.det(J)) / 2
return detJ
for tri in self.simplices:
p0, p1, p2 = self.points[tri]
# Calcul de l'aire du triangle
area = compute_area(p0, p1, p2)
# Calcul de la matrice de rigidité locale
B = np.linalg.inv(np.array([[p1[0] - p0[0], p2[0] - p0[0]],
[p1[1] - p0[1], p2[1] - p0[1]]])).T @ np.array([[-1, 1, 0], [-1, 0, 1]])
stiffness_local = (B.T @ B) * area + np.eye(3) * area / 12
# Calcul de la force locale
f_avg = f_func(*np.mean(self.points[tri], axis=0))
b_local = f_avg * area / 6
# Assemblage dans la matrice globale
for i in range(3):
for j in range(3):
A[tri[i], tri[j]] += stiffness_local[i, j]
b[tri[i]] += b_local
# Application des conditions aux bords
for (x, y), u in self.bc.items():
x = float(x)
y = float(y)
idx = np.where((np.isclose(self.points[:, 0], x)) & (np.isclose(self.points[:, 1], y)))[0]
if idx.size > 0:
idx = idx[0]
A[idx, :] = 0
A[idx, idx] = 1
b[idx] = u
return A, b
def solve(self):
"""
Résout le système linéaire pour obtenir la solution.
"""
self.generate_mesh()
A, b = self.assemble_system()
self.solution = spsolve(A.tocsr(), b)
def time_step_solution(self):
"""
Simule l'évolution de la température avec une discrétisation temporelle.
Utilisation de la méthode explicite d'Euler pour l'équation de la chaleur 1D.
"""
num_points = len(self.points)
alpha = 0.01 # Diffusivité thermique, à ajuster
u = np.zeros(num_points)
# Conditions initiales (f_expr à t=0)
f_func = sp.lambdify('x', self.f_expr, 'numpy')
u = f_func(self.points)
times = np.arange(0, self.T_max, self.dt)
# Boucle temporelle
u_values = []
for t in times:
u_values.append(u)
# Schéma explicite d'Euler pour la diffusion de la chaleur
u_new = np.copy(u)
for i in range(1, num_points - 1):
u_new[i] = u[i] + alpha * self.dt / (self.points[i+1] - self.points[i-1])**2 * (u[i-1] - 2 * u[i] + u[i+1])
u = u_new
return times, np.array(u_values)
def plot_solution(self):
"""
Visualise la solution pour différents instants.
"""
if self.dimension == 1:
times, u_values = self.time_step_solution()
# plt.figure(figsize=(10, 6))
for i, t in enumerate(times[::int(len(times)/5)]): # Affiche la solution à différents t
plt.plot(self.points, u_values[i], label=f"t={t:.2f}")
plt.title(r"Chaleur: $\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}$")
plt.xlabel("Position (x)")
plt.ylabel("Température (u)")
plt.legend()
plt.grid(True)
plt.show()
elif self.dimension == 2:
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(self.points[:, 0], self.points[:, 1], self.simplices, self.solution, cmap='viridis', edgecolor='k', linewidth=0.2)
ax.set_title(r"Chaleur: $\frac{\partial u}{\partial t} = \alpha \nabla^2 u$")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("u(x, y)")
plt.show()
In [16]:
N = 50
x, y = sp.symbols('x y')
f_2d = sp.sin(sp.pi * x) * sp.sin(sp.pi * y)
x_vals = np.linspace(0, 1, N)
y_vals = np.linspace(0, 1, N)
bc = {}
for y in y_vals:
bc[(0.0, y)] = 0
bc[(1.0, y)] = 0
for x in x_vals:
bc[(x, 0.0)] = 0
bc[(x, 1.0)] = 0
solver_2d = FEMSolver(dimension=2, domain=(0, 1, 0, 1), f_expr=f_2d, bc=bc, mesh_density=N, orientation="left")
solver_2d.solve()
solver_2d.plot_solution()
In [ ]: