Espace éléments finis (FiniteElementSpace/SubFiniteElementSpace)

Un maillage dit où sont les mailles ; un espace éléments finis dit comment on calcule dessus. Pour chaque zone du maillage, il fixe :

• une interpolation : les fonctions de forme, qui donnent une grandeur en tout point d’une maille à partir de ses valeurs aux nœuds ;
• une quadrature : les points de Gauss et leurs poids, où l’on évalue les intégrales.

C’est un agrégat de sous-espaces : un SubFiniteElementSpace par zone de maillage, dans le même ordre. Le maillage reste purement géométrique, et un même maillage peut porter plusieurs espaces.

Presque tout ce qui calcule s’appuie sur un espace : les physiques le reçoivent (model.heat_conduction(fes)), les champs aux points de Gauss posent leurs zones sur ses sous-espaces, et les opérateurs qui interpolent ou intègrent (element_field.interp_to_gauss, element_field.gradient, measure.integral…) le prennent en argument.

Construire un espace

pyrucast.FiniteElementSpace(mesh) pose sur chaque zone l’interpolation de Lagrange de son propre degré — "LAGRANGE1" sur des mailles linéaires, "LAGRANGE2" sur des quadratiques — et la quadrature "GAUSS" ; FiniteElementSpace.lagrange1(mesh) impose "LAGRANGE1" partout. Un argument interpolation= donné s’applique à toutes les zones à la fois, quadrature= de même, et FiniteElementSpace.with_choices(mesh, [(interpolation, quadrature), …]) en donne un couple par zone.

import pyrucast

# Un carré de côté 2 en deux triangles (zone 0), et son côté bas (zone 1).
c = pyrucast.Coords(2)
a = c.add_node([0.0, 0.0])
b = c.add_node([2.0, 0.0])
d = c.add_node([2.0, 2.0])
e = c.add_node([0.0, 2.0])
carre = pyrucast.Mesh(c, "TRI3")
carre.unit().add_cell([a, b, d])
carre.unit().add_cell([a, d, e])
bas = pyrucast.Mesh(c, "SEG2")
bas.unit().add_cell([a, b])
mesh = carre | bas

# Lagrange-1 et la quadrature de Gauss, sur chaque zone.
fes = pyrucast.FiniteElementSpace(mesh)
print(len(fes))  # 2
print(fes[0])  # SubFiniteElementSpace<TRI3, LAGRANGE1, GAUSS>: 3 Gauss point(s)
print(fes[1])  # SubFiniteElementSpace<SEG2, LAGRANGE1, GAUSS>: 2 Gauss point(s)

# Le même maillage, un point par maille : un autre espace.
centres = pyrucast.FiniteElementSpace(mesh, quadrature="REDUCED")
print(centres[0].gauss_count, centres[0].gauss_xi(0))  # 1 [0.333…, 0.333…]

# Un couple (interpolation, quadrature) par zone.
mixte = pyrucast.FiniteElementSpace.with_choices(
    mesh, [("LAGRANGE1", "GAUSS"), ("LAGRANGE1", "REDUCED")]
)
print([sub.quadrature for sub in mixte])  # ['GAUSS', 'REDUCED']

# Les zones portent maintenant un espace : leurs mailles sont figées.
print(carre.unit().is_sealed)  # True

Une interpolation de Lagrange doit avoir le degré des mailles : c’est ce que fait le choix par défaut, zone par zone, si bien qu’un maillage qui mêle zones linéaires et quadratiques se pose sans argument. Une interpolation donnée explicitement, elle, vaut pour toutes les zones. "REDUCED" sert à porter une valeur par maille, par exemple l’intégrale que rend element_field.cell_integral ; les coques l’emploient aussi, en interne, pour leur cisaillement transverse.

La construction est refusée pour :

• un maillage sans zone, ou une liste with_choices qui ne compte pas un couple par zone ;
• une zone POI1 : un point n’a pas d’élément de référence ;
• une interpolation que le type de maille n’admet pas ("LAGRANGE2" sur des TRI3, "HERMITE3" ailleurs que sur des SEG2) ;
• des mailles de dimension supérieure à celle des Coords (des TRI3 dans des Coords 1-D) ;
• un nom d’interpolation ou de quadrature inconnu.

