Eléments finis de Lagrange

On va présenter dans ce chapitre le type le plus simple et le plus classique d’éléments finis.

1. Unisolvance

Définition

Soit Σ={a1,,aN} un ensemble de N points distincts de Rn. Soit P un espace vectoriel de dimension finie de fonctions de Rn à valeurs dans R.On dit que Σ est P-unisolvant ssi pour tous réels α1,,αN, il existe un unique élément p de P tel que p(ai)=αi,i=1,,N.

Ceci revient à dire que la fonction :

L:PRNp(p(a1),,p(aN))

est bijective.

Preuve:

En pratique, on montrera que Σ est P-unisolvant en vérifiant que dimP=CardΣ, puis en montrant l’injectivité ou la surjectivité de L.

L’injectivité de L se démontre en établissant que la seule fonction de P s’annulant sur tous les points de Σ est la fonction nulle.

La surjectivité de L se démontre en exhibant une famille p1,,pN d’éléments de P tels que pi(aj)=δij, c’est à dire un antécédent pour L de la base canonique de RN.

En effet, étant donnés des réels α1,,αN, la fonction p=Ni=1αipi vérifie alors p(aj)=αj,j=1,,N.

2. Elément fini de Lagrange

Élément fini de Lagange

Un élément fini de Lagrange est un triplet (K,Σ,P) tel que

  • K est un élément géométrique de Rn (n= 1, 2 ou 3), compact, connexe, et d’intérieur non vide.

  • Σ={a1,,aN} est un ensemble fini de N points distincts de K.

  • P est un espace vectoriel de dimension finie de fonctions réelles définies sur K, et tel que Σ soit P-unisolvant (donc dim P=N).

Fonctions de base locales

Soit (K,Σ,P) un élément fini de Lagrange. On appelle fonctions de base locales de l’élément les N fonctions pi (i=1,,N) de P telles que

pi(aj)=δij1i,jN.

On vérifie aisément que (p1,,pN) ainsi définie forme bien une base de P.

Définition

On appelle opérateur de interpolation (ou encore Pinterpolation) sur Σ l’opérateur πK qui, à toute fonction v définie sur K, associe la fonction πKv de P définie par πKv=Ni=1v(ai)pi.

πKv est donc l’unique élément de P qui prend les mêmes valeurs que v sur les points de Σ.

3. Exemples d’éléments finis de Lagrange

3.1. Espaces de polynômes

On notera Pk l’espace vectoriel des polynômes de degré total inférieur ou égal à k.

  • Sur , Pk=Vect{1,X,,Xk} et dim Pk=k+1.

  • Sur 2, Pk=Vect{XiYj;0i+jk} et dim Pk=(k+1)(k+2)2.

  • Sur 3, Pk=Vect{XiYjZl;0i+j+lk} et dim Pk=(k+1)(k+2)(k+3)6.

On notera Qk l’espace vectoriel des polynômes de degré inférieur ou égal à k par rapport à chaque variable.

  • Sur R, Qk=Pk.

  • Sur R2, Qk=Vect{XiYj;0i,jk} et dim Qk=(k+1)2.

  • Sur R3, Qk=Vect{XiYjZl;0i,j,lk} et dim Qk=(k+1)3.

3.2. Exemples 1-D

Elément P1
  • K=[a,b]

  • Σ={a,b}

  • P=P1

Elément P2
  • K=[a,b]

  • Σ={a,a+b2,b}

  • P=P2

Elément Pk
  • K=[a,b]

  • Σ={a+ibam,i=0,,k}

  • P=Pk

3.3. Exemples 2-D triangulaires

Elément P1
  • K=triangle de sommets a1,a2,a3

  • Σ={a1,a2,a3}

  • P=P1 Les fonctions de base sont définies par pi(aj)=δij. Ce sont donc les coordonnées barycentriques : pi=λi (cf annexe [ch:bary]).

Elément P2
  • K=triangle de sommets a1,a2,a3

  • Σ={a1,a2,a3,a12,a13,a23}, où aij=ai+aj2.

  • P=P2

Les fonctions de base sont pi=λi(2λi1) et pij=4λiλj. Un exemple de calcul de ces fonctions de base est donné en annexe [ch:bary].

Eléments finis triangulaire P1, triangulaire P2 et rectangulaire Q1

elements 2D

