Méthode de Galerkin et éléments finis

Méthode de Galerkin

Nous avons vu dans la section précédente que la forme faible de chaque problème considéré pouvait se mettre sous la forme :

$V$ étant de dimension infinie, un tel problème est particulièrement difficile à résoudre. La méthode de Galerkin consiste à résoudre la formulation faible dans un sous-espace $V_h \subset V$ de dimension finie $n_h$ possédant les mêmes propriétés, soit :

où $u_h$ est une solution approchée de la solution exacte $u$.

L’idée est de résoudre cette formulation à partir d’une base $(\varphi_1,\varphi_2,\dots,\varphi_{n_h})$ de $V_h$. En décomposant $u_h$ dans cette base, on a :

Le principe de la méthode de Galerkin est alors d’utiliser les fonctions de base $\varphi_i$ comme fonctions test dans la formulation faible. Ainsi :

Soit :

En définissant les vecteurs ${\bf U_h} = \left( u_i \right)_{1\leq i \leq n_h}$ et ${\bf B_h} = \big( L(\varphi_i) \big)_{1\leq i \leq n_h}$, et la matrice ${{\bf A_h} = \big( a(\varphi_i,\varphi_j) \big)_{1\leq i,j\leq n_h}}$, le problème peut se mettre sous la forme matricielle :

En résolvant ce système linéaire, on obtient alors les composantes $u_i$ de la solution approchée $u_h$ dans la base $(\varphi_i)_{1\leq i\leq n_h}$.

$\rightarrow$ La méthode des éléments finis permet d’appliquer ce principe de résolution.

Qu’est-ce qu’un élément fini ?

On appelle élément fini le triplet $(K,P_K,\Sigma_K)$, où :

• $K$ est un domaine de l’espace appelé élément géométrique (généralement de forme simple : un segment, un triangle, un quadrangle, un tétraèdre, un hexaèdre, une pyramide, un prisme) ;
• $P_K$ est un espace fonctionnel de dimension finie $n_K$, défini sur $K$ engendré par ce qu’on appelle les fonctions de base $p_i, 1 \leq i \leq n_K$ ;
• $\Sigma_K$ est un ensemble de $n_K$ degrés de liberté représentés par $n_K$ fonctionnelles (formes) linéaires $\sigma_i,~ 1 \leq i \leq n_K,$ définies sur l’espace $P_K$ et à valeurs dans $\mathbb{R}$.

De plus, il faut qu’il existe une unique fonction de $P_K$ prenant des valeurs données aux degrés de liberté de $\Sigma_K$ : cette dernière condition définit l’unisolvance de l’élément fini $(K,P_K,\Sigma_K)$.

Nous allons maintenant voir comment mettre en place la méthode basée sur ces éléments.

Maillage

Intéressons nous tout d’abords aux éléments géométriques $K$ obtenus par discrétisation du domaine d’étude $\Omega$, ce qu’on appelle le maillage. Le domaine (polyédrique) est alors partitionné avec le type d’entités élémentaires choisi (triangle, quadrangle, tétraèdre, hexaèdre, etc…).

Le maillage doit être conforme, c’est-à-dire :

• Les éléments doivent recouvrir la fermeture $\overline{\Omega}$ de $\Omega$.
• L’intersection de deux éléments est soit l’ensemble vide, soit un sommet, soit une arête partagée.

Nous pouvons illustrer cela sur un exemple :

Maillage conforme	Maillage non conforme

Dans la suite, nous utiliserons le logiciel libre de maillage par éléments finis Gmsh développé par C. Geuzaine et J-F. Remacle.

Génération d’un maillage

Pour pouvoir discrétiser un problème donné, nous devons tout d’abord définir ses aspects purement géométriques. C’est-à-dire créer toutes les entités élémentaires :

• les points, qui permettent de définir :
• les lignes (segments, arcs de cercles, voire courbes), qui permettent de définir :
• les surfaces (planes ou courbes), qui permettent de définir :
• les volumes.

Sur notre cas d’étude, on obtient par exemple :

Après avoir défini des tailles caractéristiques d’éléments et/ou discrétiser les lignes, on peut mailler notre géométrie. Ci-dessous, deux exemples :

Triangulaire	Quadrangulaire

Pour information, je vous donne en exemple les scripts permettant de générer ces maillages. Par habitude, j’utilise deux fichiers :

• data.geo : servant à définir les paramètres qui pourront être partagés avec le solveur,
• poisson.geo : contenant la création de la géométrie.

/* 
	Fichier data.geo contenant les parametres 
*/

// Parametres geometriques
Rd = 1;
Rs = 0.15*Rd;
xs = 0.4*Rd;  
ys = 0.4*Rd;
alphaN = Pi/6;

// Parametres maillage 
finesse = 2;
nec = 4*finesse;
lcd = Rd/nec;
lcs = Rs/nec*2;
neN = 2*nec;
typele = 1;  // <=1 : triangles
             // > 1 : quadrangles
             
// Numeros regions physiques
DISQUE = 1000;
SOURCE = 1001;
DIRICHLET = 2000;
NEUMANN = 2001;