Construire l’espace scelle les zones de maillage. Leurs mailles sont figées : add_cell y lève ensuite une erreur. Tout ce qui se pose sur l’espace — les valeurs d’un champ aux points de Gauss, les lignes d’une matrice, la coloration — est rangé par rang de maille, et une maille ajoutée décalerait tout. Les positions des nœuds, elles, restent libres (voir Déplacer le maillage). Le sceau est décrit avec le maillage.

Lire un espace

fes[i], fes.unit() et for sub in fes donnent les sous-espaces. Un sous-espace donne son élément de référence (points de Gauss, poids, fonctions de forme et leurs dérivées en \( \xi \)) et, pour chaque maille, les grandeurs physiques tirées de ses nœuds : le jacobien \( J \), sa mesure \( |J| \) et les dérivées \( \partial N_i/\partial x_a \). Un point s’y désigne par le rang de la maille dans la zone et le rang du point de Gauss dans la maille.

sub = fes[0]  # le sous-espace des triangles
print(sub.element_type, sub.interpolation, sub.quadrature)  # TRI3 LAGRANGE1 GAUSS
print(sub.ref_dim, sub.space_dim, sub.nodes_per_cell)  # 2 2 3
print(sub.cell_count, sub.gauss_count)  # 2 3

# Les tables de référence : le point de Gauss 0, son poids, N_i, ∂N_i/∂ξ_k.
print(sub.gauss_xi(0), sub.gauss_weight(0))  # [0.5, 0.0] 0.16666666666666666
print(sub.n_at_g(0))  # [0.5, 0.5, 0.0]
print(sub.dn_at_g(0))  # [-1.0, -1.0, 1.0, 0.0, 0.0, 1.0]

# La maille 0 au point 0 : J, |J| et ∂N_i/∂x_a, tirés de ses nœuds.
print(sub.jacobian(0, 0))  # [2.0, 2.0, 0.0, 2.0]
print(sub.det_jacobian(0, 0))  # 4.0
print(sub.dn_dx(0, 0))  # [-0.5, 0.0, 0.5, -0.5, 0.0, 0.5]

# Une maille vue seule : son aire est la somme des w_g |J_g|.
maille = fes.element(0, 1)  # ou sub[1]
print(sum(maille.gauss_weight(g) * maille.det_jacobian(g) for g in range(3)))  # 2.0
print(len(sub.elements()), len(fes.elements(1)))  # 2 1

# Chaque sous-espace rend sa zone : la même, pas une copie.
print(len(mesh | sub.submesh()))  # 2 : l'union ne la compte pas deux fois
# mesh() rend une zone par sous-espace, dans l'ordre.
print(len(fes.mesh()), len((fes | centres).mesh()))  # 2 4
print(len(fes | fes), len(fes | pyrucast.FiniteElementSpace(mesh)))  # 2 4

Une vue de maille (Element) donne les mêmes grandeurs sans le rang de la maille — jacobian(g), det_jacobian(g), dn_dx(g)… —, les mêmes dimensions (gauss_count, nodes_per_cell, ref_dim, space_dim), plus index, cell() et nodes().

Les tableaux sont rendus à plat, ligne par ligne :

Ces méthodes lisent les fonctions de forme géométriques, celles qui placent les points de la maille. Sous "LAGRANGE1" et "LAGRANGE2", ce sont aussi celles du champ ; sous "HERMITE3" et "MODEL_EMBEDDED", elles ne décrivent que la géométrie (voir plus bas).

Un espace par formulation, construit une fois. Chaque construction crée des sous-espaces neufs : deux FiniteElementSpace(mesh) sont deux espaces, même avec les mêmes choix, et les champs posés sur l’un ne se combinent pas point à point avec ceux de l’autre (voir Champ aux points de Gauss). On construit donc l’espace une fois et on le passe à tous les opérateurs. L’union | garde chaque sous-espace une seule fois et refuse deux espaces posés sur des Coords différents.

Déplacer le maillage

Les grandeurs physiques ne sont pas stockées : chaque appel les recalcule à partir des positions actuelles des nœuds. Déplacer un nœud suffit, sans rien reconstruire ni signaler :