3.4. Exemples 2-D rectangulaires

Elément Q1
  • K=rectangle de sommets a1,a2,a3,a4, de côtés parallèles aux axes

  • Σ={a1,a2,a3,a4}

  • P=Q1 Les fonctions de base sont pi(X,Y)=(Xxj)(Yyj)(xixj)(yiyj), où (xi,yi) sont les coordonnées de ai, et où aj, de coordonnées (xj,yj) est le coin opposé à ai.

3.5. Exemples 3-D

Elément tétraèdrique P1
  • K=tétraèdre de sommets a1,a2,a3,a4

  • Σ={a1,a2,a3,a4}

  • P=P1

Elément tétraèdrique P2
  • K=tétraèdre de sommets a1,a2,a3,a4

  • Σ={ai}1i4{aij}1i<j4

  • P=P2

Les fonctions de base sont pi=λi(2λi1) et pij=4λiλj.

Elément parallélépipèdique Q1
  • K=parallélépipède de sommets a1,,a8 de côtés parallèles aux axes

  • Σ={ai}1i8

  • P=Q1

Elément prismatique
  • K=prisme droit de sommets a1,,a6

  • Σ={ai}1i6

  • P={p(X,Y,Z)=(a+bX+cY)+Z(d+eX+fY),a,b,c,d,e,fR}

Eléments finis tétraèdriques P1 et P2, parallélépipèdique Q1, et prismatique

elements 3D

4. Famille affine d’éléments finis

Équivalence affine

Deux éléments finis (ˆK,ˆΣ,ˆP) et (K,Σ,P) sont affine-équivalents ssi il existe une fonction affine F inversible (F:ˆxBˆx+b) telle que (i) K=F(ˆK) (ii) ai=F(ˆai)i=1,,N et (iii) P={ˆpF1,ˆpˆP}

Si l’on est dans Rn, B est donc une matrice n×n inversible, et b est un vecteur de Rn.

Soient (ˆK,ˆΣ,ˆP) et (K,Σ,P) deux éléments finis affine-équivalents, via une transformation F.

On note ˆpi(i=1,,N) les fonctions de base locales de ˆK.

Alors les fonctions de base locales de K sont les pi=ˆpiF1.

Définition

On appelle famille affine d’éléments finis une famille d’éléments finis tous affine-équivalents à un même élément fini (ˆK,ˆΣ,ˆP), appelé élément de référence.

D’un point de vue pratique, le fait de travailler avec une famille affine d’éléments finis permet de ramener tous les calculs d’intégrales à des calculs sur l’élément de référence.

Les éléments de référence sont :

  • En 1-D : le segment [0,1]

  • En 2-D triangulaire : le triangle unité, de sommets (0,0), (0,1) et (1,0).

  • En 2-D rectangulaire : le carré unité [0,1]×[0,1].

  • En 3-D tétraèdrique : le tétraèdre unité, de sommets (0,0,0), (1,0,0), (0,1,0) et (0,0,1).

  • En 3-D parallélépipèdique : le cube unité [0,1]×[0,1]×[0,1].

  • En 3-D prismatique : le prisme unité de sommets (0,0,0), (0,1,0), (1,0,0), et (0,0,1), (0,1,1), (1,0,1).

5. Maillages

Nous étendons ici aux dimension 2 et 3 les notions élémentaires de maillage vues en 1D, voir la figure Maillage Feel++.

Maillage Feel++

feelpp mesh

Définition

Un maillage est constituée d’une famille d’éléments(ou mailles ou cellules) {Ke}e=1,...,NeNe est le nombre d’éléments, nous noterons

Th={Km}m=1,...,Ne

avec

h=max

et

h_{K_e} = \diam(K_e) = \max_{x_1,x_2 \in K_e} \|x_1-x_2\|,\, e \in \{1,...,\Ne\}

On travaille par la suite avec des familles de maillage et on les note \set{\mathcal{T}_h}_{h > 0}.

Famille de maillage quasi-uniforme

On dira qu’une famille de maillage \set{\mathcal{T}_h}_{h > 0} est quasi-uniforme s’il existe une constante c telle que

\forall h,\ \forall K \in \calTh,\ h_K \geq c h
Cela veut dire que les élements sont tous de la même taille pour h donné.

5.1. Transformation géométrique

