III.A Mise en uvre informatique de l'espace V_h choisi pour Helmholtz.

Nous calculons analytiquement les termes des matrices et du second membre de la formulation discrète (I.1.27) du problème modèle de Helmholtz bidimensionnel sans coefficient. L'espace de discrétisation est l'espace V_h construit à l'aide des ondes planes section I.2.2.1. Nous donnons d'abord les notations nécessaires à cette section.

i)

Notations relatives à la géométrie (sommets, arêtes, normales, longueurs).

1. Soient (x₁^k,x₂^k,x₃^k) les positions des trois sommets du triangle $\Omega _k$.
2. Soient (x₁^kj,x₂^kj) les extrémités de l'arête $\Sigma_{kj}$ commune aux éléments $\Omega _k$ et $\Omega_j$.
3. Soient $( \nu_1^k,\nu_2^k,\nu_3^k)$ les trois normales sortantes du triangle $\Omega _k$, respectivement des arêtes (x₁^k,x₂^k), (x₂^k,x₃^k) et (x₃^k,x₁^k).
4. Soit $(\nu_{kj})$ la normale sortante du triangle $\Omega _k$ vers le triangle $\Omega_j$. C'est donc aussi la normale à l'arête (x₁^kj,x₂^kj).
5. Soient (L₁^k,L₂^k,L₃^k) les longueurs des trois arêtes du triangle $\Omega _k$: |x₁^k-x₂^k|, |x₂^k-x₃^k|, |x₃^k-x₁^k|.
6. Soit L_kj la longueur de l'arête $\Sigma_{kj}$: |x₁^kj-x₂^kj|.

ii)

Notations relatives aux termes exponentiels des fonctions de base.

1. Pour n=1 à 3, $h_n^k = \omega (\mathbf{v}_{km}-\mathbf{v}_{kl}) \, \frac{\displaystyle (\vec{x}_{n+1}^k-\vec{x}_n^k)}{\displaystyle 2}$ avec $\vec{x}_4^k =\vec{x}_1^k$.
2. Pour n=1 à 3, $Z_n^k = e^{(i \omega (\mathbf{v}_{km}-\mathbf{v}_{kl}) \, \vec{x}_n^k )}$.
3. Soit $h_{kj} = \omega (\mathbf{v}_{jm}-\mathbf{v}_{kl}) \, \frac{\displaystyle (\vec{x}_2^{kj}-\vec{x}_1^{kj})}{\displaystyle 2}$.
4. Soit $Z_{kj} = e^{(i \omega (\mathbf{v}_{jm}-\mathbf{v}_{kl}) \, \vec{x}_1^k )}$.

Les matrices D (I.2.18) et C (I.2.19 et I.2.20) et le second membre b (I.2.21) donnés section I.2.1.4 sont calculés à l'aide des ondes planes définissant l'espace de discrétisation V_h.

i)

Le calcul de D (I.2.18) pour des ondes planes fait apparaître un terme constant à multiplier par l'intégrale d'une onde plane sur un segment. A l'aide de l'annexe III.D.1.1 qui donne l'intégrale sur un segment de la fonction e^i<I>k<I>X, on calcule facilement que les termes non nuls de D, notés D_k^l,m, sont explicités par

$$D_{k}^{l,m} = \omega^2 \sum_{n=1}^{3}{L_n^k Z_n^k (1-\vec{\... ...isplaystyle \sin{h_n^k}}{\displaystyle h_n^k} e^{ih_n^k} } \ .$$ (III.A.1)

Le terme D_k^l,m est la somme des contributions des arêtes de l'élément $\Omega _k$ pour les fonctions de base de directions v_km et v_kl.

ii)

La matrice C, constituée de termes de couplage entre deux éléments voisins (I.2.19) et de termes diagonaux de couplage au bord du domaine (I.2.20), est calculée de la même façon que D. On obtient donc après calculs (qui utilisent toujours l'intégrale sur un segment de la fonction e^i<I>k<I>X de l'annexe III.D.1.1) les deux formules de calcul des termes non nuls de C (III.A.2) et (III.A.3) ci-dessous.

1. Dans le cas où $\Omega _k$ est voisin de $\Omega_j$, la contribution C_k,j^l,m de l'interface $\Sigma_{kj}$ pour les fonctions de base de directions v_jm et v_kl est donnée par
   

   $$C_{k,j}^{l,m} = \omega^2 L_{kj} Z_{kj} (1+\vec{\nu}_{kj}.\ma... ... } \frac{\displaystyle \sin{h_{kj}}}{\displaystyle h_{kj}} \ .$$ (III.A.2)

   
   

2. Dans le cas où $\Gamma_{k}$ est non vide, la contribution C_k,k^l,m du bord $\Gamma_{k}$ pour les fonctions de base sur $\Omega _k$ de directions v_km et v_kl est donnée par
   

   $$C_{k,k}^{l,m} = \sum_{n / [x_n^k,x_{n+1}^k] \in \Gamma_k}{{{... ...h_n^k } \frac{\displaystyle \sin{h_n^k}}{\displaystyle h_n^k}}$$ (III.A.3)

   
   
   où l'on utilise l'hypothèse ${{\mathit Q}}_k = {{\mathit Q}}_{\vert\Gamma_k}$ constant sur $\Gamma_{k}$.

iii)