# Allonger le côté bas : les deux zones voient la nouvelle position de b.
b.set_position([4.0, 0.0])
print(fes[0].det_jacobian(0, 0))  # 8.0, au lieu de 4.0
print(fes[1].det_jacobian(0, 0))  # 2.0, au lieu de 1.0

L’espace ne contrôle pas la forme des mailles, ni à la construction ni après un déplacement. det_jacobian rend une valeur absolue : une maille retournée (ses nœuds parcourus dans l’autre sens) mesure autant que la même à l’endroit. Une maille aplatie mesure 0, et dn_dx la refuse, faute de pouvoir inverser son jacobien.

Ce qui est mis en cache

Deux choses, gardées dans chaque sous-espace pour toute sa durée de vie :

• Les tables de référence, calculées à la construction : positions et poids des points de Gauss, \( N_i(\xi_g) \) et \( \partial N_i/\partial\xi_k(\xi_g) \) — sous "HERMITE3", en plus, les fonctions du champ et leurs dérivées secondes. Elles ne dépendent que du type de maille, de l’interpolation et de la quadrature : calculées une fois, elles servent ensuite à toutes les mailles et à tous les opérateurs. Leur taille ne dépend pas du nombre de mailles, quelques milliers de nombres au plus (un HEX27 à 27 points).
• La coloration des mailles, calculée au premier usage : un partage des mailles en couleurs, deux mailles de même couleur n’ayant aucun nœud commun. Les mailles d’une couleur écrivent alors en parallèle sans se gêner, dans l’assemblage des matrices (matrix.stiffness, matrix.mass, matrix.tangent…) et dans les opérateurs qui accumulent aux nœuds (node_field.internal_forces, node_field.external_forces, node_field.divergence, node_field.sum_to_nodes). La calculer demande un passage sur toute la connectivité ; la garder évite de le refaire à chaque appel, alors qu’un calcul non linéaire réassemble à chaque itération. Elle ne dépend que de la connectivité, que le sceau fige.

Ce qui n’est pas gardé : \( J \), \( |J| \), \( \partial N_i/\partial x_a \) et la position réelle des points de Gauss. Les garder coûterait une matrice par point de Gauss — de la mémoire en proportion du maillage — et les rendrait fausses au premier nœud déplacé. Les recalculer coûte une somme sur les nœuds de la maille, que chaque assemblage fait une fois par point de Gauss avant de s’en servir pour tous les termes de la matrice élémentaire.

Sauvegarde

Un espace se sauvegarde avec ses zones de maillage et ses Coords. À la lecture, le sceau des zones est conservé, et les tables de référence sont relues telles quelles. La coloration n’est pas écrite : elle se recalcule au premier usage. Un espace et les champs posés dessus, sauvegardés dans un même fichier, restent liés à la lecture : les champs relus sont posés sur l’espace relu.

Théorie

Élément de référence

Chaque type de maille a un élément de référence \( \hat K \), décrit dans un repère \( \xi \), et une numérotation locale de ses nœuds :

Les types quadratiques reprennent le repère et les sommets de leur parent linéaire, puis ajoutent les nœuds de milieu d’arête, dans l’ordre des arêtes de la convention VTK :

QUA8, HEX20 et PENTA15 sont sérendipité (nœuds d’arête seulement) ; SEG3, TRI6, TET10, QUA9 et HEX27 sont des Lagrange complets (QUA9/HEX27 : produits tensoriels \( Q_2 \), avec nœuds de face et central).

Les mailleurs produisent leurs mailles dans ces conventions, l’export VTK les écrit telles quelles et la lecture gmsh y remet ses nœuds milieux. La fiche de chaque type, avec ses fonctions de forme et sa quadrature, est dans le catalogue des éléments ; cette section décrit ce qui leur est commun.

Transformation isoparamétrique

Soit une maille à \( n \) nœuds, de coordonnées physiques \( \mathbf{x}_1, \dots, \mathbf{x}_n \in \mathbb{R}^{d_s} \), avec \( d_s \) la dimension des Coords. La transformation géométrique \( \chi : \hat{K} \to K \), qui envoie l’élément de référence sur la maille, est interpolée par les fonctions de forme :

\[ \mathbf{x}(\xi) = \chi(\xi) = \sum_{i=1}^{n} N_i(\xi)\, \mathbf{x}_i \]

