2.4. Quadratures#

2.4.1. Rappel#

Certaines intégrales ne peuvent être calculées analytiquement et devront être approchées numériquement via des règles de quadrature. Prenons pour exemple d’une fonction \(f\) quelconque mais connue, le second membre \(B\) (un vecteur) sera alors de la forme suivante, où \(\xx=(x,y)\),

(2.7)#\[\begin{split}\begin{aligned} B[I] & = \int_{\Omega} f(\mathbf{x}) \varphi_I(\mathbf{x})\;\mathrm{d}\mathbf{x} = \sum_{p}\sum_{i} \int_{K_p} f(\mathbf{x}) \varphi_i^p(\mathbf{x})\;\mathrm{d}\mathbf{x}\\ & = \sum_{p} \abs{\det(J_p)} \sum_{i} \int_{\hat{K}} f(\mathbf{x}(\xi,\eta)) \mphih_i(\xi,\eta)\;\mathrm{d}(\xi,\eta)\\ & \approx \sum_{K_p} \abs{\det(J_p)} \sum_{i}\sum_m \omega_m f(\mathbf{x}(\xi_m, \eta_m)) \mphih_i(\xi_m, \eta_m) \end{aligned}\end{split}\]

Les poids \(\omega_m\) et les points de quadrature \((\xi_m, \eta_m)\) dépendent de la précision recherchée comme expliqué dans le cours. Rappelons aussi que \(\mathbf{x}(\xi_m, \eta_m)\) s’obtient par les fonctions d’interpolation géométrique, qui dans le cas d’éléments finis isoparamétriques, sont les mêmes que les fonctions éléments finis \(\mathbb{P}^1\), c’est à dire que pour \(\mathbf{x}\) appartenant à un élément de sommets \((\mathbf{s}_i)_i\) :

\[\mathbf{x}(\xi_m, \eta_m) = \sum_{i=1} \hat{\psi}_i(\xi_m, \eta_m)\mathbf{s}_i,\]

Les fonctions \((\hat{\psi}_i)_i\) sont égales au fonctions de forme \(\Pb^1-\) Lagrange \((\mphih_i)_i\). Cependant, les fonction \((\hat{\psi}_i)_i\) sont des fonctions d’interpolation géométriques tandis que les fonctions de forme \((\hat{\varphi}_i)_i\) sont des fonctions d’interpolation de la solution.

2.4.2. Méthodes et fonctions utiles#

Pour chaque type d’élément, Triangle ou Segment, nous avons besoin de méthodes permettant d’obtenir les quantités suivantes: :

  1. Les poids des points de quadrature

  2. Les coordonnées paramétriques des points de quadrature

  3. Les coordonnées physiques des points de quadrature

  4. Les valeurs des fonctions de forme de référence \(\hat{\varphi}\) sur des coordonnées paramétriques

Remarquez que les points 1, 2 ne dépendent que du type de l’élément considéré.

Exercise 2.5

Afin de bien compartimenter chaque fonctionnalité, nous proposons :

  • Ajouter aux classes Triangle et Segment la méthode def gaussPoint(self,order=2): qui retourne, dans le format de votre choix, les poids, les coordonnées paramétriques et les coordonnées physiques des points de Gauss de l’élement considéré et pour une précision order. Vous aurez sans doute besoin de méthodes intermédiaires pour calculer, par exemple les \(\hat{\psi}_i(\xi,\eta)\).

  • Ajouter une fonction def phiRef(element, i:int, param:[float]): qui calcule \(\hat{\varphi}_i(\xi,\eta)\) sur un élément Segment ou Triangle. L’argument param est une liste des coordonnées paramétriques (\((\xi,\eta)\) pour un triangle, \(s\) pour un segment))

2.4.3. Intégrale#

Construisez maintenant une fonction de prototype suivant

def Integrale(msh:Mesh, dim:int, physical_tag:int, f, B:np.array, order=2):

Cette fonction calcul l’intégrale \(\int f \varphi_I\) sur le domaine de tag physique physical_tag et de dimension dim. Le résultat est alors ajouté dans B[I] (voir (2.7)). L’argument f sera une fonction décrite par l’utilisatrice/utilisateur, elle prendra 2 arguments, x et y et retournera un scalaire correspondant à \(f(x,y)\).