/*
	Fichier poisson.geo contenant la creation de 
	la geometrie de notre petit probleme de Poisson
*/

// Lecture des parametres
Include "data.geo" ;

// Creation des points
Point(1) = {0,0,0,lcd};  // centre domaine
Point(2) = {Rd/Sqrt[2],Rd/Sqrt[2],0,lcd};
Point(3) = {-Rd*Cos[alphaN/2],Rd*Sin[alphaN/2],0,lcd};
Point(4) = {-Rd*Cos[alphaN/2],-Rd*Sin[alphaN/2],0,lcd};
Point(5) = {Rd/Sqrt[2],-Rd/Sqrt[2],0,lcd};
Point(6) = {xs,ys,0,lcs}; // centre Omega_s
Point(7) = {xs+Rs,ys,0,lcs};
Point(8) = {xs,ys+Rs,0,lcs};
Point(9) = {xs-Rs,ys,0,lcs};
Point(10) = {xs,ys-Rs,0,lcs};

// Creation des lignes
Circle(1) = {2,1,3}; 
Circle(2) = {3,1,4};
Circle(3) = {4,1,5}; 
Circle(4) = {5,1,2};
Circle(5) = {7,6,8}; 
Circle(6) = {8,6,9};
Circle(7) = {9,6,10}; 
Circle(8) = {10,6,7};

// Discretisation de la ligne Gamma_n
Transfinite Line {2}  = neN+1 Using Progression 1.0;

// Surfaces
Line Loop(1) = {1,2,3,4};
Line Loop(2) = {5,6,7,8};
Plane Surface(1) = {2};     // Omega_s
Plane Surface(2) = {1,-2};  // Omega

// Maillage en quadrangles (si voulu)
If (typele > 1)  
  Recombine Surface '*';
EndIf

Physical Surface(DISQUE) = {2};
Physical Surface(SOURCE) = {1};
Physical Line(DIRICHLET) = {1,3,4};
Physical Line(NEUMANN) = {2};

// Un peu de cosmetisme
Color Red { Surface{1}; }
Color White { Surface{2}; }
Color Blue { Line{2};}
Color Black { Line{1,3,4,5,6,7,8};}

En plus des fonctions basiques de création des entités élémentaires, Gmsh possède des fonctions plus évoluées : transformations (symétries, translations, rotations, dilatations), extrusions. De plus, il intègre depuis quelques temps un noyau supplémentaire basé sur le moteur CAO OpenCascade permettant de créer directement des objets 2D ou 3D, d’effectuer des opérations booléennes, ou de manipuler des splines.

Exemples

Ci-après, vous trouverez des exemples de maillages utilisés dans le cadre de travaux de mon équipe de recherche au laboratoire GREEN.

En 2D

Machine à courant continu	Ventouse magnétique

En 3D

Accouplements magnétiques (maillages structurés)	Machines synchrones à griffes (mix structuré / non-structuré)

$~$

Les exemples ci-dessus sont issus des thèses de :

• G. Devornique (ENSEM, promo 2013), pour l’alterno-démarreur synchrone à griffes ;
• B. Ristagno (ENSEM, promo 2016), pour la machine à courant continu ;
• D. Giraud (ENSEM, promo 2017), pour la machine synchrone axiale à griffes.

Autres exemples 3D

On peut également utiliser des maillages beaucoup plus complexes, comme par exemple les deux ci-dessous qui me permettent de calculer des courants induits dans des têtes humaines (par exemple pour des travailleurs opérant à proximité de sources basse fréquence comme des pinces à souder) :

John Doe	Jane Doe

$\boldsymbol{\rightarrow}$ Voir même un corps humain complet (oui, j’ai un peu galéré…) :

John en entier (avec tous les organes)

Fonctions de base sur éléments nodaux

Maintenant que nous savons créer le maillage (noté $\mathcal{M}_h$), nous allons nous intéresser à la définition des fonctions qui permettront de constituer la base de l’approximation discrète $V_h$ de $V$ (fonctions de base).

Éléments de référence

En fait, nous allons définir les fonctions de base sur ce qu’on appelle un élément de référence dont la forme géométrique est fixe et prédéfinie.

En 1D

En dimension 1, les éléments sont des segments et l’élément de référence est le segment $[-1,1]$.

En 2D

En dimension 2, nous pouvons donner le triangle et le quadrangle de référence :

• pour le triangle, c’est celui de sommets : ${\bf n_1}=(0,0)$, ${\bf n_2}=(1,0)$ et ${\bf n_3} = (0,1)$. On peut le paramétrer par : $$\Delta_r = \left\{ (u,v) \in \mathbb{R}^2 : 0 \le u \le 1, 0 \le v\le 1-u\right\}$$
• pour le quadrangle, c’est celui de sommets : ${\bf n_1} = (-1,-1)$, ${\bf n_2} = (1,-1)$, ${\bf n_3} = (1,1)$ et ${{\bf n_4} =(-1,1)}$. On peut le paramétrer par : $$\text{Q}_r = \left\{ (u,v) \in \mathbb{R}^2 : -1 \le u \le 1, -1 \le v \le 1\right\}$$

Leur représentation graphique est donnée ci-dessous.

En 3D