Un maillage est généré par

  1. un élément de reference noté \hat{K}

  2. une famille de transformations géométriques mappant \hat{K} vers les éléments K_e, e=1,\ldots,\Ne dans le maillage

Nous supposerons que les transformations sont des \mathcal{C}^1- diffeomorphismes [1].

Définition

Pour une cellule K \in \mathcal{T}_h, on note T_K la transformation géométrique

T_K: \hat{K} \mapsto K

Afin de spécifier la transformation géométrique, on considère l’élément fini de Lagrange, noté (\hat{K},\hat{P}_{\mathrm{geo}}, \hat{\Sigma}_{\mathrm{geo}}), tel que

  • \ngeo = \card{\hat{\Sigma}_{\mathrm{geo}}}

  • \set{\hat{g}_1,\ldots,\hat{g}_{\ngeo}} les noeuds de \hat{K}

  • \set{\hat{\psi}_1,\ldots,\hat{\psi}_{\ngeo}} les fonctions de forme

Élément fini géométrique

On dit que (\hat{K},\hat{P}_{\mathrm{geo}}, \hat{\Sigma}_{\geo}) est l’élément fini géométrique, \set{\hat{g}_1,\ldots,\hat{g}_{\ngeo}} sont les noeuds géométriques et \set{\hat{\psi}_1,\ldots,\hat{\psi}_{\ngeo}} sont les fonctions de formes géométriques

Transformation géométrique associée à un triangle

geomap triangle

Pour chaque K \in \mathcal{T}_h, on a un \ngeo-uplet \set{g^K_1,\ldots,g^K_\ngeo}.

La transformation géométrique est définie comme suit

T_K: \hat{x} \in \hat{K} \mapsto \sum_{i=1}^\ngeo\ g^K_i \hat{\psi}_i(\hat{x})

et en particulier on a

T_K(\hat{g}_i) = g^K_i, \quad \forall i \in \set{1,\ldots,\ngeo}
On a T_K \in [\hat{P}_\geo(\hat{K})^d] et que \set{g^K_1,\ldots,g^K_\ngeo} sont les noeuds géométriques de K.

T_K est un \mathcal{C}^1-diffeomorphism donc la numérotation des noeuds \set{g^K_1,\ldots,g^K_\ngeo} doit être compatible avec les noeuds de l’élément finit géométrique.

La numérotation locale des entités géométriques dans doit être consistente avec la numérotation locale des générateurs de maillage. voir \triangleright format de fichier Gmsh pour une numérotation locale.

Un cas particulier est la transformation géométrique affine.

Définition: Maillage affine

Quand toutes les transformations géométriques \set{T_K}_{K \in \mathcal{T}_h} sont affines, cela veut dire que pour tout K \in \mathcal{T}_h, il existe un vecteur b_K \in \RR^d et une matrice J_K \in \RR^{d\times d} tels que

T_K : \hat{x} \in \hat{K} \mapsto b_K + J_K \hat{x} \in K

On dit que le maillage est affine.

Si l’élément fini géométrique est (\hat{K},\poly{P}_1,\Sigma_\ngeo) alors les éléments K sont soit des triangles soit des tétrahèdres.
Mesh<Simplex<d,1> > ou Mesh<Simplex<d> > est le type pour les maillages affines formés de simplexes dans \RR^d. 1 indique l’ordre de l’élément fini géométrique et est la valeur par défaut.

5.2. Quelques calculs avec la transformation géométrique

5.2.1. Gradient, Inverse et Jacobien

On note \xi un ensemble de n points dans \hat{K} et on note \nabla T_K(\xi) le gradient de T_K aux points \xi

\nabla T_K( \xi )\ =\ \sum_{i=0}^{\ngeo}\ g^K_i\ \nabla \psi_i (\xi)

et B_K(\xi) = \nabla T_K^{-1}(\xi) l’inverse \xi et finalement J_K(\xi) le jacobien de T_K en \xi

J_K(\xi)\ =\ |\det( \nabla T_K(\xi) )|

5.2.2. Dérivation dans l’élément de référence

Afin de dériver un polynome dans l’élément réel K, grâce à la transformation géométrique et la règle de différentiation des fonctions composées, nous dérivons seulement dans l’élément de référence \hat{K}.

Soit f: K \mapsto \RR et \hat{f}: \hat{K} \mapsto \RR telle que \hat{f} = f \circ T_K

