Équation de Poisson 2D : méthode des éléments finis

1. Introduction

2. Méthode des éléments finis

2.a. Formulation variationnelle

2.b. Conditions limites

2.c. Maillage triangulaire - Fonctions affines par morceaux

2.d. Équations algébriques

2.e. Matrice du système

1. Introduction

Ce document traite de la résolution de l'équation de Poisson bidimensionnelle par la méthode des éléments finis. La page Équation de Poisson 1D : méthode des éléments finis développe le cas unidimensionel.

Soit $Ω$ un domaine de $ℝ^{2}$ de frontière $\partial Ω$ . Soient deux fonctions $a : Ω \to ℝ$ et $s : Ω \to ℝ$ . On recherche une fonction $u : Ω \to ℝ$ deux fois dérivable vérifiant l'équation différentielle

d i v (a \vec{g r a d} u) = - s (1)

Par exemple dans un problème de diffusion thermique, u désigne la température, a est la conductivité thermique et s est la source thermique. Dans un problème d'électrostatique u est le potentiel, a la permittivité électrique et s la densité volumique de charge. Cette équation aux dérivées partielles est plus générale que l'équation de Poisson mais se ramène à celle-ci dans un milieu homogène, c'est-à-dire lorsque a est uniforme :

\nabla^{2} u = - \frac{s}{a} (2)

On se place dans le cas de coordonnées cartésiennes, qui correspond à un problème invariant par translation dans la direction Z. On a donc :

\vec{g r a d} u = \frac{\partial u}{\partial x} \vec{u_{x}} + \frac{\partial u}{\partial y} \vec{u_{y}} (3)

d i v \vec{v} = \frac{\partial v_{x}}{\partial x} + \frac{\partial v_{y}}{\partial y} (4)

Pour que l'équation $(1)$ admette une unique solution, il faut définir des conditions limites sur $\partial Ω$ .

2. Méthode des éléments finis

2.a. Formulation variationnelle

Pour une fonction $f : Ω \to ℝ$ de classe C¹, on a :

\int_{Ω} \frac{\partial f}{\partial x} d σ = \int_{\partial Ω} f \vec{n} \cdot \vec{u_{x}} d ℓ (5)

où $\vec{n}$ désigne le vecteur unitaire orthogonal à $\partial Ω$ et dirigé dans le sens sortant par rapport à $Ω$ . La même égalité est évidemment valable pour la coordonnée y. Appliquons ce théorème à :

f = v a \frac{\partial u}{\partial x} (6)

On obtient :

\int_{Ω} a \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} d σ + \int_{Ω} v \frac{\partial}{\partial x} (a \frac{\partial u}{\partial x}) d σ = \int_{\partial Ω} v a \frac{\partial u}{\partial x} \vec{n} \cdot \vec{u_{x}} d ℓ (7)

et de même :

\int_{Ω} a \frac{\partial u}{\partial y} \frac{\partial v}{\partial y} d σ + \int_{Ω} v \frac{\partial}{\partial x} (a \frac{\partial u}{\partial y}) d σ = \int_{\partial Ω} v a \frac{\partial u}{\partial y} \vec{n} \cdot \vec{u_{y}} d ℓ (8)

La somme de ces deux égalités conduit à :

\int_{Ω} a (\vec{g r a d} u) \cdot (\vec{g r a d} v) d σ + \int_{Ω} v d i v (a \vec{g r a d} u) d σ = \int_{\partial Ω} v a (\vec{g r a d} u) \cdot \vec{n} d ℓ (9)

qui devient, en utilisant l'équation $(1)$ :

\int_{Ω} a (\vec{g r a d} u) \cdot (\vec{g r a d} v) d σ = \int_{Ω} v s d σ + \int_{\partial Ω} v a (\vec{g r a d} u) \cdot \vec{n} d ℓ (10)

Cette égalité est écrite dans le cas bidimensionnel, où $Ω$ est une surface plane et $\partial Ω$ le contour de cette surface. $d σ$ est donc l'élément de surface et $d ℓ$ l'élément du contour. Le même résultat est bien sûr valable en tridimensionnel ( $Ω$ est alors un volume et $\partial Ω$ la surface qui le délimite).

Soit V l'ensemble des fonctions $Ω : \to ℝ$ de carré sommable et dont le gradient est aussi de carré sommable sur $Ω$ .

La formulation variationnelle du problème (la recherche de u) s'énonce ainsi [1] : on recherche une fonction $u \in V$ telle que pour toute fonction $v \in V$ on ait l'égalité $(10)$ . Cet énoncé est souvent appelé forme faible du problème de Poisson car les exigences sur la fonction u sont moindres que pour l'équation $(1)$ (elle n'est plus nécessairement deux fois dérivable).

2.b. Conditions limites

Séparons la frontière du domaine en deux parties : la partie $\partial Ω_{d}$ sur laquelle une condition limite de Dirichlet est imposée, la partie $\partial Ω_{n}$ sur laquelle une condition limite de Neumann est imposée. L'égalité $(10)$ s'écrit :

\int_{Ω} a (\vec{g r a d} u) \cdot (\vec{g r a d} v) d σ = \int_{Ω} v s d σ + \int_{\partial Ω_{d}} v a (\vec{g r a d} u) \cdot \vec{n} d ℓ + \int_{\partial Ω_{n}} v a (\vec{g r a d} u) \cdot \vec{n} d ℓ (11)

Comme il est montré dans Équation de Poisson 1D : méthode des éléments finis, la formulation variationnelle doit être modifiée en considérant les fonctions de V qui s'annulent sur $\partial Ω_{d}$ . Notons V₀ l'ensemble de ces fonctions. La formulation variationnelle devient : obtenir une fonction u telle que pour toute fonction $v \in V_{0}$ :

\int_{Ω} a (\vec{g r a d} u) \cdot (\vec{g r a d} v) d σ = \int_{Ω} v s d σ + \int_{\partial Ω_{n}} v a (\vec{g r a d} u) \cdot \vec{n} d ℓ (12)

La condition de Neumann consiste à fixer la valeur de $(\vec{g r a d} u) \cdot \vec{n}$ sur $\partial Ω_{n}$ . Une variante de la condition de Neumann est la condition de Robin, qui consiste à écrire sur la frontière :

(\vec{g r a d} u) \cdot \vec{n} = K_{0} (u - u_{0}) + g_{0} (13)

