Maillage et Éléments Finis

Bertrand Thierry

avr. 24, 2023

 Table des matières :

I Notes de cours 3
1 Simulation Numérique 7
1.1 Quelques exemples d’EDP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Industrie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2 Formulations Faibles 11
2.1 Espaces de Hilbert : rappels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Formulation Faible . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 Théorème de Lax-Milgram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4 Espaces de Sobolev . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.5 Application au problème modèle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3 Éléments Finis Triangulaires 25

3.1 Méthode de Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Espace P1 -Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3 Assemblage des Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.4 Calcul des Matrices Élémentaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.5 Matrice Creuse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4 Conditions aux bords 47

4.1 Conditions de Neumann hétérogène . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2 Condition de Dirichlet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.3 Condition de Fourier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

5 Avancée 57
5.1 Erreur comise et convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2 Éléments Finis P2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

II Implémentation 61
6 Maillage avec GMSH 63
6.1 Prise en main de GMSH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
6.2 API GMSH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

7 Solveur FEM Python 65

7.1 Matrices Creuses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

i
 7.2 Gestion du maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
7.3 Matrices de Masse et de Rigidité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
7.4 Quadratures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
7.5 Conditions de Dirichlet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
7.6 Résolution et Analyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

III Projet 75
8 2017 - 2018 : Wi-Fi 77
8.1 Équation de Helmholtz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
8.2 Modèle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
8.3 Implémentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
8.4 Étude et examen oral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
8.5 Résultat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

9 2020 - 2021 87
9.1 Problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
9.2 Travail demandé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
9.3 Consignes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

IV Download 91
Proof Index 93

ii
 Maillage et Éléments Finis

Le code source du cours est disponible librement sur github, n’hésitez pas à envoyer des remarques, typos, corrections,
…!
La première section présente la théorie et la mise en oeuvre des éléments finis. La deuxième partie propose un exemple
d’implémentation de la méthode des éléments en utilisant le language Python et le logiciel GMSH pour la partie maillage.

Table des matières : 1

Maillage et Éléments Finis

2 Table des matières :

Première partie

Notes de cours

3
 Maillage et Éléments Finis

Remarques préliminaires :
— Si la théorie mathématique ne vous branche pas, vous pouvez passer directement à la section pratique et revenir
plus tard sur la théorie, pour comprendre pourquoi tout cela fonctionne.
— Dans ce cours, un scalaire sera écrit normalement (𝑥) tandis qu’un vecteur sera noté en gras (x = (𝑥, 𝑦)). Nous
commençons la numérotation par 1 et non par 0.

5
Maillage et Éléments Finis

6
 CHAPITRE 1

Simulation Numérique

La dictature s’épanouit sur le terreau de l’ignorance.

---G. Orwell - 1984

1.1 Quelques exemples d’EDP

1.1.1 Thermique

Prenons un domaine borné et connexe Ω, par exemple le carré unité, qui représente un studio. En supposant le milieu
homogène, la température 𝑇 au sein du studio vérifie l’équation de Laplace :

(−Δ𝑇 ) = 0, dans Ω,
∑︀𝑑 2
où l’opérateur Δ := 𝜕
𝑖=1 𝜕𝑥2𝑖 est le Laplacien ou Opérateur de Laplace.
Supposons maintenant que le studio comporte une source de chaleur, par exemple un radiateur. Nous le modéli-
sons par une fonction 𝑞, continue, valant 0 partout sauf sur un petit domaine Ω𝑅 . En notant la conductivité 𝐾`(=
𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑒)𝑑𝑢𝑚𝑖𝑙𝑖𝑒𝑢𝑎𝑙𝑜𝑟𝑠 : 𝑚𝑎𝑡ℎ : `𝑇 vérifie l’équation de Poisson :
𝑞
(−Δ𝑇 ) = , dans Ω.
𝐾
Ce problème reste pour l’instant incomplet car il manque des conditions. Contrairement aux problèmes de Cauchy, il n’y
a pas de condition initiale, car le régime est stationnaire, mais nous avons besoin de conditions sur le bord 𝜕Ω du domaine
Ω. Nous parlons alors de conditions aux limites et par suite de problème aux limites.
Ajoutons une fenêtre à notre studio portée à la douce température de 𝑇0 = 10∘ 𝐶 grâce à l’automne frisquet. Le reste
des murs est supposé parfaitement isolants, autrement dit le flux à travers les parois est nul. Le flux d’une quantité sur une
interface étant donné par

flux := ∇𝑇 · n = 𝜕n 𝑇,

où n est le vecteur unitaire normale sortant au domaine. Le bord de notre appartement, noté Γ := 𝜕Ω, est alors divisé en
deux parties :

7
Maillage et Éléments Finis

— Γ𝐷 : la fenêtre sur laquelle la température est imposée : condition de Dirichlet.

— Γ𝑁 : les murs supposés isolants sur lesquels le flux est imposé : condition de Neumann.
Notre problème s’écrit alors :
⎧ 𝑞
⎨ −Δ𝑇 =
⎪ , dans Ω,
𝐾
𝑇 = 𝑇0 , sur Γ𝐷 ,
sur Γ𝑁 .
⎪
𝜕n 𝑇 = 0,
⎩

Nous verrons dans ce cours que ce problème admet une solution (ouf) et qui est, de plus, unique (re-ouf). Résoudre le
problème analytiquement (i.e. « à la main ») peut s’avérérer délicat, notamment si la géométrie est complexe : c’est ici que
la simulation numérique rentre en jeu et notamment les éléments finis.

1.1.2 Diffusion d’une onde Wifi dans un appartement (Projet 2017-2018)

Cet exemple est tiré du projet proposé en 2017 - 2018.

Modèle

La définition Wikipédia d’une onde est la suivante :