\nabla f\ =\ \hat{\nabla} \underbrace{\hat{f}(\xi)}_{f \circ T_K(\xi)} B_K(\xi)

en 2D, on a, en notant x=T_K(\xi),

\nabla f(x) = \begin{pmatrix} \frac{\hat{\partial} \hat{f} (\xi)}{\partial \xi_1} & \frac{\hat{\partial} \hat{f} (\xi)}{\partial \xi_2} \end{pmatrix} \begin{pmatrix} B_{K_{11}}(\xi) & B_{K_{12}}(\xi)\\ B_{K_{21}}(\xi) & B_{K_{22}}(\xi)\\ \end{pmatrix}

5.2.3. Intégration dans l’élément de référence

De manière similaire, au lieu de calculer les intégrales sur l’élément réel K, nous appliquons un changement de variables et calculons les intégrales sur l’élément de réference \hat{K}.

Soit f: K \mapsto \RR et \hat{f}: \hat{K} \mapsto \RR telle que \hat{f} = f \circ T_K, et {\mathbf{F}}: K \mapsto \RR^d et {\hat{\mathbf{F}}}: \hat{K} \mapsto \RR^d telle que \hat{{\mathbf{F}}} = {\mathbf{F}} \circ T_K,

On a alors les relations suivantes:

\int_{K} \ f\ dx\ =\ \int_{\hat{K}} f(T_K(\xi) ) J_K( \xi )\ d \xi \ =\ \int_{\hat{K}} \hat{f}(\xi) J_K( \xi )\ d \xi
\int_{K}\ \nabla f\ dx\ =\ \int_{\hat{K}} \Big(\hat{\nabla} \underbrace{\hat{f}(\xi)}_{f \circ T_K(\xi)} B_K(\xi)\Big) J_K( \xi )\ d \xi
\int_{\partial K} f( x ) dx = \int_{\partial \hat{K}} \hat{f}(\xi) \| B_K(\xi) {\hat{\mathbf{n}}}(\xi) \| J_K( \xi ) d \xi
\int_{\partial K}\ {\mathbf{F}}( x )\ \cdot\ {\mathbf{n}}(x) dx = \int_{\partial \hat{K}} {\hat{\mathbf{F}}}( \xi )\ \cdot \Big(B_K(\xi) {\hat{\mathbf{n}}}(\xi) \Big) J_K( \xi )\ d \xi

{\mathbf{n}}(x) est la normale extérieure unitaire à \partial K évaluée en x \in \partial K, et {\hat{\mathbf{n}}}(\xi) la normale unitaire extérieure à \hat{K} évaluée en \xi \in \partial \hat{K}.

Feel++ effectue automatiquement pour vous les changements de variables dans les intégrales.
Définition: Maillage conforme

Un maillage est dit conforme si l’intersection de deux éléments est soit vide, un sommet, une arête ou une face.

On ne manipule que des maillages conformes dans le cours mais Feel++ peut traiter des maillages non-conformes par exemple dans le contexte de méthode de décomposition de domaines.

6. Espaces élément fini de Lagrange

Soit \mathcal{T}_h un maillage généré par (\hat{K}, \hat{P}_{\mathrm{geo}}, \hat{\Sigma}_{\mathrm{geo}}), une cellule K \in \mathcal{T}_h est alors l’image de \hat{K} par la transformation géométrique T_K défini par [eq:23].

L’objectif à présent est de générer la famille d’éléments finis de Lagrange grâce à l’élément fini de référence (\hat{K},\hat{P}, \hat{\Sigma})

\{(K,P_K,\Sigma_K)\}_{K \in \mathcal{T}_h}
Notations

On note \{\hat{x}_1,...,\hat{x}_{\nf}\} les noeuds de l’élément fini.

On note \{\hat{\Psi}_1,...,\hat{\Psi}_{\nf}\} les fonctions de forme élément fini.

Proposition

Soit K \in \mathcal{T}_h et P_K=\{\hat{p} \circ T_K^{-1}; \hat{p} \in \hat{P}\}.

Pour tout i \in \{1,...,\nf\}, on pose x_{K,i} = T_K(\hat{x}_i)

On note \Sigma_K l’ensemble des degrés de liberté associé à \{x_{K,1},...,x_{K,\nf}\}