C’est l’hypothèse isoparamétrique : la géométrie s’interpole comme un champ. En Lagrange-1, \( \chi \) est affine sur les simplexes (SEG2, TRI3, TET4), bilinéaire sur QUA4 et trilinéaire sur HEX8.

Fonctions de forme de Lagrange

Une fonction de forme de Lagrange vérifie la propriété de Kronecker : \( N_i(\xi_j) = \delta_{ij} \) aux nœuds de référence. En Lagrange-1 :

SEG2 (\( \xi \in [-1, +1] \)) : \[ N_1(\xi) = \frac{1 - \xi}{2}, \qquad N_2(\xi) = \frac{1 + \xi}{2} \]

TRI3 (coordonnées barycentriques sur le simplexe \( \xi + \eta \le 1 \)) : \[ N_1 = 1 - \xi - \eta, \qquad N_2 = \xi, \qquad N_3 = \eta \]

QUA4 (bilinéaire sur \( [-1, +1]^2 \)) : \[ N_1 = \tfrac{1}{4}(1-\xi)(1-\eta), \quad N_2 = \tfrac{1}{4}(1+\xi)(1-\eta), \quad N_3 = \tfrac{1}{4}(1+\xi)(1+\eta), \quad N_4 = \tfrac{1}{4}(1-\xi)(1+\eta) \]

TET4 (coordonnées barycentriques sur le simplexe \( \xi + \eta + \zeta \le 1 \)) : \[ N_1 = 1 - \xi - \eta - \zeta, \quad N_2 = \xi, \quad N_3 = \eta, \quad N_4 = \zeta \]

HEX8 (trilinéaire sur \( [-1, +1]^3 \)). Pour le nœud \( i \) de coordonnées de référence \( (\xi_i, \eta_i, \zeta_i) \in \{-1, +1\}^3 \), \[ N_i(\xi, \eta, \zeta) = \tfrac{1}{8}\, (1 + \xi_i\, \xi)\, (1 + \eta_i\, \eta)\, (1 + \zeta_i\, \zeta) \]

PENTA6 (prisme, produit d’un TRI3 par un SEG2 linéaire sur \( \zeta \in [0, 1] \)). Avec les coordonnées barycentriques du triangle \( L_1 = 1 - \xi - \eta \), \( L_2 = \xi \), \( L_3 = \eta \), \[ N_j(\xi, \eta, \zeta) = L_j\,(1 - \zeta) \quad (j = 1, 2, 3), \qquad N_{j+3}(\xi, \eta, \zeta) = L_j\,\zeta \quad (j = 1, 2, 3) \]

PYRA5 a des fonctions rationnelles, qui recollent la base carrée au sommet : voir sa fiche.

En Lagrange-2, avec les coordonnées barycentriques \( L_i \) pour les simplexes :

• SEG3 : \( N_0 = \tfrac12\xi(\xi-1),\; N_1 = \tfrac12\xi(\xi+1),\; N_2 = 1-\xi^2 \) ;
• TRI6 / TET10 : sommets \( L_i(2L_i-1) \), milieux d’arête \( 4 L_a L_b \) ;
• QUA9 : produit tensoriel \( N_i(\xi,\eta) = \ell_{a}(\xi)\,\ell_{b}(\eta) \) des trois fonctions de Lagrange 1-D \( \ell_{-}(t)=\tfrac12 t(t-1),\; \ell_0(t)=1-t^2,\; \ell_{+}(t)=\tfrac12 t(t+1) \) ;
• HEX27 : produit tensoriel 3-D \( N_i(\xi,\eta,\zeta) = \ell_{a}(\xi)\,\ell_{b}(\eta)\,\ell_{c}(\zeta) \) des mêmes fonctions ;
• QUA8 (sérendipité) : sommet \( \tfrac14(1+\xi_i\xi)(1+\eta_i\eta)(\xi_i\xi+\eta_i\eta-1) \) ; milieu \( \tfrac12(1-\xi^2)(1+\eta_i\eta) \) ou \( \tfrac12(1+\xi_i\xi)(1-\eta^2) \) selon l’arête ;
• HEX20 (sérendipité) : sommet \( \tfrac18(1+\xi_i\xi)(1+\eta_i\eta)(1+\zeta_i\zeta)(\xi_i\xi+\eta_i\eta+\zeta_i\zeta-2) \) ; milieu \( \tfrac14(1-\xi^2)(1+\eta_i\eta)(1+\zeta_i\zeta) \) et ses permutations ;
• PENTA15 (sérendipité) : produit du TRI6 par un facteur quadratique en \( \zeta \), corrigé aux sommets.