Une onde est la propagation d’une perturbation produisant sur son passage une variation réversible des pro-
priétés physiques locales du milieu. Elle se déplace avec une vitesse déterminée qui dépend des caractéris-
tiques du milieu de propagation. Une onde transporte de l’énergie sans transporter de matière.
Mathématiquement, une onde E (x, 𝑡) dépend du temps 𝑡 et de l’espace x, et vérifie l’équation des ondes :
1 𝜕2E
ΔE (x, 𝑡) = (x, 𝑡),
𝑐2 𝜕𝑡2
où 𝑐 est la célérité de l’onde dans le milieu (qui peut dépendre de la position x !). Par exemple, dans le cas d’une onde
électromagnétique et dans le vide, 𝑐 est la célérité de la lumière, soit 299792458m.s−1 . La quantité x est un vecteur de
dimension 𝑑 = 2 ou 𝑑 = 3 selon le problème considéré : dans notre cas 𝑑 = 2.
Lors d’une excitation périodique,√c’est-à-dire lorsque la pulsation 𝜔 (en rad.s−1 ) de l’onde est fixée, l’onde s’écrit alors
E (x, 𝑡) = ℜ (𝑢(x)𝑒 −𝚤𝜔𝑡
) où 𝚤 = −1 et 𝐸 est une onde spatiale satisfaisant l’équation de Helmholtz :
𝜔2
Δ𝐸 + 𝐸 = 𝑓.
𝑐2
Cette équation s’obtient en remplaçant E (x, 𝑡) par 𝐸(x)𝑒−𝚤𝜔𝑡 dans l’équation des ondes. Nous notons en général 𝑘 = 𝜔𝑐
(en rad.m−1 ) le nombre d’onde et 𝜆 = 2𝜋 𝑘 (en m) la longueur d’onde, autrement dit, la distance entre deux amplitudes,
de sorte que l’équation de Helmholtz s’écrit
Δ𝐸 + 𝑘 2 𝐸 = 𝑓.
La source 𝑓 est ici spatiale, dans le cas d’une source ponctuelle de centre s la source est alors un Dirac :
Δ𝐸 + 𝑘 2 𝐸 = −𝛿s .
Les ondes Wi-Fi qui suivent la norme IEEE 802.11g sont émises à une fréquence variant de 2.4GHz à 2.5GHz. L’ap-
partement tout entier dans lequel est situé notre routeur est noté Ω. Les murs sont supposés être du même matériau : du
placo-plâtre. Le domaine Ω = Ω𝑎 ∪ Ωmur est décomposé en deux domaines, Ω𝑎 pour l’air et Ωmur pour les murs.
Une modélisation possible de ce problème est le système d’équations suivant :
En supposant que l’air a les mêmes propriétés électromagnétiques que le vide
Δ𝐸(x) + 𝑘 2 𝑛(x)2 𝐸(x) = −𝛿routeur (x) dans Ω,
{︂
𝜕n 𝐸(x) − 𝚤𝑘𝑛(x)𝐸(x) = 0 sur 𝜕Ω,
où nous avons :

8 Chapitre 1. Simulation Numérique

 Maillage et Éléments Finis

— 𝛿routeur : position du routeur. Nous l’avons placé dans le salon.

— 𝑛 : fonction de contraste qui prend en compte les caractéristiques électromagnétiques du mur et de l’air :

si x ∈ Ω𝑎 (i.e. x est dans l'air),

{︂
1
𝑛(x) =
2.4 si x ∈ Ωmur (i.e. x est dans le mur).

Notez que ces valeurs sont des valeurs physiques et ne sont pas une lubie mathématique.
— La dernière équation, 𝜕n 𝐸 − 𝚤𝑘𝑛𝐸 = 0 est une condition de Fourier-Robin (ou de Fourier ou de Robin ou
même d’impédance). Elle a pour but d’absorber (avec un succès mitigé) les ondes sortantes, mimant un mur
« transparent » (sans réflexion d’ondes).

Résolution numérique

La résolution d’un tel problème dans un appartement deux pièces avec cuisine séparée (grand luxe Parisien) avec la
méthode des éléments finis donne alors ce résultat :

Fıg. 1 – Propagation d’une onde Wi-Fi dans un appartement. Après avoir traversé 2 murs, l’onde Wi-Fi semble très
amortie. Sous le résultat est affiché le plan de l’appartement et la position du routeur (petit disque à gauche)

Vous voulez tester ?

— Téléchargez le bundle Onelab. Il contient GMSH et GetDP (un solveur éléments finis)
— Téléchargez le code, soit directement soit via Git :

git clone [Link] wifi

— Dans le dossier et dans un terminal, lancer

gmsh [Link]

— Vous pouvez modifier un peu la géométrie et la fréquence de l’onde, mise à 1GHZ. Attention, cette simulation est
très gourmande : testez d’abord avec 1GHz avant de lancer la simulation pour 2.5GHz (au risque de faire crasher
votre ordinateur) !

1.1. Quelques exemples d’EDP 9

Maillage et Éléments Finis

1.1.3 Objectifs du cours

Vous serez capable de résoudre ce genre de problème (et bien d’autres) et, ainsi, d’épater votre famille lors de ces inter-
minables dîners.

1.2 Industrie

À venir…

1.2.1 Combustion

Cerfacs, logiciel AVBP

10 Chapitre 1. Simulation Numérique

 CHAPITRE 2

Formulations Faibles

This is the problem of all great revelations : their significance so often exceeds the frame of our comprehen-
sion. We understand only after, always after. Not simply when it is too late, but precisely because it is too
late.

---R. Scott Bakker - The Darkness That Comes Before : Book 1 of the Prince of Nothing (Titre français :
Autrefois, les Ténèbres),

2.1 Espaces de Hilbert : rappels

Definition 2.1
Soit 𝑉 un R− espace vectoriel, alors l’application (·, ·) : 𝑉 × 𝑉 → R est un produit scalaire si et seulement si elle vérifie,
pour tout x, y, z ∈ 𝑉 et tout scalaire 𝛼 ∈ R :
1. (x, y) = (y, x)
2. (x + y, z) = (x, z) + (y, z)
3. (𝛼x, y) = 𝛼 (x, y)
4. (x, x) ∈ R+
5. (x, x) = 0 =⇒ x = 0

Definition 2.2
Un R− espace vectoriel 𝑉 est dit pré-Hilbertien si il est muni d’un produit scalaire.

Definition 2.3
Soit 𝑉 un R− espace vectoriel, alors l’application ‖·‖ : 𝑉 → R est une norme si et seulement si elle vérifie, pour tout
x, y ∈ 𝑉 et tout scalaire 𝛼 ∈ R :

11
Maillage et Éléments Finis

1. Séparation : ‖x‖ = 0 =⇒ 𝑥 = 0
2. Absolue homogénéité : ‖𝛼x‖ = |𝛼| ‖x‖
3. Inégalité triangulaire : ‖x + y‖ ≤ ‖x‖ + ‖y‖

Remark 2.1
Un produit scalaire induit une norme sur un espace de Hilbert :
√︀
‖x‖ := (x, x).

Nous rappelons l’inégalité de Cauchy Schwarz :

Proposition 2.1 (Inégalité de Cauchy Schwarz)

Pour tout x et y appartenant à un espace pré-Hilbertien 𝑉 :

|(x, y)| ≤ ‖x‖ ‖y‖ .

Definition 2.4
Un espace pré-Hilbertien 𝑉 est un espace de Hilbert si et seulement si il est complet pour la norme ‖·‖ induite par son
produit scalaire.

Definition 2.5
Soit 𝑉 un espace de Hilbert. L’application 𝑓 : 𝑉 × 𝑉 → R est une forme bilinéaire sur 𝑉 si et seulement si, pour tout
x, y, z de 𝑉 et 𝛼 de R :
1. 𝑓 (x, y + 𝛼z) = 𝑓 (x, y) + 𝛼𝑓 (x, z)
2. 𝑓 (𝛼x + y, z) = 𝛼𝑓 (x, z) + 𝑓 (y, z)

Theorem 2.1 (Représentation de Riesz)

Soit 𝑉 un espace de Hilbert de produit scalaire (·, ·) et de norme induite ‖·‖. Pour toute forme anti-linéaire continue ℓ, il
existe un unique 𝑤 ∈ 𝑉 tel que

ℓ(𝑣) = (𝑤, 𝑣) , ∀𝑣 ∈ 𝑉.

De plus, nous avons

|ℓ(𝑣)|
‖𝑤‖ = sup .
𝑣∈𝑉 ∖{0} ‖𝑣‖

Remark 2.2
Ce théorème montre que la forme ℓ peut être représentée par un vecteur 𝑤 qui est unique. Autrement dit, peu importe
𝑣, la quantité ℓ(𝑣) peut se calculer par la seule connaissance du vecteur 𝑤 et d’un « simple » produit scalaire.

12 Chapitre 2. Formulations Faibles

 Maillage et Éléments Finis

2.2 Formulation Faible

2.2.1 Domaine Physique

ΓD
ΓN

Fıg. 1 – Exemple de domaine de calcul avec sa normal unitaire sortante

Dans ce cours, nous considérons un ouvert polygonal Ω de R2 . Nous restons en dimension 2 pour plus de facilité mais
l’extension à la dimension 3 est relativement directe et tous les résultats énoncés sont aussi valable en dimension 3. Sur
chaque segment du bord Γ := 𝜕Ω du domaine, on définit le vecteur unitaire normale n sortant à Ω. Nous noterons que
ce vecteur n’existe pas aux intersections entre les segments. Le domaine Ω est supposé ne pas comporter de fissure ni
de point de rebroussement. Son bord est divisé en deux parties distinctes : Γ𝐷 et Γ𝑁 , potentiellement non connexe mais
d’intersection vide : Γ = Γ𝐷 ∪ Γ𝑁 et Γ𝐷 ∩ Γ𝑁 = ∅. Selon la partie du bord, une condition sera imposée à la solution :
— Sur Γ𝐷 : condition de Dirichlet, c’est à dire que la valeur de la solution y est imposée (eg 𝑢 = 0). En mécanique
on dirait que le déplacement est imposée.
— Sur Γ𝑁 : condition de Neumann, c’est à dire que le flux de la solution y est imposée (eg 𝜕n 𝑢 = 0). En mécanique,
on dirait que la force normale est imposée.

Remark
En général, on préfère travailler avec des ouverts réguliers, de classe au moins 𝒞 1 . Un tel ouvert présente l’avantage de
pouvoir clairement définir le vecteur unitaire normale n sortante à Ω. Cependant, après maillage, on se retrouve avec…
un polygone ! Alors plutôt que de travailler dans un domaine régulier pour après le casser en (petits) morceaux, nous
préférons ici commencer directement avec un polygone et mettre l’accent sur les algorithmes et la mise en oeuvre de la
méthode que les spécificités mathématiques.

Remark
Un point x ∈ R2 est parfois noté x = (𝑥, 𝑦) ou x = (𝑥1 , 𝑥2 ) selon les besoins. Nous commencerons les indiçages par 1
bien qu’en informatique cela commence souvent par 0.

2.2. Formulation Faible 13

Maillage et Éléments Finis

2.2.2 EDP (problème fort)

Ce cours se concentre sur les équations aux dérivées partielles (EDP) elliptiques du second ordre, qui font appel à l’opé-
rateur de Laplace 1 (ou Laplacien) :

𝜕2 𝜕2 𝜕2 𝜕2
Δ := + = + .
𝜕𝑥 𝜕𝑦 𝜕𝑥1 𝜕𝑥2
Nous considérons le problème générique suivant, appelé aussi problème de réaction-diffusion :
{︂
−Δ𝑢 + 𝑢 = 𝑓 (Ω),
(1)
𝜕n 𝑢 = 0 (Γ),

où nous avons défini :

— le terme 𝜕n 𝑢 désigne la dérivée normale de 𝑢 sur le bord, c’est à dire la dérivée de 𝑢 dans la direction n : 𝜕n 𝑢 =
(∇𝑢) · n, avec ∇𝑢 = [𝜕𝑥1 𝑢, 𝜕𝑥2 𝑢]𝑇 son vecteur gradient. Vous aurez sans doute remarqué que, entre deux arêtes,
le vecteur normal n n’est pas défini et donc la dérivée normale non plus. Ce « problème » n’en est pas vraiment
un et pour l’instant mettez cela de côté, nous y reviendrons !
— (−Δ𝑢) : Terme de diffusion (notez le signe négatif)
— 𝜕n 𝑢 = 0 : Condition de Neumann homogène
— 𝑓 : une fonction donnée définie sur Ω. Elle joue le rôle de terme source, c’est à dire d’apport (positif ou négatif),
par exemple de chaleur ou de force surfacique.
Pour le moment, nous ne nous intéressons pas à la régularité de la solution ni même à l’existence et l’unicité de celle-ci : nous
supposons que le problème (1) est bien posé. Une fois la méthode des éléments finis apréhendée, nous nous intéresserons
à ces questions. Cela va à l’encontre de l’habitude en mathématiques où l’on démontre le caractère bien posé avant de s’y
attaquer. Faites moi confiance et tout s’éclairera !

Remark
Pour l’instant, nous imposons une condition aux bords que nous imposons, de type Neumann homogène. Plus tard nous
verrons d’autres types de conditions : Dirichlet, où la valeur de solution est imposée, et Fourier, un mélange entre Dirichlet
et Neumann.

2.2.3 Théorème de Green

Pour une géométrie arbitraire, nous ne savons pas, en général, obtenir la solution forte (ou classique) du problème (1).
La méthode des éléments finis se base sur l’approximation numérique de la solution au sens faible du problème (1).
Nous verrons qu’une solution faible est, en fait et en général, forte. Commençons par réécrire le problème d’origine sous
sa formulation faible ou formulation variationnelle.
Un théorème central dans l’analyse des EDP et qui permet d’obtenir ces formulations faibles est celui de Green 2
∫︁ ∫︁ ∫︁
∀𝑢, 𝑣, (Δ𝑢)(x)𝑣(x)dx = − ∇𝑢(x) · ∇𝑣(x)dx + (𝜕n 𝑢)(x)𝑣(x)d𝑠(x).
Ω Ω Γ

Ce résultat est également valable en dimension 3 pour des domaines polygonaux (ouf). Le produit ∇𝑢 · ∇𝑣 est le produit
scalaire euclidien standard. La quantité 𝑣 est ici laissé non définie, c’est normal : supposez que c’est une fonction de même
régularité que 𝑢. Pour compacter les équations, nous n’indiquerons plus les quantités d’intégrations :
∫︁ ∫︁ ∫︁
∀𝑢, 𝑣, (Δ𝑢)𝑣 = − ∇𝑢 · ∇𝑣 + (𝜕n 𝑢)𝑣.
Ω Ω Γ

1. Pierre-Simon Laplace (1749 – 1827).

2. George Green (1793 – 1841). Mathématicien britannique (quasi) auto-didacte.

14 Chapitre 2. Formulations Faibles

 Maillage et Éléments Finis

Remark
Ce résultat est en quelque sorte une extension multi-dimensionnel de l’intégration par partie sur un segment Ω = [𝑎, 𝑏].
En effet, en dimension 1, l’opérateur Δ devient la dérivée seconde. La normale sortante au segment devient un scalaire
valant -1 « à gauche » (en 𝑎) et 1 « à droite » (en 𝑏) et la dérivée normale devient 𝜕n 𝑢 = ±𝑢′ :
∫︁ 𝑏 ∫︁ 𝑏 ∫︁ 𝑏
′′ ′ ′ ′ ′
𝑢 𝑣=− 𝑢 𝑣 + 𝑢 (𝑏)𝑣(𝑏) − 𝑢 (𝑎)𝑣(𝑎) = − 𝑢′ 𝑣 ′ + 𝜕n 𝑢(𝑏)𝑣(𝑏) + 𝜕n 𝑢(𝑎)𝑣(𝑎).
𝑎 𝑎 𝑎

2.2.4 Formulation faible

Le point de départ de notre analyse est la réécriture sous forme faible du problème (1). Pour cela, la méthode consiste à :
— Multiplier l’EDP par une fonction test 𝑣
— Intéger le tout sur Ω
— Appliquer le Théorème de Green
— Appliquer les conditions aux bords ∫︁ ∫︁ ∫︁
−Δ𝑢 + 𝑢 = 𝑓 =⇒ − Δ𝑢𝑣 + 𝑢𝑣 = 𝑓𝑣
∫︁ Ω ∫︁Ω Ω
∫︁ ∫︁
=⇒ ∇𝑢 · ∇ − (𝜕n 𝑢) 𝑣 + 𝑐 𝑢𝑣 = 𝑓𝑣
Ω Γ⏟ ⏞ Ω Ω
=0

Nous obtenons alors :

∫︁ ∫︁ ∫︁
−Δ𝑢 + 𝑢 = 𝑓 (Ω) =⇒ ∀𝑣, ∇𝑢 · ∇𝑣 + 𝑢𝑣 = 𝑓 𝑣. (2)
Ω Ω Ω

Ainsi, et toujours sans regarder la régularité de 𝑢 (ni de 𝑣), nous avons que : si 𝑢 est solution de l’EDP (1) alors 𝑢 est aussi
solution de la formulation faible :

⎨Trouver 𝑢 tel que

⎧
∫︁ ∫︁ ∫︁
⎩∀𝑣, ∇𝑢 · ∇𝑣 + 𝑢𝑣 = 𝑓 𝑣.
Ω Ω Ω

À gauche du signe égal se trouve l’inconnue (𝑢) et à droite la donnée (𝑓 ), c’est une convention et plus tard cette équation
s’écrira sous la forme d’un système linéaire 𝐴𝑈 = 𝐵 où le vecteur 𝐵 correspondra au membre de droite de (2) et la
matrice 𝐴 à la partie de gauche.

Remark
Attention, sur le bord Γ, 𝜕n 𝑢 = 0 n’implique pas 𝑢 = 0 !

Nous pouvons maintenant définir plus proprement la quantité 𝑣. Appelée fonction test elle n’a d’autre rôle que de « tester »
la solution. L’idée de la formulation faible est de chercher une solution qui vérifie l’EDP, non pas point à point (au sens
fort, donc) mais « en moyenne », via l’intégrale. En mécanique, 𝑣 est appelé « travaux virtuels » (avec la méthode éponyme
qui est, en fait, la formulation faible) : cette quantité est arbitraire et n’est utile que pour écrire le problème faible (2).

2.2. Formulation Faible 15

Maillage et Éléments Finis

2.3 Théorème de Lax-Milgram

Nous pouvons maintenant énoncer le théorème de Lax-Milgram (à connaître par cœur).

Theorem 2.2 (de Lax-Milgram)

Soit 𝑉 un espace de Hilbert de produit scalaire (·, ·) et de norme ‖·‖, et soit la formulation faible suivante

Trouver 𝑢 ∈ 𝑉 tel que,

{︂
(3)
∀𝑣 ∈ 𝑉, 𝑎(𝑢, 𝑣) = ℓ(𝑣).

Sous réserve des quatre hypothèses suivantes :

1. ℓ est une forme linéaire continue sur 𝑉 :
∃𝐶 > 0/∀𝑣 ∈ 𝑉, |ℓ(𝑣)| ≤ 𝐶 ‖𝑣‖ .
2. 𝑎(·, ·) est une forme biquilinéaire sur 𝑉 × 𝑉 .
3. 𝑎(·, ·) est continue :
∃𝑀 > 0/∀(𝑢, 𝑣) ∈ 𝑉 × 𝑉, |𝑎(𝑢, 𝑣)| ≤ 𝑀 ‖𝑢‖ ‖𝑣‖ .
4. 𝑎(·, ·) est coercive (ou elliptique) :
2
∃𝛼 > 0/∀𝑢 ∈ 𝑉, 𝑎(𝑢, 𝑢) ≥ 𝛼 ‖𝑢‖ .
Alors la formulation faible (3) admet une unique solution. De plus cette solution dépend continûment de la forme linéaire
ℓ:
𝑀
‖𝑢‖ ≤ 𝐶
𝛼

Proof. Comme il est question de forme linéaire, nous allons utiliser le Théorème de représentation de Riesz. En effet,
pour tout 𝑤 de 𝑉 , l’application 𝑣 → 𝑎(𝑤, 𝑣) est anti-linéaire et continue de 𝑉 dans R. Il existe donc un unique élément
de 𝑉 , noté 𝐴(𝑤) (Théorème de Riesz), tel que

∀𝑣 ∈ 𝑉, 𝑎(𝑤, 𝑣) = (𝐴(𝑤), 𝑣) .

Nous allons montrer que l’opérateur 𝐴 : 𝑉 → 𝑉 est continue, inversible et d’inverse continu. L’opérateur 𝐴 est clairement
linéaire. En prenant 𝑣 = 𝐴(𝑤) et en utilisant la continuité de 𝑎(·, ·), nous obtenons :

‖𝐴(𝑤)‖2 = (𝐴(𝑤), 𝐴(𝑤)) = 𝑎(𝑤, 𝐴(𝑤)) ≤ 𝑀 ‖𝑤‖‖𝐴(𝑤)‖.

Cette relation étant valable pour tout 𝑤, elle signifie que 𝐴 est continue, puisque :

∀𝑤 ∈ 𝑉, ‖𝐴(𝑤)‖ ≤ 𝑀 ‖𝑤‖.

Appliquons de nouveau le Théorème de représentation de Riesz au membre de droite, puisque ℓ est une forme anti-linéaire
continue :

∃!𝑓 ∈ 𝑉 tel que ‖𝑓 ‖ = ‖ℓ‖𝑉 ′ et ∀𝑣 ∈ 𝑉, ℓ(𝑣) = (𝑓, 𝑣) .

Comme 𝐴(𝑢) = 𝑓 est équivalent à ∀𝑣 ∈ 𝑉, (𝐴(𝑢), 𝑣) = (𝑓, 𝑣), alors notre formulation faible (3) devient équivalent au
problème linéaire :

Trouver 𝑢 ∈ 𝑉 tel que

{︂

𝐴(𝑢) = 𝑓.

16 Chapitre 2. Formulations Faibles

 Maillage et Éléments Finis

La question est : 𝐴 est-elle bijective ? Utilisons la coercivité de l’application 𝑎(·, ·) :

𝛼‖𝑤‖2 ≤ |(𝐴(𝑤), 𝑤)| ≤ ‖𝐴(𝑤)‖‖𝑤‖,

ce qui implique que

𝛼‖𝑤‖ ≤ ‖𝐴(𝑤)‖. (4)

Comme 𝛼 > 0, alors 𝐴 est injective. En dimension finie et comme 𝐴 est un endomorphisme, nous pourrions en déduire la
surjectivité de 𝐴. Mais nous sommes malheureusement en dimension infinie, nous devons donc montrer que Im(𝐴) = 𝑉 ,
pour cela nous montrons que Im(𝐴) est fermé dans 𝑉 et que son orthogonal (dans 𝑉 ) est réduit au singleton nul. Prenons
une suite (𝐴(𝑤𝑛 ))𝑛 de Im(𝐴) qui converge dans 𝑉 . Nous avons, pour tout 𝑛, 𝑝 ∈ N et grâce à (4),

𝛼‖𝑤𝑛 − 𝑤𝑝 ‖ ≤ ‖𝐴(𝑤𝑛 ) − 𝐴(𝑤𝑝 )‖.

Quand 𝑛 et 𝑝 tendent vers l’infini, alors ‖𝑤𝑛 − 𝑤𝑝 ‖ → 0. La suite (𝑤𝑛 )𝑛 est donc une suite de Cauchy dans 𝑉 , qui est
complet (Hilbert), elle est donc convergente et converge vers un élément 𝑤 de 𝑉 . Par continuité de 𝐴, la suite (𝐴(𝑤𝑛 ))𝑛
converge vers 𝐴(𝑤), élément de Im(𝐴). Ce qui implique que Im(𝐴) est fermé. Prenons maintenant 𝑣 ∈ Im(𝐴)⊥ , par la
coercivité de 𝑎(·, ·), nous avons
2
𝛼 |𝑣| ≤ |𝑎(𝑣, 𝑣)| = |(𝐴(𝑣), 𝑣)| = 0.

Autrement dit, 𝑣 = 0 et donc Im(𝐴)⊥ = {0} et nous avons

)︀⊥
Im(𝐴) = Im(𝐴) = Im(𝐴)⊥ = {0}⊥ = 𝑉.
(︀

L’application 𝐴 est donc bijective. Son inverse 𝐴−1 existe, et, avec (4), nous obtenons sa continuité :

⃦𝐴 (𝑤)⃦ ≤ 1 ‖𝑤‖ .
⃦ −1 ⃦
∀𝑤 ∈ 𝑉,
𝛼
Ceci prouve que 𝑢 dépend continûment du membre de droite 𝑓 (qui dépend de ℓ).

Remark 2.8
À quoi sert ce théorème ? Sous réserve de 4 hypothèses, nous aurons la garantie que la formulation faible obtenue précé-
demment admet une solution (ce qui est bien) et que cette solution est unique (encore mieux !). Il est donc d’une importance
capitale.

Avant de pouvoir appliquer ce théorème proprement dit, nous devons connaître un peu mieux les espaces de Sobolev :
Hilbert ? Norme ? Tant de questions.

2.4 Espaces de Sobolev

2.4.1 Espace 𝐿2 (Ω)

Rappelons que l’espace 𝐿2 (Ω) est l’espace des fonctions de carré mesurable (au sens de Lebesgue). Muni du produit
scalaire
∫︁
(𝑓, 𝑔)𝐿2 (Ω) = 𝑓 (x)𝑔(x)dx,
Ω

2.4. Espaces de Sobolev 17

Maillage et Éléments Finis

l’espace 𝐿2 (Ω) est un espace de Hilbert, de norme induite :

(︂∫︁ )︂ 12
2
‖𝑓 ‖𝐿2 (Ω) = |𝑓 (x)| dx .
Ω

Il est important de remarquer qu’une fonction de 𝐿2 (Ω) est définie presque partout. Autrement dit, deux fonctions 𝑓 et
𝑔 de 𝐿2 (Ω) peuvent être égales tout en ayant des valeurs différentes sur un sous-ensemble 𝜔 de Ω, de mesure nulle. Une
fonction de 𝐿2 (Ω) désigne en réalité une classe de fonctions.
Introduisons l’espace C𝑐∞ (Ω) des fonctions C ∞ sur Ω à support compact dans Ω :

C𝑐∞ (Ω) := {𝑓 ∈ C ∞ (Ω) | supp(𝑓 ) est compact dans Ω} .

Remark 2.9
— Ces fonctions (et toutes leurs dérivées) s’annulent nécesairement sur le bord de Ω (qui est ouvert)
— Dans notre cas, Ω ∈ R2 (ou même R3 ), ce qui impliqué qu’un compact de Ω est donc un fermé borné - Un
exemple d’une telle fonction est la fonction « blob » comme illustré par l’article Wikipédia

Nous rappelons/admettons le théorème de densité suivant.

Theorem 2.3 (Densité dans 𝐿2 (Ω))

L’ensemble C𝑐∞ (Ω) est dense dans 𝐿2 (Ω).

Autrement dit, pour tout élément 𝑓 de 𝐿2 (Ω), il existe une suite (𝑓𝑛 )𝑛 de fonctions de C𝑐∞ (Ω) qui converge vers 𝑓 pour
la norme de 𝐿2 (Ω). Ce théorème est extrêmement important : pour démontrer des propriétés de 𝐿2 (Ω), nous utiliserons
des propriétés de C𝑐∞ (Ω) et passerons à la limite dans 𝐿2 (Ω).

Corollary 2.1
Soit 𝑓 une fonction de 𝐿2 (Ω) telle que
∫︁
∀𝜑 ∈ C𝑐∞ (Ω), 𝑓 (x)𝜑(x)dx = (𝑓, 𝜑)𝐿2 (Ω) = 0,
Ω

alors 𝑓 (x) = 0 presque partout dans Ω.

Proof. D’après le théorème ??, il existe une suite (𝑓𝑛 )𝑛 de C𝑐∞ (Ω) qui converge vers 𝑓 . Nous avons alors
∫︁ ∫︁
2 2
0 = lim 𝑓 (x)𝑓𝑛 (x)dx = |𝑓 (x)| dx = ‖𝑓 ‖𝐿2 (Ω) ,
𝑛→∞ Ω Ω

d’où 𝑓 est nulle « au sens de » 𝐿2 (Ω), c’est-à-dire que 𝑓 (x) = 0 presque partout.

L’espace 𝐿2 (Ω) est un « petit » espace de Hilbert qui contient C 1 (Ω). Nous nous rapprochons du but… Cependant les
fonctions de 𝐿2 (Ω) ne sont pas dérivables ! Elles ne sont donc pas utilisables en pratique dans les formulations faibles.
C’est tout l’objet de la section suivante : proposer une nouvelle forme de dérivation plus faible, c’est-à-dire ici, qui ne
requiert pas de continuité.

18 Chapitre 2. Formulations Faibles

 Maillage et Éléments Finis

2.4.2 Dérivée faible

Definition 2.6
Une fonction de 𝐿2 (Ω) est dérivable au sens faible par rapport à la direction 𝑥𝑖 si et seulement si il existe un élément 𝑔𝑖
de 𝐿2 (Ω) tel que
∫︁ ∫︁
∀𝜑 ∈ C𝑐∞ (Ω), 𝑓 (x)𝜕𝑥𝑖 𝜑(x)dx = − 𝑔𝑖 (x)𝜑(x)dx.
Ω Ω

Nous notons alors 𝑔𝑖 = 𝜕𝑥𝑖 𝑓 = 𝜕𝑖 𝑓 , qui est unique en vertu du Corollaire 2.4.1.

Nous noterons maintenant 𝜕𝑥𝑖 𝑓 ∈ 𝐿2 (Ω) ou 𝜕𝑖 𝑓 ∈ 𝐿2 (Ω) pour signifier que 𝑓 est dérivable au sens faible par rapport à
𝑥𝑖 . De la même manière, nous pouvons définir le gradient faible :

Definition 2.7
Une fonction 𝑓 ∈ 𝐿2 (Ω) admet un gradient faible, noté ∇𝑓 , si et seulement si 𝑓 est dérivable au sens faible par rapport
à toutes ses variables, et nous avons alors
𝑇
∇𝑓 = (𝜕𝑥1 𝑓, 𝜕𝑥2 𝑓, . . . , 𝜕𝑥𝑑 𝑓 ) .

Proposition 2.2 (Unicité de la dérivée faible)

Si la dérivée faible dans la direction 𝑥𝑖 d’une fonction 𝑢 ∈ 𝐿2 (Ω) existe alors elle est unique.

Proof. Soit 𝑢 ∈ 𝐿2 (Ω) et supposons que 𝑢 admette deux dérivées faibles dans la direction 𝑥𝑖 : 𝑓𝑖 et 𝑔𝑖 , toutes deux dans
𝐿2 (Ω). Nous avons alors, par définition, les deux relations suivantes :
∫︁ ∫︁ ∫︁
∀𝜑 ∈ C𝑐∞ (Ω), − 𝑢(x)𝜕𝑖 𝜑(x) = 𝑓𝑖 (x)𝜑(x)dx = 𝑔𝑖 (x)𝜑(x)dx.
Ω Ω Ω

Autrement dit, nous avons

∫︁
∀𝜑 ∈ C𝑐∞ (Ω), (𝑓𝑖 − 𝑔𝑖 )(x)𝜑(x)dx = 0,
Ω

et le Corollaire 2.4.1 implique que 𝑓𝑖 = 𝑔𝑖 .

Le lien entre dérivée faible et dérivée forte (ou classique) est maintenant présenté :

Proposition 2.3
Soit 𝑢 ∈ C 1 (Ω) tel que son gradient, au sens classique, ∇𝑢 soit dans C 0 (Ω), alors 𝑢 admet un gradient au sens faible
̃︀ et l’on a ∇𝑢 = ∇𝑢.
∇𝑢 ̃︀

Proof. Il suffit de montrer ce résultat pour une direction uniquement, c’est-à-dire montrer que 𝜕̃︀𝑖 𝑢 = 𝜕𝑖 𝑢, si 𝜕̃︀𝑖 est la
dérivée partielle au sens faible. Par intégration par partie, nous avons :
∫︁ ∫︁ ∫︁
∀𝜑 ∈ C𝑐∞ (Ω), 𝜕𝑖 𝑢(x)𝜑(x)dx = − 𝑢(x)𝜕𝑖 𝜑(x)dx + 𝑢(x)𝜑(x)𝑛𝑖 (x)d𝑠(x),
Ω Ω 𝜕Ω

2.4. Espaces de Sobolev 19

Maillage et Éléments Finis

où 𝑛𝑖 est la 𝑖ème composante du vecteur normale n. Comme 𝜑 est à support compact dans Ω, nous savons que 𝜑 s’annule
sur le bord de Ω. Il vient donc
∫︁ ∫︁
∀𝜑 ∈ C𝑐∞ (Ω), 𝜕𝑖 𝑢(x)𝜑(x)dx = − 𝑢(x)𝜕𝑖 𝜑(x)dx,
Ω Ω

Autrement dit 𝑢 admet une dérivée faible. Celle-ci étant unique, nous avons bien 𝜕𝑖 𝑢 = 𝜕̃︀𝑖 𝑢.

Remark 2.10
Dans la suite, puisque nous ne travaillerons qu’avec des dérivées partielles faibles, nous omettrons le tilde.

2.4.3 Espace de Sobolev 𝐻 1 (Ω)

Nous disposons maintenant des outils nécessaires pour introduire l’espace de Sobolev 𝐻 1 (Ω) des fonctions de carré
intégrable et dérivables au sens faible dans chaque direction (𝑑 = 2, 3 est la dimension) :

𝐻 1 (Ω) = 𝑢 ∈ 𝐿2 (Ω) ⃒ ∇𝑢 ∈ (𝐿2 (Ω))𝑑 .

{︀ ⃒ }︀

Nous munissons cet espace du produit scalaire suivant (pour 𝑢 et 𝑣 dans 𝐻 1 (Ω))
∫︁ ∫︁
(𝑢, 𝑣)𝐻 1 (Ω) = 𝑢(x)𝑣(x) dx + ∇𝑢(x) · ∇𝑣(x) dx,
Ω Ω

et de la norme induite, pour 𝑢 ∈ 𝐻 1 (Ω) :

(︂∫︁ ∫︁ )︂ 12
‖𝑢‖𝐻 1 (Ω) = |𝑢(x)| dx +
2
|∇𝑢(x)| dx2
.
Ω Ω

Remark 2.11
Nous pouvons montrer que c’est effectivement un produit scalaire avec les arguments similaires à ceux utilisés pour montrer
que la « même » application est un produit scalaire sur C 1 (Ω).

Remark 2.12
Pour 𝑢 de 𝐻 1 (Ω), nous avons clairement
𝑑
∑︁
2 2 2
‖𝑢‖𝐻 1 (Ω) = ‖𝑢‖𝐿2 (Ω) + ‖𝜕𝑖 𝑢‖𝐿2 (Ω)
𝑖=1
2 2
= ‖𝑢‖𝐿2 (Ω) + ‖∇𝑢‖(𝐿2 (Ω))𝑑 ,

et donc les inégalités suivantes :

1. ‖𝑢‖𝐻 1 (Ω) ≥ ‖𝑢‖𝐿2 (Ω)
(︁∑︀ )︁1/2
𝑑 2
2. ‖𝑢‖𝐻 1 (Ω) ≥ ‖∇𝑢‖(𝐿2 (Ω))𝑑 = 𝑖=1 ‖𝜕𝑖 𝑢‖𝐿2 (Ω)
3. ‖𝑢‖𝐻 1 (Ω) ≥ ‖𝜕𝑖 𝑢‖𝐿2 (Ω) ∀𝑖 = 1, 2, . . . , 𝑑

20 Chapitre 2. Formulations Faibles

 Maillage et Éléments Finis

Nous montrons maintenant que 𝐻 1 (Ω) muni de cette norme est complet.

Theorem 2.4 (Complétude de 𝐻 1 (Ω))

L’espace 𝐻 1 (Ω) est complet pour la norme ‖·‖𝐻 1 (Ω) .

Proof. Prenons une suite de Cauchy (𝑢𝑛 )𝑛 de 𝐻 1 (Ω) et montrons qu’elle converge dans 𝐻 1 (Ω). Par définition de la suite
de Cauchy, nous avons

∀𝜀 > 0, ∃𝑁 > 0 tel que ∀𝑛 > 𝑁, ∀𝑝 > 𝑁, ‖𝑢𝑛 − 𝑢𝑝 ‖𝐻 1 (Ω) ≤ 𝜀.

Par ailleurs, pour 𝑛, 𝑝 de N l’inégalité suivante est vérifiée :

‖𝑢𝑛 − 𝑢𝑝 ‖𝐿2 (Ω) ≤ ‖𝑢𝑛 − 𝑢𝑝 ‖𝐻 1 (Ω) ,

ce qui fait de la suite (𝑢𝑛 )𝑛 une suite de Cauchy dans 𝐿2 (Ω), puisque :

∀𝜀 > 0, ∃𝑁 > 0 tel que ∀𝑛 > 𝑁, ∀𝑝 > 𝑁, ‖𝑢𝑛 − 𝑢𝑝 ‖𝐿2 (Ω) ≤ ‖𝑢𝑛 − 𝑢𝑝 ‖𝐻 1 (Ω) ≤ 𝜀.

L’espace 𝐿2 (Ω) étant complet, la suite (𝑢𝑛 )𝑛 converge dans 𝐿2 (Ω) vers 𝑢 ∈ 𝐿2 (Ω). Nous appliquons le même raison-
nement aux dérivées partielles : pour 𝑖 = 1, . . . , 𝑑, nous avons aussi

‖𝜕𝑖 𝑢𝑛 − 𝜕𝑖 𝑢𝑝 ‖𝐿2 (Ω) ≤ ‖𝑢𝑛 − 𝑢𝑝 ‖𝐻 1 (Ω) .

Ainsi, pour tout 𝑖, la suite (𝜕𝑖 𝑢𝑛 )𝑛 est aussi de Cauchy dans 𝐿2 (Ω) et converge donc vers un élément 𝑓𝑖 ∈ 𝐿2 (Ω). Il nous
faut donc montrer que 𝑢 est dérivable (au sens faible) et que 𝑓𝑖 = 𝜕𝑖 𝑢. Remarquons pour cela que, par définition,
∫︁ ∫︁
∀𝜑 ∈ C𝑐 (Ω),
∞
𝜕𝑖 𝑢𝑛 (x)𝜑(x)dx = − 𝑢𝑛 (x)𝜕𝑖 𝜑(x)dx.
Ω Ω

En passant à la limite dans 𝐿2 (Ω) dans cette expression, il vient que :

∫︁ ∫︁
∀𝜑 ∈ C𝑐∞ (Ω), 𝑓𝑖 (x)𝜑(x)dx = − 𝑢(x)𝜕𝑖 𝜑(x)dx.
Ω Ω

Autrement dit, 𝑢 est dérivable par rapport à toutes ses variables et 𝜕𝑖 𝑢 = 𝑓𝑖 , ce qui implique que 𝑢 est bien dans 𝐻 1 (Ω).
Nous avons donc montré que la suite (𝑢𝑛 )𝑛 converge dans 𝐿2 (Ω) vers un élément 𝑢 de 𝐻 1 (Ω). Il nous reste à montrer
que cette convergence est toujours valable pour la norme de 𝐻 1 (Ω). Utilisons la remarque précédente pour décomposer
la norme dans 𝐻 1 (Ω) :
𝑑
∑︁
2 2 2
‖𝑢𝑛 − 𝑢‖𝐻 1 (Ω) = ‖𝑢𝑛 − 𝑢‖𝐿2 (Ω) + ‖𝜕𝑗 𝑢𝑛 − 𝜕𝑗 𝑢‖𝐿2 (Ω) → 0(𝑛 → +∞).
𝑗=1

La suite de Cauchy (𝑢𝑛 )𝑛 est donc convergente dans 𝐻 1 (Ω), ce dernier est donc complet.

Nous en déduisons le corollaire suivant :

Corollary 2.2
𝐻 1 (Ω) est un espace de Hilbert pour le produit scalaire (·, ·)𝐻 1 (Ω) .

Nous avons également le résultat de densité suivant

2.4. Espaces de Sobolev 21

Maillage et Éléments Finis

Proposition 2.4 (Admis)

L’espace C𝑐∞ (Ω) est dense dans 𝐻 1 (Ω) pour la norme ‖·‖𝐻 1 (Ω) .

En particulier, l’espace C 1 (Ω), qui contient C𝑐∞ (Ω), est dense dans 𝐻 1 (Ω) pour la norme ‖·‖𝐻 1 (Ω) . Ce résultat nous dit
que 𝐻 1 (Ω) est le « plus petit » espace complet contenant C 1 (Ω) : c’est ce que nous cherchions !

2.5 Application au problème modèle

2.5.1 Formulation faible

Considérons un ouvert polygonal connexe Ω et le problème suivant

{︂
−Δ𝑢 + 𝑐𝑢 = 𝑓 (Ω)
(5)
𝜕n 𝑢 = 0 (Γ𝑁 = Γ)

Après mutiplication par des fonctions tests et intégrations par partie, nous obtenons la formulation faible de ce problème.

Trouver 𝑢 ∈ 𝐻 1 (Ω) tel que

{︂
(6)
∀𝑣 ∈ 𝐻 1 (Ω), 𝑎(𝑢, 𝑣) = ℓ(𝑣)

avec 𝑎(·, ·) : 𝐻 1 (Ω) × 𝐻 1 (Ω) → R et ℓ(·) : 𝐻 1 (Ω) → R définies pas

⎧ ∫︁ ∫︁
⎪
⎪
⎨ 𝑎(𝑢, 𝑣) = ∇𝑢 · ∇𝑣 + 𝑐 𝑢𝑣
∫︁Ω Ω
⎪
⎩ ℓ(𝑣) =
⎪ 𝑓𝑣
Ω

2.5.2 Existence et unicité

Tentons d’appliquer le théorème de Lax-Milgram à cette formulation faible

1. 𝐻 1 (Ω) est un espace de Hilbert
2. ℓ(·) est clairement linéaire (du fait de l’intégrale)
3. 𝑎(·, ·) est bilinéaire, pour la même raison
4. Continuité de ℓ(·) : prenons une fonction 𝑣 ∈ 𝐻 1 (Ω) :
⃒∫︁ ⃒
⃒ ⃒
|ℓ(𝑣)| = ⃒ 𝑓 𝑣 ⃒⃒
⃒
⏟ Ω⏞
(𝑓,𝑣)𝐿2 (Ω)

≤ ‖𝑓 ‖𝐿2 (Ω) ‖𝑣‖𝐿2 (Ω) Cauchy-Schwarz

≤ ‖𝑓 ‖𝐿2 (Ω) ‖𝑣‖𝐻 1 (Ω) inégalité des normes
⏟ ⏞
Constant

5. Continuité de 𝑎(·, ·) : prenons deux fonctions 𝑢 et 𝑣 de 𝐻 1 (Ω) :

22 Chapitre 2. Formulations Faibles

 Maillage et Éléments Finis

⃒∫︁ ∫︁ ⃒
⃒ ⃒
|𝑎(𝑢, 𝑣)| = ⃒⃒ ∇𝑢 · ∇𝑣 + 𝑐 𝑢𝑣 ⃒⃒
⃒∫︁Ω ⃒ Ω
⃒∫︁ ⃒
inégalité classique
⃒ ⃒ ⃒ ⃒
≤ ⃒⃒ ∇𝑢 · ∇𝑣 ⃒⃒ + |𝑐| ⃒⃒ 𝑢𝑣 ⃒⃒
⏟ Ω ⏞ ⏟ Ω⏞
(∇𝑢,∇𝑣)𝐿2 (Ω)𝑑 (𝑢,𝑣)𝐿2 (Ω)

≤ ‖∇𝑢‖(𝐿2 (Ω))𝑑 ‖∇𝑣‖(𝐿2 (Ω))𝑑 + |𝑐| ‖𝑢‖𝐿2 (Ω) ‖𝑣‖𝐿2 (Ω) inégalité triangulaire dans𝐿2 (Ω)
≤ ‖𝑢‖𝐻 1 (Ω) ‖𝑣‖𝐻 1 (Ω) + |𝑐| ‖𝑢‖𝐻 1 (Ω) ‖𝑣‖𝐻 1 (Ω) inégalité des normes
≤ (1 + 𝑐) ‖𝑢‖𝐻 1 (Ω) ‖𝑣‖𝐻 1 (Ω)
6. Coercivité de 𝑎(·, ·) : prenons une fonction 𝑢 ∈ 𝐻 1 (Ω) :
∫︁ ∫︁ ∫︁ ∫︁
𝑎(𝑢, 𝑢) = ∇𝑢 · ∇𝑢 + 𝑐 𝑢𝑢 = ‖∇𝑢‖2 + 𝑐 |𝑢|2
Ω Ω Ω Ω
(︂∫︁ ∫︁ )︂
≥ min(1, 𝑐) ‖∇𝑢‖2 + |𝑢|2
Ω Ω
2
≥ min(1, 𝑐) ‖𝑢‖𝐻 1 (Ω)
Toutes les conditions sont réunies : le problème (6) admet une unique solution d’après le théorème de Lax-Milgram.

Remark 2.13
Dans la démonstration de la continuité de ℓ, n’écrivez pas ‖𝑓 ‖𝐿2 (Ω) ≤ ‖𝑓 ‖𝐻 1 (Ω) car, d’une part nous n’en avons pas
besoin, d’autre part, nous ne savons pas si 𝑓 ∈ 𝐻 1 (Ω) !

2.5.3 Conclusion

Schématiquement, nous avons :

— Si 𝑢 est solution de (5) alors 𝑢 est solution de (6)
— Le problème (6) admet une unique solution qui appartient (au moins) à 𝐻 1 (Ω)

Remark 2.14
Pourquoi travailler dans 𝐻 1 (Ω) et non dans C 1 (Ω) ? La question est légitime, d’autant que ‖·‖𝐻 1 (Ω) est une norme de
C 1 (Ω) ! Mais… C 1 (Ω) n’est pas complet pour cette norme et n’est donc pas un espace de Hilbert si on lui adjoint cette
norme : le théorème de Lax-Milgram ne pourra alors pas s’y appliquer. Il existe des normes qui complètent C 1 (Ω), mais
les hypothèses du théorème de Lax-Milgram sont elles toujours validées avec ces normes ?

2.5. Application au problème modèle 23

Maillage et Éléments Finis

24 Chapitre 2. Formulations Faibles

 CHAPITRE 3

Éléments Finis Triangulaires

And the word was death

And the word was without light
The new beatitude
Good luck, you’re on your own !

---A Perfect Circle - The Doomed

3.1 Méthode de Galerkin

3.1.1 Contexte

Dans ce chapitre nous considérons un espace de Hilbert 𝑉 muni du produit scalaire (·, ·)𝑉 et de sa norme associée ‖·‖𝑉 .
Nous considérons la formulation variationnelle suivante
Trouver 𝑢 ∈ 𝑉 tel que
{︂
(1)
∀𝑣 ∈ 𝑉, 𝑎(𝑢, 𝑣) = ℓ(𝑣).

Les formes continues 𝑎(·, ·) et ℓ(·) sont respectivement bilinéaire et linéaire, et 𝑎(·, ·) est de plus coercive. De cette
manière, le Théorème de Lax-Milgram s’applique et le problème (1) admet une unique solution.
Nous noterons (·, ·)𝑉 et ‖·‖𝑉 respectivement le produit scalaire et la norme sur 𝑉 .

25
Maillage et Éléments Finis

3.1.2 Dimension finie

Obtenir une solution de (1) est compliqué car 𝑉 est (a priori) de dimension infinie. La méthode de Galerkin consiste à
« approcher » l’espace fonctionnel 𝑉 par un espace 𝑉ℎ ⊂ 𝑉 , de dimension finie, mais toujours de Hilbert, et ce pour le
même produit scalaire ! La formulation faible (1) est alors résolue dans 𝑉ℎ uniquement, avec pour solution 𝑢ℎ :

Trouver 𝑢ℎ ∈ 𝑉ℎ tel que

{︂
(2)
∀𝑣ℎ ∈ 𝑉ℎ , 𝑎(𝑢ℎ , 𝑣ℎ ) = ℓ(𝑣ℎ ).

On espère alors que cette solution approchée 𝑢ℎ soit une bonne estimation de la solution exacte 𝑢, c’est-à-dire que

lim ‖𝑢ℎ − 𝑢‖𝑉 = 0.

ℎ→0

Remarquons tout d’abord que la formulation faible (2) admet une unique solution.

Lemma 3.1
Le problème « approché » (2) admet une unique solution.

Proof. L’espace 𝑉ℎ ⊂ 𝑉 est un sous-espace de Hilbert de 𝑉 , nous pouvons donc appliquer le Théorème de Lax-Milgram,
dont les hypothèses sur 𝑎(·, ·) et ℓ(·) sont toujours vérifiées sur 𝑉ℎ .

Travailler dans un espace de dimension finie présente un très grand avantage : on peut en extraire une base de taille finie
et ramener le calcul de 𝑢ℎ à la résolution d’un système linéaire, pour lequel les outils (numériques) ne manquent pas.
Citons par exemple les bibliothèques suivantes :
— MUMPS : solveur open-source direct parallèle
— Pardiso : solveur direct parallèle privatif d’Intel
— PETSc : Bibliothèque contenant entres autres de nombreux solveurs directs (dont MUMPS) ou itératifs (GMRES,
…)

Lemma 3.2
Soit 𝑉 un espace de Hilbert et 𝑉ℎ un sous espace de dimension finie. Soit 𝑎(·, ·) une forme bilinéaire continue et coercive
sur 𝑉 , ℓ(·) une forme linéaire continue sur 𝑉 . Le problème approché (2) admet une unique solution. De plus, cette solution
s’obtient par la résolution d’un système linéaire de matrice définie positive.

Proof. Le problème (2) admet toujours une unique solution d’après le Théorème de Lax-Milgram. Comme 𝑉ℎ est de
dimension finie, notée 𝑁ℎ , nous pouvons en extraire une base (𝜙1 , 𝜙2 , . . . , 𝜙𝑁ℎ ) et écrire
𝑁ℎ
∑︁
𝑢ℎ = 𝑢𝐼 𝜙𝐼 .
𝐼=1

La formulation faible peut alors se réécrire sur les fonctions de cette base uniquement :
𝑁ℎ
∑︁
∀𝐼, 𝑎(𝜙𝐽 , 𝜙𝐼 )𝑢𝐽 = ℓ(𝜙𝐼 ),
𝐽=1

ou encore

𝐴ℎ 𝑈ℎ = 𝐵ℎ ,

26 Chapitre 3. Éléments Finis Triangulaires

 Maillage et Éléments Finis

avec 𝐴ℎ = (𝑎(𝜙𝐽 , 𝜙𝐼 ))1≤𝐼,𝐽≤𝑁ℎ , 𝑈ℎ = (𝑢𝐼 )1≤𝐽≤𝑁ℎ et 𝐵ℎ = (ℓ(𝜙𝐼 ))1≤𝐼≤𝑁ℎ . Montrons maintenant que la matrice
𝐴ℎ est définie positive :

∀𝑊ℎ ∈ R𝑁ℎ , 𝑊ℎ = (𝑤𝐼 )1≤𝐼≤𝑁ℎ ,

𝑁ℎ ∑︁
∑︁ 𝑁ℎ
(𝑊ℎ , 𝐴ℎ 𝑊ℎ ) = 𝑊ℎ𝑇 𝐴ℎ 𝑊ℎ = 𝑤𝑗 𝑎(𝜙𝐼 , 𝜙𝐽 )𝑤𝐼
𝐼=1 𝐽=1
𝑁ℎ ∑︁
∑︁ 𝑁ℎ
= 𝑎(𝑤𝐼 𝜙𝐼 , 𝑤𝑗 𝜙𝐽 )
𝐼=1 𝐽=1
(︃ 𝑁 𝑁ℎ
)︃
∑︁ℎ ∑︁
=𝑎 𝑤𝐼 𝜙𝐼 , 𝑤𝑗 𝜙𝐽
𝐼=1 𝐽=1

L’indice 𝐽 étant muet, nous pouvons changer son intitulé :

(︃ 𝑁 𝑁ℎ
)︃
∑︁ℎ ∑︁
(𝑊ℎ , 𝐴ℎ 𝑊ℎ ) = 𝑎 𝑤𝐼 𝜙𝐼 , 𝑤𝐼 𝜙𝐼
𝐼=1 𝐼=1

Nous utilisons maintenant la coercivité de 𝑎(·, ·) :

⃦𝑁 ⃦2
⃦∑︁ℎ ⃦
(𝑊ℎ , 𝐴ℎ 𝑊ℎ ) ≥ 𝛼 ⃦ 𝑤𝐼 𝜙𝐼 ⃦ .
⃦ ⃦
⃦ ⃦
𝐼=1 𝑉
⃦∑︀ ⃦
⃦ 𝑁ℎ
Comme 𝛼 > 0, alors le terme (𝑊ℎ , 𝐴ℎ 𝑊ℎ ) est nul si et seulement si ⃦ 𝐼=1 𝑤𝐼 𝜙𝐼 ⃦ est nulle et donc si et seulement si
⃦
∑︀𝑁ℎ 𝑉
𝐼=1 𝑤𝐼 𝜙 𝐼 est la fonction nulle. Comme la famille (𝜙𝐼 ) 1≤𝐼≤𝑁 ℎ
forme une base de 𝑉ℎ , cela revient à dire que 𝑤𝐼 = 0
pour tout 𝐼 et donc que 𝑊ℎ est le vecteur nul. Nous avons donc montré que

∀𝑊ℎ ∈ R𝑁ℎ ∖ {0}, (𝑊ℎ , 𝐴ℎ 𝑊ℎ ) > 0.

Remark 3.1
Quelques remarques :
— La matrice 𝐴ℎ discrétise l’opérateur 𝑎(·, ·) au sens où elle est de taille finie.
— La coercivité d’une forme 𝑎(·, ·) est, en quelque sorte, l’équivalent de la définie positivité de sa matrice. La
coercivité s’applique au domaine « continu » (les fonctions ou opérateurs) tandis que la définie positivité est un
terme appliqué au domaine « algébrique » (les matrices (infinies ou non)).
— L’hypothèse de Lax-Milgram sur la coercivité de 𝑎(·, ·) est une hypothèse forte puisque la matrice 𝐴ℎ discrétisant
𝑎(·, ·) doit être définie positive !

Remark 3.2
La méthode des différences finies discrétise l’opérateur différentiel (Δ) tandis que les éléments finis (issue de la méthode
de Galerkin) approche l’espace fonctionnel. C’est une différence majeure !

3.1. Méthode de Galerkin 27

Maillage et Éléments Finis

3.2 Espace P1 -Lagrange

La méthode des éléments finis est basée sur la méthode de Galerkin, ou d’approximation interne. L’idée est d’approcher
l’espace fonctionnel 𝐻 1 (Ω) par un espace de dimension finie : l’espace éléments finis. Nous nous intéressons à un tel
premier espace : P1 − Lagrange ou plus simplement P1 , composés des fonctions linéaires par triangles.

3.2.1 Maillage triangulaire (ou triangulation)

Nous découpons maintenant le domaine en triangles pour obtenir un maillage triangulaire (ou triangulation) conforme de
Ω. Un maillage est conforme s’il suit les quelques règles simples suivantes.
⋃︀𝑁𝑡Une illustration est proposée sur la figure 1.
— L’union des 𝑁𝑡 triangles doit couvrir Ω sans le dépasser : Ω = 𝑝=1 𝐾𝑝 .
— L’intersection de deux triangles est soit vide, soit une arête commune complète à chacun des deux triangles, soit
un sommet de chacun des deux triangles.
— Une arête d’un triangle est soit une arête (complète) d’un autre triangle, soit une partie de Γ, auquel cas ce segment
est complètement inclus soit dans Γ𝐷 soit dans Γ𝑁 (il n’y a pas d’arête appartenant à la fois à Γ𝐷 et à Γ𝑁 ).

Fıg. 1 – Différents maillages, conforme ou non.

Une telle triangulation sera noté Tℎ = {𝐾𝑝 , 𝑝 = 1, . . . , 𝑁𝑡 }, l’indice ℎ faisant référence à la finesse de maillage, que
l’on définit par le grand diamètre des triangles :

ℎ := max (diam(𝐾)) = max (diam(𝐾𝑝 )).

𝐾∈Tℎ 𝑝=1,...,𝑁𝑡

Le diamètre d’un triangle est la distance maximale entre deux points du triangle. Nous notons de plus Sℎ et Aℎ les
ensembles respectivement des sommets et des arêtes de Tℎ . Pour un triangle arbitraire 𝐾, nous noterons [s0 , s1 , s2 ] ses
sommets ordonnés. De même, pour un triangle 𝐾𝑝 du maillage, ses sommets ordonnées seront notés [s𝑝0 , s𝑝1 , s𝑝2 ].

Remark 3.3
Il existe aujourd’hui des mailleurs automatiques open-source, un des plus connu et que nous utiliserons est GMSH (un
tutoriel est fourni par moi même). Le maillage automatique reste un métier à part entière tant la complexité est impor-
tante notamment en 3D et avec des géométries complexes, non forcément polygonales. D’autres parts, de nombreuses
questions sont encore ouvertes aujourd’hui dans ce domaine comme la construction automatique d’un maillage composé
de quadrangles. Nous n’entrerons pas dans les détails dans ce sujet, nous serons de simple « utilisateurs et utilisatrices ».

28 Chapitre 3. Éléments Finis Triangulaires

 Maillage et Éléments Finis

Fıg. 2 – Deux orientations possibles pour un triangle. Dans les maillages considérés, tous les triangles ont la même orien-
tation.

3.2.2 Fonction linéaire sur un triangle

Cas du segment (1D). Regardons tout d’abord le cas 1𝐷 d’un segment [𝛼, 𝛽] et d’une fonction 𝑝 linéaire sur ce segment :
𝑝(𝑥) = 𝑎𝑥 + 𝑏. Les coefficients 𝑎 et 𝑏 caractérisent la fonction 𝑝 et sont, de plus, définis de manière unique dès lors
que l’on connait la valeur de 𝑝 en 𝛼 et en 𝛽 (2 équations à 2 inconnues, linéairement indépendantes). Cette propriété
reste naturellement vraie pour un segment [⃗ ⃗ « plongé » en dimension 2. Un point x de ce segment est décrit par ses
𝛼, 𝛽]
coordonées curvilignes : x(𝑠) = (1 − 𝑠)⃗ 𝛼 + 𝑠𝛽,⃗ pour 𝑠 ∈ [0, 1], et un polynôme 𝑝 de degré 1 sur [⃗ ⃗ s’écrit alors
𝛼, 𝛽]
𝑝(x(𝑠)) = (1 − 𝑠)𝑝(⃗ ⃗
𝛼) + 𝑠𝑝(𝛽) pour 𝑠 ∈ [0, 1]. On voit clairement qu’un polynôme de degré 1 sur un segment est défini
de manière unique par ses valeurs aux extrémités 1 .
Cas du triangle (2D). Revenons maintenant dans un triangle 𝐾 non plat et notons P1 l’espace des polynômes réels de
degré 1 sur R2 , de dimension 3 :

P1 (R) = 𝑝 : R2 → R ⃒ ∃!𝑎, 𝑏, 𝑐 ∈ R tels que ∀(𝑥, 𝑦) ∈ R2 , 𝑝(𝑥, 𝑦) = 𝑎 + 𝑏𝑥 + 𝑐𝑦

{︀ ⃒ }︀

L’espace P1 (𝐾) des fonctions linéaires (ou des polynômes de degré 1) sur 𝐾 est lui aussi de dimension 3 (car 𝐾 n’est pas
plat) :

P1 (𝐾) = {𝑝 : 𝐾 → R | ∃!𝑎, 𝑏, 𝑐 ∈ R tels que ∀(𝑥, 𝑦) ∈ 𝐾, 𝑝(𝑥, 𝑦) = 𝑎 + 𝑏𝑥 + 𝑐𝑦}

Une fonction 𝑝 de P1 (𝐾) est définie de manière unique par ses 3 coefficients 𝑎, 𝑏, 𝑐. Inversement, ces trois coefficients
sont calculables dès lors que l’on connait la valeur de 𝑝 sur trois points non alignés, comme les 3 sommets du triangle (voir
la proposition suivant). Une fonction 𝑝 ∈ P1 (𝐾) est donc définie de manière unique soit par la connaissance de ses trois
coefficients, soit par la connaissance de sa valeur sur les trois sommets du triangle.

Proposition 3.1
Soit 𝐾 un triangle non dégénéré de R2 de sommets s1 , s2 , s3 . Alors, pour tout jeu de données 𝛼1 , 𝛼2 , 𝛼3 ∈ R, il existe
un unique polynôme de 𝑝 ∈ P1 (𝐾) tels que 𝑝(s𝑖 ) = 𝛼𝑖 pour 𝑖 = 1, 2, 3.

1. Au lycée on disait « entre deux points ne passe qu’une et une seule droite ».

3.2. Espace P1 -Lagrange 29

Maillage et Éléments Finis

Proof. En notant s𝑖 = (𝑥𝑖 , 𝑦𝑖 ) et 𝑝(𝑥, 𝑦) = 𝑎𝑥 + 𝑏𝑦 + 𝑐 avec 𝑎, 𝑏, 𝑐 ∈ R, alors le problème revient à résoudre le système
linéaire
⎧ ⎛ ⎞⎛ ⎞ ⎛ ⎞
⎨ 𝑎𝑥1 + 𝑏𝑦1 + 𝑐 = 𝛼1 𝑥1 𝑦1 1 𝑎 𝛼1
𝑎𝑥2 + 𝑏𝑦2 + 𝑐 = 𝛼2 ⇐⇒ ⎝ 𝑥2 𝑦2 1 ⎠ ⎝ 𝑏 ⎠ = ⎝ 𝛼2 ⎠
𝑎𝑥3 + 𝑏𝑦3 + 𝑐 = 𝛼3 𝑥3 𝑦3 1 𝑐 𝛼3
⎩

Le déterminant d’un tel système n’est autre que deux fois l’aire du triangle 𝐾 qui n’est pas dégénéré :
⃒ ⃒
⃒ 𝑥1 𝑦1 1 ⃒
⃒ ⃒
Δ = ⃒⃒ 𝑥2 𝑦2 1 ⃒⃒ = 2Aire(𝐾) ̸= 0
⃒ 𝑥3 𝑦3 1 ⃒

Le système est donc bien inversible et admet une unique solution (𝑎, 𝑏, 𝑐).

Remark 3.4
Soit une fonction 𝑣 ∈ P1 (𝐾), linéaire sur le triangle 𝐾. Sa restriction 𝑣|𝜎 à une arête 𝜎 de 𝐾 est elle même une fonction
linéaire sur 𝜎. Elle est donc complètement caractérisée par sa valeur aux sommets de l’arête, qui sont aussi des sommets
de 𝐾.

3.2.3 Fonctions linéaires par éléments (= par triangles)

Nous pouvons maintenant introduire l’espace fonctionnel P1 − Lagrange, souvent abrégé P1 et noté dans ce cours 𝑉ℎ ,
contient les fonctions continues sur Ω (le fermé de Ω) et linéaires sur chaque triangle :

𝑉ℎ := 𝑣ℎ ∈ C 0 (Ω) ⃒ ∀𝐾 ∈ Tℎ , 𝑣ℎ |𝐾 ∈ P1 (𝐾) .
{︀ ⃒ }︀

Caractérisons maintenant les fonctions de cet espace. Le premier résultat montre que deux fonctions de 𝑉ℎ sont égales si
et seulement si elles coïncident sur tous les sommets de la triangulation Tℎ .

Lemma 3.3
Si 𝑢ℎ , 𝑣ℎ ∈ 𝑉ℎ vérifient 𝑢ℎ (s) = 𝑣ℎ (s) pour tout sommet s de Tℎ , alors 𝑢ℎ = 𝑣ℎ sur Ω.

Proof. En se plaçant sur le triangle 𝐾 = (s1 , s2 , s3 ) de Tℎ , nous avons 𝑢ℎ (s𝑖 ) = 𝑣ℎ (s𝑖 ) pour 𝑖 = 1, 2, 3. La proposition
3.2.2 implique que 𝑢ℎ |𝐾 = 𝑣ℎ |𝐾 . Le triangle 𝐾 étant arbitraire, cette relation vaut sur tous les éléments de la triangulation.
Le même raisonnement peut être effectué sur chaque arête pour obtenir que 𝑢ℎ − 𝑣ℎ est nulle sur Ω tout entier.

Proposition 3.2
Pour tout jeu de données réelles (𝛼𝑖 )𝑖=1,...,𝑁𝑠 , il existe une unique fonction 𝑣ℎ ∈ 𝑉ℎ vérifiant 𝑣ℎ (s𝐼 ) = 𝛼𝑖 pour tout
𝑖 = 1, . . . , 𝑁𝑠 .

Proof. L’unicité est démontrée par le lemme 3.2.3, il manque donc l’existence. Prenons un triangle 𝐾𝑝 = (s𝑝1 , s𝑝2 , s𝑝3 ) de
Tℎ et le jeu de valeurs associé (𝛼1𝑝 , 𝛼2𝑝 , 𝛼3𝑝 ) ∈ R. La proposition 3.2.2 montre qu’il existe un unique polynôme 𝑝𝐾𝑝 de

30 Chapitre 3. Éléments Finis Triangulaires

 Maillage et Éléments Finis

P1 (𝐾𝑝 ) tel que 𝑝𝐾𝑝 (s𝑝𝑖 ) = 𝛼𝑖𝑝 pour 𝑖 = 1, 2, 3. Nous pouvons répéter cette opération pour tous les triangles 𝐾𝑝 et nous
introduisons 𝑢ℎ tel que

∀𝑝 = 1, . . . , 𝑁𝑡 , 𝑢ℎ |𝐾𝑝 = 𝑝𝐾𝑝 .

La fonction 𝑣ℎ est affine sur chaque triangle, il nous faut montrer que 𝑢ℎ ∈ C 0 (Ω) pour conclure sur son appartenance à
𝑉ℎ . Comme 𝑣ℎ est continue en chaque sommet s, il reste à montrer la continuité sur les arêtes.
Prenons 2 triangles 𝐾𝑝 et 𝐾𝑞 de Tℎ ayant une arête Σ en commun. Quitte à renuméroter, notons s1 = (𝑥1 , 𝑦1 ) et
s2 = (𝑥2 , 𝑦2 ) les deux sommets de l’arête Σ et notons

𝜎(𝑡) = s1 + 𝑡(s2 − s1 ) = (𝑥1 + 𝑡(𝑥2 − 𝑥1 ), 𝑦1 + 𝑡(𝑦2 − 𝑦1 ))

une paramétrisation de Σ. Si 𝑝𝐾𝑝 (𝑥, 𝑦) = 𝑎𝑥 + 𝑏𝑦 + 𝑐, nous avons alors, pour tout 𝑡 ∈ [0, 1] :

𝑝𝐾𝑝 (𝜎(𝑡)) = 𝑎(𝑥1 + 𝑡(𝑥2 − 𝑥1 )) + 𝑏(𝑦1 + 𝑡(𝑦2 − 𝑦1 )) + 𝑐

= 𝑎(𝑥1 + 𝑡(𝑥2 − 𝑥1 )) + 𝑏(𝑦1 + 𝑡(𝑦2 − 𝑦1 )) + 𝑐 + 𝑡(𝑐 − 𝑐)
= [𝑎𝑥1 + 𝑏𝑦1 + 𝑐] + 𝑡([𝑎𝑥2 + 𝑏𝑦2 + 𝑐] + [𝑎𝑥0 + 𝑏𝑦0 + 𝑐])
= 𝑝𝐾𝑝 (s1 ) + 𝑡(𝑝𝐾𝑝 (s2 ) − 𝑝𝐾𝑝 (s1 ))
= 𝑝𝐾𝑞 (s1 ) + 𝑡(𝑝𝐾𝑞 (s2 ) − 𝑝𝐾𝑞 (s1 ))
= 𝑝𝐾𝑞 (𝜎(𝑡)).

Autrement dit, les deux polynômes 𝑝𝐾𝑝 et 𝑝𝐾𝑞 sont égaux sur l’arête Σ. La fonction 𝑣ℎ est donc continue sur toutes les
arêtes de Tℎ en plus de l’être sur tous les triangles et tous les sommets : 𝑣ℎ est donc bien continue sur tout Ω.

3.2.4 Base de 𝑉ℎ : les fonctions de forme

Au vue de ce qui précède, deux fonctions de 𝑉ℎ sont identiques si et seulement si elles possèdent la même valeur sur chaque
sommet de Tℎ . En notant 𝑁𝑠 = card(Sℎ ) le nombre de sommets du maillage, introduisons la famille des fonctions de
forme (𝜙𝐼 )1≤𝐼≤𝑁𝑠 de 𝑉ℎ , qui sont nulles sur chaque sommet sauf un :

1 si 𝐼 = 𝐽
{︂
∀𝐼, 𝐽 = 1, ..., 𝑁𝑠 , 𝜙𝐼 (s𝐽 ) = 𝛿𝐼,𝐽 =
0 sinon.

Ces fonctions sont la généralisation en 2D des fonctions chapeau unidimensionnelles (elles ressemblent d’ailleurs encore
plus à un « chapeau » !).

Proposition 3.3
La famille (𝜙𝐼 )1≤𝐼≤𝑁𝑠 est une base de 𝑉ℎ , qui est alors de dimension 𝑁𝑠 , le nombre de sommets de la triangulation Tℎ .

Proof. Montrons que la famille des fonctions de forme est une base de 𝑉ℎ . Commençon par le caractère libre de cette
famille en prenant 𝑁𝑠 données (𝛼𝑖 )1≤𝑖≤𝑁𝑠 ,
𝑁𝑠
∑︁ 𝑁𝑠
∑︁
𝛼𝐼 𝜙𝐼 = 0 =⇒ ∀𝑗 = 1, . . . , 𝑁𝑠 , 𝛼𝐼 𝜙𝐼 (s𝑗 ) = 0
𝐼=1 𝐼=1
𝑁𝑠
∑︁
=⇒ ∀𝐽 = 1, . . . , 𝑁𝑠 , 𝛼𝐽 × 1 + (𝛼𝐽 × 0) = 0
𝐼=1,𝐼̸=𝐽

=⇒ ∀𝐽 = 1, . . . , 𝑁𝑠 , 𝛼𝐽 = 0

3.2. Espace P1 -Lagrange 31

Maillage et Éléments Finis

La famille de fonctions (𝜙𝐼 )1≤𝐼≤𝑁𝑠 est libre. Pour montrer qu’elle est génératrice, prenons une fonction
(︁∑︀ )︁ ⃒ 𝑢ℎ ∈ 𝑉ℎ et
3
plaçons nous sur le triangle 𝐾 = (s1 , s2 , s3 ) (quitte à renuméroter). Le polynôme 𝐼=1 𝑢ℎ (s𝐼 )𝜙𝐼 ⃒ coïncide avec
⃒
𝐾
le polynôme 𝑢ℎ |𝐾 sur les sommets du triangle 𝐾. Les deux étant de degré 1, nous avons alors l’égalité de ces polynômes
sur tout le triangle :

(︃ 3 )︃ ⃒
∑︁ ⃒
𝑢ℎ |𝐾 = 𝑢ℎ (s𝐼 )𝜙𝐼 ⃒ .
⃒
⃒
𝐼=1 𝐾

Cette relation étant valable sur un triangle arbitraire, elle est vraie sur Ω. La famille de fonctions (𝜙𝐼 )𝐼 est donc une base
de 𝑉ℎ .

Proposition 3.4 (Admis pour le moment)

L’espace 𝑉ℎ est inclus dans 𝐻 1 (Ω).

Lemma 3.4
Le support d’une fonction de forme 𝜙𝐼 est l’union des triangles ayant pour sommet s𝐼 :

supp(𝜙𝐼 ) = {𝐾 ∈ Tℎ | s𝐼 est un sommet de 𝐾} .

Autrement dit, en dehors de ces triangles, la fonction 𝜙𝐼 est nulle.

Proof. Prenons une fonction de forme 𝜙𝐼 associée au sommet s𝐼 , et un triangle 𝐾 tel que s𝐼 n’est pas un sommet de 𝐾.
Dans ce cas, 𝜙𝐼 est nulle sur les trois sommets de 𝐾, et est donc nulle sur le triangle tout entier.

Une illustration du support des fonctions de forme est donnée sur la Figure 3.

3.2.5 Conclusion

Pour une fonction 𝑢ℎ de 𝑉ℎ , retenons que :

— 𝑢ℎ est (par définition) continue et linéaire sur chaque triangle
— La dimension de 𝑉ℎ est égale au nombre de sommets 𝑁𝑠 du maillage. Plus le maillage est fin, plus la dimension
est grande.
— La famille (𝜙𝐼 )𝐼 des fonctions de forme est une base de 𝑉ℎ . Autrement dit, il existe 𝑁𝑠 uniques coefficients (𝑢𝐼 )𝐼 ,
tels que
𝑁𝑠
∑︁
𝑢ℎ = 𝑢𝐼 𝜙𝐼
𝐼=1

— Une fonction 𝑢ℎ de 𝑉ℎ est caractérisée par sa valeur aux 𝑁𝑠 sommets

— Les coefficients sont en fait la valeur de 𝑢ℎ aux sommets : 𝑢𝐼 = 𝑢ℎ (s𝐼 ) :
𝑁𝑠
∑︁
𝑢ℎ (s𝐼 )𝜙𝐼 (s𝐽 ) = 𝑢ℎ (s𝐽 )𝜙𝐽 (s𝐽 ) = 𝑢ℎ (s𝐽 ).
𝐼=1

— Le support d’une fonction de forme 𝜙𝐼 est l’union des triangles ayant pour sommets s𝐼 . Il est donc très petit par
rapport à Ω.

32 Chapitre 3. Éléments Finis Triangulaires

 Maillage et Éléments Finis

Fıg. 3 – Support des fonctions de forme pour le maillage d’un carré.

3.2. Espace P1 -Lagrange 33

Maillage et Éléments Finis

3.3 Assemblage des Matrices

Nous devons maintenant calculer effectivement les coefficients 𝐴𝐼,𝐽 de la matrice 𝐴 et 𝐵𝐼 du vecteur 𝐵. Nous nous
intéressons pour l’instant uniquement à la matrice 𝐴.

3.3.1 Algorithme « brut-force »

Prenons deux indices de sommets 𝐼 et 𝐽 et rappelons la valeur du coefficient 𝐴𝐼,𝐽 :

∫︁ ∫︁
𝐴𝐼,𝐽 = 𝑎(𝜙𝐽 , 𝜙𝐼 ) = ∇𝜙𝐽 · ∇𝜙𝐼 + 𝑐 𝜙𝐽 𝜙𝐼
Ω Ω

Chaque intégrale sur Ω peut être décomposée comme une somme sur les triangles 𝐾𝑝 :
𝑁𝑡 ∫︁
∑︁ 𝑁
∑︁𝑡 −1 ∫︁

𝐴𝐼,𝐽 = ∇𝜙𝐽 · ∇𝜙𝐼 + 𝑐 𝜙𝐽 𝜙𝐼

𝑝=1 𝐾𝑝 𝑝=0 𝐾𝑝

𝑁𝑡 ∫︁
∑︁
𝐵𝐼 = 𝑓 (𝑥)𝜙𝐼 (𝑥)d𝑥.
𝑝=1 𝐾𝑝

Soit deux sommets s𝐼 et s𝐽 n’appartenant pas un même triangle, alors supp(𝜙𝐼 ) ∩ supp(𝜙𝐽 ) = ∅. Autrement dit, 𝜙𝐼 𝜙𝐽
est toujours nul et donc le coefficient 𝐴𝐼,𝐽 est nul ! Vue autrement, si deux sommets s𝐼 et s𝐽 ne sont pas connectés par
une arête, alors 𝐴𝐼,𝐽=0 .
Les coefficients de 𝐴 sont donc majoritairement nuls car deux sommets pris au hasard (dans le million d’un maillage) ne
sont pas connectés. En moyenne de manière empirique, un nœud (ou sommet) est connecté au maximum à 6 à 8 autres
nœuds (en 2D). Une conséquence directe est que la matrice 𝐴 est creuse, c’est-à-dire qu’un nombre important de ses
coefficients sont nuls. Une stratégie de stockage creux est donc à utiliser, ce que nous verrons plus loin. Une manière
pratique est d’utiliser le format COO pour l’assemblage puis le format CSR pour l’algèbre linéaire et la résolution du
système.
Nous devons bien entendu construire cette matrice : calculer chacun de ses coefficients et les stocker. Un algorithme naïf
ou brut-force (mais naturel) pour calculer chaque coefficient est de boucler sur les sommets et et de remplir la matrice au
fur et à mesure, c’est-à-dire de remplir les coefficients les uns après les autres. Il est présenté dans l’algorithme brut-force.
Il est à noter que la boucle sur les triangles pourraient être simplifiée(︀en ne)︀ bouclant que sur les triangles ayant pour sommet
s𝐼 et s𝐽 . Cependant, cet algorithme a tout de même un coût en 𝑂 𝑁𝑠2 ce qui est trop important pour être utilisable en
pratique.

Code source 1 – Algorithme brut-force

For I = 1:N_s
For J = 1:N_s
A(I,J) = 0
For p = 1:N_t
A(I,J) += ∫_{K_p} (∇ ϕ_J·∇ ϕ_I) +∫_{K_p}(ϕ_J ϕ_I)
EndFor
EndFor
B(I) = 0
For p = 1:N_t
B[I] += ∫_{K_p} (f ϕ_I)
EndFor
EndFor

34 Chapitre 3. Éléments Finis Triangulaires

 Maillage et Éléments Finis

3.3.2 Algorithme d’assemblage

Une autre manière de procéder, que l’on appelle assemblage, se base sur une boucle sur les triangles plutôt que sur les
sommets. Le principe est de parcourir les triangles et de calculer des contributions élémentaires, qui vont s’ajouter petit
à petit dans la matrice 𝐴. Reprenons l’expression du coefficient 𝐴𝐼,𝐽 :
𝑁𝑡 ∫︁
∑︁ 𝑁
∑︁𝑡 −1 ∫︁
𝐴𝐼,𝐽 = ∇𝜙𝐽 · ∇𝜙𝐼 +𝑐 𝜙𝐽 𝜙𝐼
𝑝=1 𝐾𝑝 𝑝=0 𝐾𝑝
⏟ ⏞ ⏟ ⏞
Contrib. élémentaire Contrib. élémentaire

Introduisons 𝑎𝑝 (·, ·) la famille de forme bilinéaire suivante, pour 𝑝 = 1, . . . , 𝑁𝑡 :

∫︁ ∫︁
𝑎𝑝 (𝜙𝐽 , 𝜙𝐼 ) = ∇𝜙𝐽 (x) · ∇𝜙𝐼 (x)dx + 𝑐 𝜙𝐽 (x)𝜙𝐼 (x)dx
𝐾𝑝 𝐾𝑝

Ensuite, nous réécrivons la matrice 𝐴 sous la forme suivante

𝑠 −1
𝑁𝑠 𝑁∑︁
∑︁
𝐴= 𝑎(𝜙𝐽 , 𝜙𝐼 )e𝑇𝐼 e𝐽 ,
𝐼=1 𝑗=0

où e𝐼 est le vecteur de la base canonique de R𝑁𝑠 . Nous avons alors

𝑁𝑠 ∑︁
∑︁ 𝑁𝑠
𝐴= 𝑎(𝜙𝐽 , 𝜙𝐼 )e𝑇𝐼 e𝐽
𝐼=1 𝐽=1
𝑁𝑠 ∑︁
∑︁ 𝑁𝑠 ∑︁
𝑁𝑡
= 𝑎𝑝 (𝜙𝐽 , 𝜙𝐼 )e𝑇𝐼 e𝐽 (3)
𝐼=1 𝐽=1 𝑝=1
𝑁𝑡 ∑︁
∑︁ 𝑁𝑠 ∑︁
𝑁𝑠
= 𝑎𝑝 (𝜙𝐽 , 𝜙𝐼 )e𝑇𝐼 e𝐽
𝑝=1 𝐼=1 𝐽=1

Nous remarquons maintenant que 𝑎𝑝 (𝜙𝐽 , 𝜙𝐼 ) est nul dès lors que s𝐼 ou s𝐽 ne sont pas des sommets de 𝐾𝑝 (car 𝜙𝐼 𝜙𝐽 = 0
sur 𝐾𝑝 ). Finalement, la somme sur tous les sommets du maillage se réduit à une somme sur les 3 sommets du triangle
𝐾𝑝 considéré.
Nous comprenons que nous devons maintenant travailler localement dans chaque triangle. Pour cela, nous avons besoin
d’introduire une numérotation locale de chaque sommet une fonction L2G (Local To Global)permettant de basculer du
local vers le global une fonction telle que, pour 𝑝 = 1, . . . , 𝑁𝑡 et 𝑖 = 1, 2, 3 :

L2G(𝑝, 𝑖) = 𝐼 ⇐⇒ s𝑝𝑖 = s𝐼

Ainsi, pour un triangle 𝐾𝑝 , ses sommets sont numérotés [s𝑝1 , s𝑝2 , s𝑝3 ] en numérotation locale ou
[sL2G(𝑝,1) , sL2G(𝑝,2) , sL2G(𝑝,3) ] en numérotation globale, comme le montre la figure 4. Nous distinguerons la numé-
rotation globale par des lettres capitales (𝐼, 𝐽) et la numérotation locale par des minuscules (𝑖, 𝑗). Nous introduisons
aussi les fonctions de forme locales :

𝜙𝑝𝑖 = 𝜙L2G(𝑝,𝑖) |𝐾𝑝 .

Remark 3.5
Pour mieux comprendre la différence entre numérotation locale et globale, une application est disponible en ligne.

3.3. Assemblage des Matrices 35

Maillage et Éléments Finis

6 7 8

4 5 1

1 2 3
Fıg. 4 – Numérotation locale et globale

Utilisons ces nouvelles notations dans l’équation (3), en ramenant la somme sur les sommets à uniquement les sommets
du triangle considéré :
𝑁𝑡 ∑︁
∑︁ 3 ∑︁
3
𝐴= 𝑎𝑝 (𝜙𝑝𝑗 , 𝜙𝑝𝑖 )e𝑇L2G(𝑝,𝑖) eL2G(𝑝,𝑗)
𝑝=1 𝑖=1 𝑗=1

L’algorithme d’assemblage est alors

(︀ )︀complet ! Une version pseudo-code est présenté par l’algorithme d’assemblage. Sa
complexité est en 𝑂 (𝑁𝑡 ) ≪ 𝑂 𝑁𝑠2 . Comme le premier algorithme, il possède en plus l’avantage d’être parallélisable.

Code source 2 – Algorithme d’assemblage

A = 0
B = 0
For p = 1:N_t
For i = 1:3
I = L2G(p,i)
For j = 1:3
J = L2G(p,j)
A(I,J) += a_p(ϕ_j^p,ϕ_i^p)
EndFor
B(I) += l_p(ϕ_i^p)
EndFor
EndFor

Remark 3.6
Cet algorithme n’est pas encore utilisable, nous devons calculer la valeur de 𝑎𝑝 (𝜙𝑝𝑗 , 𝜙𝑝𝑖 ) et ℓ𝑝 (𝜙𝑝𝑖 ). De plus, il manque
encore les conditions de Dirichlet.

36 Chapitre 3. Éléments Finis Triangulaires

 Maillage et Éléments Finis

3.4 Calcul des Matrices Élémentaires

3.4.1 Matrices de Masse et de Rigidité

La matrice 𝐴 peut être décomposée en deux matrices : la masse et la rigidité :

𝐴 = 𝐷 + 𝑐𝑀,
— 𝑀 : la matrice de masse (ou de volume), de coefficient
∫︁
𝑀𝐼,𝐽 = 𝜙𝐽 𝜙𝐼 .
Ω

— 𝐷 : la matrice de rigidité, de coefficient

∫︁
𝐷𝐼,𝐽 = ∇𝜙𝐽 ∇𝜙𝐼 .
𝜔

Remark 3.7
Dans la littérature, cette matrice est souvent notée 𝐾, mais nous l’appelons 𝐷 pour éviter toute confusion avec les triangles,
nommés 𝐾 également.

Remark 3.8
La matrice de masse 𝑀 représente l’opérateur Identité dans la base des fonctions de forme (qui n’est pas orthogonale ni
normée !). Pour s’en convaincre, il faut regarder « l’équation » 𝑢 = 𝑓 (ou 𝐼𝑑.𝑢 = 𝑓 ) et appliquer la méthode des éléments
finis pour obenir la « formulation faible »
∫︁ ∫︁
∀𝑣ℎ , 𝑈 𝑣ℎ = 𝑓 𝑣ℎ ,
Ω Ω

qui aboutit au système linéaire suivant : 𝑀 𝑈 = 𝐵. L’opérateur Identité, appliqué à 𝑢, est bien discrétisé en 𝑀 .

Les contributions élémentaires, c’est à dire les quantités 𝑎𝑝 (𝜙𝑝𝑗 , 𝜙𝑝𝑖 ) et ℓ𝑝 (𝜙𝑝𝑖 ), peuvent elles aussi être décomposées en
deux parties. Pour rappel, les sommets d’un triangle 𝐾𝑝 seront notés [s𝑝0 , s𝑝1 , s𝑝2 ] et ordonnés dans le sens trigonométrique.
Nous noterons s𝑝𝑖 = (𝑥𝑝𝑖 , 𝑦𝑖𝑝 ) un sommet de 𝐾𝑝 et 𝜙𝑝𝑖 la fonction de forme locale associée. Nous notons 𝑀𝑝𝑒 et 𝐷𝑝𝑒 les
matrices de masse et de rigidité élémentaire du triangle 𝐾𝑝 , de coefficient respectif (𝑀𝑝𝑒 )𝑖,𝑗 et (𝐷𝑝𝑒 )𝑖,𝑗 donnés par
∫︁
𝑒
(𝑀𝑝 )𝑖,𝑗 = 𝜙𝑝𝑗 𝜙𝑝𝑖
𝐾𝑝
∫︁
(𝐷𝑝𝑒 )𝑖,𝑗 = ∇𝜙𝑝𝑗 · ∇𝜙𝑝𝑖 .
𝐾𝑝

3.4.2 Matrice de masse élémentaire

Nous nous focalisons sur la matrice de masse, le principe est similaire pour la matrice 𝐷 et est détaillé juste après.
Pour construire la matrice 𝑀 , nous avons vu qu’il était préférable de parcourir les triangles plutôt que les sommets, autre-
ment dit, plutôt que de calculer 𝑀𝐼,𝐽 directement, mieux vaut calculer, pour tout triangle 𝑝, la contribution élémentaire
(𝑀𝑝𝑒 )𝑖,𝑗 pour 𝑖, 𝑗 = 1, 2, 3, définie par :
∫︁
(𝑀𝑝𝑒 )𝑖,𝑗 = 𝜙𝑝𝑗 (x) 𝜙𝑝𝑖 (x)dx. (4)
𝐾𝑝

Chaque contribution élémentaire (𝑀𝑝𝑒 )𝑖,𝑗 est ensuite ajoutée à 𝑀𝐼,𝐽 , avec 𝐼 = L2G(𝑝, 𝑖) et 𝐽 = L2G(𝑝, 𝑗).

3.4. Calcul des Matrices Élémentaires 37

Maillage et Éléments Finis

Triangle de référence

Pour calculer la quantité élémentaire (4), plaçons nous tout d’abord dans un triangle « simple » 𝐾, ̂︀ appelé triangle
de référence. Celui-ci est souvent choisi comme étant le triangle rectangle de sommets ̂︀s1 = (0, 0), ̂︀s2 = (1, 0) et
̂︀s3 = (0, 1), ordonnés dans le sens trigonométrique. Pour différencier ce triangle d’un triangle du maillage, nous lui
adjoignons un repère (𝜉, 𝜂) dit repère paramétrique.

Fıg. 5 – Triangle de référence 𝐾

̂︀ et son repère paramétrique (𝜉, 𝜂).

Nous notons 𝜙 ̂︀𝑖 ∈ P1 (𝐾)̂︀ les trois fonctions de forme associées aux sommets ̂︀s𝑖 , pour 𝑖 = 1, 1, 3, définies par 𝜙
̂︀𝑖 (̂︀s𝑗 ) =
𝛿𝑖𝑗 . Ces fonctions 𝜙̂︀𝑖 étant des polynômes de degré un, nous pouvons les calculer analytiquement :
⎧
⎨ 𝜙 ̂︀1 (𝜉, 𝜂) = 1 − 𝜉 − 𝜂
𝜙
̂︀2 (𝜉, 𝜂) = 𝜉
𝜙
̂︀3 (𝜉, 𝜂) = 𝜂
⎩

Lemma 3.5
Dans le triangle 𝐾,
̂︀ la matrice de masse élémentaire 𝑀̂︁𝑒 = (𝑀̂︁𝑒 )1≤𝑖,𝑗≤3 de coefficient
𝑖,𝑗
∫︁
̂︁𝑒 =
𝑀 𝑖,𝑗 𝜙 ̂︀𝑖 d(𝜉, 𝜂),
̂︀𝑗 𝜙
𝐾
̂︀

est donnée par

⎛ ⎞
2 1 1
̂︁𝑒 1 ⎝
𝑀 = 1 2 1 ⎠.
24
1 1 2

38 Chapitre 3. Éléments Finis Triangulaires

 Maillage et Éléments Finis

Proof. Prenons tout d’abord le cas 𝑖 = 𝑗 = 2, soit 𝜙

̂︀𝑖 (𝜉, 𝜂) = 𝜙̂︀𝑗 (𝜉, 𝜂) = 𝜉. Dans ce cas :
∫︁ 1 ∫︁ 1−𝜉 ∫︁ 1 [︂ 3 ]︂1
𝜉4
∫︁
𝜉 1 1 1
𝜉 d(𝜉, 𝜂) =
2
𝜉 d𝜂d𝜉 =
2
(1 − 𝜉)𝜉 d𝜉 =
2
− = − = .
𝐾̂︀ 0 0 0 3 4 0 3 4 12
Les calculs sont similaires pour 𝑖 = 1 et 𝑖 = 3. Prenons maintenant 𝑖 ̸= 𝑗, par exemple 𝑖 = 3 et 𝑗 = 2 :
∫︁ 1 (︃∫︁ 1−𝜉 )︃
1 1
∫︁ ∫︁ [︂ ]︂
2 1 1 2 1 1
𝜉𝜂d(𝜉, 𝜂) = 𝜂d𝜂 𝜉d𝜉 = (1 − 𝜉) 𝜉d𝜉 = − + = .
𝐾̂︀ 0 0 2 0 2 2 3 4 24
Les calculs sont similaires pour les autres combinaisons.

Triangle quelconque

Changement de coordonnées. Soit un triangle 𝐾𝑝 du maillage et supposons que nous disposions d’une transformation
̂︀ en 𝐾𝑝 avec en plus 𝑇𝑝 (̂︀s𝑖 ) = s𝑝 (conservation
bijective et linéaire 𝑇𝑝 permetteant de transformer le triangle de référence 𝐾 𝑖
de l’ordre des sommets). Cette fonction 𝑇𝑝 transforme les coordonnées paramétriques (𝜉, 𝜂) en coordonnées physiques
(𝑥, 𝑦) avec (𝑥, 𝑦) = 𝑇𝑝 (𝜉, 𝜂) ∈ 𝐾𝑝 , et conserve « l’ordre des sommets ».

Fıg. 6 – Transformation entre le triangle de référence 𝐾

̂︀ et un triangle quelconque 𝐾𝑝 .

Nous avons 𝜙𝑝𝑗 (𝑥, 𝑦) = 𝜙𝑝𝑗 (𝑇𝑝 (𝜉, 𝜂)) avec 𝜙𝑝𝑗 ∘ 𝑇𝑝 ∈ P1 (𝐾)
̂︀ et 𝜙𝑝 ∘ 𝑇𝑝 (̂︀s𝑖 ) = 𝛿𝑖𝑗 , soit exactement les mêmes propriétés
𝑗
𝑝
que les 𝜙
̂︀𝑖 . Par unicité, nous avons 𝜙𝑗 ∘ 𝑇𝑝 = 𝜙̂︀𝑗 .
En notant 𝐽𝑝 la matrice Jacobienne de 𝑇𝑝 , alors la quantité (𝑀𝑝𝑒 )𝑖,𝑗 peut alors s’écrire, par changement de variables :
∫︁ ∫︁
(𝑀𝑝𝑒 )𝑖,𝑗 = 𝜙𝑝𝑗 (𝑥, 𝑦)𝜙𝑝𝑖 (𝑥, 𝑦)d(𝑥, 𝑦) = |det(𝐽𝑝 )| 𝜙
̂︀𝑗 (𝜉, 𝜂)𝜙
̂︀𝑖 (𝜉, 𝜂)d(𝜉, 𝜂)
𝐾𝑝
⏟𝐾
̂︀
⏞
Déjà calculé !

Ainsi, pour calculer la matrice élémentaire d’un triangle 𝐾𝑝 quelconque, nous n’avons besoin que du déterminant de la
Jacobienne : det(𝐽𝑝 ).
Expression et Jacobienne de la transformation. La transformation que nous cherchons, 𝑇𝑝 , est linéaire et « conserve »
les sommets et leur ordre. Pour obtenir son expression, nous construisons des fonctions d’interpolation géométrique,
(𝜓̂︀𝑖 )1≤𝑖≤3 , linéaires sur 𝐾
̂︀ et telles que :

∀𝑖, 𝑗 = 1, 2, 3 𝜓̂︀𝑖 (̂︀s𝑗 ) = 𝛿𝑖𝑗 .

La transformation aura alors pour expression :
𝑇𝑝 : 𝐾
̂︀ → 𝐾𝑝
(𝜉, 𝜂) ↦→ 𝑇𝑝 (𝜉, 𝜂) = (𝑥, 𝑦) = 𝜓̂︀1 (𝜉, 𝜂)s𝑝1 + 𝜓̂︀2 (𝜉, 𝜂)s𝑝2 + 𝜓̂︀3 (𝜉, 𝜂)s𝑝3 .

3.4. Calcul des Matrices Élémentaires 39

Maillage et Éléments Finis

En d’autres termes, les fonctions d’interpolation géométrique 𝜓̂︀𝑖 sont ici identiques aux fonctions de forme 𝜙
̂︀𝑖 :
⎧
⎨ 𝜓̂︀1 (𝜉, 𝜂) = 1 − 𝜉 − 𝜂
⎪
𝜓̂︀2 (𝜉, 𝜂) = 𝜉
⎪
⎩ 𝜓̂︀ (𝜉, 𝜂) = 𝜂
3

La matrice Jacobienne de la transformation est alors donnée par

𝜕𝑥 𝜕𝑥
⎛ ⎞
(︂ 𝑝
𝑥2 − 𝑥𝑝1 𝑥𝑝3 − 𝑥𝑝1
)︂
⎜ 𝜕𝜉 𝜕𝜂 ⎟ =
𝐽𝑝 = ⎝ 𝜕𝑦 ,
𝜕𝑦 ⎠ 𝑦2𝑝 − 𝑦1𝑝 𝑦3𝑝 − 𝑦1𝑝
𝜕𝜉 𝜕𝜂
et son déterminant vaut
|det(𝐽𝑝 )| = |(𝑥𝑝2 − 𝑥𝑝1 )(𝑦3𝑝 − 𝑦1𝑝 ) − (𝑥𝑝3 − 𝑥𝑝1 )(𝑦2𝑝 − 𝑦1𝑝 )|
= 2|𝐾𝑝 | ̸= 0,

ce qui implique que le déterminant est non nul puisque le triangle n’est pas dégénéré : la transformation 𝑇𝑝 est bien
inversible.

Remark 3.9
Quand 𝜓̂︀𝑖 = 𝜙 ̂︀𝑖 , nous parlons d’éléments finis isoparamétriques. Il convient de retenir que ce choix n’est pas obligatoire
et les fonctions 𝜓̂︀𝑖 et 𝜙
̂︀𝑖 sont indépendantes. En particulier, pour obtenir des éléments courbes, les fonctions 𝜓̂︀𝑖 pourraient
être quadratiques par exemple.

Expression finale de la matrice élémentaire.

Lemma 3.6
La matrice de masse élémentaire 𝑀𝑝𝑒 = ((𝑀𝑝𝑒 )𝑖,𝑗 )0≤𝑖,𝑗≤2 du triangle 𝐾𝑝 a pour expression
⎛ ⎞
2 1 1
|𝐾𝑝 | ⎝
𝑀𝑝𝑒 = 1 2 1 ⎠.
12
1 1 2

3.4.3 Matrice de rigidité élémentaire

Nous appliquons la même procédure pour la matrice de rigidité 𝐷, autrement dit, nous calculons les matrices de rigidité
élémentaire 𝐷𝑝𝑒 définies par
∫︁
(𝐷𝑝𝑒 )𝑖,𝑗 = ∇𝜙𝑝𝑗 (𝑥, 𝑦) · ∇𝜙𝑝𝑖 (𝑥, 𝑦)d(𝑥, 𝑦).
𝐾𝑝

40 Chapitre 3. Éléments Finis Triangulaires

 Maillage et Éléments Finis

Triangle de référence

Lemma 3.7
Dans le triangle de référence 𝐾,
̂︀ la matrice de rigidité élémentaire 𝐷 ̂︀ = (𝐷 ̂︀ 𝑖,𝑗 )1≤𝑖,𝑗≤3 de coefficient
∫︁
𝐷
̂︀ 𝑖,𝑗 = ∇𝜙̂︀𝑗 (𝜉, 𝜂) · ∇𝜙
̂︀𝑖 (𝜉, 𝜂)d(𝜉, 𝜂),
𝐾
̂︀

a pour expression
⎛ ⎞
2 −1 −1
1
̂︀ = ⎝ −1
𝐷 1 0 ⎠
2
−1 0 1

Proof. Les gradients des fonctions de forme 𝜙

̂︀𝑗 sont donnés par :
(︂ )︂ (︂ )︂ (︂ )︂
−1 1 0
∇𝜉,𝜂 𝜙̂︀0 = , ∇𝜉,𝜂 𝜙 ̂︀1 = , ∇𝜉,𝜂 𝜙
̂︀2 = .
−1 0 1
La matrice étant symétrique, nous pouvons limiter les calculs à la partie triangulaire supérieure :
∫︁ ∫︁ (︂ )︂ ∫︁
−1
𝐷1,1 =
̂︀ ∇𝜙 ̂︀1 d(𝜉, 𝜂) =
̂︀1 · ∇𝜙 (−1, −1) d(𝜉, 𝜂) = 2 d(𝜉, 𝜂) = 1
𝐾̂︀ 𝐾̂︀ −1 𝐾 ̂︀
∫︁ ∫︁ (︂ )︂ ∫︁
1 1
𝐷̂︀ 2,2 = ∇𝜙 ̂︀2 d(𝜉, 𝜂) =
̂︀2 · ∇𝜙 (1, 0) d(𝜉, 𝜂) = d(𝜉, 𝜂) = =𝐷 ̂︀ 3,3
𝐾̂︀ 𝐾̂︀ 0 𝐾̂︀ 2
∫︁ ∫︁ (︂ )︂ ∫︁
1 1
𝐷1,2 =
̂︀ ∇𝜙 ̂︀2 d(𝜉, 𝜂) =
̂︀1 · ∇𝜙 (−1, −1) d(𝜉, 𝜂) = − d(𝜉, 𝜂) =−
𝐾̂︀ 𝐾̂︀ 0 𝐾̂︀ 2
∫︁ ∫︁ (︂ )︂ ∫︁
0 1
𝐷̂︀ 1,3 = ∇𝜙 ̂︀3 d(𝜉, 𝜂) =
̂︀1 · ∇𝜙 (−1, −1) d(𝜉, 𝜂) = − d(𝜉, 𝜂) =−
𝐾̂︀ 𝐾̂︀ 1 𝐾̂︀ 2
∫︁ ∫︁ (︂ )︂
0
𝐷̂︀ 2,3 = ∇𝜙 ̂︀3 d(𝜉, 𝜂) =
̂︀2 · ∇𝜙 (1, 0) d(𝜉, 𝜂) = 0.
𝐾̂︀ 𝐾̂︀ 1

Triangle quelconque

Pour calculer les dérivées partielles selon 𝑥 et 𝑦 de 𝜙

̂︀𝑗 , nous utilisons la dérivée de fonction composée :
⎛ 𝑝⎞ ⎛
𝜕𝜙𝑗 𝜕𝜉 𝜕𝜂 ⎞ ⎛ 𝜕 𝜙 ̂︀𝑗 ⎞
⎜ 𝜕𝑥 ⎟ ⎜ 𝜕𝑥 𝜕𝑥 ⎟ ⎜ 𝜕𝜉 ⎟
⎜ 𝑝⎟ = ⎝
⎝ 𝜕𝜙𝑗 ⎠ 𝜕𝜉 𝜕𝜂 ⎠ ⎝ 𝜕 𝜙 ̂︀𝑗 ⎠
𝜕𝑦 𝜕𝑦 𝜕𝑦 𝜕𝜂
En notant 𝐵𝑝 la matrice de passage, nous avons

∇𝑥,𝑦 𝜙𝑝𝑗 (𝑥, 𝑦) = 𝐵𝑝 ∇𝜉,𝜂 𝜙

̂︀𝑗 (𝜉, 𝜂).

L’opération « inverse » nous donne :

⎛ 𝑝⎞
⎛ 𝜕𝜙
̂︀ ⎞
𝑗
⎛ 𝜕𝑥 𝜕𝑦 ⎞ 𝜕𝜙𝑗
⎜ 𝜕𝜉 ⎟ ⎜ 𝜕𝜉 𝜕𝜉 ⎟ ⎜
⎜ 𝜕𝑥 ⎟ ̂︀𝑗 (𝜉, 𝜂) = (𝐽𝑝 )𝑇 ∇𝑥,𝑦 𝜙𝑝𝑗 (𝑥, 𝑦).
⎟
̂︀ ⎠ = ⎝ 𝜕𝑥
⎝ 𝜕𝜙
𝑗 𝜕𝑦 ⎠ ⎝ 𝜕𝜙𝑝𝑗 ⎠ ⇐⇒ ∇𝜉,𝜂 𝜙
𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝑦

3.4. Calcul des Matrices Élémentaires 41

Maillage et Éléments Finis

Nous en déduisons que 𝐵𝑝 = (𝐽𝑝𝑇 )−1 , en particulier, dans le cas d’une transformation linéaire de triangle, nous obtenons :

𝑦3𝑝 − 𝑦1𝑝 𝑦1𝑝 − 𝑦2𝑝

(︂ )︂
1
𝐵𝑝 = .
det(𝐽𝑝 ) 𝑥𝑝1 − 𝑥𝑝3 𝑥𝑝2 − 𝑥𝑝1

Au final, comme 𝑋 · 𝑌 = 𝑋 𝑇 𝑌 , nous obtenons

∫︁ ∫︁
(∇𝜙𝑝𝑗 )𝑇 ∇𝜙𝑝𝑖 d(𝑥, 𝑦) = |det(𝐽𝑝 )| (∇𝜙 ̂︀𝑖 d(𝜉, 𝜂).
̂︀𝑗 )𝑇 (𝐵𝑝𝑇 𝐵𝑝 )∇𝜙 (5)
𝐾𝑝 𝐾
̂︀

En éléments finis P1 , les fonctions de forme sont linéaires et leur gradient est donc constant. Nous pouvons alors sortir les
termes ∇𝜙
̂︀𝑖 et ∇𝜙 ̂︀𝑗 de l’intégral pour obtenir le lemme suivant.

Lemma 3.8
Les coefficients a matrice de rigidité élémentaire 𝐷𝑝𝑒 = ((𝐷𝑝𝑒 )𝑖,𝑗 )1≤𝑖,𝑗≤3 sont obtenus pas la relation suivante
∫︁
𝑒
(𝐷𝑝 )𝑖,𝑗 = ∇𝜙𝑝𝑗 (𝑥, 𝑦) · ∇𝜙𝑝𝑖 (𝑥, 𝑦)d(𝑥, 𝑦),
𝐾𝑝

̂︀𝑗 )𝑇 (𝐵𝑝𝑇 𝐵𝑝 )∇𝜙

= |𝐾𝑝 | (∇𝜙 ̂︀𝑖 .

Proof. Pour les éléments

⃒ ⃒finis P , les gradients ∇𝜙
1
̂︀𝑗 sont constants et peuvent être sortis de l’intégrale. De plus, comme
|det(𝐽𝑝 )| = 2 |𝐾𝑝 | et ⃒𝐾 ⃒ = , nous avons
⃒ ̂︀ ⃒ 1
2

∫︁
∇𝜙𝑝𝑗 · ∇𝜙𝑝𝑖 dx = |𝐾𝑝 | (∇𝜙
̂︀𝑗 )𝑇 (𝐵𝑝𝑇 𝐵𝑝 )∇𝜙
̂︀𝑖 .
𝐾𝑝

3.4.4 Second membre (ou RHS ou Membre de droite)

Étudions maintenant les termes du membre de droite comme

∫︁
𝑓 (x)𝜙𝑝𝑖 (x)dx.
𝐾𝑝

Sauf pour certaines fonctions 𝑓 particulières, nous ne pourrons certainement pas calculer explicitement ce terme, nous
devons approcher cette intégrale à l’aide d’une formule de quadrature en passant à l’éléments de référence :
∫︁ ∫︁
𝑝
𝑓 (x)𝜙𝑖 (x)dx = |det(𝐽𝑝 )| 𝑓 (x(𝜉, 𝜂))𝜙
̂︀𝑖 (𝜉, 𝜂)d(𝜉, 𝜂)
𝐾𝑝 𝐾
̂︀
𝑀
∑︁−1
≃ |det(𝐽𝑝 )| 𝜔𝑚 𝑓 (x(𝜉𝑚 , 𝜂𝑚 ))𝜙(𝜉
̂︀ 𝑚 , 𝜂𝑚 ).
𝑚=0

Les points (𝜉𝑚 , 𝜂𝑚 ) sont appelés points de quadrature (parfois points de Gauss, même si la règle de quadrature utilisée
n’est pas de Gauss) et les quantités 𝜔𝑚 ∈ R les poids associés. Notons que le point x(𝜉𝑚 , 𝜂𝑚 ) s’obtient par l’expression
vue précédemment :
2
∑︁
x(𝜉𝑚 , 𝜂𝑚 ) = s𝑝𝑖 𝜓̂︀𝑖 (𝜉𝑚 , 𝜂𝑚 ).
𝑖=0

42 Chapitre 3. Éléments Finis Triangulaires

 Maillage et Éléments Finis

Nous présentons ici deux règles de quadrature pour l’intégrale 𝐾̂︀ 𝑔(x)dx sur 𝐾
̂︀ d’une fonction 𝑔 quelconque. La première
∫︀
̂︀
règle est exacte pour des polynômes de degré 1, la deuxième pour des polynômes de degré 2 (règles de Hammer) :

𝜉𝑚 𝜂𝑚 𝜔𝑚 Degré de précision
1/3 1/3 1/6 1
1/6 1/6 1/6 2
4/6 1/6 1/6
1/6 4/6 1/6

Remark 3.10
Les formules de quadrature ont évidemment un impact sur la qualité de l’approximation, toutefois, elles jouent un rôle
relativement mineur par rapport aux autres approximations (et l’on peut choisir plus de points d’intégration !).

3.5 Matrice Creuse

La méthde des éléments finis mène à la résolution d’un problème linéaire du type :

𝐴𝑥 = 𝑏,

où la matrice 𝐴 est creuse, c’est-à-dire que, majoritairement, les coefficients de 𝐴 sont nuls. Pour minimiser la mémoire
occupée par la matrice, seuls les coefficients non-nuls sont stockés. Ceci permet également d’améliorer notablement les
performances du produit matrice-vecteur, en passant d’une complexité de 𝑂(𝑁 2 ) à 𝑂(nnz) où nnz est le nombre de
coefficients non-nuls (nnz = number of non-zeros).
Il existe plusieurs formats de matrices creuses. Parmi les plus connus et utilisés : les formats COO (COOrdinates) et CSR
(Compressed Sparse Row ou CRS pour Compressed Row Storage).

3.5.1 Format COO

Principe

Relativement naturel et simple à comprendre et utiliser. La matrice est stockée sous la forme de trois tableaux row, col
et val, tous trois de taille nnz et contenant respectivement l’indice ligne, colonne et le coefficient non nuls de la matrice.
En d’autre termes, pour i = 0, …, (nnz-1),

𝐴(row[𝑖], col[𝑖]) = val[𝑖]. (6)

L’avantage de ce format est la facilité d’implémentation et la possibilité d’ajouter des coefficients « à la volée ». En effet,
les tableaux row, col et val n’ont pas besoin d’être triés selon l’ordre indices.
Prenons la matrice exemple suivante avec nnz =10 :
⎛ ⎞
3 0 0 2 1
⎜0 0 5 8 0⎟
(7)
⎜ ⎟
⎜0
𝐴=⎜ 1 2 0 0⎟
⎟.
⎝0 0 9 0 0⎠
0 0 10 4 0

Le stockage COO de cette matrice prendra alors la forme suivante :

3.5. Matrice Creuse 43

Maillage et Éléments Finis

Indice 0 1 2 3 4 5 6 7 8 9
row 0 0 0 1 1 2 2 3 4 4
col 0 3 4 2 3 1 2 2 2 3
val 3 2 1 5 8 1 2 9 10 4

Doublons

Le format COO peut autoriser les doublons, c’est-à-dire des coefficients ayant les mêmes indices ligne et colonne qu’un
autre. En reprenant l’exemple ci-dessus et en divisant le dernier coefficient en deux, nous pourrions obtenir le stockage
suivant :

Indice 0 … 8 9 10
row 0 … 4 4 4
col 0 … 2 3 3
val 3 … 10 1 3

Cette propriété est extrêmement pratique pour les éléments finis et l’algorithme d’assemblage ! En effet, chaque contribu-
tions élémentaires peut être ajouté aux tableaux row, col et val.

Du COO au Dense

Pour reconstuire la matrice sous format dense, le pseudo-code ci-dessous fonctionnerait et autorise d’avoir une dupplica-
tion de coefficients (du fait du += ) :

A = zeros(N,N)
for (i = 0; i < [Link](); i++)
A(row[i], col[i]) = val[i]
end

Produit Matrice-Vecteur

Un pseudo code serait le suivant :

// y = A*x
y = zeros(n) // vecteur nul
for (i = 0; i < [Link](); i++)
y[row[i]] += val[i] * x[col[i]]
end

44 Chapitre 3. Éléments Finis Triangulaires

 Maillage et Éléments Finis

Triplets

Plutôt que 3 tableaux, une matrice au format COO peut aussi être stockée sous forme d’un tableau de triplets (i,j,val), ce
qui donnerait pour la matrice (7) :

Indice 0 1 2 3 4 5 6 7 8 9
Triplets [0,0,3] [0,3,2] [0,4,1] [1,2,5] [1,3,8] [2,1,1] [2,2,2] [3,2,9] [4,2,10] [4,3,4]

Conclusion

Le format COO est très souple et permet de construire une matrice aisément, cependant il présente les défauts suivants :
— Deux adressages indirects sont nécessaires pour effectuer le produit matrice vecteur
— Les accès aux données ne sont pas a priori connus
— Absence de méthode rapide pour obtenir un terme de la matrice connaissant ses indices ligne et colonne
Dans la pratique, le format COO est souvent utilisée comme format « tampon » pour stocker la matrice au format CSR,
bien plus efficace pour les opérations d’algèbre linéaire. Le stockage sous forme de triplets est alors le plus pratique.

3.5.2 Format CSR

Principe

Le format CSR est spécialisé dans les opérations d’algèbres linéaires et pallie les défauts du COO. Son nom vient du fait
que le tableau row est compressé. Une matrice au format CSR est composée des deux tableaux col et val, comme pour
le COO et ordonnés par « lignes », et le tableau row est défini ainsi :
— Sa taille est fixée à n+1 (n=nombre de lignes de la matrice)
— row[i] est maintenant l’indice du premier élément non nul de la ligne i dans les tableaux col et val
Par exemple, le stockage CSR de la matrice (7) est :

Indice 0 1 2 3 4 5 6 7 8 9
row 0 3 5 7 8 10
col 0 3 4 2 3 1 2 2 2 3
val 3 2 1 5 8 1 2 9 10 4

Le tableau row est compressé par rapport au format COO puisque sa taille est maintenant de n+1, bien inférieure à
nnz ! Sur une petite matrice, le gain mémoire est très faible, mais sur une matrice à plusieurs millions d’entrée, cette
stratégie devient payante. D’autre part, l’absence de doublon de coefficients et le fait que les tableaux sont triés permettent
d’améliorer significativement les opérations d’algèbres linéaires.

Du CSR au Dense

Le pseudo code pour reconstruire la matrice dense associé ressemblerait à ceci :

A = zeros(N,N)
for (i = 0; i < [Link]() - 1; i++)
for (j = row[i]; j < row[i+1]; j++)
A(i, col[j]) = val[j]
end
end

3.5. Matrice Creuse 45

Maillage et Éléments Finis

Produit Matrice - Vecteur

Le pseudo-code est le suivant

// y = A*x
y = zeros([Link]() - 1)
for (i = 0; i < [Link]()-1; i++)
for (j = row[i]; j < row[i+1]; j++)
// Parcours des indices colonnes de la ligne i
y[i] += val[j]*x[col[j]];
end
end

Nous noterons que, cette fois-ci, les coefficients des vecteurs sont parcourus consécutivement.

Conclusion

Le format CSR est rigide : il est très coûteux d’ajouter des éléments dans la matrice. Ainsi et afin de ne pas perdre en
efficacité, il est nécessaire de connaître à l’avance l’emplacement des coefficients non nuls de la matrice avant de la
construire. En revanche, une fois construite, cette forme de stockage est très efficace.

3.5.3 Du COO au CSR

Principe

La souplesse du format COO permet de construire la matrice en ajoutant les triplets des coefficients (i,j,val) au fur et à
mesure. Ensuite, une fois tous les triplets sauvegardés, ils sont triés (ou assemblés) et les doublons fusionnés. Il ne reste
alors plus qu’à extraire les tableaux row, col et val du tableau de triplets et à compresser le vecteur row pour obtenir
une matrice CSR.

Utilisation

En supposant les fonctions existantes, le pseudo-code suivant permet de passer d’une matrice A au format COO à une
matrice B au format CSR :

MatriceCOO A(n) // COO

MatriceCSR B(n) // CSR
// Ajout des triplets
[Link](0,0,2.);
[Link](0,1,-1.1);
[...]
// Convertisseur en CSR
B = A.to_csr();

46 Chapitre 3. Éléments Finis Triangulaires

 CHAPITRE 4

Conditions aux bords

Thug : Where is the money Lebowski ?

The Dude : It’s uh… uh… it’s down there somewhere, let me take another look.

---The Big Lebowski (movie)

4.1 Conditions de Neumann hétérogène

4.1.1 Théorie

Rajoutons maintenant la condition de Neumann hétérogène à notre problème (𝑔𝑁 ̸= 0) :

{︂
−Δ𝑢 + 𝑢 = 𝑓 (Ω),
𝜕n 𝑢 = 𝑔𝑁 (Γ := 𝜕Ω).

Après multiplication par des fonctions test et intégration par partie, nous obtenons la formulation variationnelle
∫︁ ∫︁ ∫︁ ∫︁
∇𝑢 · ∇𝑣 + 𝑢𝑣 − (𝜕n 𝑢)𝑣 = 𝑓 𝑣.
Ω Ω 𝜕Ω=Γ Ω

En utilisant la condition 𝜕n 𝑢 = 𝑔𝑁 sur Γ, nous obtenons la formulation variationnelle suivante :

Trouver 𝑢 ∈ 𝐻 1 (Ω) tel que

{︂
(1)
∀𝑣 ∈ 𝐻 1 (Ω), 𝑎(𝑢, 𝑣) = ℓ(𝑣),
avec
∫︁ ∫︁
𝑎(𝑢, 𝑣) := ∇𝑢 · ∇𝑣 + 𝑢𝑣
Ω Ω
∫︁ ∫︁
ℓ(𝑣) := 𝑔𝑁 𝑣 + 𝑓 𝑣.
Γ Ω

Pour pouvoir appliquer le Théorème de Lax-Milgram, nous savons par le cas de Neumann homogène que l’application
𝑎(·, ·) est continue et coercive. Rien de neuf sous le soleil me direz-vous ? Oui mais non :

47
Maillage et Éléments Finis

— Nous ne savons pas si ℓ(·) est continue

— Pire encore, est-ce que le terme Γ 𝑔(x)𝑣(x)d𝑠(x) a un sens quand 𝑣 est dans 𝐻 1 (Ω) ?
∫︀

Nous n’avons en effet pas (encore) donné de sens à la trace (= la « restriction ») sur Γ d’une fonction de 𝐻 1 (Ω), c’est-à-dire
à 𝑣|Γ . C’est l’objet du théorème ci-dessous (admis).

Theorem 4.1 (Continuité de la Trace (admis))

Soit Γ ⊂ 𝜕Ω une partie du bord de mesure non nulle au sens de la mesure de surface. Alors il existe une unique application
𝛾Γ : 𝐻 1 (Ω) → 𝐿2 (Γ) qui est continue au sens de ‖·‖𝐻 1 (Ω) :

∃𝐶 > 0 tel que ∀𝑣 ∈ 𝐻 1 (Ω), ‖𝛾Γ (𝑣)‖𝐿2 (𝜕Ω) ≤ 𝐶 ‖𝑣‖𝐻 1 (Ω) .

Cette application est de plus caractérisée par

∀𝜙 ∈ C 1 (Ω), 𝛾Γ (𝜙) = 𝜙|Γ .

Ce théorème nous permet de montrer que la forme ℓ a un sens (chaque quantité existe) et est bien continue puisque, pour
tout 𝑣 de 𝐻 1 (Ω) :
⃒∫︁ ⃒ ⃒∫︁ ⃒
Inégalité Triangulaire
⃒ ⃒ ⃒ ⃒
|ℓ(𝑣)| ≤ ⃒⃒ 𝑓 (x)𝑣(x)dx⃒⃒ + ⃒⃒ 𝑔𝑁 (x)𝑣(x)dx⃒⃒
⃒∫︁Ω ⃒ ⃒∫︁𝜕Ω ⃒
Réécriture
⃒ ⃒ ⃒ ⃒
≤ ⃒⃒ 𝑓 (x)𝑣(x)dx⃒⃒ + ⃒⃒ 𝑔𝑁 (x)𝛾𝜕Ω (𝑣(x))⃒⃒
Ω 𝜕Ω
≤ ‖𝑓 ‖𝐿2 (Ω) ‖𝑣‖𝐿2 (Ω) + ‖𝑔‖𝐿2 (𝜕Ω) ‖𝛾𝜕Ω (𝑣)‖𝐿2 (𝜕Ω) Cauchy-Schwarz
(︁ )︁
≤ ‖𝑓 ‖𝐿2 (Ω) + 𝐶 ‖𝑔‖𝐿2 (𝜕Ω) ‖𝑣‖𝐻 1 (Ω) Cont. Trace.

4.1.2 Implémentation dans le cas P1

Nous discrétisons la formulation faible (1) en remplaçant formellement 𝐻 1 (Ω) par 𝑉ℎ pour obtenir

Trouver 𝑢ℎ ∈ 𝑉ℎ tel que

{︂
(2)
∀𝑣ℎ ∈ 𝐻 1 (Ω), 𝑎(𝑢ℎ , 𝑣ℎ ) = ℓ(𝑣ℎ ),

Nous appliquons la méthode vue précédemment pour obtenir un système linéaire équivalent à (2) :

𝐴𝑈 = 𝐵.

Les coefficients de la matrice 𝐴 et du vecteur 𝐵 sont donnés par

∫︁ ∫︁
𝐴 = (𝐴𝐼,𝐽 )𝐼,𝐽 , 𝐴𝐼,𝐽 = 𝑎(𝜙𝐽 , 𝜙𝐽 ) = ∇𝜙𝐽 · ∇𝜙𝐼 + 𝜙𝐽 𝜙𝐼
Ω Ω
∫︁ ∫︁
𝐵 = (𝐵𝐼 )𝐼 , 𝐵𝐼 = ℓ(𝜙𝐼 ) = 𝑓 𝜙𝐼 + 𝑔𝑁 𝜙𝐼
Ω
⏟ Γ⏞
Nouveau !

, Au final, seule le membre de droite est modifié par rapport au∫︀ cas de la condition de Neumann homogène. Autrement
dit, la matrice 𝐴 est identique et il nous suffit de savoir calculer Γ 𝑔𝑁 𝜙𝐼 pour obtenir le second membre
∫︀ : pour cela, nous
utilisons une règle de quadrature sur des segments. La technique est la même que pour le calcul de Ω 𝑓 𝜙𝐼 .
Décomposons Γ en 𝑁Γ segments (qui sont aussi des arêtes de triangles !) 𝜎𝑞 , 𝑞 = 1, . . . , 𝑁Γ . Chaque segment a deux
𝜎 𝜎 𝜎
sommets indicés [s1 𝑞 , s2 𝑞 ]. Nous notons 𝜙𝑖 𝑞 = 𝜙𝐼 |𝜎𝑞 la restriction de la fonction de forme 𝜙𝐼 au segment 𝑞, tel que

48 Chapitre 4. Conditions aux bords

 Maillage et Éléments Finis

𝜎
s𝑖 𝑞 = s𝐼 . Nous pouvons alors décomposer 𝐵 comme une somme de contributions élémentaires sur les triangles et les
segments.
𝑁𝑡 ∑︁
3 ∫︁ 𝑁Γ ∑︁
2 ∫︁
𝜎
∑︁ ∑︁
𝐵= 𝑓 𝜙𝑝𝑖 + 𝑔𝑁 𝜙𝑖 𝑞 .
𝑝=1 𝑖=1 𝐾𝑝 𝑞=1 𝑖=1 𝜎𝑞

Nous savons comment approcher la quantité de gauche avec une formule de quadrature adaptée. Le terme de droite
s’approche lui aussi avec une formule de quadrature 1D, par exemple la formule « 1/3 Simpson », qui est exacte pour des
𝜎 𝜎
⃦ 𝜎 𝜎 ⃦ s 𝑞 +s 𝑞
polynômes de degré 2. Nous notons |𝜎| = ⃦s1 𝑞 − s2 𝑞 ⃦ la taille du segment et s12 = 1 2 2 le milieu du segment, alors
la formule est la suivante :
|𝜎| (︀ 𝜎𝑞
∫︁
𝜎 𝜎 )︀
𝑔≈ 𝑔(s1 ) + 4𝑔(s12𝑞 ) + 𝑔(s2 𝑞 )
𝜎 6

4.2 Condition de Dirichlet

4.2.1 Homogène

Formulation Faible

Soit le problème suivant (notez l’absence du terme en 𝑢)

{︂
−Δ𝑢 = 𝑓 (Ω),
(3)
𝑢 = 0 (Γ := 𝜕Ω).

Multiplions l’EDP (3) par des fonctions tests 𝑣, intégrons sur Ω et appliquons le Théorème de Green :

−Δ𝑢 = 𝑓 =⇒ ∀𝑣, (−Δ𝑢)𝑣 = 𝑓 𝑣

∫︁ ∫︁
=⇒ ∀𝑣, − (Δ𝑢)𝑣 = 𝑓𝑣
∫︁ Ω ∫︁ Ω ∫︁
=⇒ ∀𝑣, ∇𝑢 · ∇𝑣 − (𝜕n 𝑢)𝑣 = 𝑓𝑣
Ω Γ Ω

Nous sommes théoriquement bloqué. Nous sommes tentés de dire que 𝜕n 𝑢 = 0 mais non seulement nous ne le savons pas,
mais en plus c’est très probablement faux ! Utiliser les conditions aux bords est en revanche la bonne idée. Nous savons
que 𝑢 est nul sur le bord Γ. Autrement dit, nous ne cherchons pas la valeur de la solution sur ce bord, nous la connaissons
déjà. Afin de conserver la symétrie entre 𝑢 et 𝑣, imposons à 𝑣 d’être aussi nul sur le bord et regardons ce que l’on obtient :
∫︁ ∫︁
−Δ𝑢 = 𝑓 =⇒ ∀𝑣, 𝑣|Γ = 0 ∇𝑢 · ∇𝑣 = 𝑓𝑣
Ω Ω

En terme de dérivabilité, l’espace 𝐻 1 (Ω) est suffisant pour la solution 𝑢 et les fonctions tests 𝑣. Il manque toutefois la
condition de Dirichlet, dite essentielle, qui doit être incluse dans l’espace fonctionnel. Pour cela, en rappelant que 𝛾 est
l’application trace de 𝐻 1 (Ω) sur 𝐿2 (Γ), nous introduisons un espace de Sobolev qui prend en compte cette condition de
Dirichlet.

𝐻01 (Ω) = 𝑢 ∈ 𝐻 1 (Ω) ⃒ 𝛾𝑢 = 0 ,

{︀ ⃒ }︀

La formulation faible s’écrit alors

{︃
Trouver 𝑢 ∈ 𝐻01 (Ω) tel que
∀𝑣 ∈ 𝐻01 (Ω), 𝑎(𝑢, 𝑣) = ℓ(𝑣),

4.2. Condition de Dirichlet 49

Maillage et Éléments Finis

avec
𝑎 : 𝐻01 (Ω) × 𝐻01 (Ω) −→ R∫︁
(𝑢, 𝑣) ↦−→ ∇𝑢 · ∇𝑣
Ω
ℓ: 𝐻01 (Ω) −→ R∫︁
𝑣 ↦−→ 𝑓𝑣
Ω

Remark 4.1
Attention, c’est parce que 𝑣 est nul sur Γ que l’intégrale sur Γ s’annule. Ce n’est pas parce que 𝜕n 𝑢 = 0 ! D’ailleurs, sauf
si 𝑢 = 0 partout, il y a fort à parier que 𝜕n 𝑢 ̸= 0 !

Démontrons maintenant que cette formulation faible admet une unique solution. Commençons tout d’abord par montrer
que 𝐻01 (Ω) est un espace de Hilbert.

Lemma 4.1
L’espace 𝐻 1 (Ω) est de Hilbert

Ensuite, la continuité de ℓ a déjà été démontrée dans 𝐻 1 (Ω) et donc dans 𝐻01 (Ω). Occupons nous de 𝑎(·, ·).
— Continuité de 𝑎(·, ·).
⃒∫︁ ⃒
∀𝑢, 𝑣 ∈ 𝐻01 (Ω), |𝑎(𝑢, 𝑣)| = ⃒⃒ ∇𝑢 · ∇𝑣 ⃒⃒
⃒ ⃒
Ω
≤ ‖∇𝑢‖𝐿2 (Ω) ‖∇𝑣‖𝐿2 (Ω) Cauchy Schwarz
≤ ‖𝑢‖𝐻 1 (Ω) ‖𝑣‖𝐻 1 (Ω) Inégalité des normes

— Coercivité de 𝑎(·, ·) :
∫︁
2
∀𝑢 ∈ 𝐻01 (Ω), 𝑎(𝑢, 𝑢) = ∇𝑢 · ∇𝑢 = ‖∇𝑢‖𝐿2 (Ω) ≥ . . .
Ω

La coercivité est en réalité compliquée à obtenir puisque nous aimerions avoir :

(︁ )︁
2 2 2 2
‖∇𝑢‖𝐿2 (Ω) ≥ 𝐶 ‖𝑢‖𝐻 1 (Ω) = 𝐶 ‖𝑢‖𝐿2 (Ω) + ‖∇𝑢‖𝐿2 (Ω)

L’inégalité de Poincaré vient alors à notre rescousse !

Proposition 4.1 (Inégalité de Poincaré (admise))

Il existe une constante 𝐶 ne dépendant que de Ω telle que

∀𝑢 ∈ 𝐻01 (Ω), ‖∇𝑢‖𝐿2 (Ω) ≥ 𝐶 ‖𝑢‖𝐻 1 (Ω)

Remark 4.2
L’inégalité de Poincaré est également valable si{︀la condition ⃒de Dirichlet }︀n’est posée que sur une partie Γ𝐷 du bord Γ.
Dans ce cas, l’espace considéré est 𝐻Γ1𝐷 (Ω) := 𝑣 ∈ 𝐻 1 (Ω) ⃒ 𝛾Γ𝐷 𝑣 = 0 où 𝛾Γ𝐷 : 𝐻 1 (Ω) → 𝐿2 (Γ𝐷 ) est l’application
trace sur Γ𝐷 . À noter que 𝐻Γ1𝐷 (Ω) est un Hilbert pour les mêmes raisons 𝐻01 (Ω) l’est.

50 Chapitre 4. Conditions aux bords

 Maillage et Éléments Finis

Remark 4.3
L’inégalité de Poincaré montre que la semi-norme 𝑣 ↦→ ‖∇𝑣‖𝐿2 (Ω) est une norme sur 𝐻 1 (Ω) et est équivalente à la
norme usuelle ‖·‖𝐻 1 (Ω) , puisque l’on a ‖∇𝑣‖𝐿2 (Ω) ≥ 𝐶 ‖𝑣‖𝐻 1 (Ω) ≥ 𝐶 ‖∇𝑣‖𝐿2 (Ω) .

L’inégalité de Poincaré implique la coercivité de 𝑎(·, ·). Toutes les hypothèses du théorème de Lax-Milgram sont vérifiées
et la formulation faible du problème de Dirichlet homogène admet bien une unique solution.

Implémentation

Si 𝑉ℎ est l’espace des éléments finis P1 sur Ω, alors une discrétisation naturelle de 𝐻01 (Ω) est l’espace 𝑉ℎ,0 défini par

𝑉ℎ,0 = {𝑢 ∈ 𝑉ℎ | 𝑢|Γ = 0}

Nous pouvons aussi raisonner sur le système linéaire directement. Nous séparons les degrés de liberté en deux
sous-ensembles :
1. Ceux qui appartiennent à Ω ou à Γ𝑁 : nous les noterons avec un indice 𝐼 (pour Intérieur) : 𝑢𝐼
2. Ceux qui appartiennent à Γ𝐷 , ils seront notés avec un indice 𝐷 : 𝑢𝐷
Quitte à renuméroter, le vecteur 𝑈 de degrés de liberté se réécrit
(︂ )︂
𝑢𝐼
𝑈= ,
𝑢𝐷
et le système linéaire 𝐴𝑈 = 𝐵 devient :
(︂ )︂ (︂ )︂ (︂ )︂
𝐴𝐼,𝐼 𝐴𝐼,𝐷 𝑢𝐼 𝐵𝐼
𝐴𝑈 = 𝐵 ⇐⇒ =
𝐴𝐷,𝐼 𝐴𝐷,𝐷 𝑢𝐷 𝐵𝐷
Les degrés de liberté 𝑢𝐷 sont en réalité fixés à 0 du fait de la condition de Dirichlet, autrement dit, le système à résoudre
se résume à (𝐼𝐷,𝐷 étant la matrice identité) :
(︂ )︂ (︂ )︂ (︂ )︂
𝐴𝐼,𝐼 𝐴𝐼,𝐷 𝑢𝐼 𝐵𝐼
𝐴𝑈 = 𝐵 ⇐⇒ = (4)
0 𝐼𝐷,𝐷 𝑢𝐷 0
Informatiquement, nous devons donc rendre les lignes et colonnes associées aux degrés de liberté de Dirichlet, nulles,
sauf sur la diagonale avec la valeur 1. Cette opération peut être effectuée après l’assemblage de la matrice ou lors de
l’algorithme directement.

Remark 4.4
La valeur de 1 sur la diagonale est finalement arbitraire : nous pouvons choisir n’importe quelle valeur. Pour des raisons
de précision numérique, il peut être plus pertinent de choisir comme valeur la moyenne de la somme de la diagonale de
𝐴𝐼,𝐼 (sa trace). Cette technique peu coûteuse permet d’éviter de polluer le conditionnement de la matrice par des valeurs
potentiellement trop grande ou trop petite par rapport à la « moyenne ».

Remark 4.5
Dans le cas de condition de Dirichlet homogène, ce système ce simplifie :
(︂ )︂ (︂ )︂ (︂ )︂
𝐴𝐼,𝐼 0 𝑢𝐼 𝐵𝐼
𝐴𝑈 = 𝐵 ⇐⇒ = ,
0 𝐼𝐷,𝐷 𝑢𝐷 0
ou encore, plus simplement : 𝐴𝐼,𝐼 𝑈𝐼 = 𝐵𝐼 . Le système obtenu est de plus petite taille : c’est logique, l’espace 𝑉ℎ,0 est
de dimension le nombre de sommets du maillage moins le nombre de sommets sur le bord Γ.

4.2. Condition de Dirichlet 51

Maillage et Éléments Finis

4.2.2 Condition hétérogène

Notion de relèvement

Nous considérons maintenant le cas d’une condition de Dirichlet non homogène, autrement dit, si 𝑔 ̸= 0 :
{︂
−Δ𝑢 = 𝑓 (Ω)
(5)
𝑢 = 𝑔 (Γ)
Nous pouvons introduire l’ensemble suivant

𝐻𝑔1𝐷 (Ω) = 𝑢 ∈ 𝐻 1 (Ω) ⃒ 𝛾𝑢 = 𝑔 ,

{︀ ⃒ }︀

mais ce n’est pas un espace vectoriel ! Pour remédier à ce problème, nous nous ramenons au cas d’une condition de
Dirichlet homogène en introduisant un relèvement (= une « extension », l’inverse d’une « restriction ») 𝑢𝑔 de 𝑔 : une
fonction de 𝐻 1 (Ω) telle que 𝛾𝑢𝑔 = 𝑔. Nous ne nous préoccuperons pas de savoir si une telle fonction existe et supposons
que tel est le cas. Le problème devient alors de chercher 𝑢𝑡 = 𝑢 − 𝑢𝑔 satisfaisant :
{︂
−Δ𝑢𝑡 = 𝑓 + Δ𝑢𝑔 (Ω)
(6)
𝑢𝑡 = 0 (Γ)
Nous avons vu plus haut que ce problème admet une unique solution, ce qui implique que (5) admet également une unique
solution.

Remark 4.6
Le relèvement n’est pas unique, puisque si 𝑢0 ∈ 𝐻01 (Ω) alors 𝑢𝑔 + 𝑢0 est aussi un relèvement acceptable.

Remark 4.7
Pour que le relèvement existe, il suffit que 𝑔 ∈ 𝐻 1/2 (Γ). Cet espace est composé des traces sur Γ des fonctions de 𝐻 1 (Ω) :

𝐻 1/2 (Γ) = 𝛾𝑣 ⃒ 𝑣 ∈ 𝐻 1 (Ω)

{︀ ⃒ }︀

Il contient naturellement 𝐿2 (Γ) puisque 𝛾𝑣 ∈ 𝐿2 (Γ).

Relèvement en P1

En éléments finis P1 , un relèvement naturel est la fonction 𝑢ℎ,𝑔 de 𝑉ℎ telle que

𝑔(s𝑗 ) si s𝑗 ∈ Γ𝐷 ,
{︂
𝑢ℎ,𝑔 (s𝑗 ) =
0 sinon.
Cette fonction n’est pas un relèvement de 𝑔 puisqu’elle ne coïncide avec 𝑔 que sur les sommets, mais pas nécessairement
entre ceux-ci. Toutefois, au niveau discret, elle remplit ce rôle : c’est un relèvement de l’interpolée de 𝑔 dans 𝑉ℎ (voir
la section suivante). Nous notons 𝑔ℎ le vecteur de même taille que 𝐵𝐷 et de coefficient 𝑔(s𝐼 ) avec s𝐼 ∈ Γ. En pratique,
appliquer la condition de Dirichlet hétérogène se traduit par la décomposition de la matrice ainsi :
(︂ )︂ (︂ )︂ (︂ )︂
𝐴𝐼,𝐼 𝐴𝐼,𝐷 𝑢𝐼 𝐵𝐼
= .
0 𝐼𝐷,𝐷 𝑢𝐷 𝑔ℎ
La quantité 𝑔ℎ est ici à voir comme un vecteur de coefficient 𝑔(s). La matrice obtenue est non symétrique, ce qui peut
poser des problèmes (par ex. augmentation du coût de stockage mémoire). Une astuce simple consiste à réécrire sous la
forme suivante :
(︂ )︂ (︂ )︂ (︂ )︂
𝐴𝐼,𝐼 0 𝑢𝐼 𝐵𝐼 − 𝐴𝐼,𝐷 𝑔ℎ
= .
0 𝐼𝐷,𝐷 𝑢𝐷 𝑔ℎ

52 Chapitre 4. Conditions aux bords

 Maillage et Éléments Finis

Remark 4.8
Comme pour Dirichlet homogène, nous pouvons aussi nous contenter de résoudre un système plus petit : 𝐴𝐼,𝐼 𝑢𝐼 =
𝐵𝐼 − 𝐴𝐼,𝐷 𝑔ℎ .

Remark 4.9
Le terme 𝐴𝐼,𝐷 𝑔ℎ est la version discrète du terme Δ𝑢𝑔 qui∫︀ apparait dans (6). En effet, la matrice 𝐴 discrétise l’opérateur
𝑎(·, ·) qui, ici, représente le laplacien sous sa forme faible Ω ∇𝑢 · ∇𝑣. Gardez à l’esprit que 𝐴𝐼,𝐷 n’est pas carré et prend
en argument un vecteur de la taille le nombre de sommets de Γ pour retourner un vecteur de taille le nombre de sommets
du maillage.

4.3 Condition de Fourier

4.3.1 Problème

Étudions le problème suivant pour 𝑓 et 𝑔 suffisamment régulières :

{︂
−Δ𝑢 + 𝑢 = 𝑓 (Ω),
(7)
𝜕n 𝑢 + 𝑢 = 𝑔 (Γ := 𝜕Ω).

La condition de Fourier (ou Robin ou Fourier-Robin) s’écrit aussi 𝜕n 𝑢 = 𝑔 − 𝑢 sur Γ. Après calcul, la formulation
variationnelle s’écrit
Trouver 𝑢 ∈ 𝐻 1 (Ω) tel que
{︂

∀𝑣 ∈ 𝐻 1 (Ω), 𝑎(𝑢, 𝑣) = ℓ(𝑣),

avec (𝛾 est l’application trace sur Γ) :

𝑎 : 𝐻 1 (Ω) × 𝐻 1 (Ω) → 𝐻 1
∫︁ (Ω) ∫︁ ∫︁
(𝑢, 𝑣) ↦→ ∇𝑢 · ∇𝑣 + 𝑢𝑣 + 𝛾(𝑢)𝛾(𝑣)
Ω Ω Γ
ℓ: 𝐻 1 (Ω) → 𝐻 1
∫︁ (Ω) ∫︁
𝑣 ↦→ 𝑓 𝑣 + 𝑔𝛾(𝑣)
Ω Γ

4.3.2 Existence et unicité

Nous avons vu que dans la section 4.1 consacrée à Neumann hétérogème que les intégrales sur le bord ont un sens du
fait de l’existence de l’opérateur trace. Pour une condition de Neumann, l’opérateur ℓ est le même et nous avons déjà vu
dans la section 4.1 qu’il vérifie les hypothèses du Théorème de Lax-Milgram. Il ne nous reste qu’à vérifier que 𝑎(·, ·) est
bilinéaire (trivial), continue et coercive.
— Continuité de 𝑎(·, ·) pour tout 𝑢, 𝑣 ∈ 𝐻 1 (Ω) :

4.3. Condition de Fourier 53

Maillage et Éléments Finis

⃒∫︁ ∫︁ ∫︁ ⃒
⃒ ⃒
|𝑎(𝑢, 𝑣)| = ⃒⃒ ∇𝑢 · ∇𝑣 + 𝑢𝑣 + 𝛾(𝑢)𝛾(𝑣)⃒⃒
⃒∫︁Ω ∫︁
Ω
⃒ ⃒Γ∫︁ ⃒
Inégalité Triang.
⃒ ⃒ ⃒ ⃒
≤ ⃒⃒ ∇𝑢 · ∇𝑣 + 𝑢𝑣 ⃒⃒ + ⃒⃒ 𝛾(𝑢)𝛾(𝑣)⃒⃒
Ω Ω Γ
≤ ‖𝑢‖𝐻 1 (Ω) ‖𝑣‖𝐻 1 (Ω) + ‖𝛾(𝑢)‖𝐿2 (Γ) ‖𝛾(𝑣)‖𝐿2 (Γ) Cauchy-Schwarz
2
≤ ‖𝑢‖𝐻 1 (Ω) ‖𝑣‖𝐻 1 (Ω) + 𝐶 ‖𝑢‖𝐻 1 (Ω) ‖𝑣‖𝐻 1 (Ω) Cont. Trace
2
≤ (1 + 𝐶 ) ‖𝑢‖𝐻 1 (Ω) ‖𝑣‖𝐻 1 (Ω)
La constante 𝐶 est la constante de continuité de l’opérateur Trace sur Γ définie dans le théorème ??.
— Coercivité de 𝑎(·, ·), avec 𝑢 ∈ 𝐻 1 (Ω) ∫︁: ∫︁ ∫︁
𝑎(𝑢, 𝑢) = ∇𝑢 · ∇𝑢 + 𝑢𝑢 + 𝛾(𝑢)𝛾(𝑢)
Ω
∫︁ Ω Γ
2 2
= ‖𝑢‖𝐻 1 (Ω) + |𝛾(𝑢)|
⏟Γ ⏞
≥0
2
≥ ‖𝑢‖𝐻 1 (Ω)
Le problème admet donc une unique solution.

4.3.3 Matrice de masse sur le bord

Après discrétisation dans la base éléments finis, nous sommes ramenés à la résolution du système linéaire

𝐴𝑈 = 𝑏,

où la matrice 𝐴 et le vecteur 𝑏 sont donnés par

∫︁ ∫︁ ∫︁
𝐴(𝐼, 𝐽) = ∇𝜙𝐽 · ∇𝜙𝐼 + 𝜙𝐽 𝜙𝐼 + 𝜙𝐽 |Γ 𝜙𝐼 |Γ
∫︁Ω ∫︁ Ω Γ

𝐵(𝐼) = 𝑓 𝜙𝐼 + 𝑔𝜙𝐼 |Γ
Ω Γ

Le vecteur 𝐵 se calcule grâce aux formules de quadratures vues dans les paragraphes 3.4.4 pour les triangles et dans 4.1.2
pour les segments. La matrice 𝐴 est obtenue par la somme de la matrice de rigidité 𝐷, de masse 𝑀 et d’une dernière
matrice 𝑀Γ de coefficients :
∫︁
𝑀Γ (𝐼, 𝐽) = 𝜙𝐽 |Γ 𝜙𝐼 |Γ .
Γ

Cette matrice correspond à une matrice de masse sur le bord Γ. Nous pouvons tout d’abord remarquer que 𝜙𝐼 |Γ = 0 dès
que s𝐼 n’est pas sur Γ. Comme toujours, nous préférons la décomposer en contributions élémentaires où, ici, un élément
sera un segment :
∑︁ ∫︁
𝑀Γ (𝐼, 𝐽) = 𝜙𝐽 |𝜎 𝜙𝐼 |𝜎 .
𝜎∈Γ 𝜎

Nous pouvons maintenant remarquer que la somme sur les arêtes n’en est pas une puisque l’intégrale sur 𝜎 est nulle dès
que s𝐼 ou s𝐽 n’est pas un sommet de l’arête. Cependant, n’oublions pas que nous ne calculons pas les coefficients un à un
mais que nous assemblons la matrice, autrement dit, nous parcourons chaque segment, calculons toutes les contributions
élémentaires associées à ce dernier, et additionnons le tout dans la grande matrice du système.
Autrement dit et quitte à renuméroter, nous considérons une arête 𝜎 = [s𝜎1 , s𝜎2 ], nous cherchons à calculer :
∫︁
𝑀𝜎𝑒 (𝑖, 𝑗) = 𝜙𝜎𝑗 𝜙𝜎𝑖 ,
𝜎

avec 𝜙𝜎𝑖 = 𝜙𝐼 |𝜎 et s𝜎𝑖 = s𝐼 . La matrice 𝑀𝜎𝑒 est de dimension 2x2.

54 Chapitre 4. Conditions aux bords

 Maillage et Éléments Finis

4.3.4 Calcul de la matrice

Nous introduisons la coordonnée curviligne 𝑡

‖x − s𝜎1 ‖
∀x ∈ 𝜎, 𝑡(x) = ∈ [0, 1].
‖s𝜎1 − s𝜎2 ‖

Quand x = s𝜎1 alors 𝑡 = 0 et x = s𝜎2 alors 𝑡 = 1.

La trace d’une fonction de forme P1 sur 𝜎 est la « fonction chapeau » 1D classique. Plus précisément :

𝜙𝜎1 (x) = 𝜙𝜎1 (x(𝑡)) = 𝜑^1 (𝑡) = 1 − 𝑡

𝜙𝜎2 (x) = 𝜙𝜎2 (x(𝑡)) = 𝜑^2 (𝑡) = 𝑡

À l’inverse, connaissant 𝑡 on peut retrouver le point x :

x(𝑡) = (1 − 𝑡)s1 + 𝑡s2 .

Nous avons une transformation bijective entre 𝜎 et le segment [0, 1] dit de référence et noté 𝜎
^ . Nous pouvons opérer un
changement de variable (𝑖, 𝑗 = 1, 2) :
∫︁ ∫︁ 1 ∫︁ 1
𝜙𝜎𝑗 𝜙𝜎𝑖 dx = |𝜎| 𝜙𝜎𝑗 (x(𝑡))𝜙𝜎𝑖 (x(𝑡))d𝑡 = |𝜎| 𝜑^𝑗 (𝑡)𝜑^𝑖 (𝑡)d𝑡.
𝜎 0 0

Les coefficients de masse de bord se calculent alors aisément et on obtient :

(︂ )︂
|𝜎| 2 1
𝑀𝜎𝑒 (𝑖, 𝑗) = .
6 1 2

4.3. Condition de Fourier 55

Maillage et Éléments Finis

56 Chapitre 4. Conditions aux bords

 CHAPITRE 5

Avancée

Und der Haifisch der hat Tränen

Und die laufen vom Gesicht
Doch der Haifisch lebt im Wasser
So die Tränen sieht man nicht
In der Tiefe ist es einsam
Und so manche Zähre fließt
Und so kommt es dass das Wasser
In den Meeren salzig ist

---Rammstein - Haifisch

5.1 Erreur comise et convergence

5.1.1 Méthode de Galerkin : Erreur commise

Nous considérons ici une formulation variationnelle abstraite :

Trouver 𝑢 ∈ 𝑉 tel que
{︂
(1)
∀𝑣 ∈ 𝑉, 𝑎(𝑢, 𝑣) = ℓ(𝑣).

Les formes continues 𝑎(·, ·) et ℓ(·) sont respectivement bilinéaire et linéaire, et 𝑎(·, ·) est de plus coercive. De cette
manière, le Théorème de Lax-Milgram s’applique et le problème (1) admet une unique solution. Nous noterons (·, ·)𝑉 et
‖·‖𝑉 respectivement le produit scalaire et la norme sur 𝑉 .
Nous restons dans un cadre abstrait et introduisons 𝑉ℎ , un sous-espace fonctionnel de 𝑉 , de Hilbert et de dimension finie.
Nous appliquons la méthode de Galerkin pour obtenir la formulation faible « approchée » :
Trouver 𝑢ℎ ∈ 𝑉ℎ tel que
{︂
(2)
∀𝑣ℎ ∈ 𝑉ℎ , 𝑎(𝑢ℎ , 𝑣ℎ ) = ℓ(𝑣ℎ ).

57
Maillage et Éléments Finis

Nous quantifions maintenant l’erreur commise en approchant 𝑢 par 𝑢ℎ , ou plus exactement, 𝑉 par 𝑉ℎ . Notons une pro-
priété très intéressante de la solution approchée 𝑢ℎ :

Lemma 5.1
Soit 𝑢 la solution exacte (i.e. solution de eqref{eq3-pbmodel}) et 𝑢ℎ la solution approchée (i.e. solution de (2)). Soit
𝑒ℎ = 𝑢 − 𝑢ℎ est l’erreur d’approximation, alors nous avons l’égalité suivante

∀𝑣ℎ ∈ 𝑉ℎ , 𝑎(𝑒ℎ , 𝑣ℎ ) = 0.

Proof. Comme 𝑉ℎ ⊂ 𝑉 , nous pouvons choisir 𝑣 = 𝑣ℎ dans la formulation variationnelle eqref{eq3-pbmodel} :

∀𝑣ℎ ∈ 𝑉ℎ , 𝑎(𝑢 − 𝑢ℎ , 𝑣ℎ ) = 𝑎(𝑢, 𝑣ℎ ) − 𝑎(𝑢ℎ , 𝑣ℎ )

= ℓ(𝑣ℎ ) − ℓ(𝑣ℎ )
=0

Remark 5.1
Si 𝑎(·, ·) est symmétrique, le lemme précédent implique que l’erreur est orthogonale à 𝑉ℎ par rapport au produit scalaire
𝑎(·, ·).

Nous pouvons maintenant montrer que l’erreur d’approximation 𝑢ℎ de 𝑢 est uniformément bornée par la distance entre 𝑢
et l’espace 𝑉ℎ . Ce résultat est connu comme étant le Lemme de Céa, démontré par Jean Céa durant sa thèse, en 1964.

Lemma 5.2 (de Céa)

Soit 𝑢 la solution exacte (i.e. solution de eqref{cv-pbmodel}) et 𝑢ℎ la solution approchée (i.e. solution de (2)). Nous avons
𝑀
‖𝑢 − 𝑢ℎ ‖𝑉 ≤ inf ‖𝑢 − 𝑣ℎ ‖𝑉 ,
𝛼 𝑣ℎ ∈𝑉ℎ
où 𝑀 et 𝛼 sont respectivement les constantes de continuité et de coercivité de 𝑎(·, ·) qui apparaissent dans le Théorème
de Lax-Milgram.

Proof. Pour 𝑣ℎ ∈ 𝑉ℎ , la quantité 𝑣ℎ − 𝑢ℎ est aussi un élément de 𝑉ℎ , ce qui implique d’après le lemme précédent que

𝑎(𝑢 − 𝑢ℎ , 𝑢 − 𝑢ℎ ) = 𝑎(𝑢 − 𝑢ℎ , 𝑢 − 𝑣ℎ + 𝑣ℎ − 𝑢ℎ )
= 𝑎(𝑢 − 𝑢ℎ , 𝑢 − 𝑣ℎ ) + 𝑎(𝑢 − 𝑢ℎ , 𝑣ℎ − 𝑢ℎ )
= 𝑎(𝑢 − 𝑢ℎ , 𝑢 − 𝑣ℎ ).

La coercivité et la continuité de 𝑎(·, ·) impliquent que

2
∀𝑣ℎ ∈ 𝑉ℎ , 𝛼 ‖𝑢 − 𝑢ℎ ‖𝑉 ≤ |𝑎(𝑢 − 𝑢ℎ , 𝑢 − 𝑢ℎ )|
≤ |𝑎(𝑢 − 𝑢ℎ , 𝑢 − 𝑣ℎ )|
≤ 𝑀 ‖𝑢 − 𝑢ℎ ‖𝑉 ‖𝑢 − 𝑣ℎ ‖𝑉 .

Nous en déduisons le résultat cherché :

𝑀
∀𝑣ℎ ∈ 𝑉ℎ , ‖𝑢 − 𝑢ℎ ‖𝑉 ≤ ‖𝑢 − 𝑣ℎ ‖𝑉 .
𝛼

58 Chapitre 5. Avancée
 Maillage et Éléments Finis

Remark 5.2
Le point important du Lemme de Céa est de remplacer le problème d’estimation de l’erreur par un problème d’approxi-
mation. En effet, il nous suffit de montrer que la solution est « bien approchée » par les fonctions de 𝑉ℎ pour savoir que
l’erreur ne sera qu’une constante fois plus grande que cette erreur d’approximation.

Nous pouvons maintenant donner une condition pour que la méthode de Galerkin converge.

Lemma 5.3
Soit Πℎ : 𝑉 → 𝑉ℎ un opérateur d’interpolation tel que

∀𝑣 ∈ 𝑉, lim ‖𝑣 − Πℎ 𝑣‖𝑉 = 0,
ℎ→0

alors la méthode de Galerkin converge, c’est-à-dire :

lim ‖𝑢 − 𝑢ℎ ‖𝑉 = 0.
ℎ→0

Proof. C’est une conséquence directe du lemme de Céa, puisque :

𝑀
0 ≤ ‖𝑢 − 𝑢ℎ ‖𝑉 ≤ ‖𝑢 − Πℎ 𝑢‖𝑉 → 0 (ℎ → 0).
𝛼

La propriété demandée à l’opérateur de projection Φℎ : 𝑉 → 𝑉ℎ est assez naturel : plus ℎ est petit et plus le projeté d’une
fonction 𝑣 ∈ 𝑉 doit être proche de 𝑣. Nous pouvons voir cela comme l’espace 𝑉ℎ est « proche » de 𝑉 .

5.1.2 Converge des éléments finis P1

Afin de montrer que la méthode des éléments finis P1 converge, nous devons obtenir un opérateur d’interpolation et
montrer qu’il vérifie l’hypothèse nécessaire du lemme 5.1.1.
TODO : !

5.2 Éléments Finis P2

5.2. Éléments Finis P2 59

Maillage et Éléments Finis

60 Chapitre 5. Avancée
Deuxième partie

Implémentation

61
 CHAPITRE 6

Maillage avec GMSH

6.1 Prise en main de GMSH

Suivez ce tutoriel GMSH pour apprendre à :

— Générer de géométries simples (e.g. stade)
— Générer de géométries plus compliquées avec OpenCascade (e.g. tasse)
— Gerer des labels Physical
— Débuter avec l’API Python

6.2 API GMSH

Déjà abordé dans le tutoriel GMSH, nous nous intéressons à l’API Python de GMSH.

6.2.1 Fonction de Forme

Le but est d’afficher une fonction de forme P1 − Lagrange, c’est à dire une fonction 𝜙𝐼𝐽 de 𝑉ℎ qui vaut 0 sur tous les
sommets 𝐽 ̸= 𝐼 du maillage sauf sur le sommet 𝐼, pour lequelle la fonction prend la valeur 1 :

∀𝐼, 𝐽 = 0, . . . , 𝑁𝑠 − 1, 𝜙𝐼 (s𝐽 ) = 𝛿𝐼𝐽

Exercise 1.1
À l’aide de l’API Python de GMSH :
— Générez un carré unitaire avec un pas de maillage de 0.25
— Appliquez des labels Physical : un pour la surface et un pour son bord
— Générez le maillage 2D
— Construisez 3 numpy array Phi, X et Y unidimensionnel et de taille le nombre de sommets 𝑁𝑠 du maillage tels
que :

63
Maillage et Éléments Finis

— Phi est le vecteur nul sauf en un coefficient où Phi[I] = 1 (choisissez le I)

— X et Y sont respectivement les coordonnées x et y des points du maillage (voir ci-dessous)
— Affichez le tout à l’aide de Matplotlib et de la projection 3D

Pour obtenir les coordonnées des points d’un groupe Physical donné, vous pouvez utilisez [Link].
getNodesForPhysicalGroup(dim, tag) (voir [Link]) :

def getNodesForPhysicalGroup(dim, tag):

"""
Get the nodes from all the elements belonging to the physical group of
dimension `dim' and tag `tag'. `nodeTags' contains the node tags; `coord
˓→'

is a vector of length 3 times the length of `nodeTags' that contains␣

˓→the x,

y, z coordinates of the nodes, concatenated: [n1x, n1y, n1z, n2x, ...].

Return `nodeTags', `coord'.

"""

Remark 1.1
— GMSH commence la numérotation des sommets à 1
— La liste retournée par getNodesForPhysicalGroup n’est pas triée

6.2.2 Interpolation P1

Prenons une fonction f définie sur Ω. Une interpolation possible de f sur l’espace P1 est la fonction Πℎ 𝑓 telle que Πℎ 𝑓 (s) =
𝑓 (s) pour chaque sommet s du maillage.

Exercise 1.2
Construisez l’interpollée Πℎ 𝑓 de la fonction 𝑓 (𝑥, 𝑦) = sin(𝜋𝑥) sin(𝜋𝑦)

64 Chapitre 6. Maillage avec GMSH

 CHAPITRE 7

Solveur FEM Python

7.1 Matrices Creuses

Le but maintenant est d’implémenter la méthode des éléments finis P1 en 2D, autrement dit, à un maillage donné, de
calculer :
— Les matrices de masse 2D (union de triangles) et 1D (union de segments)
— Les matrices de rigité (union de triangles) et 1D (union de segments)
— Appliquer les conditions de Dirichlet éventuelles
— Approcher des intégrales par des quadratures adaptées
La matrice du système finale sera stockée sous le format COOrdinate. L’algorithme d’assemblage se prête particulièrement
bien au stockage COO puisque la redondance des coefficients est autorisée par Scipy. Autrement dit, les fonctions calculant
les matrices de masse et de rigidité retourneront des triplets de type [I, J, Valeur] où I est l’indice ligne, J l’indice
colonne et Valeur le coefficient (potentiellement partiel).
Avant de résoudre le système, la matrice sera transformée au format Compressed Storage Raw (CSR).

7.1.1 Matrice COO

Rappels

Le format COO (COOrdinates) propose de ne stocker que les coefficients non nuls d’une matrice A sous la forme de 2
listes d’entiers row et col et une liste de réels val, toutes trois de même taille, telles que

A[row[i], col[i]] = val[i]

Les trois listes peuvent aussi être combinées en une seule liste data de triplets de type [int, int, double] et
telle que :

data[i] = [row[i], col[i], val[i]]

Le format COO est assez permissif :

65
Maillage et Éléments Finis

— Redondance : si deux triplets possèdent les mêmes indices ligne et colonne alors le coefficient de la matrice, associé
à ces indices ligne et colonne, sera obtenu en sommant les valeurs de ces triplets.
— Pas d’ordre : aucune nécessité d’ordonner les triplets selon les lignes ou les colonnes
Par exemple, les deux jeux de listes ci-dessous décrivent la même maitrce, la seule différence est que le coefficient en
(2,2) est scindé en deux et le premier et le dernier triplet ont été permutés :

i 0 1 2 3 4 5 6
row[i] 0 0 1 2 3 3 3
col[i] 0 3 1 2 0 1 3
val[i] 1.1 2 1 2.3 0.5 2 2

i 0 1 2 3 4 5 6 7
row[i] 3 0 1 2 2 3 3 0
col[i] 3 3 1 2 2 0 1 0
val[i] 2 2 1 1.3 1 0.5 2 1.1

Scipy

Pour construire une matrice COO dans Scipy à partir d’un jeu de données data, ce dernier doit être de type ([],
([],[])) : un Tuple contenant une List (val) ainsi qu’un Tuple contenant deux List (row et col) (list vs.
tuple ?) :

data = (val, (row, col))

Une matrice COO se construit alors ainsi

A = coo_matrix(data) # si data dans le format ci-dessus

A = coo_matrix((val, (row, col))) # si row, col et val sont séparées

Remark 2.15
Une matrice COO peut être visualisée en se transformant en array avec la méthode toarray() :

print([Link]())

Triplets

Nous proposons de construire notre future matrice par concaténation de triplets de type (I, J, valeur). une classe Tri-
plets qui encapsule cette structure de données. Nous lui adjoignons une méthode append permettant d’ajouter un
triplet au bout des autres :

Triplets t; print([Link]) # ([], ([], []))

[Link](0, 1 ,2.); print([Link]) # ([2.], ([0], [1]))
[Link](3, 4 ,5.2); print([Link]) # ([2., 5.2], ([0, 3], [1, 4]))

La classe ressemble alors à cela :

66 Chapitre 7. Solveur FEM Python

 Maillage et Éléments Finis

def class Triplet:

def __init__():
[Link] = ([], ([], []))
def __str__():
return str([Link])
def append(self, I, J, val):
# Ajoute le triplet [I, J, val] dans [Link]
# ...

Exercise 2.1
Construisez la classe Triplet et implémentez la méthode append. N’oubliez pas de tester votre classe.

Exercise 2.2
Testez votre classe Triplet en construisant la matrice suivante (au format COO évidemment) :

⎛ ⎞
1.1 0 0 2
⎜ 0 1 0 0⎟
𝐴=⎜ ⎟
⎝ 0 0 2.3 0⎠
0.5 2 0 2

7.1.2 Format CSR

Une fois la matrice au format COO construite, nous la transformerons au format CSR par la méthode tocsr() :

A = coo_matrix((val, (row, col))).tocsr()

7.2 Gestion du maillage

7.2.1 Une classe par élément

Nous proposons de construire 3 classes : Point, Segment et Triangle représentant un élément de type point,
segment et triangle. La première étape consiste à transcrire les informations du maillage dans notre structure de données.
Nous pouvons aussi ne pas passer par cette étape et n’utiliser que GMSH, cependant, cela nous donnera un meilleur
contrôle par la suite.
La classe Point contient les coordonnées du point et un identifiant (0,1,2,3, …), son indice globale, qui correspondra
en P1 à la ligne dans la matrice du degré de liberté (DOF). Les classes Segment et Triangle ont pour paramètres :
— Un identifiant (0,1,2,3…)
— Une liste orientée de Point qui le définissent
— Un tag Physical (ou -1 sinon)
L’identifiant n’a pas besoin d’être le même que dans GMSH, il est même plus logique que les identifiants soient consécutifs
et incrémenté à chaque nouvelle instance (0, 1, 2, …). Le tag Physical doit en revanche être le même que GMSH pour
éviter toute confusion. Vous êtes évidemment libre d’ajouter des paramètres et des méthodes à ces classes !
Par ailleurs, les 2 classes d’élément Segment et Triangle, disposeront au moins des méthodes suivantes :

7.2. Gestion du maillage 67

Maillage et Éléments Finis

— area() : calcule l’aire de l’élément (pour un segment, sa taille ; peut être stockée dans un paramètre et être
calculé une fois pour toute)
— jac() : calcule son jacobien (son aire pour un segment, 2 fois l’aire pour un triangle)

7.2.2 Une classe pour les gouverner toutes…

Il peut être assez malin de construire une classe Mesh qui représente le maillage et permettra d’effectuer des recherches
d’éléments. Cette classe aura pour paramètre trois listes : une liste de Point, une de Segment et une de Triangle.
De plus, nous lui ajoutons deux méthodes :
— [Link](dim, physical_tag) : retourne une liste de tous les éléments ayant la dimension
dim (=1 pour segment, 2 pour triangle) et le tag physique physical_tag
— [Link](dim, physical_tag) : retourne uniquement les points du domaine de dimension
dim et de label physique physical_tag

7.2.3 …Et dans les ténébres les lier (à GMSH)

Le maillage étant construit avec GMSH, il s’agit maintenant de convertir les données issues de GMSH dans notre structure.
Nous supposons ici que le maillage est déjà chargé en mémoire, soit parce qu’il vient d’être construit par GMSH (via
[Link]) soit parce qu’il a été lu sur disque (via [Link]("[Link]")). Nous proposons
que le constructeur de Mesh soit vide (ne construit rien) et que l’instance de Mesh soit construite à l’aide d’une méthode
GmshToMesh(string filename), qui lit les données GMSH et construit les Point, Segment et Triangle
et les ajoute dans les listes correspondent à Mesh. L’argument filename peut être optionnel si le maillage est déjà
construit.

Remark 2.16
La méthode la plus délicate à construire est GmshToMesh. Pour vous aider un petit peu, n’hésitez pas à fouiller dans
l’API de GMSH) :
— [Link]() : retourne tous les noeuds
— [Link]() : retourne tous les groupes physiques avec leurs dimension et tag
— [Link](dim, physical_tag) : retourne toutes les enti-
tés d’un groupe physique
— [Link](dim, tag) : retourne tous les éléments de dimension dim (seg-
ments (dim=1), triangles (dim=2), …) appartenant à l’entity de label tag.
L’algorithme ressemble surement à ceci :

// Création des points

For every Nodes
Point(...)
End
//Création des éléments
For every Physical Entity
For every Entity
For every Element

68 Chapitre 7. Solveur FEM Python

 Maillage et Éléments Finis

7.3 Matrices de Masse et de Rigidité

7.3.1 Rappel : Algorithme d’assemblage

Plutôt que de calculer les coefficients un à un de la matrice 𝐴, l’algorithme d’assemblage propose de parcours chaque
élément et d’ajouter leur contribution élémentaire à la grande matrice. Avec notre notation en Triplets cela donne :

Triplets t;
For p = 1, ... Nt // Parcours des Triangles
Mp = MatElem(p); // Matrice Elementaire du triangle p
For i = 1,2,3
I = Loc2Glob(p, i);
For j = 1,2,3
J = Loc2Glob(p,j);
[Link](I, J, Mp(i,j)); // contribution élémentaire
End
End

Si les élements parcourus sont ne sont pas des triangles à 3 points (segments, tétrahèdres, …), il suffit d’adapter le
pseudo-code ci-desssus. Nous devons donc implémenter les calculs des matrices de masse et de rigidité élémentaire pour
chaque élément.

7.3.2 Matrices de masse élémentaires

Construisez une fonction mass_elem prenant en argument un Segment ou un Triangle, un Triplets et un

scalaire optionnel :

# element = Segment ou Triangle ; triplets = Triplets ; alpha un scalaire optionnel

def mass_elem(element, triplets, alpha =1.):
# ...
# return triplets

Cette fonction calcule [les coefficients de la matrice élémentaire]([Link]

implementation_matrices_elementaires/) de l’élément (selon son type) et les ajoute à triplets.

Remark 2.17
Pour un élément donné, son type (Segment ou Triangle) est donné par son paramètre name.

7.3.3 Matrice de masse globale

Nous proposons de construire une fonction qui calcule toutes les contributions élémentaires de la matrice de masse d’un
domaine de tag Physical et de dimension dim issue d’un maillage msh. Les coefficients partiels seront ajoutés sous
forme de Triplet dans une liste envoyée en argument. Nous séparons pour le moment les calculs de la matrice de
masse de ceux de la matrice de rigidité :

# msh = Mesh, dim = int, physical_tag = int, triplets = Triplets

def Mass(msh, dim, physical_tag, triplets):
#...

Exercise 2.3

7.3. Matrices de Masse et de Rigidité 69

Maillage et Éléments Finis

Au boulot ! Assurez vous que la matrice de masse globale 𝑀 associée au domaine Ω vérifie la relation suivante

𝑈 𝑇 𝑀.𝑈 = |Ω|, 𝑈 = [1, 1, 1, . . . , 1]𝑇 .

7.3.4 Matrice de Rigidité

Pour les matrices de rigidité, il faut calculer des quantités supplémentaires, comme la matrice 𝐵_𝑝 ou les gradients des
fonctions de forme, par exemple :

def gradPhi(element, i:int):

# ...

Exercise 2.4
Ajoutez les fonctionnalités dans votre code permettant de calculer les contributions élémentaires des matrices de rigidité
puis la matrice globale.
Vérifiez que votre matrice de rigidité 𝐷 satisfait bien la relation suivante :

𝐷𝑈 = 0, 𝑈 = [1, 1, 1, . . . , 1]𝑇 .

7.4 Quadratures

7.4.1 Rappel

Certaines intégrales ne peuvent être calculées analytiquement et devront être approchées numériquement via des règles
de quadrature. Prenons pour exemple d’une fonction 𝑓 quelconque mais connue, le second membre 𝐵 (un vecteur) sera
alors de la forme suivante, où x = (𝑥, 𝑦),
∫︁ ∑︁ ∑︁ ∫︁
𝐵[𝐼] = 𝑓 (x)𝜙𝐼 (x) dx = 𝑓 (x)𝜙𝑝𝑖 (x) dx
Ω 𝑝 𝑖 𝐾𝑝
∑︁ ∑︁ ∫︁
= |det(𝐽𝑝 )| ̂︀𝑖 (𝜉, 𝜂) d(𝜉, 𝜂)
𝑓 (x(𝜉, 𝜂))𝜙 (1)
^
𝐾
𝑝 𝑖
∑︁ ∑︁ ∑︁
≈ |det(𝐽𝑝 )| 𝜔𝑚 𝑓 (x(𝜉𝑚 , 𝜂𝑚 ))𝜙
̂︀𝑖 (𝜉𝑚 , 𝜂𝑚 )
𝐾𝑝 𝑖 𝑚

Les poids 𝜔𝑚 et les points de quadrature (𝜉𝑚 , 𝜂𝑚 ) dépendent de la précision recherchée comme expliqué dans le cours.
Rappelons aussi que x(𝜉𝑚 , 𝜂𝑚 ) s’obtient par les fonctions d’interpolation géométrique, qui dans le cas d’éléments finis
isoparamétriques, sont les mêmes que les fonctions éléments finis P1 , c’est à dire que pour x appartenant à un élément de
sommets (s𝑖 )𝑖 :
∑︁
x(𝜉𝑚 , 𝜂𝑚 ) = 𝜓^𝑖 (𝜉𝑚 , 𝜂𝑚 )s𝑖 ,
𝑖=1

Les fonctions (𝜓^𝑖 )𝑖 sont égales au fonctions de forme P1 − Lagrange (𝜙

̂︀𝑖 )𝑖 . Cependant, les fonction (𝜓^𝑖 )𝑖 sont des fonctions
d’interpolation géométriques tandis que les fonctions de forme (𝜙^𝑖 )𝑖 sont des fonctions d’interpolation de la solution.

70 Chapitre 7. Solveur FEM Python

 Maillage et Éléments Finis

7.4.2 Méthodes et fonctions utiles

Pour chaque type d’élément, Triangle ou Segment, nous avons besoin de méthodes permettant d’obtenir les quantités
suivantes : :
1. Les poids des points de quadrature
2. Les coordonnées paramétriques des points de quadrature
3. Les coordonnées physiques des points de quadrature
4. Les valeurs des fonctions de forme de référence 𝜙^ sur des coordonnées paramétriques
Remarquez que les points 1, 2 ne dépendent que du type de l’élément considéré.

Exercise 2.5
Afin de bien compartimenter chaque fonctionnalité, nous proposons :
— Ajouter aux classes Triangle et Segment la méthode def gaussPoint(self,order=2): qui re-
tourne, dans le format de votre choix, les poids, les coordonnées paramétriques et les coordonnées physiques des
points de Gauss de l’élement considéré et pour une précision order. Vous aurez sans doute besoin de méthodes
intermédiaires pour calculer, par exemple les 𝜓^𝑖 (𝜉, 𝜂).
— Ajouter une fonction def phiRef(element, i:int, param:[float]): qui calcule 𝜙^𝑖 (𝜉, 𝜂) sur un
élément Segment ou Triangle. L’argument param est une liste des coordonnées paramétriques ((𝜉, 𝜂) pour
un triangle, 𝑠 pour un segment))

7.4.3 Intégrale

Construisez maintenant une fonction de prototype suivant

def Integrale(msh:Mesh, dim:int, physical_tag:int, f, B:[Link], order=2):

Cette fonction calcul l’intégrale 𝑓 𝜙𝐼 sur le domaine de tag physique physical_tag et de dimension dim. Le résultat
∫︀

est alors ajouté dans B[I] (voir (1)). L’argument f sera une fonction décrite par l’utilisatrice/utilisateur, elle prendra 2
arguments, x et y et retournera un scalaire correspondant à 𝑓 (𝑥, 𝑦).

7.5 Conditions de Dirichlet

Pour prendre en compte les éventuelles condition de Dirichlet, nous avons besoin d’une fonction de prototype suivant

def Dirichlet(msh, dim, physical_tag, g, triplets, B):

Cette fonction prend comme argument le Triplets triplets et le vecteur B du système linéaire et les modifie pour
prendre en compte la condition de dirichlet 𝑢 =g sur le domaine de dimension dim et de tag physique physical_tag.
La technique utilisée pour forcer cette condition est celle vue en cours.
Pour cela, nous parcourons les noeuds I du domaine de Dirichlet. Puis, dans la liste des indices ligne de triplets,
dès qu’un occurence à I est obtenu, la valeur de ce triplet est mise à 0. Il ne faut pas oublier, à la fin, d’ajouter un triplet
(I,I,1) correspondant au terme diagonal et de modifier le coefficient B[I] = g(x,y).

Remark 2.18
Cette technique n’est peut être pas la plus optimale ! Mais elle a le mérite de fonctionner…

7.5. Conditions de Dirichlet 71

Maillage et Éléments Finis

7.6 Résolution et Analyse

7.6.1 Problème de référence

Résumons ici l’utilisation de notre programme éléments finis sur le problème suivant :
{︂
−Δ𝑢 + 𝑢 = 𝑓 (Ω)
(2)
𝑢 = 0 (𝜕Ω)

La formulation variationnelle est donnée par

⎧
⎨ Trouver 𝑢 ∈ 𝐻01 (Ω)∫︁ tel que ∫︁ ∫︁
⎩ ∀𝑣 ∈ 𝐻01 (Ω), ∇𝑢 · ∇𝑣 + 𝑢𝑣 = 𝑓𝑣
Ω Ω Ω

Pour simplifier nous prenons Ω =]0, 1[×]0, 1[ le carré unitaire et 𝑓 (𝑥, 𝑦) = (1 + 2𝜋 2 ) sin(𝜋𝑥) sin(𝜋𝑦) de sorte que la
solution exacte est connue et vaut

𝑢(𝑥, 𝑦) = 𝑔(𝑥, 𝑦) = sin(𝜋𝑥) sin(𝜋𝑦).

Fıg. 1 – Solution

72 Chapitre 7. Solveur FEM Python

 Maillage et Éléments Finis

7.6.2 Résolution

Dans notre programme, cela reviendra à écrire quelque chose comme

#import ...

#Données
def g(x,y):
return [Link]([Link]*x)*[Link]([Link]*y)
def f(x,y):
return g(x,y)*(2*[Link]*[Link] +1 )
def diri(x,y):
return 0.
#Maillage
msh = [Link]("[Link]")
# Triplets
t = [Link]()
fem_p1.Mass(msh, 2,10, t)
fem_p1.Stiffness(msh, 2,10, t)
b = [Link](([Link],))
fem_p1.Integral(msh, 2, 10, f, b, 2)
fem_p1.Dirichlet(msh, t, b, 1, 1, diri)
# Résolution
A = (sparse.coo_matrix([Link])).tocsr()
U = [Link](A, b)

# Visualisation
x= [pt.x for pt in [Link]]
y= [pt.y for pt in [Link]]
connectivity=[]
for tri in [Link]:
[Link]([ [Link] for p in tri.p])

[Link](x, y, connectivity, U, 12)

[Link]()
[Link]()

### U de référence
Uref = [Link](([Link],))
for pt in [Link]:
I = int([Link])
Uref[I] = g(pt.x, pt.y)
[Link](x, y, connectivity, Uref, 12)
[Link]()
[Link]()

7.6.3 Convergence

Exercise 2.6

1. Pour différents pas de maillage, calculez l’erreur en norme 𝐿2 entre la solution exacte et la solution approchée
pour le problème (2).
2. Affichez la courbe de l’erreur en fonction de ℎ en échelle log-log.
3. Calculez la pente de la courbe et déduisez-en la vitesse de convergence par rapport au pas de maillage (ℎ).
Sauvegardez par ailleurs une copie de la courbe en format données (JSON ou autre) ou image (PNG par exemple,

7.6. Résolution et Analyse 73

Maillage et Éléments Finis

pas de JPG nous ne sommes pas des sauvages !).

74 Chapitre 7. Solveur FEM Python

Troisième partie

Projet

75
 CHAPITRE 8

2017 - 2018 : Wi-Fi

Ce projet propose d’étudier la propagation des ondes Wi-Fi dans un appartement de 35m2 . Nous ferons appels aux logiciels
GMSH, pour la gestion du maillage et la visualisation, et FreeFem++ pour la résolution par éléments finis. Le sujet est
décomposé en trois parties : la théorie (le modèle), l’implémentation et l’étude des résultats.

8.1 Équation de Helmholtz

La définition Wikipédia d’une onde est la suivante :

Une onde est la propagation d’une perturbation produisant sur son passage une variation réversible des pro-
priétés physiques locales du milieu. Elle se déplace avec une vitesse déterminée qui dépend des caractéris-
tiques du milieu de propagation. Une onde transporte de l’énergie sans transporter de matière.
Mathématiquement, une onde E (x, 𝑡) dépend du temps 𝑡 et de l’espace x, et vérifie l’équation des ondes :

1 𝜕2E
ΔE (x, 𝑡) = (x, 𝑡),
𝑐2 𝜕𝑡2
𝜕2
∑︀3
où Δ = 𝑗=1 𝜕𝑥2𝑗 est l’opérateur Laplacien (spatiale) et 𝑐 est la célérité de l’onde dans le milieu (peut dépendre de x !).
Par exemple, dans le cas d’une onde électromagnétique et dans le vide, 𝑐 est la célérité de la lumière, soit 299792458m.s−1 .
La quantité x est un vecteur de dimension 𝑑 = 2 ou 𝑑 = 3 selon le problème considéré : dans notre cas 𝑑 = 2.
c’est-à-dire lorsque la pulsation 𝜔 (en rad.s−1 ) de l’onde est fixée, l’onde s’écrit alors
Lors d’une excitation périodique,√
E (x, 𝑡) = ℜ (𝑢(x)𝑒−𝚤𝜔𝑡 ) où 𝚤 = −1 et 𝐸 est une onde spatiale satisfaisant l’équation de Helmholtz :

𝜔2
Δ𝐸 + 𝐸 = 𝑓.
𝑐2
Cette équation s’obtient en remplaçant E (x, 𝑡) par 𝐸(x)𝑒−𝚤𝜔𝑡 dans l’équation des ondes. Nous notons en général 𝑘 = 𝜔𝑐
(en rad.m−1 ) le nombre d’onde et 𝜆 = 2𝜋 𝑘 (en m) la longueur d’onde, autrement dit, la distance entre deux amplitudes,
de sorte que l’équation de Helmholtz s’écrit

Δ𝐸 + 𝑘 2 𝐸 = 𝑓.

77
Maillage et Éléments Finis

Les ondes Wi-Fi qui suivent la norme IEEE 802.11g sont émises à une fréquence variant de 2.4GHz à 2.5GHz. Dans
notre projet, nous nous limiterons à des ondes de fréquence 2.4GHz textbf{si votre machine vous le permet !} En effet,
les simulations de propagation d’ondes sont parmis les plus coûteuses en terme de CPU. Si votre machine n’est pas assez
puissante, nous prendrons une onde de fréquence plus faible, comme 1GHz voire moins.

Exercise 1.3
Sachant que 𝜔 = 2𝜋𝐹 où 𝐹 est la fréquence, en Hertz (Hz), calculez le nombre d’onde 𝑘 et la longueur d’onde 𝜆 dans le
vide, pour une onde électromagnétique et pour 𝐹 = 2.4GHz et 𝐹 = 1GHz.

8.2 Modèle

8.2.1 Entre les murs

Nous notons Ω l’appartement tout entier dans lequel est situé notre routeur. Les murs sont supposés être du même maté-
riau : du placo-plâtre. Le domaine Ω est décomposé en deux domaines :

Ω = Ωair ∪ Ωmur ,

où Ωair est l’intérieur de l’appartement, composé d’air (nous le supposerons vide de meubles), et Ωmur contient les murs
en placo-plâtre. Nous ne prendrons pas en compte les appartements voisins au nôtre. De plus, nous nous limiterons à la
dimension 2 de l’espace.
L’air est modélisé comme étant le vide : 𝑐air = 𝑐, tandis que, pour les murs en placo-plâtre, nous avons 𝑐placo = 2.4
𝑐
. Plutôt
que de rendre le nombre d’onde dépendant de l’espace, nous introduisons la fonction de contraste 𝑛 définie par :

si x ∈ Ωair (i.e. x est dans l'air),

{︂
1
𝑛(x) =
2.4 si x ∈ Ωmur (i.e. x est dans le mur).

Nous modifions alors l’équation de Helholtz ainsi :

Δ𝐸(x) + 𝑘 2 𝑛(x)2 𝐸(x) = 𝑓 (x) ∀x ∈ Ω.

Le nombre d’onde 𝑘 est donc ici celui de l’air (ou du vide).

8.2.2 Le routeur : la fonction 𝑓

Idéalement, le routeur devrait être modélisé comme un point source et donc mathématiquement par la distribution de
Dirac 𝛿x𝑠 centré sur la position x𝑠 du routeur. L’équation à résoudre serait alors :

Δ𝐸(x) + 𝑘 2 𝑛(x)2 𝐸(x) = −𝛿x𝑠 .

Cependant nous ne pouvons pas utiliser de distribution de Dirac avec la méthode des éléments finis (car 𝛿x𝑠 ̸∈ 𝐿2 (Ω)).
Nous modélisons alors le routeur par un disque Ω𝑠 de rayon 𝜀𝑠 = 0.1 et de centre x𝑠 = (𝑥𝑠 , 𝑦𝑠 ) et la distribution 𝛿x𝑠 est
approchée par la fonction chapeau 𝑓𝑠 définie par

‖x𝑠 − x‖
⎧ (︂ )︂
⎨ 3
1− si x ∈ Ω𝑠 ,
∀x ∈ Ω, 𝑓𝑠 (x) = 𝜋𝜀2𝑠 𝜀𝑠
0 sinon.
⎩

La fonction 𝑓𝑠 est d’intégrale totale égale à 1 et quand 𝜀𝑠 tend vers 0, la suite de fonction (𝑓𝑠 )𝜀𝑠 tend vers la distribution
de Dirac 𝛿x𝑠 .

78 Chapitre 8. 2017 - 2018 : Wi-Fi

 Maillage et Éléments Finis

Exercise 1.4
Montrez que :
∫︁
𝑓𝑠 (x)dx = 1.
Ω

Pour simplifier, on pourra se placer dans le cas où Ω𝑠 est centré en (0, 0) et utiliser un changement de coordonnées
adéquat.

8.2.3 Réflexions parasites

Sur le bord extérieur de notre appartement, nous souhaitons que l’onde ne se réfléchisse pas, pour cela nous imposons la
condition aux limites suivantes :

𝜕n 𝐸(x) − 𝚤𝑘𝑛(x)𝐸(x) = 0, sur 𝜕Ω.

Cette condition est une approximation de la condition non réfléchissante exacte : l’onde sera légèrement réfléchie, toutefois
la mise en oe{}uvre d’une telle condition reste très simple et peu coûteuse, ce qui explique que nous l’utilisions. Notez
que n est le vecteur normale unitaire sortant à Ω (ne pas confondre avec 𝑛 la fonction de constraste).

8.2.4 Le système à résoudre

Au final, le système que nous devons résoudre est le suivant. Nous cherchons l’onde spatiale textbf{complexe} 𝐸, telle
que, pour 𝑘 = 𝜔/𝑐 donné, elle vérifie :

Δ𝐸(x) + 𝑘 2 𝑛(x)2 𝐸(x) = −𝑓𝑠 (x) dans Ω,

{︂
(1)
𝜕n 𝐸(x) − 𝚤𝑘𝑛(x)𝐸(x) = 0 sur 𝜕Ω.

8.2.5 Formulation variationnelle

Exercise 1.5
Montrez que la formulation variationnelle du problème (ref{eq:helmholtz}) s’écrit

Trouver 𝑢 ∈ 𝐻 1 (Ω) tel que

{︂
(2)
∀𝑣 ∈ 𝐻 1 (Ω), 𝑎(𝑢, 𝑣) = ℓ(𝑣),
avec
∫︁ ∫︁ ∫︁
𝑎(𝑢, 𝑣) = − ∇𝐸(x) · ∇𝑣(x)dx + 𝑘 2 𝑛(x)2 𝐸(x)𝑣(x)dx + 𝚤𝑘𝑛(x)𝐸(x)𝑣(x)dx,
Ω Ω 𝜕Ω
∫︁
ℓ(𝑣) = − 𝑓𝑠 (x)𝑣(x)dx.
Ω

Exercise 1.6
Montrez que les applications 𝑎(·, ·) et ℓ(·) sont continues sur respectivement 𝐻 1 (Ω) × 𝐻 1 (Ω) et 𝐻 1 (Ω).

8.2. Modèle 79
Maillage et Éléments Finis

Malheureusement pour nous, 𝑎(·, ·) n’est pas coercive : nous ne pouvons pas appliquer le Théorème de Lax-Milgram et
de ses corollaires. Nous admettrons le Théorème suivant

Theorem 1.1
Les problèmes (1) et (2) admettent une unique solution.

8.3 Implémentation

8.3.1 Le pas de maillage

Pour résoudre des problèmes de propagation d’ondes, le pas de maillage (ou la finesse de maillage), noté ℎ (diamètre du
plus grand élément), dépend du nombre d’onde 𝑘, ou plutôt, de la longueur d’onde 𝜆 = 2𝜋/𝑘. En effet, si nous ne prenons
pas assez de points de discrétisation, l’onde ne sera pas suffisamment approchée comme le montre la Figure 1. Cela rejoint
le Théorème de Shanon d’échantillonage.
En pratique, nous prenons un nombre de points 10 ≤ 𝑛𝜆 ≤ 20 par longueur d’onde. Autrement dit, nous avons, si ℎ est
la taille caractéristique d’un élement :

𝜆 2𝜋
ℎ= = .
𝑛𝜆 𝑘𝑛𝜆
Il faudra veiller, lors de nos simulations, à ce que 𝑛𝜆 = 10 textbf{au minimum}. Pour les simulations finales à haute
fréquence, il vaut mieux privilégier 𝑛𝜆 = 15 ou 𝑛𝜆 = 20, si la machine le permet.

Remark 1.2
Avant de lancer un calcul et/ou un maillage pour tester (i.e mon code plante-il ?) : choisissez un nombre d’onde faible !
Par exemple 𝑘 = 10 ou 𝑘 = 5 (mais pas 𝑘 = 0 !).

8.3.2 L’appartement : GMSH

Le plan de l’appartement que nous considérons est donné par la Figure 2 où chaque dimension est exprimée en mètre.
Nous proposons les contraintes suivantes sur l’appartement : sa longueur 𝐿 est fixée à 7, sa largeur ℓ à 5 et l’épaisseur 𝑑
des murs est égale à 0.2.

Exercise 1.7
Implémentez un code GMSH qui reproduit le plan de la figure 2, en respectant las consignes suivantes :
1. Nous devons être en mesure de pouvoir modifier via l’interface graphique les quantités suivantes :
— Le nombre de points 𝑛𝜆 de discrétisation par longueur d’onde 𝜆. N’oubliez pas que 𝑛𝜆 est un entier.
— La position du routeur x𝑠 (Il n’est pas nécessaire de contraindre/vérifier sa position par rapport aux murs
(mais par rapport aux dimensions de l’appartement, oui !), nous supposerons l’utilisateur suffisamment
malin pour cela)
— Remarque : dans la version originale du projet, il était demandé de pouvoir modifier les dimensions des pièces et
donc de pouvoir déplacer les murs.

2. La largeur des portes est fixées à 80cm.

3. Le rayon de Ω𝑠 est fixé à 0.1cm.

80 Chapitre 8. 2017 - 2018 : Wi-Fi

 Maillage et Éléments Finis

Fıg. 1 – Illustration d’un échantillonage trop faible en 1D : la solution discrète est linéaire et non oscillante ! Elle est trop
éloignée de la solution exacte. Nous avons besoin de plus de points de discrétisations.

8.3. Implémentation 81
Maillage et Éléments Finis

Fıg. 2 – Plan de l’appartement (en mètre). Rappelons que Ω = Ωair ∪ Ωmur et 𝜕Ω est le bord « externe » de l’appartement.

82 Chapitre 8. 2017 - 2018 : Wi-Fi

 Maillage et Éléments Finis

4. Tous les murs sont de même épaisseur 𝑑.

Remark 1.3
À vous de définir les entités Physical qui vous semblent d’importance ainsi que les dimensions des pièces.

8.3.3 La formulation variationnelle : FreeFem++

Exercise 1.8
Rédigez le code freefem permettant de résoudre le problème (2) dans le maillage réalisé par gmsh. Les quantités en sorties
qui nous intéressent sont :
— Le champ 𝑢 : partie réelle et partie imaginaire, dans tout l’appartement.
— La valeur absolue du champ : |𝑢|, dans tout l’appartement.

Remark 1.4
Plusieurs remarques :
— Vous pouvez rédiger un script qui modifie le nombre d’onde 𝑘 dans le(s) fichier(s) freefem et gmsh, puis lance
gmsh et freefem… En bref : vous êtes libre de geeker comme bon vous semble.
— En phase de test, rappelez-vous de textbf{ne pas lancer} de résolution du problème pour un 𝑘 élevé.
— Mettez la solution à l’échelle : faites en sorte que le maximum de la solution soit égale à 1.
— La solution que l’on calcule présente très certainement un pic au niveau du routeur, qui empêche de visionner
correctement la solution en dehors de celui-ci. Il peut, dès lors, être intéressant de regarder la solution obtenue
partout en dehors du routeur et/ou de modifier la range de couleur.

8.4 Étude et examen oral

Maintenant que nous disposons d’un code qui permet de résoudre le problème de propagation d’ondes Wi-Fi dans un
appartement, nous pouvons le tester.
Il n’y a ici plus d’exercice à proprement parler : à vous d’effectuer les simulations qui vous paraissent intéressantes et de
sauvegarder les résultats issus de ces résolutions en vue de les présenter à l’oral. La qualité des réponses et le nombre de
questions traitées seront pris en compte. Soyez rigoureux/rigoureuses et précis(es) : notez scrupuleusement les paramètres
utilisés, ne vous embrouillez pas entre les différentes simulations, sauvegardez les fichiers images et les fichiers de données,
…
Nous donnons ici quelques idées de questions à se poser :
— Si le routeur est situé obligatoirement dans le salon et le long du mur « gauche », peut-on lui trouver une place
pour pouvoir toujours être connecté au Wi-Fi tout en étant aux toilettes ? (nous pouvons par exemple supposer
qu’en deça de 25% de la puissance max, la réception n’est pas suffisante pour glander sur youtube).
— À défaut de trouver l’optimal, quel serait un bon emplacement pour le routeur, pour obtenir du réseau Wi-Fi
partout ?
— Sachant que l’indice 𝑛 du béton est de l’ordre de 7. Que se passe-t-il si un des murs de l’appartement est en béton
(mur porteur) ? Pouvons-nous obtenir une cartographie des différences entre la solution avec mur en placo et avec
un mur porteur ?
— Influence du pas de maillage sur la solution ?

8.4. Étude et examen oral 83

Maillage et Éléments Finis

— Que se passe-t-il si les murs deviennent comme de l’air (𝑛(Ωmur ) = 1)) ? Pouvons nous obtenir une cartographie
des différences entre la solution avec mur en placo et mur « en air » ?
— Et si les murs étaient en béton (𝑛 = 7) ?
— Et si c’était votre appartement ?
— …
Réfléchissez également à d’autres questions plus théoriques et essayez d’apporter des réponses, comme par exemple
— Quelle sont les limites du modèle ?
— Que et comment pourrions nous améliorer notre modèle et notre code si nous disposions d’une machine de puis-
sance infinie ?
— Quelles améliorations, du point de vue géométrique, pouvons-nous apporter à notre appartement ?
Plus que l’obtention de réponse « parfaite », ce sera votre capacité et votre volonté de recherche qui seront pris en compte.

8.5 Résultat

8.5.1 Exemple

La figure 3 illustre un exemple de ce que l’on peut obtenir.

Fıg. 3 – Propagation d’une onde Wi-Fi dans un appartement. Après avoir traversé 2 murs, l’onde Wi-Fi semble très
amortie. Sous le résultat est affiché le plan de l’appartement et la position du routeur (petit disque à gauche)

8.5.2 Vous voulez tester ?

— Téléchargez le bundle Onelab. Il contient GMSH et GetDP (un solveur éléments finis)
— Téléchargez le code, soit directement soit via Git :

git clone [Link] wifi

— Dans le dossier et dans un terminal, lancer

gmsh [Link]

84 Chapitre 8. 2017 - 2018 : Wi-Fi

 Maillage et Éléments Finis

— Vous pouvez modifier un peu la géométrie et la fréquence de l’onde, mise à 1GHZ. Attention, cette simulation est
très gourmande : testez d’aborg avec 1GHz avant de lancer la simulation pour 2.5GHz (au risque de faire crasher
votre ordinateur) !

8.5. Résultat 85
Maillage et Éléments Finis

86 Chapitre 8. 2017 - 2018 : Wi-Fi

 CHAPITRE 9

2020 - 2021

9.1 Problème

Soit le domaine suivant modélisant un appartement :

Les murs sont supposés parfaitement isolant (nous pouvons toujours rêver !) , ce qui explique la condition de Neumann
homogène que nous imposons. Nous remarquons que le bord de Ω est séparé en plusieurs parties : les radiateurs (ΓRad ),
les fenêtres (ΓFen ) et les murs (ΓMur ). Nous cherchons à calculer la température 𝑢 dans la pièce, qui vérifie le système
suivant
⎧
⎪
⎪ −Δ𝑢 = 0 (Ω)
𝑢 = 𝑇𝑐 (ΓRad )
⎨
(1)
⎪
⎪ 𝑢 = 𝑇𝑓 (ΓFen )
𝜕n 𝑢 = 0 (ΓMur )
⎩

Les paramètres sont les suivants :

— La longueur 𝐿 ≃ 10
— la largeur ℓ ≃ 10
— l’épaisseur des murs 𝑑 ≃ 0.5
— la longueur d’une fenêtre est d’environ 1
— la longueur d’un radiateur est d’environ 1
— Les températures 𝑇𝑐 = 25 et 𝑇𝑓 = −10 sont les températures respectivement des radiateurs (ça chauffe) et de
dehors (ça caille)
Tous ces paramètres peuvent être librement modifiés par vous même et les paramètres de la géométrie doivent même être
choisis par vous !

87
Maillage et Éléments Finis

Fıg. 1 – Domaine Ω et ses bords : un petit appartement tout mignon

88 Chapitre 9. 2020 - 2021

 Maillage et Éléments Finis

9.2 Travail demandé

Nous souhaitons résoudre ce problème à l’aide de la méthode des éléments finis P1 − Lagrange.
1. Construisez la géométrie. Il ne s’agit pas de reproduire exactement l’appartement décrit plus haut mais de
construire un appartement : libre à vous d’ajouter des pièces, fenêtres, des radiateurs ou un poster de Justin Bieber.
2. Programmez un code éléments finis P1 qui résolve le problème (1)

9.3 Consignes

1. Vous pouvez réaliser ce projet en binôme

2. Vous devez rendre ce projet sous la forme d’un dépôt git : envoyez moi l’URL du dépôt uniquement, pas de
fichier zip !
3. Date limite de rendu : 28 février à 23h59. Pas de blague, tout projet rendu après le 28 février ne sera pas
noté.

En outre, votre dépôt doit contenir :

1. Un script qui résout le problème (1) et affiche la solution.
2. Un court fichier [Link] facilitant sa compréhension, répondant notamment aux questions : « comment
lance-t-on vos programmes ? » et « que doit-on obtenir ? » (exemple : « Exécutez “[Link]” et vous devez obtenir
la même image que “[Link]” qui résout le problème méga compliqué »)
3. Ajoutez à votre dépôt une image de la solution que vous avez obtenue (en PNG (pas très grosse svp !) et surtout
pas en JPG). Vous pouvez même afficher l’image dans le fichier [Link].

Remark 2.19
Pour faire simple, je dois pouvoir télécharger votre dépôt, lancer un fichier, et voir la solution, le tout sans avoir à réfléchir
de mon côté :-)

9.2. Travail demandé 89

Maillage et Éléments Finis

90 Chapitre 9. 2020 - 2021

Quatrième partie

Download

91
 Proof Index

corollary-11 lemma-2
corollary-11 (lecture/weak-solution/sobolev-spaces), lemma-2 (lecture/elements-finis-triangulaires/contributions-elementaires),
21 38

corollary-2 lemma-3
corollary-2 (lecture/weak-solution/sobolev-spaces), lemma-3 (lecture/elements-finis-triangulaires/espace-p1),
18 30

definition-0 lemma-4
definition-0 (lecture/weak-solution/rappels), 11 lemma-4 (lecture/elements-finis-triangulaires/contributions-elementaires),
40
definition-1
definition-1 (lecture/weak-solution/rappels), 11 lemma-5
lemma-5 (lecture/elements-finis-triangulaires/contributions-elementaires),
definition-2 41
definition-2 (lecture/weak-solution/rappels), 11
lemma-6
definition-3 lemma-6 (lecture/elements-finis-triangulaires/contributions-elementaires),
42
definition-3 (lecture/weak-solution/sobolev-spaces),
19 lemma-7
definition-4 lemma-7 (lecture/elements-finis-triangulaires/espace-p1),
32
definition-4 (lecture/weak-solution/sobolev-spaces),
19 proposition-1
definition-5 proposition-1 (lec-
ture/elements-finis-triangulaires/espace-p1),
definition-5 (lecture/weak-solution/rappels), 12
29
definition-6 proposition-12
definition-6 (lecture/weak-solution/rappels), 12 proposition-12 (lec-
ture/weak-solution/sobolev-spaces), 21
lemma-0
lemma-0 (lecture/elements-finis-triangulaires/galerkin), proposition-2
26 proposition-2 (lec-
ture/boundary-conditions/dirichlet), 50
lemma-1
lemma-1 (lecture/elements-finis-triangulaires/galerkin), proposition-4
26
proposition-4 (lecture/weak-solution/rappels), 12

93
Maillage et Éléments Finis

proposition-5 theorem-7
proposition-5 (lec- theorem-7 (lecture/weak-solution/rappels), 12
ture/weak-solution/sobolev-spaces), 19

proposition-6
proposition-6 (lec-
ture/weak-solution/sobolev-spaces), 19

remark-0
remark-0 (projet/2020-2021/index), 89

remark-1
remark-1 (numeric/gmsh/api), 64

remark-10
remark-10 (lecture/boundary-conditions/dirichlet), 53

remark-2
remark-2 (lecture/weak-solution/weak-formulation), ??

remark-3
remark-3 (lecture/weak-solution/weak-formulation), ??

remark-4
remark-4 (lecture/weak-solution/weak-formulation), ??

remark-5
remark-5 (projet/2017-2018/index), 80

remark-6
remark-6 (lecture/boundary-conditions/dirichlet), 51

remark-7
remark-7 (projet/2017-2018/index), 83

remark-8
remark-8 (lecture/weak-solution/sobolev-spaces), 20

remark-9
remark-9 (projet/2017-2018/index), 83

theorem-0
theorem-0 (lecture/weak-solution/lax-milgram), 16

theorem-1
theorem-1 (lecture/weak-solution/sobolev-spaces), 18

theorem-10
theorem-10 (lecture/weak-solution/sobolev-spaces), 21

theorem-4
theorem-4 (projet/2017-2018/index), 80

94 Proof Index