où K₀, u₀ et g₀ sont des constantes. Pour K₀=0, on retrouve la condition de Neumann avec g₀ la composante normale du gradient sur la frontière. Plus généralement, ces trois constantes peuvent dépendre du point de la frontière (auquel cas il ne s'agit plus de constantes). Dans ce document, on nommera condition de Neumann la condition définie par l'équation $(13)$ .

Finalement, l'égalité de la formulation variationnelle s'écrit :

\int_{Ω} a (\vec{g r a d} u) \cdot (\vec{g r a d} v) d σ = \int_{Ω} v s d σ + \int_{\partial Ω_{n}} v a (K_{0} (u - u_{0}) + g_{0}) d ℓ (14)

2.c. Maillage triangulaire - Fonctions affines par morceaux

Le domaine $Ω$ est subdivisé en triangles : chaque maille est un triangle et comporte donc trois nœuds. Une maille partage avec une maille voisine deux nœuds. La figure suivante montre un maillage triangulaire très simple pour un domaine rectangulaire (celui que nous utiliserons pour les tests) :

Figure pleine page

Le maillage comporte P nœuds. Le nœud N_i a pour coordonnées (x_i,y_i). En général, les nœuds ne sont évidemment pas disposés sur une grille.

Les fonctions u et v du problème variationnel sont considérées dans l'ensemble, noté V_h, des fonctions continues et affines par morceaux sur le maillage (l'indice h désigne le maillage). Pour $u_{h} \in V_{h}$ on a sur chaque maille : $u_{h} (x, y) = c_{0} x + c_{1} y + c_{2}$ . Si les mailles sont assez petites, les valeurs de u_h aux nœuds sont censées être de bonnes approximations de u, la solution de l'équation $(1)$ avec les conditions limites. La fonction u_h est une solution du problème variationnel mais elle n'est évidemment pas une solution de $(1)$ .

Si l'on note $c_{0}, c_{1}, c_{2}$ et $c'_{0}, c'_{1}, c'_{2}$ les coefficients de u_h pour deux mailles voisines, on a sur le côté commun à ces deux mailles, par continuité :

c_{0} x + c_{1} y + c_{2} = c'_{0} x + c'_{1} y + c'_{2} (15)

ce qui est l'équation de la droite passant par les deux nœuds communs. On a bien évidemment :

u_{h} (N_{i}) = c_{0} x_{i} + c_{1} y_{i} + c_{2} (16)

où $c_{0}, c_{1}, c_{2}$ sont les coefficients d'un des trois triangles qui partagent le nœud N_i.

On définit les fonctions affines par morceaux suivantes, pour $i = 1, \dots P$ :

\begin{matrix} ϕ_{i} (N_{j}) = 1, j = i & (17) \\ ϕ_{i} (N_{j}) = 0, j \neq i & (18) \end{matrix}

Autrement dit, la fonction $ϕ_{i}$ vaut 1 sur le nœud N_i et 0 sur tous les autres nœuds. Ces fonctions constituent une base V_h. Pour toute fonction $u_{h} \in V_{h}$ , on a en effet :

u_{h} = \sum_{i = 1}^{P} u_{h} (N_{i}) ϕ_{i} (19)

Le support de $ϕ_{i}$ est constitué des triangles qui ont le nœud N_i en commun. La figure suivante montre le support de $ϕ_{13}$ et de $ϕ_{1}$ :

Figure pleine page

Pour un noeud N_i, la fonction $ϕ_{i}$ s'annule sur le côté opposé de chaque triangle auxquels ce nœud appartient.

Considérons un nœud N_i et ses nœuds voisins. Le support de $ϕ_{i}$ est constitué des triangles définis avec le nœud N_i et ces voisins. La figure suivante montre ces triangles :

Figure pleine page

Considérons un de ces triangles, dont les sommets sont notés 1,2,3, le sommet 1 étant le nœud N_i. Les coefficients $c_{0}, c_{1}, c_{2}$ de $ϕ_{i}$ sur ce triangle sont solutions des trois équations suivantes :

\begin{matrix} c_{0} x_{1} + c_{1} y_{1} + c_{2} = 1 \\ c_{0} x_{2} + c_{1} y_{2} + c_{2} = 0 \\ c_{0} x_{3} + c_{1} y_{3} + c_{2} = 0 \end{matrix}

Si detT désigne le déterminant de ce système, il est facile de démontrer que, si les sommets N₁,N₂,N₃ sont ordonnés dans le sens direct, l'aire du triangle s'écrit :

S_{T} = \frac{1}{2} (\vec{N_{1} N_{2}} \land \vec{N_{1} N_{3}}) \cdot \vec{u_{z}} = \frac{1}{2} d e t T (20)

Si les trois points ne sont pas alignés, ce système a un déterminant non nul. L'unique solulion peut être obtenue sans difficulté avec les règles de Cramer :

c_{0} = \frac{y_{2} - y_{3}}{d e t T}, c_{1} = \frac{x_{3} - x_{2}}{d e t T}, c_{2} = \frac{x_{2} y_{3} - x_{3} y_{2}}{d e t T} (21)

Soit V_h,0 l'ensemble des fonctions affines par morceaux sur le maillage et qui s'annulent sur les nœuds qui appartiennent à la frontière $\partial Ω_{d}$ (la partie de la frontière où une condition de Dirichlet est appliquée). Notons $N i n t$ les indices de tous les nœuds à l'exception de ces nœuds (indices des nœuds intérieurs). Pour toute fonction $v \in V_{h, 0}$ on a :

v = \sum_{i \in N i n t} v (N_{i}) ϕ_{i} (22)

La formulation variationnelle s'énonce comme suit. Trouver une fonction $u_{h} \in V_{h}$ telle que pour toute fonction $v \in V_{h, 0}$ on ait :

\int_{Ω} a (\vec{g r a d} u_{h}) \cdot (\vec{g r a d} v) d σ = \int_{Ω} v s d σ + \int_{\partial Ω_{n}} v a (K_{0} (u_{h} - u_{0}) + g_{0}) d ℓ (23)

Un énoncé équivalent est : trouver une fonction $u_{h} \in V_{h}$ telle que :