Toutes forment une partition de l’unité, \( \sum_i N_i(\xi) = 1 \) en tout \( \xi \) : un champ constant est reproduit exactement. Les tests le vérifient, avec Kronecker, et recoupent les dérivées des types quadratiques par différences finies.

Dérivées de référence

Les dérivées de référence \( \partial N_i / \partial \xi_k \) s’obtiennent en dérivant les formules. En Lagrange-1, elles sont constantes sur SEG2, TRI3 et TET4 ; sur les autres types, et sur tous les quadratiques, elles sont polynomiales en \( \xi \). La partition de l’unité dérivée donne \( \sum_i \partial N_i / \partial \xi_k = 0 \) dans chaque direction \( k \), identité elle aussi testée.

Fonctions de forme d’Hermite (HERMITE3)

Les fonctions de Lagrange sont C⁰ : le champ est continu d’une maille à la suivante, sa dérivée ne l’est pas. Cela suffit aux équations du second ordre. Une poutre d’Euler-Bernoulli obéit à \( (EIw'')'' = q \), du quatrième ordre, dont la forme faible exige un champ C¹ : flèche et pente continues.

C’est ce que fournit "HERMITE3", sur SEG2 :

\[ \begin{aligned} H_1 &= \tfrac14(2 - 3\xi + \xi^3), &\qquad H_2 &= \tfrac14(1 - \xi - \xi^2 + \xi^3),\\ H_3 &= \tfrac14(2 + 3\xi - \xi^3), &\qquad H_4 &= \tfrac14(-1 - \xi + \xi^2 + \xi^3). \end{aligned} \]

La propriété de Kronecker y porte sur deux grandeurs : à chaque extrémité, une fonction vaut 1 et a une pente nulle, l’autre vaut 0 et a une pente 1, et les deux fonctions de l’extrémité opposée s’y annulent avec leur pente. D’où l’ordre des degrés de liberté, \( [w_A,\ w'_A,\ w_B,\ w'_B] \) : une valeur et une pente par nœud.

Quatre fonctions pour deux nœuds. Le nombre de fonctions de forme n’est plus le nombre de nœuds : les tables du champ sont dimensionnées sur le nombre de fonctions, celles de la géométrie sur le nombre de nœuds.

L’élément est sous-paramétrique. Un SEG2 reste un segment droit, interpolé en Lagrange-1 quel que soit le champ qu’il porte. Le sous-espace tabule donc deux bases : la géométrique, qui donne le jacobien \( \partial x/\partial\xi \) et que lisent n_at_g, dn_at_g et dn_dx, et celle du champ, qui interpole l’inconnue. Sous Lagrange, les deux coïncident et seule la géométrique est stockée.

Dérivées secondes. "HERMITE3" est la seule interpolation à tabuler \( \partial^2 N_i/\partial\xi^2 \) : la courbure est une grandeur primaire d’un élément C¹. Sur SEG2, ces dérivées sont linéaires en \( \xi \) : la courbure varie dans la maille. Le passage au physique est une règle de chaîne, sans terme en \( \partial J/\partial\xi \) puisqu’un segment a un jacobien constant :

\[ \frac{\partial^2 N_i}{\partial x^2} = \frac{1}{J^2}\, \frac{\partial^2 N_i}{\partial \xi^2}, \qquad J = \frac{L}{2}. \]

La pente de référence n’est pas la rotation. \( H_2 \) et \( H_4 \) portent une pente \( \partial w/\partial\xi \), alors que le degré de liberté d’une poutre est \( \theta = \partial w/\partial x \). Les deux diffèrent du jacobien, \( \partial w/\partial\xi = J\,\partial w/\partial x \). Ce facteur n’est pas dans les fonctions de forme : c’est un passage référence → physique, appliqué là où sont appliqués tous les autres jacobiens. Avec un \( L \) dans la fonction de forme, la base dépendrait de la maille et ne pourrait plus être tabulée une fois par type.