Alors (K,P_K,\Sigma_K) est un élément fini de Lagrange.

Les fonctions de forme sont définies de la façon suivante

\Psi_{K,i} = \hat{\Psi}_i \circ T_K^{-1}, \quad i=1,...,\nf

et l’opérateur d’interpolation local comme

\Ilag{K}: v \in \mathcal{C}^0(K) \mapsto \sum_{i=1}^\nf\ v(x_{K,i})\ \Psi_{K,i}\ \in P_K

Une propriété importante de \Ilag{K} est que

\forall v \in \mathcal{C}^0(K),\quad \Ilag{K}( v \circ T_K ) = \Ilag{K}( v ) \circ T_K.
Théorème

Soit T_K une transformation affine.

Soit \mathbb{P}_k \subset \hat{P} et k+1 > \frac{d}{2}.

Soit h_K le diamètre de K et \rho_K le diamètre de la plus grande boule inscrite dans K et \omega_K = \frac{h_K}{\rho_K}.

Alors il existe une constante c independente de K telle que \forall v \in H^{k+1}(K) et pour tout m \in \{0,...,k+1\},

|v - \Ilag{K}(v)|_{m,K} \leq c h^{k+1-m} \omega_K^m |v|_{k+1,K}
  • \omega_K devrait être aussi proche de 1 que possible

  • La deuxième hypothèse technique permet d’assurer que H^{k+1}(K) \subset C^0(K).

  • on obtient des résultats similaires si v n’est pas suffisamment régulière

Définition: Espace H^1-conforme

Un espace vectoriel V_h de fonctions définies sur un domaine \Omega_h est dit être H^1-conforme si V_h \subset H^1(\Omega_h).

Afin de construire un tel espace on introduit tout d’abord

W_h = \{v_h \in L^2(\Omega_h); \forall K \in \mathcal{T}_h, v_h|_K \in P_K\}

mais ce n’est pas suffisant: les fonctions W_h peuvent avoir des sauts entre les éléments du maillage. Nous avons donc besoin d’assurer la continuité de ces fonctions

V_h = W_h \cap C^0(\Omega_h) = \{ v_h \in W_h; \forall F \in \mathcal{F}^i_h. \jump{v_h}_F = 0\}

Concernant l’implémentation, nous avons besoin de d’indentifier les degrés de liberté communs entre les éléments quand nous construisons la tables des degrés de liberté.

Voici deux exemples d’espace H^1-conforme

Espace H^1-conforme sur des simplexes (e.g. triangle)
P^k_{c,h} = \{ v_h \in C^0(\Omega_h); \forall K \in \mathcal{T}_h, v_h \circ T_K \in \mathbb{P}_k\}
Espace H^1-conforme sur des hypercubes (e.g quadrilatère)
Q^k_{c,h} = \{ v_h \in C^0(\Omega_h); \forall K \in \mathcal{T}_h, v_h \circ T_K \in \mathbb{Q}_k\}
Table 1. Dimension des espaces P_{\mathrm{c}, h}^k et Q_{\mathrm{c}, h}^k en dimension 2 et 3 pour k \in\{1,2,3\}.
espace d k=1 k=2 k=3

P_{c,h}^k

2

N_{so}

N_{so}+N_{ar}

N_{so}+2N_{ar}+N_{ma}

Q_{c,h}^k

2

N_{so}

N_{so}+N_{ar}+N_{ma}

N_{so}+2N_{ar}+4N_{ma}

P_{c,h}^k

3

N_{so}

N_{so}+N_{ar}

N_{so}+2N_{ar}+N_{fa}

Q_{c,h}^k

3

N_{so}

N_{so}+N_{ar}+N_{fa}+N_{ma}

N_{so}+2N_{ar}+4N_{fa}+8N_{ma}

6.1. Projections orthogonales

On note

\begin{aligned} \Pi^{0,k}_{c,h} : L^2(\Omega) \rightarrow P^k_{c,h}\\ \Pi^{1,k}_{c,h} : H^1(\Omega) \rightarrow P^k_{c,h} \end{aligned}

les projections orthogonales associées respectivement aux produits scalaires (u,v)_{0,\Omega} = \int_\Omega u v and (u,v)_{1,\Omega} = \int_\Omega u v + \int_\Omega \nabla u \cdot \nabla v