Nous nous contenterons de donner le tétraèdre et l’hexaèdre de référence :

• pour le tétraèdre, c’est celui de sommets : ${\bf n_1} = (0,0,0)$, ${\bf n_2} = (1,0,0)$, ${\bf n_3} = (0,1,0)$ et ${{\bf n_4} =(0,0,1)}$. On peut le paramétrer par : $$\text{T}_r = \left\{ (u,v,w) \in \mathbb{R}^3 : 0 \le u \le 1, 0 \le v\le 1-u, 0 \le w\le 1-u-v\right\}$$
• pour l’hexaèdre, c’est celui de sommets : ${\bf n_1} = (-1,-1,-1)$, ${\bf n_2} = (1,-1,-1)$, ${\bf n_3} = (1,1,-1)$, ${\bf n_4} = (-1,1,-1)$, ${\bf n_5} = (-1,-1,1)$, ${\bf n_6} = (1,-1,1)$, ${\bf n_7} = (1,1,1)$ et ${\bf n_8} = (-1,1,1)$. On peut le paramétrer par : $$\text{H}_r = \left\{ (u,v,w) \in \mathbb{R}^3 : -1 \le u \le 1, -1 \le v \le 1, -1 \le w \le 1\right\}$$

Et leur représentation graphique est donnée ci-dessous.

Fonctions de base sur éléments de référence

Sur ces élément de références, nous pouvons définir des fonctions linéaires valant 1 sur un sommet (un nœud) et 0 sur les autres. Dans ce cas, on fera une approximation linéaire par morceaux de notre fonction inconnue. Ces fonctions seront notées $p_i$ pour $i \in [\![1,n_K]\!]$.

Nous aurons donc autant de fonctions de base que de sommets $n_K$ de l’élément de référence.

Cas du triangle

Dans le cas du triangle de référence, ces fonctions sont données par la table ci-dessous :

Une représentation graphique est donnée par :

On peut également les tracer en surélevé :

$p_1$	$p_2$	$p_3$

Cas du quadrangle

Dans le cas du triangle de référence, ces fonctions sont données par la table ci-dessous :

Une représentation graphique est donnée par :

On peut également les tracer en surélevé :

$p_1$	$p_2$	$p_3$	$p_4$

Cas du tétraèdre

Dans le cas du tétraèdre de référence, ces fonctions sont données par la table ci-dessous :

Et une représentation graphique de $p_1$ est donnée ci-dessous :

Cas de l’hexaèdre

Dans le cas de l’hexaèdre de référence, ces fonctions sont données par la table ci-dessous :

Et une représentation graphique de $p_1$ est donnée ci-dessous :

Éléments quelconques

Considérons maintenant un élément réel $K$ de $\mathcal{M}_h$, et intéressons nous à la transformation permettant de passer de l’élément de référence associé $K_r$ à cet élément, ainsi que sa réciproque.

Transformations géométriques

Soit ${\bf x}$ un point quelconque de $K$ de coordonnées $(x,y,z)$. Une paramétrisation possible de ${\bf x}$ est :

où ${\bf u} = (u,v,w) \in K_r$ permet de faire correspondre tout point $(x,y,z)$ de $K$ avec un point $(u,v,w)$ de $K_r$ de façon biunivoque (bijective). La transformation inverse est alors notée :

Une représentation schématique est donnée ci-dessous :

Les formules de différentiation des fonctions composées nous donnent alors :

La matrice ${\bf J}$ est appelée matrice jacobienne de transformation. Les relations ci-dessus peuvent s’écrire sous forme plus compacte :

On rappelle au passage la formule de l’inverse :

En pratique, on ne la calcule pas de cette façon. En effet, il est intéressant de remarquer que les colonnes de ${\bf J^{-1}}$ sont les gradients selon les coordonnées réelles ${\bf x}$ des fonctions $u({\bf x}),$ $v({\bf x}),$ et $w({\bf x})$. On les note : ${\bf grad_x}\,u({\bf x}),$ ${\bf grad_x}\,v({\bf x})$ et ${\bf grad_x}\,w({\bf x})$.

Nous pouvons calculer ces termes à partir des lignes de la matrice Jacobienne ${\bf J}$ qui sont ${}^{\operatorname t}\partial_u\,{\bf x},$ ${}^{\operatorname t}\partial_v\,{\bf x}$ et ${}^{\operatorname t}\partial_w\,{\bf x},$ et qui sont plus faciles à calculer.

Rappelons les deux propriétés de la comatrice suivantes :

Ainsi, on a :

On peut faire de même pour les deux autres colonnes, au final, on obtient :

Le déterminant de la matrice jacobienne de transformation ($\det\,{\bf J},$ appelé Jacobien) peut se calculer par :

Ce Jacobien est très important car il intervient directement dans les formules d’intégration sur l’élément $K$ qui sont (par la formule de changement de variable) :

Ci-dessous, vous trouverez une petite application permettant de calculer le Jacobien d’un triangle quelconque (avec l’aimable autorisation de son auteur B. Thierry) :

Éléments isoparamétriques

