3.3. Assemblage des Matrices#

Nous devons maintenant calculer effectivement les coefficients \(A_{I,J}\) de la matrice \(A\) et \(\Bh_{I}\) du vecteur \(B\). Nous nous intéressons pour l’instant uniquement à la matrice \(A\).

3.3.1. Algorithme « brut-force »#

Prenons deux indices de sommets \(I\) et \(J\) et rappelons la valeur du coefficient \(A_{I,J}\) :

\[A_{I,J} = a(\mphi_J, \mphi_I) = \int_{\Omega}\nabla \mphi_J \cdot\nabla \mphi_I+ c\int_{\Omega}\mphi_J\mphi_I\]

Chaque intégrale sur \(\Omega\) peut être décomposée comme une somme sur les triangles \(\tri_p\) :

\[\begin{split}\begin{aligned} A_{I,J} &= \sum_{p=1}^{\Nt} \int_{\tri_p}\nabla \mphi_J \cdot\nabla \mphi_I+ c\sum_{p=0}^{\Nt-1} \int_{\tri_p}\mphi_J\mphi_I\\ \Bh_{I} &= \sum_{p=1}^{\Nt}\int_{\tri_p}f(x)\mphi_I(x)\diff x. \end{aligned}\end{split}\]

Soit deux sommets \(\vertice_I\) et \(\vertice_J\) n’appartenant pas un même triangle, alors \(\supp(\mphi_I)\cap\supp(\mphi_J) =\emptyset\). Autrement dit, \(\mphi_I\mphi_J\) est toujours nul et donc le coefficient \(A_{I,J}\) est nul ! Vue autrement, si deux sommets \(\vertice_I\) et \(\vertice_J\) ne sont pas connectés par une arête, alors \(A_{I,J=0}\).

Les coefficients de \(A\) 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 \(A\) 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 \(\vertice_I\) et \(\vertice_J\). Cependant, cet algorithme a tout de même un coût en \(\grandO{\Ns^2}\) ce qui est trop important pour être utilisable en pratique.

Code source 3.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

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 \(A\). Reprenons l’expression du coefficient \(A_{I,J}\):

\[A_{I,J} = \sum_{p=1}^{\Nt} \underbrace{\int_{\tri_p}\nabla \mphi_J \cdot\nabla \mphi_I}_{\text{Contrib. élémentaire}}+ c\sum_{p=0}^{\Nt-1} \underbrace{\int_{\tri_p}\mphi_J\mphi_I}_{\text{Contrib. élémentaire}}\]

Introduisons \(a_p(\cdot,\cdot)\) la famille de forme bilinéaire suivante, pour \(p=1,\ldots,\Nt\) :

\[a_p(\mphi_J,\mphi_I) = \int_{\tri_p}\nabla \mphi_J(\xx) \cdot\nabla \mphi_I(\xx)\diff \xx +c\int_{\tri_p}\mphi_J(\xx)\mphi_I(\xx)\diff \xx\]

Ensuite, nous réécrivons la matrice \(A\) sous la forme suivante

\[A = \sum_{I=1}^{\Ns}\sum_{j=0}^{\Ns-1}a(\mphi_J,\mphi_I) \ee_I^T\ee_J,\]

où \(\ee_I\) est le vecteur de la base canonique de \(\Rb^{\Ns}\). Nous avons alors

(3.3)#\[\begin{split}\begin{aligned} A &= \sum_{I=1}^{\Ns}\sum_{J=1}^{\Ns}a(\mphi_J,\mphi_I) \ee_I^T\ee_J\\ &= \sum_{I=1}^{\Ns}\sum_{J=1}^{\Ns}\sum_{p=1}^{\Nt}a_{p}(\mphi_J,\mphi_I) \ee_I^T\ee_J\\ &= \sum_{p=1}^{\Nt}\sum_{I=1}^{\Ns}\sum_{J=1}^{\Ns}a_{p}(\mphi_J,\mphi_I) \ee_I^T\ee_J\\ \end{aligned}\end{split}\]

Nous remarquons maintenant que \(a_{p}(\mphi_J,\mphi_I)\) est nul dès lors que \(\vertice_I\) ou \(\vertice_J\) ne sont pas des sommets de \(\tri_p\) (car \(\mphi_I\mphi_J = 0\) sur \(\tri_p\)). Finalement, la somme sur tous les sommets du maillage se réduit à une somme sur les 3 sommets du triangle \(\tri_p\) 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 \(\locToGlob\) (Local To Global)permettant de basculer du local vers le global une fonction telle que, pour \(p=1,\ldots,\Nt\) et \(i=1,2,3\) :

\[\locToGlob(p,i) = I \iff \vertice_i^p = \vertice_I\]

Ainsi, pour un triangle \(\tri_p\), ses sommets sont numérotés \([\vertice_{1}^{p},\vertice_{2}^{p},\vertice_{3}^{p}]\) en numérotation locale ou \([\vertice_{\locToGlob(p,1)},\vertice_{\locToGlob(p,2)},\vertice_{\locToGlob(p,3)}]\) en numérotation globale, comme le montre la figure 3.4. Nous distinguerons la numérotation globale par des lettres capitales (\(I\), \(J\)) et la numérotation locale par des minuscules (\(i\), \(j\)). Nous introduisons aussi les fonctions de forme locales :

\[\mphi_i^p = \mphi_{\locToGlob(p,i)}|_{\tri_p}.\]

Remark 3.5

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

Cliquez sur un triangle pour faire apparaitre la numérotation locale des sommets du triangle. Recliquez dessus pour revenir en numérotation globale

Utilisons ces nouvelles notations dans l’équation (3.3), en ramenant la somme sur les sommets à uniquement les sommets du triangle considéré :

\[A = \sum_{p=1}^{\Nt}\sum_{i=1}^{3}\sum_{j=1}^{3}a_{p}(\mphi_j^p,\mphi_i^p) \ee_{\locToGlob(p,i)}^T\ee_{\locToGlob(p,j)}\]

L’algorithme d’assemblage est alors complet ! Une version pseudo-code est présenté par l’algorithme d’assemblage. Sa complexité est en \(\grandO{\Nt} \ll \grandO{\Ns^2}\). Comme le premier algorithme, il possède en plus l’avantage d’être parallélisable.

Code source 3.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 \(a_p(\mphi_j^p,\mphi_i^p)\) et \(\ell_p(\mphi_i^p)\). De plus, il manque encore les conditions de Dirichlet.

Assemblage des Matrices

Contenu

3.3. Assemblage des Matrices#

3.3.1. Algorithme « brut-force »#

3.3.2. Algorithme d’assemblage#