Le second membre, donné par l'équation (I.2.21), est calculé pour des fonctions de base issues d'ondes planes par

$$b_{k,l} = { -2i \omega \displaystyle \int _{ \Omega_k} f. \... ..._k} + i \omega ) e^{ (i\omega\mathbf{v}_{kl}.\vec{x}) }} } \ .$$ (III.A.4)

Nous explicitons le calcul analytique de (III.A.4) dans les cas suivants de valeurs des fonctions f et g.

1. La contribution d'un terme source localisé en un point x₀ de $\Omega$, soit pour $f = \delta_{x_0}$, est donnée par
   

   $$% b_{k,l} = \left\vert \null\,\vcenter{\openup\jot \let\\ =\@... ...0 \in \Omega_k \cr & 0 & \mbox{ sinon. } \cr \crcr}}\, \right.$$ (III.A.5)

   
   

2. La contribution d'une source surfacique sur un bord $B = \Gamma \cap \partial \Omega_k$ correspondant à l'excitation par une onde plane incidente de direction v₀, modélisée par
   
   $$g_{\vert\Gamma_k \subset B} = [(1+{{\mathit Q}}_k)({\partial... ...omega(1-{{\mathit Q}}_k)] e^{i\omega (\mathbf{v}_{0}.\vec{x})}$$
   
   
   où l'on suppose ${{\mathit Q}}_k = {{\mathit Q}}_{\vert\Gamma_k}$ constant, est
   

   $$% \null\,\vcenter{\openup\jot \let\\ =\@ \ialign{\strut\hfil... ...)\vec{\nu}_n\mathbf{v}_{0}+1-{{\mathit Q}}_k \ . \cr \crcr}}\,$$ (III.A.6)

   
   

3. On calcule la contribution d'une source volumique $f = \mu \omega^2 e^{i\omega(\mathbf{v}_{0}.\vec{x})}$. Sur un triangle de surface S et de barycentre $\vec{G}$, le calcul de l'intégrale de l'onde plane e^i<I>k<I>X est effectué annexe III.D.1.2. Alors, pour ${{\mathbf k}}= \omega(\mathbf{v}_0-\mathbf{v}_{kl})$, la contribution de la source volumique f est, pour un triangle $(\Omega_k)$,
   

   $$\null\,\vcenter{\openup1\jot \let\\ =\@ \ialign{\strut\hfil... ...ec{x}_3^k-\vec{x}_2^k)}{2} \ . \cr \crcr}}\, \right. \crcr}}\,$$ (III.A.7)

   
   

La mise en uvre informatique du calcul des termes explicités ci-dessus n'est pas directe. En effet, les intégrations des ondes planes sur un segment ou sur un triangle donnent des fonctions dont le domaine de définition pose problème.

1. Les termes matriciels D et C ne sont pas définis pour $\alpha=h_n^k$ ou $\alpha=h_{kj}$ nuls. En effet, la fonction $I_1(\alpha)$, intégrale sur un segment [x₁,x₂] de e^i<I>k<I>X, vaut
   
   $$I_1(\alpha) = L e^{i{{\mathbf k}}\vec{G}} \frac{\displaystyle \sin{\alpha}}{\displaystyle \alpha} \ ,$$
   
   
   où $\vec{G}$ est la position du barycentre du segment et L sa longueur, $\alpha$ est donné par $\alpha={{\mathbf k}}\frac{(x_1-x_2)}{2}$. Cette fonction, qui permet de calculer les termes matriciels, est définie sur ${\mathbb{R}}-\{0\}$. Elle se prolonge mathématiquement par continuité en 0, mais cette notion est ignorée par les calculateurs. Nous expliquons annexe III.D.2.1 comment calculer la fonction $I_1(\alpha)$ sur ${\mathbb{R}}$ ou ${\mathbb{C}}$ sur tout calculateur de façon à obtenir la meilleure précision possible. Le lecteur intéressé par les problèmes d'optimisation sur ``super-calculateurs'' consultera l'annexe III.D.3.1 pour l'implémentation sur machine à architecture vectorielle.
2. Le calcul du second membre dans le cas d'une source volumique de la forme $f = \mu \omega^2 e^{i\omega(\mathbf{v}_{0}.\vec{x})}$ n'est pas défini sur les droites $\alpha=0$, $\beta=0$ et $\alpha=\beta$ dans l'équation (III.A.7). Comme pour l'intégrale sur un segment $I_1(\alpha)$, nous expliquons, annexe III.D.2.2, comment calculer par ordinateur la fonction intégrale sur un triangle [x₁,x₂,x₃] de surface S de e^i<I>k<I>X, soit, pour $\vec{G}$ la position du barycentre du triangle et ${{\mathbf k}}= \omega(\mathbf{v}_0-\mathbf{v}_{kl})$,
   
   $$I_2(\alpha,\beta) = 2S e^{i{{\mathbf k}}\vec{G}} e^{ -\frac{... ...{\displaystyle \beta }}{\displaystyle 2i( \alpha-\beta )}) \ ,$$
   
   
   et comment optimiser le calcul sur machine vectorielle annexe III.D.3.2.

Le lecteur trouvera comment nous avons optimisé le problème suivant:

• garder une excellente précision de calcul,
• effectuer le calcul le plus rapidement possible.