Un élément isoparamétrique est un élément fini dont les fonctions de base nodales (qui permettent l’interpolation des champs scalaires) sont aussi utilisées pour paramétrer l’élément géométrique associé (fonctions d’interpolation géométrique) [Dula94].

Les fonctions de base vues ci-dessus sont définies par morceaux, sur chacun des éléments géométriques qui remplissent le domaine étudié, et permettent alors d’assurer leur continuité aux interfaces entre éléments. Ainsi, il n’y aura pas de discontinuité ni des champs scalaires interpolés, ni des coordonnées après transformation des éléments de référence vers les éléments réels. De telles fonctions de base sont dites conformes.

Considérons un élément fini nodal $(K,P_K,\Sigma_K)$. Si $n_K$ est l’ensemble des nœuds de $K$, de coordonnées ${\bf x_i},~i \in \mathbb{N},$ et si les $p_i({\bf u}),~i \in [\![1,n_K]\!],$ sont ses fonctions de base exprimées dans les coordonnées ${\bf u}$ de l’élément de référence $K_r$ associé à $K$, alors la paramétrisation de $K$ $\left({\bf x} = {\bf x}({\bf u})\right)$ est donnée par :

et cet élément est isoparamétrique.

Fonctions de base globales

À l’aide de la transformation géométrique précédente et des fonctions de bases $p_i,~i \in [\![1, n_K]\!]$ définies sur l’élément de référence, nous pouvons définir les fonctions de base globales $\varphi_i,~i \in [\![1, n_h]\!]$ définies sur le maillage complet $\mathcal{M}_h$. Il y en aura autant que de nœuds et chacune sera liée à un nœud particulier. Elle seront non-nulles uniquement sur les éléments contenant le nœud considéré. Cet ensemble d’éléments est appelé le support de la fonction de base.

Finalement, d’après ce qui précède, la fonction de base $\varphi_k$ liée au nœud $k \in [\![1, n_h]\!]$ sera nulle sur tous les nœuds du maillage sauf le nœud $k$ où elle prendra la valeur $1$. Soit, en notant ${\bf x_j}$ les coordonnées du nœud $j$ :

Une illustration interactive est donnée par l’application suivante (toujours avec l’aimable autorisation de son auteur) :

La famille $(\varphi_i)_{0\le i \le n_h}$ ainsi construite forme une base de $V_h$ qui est donc de dimension $n_h$ égale aux nombre de nœuds du maillage $\mathcal{M}_h$.

Soit $u_h$ une fonction de $V_h$. En la décomposant dans la base obtenue :

Les coefficients $u_i$ sont en fait les valeurs de la fonctions au nœuds. En effet en appliquant en ${\bf x_j}$ :

Finalement, une fonction $u_h$ de $V_h$ se décompose donc dans la base $(\varphi_i)_{0\le i \le n_h}$ selon :

et elle est entièrement caractérisée par ses $n_h$ valeurs aux nœuds.

Remarques complémentaires importantes

1. Dans le cas présenté, les fonctions de base de l’élément de référence servent aussi à paramétrer les éléments réels (fonctions d’interpolation géométrique). Cela correspond au cas particulier des éléments isoparamétriques, mais ces deux types de fonctions sont indépendants l’un de l’autre :

   • On pourrait tout à fait utiliser des fonctions de base d’ordre 2 (interpolation quadratique) ou d’ordre 3 (interpolation cubique) avec la transformation géométrique précédente (et les éléments rectilignes précédents). Un exemple fonctionnant sous Onelab vous sera fourni page suivante.
   • On pourrait utiliser des fonctions d’interpolation géométrique d’ordre 2 pour définir des éléments courbes. Le cas isoparamétrique (ou les deux types de fonctions sont d’ordre 2) correspond aux éléments dits $\mathbf{P_2}$ (ou $\mathbf{Q_2}$ si quadrangles).

2. Les éléments présentés ci-dessus conviennent parfaitement pour représenter des champs scalaires, nous verrons plus tard lesquels utiliser pour les champs vectoriels.

3. Plus que les fonctions de bases globales, ce sont les fonctions de base sur les éléments de référence et la matrice Jacobienne qui sont utiles pour créer le système à résoudre, c’est ce que nous allons voir de suite.

Assemblage et résolution

Assemblage du système à résoudre

On rappelle l’approximation de la formulation faible à résoudre obtenue par la méthode de Galerkin :

Avec :

Et :

Cette formulation est équivalente à la forme matricielle :

avec :

Calcul de ${\bf A_h}$ (matrice de rigidité ou de raideur)

Nous pourrions calculer les $n_h^2$ termes $A_{ij}$ séparément en effectuant une double boucle sur les nœuds. Cette façon de procéder n’est pas très pertinente car de nombreux termes sont nuls : typiquement, lorsque le nœud $i$ n’appartient pas à un élément du support de $\varphi_j$.

Une méthode plus efficace consiste à boucler sur tous les éléments de $\mathcal{M}_h$ et pour chaque élément $K_q$ calculer sa contribution $a_q(\varphi_k,\varphi_l),$ $(k,l) \in [\![1,n_{K_q}]\!]$ à la matrice globale ${\bf A_h}$ :

D’après notre formule de changement de base :

En utilisant les relations entre les gradients et la matrice jacobienne de transformation géométrique :