Ce qui le vérifie. La raideur d’Euler-Bernoulli est l’oracle capable de mettre cette base en défaut :

\[ K = \int_0^L EI \left(\frac{\partial^2 N}{\partial x^2}\right)^{\!\top} \frac{\partial^2 N}{\partial x^2}\, dx = \frac{EI}{L^3} \begin{bmatrix} 12 & 6L & -12 & 6L \\ 6L & 4L^2 & -6L & 2L^2 \\ -12 & -6L & 12 & -6L \\ 6L & 2L^2 & -6L & 4L^2 \end{bmatrix}. \]

tests/hermite.rs intègre le membre de gauche depuis la base tabulée et le compare à la forme fermée, à la précision machine. Une erreur dans les fonctions, dans leurs dérivées secondes ou dans le facteur jacobien des pentes apparaît dans cette comparaison. L’intégrande étant quadratique en \( \xi \), les deux points de Gauss du SEG2 l’intègrent exactement.

Pas de base du tout : MODEL_EMBEDDED

Un élément dont la matrice élémentaire est une forme fermée n’évalue jamais de fonction de forme : son intégrale a été faite une fois, à la main, et seul le résultat figure dans le code. La poutre de Timoshenko exacte en est un : ses fonctions de forme dépendent du matériau, à travers \( \Phi = 12EI/(G A_s L^2) \), et aucun espace ne peut les tabuler.

Y déclarer une interpolation de Lagrange affirmerait quelque chose que la physique n’utilise pas, et qui serait faux : la forme fermée vient de fonctions cubiques et quadratiques, pas linéaires. "MODEL_EMBEDDED" déclare qu’il n’y a pas de base du champ :

La géométrie reste entièrement définie — coordonnées, jacobien, mesure — et tout le reste fonctionne comme sur un autre espace : champ matériau, champs aux points de Gauss, assemblage, coloration. Seule l’interpolation de l’inconnue n’appartient pas à l’espace, et les opérateurs qui en ont besoin lèvent une erreur qui le dit, plutôt que de se rabattre sur la base géométrique :

Le repli serait plausible, et c’est ce qui le rendrait dangereux : sur une poutre, la charge répartie cohérente vaut \( qL/2 \) et \( qL^2/12 \) en moment ; une base linéaire donnerait \( qL/2 \) sans moment.

La visualisation, elle, reste sur la base géométrique : ce qu’elle colorie est ce qu’elle dessine, un segment droit ou une facette plane. C’est une image, pas une valeur calculée.

La formulation fournit alors elle-même ce que l’espace ne peut plus interpoler : element_field.beam_deformation évalue les déformations de la poutre depuis les fonctions de forme de l’élément, et reçoit pour cela le champ matériau, dont \( \Phi \) dépend.

Quadrature de Gauss

Pour intégrer une fonction \( f \) sur une maille, on se ramène à l’élément de référence par le changement de variables \( \mathbf{x} = \chi(\xi) \) :

\[ \int_K f(\mathbf{x})\, d\mathbf{x} = \int_{\hat{K}} f(\chi(\xi))\, |J(\xi)|\, d\xi \approx \sum_{g=1}^{n_g} w_g\, f(\chi(\xi_g))\, |J(\xi_g)| \]

avec \( |J| \) la mesure du jacobien, définie plus bas. Le couple \( (\xi_g, w_g) \) est la règle de quadrature. La règle "GAUSS" de chaque type intègre exactement sa matrice de masse sur une maille à géométrie droite :

La règle "REDUCED" n’a qu’un point, au centre de gravité de l’élément de référence, de poids égal à sa mesure : elle intègre exactement les fonctions affines. D’où viennent ces règles, et ce que chacune intègre exactement, est expliqué avec les éléments finis.

La somme des poids vaut la mesure de l’élément de référence : 2 pour SEG2/SEG3, 1/2 pour TRI3/TRI6, 4 pour QUA4/QUA8/QUA9, 1/6 pour TET4/TET10, 4/3 pour PYRA5, 1/2 pour PENTA6/PENTA15, 8 pour HEX8/HEX20/HEX27. Les tests le vérifient pour tous les types, et vérifient l’exactitude des règles TRI6 et TET10 en intégrant des monômes.

Jacobien et grandeurs physiques

