3.4. Calcul des Matrices Élémentaires#

3.4.1. Matrices de Masse et de Rigidité#

La matrice \(A\) peut être décomposée en deux matrices : la masse et la rigidité :

\[A = D + c M,\]

  • \(M\) : la matrice de masse (ou de volume), de coefficient

    \[M_{I,J} = \int_{\Omega} \mphi_J\mphi_I.\]

  • \(D\) : la matrice de rigidité, de coefficient

    \[D_{I,J}= \int_{\omega}\nabla\mphi_J\nabla\mphi_I.\]

Remark 3.7

Dans la littérature, cette matrice est souvent notée \(K\), mais nous l’appelons \(D\) pour éviter toute confusion avec les triangles, nommés \(K\) également.

Remark 3.8

La matrice de masse \(M\) 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 » \(u=f\) (ou \(Id. u = f\)) et appliquer la méthode des éléments finis pour obenir la « formulation faible »

\[\forall \vh, \quad \int_{\Omega} U\vh = \int_{\Omega}f\vh,\]

qui aboutit au système linéaire suivant : \(MU = B\). L’opérateur Identité, appliqué à \(u\), est bien discrétisé en \(M\).

Les contributions élémentaires, c’est à dire les quantités \(a_p(\mphi_j^p,\mphi_i^p)\) et \(\ell_{p}(\mphi_i^p)\), peuvent elles aussi être décomposées en deux parties. Pour rappel, les sommets d’un triangle \(\tri_p\) seront notés \([\vertice_{0}^{p}, \vertice_{1}^{p},\vertice_{2}^{p}]\) et ordonnés dans le sens trigonométrique. Nous noterons \(\vertice_i^p=(x_i^p, y_i^p)\) un sommet de \(\tri_p\) et \(\mphi_i^p\) la fonction de forme locale associée. Nous notons \(\Me{p}\) et \(\De{p}\) les matrices de masse et de rigidité élémentaire du triangle \(\tri_p\), de coefficient respectif \((\Me{p})_{i,j}\) et \((\De{p})_{i,j}\) donnés par

\[\begin{split}\begin{aligned} (\Me{p})_{i,j} &= \int_{\tri_p}\mphi_j^p\mphi_i^p\\ (\De{p})_{i,j} &=\int_{\tri_p}\nabla\mphi_j^p\cdot\nabla\mphi_i^p. \end{aligned}\end{split}\]

3.4.2. Matrice de masse élémentaire#

Nous nous focalisons sur la matrice de masse, le principe est similaire pour la matrice \(D\) et est détaillé juste après.

Pour construire la matrice \(M\), nous avons vu qu’il était préférable de parcourir les triangles plutôt que les sommets, autrement dit, plutôt que de calculer \(M_{I,J}\) directement, mieux vaut calculer, pour tout triangle \(p\), la contribution élémentaire \((\Me{p})_{i,j}\) pour \(i,j = 1,2,3\), définie par :

(3.4)#\[(\Me{p})_{i,j}= \int_{\tri_p} \mphi_j^p(\xx)\ \mphi_i^p(\xx)\diff\xx.\]

Chaque contribution élémentaire \((\Me{p})_{i,j}\) est ensuite ajoutée à \(M_{I,J}\), avec \(I=\locToGlob(p,i)\) et \(J=\locToGlob(p,j)\).

3.4.2.1. Triangle de référence#