On a pour l=1,...,k et si v \in H^{l+1}(\Omega)

\begin{array}{rl} \|v - \Pi^{0,k}_{c,h}(v)\|_{0,\Omega} & \leq c h^{l+1} |v|_{k+1,\Omega} \\ \|v - \Pi^{1,k}_{c,h}(v)\|_{1,\Omega} & \leq c h^{l} |v|_{k+1,\Omega} \end{array}

Pour calculer \Pi^{0,k}_{c,h} et \Pi^{1,k}_{c,h} on a besoin de résoudre les problèmes

Problème: Projection L^2

Soit v une fonction de L^2, calculer \Pi^{0,k}_{c,h}(v) \in P^k_{c,h} tel que \forall v_h \in P^k_{c,h} on a

(\Pi^{0,k}_{c,h}(v), v_h)_{0,\Omega} = (v, v_h)_{0,\Omega}
Projection H^1

Soit v une fonction de H^1, calculer \Pi^{1,k}_{c,h}(v) \in P^k_{c,h} tel que \forall v_h \in P^k_{c,h} on a

(\Pi^{1,k}_{c,h}(v), v_h)_{1,\Omega} = (v, v_h)_{1,\Omega}
L’opérateur \Pi_{\mathrm{c}, b}^{1, k} est appelé projecteur elliptique. On peut également considérer la projection orthogonale de H_0^1(\Omega) dans P_{\mathrm{c}, b, 0}^k avec le produit scalaire \int_{\Omega} \nabla v \cdot \nabla w. Ce dernier satisfait les mêmes propriétés que celles énoncées ci-dessus pour l’opérateur \Pi_{\mathrm{c}, b}^{1, k}.
Des résultats analogues existent pour des opérateurs de projection orthogonale sur les espaces Q_{c, b}^k.

6.2. Interpolant de Lagrange

Définition

Notons V_h un espace H^1 conforme, \{\Psi_i\}_{1,...,N} une base nodale de V_h et \{x_1,...,x_N\} les noeuds associés alors l’interpolant de Lagrange est défini par

\Ilag{h}: v \in C^0(\Omega_h) \mapsto \sum_{i=1}^N v(x_{i}) \Psi_i \in V_h
La restriction de l’interpolant de Lagrange à une cellule K coincide avec l’interpolant de Lagrange appliqué à la fonction dans la cellule K:
\Ilag{h}(v)|_K = \Ilag{h}(v|_K)
Théorème

Supposons \{\mathcal{T}_h\}_{h>0} une famille de maillage quasi-uniforme et conformes, \mathbb{P}_k \subset \hat{P} et k+1 > \frac{d}{2}.

Alors il existe une constante c telle que pour tout h et v \in H^{k+1}(\Omega_h)

\|v - \Ilag{h}(v)\|_{0,\Omega_h} + h |v - \Ilag{h}(v)|_{1,\Omega_h} \leq c h^{k+1} |v|_{k+1,\Omega_h}

De plus, pour tout m \in\{2, \ldots, k+1\}, on a

\begin{equation*} \left(\sum_{K \in \mathcal{T}_b}\left|v-\mathcal{I}_b^{\mathrm{Lag}} v\right|_{m, K}^2\right)^{\frac{1}{2}} \leqslant c h^{k+1-m}|v|_{k+1, \Omega_b} . \end{equation*}

Ce théorème est une conséquence directe du théorème d’interpolation local. On notera que l’hypothèse portant sur la régularité de la famille \left\{\mathcal{T}_h\right\}_{h>0} permet de contrôler uniformément les rapports \varpi_K. Enfin, lorsque la fonction à interpoler n’est pas suffisamment régulière, par exemple si v \in H^s\left(\Omega_h\right) et v \notin H^{s+1}\left(\Omega_h\right) pour un entier s \in\{1, \ldots, k\}, l’estimation [eq.62] ci-dessus devient

\begin{equation*} \left\|v-\mathcal{I}_b^{\mathrm{Lag}} v\right\|_{0, \Omega_b}+h\left|v-\mathcal{I}_b^{\mathrm{Lag}} v\right|_{1, \Omega_b} \leqslant c h^s|v|_{s, \Omega_b} \end{equation*}
Example 1. Vérification du théorème précédent sur un exemple