Le jacobien de la transformation \( \chi \) est la matrice des dérivées \[ J_{a,k}(\xi) = \frac{\partial x_a}{\partial \xi_k} = \sum_{i=1}^{n} \mathbf{x}_{i,a}\, \frac{\partial N_i}{\partial \xi_k}(\xi) \quad \text{de taille } d_s \times d_r \] où \( a \) parcourt les \( d_s \) directions physiques et \( k \) les \( d_r = \dim \hat K \) directions de référence.

Cas plein : \( d_s = d_r \). Quand la maille a la dimension des Coords (un TRI3 dans des Coords 2-D), \( J \) est carrée. Son déterminant mesure la dilatation locale du volume, et c’est sa valeur absolue qui entre dans l’intégration : det_jacobian rend \( |\det J| \). Les dérivées physiques des fonctions de forme passent par l’inverse de \( J \) : \[ \frac{\partial N_i}{\partial x_a} = \sum_{k=1}^{d_r} (J^{-1})_{k, a} \, \frac{\partial N_i}{\partial \xi_k} \quad \Longleftrightarrow \quad \nabla_x N_i = J^{-T}\, \nabla_\xi N_i \]

Cas plongé : \( d_s > d_r \). Une zone peut être plongée dans un espace de dimension supérieure : des SEG2 dans des Coords 2-D ou 3-D (un contour, une barre), des TRI3 dans des Coords 3-D (une surface, une coque). \( J \) est alors rectangulaire et n’a plus de déterminant, mais la métrique \[ G(\xi) = J(\xi)^T\, J(\xi) \quad \text{de taille } d_r \times d_r \] est symétrique définie positive tant que la maille n’est pas dégénérée. L’élément de mesure devient \[ d\mu = \sqrt{\det G}\, d\xi \] qui joue le rôle de \( |J| \) dans la quadrature : det_jacobian rend ce \( \sqrt{\det G} \). Le gradient tangent d’un champ (la projection du gradient sur la maille) passe par la pseudo-inverse : \[ \nabla_s N_i = J\, G^{-1}\, \nabla_\xi N_i \] dn_dx rend ces composantes, dans le repère des Coords. Pour \( d_s = d_r \), ces formules redonnent celles du cas plein (\( J G^{-1} = J^{-T} \)), et c’est la seule forme calculée.

Dimensions admises :

\( d_s < d_r \) n’a pas de sens (pas de maille 2-D dans un espace 1-D) : la construction le refuse.

Maille dégénérée. dn_dx doit inverser \( G \). Il refuse une maille dont \( \det G \) tombe sous \( 10^{-12} \) fois l’échelle de \( G \) (son plus grand terme, à la puissance \( d_r \)). Le seuil est relatif : une maille d’un demi-millimètre décrite en mètres a un \( \det G \) de l’ordre de \( 10^{-21} \) sans être dégénérée pour autant, et l’unité de longueur choisie ne change rien au verdict.

Axisymétrie. Sur des Coords axisymétriques, le plan du maillage est le plan méridien \( (r, z) \), et une maille y représente un anneau. det_jacobian reste le \( |J| \) de la section : le facteur \( 2\pi r \) entre dans le poids d’intégration, appliqué par les opérateurs qui intègrent, de sorte que raideur, masse, flux, volumes et efforts intérieurs portent sur l’anneau entier.

Côté Rust

SubFiniteElementSpace détient :

• le Handle<SubMesh> de sa zone, scellé à la construction (mesh::seal) ;
• son Interpolation et sa QuadratureRule, la dimension des Coords et leur caractère axisymétrique, lus une fois à la construction (un Coords n’en change pas) ;
• les tables de référence, en vecteurs plats : gauss_xi, gauss_w, n_at_g, dn_at_g pour la géométrie, et field_n_at_g, field_dn_at_g, field_d2n_at_g pour le champ, vides quand le champ suit la géométrie (toute interpolation de Lagrange) ;
• la coloration, un OnceLock<Vec<Vec<usize>>> marqué #[serde(skip)].

FiniteElementSpace n’est qu’un Vec<Handle<SubFiniteElementSpace>> : il ne garde pas de Mesh. mesh() en reconstruit un avec la zone de chaque sous-espace, partagée, dans l’ordre, et check_push refuse un sous-espace sur d’autres Coords que le premier.