\int_{Ω} a (\vec{g r a d} u_{h}) \cdot (\vec{g r a d} ϕ_{i}) d σ = \int_{Ω} ϕ_{i} s d σ + \int_{\partial Ω_{n}} ϕ_{i} a (K_{0} (u_{h} - u_{0}) + g_{0}) d ℓ, i \in N i n t (24)

2.d. Équations algébriques

Introduisons les composantes de u_h sur la base :

u_{h} = \sum_{j = 1}^{P} ξ_{j} ϕ_{j} (25)

Les $ξ_{i}$ sont les valeurs de u_h aux nœuds, autrement dit les approximations de la solution u. Le report de cette décomposition dans $(23)$ conduit à un système d'équations linéaires pour les composantes $ξ_{i}$ :

\sum_{j = 1}^{P} ξ_{j} \int_{Ω} a (\vec{g r a d} ϕ_{j}) \cdot (\vec{g r a d} ϕ_{i}) d σ = \int_{Ω} ϕ_{i} s d σ + \sum_{j = 1}^{P} \int_{\partial Ω_{n}} ϕ_{i} a K_{0} ξ_{j} ϕ_{j} d ℓ + \int_{\partial Ω_{n}} ϕ_{i} (a g_{0} - K_{0} u_{0}) d ℓ, i \in N i n t (26)