Nous considérons pour cela \alpha un réel et \mathbf{x}=(x_1,...,x_d) \in \mathbb{R}^d un point de \Omega = [0,1]^d, d=1,2,3 et v la fonction définie par

\begin{array}{rl} v : \Omega &\rightarrow \mathbb{R}\\ \mathbf{x} &\rightarrow ( \mathbf{x} \cdot \mathbf{x} )^{\alpha/2}\ \Pi_{i=1}^d( 1-x_i^2) \end{array}

Nous construisons l’interpolant de Lagrange de v dans P^k_{c,h} avec k=1,...,5 et d=1,2,3 et étudions l’erreur d’interpolation L^2 et H^1 du théorème [eq:62] en échelle log-log.

Nous devons obtenir des droites de pentes k (resp. k+1) pour la norme L^2 (resp. H^1.)

7. Approximation iso-paramétrique

Quand le domaine est courbe, si nous désirons obtenir des propriétés de convergence optimale nous avons besoin de discretiser le bord du domaine avec suffisamment de précision

Notons (\hat{K},\hat{P}_{\mathrm{geo}}, \hat{\Sigma}_{\geo}) l’élément fini géométrique et (\hat{K},\hat{P}, \hat{\Sigma}) l’élément fini de référence pour V_h.

Définition
  • Si \hat{P}_{\mathrm{geo}} = \hat{P} l’approximation est dite iso-parametrique.

  • Si \hat{P}_{\mathrm{geo}} \subset \hat{P} l’approximation est dite sub-parametrique.

  • Si \hat{P} \subset\hat{P}_{\mathrm{geo}} l’approximation est dite sur-parametrique.

Gmsh cite:[gmsh] un mailleur libre permet de générer des maillage d’ordre élevé jusqu’à l’ordre 5 en 2D et 4 en 3D
Théorème

Supposons que \{\mathcal{T}_h\}_{h>0} une famille de maillage quasi-uniformes et conformes, \mathbb{P}_k \subset \hat{P} et k+1 > \frac{d}{2} et k_{\mathrm{geo}} = k.

Alors il existe une constante c telle que pour tout h et v \in H^{k+1}(\Omega_h)

\|v - \Ilag{h}(v)\|_{0,\Omega_h} + h |v - \Ilag{h}(v)|_{1,\Omega_h} \leq c h^{k+1} |v|_{k+1,\Omega_h}
Les résultats sont identiques à ceux du theorème Théorème.

Nous allons vérifier sur un exemple à l’erreur d’approximation de l’élément géométrique. Considérons les cercles unité généré par une transformation affine, noté \Omega^1_h, et d’ordre 2, noté \Omega^2_h, et calculons leur aire respective.

Construisons une famille de maillage \{\calTh\}_{h >0}, par exemple h\in \{0.4, 0.2, 0.1, 0.05\} et calculons l’erreur entre le calcul exact de l’aire \pi et le calcul numérique \int_{\Omega^1_h} 1 et \int_{\Omega^2_h} 1 respectivement.

8. Feel++

En Feel++, un maillage Mesh est décomposé en un ensemble d’éléments décomposés en sous entités (volume,face,arête,point):

  • faces(decomposés en sous entités),

  • arêtes(decomposés en sous entités) et

  • points.

et à chaque élément K est associé une transformation géométrique T_K.

L’ordre polynomial de la transformation géométrique est donné par le second argument template

Type d’un object maillage en Feel++
Mesh<Simplex<d, k_geo>> (1)
Mesh<Hypercube<d, k_geo>> (2)
cpp
1 Un maillage de simplexes en dimension d avec transformation géométrique d’ordre k_{\mathrm{geo}}.
2 Un maillage d’hypercubes en dimension d avec transformation géométrique d’ordre k_{\mathrm{geo}}.

Ain de parcourir les éléments et faces du maillage, Feel++ fournit des fonctions renvoyant des itérateurs (début et fin) sur ces ensembles

elements(mesh)

retourne 2 itérateurs sur l’ensemble des éléments du maillage

markedelements(mesh,<int>)

et markedelements(mesh,<string>) retourne 2 itérateurs sur les éléments marqués par l’entier <int> et la chaîne des caractères <string> respectivement, cela correspondra typiquement à des propriétés de matériau

boundaryfaces(mesh)

retourne 2 itérateurs sur les faces au bord du maillage