Finalement :

Nous avons vu comment calculer ${\bf x}({\bf u})$, $\text{det}\,{\bf J_q}$ et ${\bf J_q^{-1}}$ à la page précédente, et les gradients des fonctions de base se déterminent facilement à partir de leurs expressions.

Ainsi, pour calculer les $n_{K_q}$ contributions de l’élément $K_q$, il nous reste plus qu’à savoir comment déterminer l’intégrale sur l’élément de référence $K_r$. Nous le verrons ci-après.

Calcul de ${\bf B_h}$ (second membre)

Intéressons nous maintenant au membre de droite en l’absence de condition de Neumann particulière.

Cas sans condition de Neumann hétérogène

Nous allons procéder de la même manière que pour ${\bf A_h}$ en calculant les contributions $L_q(\varphi_k)$ de chaque élément $K_q$ à ${\bf B_h}$ :

Finalement :

Cas avec condition de Neumann hétérogène

Dans le cas d’une condition de Neumann non nulle, nous procédons comme ci-dessus mais en plus nous devons également calculer les contributions $L_q^n(\varphi_k)$ des éléments de la frontière $\Gamma_n$. Ceux-ci sont d’une dimension moindre que les éléments constituant le domaine $\Omega$. Ce seront donc des triangles ou des quadrangles en 3D, ou des segments en 2D.

Pour chaque éléments $K_q^n$ de $\Gamma_n$, nous calculons ainsi :

où $K_r^n$ est l’élément de référence associé aux éléments du bord du domaine.

Intégration numérique

Pour achever le calcul des coefficients des deux matrices, nous devons calculer les intégrales des différents termes sur le (ou les) élément(s) de référence. Pour cela nous allons utiliser une méthode de quadrature ; et généralement, c’est une formule de quadrature de Gauss.

Ainsi, on aura :

où les ${\bf x_{g_m}}$ sont les coordonnées des $n_g$ points de Gauss choisis et $w_m$ les poids associés.

On notera que le calcul donne même la valeur exacte dans le cas de polynôme (à condition d’avoir le bon nombre de points de Gauss).

Dans la table ci-dessous, sont donnés quelques exemples :

On remarquera que dans tous les cas, la somme des poids est égale à la longueur, l’aire ou le volume de l’élément considéré.

Résolution

Après génération du système linéaire précédent, il ne reste plus qu’à résoudre. Il convient donc de choisir une méthode / un algorithme pouvant résoudre efficacement des systèmes de grande taille mais creux.

Dans le cadre de ce cours, nous utiliserons GetDP pour définir et résoudre nos problèmes. Ce dernier s’appuie sur la suite PETSc configurée pour appeler le solveur linéaire direct multithreadé MUMPS (avec factorisation LU).

Après résolution, nous avons accès aux valeurs des degrés de liberté (valeurs du champ scalaire recherché) en chaque nœud du maillage. Nous pouvons ensuite calculer ce champ ou ses dérivés à partir de ces valeurs pour les tracer ou calculer tout autre grandeur locale ou globale. C’est ce qu’on appelle le post-traitement (post-processing).

Application à notre exemple (résolution avec GetDP)

Nous allons mettre en application ce qui précède sur notre exemple fil rouge, en écrivant le script GetDP permettant de le résoudre. Ce dernier doit porter le même nom que le fichier *.geo avec l’extension *.pro, ici : poisson.pro. Le déroulement suit les différentes étapes décrites dans les sections ci-avant, soit :

1. Lien avec le maillage : Où nous chargeons les différents paramètres et groupons les éléments en fonction de leur région d’appartenance.

Include "data.geo" ;

Group{
  Domaine = Region[{DISQUE,SOURCE}];
  Source = Region[{SOURCE}];
  Dirichlet = Region[{DIRICHLET}];
  Neumann = Region[{NEUMANN}];
}

2. Définition des fonctions :