Pour calculer la quantité élémentaire (3.4), plaçons nous tout d’abord dans un triangle « simple » \(\trih\), appelé triangle de référence. Celui-ci est souvent choisi comme étant le triangle rectangle de sommets \(\verticeh_{1}=(0,0)\), \(\verticeh_{2}=(1,0)\) et \(\verticeh_{3}=(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 \((\xi,\eta)\) dit repère paramétrique.

Triangle de référence

Fig. 3.5 Triangle de référence \(\trih\) et son repère paramétrique \((\xi,\eta)\).#

Nous notons \(\mphih_i \in \Pb^1(\trih)\) les trois fonctions de forme associées aux sommets \(\verticeh_i\), pour \(i=1,1,3\), définies par \(\mphih_i(\verticeh_j) = \delta_{ij}\). Ces fonctions \(\mphih_i\) étant des polynômes de degré un, nous pouvons les calculer analytiquement :

\[\begin{split}\left\{ \begin{array}{l} \mphih_1(\xi,\eta) = 1-\xi-\eta\\ \mphih_2(\xi,\eta) = \xi\\ \mphih_3(\xi,\eta) = \eta\\ \end{array} \right.\end{split}\]

Lemma 3.5

Dans le triangle \(\trih\), la matrice de masse élémentaire \(\Meh = (\Meh_{i,j})_{1\leq i,j\leq 3}\) de coefficient

\[\Meh_{i,j} = \int_{\trih} \mphih_j\mphih_i \diff(\xi,\eta),\]

est donnée par

\[\begin{split}\Meh = \frac{1}{24}\left( \begin{array}{c c c} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{array} \right).\end{split}\]

Proof. Prenons tout d’abord le cas \(i=j=2\), soit \(\mphih_i(\xi,\eta) = \mphih_j(\xi,\eta) = \xi\). Dans ce cas :

\[\int_{\trih} \xi^2 \diff (\xi,\eta) = \int_0^1\int_0^{1-\xi} \xi^2 \diff\eta\diff\xi = \int_0^1(1-\xi)\xi^2\diff\xi = \left[\frac{\xi^3}{3} - \frac{\xi^4}{4}\right]_0^1=\frac{1}{3}-\frac{1}{4} = \frac{1}{12}.\]

Les calculs sont similaires pour \(i=1\) et \(i=3\). Prenons maintenant \(i\neq j\), par exemple \(i=3\) et \(j=2\) :

\[\int_{\trih} \xi\eta \diff (\xi,\eta) = \int_0^1\left(\int_0^{1-\xi} \eta \diff\eta\right)\xi\diff\xi = \frac{1}{2}\int_0^1(1-\xi)^2\xi\diff\xi = \frac{1}{2}\left[ \frac{1}{2} - \frac{2}{3} +\frac{1}{4}\right] =\frac{1}{24}.\]

Les calculs sont similaires pour les autres combinaisons.

3.4.2.2. Triangle quelconque#

Changement de coordonnées. Soit un triangle \(\tri_p\) du maillage et supposons que nous disposions d’une transformation bijective et linéaire \(\trihToTri{p}\) permetteant de transformer le triangle de référence \(\trih\) en \(\tri_p\) avec en plus \(\trihToTri{p}(\verticeh_i) = \vertice_i^p\) (conservation de l’ordre des sommets). Cette fonction \(\trihToTri{p}\) transforme les coordonnées paramétriques \((\xi,\eta)\) en coordonnées physiques \((x,y)\) avec \((x,y)=\trihToTri{p}(\xi,\eta)\in\tri_p\), et conserve « l’ordre des sommets ».

Transformation entre le triangle de référence et un triangle quelconque

Fig. 3.6 Transformation entre le triangle de référence \(\trih\) et un triangle quelconque \(\tri_p\).#

Nous avons \(\mphi_j^p(x,y) = \mphi_j^p(\trihToTri{p}(\xi,\eta))\) avec \(\mphi_j^p\circ\trihToTri{p}\in\Pb^1(\trih)\) et \(\mphi_j^p\circ\trihToTri{p}(\verticeh_i) = \delta_{ij}\), soit exactement les mêmes propriétés que les \(\mphih_i\). Par unicité, nous avons \(\mphi_j^p\circ\trihToTri{p} = \mphih_j\).

En notant \(\JK{p}\) la matrice Jacobienne de \(\trihToTri{p}\), alors la quantité \((\Me{p})_{i,j}\) peut alors s’écrire, par changement de variables :

\[(\Me{p})_{i,j} = \displaystyle\int_{\tri_p}\mphi_j^p(x,y)\mphi_i^p(x,y) \diff(x,y) =\displaystyle \abs{\det(\JK{p})}\underbrace{\int_{\trih}\mphih_j(\xi,\eta)\mphih_i(\xi,\eta)\diff(\xi,\eta)}_{\text{Déjà calculé !}}\]

Ainsi, pour calculer la matrice élémentaire d’un triangle \(\tri_p\) quelconque, nous n’avons besoin que du déterminant de la Jacobienne : \(\det(\JK{p})\).

Expression et Jacobienne de la transformation. La transformation que nous cherchons, \(\trihToTri{p}\), est linéaire et « conserve » les sommets et leur ordre. Pour obtenir son expression, nous construisons des fonctions d’interpolation géométrique, \((\psih_i)_{1\leq i \leq 3}\), linéaires sur \(\trih\) et telles que :

\[\forall i,j=1,2,3 \quad \psih_i(\verticeh_j) = \deltaij.\]

La transformation aura alors pour expression :

\[\begin{split}\begin{array}{r c c l} \trihToTri{p}\colon & \trih & \to & \tri_p\\ & (\xi,\eta) & \mapsto & \trihToTri{p}(\xi,\eta) = (x,y) = \psih_{1}(\xi,\eta) \vertice_{1}^{p} + \psih_{2}(\xi,\eta) \vertice_{2}^{p} + \psih_{3}(\xi,\eta) \vertice_{3}^{p}. \end{array}\end{split}\]

En d’autres termes, les fonctions d’interpolation géométrique \(\psih_i\) sont ici identiques aux fonctions de forme \(\mphih_i\) :

\[\begin{split}\left\{ \begin{array}{l} \psih_{1}(\xi,\eta) = 1 - \xi - \eta\\ \psih_{2}(\xi,\eta) = \xi\\ \psih_{3}(\xi,\eta) = \eta\\ \end{array} \right.\end{split}\]

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

\[\begin{split}\JK{p} = \left( \begin{array}{c c} \displaystyle\frac{\partial x}{\partial \xi} &\displaystyle \frac{\partial x}{\partial \eta} \\ \displaystyle\frac{\partial y}{\partial \xi} &\displaystyle \frac{\partial y}{\partial \eta} \end{array} \right) = \left( \begin{array}{c c} x_{2}^{p} - x_{1}^{p} & x_{3}^{p} - x_{1}^{p}\\ y_{2}^{p} - y_{1}^{p} & y_{3}^{p} - y_{1}^{p} \end{array} \right),\end{split}\]

et son déterminant vaut

\[\begin{split}\begin{aligned} \abs{\det(\JK{p})} &= \abs{(x_{2}^{p}-x_{1}^{p})(y_{3}^{p}-y_{1}^{p}) - (x_{3}^{p}-x_{1}^{p})(y_{2}^{p}-y_{1}^{p})}\\ &= 2|\tri_p| \neq 0, \end{aligned}\end{split}\]

ce qui implique que le déterminant est non nul puisque le triangle n’est pas dégénéré : la transformation \(\trihToTri{p}\) est bien inversible.

Remark 3.9

Quand \(\psih_i = \mphih_i\), nous parlons d’éléments finis isoparamétriques. Il convient de retenir que ce choix n’est pas obligatoire et les fonctions \(\psih_i\) et \(\mphih_i\) sont indépendantes. En particulier, pour obtenir des éléments courbes, les fonctions \(\psih_i\) pourraient être quadratiques par exemple.

Déplacez les sommets du triangle pour modifier la valeur du Jacobien. Quand il est négatif cela signifie que le triangle est "retourné" par rapport au triangle de référence

Expression finale de la matrice élémentaire.

Lemma 3.6

La matrice de masse élémentaire \(\Me{p} = ((\Me{p})_{i,j})_{0\leq i,j\leq 2}\) du triangle \(\tri_p\) a pour expression

\[\begin{split}\Me{p} = \frac{\abs{\tri_p}}{12} \left( \begin{array}{c c c} 2 & 1 & 1\\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right).\end{split}\]

Cliquez sur un triangle pour ajouter ses contributions élémentaires dans la matrice de masse. Recliquez dessus pour les soustraire.

3.4.3. Matrice de rigidité élémentaire#

Nous appliquons la même procédure pour la matrice de rigidité \(D\), autrement dit, nous calculons les matrices de rigidité élémentaire \(\De{p}\) définies par

\[(\De{p})_{i,j} = \int_{\tri_p}\nabla \mphi_j^p(x,y)\cdot \nabla\mphi_i^p(x,y)\diff(x,y).\]

3.4.3.1. Triangle de référence#

Lemma 3.7

Dans le triangle de référence \(\trih\), la matrice de rigidité élémentaire \(\widehat{D}= (\widehat{D}_{i,j})_{1\leq i,j\leq 3}\) de coefficient

\[\widehat{D}_{i,j} = \int_{\trih}\nabla \mphih_j(\xi,\eta)\cdot \nabla\mphih_i(\xi,\eta)\diff(\xi,\eta),\]

a pour expression

\[\begin{split}\widehat{D} = \frac{1}{2} \left( \begin{array}{l l c} 2 & -1 & -1 \\ -1 & 1 & 0 \\ -1 & 0 & 1 \end{array} \right)\end{split}\]

Proof. Les gradients des fonctions de forme \(\mphih_j\) sont donnés par :

\[\begin{split}\nabla_{\xi,\eta}\mphih_0 = \begin{pmatrix} -1\\ -1 \end{pmatrix} , \quad \nabla_{\xi,\eta}\mphih_1 = \begin{pmatrix} 1\\ 0 \end{pmatrix}, \quad \nabla_{\xi,\eta}\mphih_2 = \begin{pmatrix} 0\\ 1 \end{pmatrix}.\end{split}\]

La matrice étant symétrique, nous pouvons limiter les calculs à la partie triangulaire supérieure :

\[\begin{split}\begin{aligned} \widehat{D}_{1,1} &= \int_{\trih}\nabla\mphih_1\cdot\nabla\mphih_1 \diff (\xi,\eta) = \int_{\trih} (-1,-1)\begin{pmatrix}-1\\ -1\end{pmatrix}\diff (\xi,\eta) = 2 \int_{\trih} \diff(\xi,\eta) &&= 1\\ \widehat{D}_{2,2} &= \int_{\trih}\nabla\mphih_2\cdot\nabla\mphih_2 \diff (\xi,\eta) = \int_{\trih} (1,0)\begin{pmatrix}1\\ 0\end{pmatrix} \diff (\xi,\eta) = \int_{\trih} \diff(\xi,\eta) &&= \frac{1}{2} =\widehat{D}_{3,3}\\ \widehat{D}_{1,2} &= \int_{\trih}\nabla\mphih_1\cdot\nabla\mphih_2 \diff (\xi,\eta) = \int_{\trih} (-1,-1)\begin{pmatrix}1\\ 0\end{pmatrix} \diff (\xi,\eta) = -\int_{\trih} \diff(\xi,\eta) &&= -\frac{1}{2}\\ \widehat{D}_{1,3} &= \int_{\trih}\nabla\mphih_1\cdot\nabla\mphih_3 \diff (\xi,\eta) = \int_{\trih} (-1,-1)\begin{pmatrix}0\\ 1\end{pmatrix} \diff (\xi,\eta) = -\int_{\trih} \diff(\xi,\eta)&& = -\frac{1}{2}\\ \widehat{D}_{2,3} &= \int_{\trih}\nabla\mphih_2\cdot\nabla\mphih_3 \diff (\xi,\eta) = \int_{\trih} (1,0)\begin{pmatrix}0\\ 1\end{pmatrix} \diff (\xi,\eta) &&= 0. \end{aligned}\end{split}\]

3.4.3.2. Triangle quelconque#

Pour calculer les dérivées partielles selon \(x\) et \(y\) de \(\mphih_j\), nous utilisons la dérivée de fonction composée :

\[\begin{split}\begin{pmatrix} \displaystyle \frac{\partial \mphi_j^p}{\partial x}\\[0.2cm] \displaystyle \frac{\partial \mphi_j^p}{\partial y} \end{pmatrix} = \begin{pmatrix} \displaystyle \frac{\partial \xi}{\partial x} & \displaystyle \frac{\partial \eta}{\partial x}\\[0.2cm] \displaystyle \frac{\partial \xi}{\partial y} & \displaystyle \frac{\partial \eta}{\partial y} \end{pmatrix} \begin{pmatrix} \displaystyle \frac{\partial \mphih_j}{\partial \xi}\\[0.2cm] \displaystyle \frac{\partial \mphih_j}{\partial \eta} \end{pmatrix}\end{split}\]

En notant \(\BK{p}\) la matrice de passage, nous avons

\[\nabla_{x,y}\mphi_j^p(x,y) = \BK{p}\nabla_{\xi,\eta}\mphih_j(\xi,\eta).\]

L’opération « inverse » nous donne :

\[\begin{split}\begin{pmatrix} \displaystyle \frac{\partial \mphih_j}{\partial \xi}\\[0.2cm] \displaystyle \frac{\partial \mphih_j}{\partial \eta} \end{pmatrix} = \begin{pmatrix} \displaystyle \frac{\partial x}{\partial \xi} & \displaystyle \frac{\partial y}{\partial \xi}\\[0.2cm] \displaystyle \frac{\partial x}{\partial \eta} & \displaystyle \frac{\partial y}{\partial \eta} \end{pmatrix} \begin{pmatrix} \displaystyle \frac{\partial \mphi_j^p}{\partial x}\\[0.2cm] \displaystyle \frac{\partial \mphi_j^p}{\partial y} \end{pmatrix} \iff \nabla_{\xi,\eta}\mphih_j(\xi,\eta) = (\JK{p})^T\nabla_{x,y}\mphi_j^p(x,y).\end{split}\]

Nous en déduisons que \(\BK{p} = (\JK{p}^T)^{-1}\), en particulier, dans le cas d’une transformation linéaire de triangle, nous obtenons :

\[\begin{split}\BK{p} = \frac{1}{\det(\JK{p})} \left( \begin{array}{c c} y_{3}^{p}-y_{1}^{p} & y_{1}^{p}-y_{2}^{p}\\ x_{1}^{p}-x_{3}^{p} & x_{2}^{p}-x_{1}^{p} \end{array} \right).\end{split}\]

Au final, comme \(X\cdot Y = X^TY\), nous obtenons

(3.5)#\[\int_{\tri_p} (\nabla\mphi_j^p)^T\nabla\mphi_i^p \diff(x,y) = \abs{\det(\JK{p})}\int_{\trih} (\nabla\mphih_j)^T (\BK{p}^T \BK{p})\nabla\mphih_i \diff (\xi,\eta).\]

En éléments finis \(\Pb^1\), les fonctions de forme sont linéaires et leur gradient est donc constant. Nous pouvons alors sortir les termes \(\nabla\mphih_i\) et \(\nabla\mphih_j\) de l’intégral pour obtenir le lemme suivant.

Lemma 3.8

Les coefficients a matrice de rigidité élémentaire \(\De{p} = ((\De{p})_{i,j})_{1\leq i,j\leq 3}\) sont obtenus pas la relation suivante

\[\begin{split}\begin{aligned} (\De{p})_{i,j} &= \int_{\tri_p}\nabla \mphi_j^p(x, y)\cdot\nabla\mphi_i^p(x,y)\diff(x,y),\\ &= \abs{\tri_p}(\nabla\mphih_j)^T (\BK{p}^T \BK{p})\nabla\mphih_i. \end{aligned}\end{split}\]

Proof. Pour les éléments finis \(\Pb^1\), les gradients \(\nabla\mphih_j\) sont constants et peuvent être sortis de l’intégrale. De plus, comme \(\abs{\det(\JK{p})} = 2\abs{\tri_p}\) et \(\abs{\trih}= \frac{1}{2}\), nous avons

\[\int_{\tri_p} \nabla\mphi_j^p\cdot\nabla\mphi_i^p \diff\xx=\abs{\tri_p}(\nabla\mphih_j)^T (\BK{p}^T \BK{p})\nabla\mphih_i.\]

3.4.4. Second membre (ou RHS ou Membre de droite)#

Étudions maintenant les termes du membre de droite comme

\[\int_{\tri_p}f(\xx)\mphi_i^p(\xx)\diff \xx.\]

Sauf pour certaines fonctions \(f\) 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 :

\[\begin{split}\begin{aligned} \displaystyle \int_{\tri_p}f(\xx)\mphi_i^p(\xx)\diff \xx &= \displaystyle \abs{\det(\JK{p})}\int_{\trih}f(\xx(\xi,\eta))\mphih_i(\xi,\eta)\diff (\xi,\eta) \\ & \displaystyle \simeq \abs{\det(\JK{p})}\sum_{m=0}^{M-1}\omega_m f(\xx(\xi_m,\eta_m))\mphih(\xi_m,\eta_m). \end{aligned}\end{split}\]

Les points \((\xi_m,\eta_m)\) 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 \(\omega_m\in\Rb\) les poids associés. Notons que le point \(\xx(\xi_m,\eta_m)\) s’obtient par l’expression vue précédemment :

\[\xx(\xi_m,\eta_m) = \sum_{i=0}^2\vertice_i^p\psih_i(\xi_m,\eta_m).\]

Nous présentons ici deux règles de quadrature pour l’intégrale \(\int_{\trih}\gh(\xx)\diff\xx\) sur \(\trih\) d’une fonction \(g\) 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) :

\(\xi_m\)

\(\eta_m\)

\(\omega_m\)

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 !).