markedfaces(mesh,<int>) et markedelements(mesh,<string>)

retourne 2 itérateurs sur les faces marquées par l’entier <int> et la chaîne des caractères <string> respectivement, ca correspondra typiquement aux conditions aux limites

link pending to Feel++ user manual mesh iterators section.

L’espace d’approximation V_h H^1 conforme (espaces de functions continues sur \Omega polynomiales par morceaux de degré \leq k) est défini comme suit

FunctionSpace<Mesh<Simplex<d, k_geo> >, bases<Lagrange<k> > > V_h; (1)
FunctionSpace<Mesh<Hypercube<d, k_geo> >, bases<Lagrange<k> > > V_h; (2)
cpp
1 \Pch{k}
2 \Qch{k}

9. Du problème global aux éléments locaux

On va maintenant faire le lien entre la résolution d’un problème par méthode d’éléments finis et les notions qui viennent d’être introduites.

Soit une EDP à résoudre sur un domaine \Omega, et V l’espace de Hilbert dans lequel on cherche une solution de la formulation variationnelle du problème.

On réalise un maillage de \Omega par une famille affine de N_e éléments finis (K_i,\Sigma_i,P_i)_{i=1,\ldots,N_e}.

Par unisolvance, la solution approchée u_h sera entièrement définie sur chaque élément (K_i,\Sigma_i,P_i) par ses valeurs sur les points de \Sigma_i, qu’on appellera les noeuds du maillage.

Il est à noter qu’un noeud sera en général commun à plusieurs éléments adjacents.

Le nombre total de noeuds N_h est donc inférieur à N_e\times\hbox{Card} \Sigma_i, et on a dim V_h = N_h.

Notons a_1,\ldots,a_{N_h} les noeuds du maillage.

Le problème approché se ramène donc à la détermination des valeurs de u_h aux points a_i: ce sont les degrés de liberté du problème approché.

On va construire une base de V_h en associant à chaque ddl a_i un vecteur de la base. On définit ainsi les fonctions de base globales \varphi_i (i=1,\ldots,N_h) par

Fonctions de base globales
{\varphi_i}_{|K_j} \in P_j, \quad j=1,\ldots,N_e \mbox{ et } \varphi_i(a_j)=\delta_{ij}, 1\le i,j \le N_h

L’espace d’approximation interne est donc alors :

V_h = \hbox{Vect }\left\{\varphi_1,\ldots,\varphi_{N_h}\right\}
Exemple de fonction de base globale (élément triangulaire P_1)

fonction globale

On remarque qu’une telle fonction \varphi_i est nulle partout, sauf sur les éléments dont a_i est un noeud.

En effet, si a_i n’appartient pas à un élément K, \varphi_i est nulle sur tous les noeuds de K, et donc sur K tout entier par unisolvance.

De plus, sur un élément K dont a_i est un noeud, \varphi_i vaut 1 sur a_i et 0 sur les autres noeuds de K.

Donc \varphi_i\, _{|K} est une fonction de base locale de K.

la fonction de base globale \varphi_i est donc construite comme la réunion des fonctions de base locales sur les éléments du maillage dont a_i est un noeud.

C’est à ce niveau que se situe le lien entre les définitions locales introduites au Elément fini de Lagrange et le problème global approché à résoudre.

Par ailleurs, ceci implique que tous les calculs à effectuer sur les fonctions de base globales peuvent se ramener à des calculs sur les fonctions de base locales, et donc simplement à des calculs sur l’élément de référence (car on a maillé le domaine avec une famille d’éléments finis affine-équivalents).

Maillage non conforme

conforme

Ce type de définition des fonctions de base n’est possible que si le maillage est conforme, c’est à dire si l’intersection entre deux éléments est soit vide, soit réduite à un sommet ou une arête en dimension 2 (ou à un sommet, une arête ou une face en dimension 3).

On interdit ainsi les situations du type de celle de la figure.

10. Exercices

  1. Calculer les fonctions de base locales des éléments finis de Lagrange introduits dans ce chapitre.

  2. Donner l’allure des fonctions de base globales correspondantes. Sont-elles continues ? dérivables ?

  3. Pour les éléments finis de Lagrange introduits dans ce chapitre, écrire le changement de variable affine entre élément quelconque et élément de référence.


1. la transformation et son inverse sont \mathcal{C}^1 et bijectives