Ce système comporte plus d'inconnues que d'équations : le nombre d'équations est égal aux nombres de nœuds à l'exception de ceux de la frontière où une condition de Dirichlet est appliquée alors que le nombre d'inconnues est égal aux nombre total de nœuds. Il faut donc enlever du système les $ξ_{i}$ des nœuds qui sont sur $\partial Ω_{d}$ . Pour cela, on écrit u_h comme la somme d'une fonction de V_h,0 (qui s'annule sur la frontière de Dirichlet) et d'une fonction de V_h,d, l'ensemble des fonctions affines par morceaux qui vérifient précisément la condition de Dirichlet choisie :

u_{h} = u_{h, 0} + u_{h, d} (27)

Il peut être démontré que u_h est indépendant du choix de u_h,d [1]. On choisit donc la fonction la plus simple : celle qui s'annule sur tous les nœuds à l'exception de ceux où la condition de Dirichlet est appliquée. Notons Nd les indices des nœuds où une condition de Dirichlet est appliquée. On a :

u_{h} = \sum_{j \in N i n t} ξ_{j} ϕ_{j} + \sum_{k \in N d} ξ_{k} ϕ_{k} (28)

Les $ξ_{k}$ dans la seconde somme sont connues puisqu'il s'agit des valeurs de u sur les nœuds appartenant à la frontière où une condition de Dirichlet est appliquée.

On obtient un système d'équations pour les inconnues $ξ_{j}, j \in N i n t$ :

\begin{matrix} \sum_{j \in N i n t} ξ_{j} \int_{Ω} a (\vec{g r a d} ϕ_{j}) \cdot (\vec{g r a d} ϕ_{i}) d σ + \sum_{k \in N d} ξ_{k} \int_{Ω} a (\vec{g r a d} ϕ_{k}) \cdot (\vec{g r a d} ϕ_{i}) d σ \\ = \int_{Ω} ϕ_{i} s d σ + \sum_{j \in N i n t} \int_{\partial Ω_{n}} ϕ_{i} a K_{0} ξ_{j} ϕ_{j} d ℓ + \sum_{k \in N d} \int_{\partial Ω_{n}} ϕ_{i} a K_{0} ξ_{k} ϕ_{k} d ℓ + \int_{\partial Ω_{n}} ϕ_{i} a (g_{0} - K_{0} u_{0}) d ℓ, i \in N i n t \end{matrix}

Réécrivons ces équations en mettant à gauche la combinaison linéaire des inconnues :

\begin{matrix} \sum_{j \in N i n t} ξ_{j} \int_{Ω} a (\vec{g r a d} ϕ_{j}) \cdot (\vec{g r a d} ϕ_{i}) d σ - \sum_{j \in N i n t} ξ_{j} \int_{\partial Ω_{n}} ϕ_{i} a K_{0} ϕ_{j} d ℓ \\ = - \sum_{k \in N d} ξ_{k} \int_{Ω} a (\vec{g r a d} ϕ_{k}) \cdot (\vec{g r a d} ϕ_{i}) d σ + \int_{Ω} ϕ_{i} s d σ + \sum_{k \in N d} \int_{\partial Ω_{n}} ϕ_{i} a K_{0} ξ_{k} ϕ_{k} d ℓ + \int_{\partial Ω_{n}} ϕ_{i} a (g_{0} - K_{0} u_{0}) d ℓ, i \in N i n t \end{matrix}

Posons :

\begin{matrix} A_{i j} = \int_{Ω} a (\vec{g r a d} ϕ_{j}) \cdot (\vec{g r a d} ϕ_{i}) d σ, (i, j) \in N i n t & (29) \\ R_{i j} = - \int_{\partial Ω_{n}} a K_{0} ϕ_{i} ϕ_{j} d ℓ, (i, j) \in N i n t & (30) \\ B_{i} = - \sum_{k \in N d} ξ_{k} \int_{Ω} a (\vec{g r a d} ϕ_{k}) \cdot (\vec{g r a d} ϕ_{i}) d σ + \int_{Ω} ϕ_{i} s d σ + \sum_{k \in N d} ξ_{k} \int_{\partial Ω_{n}} ϕ_{i} a K_{0} ϕ_{k} d ℓ + \int_{\partial Ω_{n}} ϕ_{i} a (g_{0} - K_{0} u_{0}) d ℓ, i \in N i n t & (31) \end{matrix}

Le système d'équations linéaires s'écrit :

\sum_{j \in N i n t} A_{i j} ξ_{j} + \sum_{j \in N i n t} R_{i j} ξ_{j} = B_{i}, i \in N i n t (32)

2.e. Matrice du système

La matrice du système d'équations est M=A+R. La matrice A est indépendante des conditions limites. La matrice R est non nulle seulement si une condition de Robin ( $K_{0} \neq 0$ ) est appliquée sur une partie de la frontière $\partial Ω_{n}$ .

Pour $i \neq j$ , le coefficient A_ij est non nul si les nœuds i et j ont une partie de leur support en commun, autrement dit si ce sont deux nœuds voisins (rappelons que les nœuds où une condition de Dirichlet est appliquée sont exclus). Pour évaluer approximativement l'intégrale sur un triangle T dont les sommets sont N₁,N₂,N₃, on utilise la formule de quadrature du barycentre :

\int_{T} f d σ \approx f (\frac{N_{1} + N_{2} + N_{3}}{3}) S_{T} (33)

où S_T est l'aire du triangle. Considérons $A_{i j}^{T}$ , l'intégrale sur un des deux triangles dont un côté est constitué des nœuds i et j, le troisième étant noté N_k. Soient $c_{i 0}, c_{i 1}, c_{i 2}$ les coefficients qui définissent $ϕ_{i}$ sur ce triangle et $c_{j 0}, c_{j 1}, c_{j 2}$ les coefficients qui définissent $ϕ_{j}$ sur le même triangle.

Figure pleine page

On a :

(\vec{g r a d} ϕ_{j}) \cdot (\vec{g r a d} ϕ_{i}) = c_{j 0} c_{i 0} + c_{j 1} c_{i 1} (34)

La formule du barycentre donne :

A_{i j}^{T} = a (\frac{N_{i} + N_{j} + N_{k}}{3}) (c_{j 0} c_{i 0} + c_{j 1} c_{i 1}) S_{T} (35)

En pratique, il sera plus simple d'attribuer une valeur de a à chaque nœud et on écrira donc :

A_{i j}^{T} = \frac{1}{3} (a_{i} + a_{j} + a_{k}) (c_{j 0} c_{i 0} + c_{j 1} c_{i 1}) S_{T} (36)

Pour chacun des triangles auxquels le nœud N_i appartient, on calcule A^T_ij pour les deux autres nœuds (d'indice $j \in N i n t$ ) et on ajoute le résultat au coefficient A_ij de la matrice (pour chaque sommet voisin, il y a deux termes).

Il faut aussi déterminer le coefficient sur la diagonale :

A_{i i} = \int_{Ω} a {(\vec{g r a d} ϕ_{i})}^{2} d σ (37)

Le support de la fonction à intégrer est constitué des triangles auquel le nœud N_i appartient. Pour un de ces triangles, l'intégrale vaut approximativement :

A_{i i}^{T} = \frac{1}{3} (a_{i} + a_{j} + a_{k}) (c_{i 0}^{2} + c_{i 1}^{2}) S_{T} (38)

Pour $i \neq j$ , R_ij est non nul si les deux nœuds N_i et N_j sont voisins et s'ils appartiennent à une frontière où une condition de Robin est appliquée.

Figure pleine page

Le produit $ϕ_{i} ϕ_{j}$ est non nul sur le triangle constitué de ces deux nœuds et du nœud N_k. R_ij est l'intégrale sur le segment N_iN_j de $f = - a K_{0} ϕ_{i} ϕ_{j}$ . Les coefficients a et K₀ seront dans la plupart des cas constants sur le segment ou peuvent être considérés comme tel. Si $ℓ$ est l'abscisse sur le segment, on a :

f (ℓ) = - a K_{0} \frac{ℓ (d - ℓ)}{d^{2}} (39)

où d est la distance entre les deux nœuds. On a donc :

R_{i j} = - \frac{a K_{0}}{d^{2}} \int_{0}^{d} ℓ (d - ℓ) d ℓ = - \frac{a K_{0} d}{6} (40)

Il faut aussi déterminer le coefficient sur la diagonale :

R_{i i} = - \int_{\partial Ω_{n}} a K_{0} ϕ_{i}^{2} d ℓ (41)

La fonction à intégrer est non nulle sur les deux segments (notés 1 et 2) de la frontière auquel le nœud N_i appartient :

Figure pleine page

Notons a₁,K₁ et a₂,K₂ les valeurs des coefficients aux milieux des segments 1 et 2. On a :

R_{i i} \approx - a_{1} K_{1} \int_{0}^{d_{1}} \frac{ℓ^{2}}{d_{1}^{2}} d ℓ - a_{2} K_{2} \int_{0}^{d_{2}} \frac{ℓ^{2}}{d_{2}^{2}} d ℓ = - \frac{1}{3} a_{1} K_{1} d_{1} - \frac{1}{3} a_{2} K_{2} d_{2} (42)

Les matrices A et R sont symétriques. La ligne i de la matrice A comporte un coefficient diagonal et des coefficients non diagonaux correspondant aux nœuds voisins du nœud N_i. La ligne i de la matrice R comporte un coefficients diagonal et des coefficients non diagonaux correspondant aux nœuds voisins de N_i situés sur le bord $\partial Ω_{n}$ .

Rappelons que tous les nœuds qui interviennent dans le calcul de ces deux matrices sont des nœuds intérieurs, c'est-à-dire qu'il faut exclure les nœuds où une condition de Dirichlet et appliquée.

2.f. Matrice colonne

Il faut expliciter les coefficients de la matrice colonne B. Le premier terme de B_ij est :

B_{1, i} = - \sum_{k \in N d} ξ_{k} \int_{Ω} a (\vec{g r a d} ϕ_{k}) \cdot (\vec{g r a d} ϕ_{i}) d σ, i \in N i n t (43)

Pour $i \in N i n t$ (nœud intérieur), il faut considérer les nœuds N_k (de la frontière où une condition de Dirichlet est appliquée) tels que le support de $ϕ_{k}$ ait une partie commune avec celui de $ϕ_{i}$ . $B_{1, i}$ est non nul si et seulement si N_i est voisin d'un nœud de la frontière de Dirichlet. La figure suivante montre un nœud de ce type et les triangles qui constituent la partie commune des supports.

Figure pleine page

L'intégrale sur un triangle se calcule comme indiqué plus haut au sujet de A_ij. Pour le nœud N_k voisin de N_i, dont la valeur est $ξ_{k}$ , l'intégrale se fait sur deux triangles, les deux triangles qui partagent le côté N_iN_k. Pour chacun des triangles auxquels le nœuds N_i appartient, on repère parmi les deux autres nœuds ceux qui sont sur la frontière de Dirichlet puis on ajoute le terme correspondant dans B_i. Pour chaque triangle auquel le nœud N_i appartient, on identifie ceux pour lesquels au moins un des deux autres nœuds est sur un bord de Dirichlet (triangles grisés sur la figure ci-dessus). À chacun de ces triangles correspond un ou deux termes dans $(43)$ .

Le deuxième terme de B_ij est :

B_{2, i} = \int_{Ω} ϕ_{i} s d σ, i \in N i n t, i \in N i n t (44)

Il s'agit du terme qui contient la source de l'équation de Poisson. Pour $i \in N i n t$ , il faut évaluer l'intégrale sur le support de $ϕ_{i}$ , qui est constitué des triangles auxquels ce nœud appartient. Pour un de ces triangles, dont les sommets sont N₁,N₂,N₃, la formule de quadrature du barycentre donne :

B_{2, i} \approx s (\frac{N_{1} + N_{2} + N_{3}}{3}) ϕ_{i} (\frac{N_{1} + N_{2} + N_{3}}{3}) S_{T} (45)

La fonction $ϕ_{i}$ vaut 1/3 au barycentre. Pour la source, on considère plutôt la moyenne des valeurs aux sommets du triangle :

B_{2, i} \approx \frac{s_{1} + s_{2} + s_{3}}{3} \frac{1}{3} S_{T} (46)

Le troisième terme de B_ij est :

B_{3, i} = \sum_{k \in N d} ξ_{k} \int_{\partial Ω_{n}} a K_{0} ϕ_{i} ϕ_{k} d ℓ, i \in N i n t (47)

Ce coefficient peut être non nul s'il existe un nœud N_k de la frontière de Dirichlet tel que $ϕ_{k}$ soit non nul sur au moins un segment de la frontière de Neumann. Cela ne peut se rencontrer qu'à la jonction entre ces deux frontières, comme le montre la figure suivante :

Figure pleine page

L'intégrale se fait sur le segment N_kN_i. Le calcul est similaire à celui du coefficient R_ij. On obtient :

B_{3, i} = ξ_{k} a K_{0} \frac{d}{6} (48)

où d est la distance N_kN_i.

Le quatrième terme de B_ij est :

B_{4, i} = \int_{\partial Ω_{n}} ϕ_{i} a (g_{0} - K_{0} u_{0}) d ℓ, i \in N i n t (49)

Ce terme contient la condition de Neumann (gradient g₀) et le terme constant de la condition de Robin. Il est non nul si et seulement si le nœud N_i appartient à la frontière $\partial Ω_{n}$ . L'intégrale se fait sur le ou les deux segments de la frontière adjacents à ce nœud.

3. Implémentation

La structure de données qui contient la définition du maillage est une liste de nœuds qui est générée au moment de la création du maillage. Chaque nœud N_i est un objet qui contient les informations suivantes :

Ses coordonnées.
La liste des nœuds voisins, classés dans le sens de rotation direct. Deux nœuds voisins dans cette liste constituent avec le nœud N_i un triangle. On peut donc à partir cette liste obtenir rapidement tous les triangles auquel le nœud N_i appartient.
Le type du nœud : sur une frontière de Neumann, sur une frontière de Dirichlet ou autre.
Les coefficients relatifs à la condition limite, en particulier la valeur de $ξ_{i}$ si le nœud est sur une condition limite de Dirichlet.
Les valeurs de a et s. Si s est définie par une fonction, on fait le choix d'évaluer cette fonction sur les nœuds et non pas aux barycentres des triangles. La valeur au barycentre sera simplement évaluée comme la moyenne des valeurs aux sommets.

import numpy as np
import scipy.sparse as sparse
import scipy.sparse.linalg as linalg
from scipy.linalg import solve

DIR = 1
NEU = 2

class Noeud:
    def __init__(self,x,y,a,s,lim=0,ksi=0,K0=0,u0=0,g0=0,bord=False):
        self.xy = np.array([x,y])
        self.x = x 
        self.y = y
        self.a = a
        self.s = s
        self.lim = lim
        self.ksi = ksi
        self.K0 = K0
        self.u0 = u0
        self.g0 = g0
        self.noeuds_voisins = []
        self.indice = 0 # indice dans le maillage
        self.bord= bord
    def ajouterVoisin(self,noeud):
        self.noeuds_voisins.append(noeud)              
    def extraireTriangles(self):
        N = len(self.noeuds_voisins)
        triangles =  []
        for i in range(N-1):
            n1 = self
            n2 = self.noeuds_voisins[i]
            n3 = self.noeuds_voisins[i+1]
            triangles.append([n1,n2,n3])
        if not(self.bord):
            n1 = self
            n2 = self.noeuds_voisins[N-1]
            n3 = self.noeuds_voisins[0]
            triangles.append([n1,n2,n3])   
            
        return triangles

La fonction extraireTriangles extrait la liste des triangles auxquels le nœud appartient. Chaque triangle est une liste de trois nœuds, le premier étant le nœud central. Dans un triangle, les nœuds sont rangés dans le sens direct, c'est-à-dire que $(\vec{N_{1} N_{2}} \land \vec{N_{1} N_{3}}) \cdot \vec{u_{z}} > 0$ . Lorsqu'un nœud n'est pas sur un bord, le dernier nœud de la liste des voisins forme un triangle avec le premier.

class ElementsFinis2D:
    def __init__(self):
        self.noeuds_dirichlet = []
        self.noeuds = []
        self.indice = 0
        self.indice_dirichlet = 0
    def ajouterNoeud(self,noeud):
        if noeud.lim == DIR:
            noeud.indice = self.indice_dirichlet
            self.noeuds_dirichlet.append(noeud)
            self.indice_dirichlet += 1
        else:
            noeud.indice = self.indice
            self.noeuds.append(noeud)
            self.indice += 1
    def definirMaillage(self,maillage):
        self.noeuds_dirichlet = []
        self.noeuds = []
        self.indice = 0
        self.indice_dirichlet = 0
        for noeud in maillage:
            self.ajouterNoeud(noeud)
    def fonctionPhi(self,triangle):
        n1 = triangle[0]
        n2 = triangle[1]
        n3 = triangle[2]
        detT = n2.x*n3.y-n2.x*n1.y-n1.x*n3.y-n2.y*n3.x+n2.y*n1.x+n1.y*n3.x
        if detT <=0 : 
            print("Erreur de définition du triangle")
            print(n1.xy)
            print(n2.xy)
            print(n3.xy)
        ST = 0.5*detT
        c0 = (n2.y-n3.y)/detT
        c1 = (n3.x-n2.x)/detT 
        c2 = (n2.x*n3.y-n3.x*n2.y)/detT
        return [ST,c0,c1,c2]
    def creerMatrices(self):
        self.Nnoeuds = len(self.noeuds)
        self.M = np.zeros((self.Nnoeuds,self.Nnoeuds),dtype=np.float64)
        self.B = np.zeros(self.Nnoeuds,dtype=np.float64)
    def matriceA(self):
        for i in range(self.Nnoeuds):
            Ni = self.noeuds[i]
            triangles = Ni.extraireTriangles()
            Aii = 0
            for t in range(len(triangles)):
                [ST,ci0,ci1,ci2] = self.fonctionPhi(triangles[t])
                Aii += (triangles[t][0].a+triangles[t][1].a+triangles[t][2].a)/3*(ci0*ci0+ci1*ci1)*ST
                if triangles[t][1].lim != DIR:
                    j = triangles[t][1].indice
                    [ST,cj0,cj1,cj2] = self.fonctionPhi([triangles[t][1],triangles[t][2],triangles[t][0]])
                    self.M[i,j] += (triangles[t][0].a+triangles[t][1].a+triangles[t][2].a)/3*(ci0*cj0+ci1*cj1)*ST
                if triangles[t][2].lim != DIR:
                    j = triangles[t][2].indice
                    [ST,cj0,cj1,cj2] = self.fonctionPhi([triangles[t][2],triangles[t][0],triangles[t][1]])
                    self.M[i,j] += (triangles[t][0].a+triangles[t][1].a+triangles[t][2].a)/3*(ci0*cj0+ci1*cj1)*ST
            self.M[i,i] += Aii
    
    def matriceR(self):
        for i in range(self.Nnoeuds):
            Ni = self.noeuds[i]
            Rii = 0
            if Ni.lim==NEU:
                for n in Ni.noeuds_voisins:
                    if n.lim == NEU:
                        d = np.sqrt((Ni.x-n.x)**2+(Ni.y-n.y)**2)
                        Rii += -(Ni.a*Ni.K0+n.a*n.K0)/2*d/3
                        j = n.indice
                        self.M[i,j] += -1/6*(Ni.a*Ni.K0+n.a*n.K0)/2*d
            self.M[i,i] += Rii
    def matriceB1(self):
        for i in range(self.Nnoeuds):
            B1i = 0
            Ni = self.noeuds[i]
            triangles = Ni.extraireTriangles()
            for t in range(len(triangles)):
                [ST,ci0,ci1,ci2] = self.fonctionPhi(triangles[t])
                if triangles[t][1].lim == DIR :
                    [ST,cj0,cj1,cj2] = self.fonctionPhi([triangles[t][1],triangles[t][2],triangles[t][0]])
                    B1i += -triangles[t][1].ksi*(triangles[t][0].a+triangles[t][1].a+triangles[t][2].a)/3*(ci0*cj0+ci1*cj1)*ST
                if triangles[t][2].lim == DIR:
                    [ST,cj0,cj1,cj2] = self.fonctionPhi([triangles[t][2],triangles[t][0],triangles[t][1]])
                    B1i += -triangles[t][2].ksi*(triangles[t][0].a+triangles[t][1].a+triangles[t][2].a)/3*(ci0*cj0+ci1*cj1)*ST
            self.B[i] += B1i
    def matriceB2(self):
        for i in range(self.Nnoeuds):
            B2i = 0
            Ni = self.noeuds[i]
            triangles = Ni.extraireTriangles()
            for t in range(len(triangles)):
                [ST,ci0,ci1,ci2] = self.fonctionPhi(triangles[t])
                B2i += (triangles[t][0].s+triangles[t][1].s+triangles[t][2].s)/9*ST
            self.B[i] += B2i
    def matriceB3(self):
        for i in range(self.Nnoeuds):
            Ni = self.noeuds[i]
            if Ni.lim == NEU:
                B3i = 0
                for n in Ni.noeuds_voisins:
                    if n.lim == DIR:
                        d = np.sqrt((Ni.x-n.x)**2+(Ni.y-n.y)**2)
                        B3i += n.ksi*(Ni.a*Ni.K0+n.a*n.K0)/2*d/6
                self.B[i] += B3i
    def matriceB4(self):
        for i in range(self.Nnoeuds):
            Ni = self.noeuds[i]
            if Ni.lim == NEU:
                B4i = 0
                for n in Ni.noeuds_voisins:
                    if n.lim == NEU:
                        d = np.sqrt((Ni.x-n.x)**2+(Ni.y-n.y)**2)
                        B4i += 0.5*(Ni.a*(Ni.g0-Ni.K0*Ni.u0)+n.a*(n.g0-n.K0*n.u0))*0.5*d
                self.B[i] += B4i
    def solve(self,solver="iterative"):
        self.creerMatrices()
        self.matriceA()
        self.matriceR()
        self.matriceB1()
        self.matriceB2()
        self.matriceB3()
        self.matriceB4()
        if solver=="iterative":
            M = sparse.bsr_array(self.M)
            ksi,exitCode = linalg.cg(M, self.B, atol=1e-5)
            print("exitCode : ",exitCode)
        else:
            ksi = solve(self.M,self.B)
        for i in range(self.Nnoeuds):
            self.noeuds[i].ksi = ksi[i]

La fonction ajouterNoeud ajoute un nœud, soit dans la liste des nœuds internes, soit dans la liste des nœuds d'une frontière de Dirichlet. L'indice dans la liste est affecté au nœud. La fonction definirMaillage définit le maillage complet donné par une liste de nœuds.

La fonction fonctionPhi calcule, pour un triangle donné (liste de trois nœuds), l'aire S_T et les coefficients c₀,c₁,c2 de la fonction $ϕ_{i}$ relative au premier nœud du triangle.

La fonction creerMatrices crée les deux matrices qui définissent le système linéaire. Les fonctions matriceX remplissent ces matrices en calculant les différents coefficients A_ij,R_ij,B_1i,B_2i,B_3i,B_4i.

La fonction solve appelle les fonctions de création et de remplissage précédentes et calcule la solution du système, par défaut par la méthode itérative du gradient conjugué. Les valeurs de $ξ_{i}$ obtenues sont affectées aux nœuds.

4. Exemples

4.a. Exemple A

Définition

On considère le domaine carré $Ω : [0, 1] \times [0, 1]$ . La source est nulle (s=0) et le coefficient a est uniforme (le même partout). En x=0 et x=1 on applique des conditions limites de Dirichlet, avec respectivement les valeurs 0 et 1. En y=0 et y=1 on applique des conditions de Neumann avec g₀=0. La solution est : u(x,y)=x (il s'agit en fait d'un problème unidimensionnel).

La structure du maillage est montré sur la figure suivante :

Figure pleine page

Les 4 nœuds des coins sont affectés de la condition de Dirichlet. Chaque nœud non situé sur un bord a 6 nœuds voisins.

Génération du maillage

Il s'agit de générer la liste des nœuds puis d'attribuer les voisins de chaque nœud. Soit N le nombre de nœuds sur chaque axe.

N = 50
x = np.linspace(0,1,N)
y = np.linspace(0,1,N)
maillage = []
a = 1
s = 0
for j in range(N):
    for i in range(N):
        if i!=0 and i!=N-1 and j!=0 and j!=N-1:
            n = Noeud(x[i],y[j],a,s)
        else:  
            if i==0:
                n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
            elif i==N-1:
                n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=1,bord=True)
            if j==0 and i!=0 and i!=N-1:
                n = Noeud(x[i],y[j],a,s,lim=NEU,g0=0,bord=True)
            elif j==N-1 and i!=0 and i!=N-1:
                n = Noeud(x[i],y[j],a,s,lim=NEU,g0=0,bord=True)
        
        maillage.append(n)

L'attribution des voisins pour chaque nœud dépend seulement du maillage. On la réalise donc dans une fonction afin de pouvoir la réutiliser pour d'autres exemples :

def attribuerVoisins(N,maillage):
    indice = 0
    for j in range(N):
        for i in range(N):
            n = maillage[indice]
            if i!=0 and i!=N-1 and j!=0 and j!=N-1:
                n.ajouterVoisin(maillage[indice+1])
                n.ajouterVoisin(maillage[indice+N+1])
                n.ajouterVoisin(maillage[indice+N])
                n.ajouterVoisin(maillage[indice-1])
                n.ajouterVoisin(maillage[indice-1-N])
                n.ajouterVoisin(maillage[indice-N])
            else:
                if i==0 and j==0:
                    n.ajouterVoisin(maillage[indice+1])
                    n.ajouterVoisin(maillage[indice+N+1])
                    n.ajouterVoisin(maillage[indice+N])
                elif i==0 and j==N-1:
                    n.ajouterVoisin(maillage[indice-N])
                    n.ajouterVoisin(maillage[indice+1])
                elif i==0:
                    n.ajouterVoisin(maillage[indice-N])
                    n.ajouterVoisin(maillage[indice+1])
                    n.ajouterVoisin(maillage[indice+N+1])
                    n.ajouterVoisin(maillage[indice+N])
                if i==N-1 and j==0:
                    n.ajouterVoisin(maillage[indice+N])
                    n.ajouterVoisin(maillage[indice-1])
                elif i==N-1 and j==N-1:
                    n.ajouterVoisin(maillage[indice-1])
                    n.ajouterVoisin(maillage[indice-1-N])
                    n.ajouterVoisin(maillage[indice-N])
                elif i==N-1:
                    n.ajouterVoisin(maillage[indice+N])
                    n.ajouterVoisin(maillage[indice-1])
                    n.ajouterVoisin(maillage[indice-1-N])
                    n.ajouterVoisin(maillage[indice-N])
                if j==0 and i!=0 and i!=N-1:
                    n.ajouterVoisin(maillage[indice+1])
                    n.ajouterVoisin(maillage[indice+N+1])
                    n.ajouterVoisin(maillage[indice+N])
                    n.ajouterVoisin(maillage[indice-1])
                if j==N-1 and i!=0 and i!=N-1:
                    n.ajouterVoisin(maillage[indice-1])
                    n.ajouterVoisin(maillage[indice-N-1])
                    n.ajouterVoisin(maillage[indice-N])
                    n.ajouterVoisin(maillage[indice+1])
            indice+= 1
            
attribuerVoisins(N,maillage)

Résolution et tracé de u(x,y)

Pour le maillage considéré, il est aisé de tracer u(x,y) car les nœuds sont disposés sur une grille. Dans le cas général, il faudra procéder à une interpolation dans les triangles (par la fonction affine par morceaux) pour obtenir les valeurs de u sur une grille.

elem = ElementsFinis2D()
elem.definirMaillage(maillage)
elem.solve()
U = np.zeros((N,N),dtype=np.float64)
indice = 0
for j in range(N):
    for i in range(N):
        n = maillage[indice]
        U[j,i] = n.ksi
        indice+= 1
from matplotlib.pyplot import *
figure()
contour(U,10,extent=(0,1,0,1))

fig1.pdf

figure()
plot(x,U[N//2,:])
grid()

fig2.pdf

4.b. Exemple B

On reprend l'exemple A avec une source uniforme et des conditions de Dirichlet nulle en x=0 et x=1 :

N = 50
x = np.linspace(0,1,N)
y = np.linspace(0,1,N)
maillage = []
a = 1
s = 1
for j in range(N):
    for i in range(N):
        if i!=0 and i!=N-1 and j!=0 and j!=N-1:
            n = Noeud(x[i],y[j],a,s)
        else:  
            if i==0:
                n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
            elif i==N-1:
                n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
            if j==0 and i!=0 and i!=N-1:
                n = Noeud(x[i],y[j],a,s,lim=NEU,g0=0,bord=True)
            elif j==N-1 and i!=0 and i!=N-1:
                n = Noeud(x[i],y[j],a,s,lim=NEU,g0=0,bord=True)
        
        maillage.append(n)
attribuerVoisins(N,maillage)

elem = ElementsFinis2D()
elem.definirMaillage(maillage)
elem.solve()
U = np.zeros((N,N),dtype=np.float64)
indice = 0
for j in range(N):
    for i in range(N):
        n = maillage[indice]
        U[j,i] = n.ksi
        indice+= 1
from matplotlib.pyplot import *
figure()
contour(U,10,extent=(0,1,0,1))

fig3.pdf

figure()
plot(x,U[N//2,:])
grid()

fig4.pdf

4.c. Exemple C

On reprend l'exemple A avec une source uniforme et des conditions de Dirichlet nulle sur les quatre bords :

N = 50
x = np.linspace(0,1,N)
y = np.linspace(0,1,N)
maillage = []
a = 1
s = 1
for j in range(N):
    for i in range(N):
        if i!=0 and i!=N-1 and j!=0 and j!=N-1:
            n = Noeud(x[i],y[j],a,s)
        else:  
            if i==0:
                n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
            elif i==N-1:
                n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
            if j==0 and i!=0 and i!=N-1:
                n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
            elif j==N-1 and i!=0 and i!=N-1:
                n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
        
        maillage.append(n)
attribuerVoisins(N,maillage)

elem = ElementsFinis2D()
elem.definirMaillage(maillage)
elem.solve()
U = np.zeros((N,N),dtype=np.float64)
indice = 0
for j in range(N):
    for i in range(N):
        n = maillage[indice]
        U[j,i] = n.ksi
        indice+= 1
from matplotlib.pyplot import *
figure()
contour(U,10,extent=(0,1,0,1))

fig5.pdf

figure()
plot(x,U[N//2,:])
grid()

fig6.pdf

4.d. Exemple D

On reprend l'exemple A avec deux milieux différents : la moitié à droite du domaine a une valeur de a deux fois plus grande.

N = 50
x = np.linspace(0,1,N)
y = np.linspace(0,1,N)
maillage = []
a1 = 1
a2 = 2
s = 0
for j in range(N):
    for i in range(N):
        if i < N//2 : a=a1
        else : a=a2
        if i!=0 and i!=N-1 and j!=0 and j!=N-1:
            n = Noeud(x[i],y[j],a,s)
        else:  
            if i==0:
                n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
            elif i==N-1:
                n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=1,bord=True)
            if j==0 and i!=0 and i!=N-1:
                n = Noeud(x[i],y[j],a,s,lim=NEU,g0=0,bord=True)
            elif j==N-1 and i!=0 and i!=N-1:
                n = Noeud(x[i],y[j],a,s,lim=NEU,g0=0,bord=True)
        
        maillage.append(n)
attribuerVoisins(N,maillage)

elem = ElementsFinis2D()
elem.definirMaillage(maillage)
elem.solve()
U = np.zeros((N,N),dtype=np.float64)
indice = 0
for j in range(N):
    for i in range(N):
        n = maillage[indice]
        U[j,i] = n.ksi
        indice+= 1
from matplotlib.pyplot import *
figure()
contour(U,10,extent=(0,1,0,1))

fig7.pdf

figure()
plot(x,U[N//2,:])
grid()

fig8.pdf

4.e. Exemple E

Un objet de forme carrée est introduit dans le domaine, avec une condition de Dirichlet sur ses bords. Les quatres bords du domaine ont une condition de Dirichlet. La figure suivante montre comment l'objet est défini :

Figure pleine page

N = 50
x = np.linspace(0,1,N)
y = np.linspace(0,1,N)
maillage = []
a = 1
s = 0
# bords de l'objet
e = 5
i1 = N//2-e 
i2 = N//2+e 
j1 = N//2-e 
j2 = N//2+e
for j in range(N):
    for i in range(N):
        if (i==i1 and j1<=j<=j2) or (i==i2 and j1<=j<=j2) or (j==j1 and i1<=i<=i2) or (j==j2 and i1<=i<=i2):
            n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=1)
        else:
            if i!=0 and i!=N-1 and j!=0 and j!=N-1:
                n = Noeud(x[i],y[j],a,s)
            else:  
                if i==0:
                    n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
                elif i==N-1:
                    n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
                if j==0 and i!=0 and i!=N-1:
                    n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
                elif j==N-1 and i!=0 and i!=N-1:
                    n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
        maillage.append(n)
attribuerVoisins(N,maillage)

elem = ElementsFinis2D()
elem.definirMaillage(maillage)
elem.solve()
U = np.zeros((N,N),dtype=np.float64)
indice = 0
for j in range(N):
    for i in range(N):
        n = maillage[indice]
        U[j,i] = n.ksi
        indice+= 1
from matplotlib.pyplot import *
figure()
contour(U,10,extent=(0,1,0,1))

fig9.pdf

figure()
plot(x,U[N//2,:])
grid()

fig10.pdf

4.f. Exemple F

On reprend l'exemple A mais avec une condition de Robin sur les bords en y=0 et y=1 :

N = 50
x = np.linspace(0,1,N)
y = np.linspace(0,1,N)
maillage = []
a = 1
s = 0
for j in range(N):
    for i in range(N):
        if i!=0 and i!=N-1 and j!=0 and j!=N-1:
            n = Noeud(x[i],y[j],a,s)
        else:  
            if i==0:
                n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=0,bord=True)
            elif i==N-1:
                n = Noeud(x[i],y[j],a,s,lim=DIR,ksi=1,bord=True)
            if j==0 and i!=0 and i!=N-1:
                n = Noeud(x[i],y[j],a,s,lim=NEU,g0=0,K0=-1,u0=0,bord=True)
            elif j==N-1 and i!=0 and i!=N-1:
                n = Noeud(x[i],y[j],a,s,lim=NEU,g0=0,K0=-1,u0=0,bord=True)
        
        maillage.append(n)
attribuerVoisins(N,maillage)

elem = ElementsFinis2D()
elem.definirMaillage(maillage)
elem.solve()
U = np.zeros((N,N),dtype=np.float64)
indice = 0
for j in range(N):
    for i in range(N):
        n = maillage[indice]
        U[j,i] = n.ksi
        indice+= 1
from matplotlib.pyplot import *
figure()
contour(U,10,extent=(0,1,0,1))

fig11.pdf

figure()
plot(x,U[N//2,:])
grid()

fig12.pdf

Équation de Poisson 2D : méthode des éléments finis

1. Introduction

2. Méthode des éléments finis

2.a. Formulation variationnelle

2.b. Conditions limites

2.c. Maillage triangulaire - Fonctions affines par morceaux

2.d. Équations algébriques

2.e. Matrice du système

2.f. Matrice colonne

3. Implémentation

4. Exemples

4.a. Exemple A

Définition

Génération du maillage

Résolution et tracé de u(x,y)

4.b. Exemple B

4.c. Exemple C

4.d. Exemple D

4.e. Exemple E

4.f. Exemple F

Références