On rappelle l’approximation de la formulation faible à résoudre obtenue par la méthode de Galerkin :
$$\forall i \in [\![1,n_h]\!]~,~ \sum\limits_{j=1}^{n_h} a(\varphi_i,\varphi_j)\,u_j = L(\varphi_i)$$
Avec : $$ a(\varphi_i,\varphi_j) = \iiint_{\Omega} \alpha\,{\bf grad}\,\varphi_i \cdot {\bf grad}\,\varphi_j~\text{d}\Omega$$
Et : $$\displaystyle L(\varphi_i) = \underbrace{\iiint_{\Omega} \beta\,\varphi_i~\text{d}\Omega}_{\left(\beta,\varphi_i\right)_{\Omega}} + \underbrace{\iint_{\Gamma_n} \alpha\,\gamma\,\varphi_i~\text{d}\Gamma}_{\left<\alpha\,\gamma,\varphi_i\right>_{\Gamma_n}}$$
Cette formulation est équivalente à la forme matricielle :
$${\bf A_h}\,{\bf U_h} = {\bf B_h}$$
avec : $$\left\{\begin{aligned}&{\bf A_h} = \big(a(\varphi_i,\varphi_j)\big)_{1\le i,j\le n_h} = (A_{ij})_{1\le i,j\le n_h}\\[1ex] &{\bf B_h} = \big(L(\varphi_i)\big)_{1\le i\le n_h} = (B_i)_{1\le i\le n_h} \end{aligned}\right.$$
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}$ :
$$a_q(\varphi_k,\varphi_l) = \iiint_{K_q} \alpha({\bf x})\,{\bf grad}\,\varphi_k({\bf x}) \cdot {\bf grad}\,\varphi_l({\bf x})~\text{d}{\bf x} $$
D’après notre formule de changement de base :
$$a_q(\varphi_k,\varphi_l) = \iiint_{K_r} \alpha\left({\bf x}({\bf u})\right)\,{\bf grad}\,\varphi_k\left({\bf x}({\bf u})\right) \cdot {\bf grad}\,\varphi_l\left({\bf x}({\bf u})\right)\,\text{det}\,{\bf J_q}~\text{d}{\bf u} $$
En utilisant les relations entre les gradients et la matrice jacobienne de transformation géométrique :
$$a_q(\varphi_k,\varphi_l) = \iiint_{K_r} \alpha\left({\bf x}({\bf u})\right)~\left({\bf J_q^{-1}}\,{\bf grad}\,p_k({\bf u})\right) \cdot \left({\bf J_q^{-1}}\,{\bf grad}\,p_l({\bf u})\right)\,\text{det}\,{\bf J_q}~\text{d}{\bf u}$$
Finalement :
$$\boxed{a_q(\varphi_k,\varphi_l) = \text{det}\,{\bf J_q} \iiint_{K_r} \alpha\left({\bf x}({\bf u})\right)~{}^{\operatorname t}\left({\bf grad}\,p_k({\bf u})\right)\,{}^{\operatorname t}({\bf J_q^{-1}})\,{\bf J_q^{-1}}\,{\bf grad}\,p_l({\bf u})~\text{d}{\bf u}}$$
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.
Intéressons nous maintenant au membre de droite en l’absence de condition de Neumann particulière.
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}$ :
$$\begin{aligned}L_q(\varphi_k) &= \iiint_{K_q} \beta({\bf x})\,\varphi_k({\bf x})\,\text{d}{\bf x}\\ &= \iiint_{K_r} \beta\left({\bf x}({\bf u})\right)\,p_k({\bf u})\,\text{det}\,{\bf J_q}~\text{d}{\bf u}\end{aligned}$$
Finalement :
$$\boxed{L_q(\varphi_k) = \text{det}\,{\bf J_q} \iiint_{K_r} \beta\left({\bf x}({\bf u})\right)\,p_k({\bf u})~\text{d}{\bf u}}$$
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 :
$$L_q^n(\varphi_k) = \iint_{K_q^n} \alpha({\bf x})\,\gamma({\bf x})\,\varphi_k({\bf x})\,\text{d}{\bf x}$$
$$\boxed{L_q^n(\varphi_k) = \text{det}\,{\bf J_q^n} \iint_{K_r^n} \alpha\left({\bf x}({\bf u})\right)\,\gamma\left({\bf x}({\bf u})\right)\,p_k({\bf u})\,\text{d}{\bf u}}$$
où $K_r^n$ est l’élément de référence associé aux éléments du bord du domaine.
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 :
$$\iiint_{K_r} f({\bf u})\,\text{d}{\bf u} \approxeq \sum\limits_{m=1}^{n_g} \omega_m\,f({\bf x_{g_m}})$$
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 :
| Élément | $n_g$ | Points de Gauss | Poids | Exacte pour |
|---|---|---|---|---|
| Segment | $1$ | $x_g = 0$ | $\omega =2$ | degré 1 |
| $2$ | $x_{g_1} = -\frac{1}{\sqrt{3}}$, $x_{g_2} = \frac{1}{\sqrt{3}}$ | $\omega_1 = \omega_2 = 1$ | degré 2 | |
| Triangle | $1$ | ${\bf x_g} = (\frac{1}{3},\frac{1}{3})$ | $\omega=\frac{1}{2}$ | degré 1 |
| $3$ | ${\bf x_{g_1}} = (\frac{1}{6},\frac{1}{6})$, ${\bf x_{g_2}} = (\frac{2}{3},\frac{1}{6})$, ${\bf x_{g_3}} = (\frac{1}{6},\frac{2}{3})$ | $\omega_1 = \omega_2 = \omega_3 = \frac{1}{6}$ | degré 2 | |
| Quadrangle | $1$ | ${\bf x_g} = (0,0)$ | $\omega=4$ | degré 1 |
| $4$ | ${\bf x_{g_1}} = (\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}})$, ${\bf x_{g_2}} = (-\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}})$, ${\bf x_{g_3}} = (\frac{1}{\sqrt{3}},-\frac{1}{\sqrt{3}})$, ${\bf x_{g_4}} = (-\frac{1}{\sqrt{3}},-\frac{1}{\sqrt{3}})$ | $\omega_1 = \omega_2 = \omega_3 = \omega_4 = 1$ | degré 2 | |
| Tétraèdre | $1$ | ${\bf x_g} = (\frac{1}{4},\frac{1}{4},\frac{1}{4})$ | $\omega=\frac{1}{6}$ | degré 1 |
| $4$ | ${\bf x_{g_1}} = (\alpha,\alpha,\alpha)$, ${\bf x_{g_2}} = (\alpha,\alpha,\beta)$, ${\bf x_{g_3}} = (\alpha,\beta,\alpha)$, ${\bf x_{g_4}} = (\beta,\alpha,\alpha)$, avec : $\alpha = 0,138196601125$, $\beta = 0,585410196625$ | $\omega_1 = \omega_2 = \omega_3 = \omega_4 = \frac{1}{24}$ | degré 2 |
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é.
Contrairement à ce qu’on pourrait penser, le choix du nombre de points de Gauss n’est pas primordial car l’erreur commise par la méthode d’intégration est faible devant les autres approximations de la méthode. Les seconds choix de la table précédente sont donc tout à fait pertinents (on n’utilisera que des fonctions de bases linéaires, voire quadratiques).
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.
Pour les différentes méthodes de résolution (avec ou sans pré-conditionnement), je vous renvoie à l’EC d’« Analyse numérique », du tronc commun de S6.
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).
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 :
Include "data.geo" ;
Group{
Domaine = Region[{DISQUE,SOURCE}];
Source = Region[{SOURCE}];
Dirichlet = Region[{DIRICHLET}];
Neumann = Region[{NEUMANN}];
}
Function{
beta[#{Source}] = 1.0;
gamma[#{Neumann}] = 1.0;
}
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; }
}
}
}
}
}
Constraint {
{Name ConditionDirichlet ;
Case {
{Region Dirichlet ; Type Assign ; Value 0. ;}
}
}
}
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 ; }
}
}
}
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 ; }
}
}
}
Resolution {
{ Name analyse ;
System {
{ Name Sys ; NameOfFormulation Poisson ; }
}
Operation {
Generate[Sys] ;
Solve[Sys] ;
SaveSolution[Sys] ;
}
}
}
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 ou sur votre smartphone avec l’application Onelab (disponible sur les stores officiels).
J’ai rajouté :
Exemple fil rouge à télécharger (compatible PC et smartphone) :