Les constructeurs sont FiniteElementSpace::lagrange1(&mesh), FiniteElementSpace::new(&mesh, interpolation) (Gauss partout) et FiniteElementSpace::with(&mesh, &[(interpolation, quadrature), …]), un couple par zone ; SubFiniteElementSpace::new(submesh, interpolation, quadrature) construit un sous-espace seul.

#[test]
fn le_constructeur_par_defaut_pose_lagrange1_et_gauss() -> Result<()> {
    let coords = Handle::new(Coords::new(2)?);
    let a = Node::create_in(coords.clone(), &[0.0, 0.0])?;
    let b = Node::create_in(coords.clone(), &[2.0, 0.0])?;
    let c = Node::create_in(coords.clone(), &[0.0, 2.0])?;

    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::TRI3));
    mesh.add_cell(&[a.id(), b.id(), c.id()])?;

    let fes = FiniteElementSpace::lagrange1(&mesh)?;
    let sub = fes.get(0)?;
    let s = sub.read();
    assert_eq!(s.gauss_count(), 3);
    // Le triangle (0,0), (2,0), (0,2) a |J| = 4 partout : la transformation
    // est affine, et |J| vaut son aire (2) sur celle du triangle de référence (1/2).
    for g in 0..s.gauss_count() {
        let dj = s.det_jacobian(0, g)?;
        assert!((dj - 4.0).abs() < 1e-12);
    }
    Ok(())
}

#[test]
fn le_constructeur_explicite_choisit_par_sous_maillage() -> Result<()> {
    let (_, mesh, _) = triangle()?;

    let fes =
        FiniteElementSpace::with(&mesh, &[(Interpolation::Lagrange1, QuadratureRule::Gauss)])?;

    assert_eq!(fes.len(), 1);
    Ok(())
}

Les évaluations sont les mêmes qu’en Python, sur le guard du sous-espace ; les tables de référence sont prêtées en &[f64], les grandeurs physiques rendues en Vec<f64> :

#[test]
fn evaluer_les_grandeurs_sur_une_cellule() -> Result<()> {
    let (_, _, fes) = triangle()?;
    let sub = fes.get(0)?;

    let s = sub.read();
    for cell_idx in 0..s.cell_count() {
        for g in 0..s.gauss_count() {
            let n = s.n_at_g(g)?; // N_i(ξ_g)
            let dn = s.dn_at_g(g)?; // ∂N_i/∂ξ_k(ξ_g)
            let jac = s.jacobian(cell_idx, g)?;
            let det_j = s.det_jacobian(cell_idx, g)?;
            let dn_dx = s.dn_dx(cell_idx, g)?;
            // … de quoi écrire la matrice élémentaire …
            let _ = (n, dn, jac, det_j, dn_dx);
        }
    }
    Ok(())
}

Le reste de l’API sert les physiques :

• La base du champ. shape_count(), field_n_at_g(g), field_dn_at_g(g) et field_d2n_at_g(g) lisent la base qui interpole l’inconnue. Sous Lagrange, elles rendent les tables géométriques ; sous ModelEmbedded, elles lèvent une erreur qui nomme la situation ; field_d2n_at_g n’existe que sous Hermite3.
• La coloration. coloring(compute) rend la coloration, en la calculant au premier appel par la fermeture fournie : c’est l’assemblage (models::kernel) qui choisit les clés de conflit, ici les nœuds des mailles, et la couche conteneur n’en dépend pas.
• Les pilotes de models::kernel. Un opérateur ne passe pas par jacobian ou dn_dx, qui allouent leur résultat. Il copie une fois, par zone et par appel, les tables de référence du sous-espace, puis parcourt les mailles en parallèle. Chaque maille reçoit un CellGeom qui lit les coordonnées de ses nœuds une seule fois, au premier besoin, et calcule \( \partial N_i/\partial x_a \) et \( |J|\,w_g \) dans des tampons sur la pile (neuf nombres au plus pour \( J \), \( G \) ou leur inverse). CellGeom::det_j_w est le seul endroit où se décide la mesure, facteur \( 2\pi r \) compris, et CellGeom::field_n_at_g celui où les pentes d’Hermite sont multipliées par \( J \).