Function{
  beta[#{Source}] = 1.0;
  gamma[#{Neumann}] = 1.0;
}

3. Définition du Jacobien et de la méthode d’intégration :

Jacobian{
  { Name Jvol;
    Case {
      { Region All ; Jacobian Vol; }
    }
  }
}

Integration {
  { Name Integ;
    Case { 
      {Type Gauss;
	Case { 
	  { GeoElement Line ; NumberOfPoints 2; }
	  { GeoElement Triangle ; NumberOfPoints 6; }
	  { GeoElement Quadrangle ; NumberOfPoints 4; }
	}
      }
    }
  }
}

4. Définition des contraintes : permettant de définir les conditions aux limites dans les espaces fonctionnels

Constraint {	
  {Name ConditionDirichlet ;
    Case {      
      {Region Dirichlet ; Type Assign ; Value 0. ;}
    }
  }
}

5. Définition de l’espace fonctionnel :

FunctionSpace {
  { Name Hgrad ; Type Form0 ; 
    BasisFunction { 
      { Name sn ; NameOfCoef phi ; Function BF_Node ; 
	Support Region[{Domaine,Neumann}]  ; Entity NodesOf[All] ; }
    } 
    Constraint {
      { NameOfCoef phi  ; EntityType NodesOf ; NameOfConstraint ConditionDirichlet ; }
    }
  }
}

6. Définition de la formulation (faible) :

Formulation {
  { Name Poisson ; Type FemEquation ;
    Quantity {
      { Name u  ; Type Local  ; NameOfSpace Hgrad ; }
    }
   Equation {
      Integral { [ Dof{Grad u}  , {Grad u} ] ;
		In Domaine ; Jacobian Jvol ; Integration Integ ; }
      Integral { [  -beta[] , {u} ] ;
		In Source ; Jacobian Jvol ; Integration Integ ; }     
      Integral { [ -gamma[] , {u} ] ; 
      	In Neumann ; Jacobian Jvol ; Integration Integ ; }
    }
  }
}

7. Résolution :

Resolution {
  { Name analyse ;
    System {
      { Name Sys ; NameOfFormulation Poisson ; }
    }
    Operation {
      Generate[Sys] ; 
      Solve[Sys] ;
      SaveSolution[Sys] ;
    }
  }
}

8. Post-traitement : Calcul des grandeurs (PostProcessing) et tracé ou afffichage (PostOperation).

PostProcessing {
  { Name PostProc ; NameOfFormulation Poisson ;
    Quantity {	
      { Name u  ; Value { Term { [ {u} ] ; In Domaine ; Jacobian Jvol ; } } }
      { Name Grad_u  ; Value { Term { [ {Grad u} ] ; In Domaine ; Jacobian Jvol ; } } }
      { Name normu ; Value { Integral { [ Norm[{u}] ]  ; In Domaine ; 
	    Jacobian Jvol ; Integration Integ ;} } } 
    }
  }
}

PostOperation PostOp UsingPost PostProc {
  Print[ Grad_u, OnElementsOf Domaine, File "gradu.pos"] ;
  Print[ u, OnElementsOf Domaine, File "u.pos"] ;
  Print[ normu[Domaine], OnGlobal, File "normu.dat", 
    SendToServer "Modele/Resultats/1Norme (u)", Color "AliceBlue" ];
}

J’ai regroupé l’ensemble des fichiers permettant de résoudre notre exemple dans un petit modèle Onelab pouvant être utilisé directement dans Gmsh.
Vous pouvez même l’exécuter sur votre smartphone avec l’application Onelab (disponible sur les stores officiels) !
Edit (janv. 2026) : Pour ceux qui sont sous Android ou dérivés, l’app n’est malheureusement plus dans le Play Store. Vous pouvez néanmoins télécharger l’apk à partir de ce lien et l’installer directement (je procède comme ça sous GrapheneOS).

J’ai rajouté :

• une option pour choisir un maillage triangulaire ou quadrangulaire,
• un paramètre de contrôle du maillage (plus la valeur est grande, plus le maillage est fin),
• une option pour choisir l’ordre d’interpolation (ordre des fonctions de base) : linéaire, quadratique, ou cubique (dernier cas seulement pour les triangles).

Exemple fil rouge à télécharger (compatible PC et smartphone) :

Éléments finis de Whitney

Contexte

Nous venons de développer complètement le principe de la méthode des éléments finis pour des champs scalaires appartenant à $\text{H}^1(\Omega) \left(= \text{H}({\bf grad},\Omega)\right)$. Pour cela, nous avons utilisé des éléments de Lagrange bien adaptés pour ce type d’inconnues.

Cependant les différents champs constituant le champ électromagnétique sont également des grandeurs vectorielles pouvant appartenir à ${\bf H}({\bf rot},\Omega)$ ou ${\bf H}(\text{div},\Omega)$. Typiquement nous ne pouvons résoudre notre formulation rot-rot avec des éléments de Lagrange tels quels. Nous pourrions choisir de résoudre avec trois champs scalaires pour déterminer notre vecteur inconnu (un par composante), mais cette solution est loin d’être la plus pertinente en terme d’efficacité et d’aspects géométriques (respect de certaines conditions de continuité aux interfaces).

Parmi la grande variété d’éléments finis existants, une famille est très largement utilisée pour modéliser des problèmes électromagnétiques : les éléments finis de Whitney. En fait, ces derniers ont été spécifiquement créés dans ce but [Boss93]. Basés sur les formes différentielles du même nom, ils offrent un cadre idéal pour approximer nos champs en 3D car ils permettent de construire de manière séquentielle les espaces discrets auxquels ils appartiennent. Nous allons donc présenter succinctement ces éléments.

Éléments de Whitney

Contrairement aux éléments précédents, les éléments de Whitney sont initialement définis en 3D et s’appuient sur des degrés de libertés pouvant être reliés aux valeurs aux nœuds (éléments nodaux), mais aussi aux circulations le long des arêtes (éléments d’arête), aux flux au travers des facettes (éléments de facette), et aux intégrales sur les volumes (éléments de volume) des éléments du maillage. Compte tenu de ces différences de nature entre degrés de liberté, on les appelle parfois éléments finis mixtes.

Espaces discrets

Le principe de base est de construire successivement les espaces discrets $W^0$, $W^1$, $W^2$ et $W^3$ permettant d’approximer respectivement $\text{H}({\bf grad},\Omega)$, ${\bf H}({\bf rot},\Omega)$, ${\bf H}(\text{div},\Omega)$ et $L^2(\Omega)$, à partir d’une séquence de fonctions de base, liées respectivement aux nœuds $(\mathcal{N_h})$, aux arêtes $(\mathcal{A_h})$, aux facettes $(\mathcal{F_h})$ et aux volumes $(\mathcal{V_h})$ des éléments du maillage $\mathcal{M_h}$.

On peut synthétiser l’ensemble dans la table suivante :

En procédant ainsi, on respecte la séquence vue sur les espaces continus, à savoir :

Et l’ensemble nous permet également d’assurer les continuités aux interfaces :

• de valeur pour les $s_n$ ;
• de composante tangentielle pour les ${\bf s_a}$ ;
• de composante normale pour les ${\bf s_f}$.

Voyons maintenant comment définir ces fonctions de base (sur les éléments de références).

Fonctions de base

Fonctions de base nodales

Puisque $W_0$ est une approximation de $\text{H}({\bf grad},\Omega) = \text{H}^1$, nous pouvons réutiliser les fonctions définies dans les parties précédentes. Nous aurons donc pour tout nœud $n_i = \{i\}$ de $\mathcal{N_h}$ :

Tout ce qui précède reste donc valable. Les fonctions de base nodales des éléments de Whitney sont celles des éléments de Lagrange.

Fonctions de base d’arête

À une arête ${\bf a_{ij}} = \{i\rightarrow j\}$ de $\mathcal{A_h}$, on associe :

Où $\operatorname{N}(i,\overline{j})$ représente l’ensemble des nœuds de la facette $f$ contenant le nœud $i$ mais pas le nœud $j$.

Remarque : Dans le cas du tétraèdre, l’expression ci-dessus se ramène à :

À titre d’illustration, nous traçons un exemple pour les triangle, quadrangle, tétraèdre et hexaèdre de référence ci-dessous :

${\bf s_{a_{23}}}$ sur tétraèdre	${\bf s_{a_{23}}}$ sur hexaèdre

En considérant une facette fictive $f$ orthogonale au plan d’étude, on peut définir nos fonctions de base en 2D, et dans ce cas, on obtient :

${\bf s_{a_{23}}}$ sur triangle	${\bf s_{a_{23}}}$ sur quadrangle

Fonctions de base de facette

À une facette orientée ${\bf f_{ijk}} = \{i,j,k\} = \operatorname{N}({\bf f_{ijk}})$ (si tétrahèdre) ou ${\bf f_{ijkl}} = \{i,j,k,l\} = \operatorname{N}({\bf f_{ijkl}})$ (si hexaèdre) de $\mathcal{F_h}$, on associe :

où la liste des nœuds de la facette $\operatorname{N}({\bf f})$ est rendue circulaire, $\operatorname{N}(i,\overline{j})$ conserve sa définition précédente et $a$ est un coefficient tel que : $ a= 2$ si $\dim(\operatorname{N}({\bf f})) = 3$ et $a = 1$ si $\dim(\operatorname{N}({\bf f})) = 4$.

On représente graphiquement ces fonctions de facette en 2D (facette fictive perpendiculaire au plan d’étude) dans la table ci-dessous :

${\bf s_{f_{23\perp}}}$ sur triangle	${\bf s_{f_{23\perp}}}$ sur quadrangle

Fonctions de base de volume

À l’élément $K$, on fait correspondre la fonction de base de volume :

qui est donc constante sur l’élément.

Degrés de liberté

Avec les définitions précédentes on peut décomposer tout champ scalaire ou vectoriel selon son espace d’appartenance :

1. Soit $u \in \text{H}({\bf grad},\Omega)$, son approximation sur $\mathcal{M_h}$ est $u_h \in W^0$ tel que :

   $$u_h = \sum\limits_{n\in\mathcal{N_h}} \sigma_n(u_h)\,s_n = \sum\limits_{n\in\mathcal{N_h}} u_n\,s_n$$

   avec :

   $$u_n = u_h({\bf x_n})$$

2. Soit ${\bf u} \in \textbf{H}({\bf rot},\Omega)$, son approximation sur $\mathcal{M_h}$ est ${\bf u_h} \in W^1$ tel que :

   $${\bf u_h} = \sum\limits_{a\in\mathcal{A_h}} \sigma_a({\bf u_h})\,{\bf s_a} = \sum\limits_{a\in\mathcal{A_h}} {u_a}\,{\bf s_a}$$

   avec :

   $$u_a = \int_a {\bf u_h}\cdot{\bf dl}$$

3. Soit ${\bf u} \in \textbf{H}(\text{div},\Omega)$, son approximation sur $\mathcal{M_h}$ est ${\bf u_h} \in W^2$ tel que :

   $${\bf u_h} = \sum\limits_{f\in\mathcal{F_h}} \sigma_f({\bf u_h})\,{\bf s_f} = \sum\limits_{f\in\mathcal{F_h}} {u_f}\,{\bf s_f}$$

   avec :

   $$u_f = \iint_f {\bf u_h}\cdot{\bf n}~\text{d}S$$

4. Soit $u \in L^2(\Omega)$, son approximation sur $\mathcal{M_h}$ est $u_h \in W^3$ tel que:

   $$u_h = \sum\limits_{v\in\mathcal{V_h}} \sigma_v(u_h)\,s_v = \sum\limits_{v\in\mathcal{V_h}} u_v\,s_v$$

   avec :

   $$u_v = \iiint_v u_h~\text{d}V$$

Ainsi, en décomposant nos approximations dans les bons espaces, on peut appliquer la méthode de Galerkin en utilisant les fonctions de base ci-dessus comme fonctions tests. Par analogie avec les éléments nodaux de la section précédente, on se ramènera donc à la résolution d’un système linaire où le vecteur inconnu ${\bf U_h}$ correspondra aux différents degrés de liberté.

Application des éléments d’arête à la formulation rot-rot

Formulation

On rappelle la formulation faible continue à résoudre :

Qui devient en discret :

Où $W_0^1 \in W^1$ est obtenu en annulant les degrés de liberté sur les arêtes appartenant à la frontière $\Gamma_d$. On constate que cela implique bien notre condition de bord :

La façon de procéder est identique à celle basée sur les éléments nodaux. Nous pouvons calculer les contributions de chaque arête élément par élément et les ajouter à la matrice de raideur globale ${\bf A_h}$, ou au vecteur second membre ${\bf B_h}$. On peut finalement résoudre le système obtenu, mais seulement après avoir implanté une condition de jauge comme expliqué ci-après. En effet, une telle formulation sous forme brute conduit à une matrice ${\bf A_h}$ singulière.

Condition de jauge

Nous savons malheureusement qu’un problème général (en 3D) comme celui-ci n’a pas de solution unique (vient de l’équivalence avec la formulation forte) et que le champ ${\bf u_h}$ solution est défini à un champ de gradient près. Il est alors impossible de résoudre le système obtenu sans avoir ajouté de condition de jauge particulière. Nous allons donc voir comment procéder.

Jauge d’arbre

Dans le domaine continu, pour rendre ${\bf u}$ unique, on peut ajouter une condition du type :

Où ${\bf w}$ est un champ de vecteurs dont les lignes de champ ne sont pas fermées et sont telles qu’elles peuvent relier toute paire de points quelconque du domaine d’étude $\Omega$ (par exemple un champ de gradient divergeant).

Dans le domaine discret, nous allons construire un arbre. C’est un ensemble d’arêtes qui relient tous les nœuds du maillage entre eux sans former de chemin fermé. Les arêtes du maillage qui ne sont pas dans l’arbre constituent le co-arbre associé. Un moyen de jauger ${\bf u_h}$ est ainsi d’annuler les composantes sur les arrêtes de l’arbre (condition équivalente à celle continue ci-dessus). Afin de ne pas dénaturer les conditions imposées sur $\Gamma_d$, il convient de commencer à construire l’arbre à partir des arrêtes de cette frontière. Une représentation graphique est donnée ci-dessous :

En rajoutant la condition de jauge d’arbre dans l’espace fonctionnel auquel appartient ${\bf u_h}$ (que nous noterons $W_{0,\text{JA}}^1 \subset W_0^1$), nous pouvons alors résoudre notre formulation faible rot-rot qui s’exprime :

Cette façon de procéder est la plus efficace d’un point de vue numérique car elle permet de réduire la taille du système à résoudre (par annulation des degrés de liberté liés au arêtes de l’arbre). Cependant la solution obtenue perd son sens physique, seul son rotationnel pourra être interprété physiquement. Pour de plus amples compléments, je vous renvoie à [Dula94] .

Jauge de Coulomb

Une autre Jauge possible est d’annuler la divergence de notre inconnue (Jauge de Coulomb) :

Pour obtenir la forme faible de cette condition, considérons un champ scalaire $\mu$ quelconque de $\text{H}({\bf grad},\Omega)$. On a :

En utilisant la formule de Green en grad-div :

En appliquant le théorème de la divergence et en séparant les deux frontières :

Nous ne considérerons que le cas d’une condition de Neumann homogène sur $\Gamma_n$ $(\boldsymbol{\gamma} = 0)$, soit : $\left.{\bf rot\,u}\wedge{\bf n}\right|_{\Gamma_n} = 0$. On peut donc en déduire : $\left.{\bf u}\cdot{\bf n}\right|_{\Gamma_n} = 0$. En choisissant $\mu = 0$ sur $\Gamma_d$, on annule le second terme surfacique.

Ainsi, notre condition de Jauge est équivalente à :

Si nous ajoutons uniquement ce terme à la formulation faible initiale, nous nous trouverons face à une formulation avec une fonction inconnue et deux fonctions tests. Pour remédier à cela, nous devons introduire un multiplicateur de Lagrange $\lambda$. Ainsi la formulation variationnelle complète devient :

La transcription dans le domaine discret est immédiate :

Cette façon de procéder est plus lourde à résoudre : le système est beaucoup plus grand que le précédent car on a rajouté autant de degrés de liberté que de nœuds du maillage. Cependant, elle a le mérite de fournir une solution « plus physique » et plus facilement interprétable.