Opérateurs sur les champs

Les modules ops::node_field, ops::element_field, ops::coords, ops::measure et ops::field dérivent et transforment les champs — chacun nommé d’après ce qu’il produit : coordonnées, restriction, projection sur un autre maillage, fusion, et dérivations géométriques vers les points de Gauss. Ces opérateurs croisent des conteneurs (maillage + champ, espace EF

• champ) — ce sont donc des fonctions libres. L’arithmétique scalaire et par composante, elle, reste sur les types de champ (cf. Champ).

Les autres thèmes d’opérateurs ont leur propre page : construction du matériau (Construction), Assemblage (dont le chargement réparti flux), Comportement, Solveur.

Coordonnées et déplacement

positions est le pont géométrie → champ : on en tire un NodeField qu’on peut tracer, dériver, ou réinjecter après calcul (displace pour passer à la configuration déformée). En grandes transformations, la configuration courante vit dans un jeu à part du même Coords, posé par displace_from et rendu actif le temps d’un calcul par in_config : voir Grandes transformations.

Restriction, fusion, consolidation

La consolidation d’un NodeField est exactement la finalisation de l’union | : après déduplication par handle, les zones définies sur le même SubMesh deviennent une seule zone portant l’union de leurs composantes — une composante définie par plusieurs zones doit y avoir la même valeur partout (sinon erreur). Une vérification inter-supports finale impose qu’un nœud partagé par des zones de supports différents s’accorde sur toute composante commune. Le même consolidate accepte un ElementField (opération element_field.consolidate) : les sous-champs d’une même FiniteElementSpace fusionnent en une zone portant l’union de leurs composantes — utile pour réunir des zones matériau bâties par physique sur une fespace partagée (k thermique + E/nu/alpha mécanique) en un champ matériau unique lu par chaque physique.

Projection sur un autre maillage

project(field, mesh, target, tol=1e-6) (aussi field.project(mesh, target)) porte un champ nodal aux nœuds de target en l’interpolant par les mailles de mesh : un nœud de target situé dans une maille de mesh reçoit Σ Nᵢ(ξ) f(nœudᵢ), les fonctions de forme de cette maille en ses coordonnées de référence (une maille quadratique interpole quadratiquement), pour toutes les composantes. C’est le PROI de Cast3M sur une géométrie. Un champ nodal ne porte pas de mailles, d’où l’argument mesh, qui doit être sur le même Coords que le champ ; target peut vivre dans un autre Coords, pour porter un champ sur une pièce remaillée.

Un nœud de target hors de mesh — au-delà de tol hors de l’élément de référence de chaque maille, un jeu sans unité — reçoit 0, comme un champ nodal lit un nœud qu’il ne définit pas : on projette ainsi un patch sur un maillage plus grand. Le résultat a une zone par sous-maillage de target, comme positions. Les nœuds de target sont localisés une seule fois.

# La température de la plaque portée aux nœuds d'un autre maillage — ici d'un
# autre Coords, comme une pièce remaillée — par les fonctions de forme des
# mailles de la plaque. Un nœud hors de la plaque reçoit 0 : on projette aussi
# un patch sur un maillage plus grand.
autre = pyrucast.Coords(2)
cibles = [autre.add_node(p) for p in ([0.5, 0.5], [1.5, 0.25], [3.0, 0.5])]
nuage = pyrucast.mesh.poi1_from_nodes(cibles)
portee = pyrucast.node_field.project(temperature, mesh, nuage)  # tol=1e-6
# La même, en méthode : temperature.project(mesh, nuage)

Un champ aux points de Gauss se projette aussi. Avec un ElementField, chaque nœud de target reçoit la valeur que donnent les seuls points de Gauss de la maille qui le contient : leur ajustement aux moindres carrés sur cette maille — celui-là même qui colore la maille sur un tracé —, évalué au point. Rien n’est moyenné entre mailles : deux voisines peuvent donner deux valeurs différentes sur leur face commune, et une suite de points qui passe de l’une à l’autre montre ce saut. Un point posé exactement sur une face partagée est lu dans l’une des deux ; pour suivre un champ le long d’une ligne, on place les points entre les nœuds. mesh est alors le maillage du champ : chacun de ses sous-maillages doit porter une zone du champ.

# A field at the Gauss points read at points: each point takes the value its
# own cell's Gauss points give there — nothing averaged with the neighbour.
de_part_et_d_autre = [autre.add_node(p) for p in ([0.9, 0.5], [1.1, 0.5])]
deux_points = pyrucast.mesh.poi1_from_nodes(de_part_et_d_autre)
lue = pyrucast.node_field.project(par_maille, mesh, deux_points)

Les valeurs relevées sur un nuage de points se tracent ensuite le long de la ligne qu’ils forment : voir Le long d’une ligne.

Bande de valeurs (ge / gt / le / lt)

select et mask partagent la même bande de valeurs, fixée par quatre bornes de comparaison qui reprennent une pour une les opérateurs Python :

On donne au plus une borne basse (ge ou gt) et au plus une borne haute (le ou lt), et au moins une borne en tout ; une borne absente laisse ce côté ouvert. Erreur si aucune borne, ou si la borne basse dépasse la haute.

Sélection par valeur

select extrait, zone par zone, la partie du support d’un champ dont les valeurs tombent dans la bande. C’est un filtre par valeur qui renvoie un Mesh — un sous-maillage par zone traitée (les zones restent séparées, rien n’est moyenné ni fusionné).

• Type de champ. Sur un NodeField / SubNodeField, on sélectionne les nœuds : chaque zone donne un sous-maillage POI1 des nœuds retenus. Sur un ElementField / SubElementField, on sélectionne les cellules : chaque zone donne un sous-maillage de son propre type d’élément, et une cellule n’est retenue que si tous ses points de Gauss passent (la bande doit tenir tout le long de la cellule).
• Composantes. components=None teste toutes les composantes de chaque zone. Une liste components ne teste que ces composantes, et seulement sur les zones qui les portent toutes — une zone à laquelle il manque une composante demandée est ignorée (aucun sous-maillage produit).
• Combinaison (ET). Quand plusieurs composantes sont testées, elles sont combinées en ET : un nœud / une cellule n’est retenu que si chaque composante testée est dans la bande.

# Nodes whose temperature lies between 20 and 80 °C (inclusive bounds).
chauds = pyrucast.mesh.select(temperature, ge=20.0, le=80.0)

# Cellules dont la contrainte de von Mises dépasse un seuil (borne basse seule).
critiques = pyrucast.mesh.select(sigma, ge=250e6, components=["vm"])

Masque par valeur

mask garde la structure exacte du champ (mêmes zones, même support, mêmes composantes) et se contente de réécrire les valeurs : 1.0 là où la bande tient, 0.0 sinon — composante par composante (le MASQUE de Cast3M). Le résultat est donc du même type que l’entrée et se multiplie terme à terme avec elle. Un NodeField est masqué par nœud, un ElementField par point de Gauss.

• Pas de ET entre composantes (contrairement à select) : chaque valeur est testée pour elle-même.
• Composantes. components=None teste toutes les composantes. Une liste components ne teste que celles-ci ; les autres restent à 1.0 (neutre pour le produit), et une zone à laquelle il manque une composante demandée reste tout à 1.0.

# Resets a field's negative values to zero, component by component.
positif = champ * champ.mask(ge=0.0)

# Sugar: comparisons build a mask directly.
positif = champ * (champ >= 0.0)  # the same thing
chauds = temperature > 80.0  # NodeField 0/1

Les opérateurs >=, >, <=, < sur un champ (NodeField, SubNodeField, ElementField, SubElementField) renvoient le masque correspondant contre le scalaire de droite. == / != gardent leur sens Python habituel (identité).

Exemple complet et exécutable : examples/field_mask.py (lancer avec python examples/field_mask.py après maturin develop).

Extraction et renommage de composantes (EXCO)

Deux opérateurs travaillent sur le jeu de composantes d’un champ, sans toucher au support ni aux valeurs — l’équivalent de EXCO de Cast3M. Tous deux acceptent les quatre saveurs (NodeField, SubNodeField, ElementField, SubElementField) et renvoient la même saveur.

filter_components traite chaque zone indépendamment :

• une zone ne portant aucune des composantes demandées est abandonnée ;
• une zone ne portant que des composantes demandées (rien à retirer) voit son sous-champ partagé tel quel (handle copié, pas de duplication) ;
• une zone mixte est reconstruite sur le même support avec les seules composantes demandées, dans son propre ordre.

components peut être un sur-ensemble des composantes du champ (les noms absents sont ignorés) : passer model.primal_vars() à un résultat de solve pour en retirer les inconnues duales (multiplicateurs de Lagrange) est l’usage visé. Erreur si aucune zone ne porte l’une des composantes demandées.

rename_component laisse inchangée (handle partagé) toute zone ne portant pas old. Erreur si aucune zone ne porte old, ou si une zone concernée a déjà une composante nommée new.

# Removes the Lagrange multipliers from a solve result.
u = solution.filter_components(model.primal_vars())

# Renames a component before exporting.
export = u.rename_component("u_x", "DX")

Sucre d’indexation (façon pandas/numpy) : sur les quatre saveurs (NodeField, SubNodeField, ElementField, SubElementField), une clé chaîne ou liste de chaînes appelle filter_components et renvoie la même saveur. Les autres clés gardent leur sens : int/slice → accès aux zones sur un agrégat ; le tuple d’accès à une valeur sur un sous-champ (sub[node, "UX"], sub[cell, gauss, "E"]) est inchangé.

ux = champ["u_x"]  # == filter_components(champ, "u_x")
depl = champ[["u_x", "u_y"]]  # == filter_components(champ, ["u_x", "u_y"])
zone = champ[0]  # inchangé : la zone (SubNodeField)
val = champ[0][node, "u_x"]  # unchanged: the value at the node

L’accesseur champ.components() (présent sur les quatre saveurs) donne la liste des composantes, d’où l’idiome « reprojeter u1 sur les composantes de u2 » :

u = u1[u2.components()]  # u1 cut down to u2's set of components

Côté Rust, filter_components / select_components acceptent indifféremment un &str, un tableau ["u_x", "u_y"] ou un Vec<String> (trait IntoComponentNames) — donc field.filter_components(model.primal_vars()) passe directement.

Dérivation géométrique (vers les points de Gauss)

Ces opérateurs ne dépendent que de l’espace EF et du champ — aucune physique. Ils produisent l’ElementField que le comportement (integrate_behavior) consomme ensuite ; au niveau d’un modèle, strain choisit pour chaque loi celui qui la nourrit. Ils partagent tous le même moteur parallèle : le driver nodal_pointwise (déterministe bit-à-bit, cf. Parallélisme), pendant nodal de element_pointwise.

interp_to_gauss(field, fespace) → ElementField

Interpole un champ nodal vers les points de Gauss (valeurs, pas dérivées) :

\[ f(\xi_g) = \sum_i f_i\, N_i(\xi_g). \]

Le résultat porte les mêmes composantes que l’entrée, une valeur par (cellule, point de Gauss). C’est le pendant « valeurs » de gradient (direction nœuds → Gauss du CHAN de Cast3M) : typiquement pour porter une température nodale aux points de Gauss avant thermal_strain.

gradient(field, fespace) → ElementField

Gradient d’un champ nodal aux points de Gauss, cellule par cellule :

\[ \nabla f = \sum_i f_i\, \nabla N_i \quad \text{évalué en chaque } \xi_g. \]

Une composante de sortie <comp>_<axe> par couple (composante d’entrée, axe).

deformation(u, fespace, measure="linear") → ElementField

Déformation d’un champ de déplacement aux points de Gauss, dans la mesure choisie. u doit porter exactement space_dim composantes (déplacement selon x, y, z). On note \(H = \nabla u\) et \(F = I + H\) :

Chaque mesure a son préfixe : on ne peut pas en lire une à la place d’une autre. Une mesure symétrique produit une composante <préfixe>_<ai><aj> par entrée indépendante i ≤ j. Le gradient F, qui n’est pas symétrique, les produit toutes (9 en 3D, 4 en 2D). En 2D, la direction hors plan n’est pas déformée (\(F_{zz} = 1\)) et ses entrées ne sont pas produites.

La mesure linéarisée est en convention tenseur (eps_xy = ½(∂u_x/∂y + ∂u_y/∂x), pas le cisaillement ingénieur γ). C’est l’entrée du comportement de l’élasticité. Les mesures finies valent zéro (ou \(I\) pour C et b) sous une rotation rigide finie, contrairement à eps.

Sur un sous-espace axisymétrique, une composante <préfixe>_zz est ajoutée. C’est l’entrée orthoradiale, construite sur \(H_{\theta\theta} = u_r / r\), que le gradient méridien ne peut pas exprimer : eps_zz = u_r / r, F_zz = 1 + u_r / r… (cf. Axisymétrie).

Une transformation dégénérée (\(\det F = 0\)) donne des valeurs non finies pour Hencky et Almansi : rien n’est vérifié au point de Gauss.

"jacobian" rend le seul scalaire J, entrée orthoradiale comprise sur un sous-espace axisymétrique. Une maille retournée a J ≤ 0 en un de ses points : c’est ce que le pas à pas vérifie à chaque itéré en grands déplacements. Aucune autre mesure ne le voit : C = FᵀF reste défini positif quand F change de signe, si bien que Green–Lagrange et Hencky se calculent sans broncher sur une maille retournée.

# Linearized strain, the default: eps_xx, eps_xy, eps_yy.
eps = pyrucast.element_field.deformation(u, fes)
# Finite measures share the call; each has its own prefix.
E = pyrucast.element_field.deformation(u, fes, "green_lagrange")  # E_xx…
F = pyrucast.element_field.deformation(u, fes, "F")  # F_xx, F_xy, F_yx, F_yy
H = pyrucast.element_field.deformation(u, fes, "hencky")  # H_xx…

thermal_strain(temperature, materials, fespace, t_ref, measure="linear") → ElementField

Déformation thermique de libre dilatation (Cast3M EPTH), pour la thermomécanique non couplée :

\[ \varepsilon_{th} = \alpha\,(T - T_{ref})\,\big[\,1,1,(1),0,0,0\,\big]. \]

temperature est un champ par éléments portant "T" (p. ex. produit par interp_to_gauss) ; alpha est lu dans le champ matériau, où il voyage comme composante facultative de l’élasticité (à côté de E/nu, cf. Élasticité). La sortie a exactement la même disposition que deformation (composantes normales à α·ΔT, cisaillements nuls — y compris l’orthoradiale eps_zz en axisymétrique, un solide de révolution se dilatant aussi circonférentiellement), si bien que deformation(u, fespace) - thermal_strain(...) donne la déformation mécanique ε(u) − ε_th. Aucun couplage n’est fait ici : l’utilisateur compose la charge thermique et la contrainte réelle σ = D:(ε − ε_th) à partir des briques (integrate_behavior, internal_forces) — cf. l’exemple thermomécanique.

measure est la mesure que lit la loi : la même dilatation y est donnée exactement, ½((1 + αΔT)² − 1) en Green–Lagrange (E_*), ln(1 + αΔT) en Hencky (H_*) — cf. Grandes transformations.

elastic_stress(model, strain, materials) / elastic_strain(model, stress, materials) → ElementField

L’opérateur élastique d’un modèle aux points de Gauss, dans les deux sens (Cast3M ELAS) :

\[ \sigma = \mathbb D : \varepsilon \qquad\text{et}\qquad \varepsilon = \mathbb D^{-1} : \sigma . \]

D est l’opérateur avec lequel chaque loi prédit — pas la réponse de la loi elle-même : la matrice de Hooke d’une élasticité dans sa symétrie (isotrope, orthotrope, anisotrope), le D isotrope du prédicteur élastique d’une plasticité, le D non endommagé que dégrade une loi d’endommagement. Le matériau n’a besoin que des constantes élastiques (E, nu en isotrope) : le sigma_y d’une plasticité n’y joue aucun rôle. C’est ce dont un solveur incrémental a besoin pour transformer un incrément de déformation imposé (thermique, donné) en la contrainte qu’il causerait si rien ne plastifiait — la pseudo-force F* = ∫ Bᵀ D Δε* — ou pour ramener une contrainte à la déformation élastique qui la porte sous d’autres modules ((D⁻¹ₙ₊₁ − D⁻¹ₙ) σₙ quand le matériau varie avec la température).

L’entrée et la sortie en déformation sont en convention tensorielle, disposition de deformation et thermal_strain ; en contraintes planes, elastic_strain n’écrit pas eps_zz, que la loi plane ne lit pas. Une zone par sous-modèle continu ; un sous-modèle sans opérateur élastique (conduction, barre, poutre, contrainte) n’en produit pas. elastic_strain échoue là où D est singulier (un module nul).

Deux fonctions plutôt qu’une qui lirait le sens sur les noms des composantes : les deux sens ne font pas la même opération, et l’état que rend integrate_behavior porte à la fois des sigma_* et des eps_*.

# The stress an imposed strain would cause if nothing yielded: D:ε.
sigma = pyrucast.element_field.elastic_stress(plasticite, dilatation, materiaux)
# And back: the elastic strain carrying a stress, D⁻¹:σ.
retour = pyrucast.element_field.elastic_strain(plasticite, sigma, materiaux)

von_mises(model, stress, coefficients=[0.5, 0.5, 0.5, 1.5, 1.5, 1.5]) → ElementField

La contrainte équivalente aux points de Gauss : von Mises par défaut, la forme anisotrope de Hill avec ses six coefficients [F, G, H, L, M, N] :

\[ \sigma_{eq}^2 = F(\sigma_{yy}-\sigma_{zz})^2 + G(\sigma_{zz}-\sigma_{xx})^2 + H(\sigma_{xx}-\sigma_{yy})^2 + 2L\,\sigma_{yz}^2 + 2M\,\sigma_{xz}^2 + 2N\,\sigma_{xy}^2 \]

Avec \( F = G = H = 1/2 \) et \( L = M = N = 3/2 \), c’est von Mises, \( \sqrt{3 J_2} \) — la valeur par défaut. stress porte les contraintes qu’écrivent les lois (sigma_xx, … : la sortie d’elastic_stress ou l’état d’integrate_behavior conviennent, quoi qu’ils portent d’autre). Le modèle dit lesquelles lire, car elles dépendent de la cinématique :

Le résultat a une zone par sous-modèle continu, sur ses points de Gauss, et une seule composante, von_mises. Une erreur nomme ce qui manque : des coefficients qui laisseraient la forme devenir négative (chacun de L, M, N, F + G, G + H, H + F et FG + GH + HF doit être ≥ 0), ou une contrainte sans sigma_zz en déformations planes — les lois plastiques le portent, une contrainte élastique non. Elle n’existe pas en méthode : elle ne s’applique qu’à un champ de contraintes.

# von Mises at each Gauss point, from the stresses the law writes; the model
# says which ones to read (here a plane stress: σzz = 0).
vm = pyrucast.element_field.von_mises(plasticite, contraintes)
# Hill's anisotropic form: the coefficients [F, G, H, L, M, N].
hill = pyrucast.element_field.von_mises(
    plasticite, contraintes, [0.6, 0.4, 0.5, 1.5, 1.5, 2.0]
)

stress_measure(stress, gradient, target) → ElementField

Une contrainte convertie, point par point, d’une mesure à l’autre, avec le gradient de la transformation F (deformation(u, fespace, "gradient")) : Cauchy sigma_*, Kirchhoff tau_*, Piola–Kirchhoff 1 P_* (non symétrique) et 2 S_*, duale de Hencky T_*. La mesure d’entrée se reconnaît à ses noms, une par zone ; target nomme celle qu’on veut ("cauchy", "kirchhoff", "pk1", "pk2", "hencky"). C’est ainsi que la contrainte lagrangienne d’une loi en grandes transformations devient la P que lisent les forces internes, ou la Cauchy des résultats — la moitié « contraintes » de Cast3M CONF. Une barre se reconnaît à son effort de Piola–Kirchhoff 2 axial S_n : "pk1" en fait le vecteur effort nominal P_x, P_y(, P_z) = S_n F t₀, "cauchy" l’effort vrai n = λ S_n (t₀ sa direction sur la configuration active, la référence, comme F). Formules et exemples : Grandes transformations.

rotate_frame(field, gradient) → ElementField

Un champ de matériau dont les axes (V1X, V1Y en 2D, V1X…V2Z en 3D) tournent avec la pièce, point par point, de la rotation R de la décomposition polaire F = R U : v = R v₀, d’où k = R k₀ Rᵀ. F est lu dans gradient, que produit deformation(u, fespace, "gradient") sur la référence. Les autres composantes sont recopiées telles quelles : la matière est la même, seuls ses axes suivent la rotation. Une zone sans axes (matériau isotrope) revient partagée, sans copie ; une zone qui n’en porte qu’une partie est une erreur. C’est ce que lit une conduction ou une diffusion calculée sur la pièce déformée. Pas de méthode : elle exige des axes et un gradient nommés. Exemple : Grandes transformations.

divergence(field) → NodeField

Divergence faible (consistante) d’un champ vectoriel par éléments — l’adjoint de gradient :

\[ d_i = \int_\Omega \nabla N_i \cdot F\, d\Omega \approx \sum_{\text{cell}} \sum_g (\nabla N_i \cdot F)\big|_g\, |J|_g\, w_g, \]

accumulé par nœud. C’est l’opérateur Bᵀ, transposé du gradient : il vérifie ⟨∇f, F⟩ = ⟨f, div F⟩. Un tenseur se lit à ses noms : symétrique (A_xx, A_xy, A_yy), ou plein s’il porte aussi A_yx — comme la contrainte nominale P des grandes transformations. Le champ d’entrée doit porter exactement space_dim composantes (F_x, F_y, F_z) ; la sortie porte une composante "div", et un nœud partagé par plusieurs sous-espaces reçoit la somme de leurs contributions. Ce sont des quantités intégrées, pas les valeurs ponctuelles de ∇·F (pas de projection L²).

cell_integral(field, target) → ElementField et sum_to_nodes(field) → NodeField

Deux briques par maille (Cast3M INTG 'ELEM' et CHAN 'CHPO' 'SOMM') :

• element_field.cell_integral intègre chaque composante maille par maille, avec la quadrature du champ d’entrée : v_e = Σ_g f(g)·|J|_g·w_g (le 2πr inclus en axisymétrique). Le résultat vit sur target, un espace à un point par maille du même maillage — FiniteElementSpace(mesh, quadrature="REDUCED"), un seul point au barycentre, le support 'GRAVITE' de Cast3M ;
• node_field.sum_to_nodes rassemble un tel champ aux nœuds par somme : chaque nœud reçoit Σ_{e ∋ i} v_e, la somme des valeurs des mailles qui le touchent, zones confondues. Ce n’est pas un lissage : ni moyenne ni poids.

Intégrer 1 donne la mesure de chaque maille ; la sommer aux nœuds donne le volume que représente chaque nœud — le poids d’un solveur incrémental pour mettre ses résidus à une échelle comparable (l’accélération d’Anderson pondère par (Σ_e |Ω_e|)^{-1/2}).

# One point per cell, at its centroid: where per-cell values live.
centres = pyrucast.FiniteElementSpace(mesh, quadrature="REDUCED")
aires = pyrucast.element_field.cell_integral(un, centres)  # ∫ 1 = the cell's area
# Each node receives the sum over the cells around it.
volumes = pyrucast.node_field.sum_to_nodes(aires)

beam_deformation(field, fespace, material) → ElementField

Déformations généralisées d’une poutre — Bernoulli comme Timoshenko — à partir d’un champ de déplacements et de rotations. Un opérateur pour les trois configurations, lues sur la dimension du maillage :

C’est le produit B · u, avec le B de l’élément (models::beam::b_into) — celui-là même dont la rigidité intègre Bᵀ D B et dont les forces internes intègrent le transposé. Les déformations sont évaluées à chaque point de Gauss : la courbure varie le long d’une travée non chargée (M' = V), seul le cisaillement est constant (V' = 0).

Le matériau est exigé, et c’est la signature honnête : Φ = 12EI/(G·A_s·L²) décide de la distribution de courbure. Un matériau sans constantes de cisaillement — celui d’une poutre de Bernoulli, qui ne demande ni G ni A_s — signifie Φ = 0, si bien que le même opérateur sert les deux théories sans qu’on ait à lui dire laquelle. Le résultat se donne au comportement pour obtenir les efforts de section.

shell_deformation(field, fespace, formulation) → ElementField

Déformations généralisées d’une coque, depuis les six DDL u_x…u_z, r_x…r_z. formulation vaut "thick" ou "kirchhoff" : les lignes de membrane, de vrillage et de cisaillement sont partagées, mais celles de flexion sont toute la différence entre les deux.

Toutes dans le repère local de la facette. drill est le résidu de vrillage θ_z − ω_z, dont l’effort conjugué est un moment que la rigidité intègre comme les autres.

Le cisaillement transverse est constant par élément, échantillonné au point réduit — le point même où la rigidité l’intègre. Ce n’est pas une simplification : c’est l’intégration réduite qui empêche une coque mince de bloquer, lue de l’autre côté. L’échantillonner ailleurs rapporterait une déformation que l’élément ne porte pas.

Le B du modèle : strain(model, solution, materials) → ElementField

Les producteurs ci-dessus connaissent un espace EF, pas une loi. strain est leur pendant au niveau du modèle : il demande à chaque loi quel producteur la nourrit, l’applique à ses inconnues sur son espace, et rend l’entrée de toutes les lois du modèle — le EPSI de Cast3M, étendu à toutes les physiques. Aussi en méthode : model.strain(solution, materials).

Chaque loi lit ses inconnues de solution par leur nom : une même solution — déplacement, température, multiplicateurs — sert tout le modèle, et une inconnue absente vaut zéro. Les sous-modèles sans loi (contraintes, charges, échanges) ne produisent rien. materials n’est lu que par les poutres, dont la répartition flexion / cisaillement dépend de Φ = 12EI/(G·A_s·L²). Ce que la loi déclare, c’est son producteur (Behavior::strain_operator, méthode requise : une physique nouvelle ne peut pas l’oublier) ; l’opérateur tient la vue de la solution une fois pour toutes les lois.

strain et internal_forces sont transposés. Les éléments finis portent deux couples de variables duales : déplacement / force aux nœuds, déformation / contrainte aux points de Gauss. Le travail se mesure des deux côtés — f·u aux nœuds, Σ_g w_g |J_g| σ:ε aux points — et Bᵀ W est l’adjoint de B pour ces deux produits :

\[ \sigma \cdot W (B u) = (B^\mathsf{T} W \sigma) \cdot u \qquad \forall\, u,\ \sigma , \]

le principe des travaux virtuels écrit au niveau discret :

   nœuds                         points de Gauss
   u      ───── strain (B) ─────▶  ε
                                   │  integrate_behavior (la loi)
   f_int  ◀── internal_forces ───  σ
                (Bᵀ W)

D’où K = Bᵀ W D B, symétrique si et seulement si D l’est, et, pour une loi linéaire, internal_forces(integrate_behavior(strain(u))) = K·u à l’arrondi près — ce que vérifie l’exemple. divergence est le cas scalaire du même Bᵀ, et l’adjoint de gradient. En grandes transformations, la transposition lie la variation δE = sym(Fᵀ∇₀δu) aux efforts : internal_forces lit alors la contrainte nominale P = F S (cf. Grandes transformations), et la tangente gagne la rigidité géométrique. Elle ne tient plus pour un terme non variationnel (advection, pression suiveuse, relation non associée), ni dès que les deux côtés n’utilisent pas les mêmes points : c’est pourquoi shell_deformation lit le cisaillement au point réduit où la rigidité l’intègre.

Un modèle couplé se passe de consolidate. Deux lois sur un même maillage (thermique et mécanique) donnent deux zones sur un même espace. Chaque consommateur — integrate_behavior (entrée et état précédent), internal_forces, tangent, geometric — prend la zone qui porte ce qu’il lit, comme le matériau est déjà résolu par ses composantes.

Pourquoi strain et pas deformation. En anglais, deformation désigne la transformation elle-même (F, et les tenseurs C, b, qui valent I en mouvement rigide) ; strain en est une mesure, qui s’annule en mouvement rigide (ε, E, ln U), et le mot consacré pour les déformations généralisées des poutres et coques. L’opérateur géométrique, qui produit les deux, garde donc deformation ; ce que lit une loi est un strain. Le gradient thermique en est un abus assumé, comme les internal_forces d’un problème thermique.

# B: each law gets its own input from one solution — a strain for the solid,
# a temperature gradient for conduction.
eps = couple.strain(solution, materiaux)
repos = couple.initial_state(materiaux)  # l'état de départ des deux lois
etat = couple.integrate_behavior(eps, materiaux, repos)  # the laws
f_int = couple.internal_forces(etat, solution, materiaux)  # Bᵀ: back to the nodes

Maths élément par élément

Onze fonctions appliquent une fonction scalaire à chaque valeur d’un champ et renvoient un nouveau champ du même type (style numpy). Elles acceptent indifféremment les quatre saveurs de champ — NodeField, SubNodeField, ElementField, SubElementField — et dispatchent par type.

Les résultats sont non bornés, comme en numpy : aucune protection sur le domaine (log de ≤ 0 donne -inf/nan, sqrt d’un négatif donne nan). Ces fonctions se combinent à l’arithmétique scalaire des champs (f + s, f * s, cf. Champ) pour bâtir des expressions par composante.

# Atténuation exponentielle d'un champ de température.
attenue = pyrucast.field.exp(temperature * -0.1)

# Magnitude of a field (combined with scalar field arithmetic).
amplitude = pyrucast.field.abs(signal)

Réduction

Deux produits scalaires, à ne pas confondre — ils diffèrent par ce qui est réduit et donc par le type du résultat :

Les deux exigent la même saveur d’opérandes (NodeField, SubNodeField, ElementField, SubElementField), alignent les composantes par nom (l’ordre peut différer) et associent les valeurs comme l’arithmétique de Champ : la somme ne porte que sur les (support, composante) partagés par les deux champs ; un support ou une composante d’un seul côté n’a pas de vis-à-vis et ne contribue pas. Un nœud que plusieurs zones portent compte une fois dans xty.

• xty (produit scalaire global, dot / dot_field) → un float ;
• psca (pscal / pscal_field) → un champ (composante "psca"), une zone par support partagé.

xty(x, y) → float

Produit scalaire global des deux champs entiers :

\[ x \cdot y = \sum_i \sum_c x_{i,c}\, y_{i,c}, \]

la somme parcourant toutes les valeurs, un nœud que plusieurs zones portent compté une fois. Le résultat est un unique float : le produit scalaire qui sert au calcul d’énergie (F·u), aux normes de résidu, etc. L’addition flottante n’étant pas associative, le total dépend du nombre de threads jusqu’au dernier ULP — comme le solveur, ce n’est pas reproductible bit à bit.

# External strain energy: work of the nodal forces in the displacement field
# (same components, same mesh).
energie = pyrucast.measure.xty(forces, deplacements)

psca(x, y) → champ (même saveur que les entrées)

Produit scalaire nœud par nœud (ou point par point) — réduction sur les composantes seules, le support est conservé :

\[ p_i = \sum_c x_{i,c}\, y_{i,c}. \]

Le résultat est un nouveau champ de la même saveur que les entrées, portant une seule composante "psca" : la valeur du produit scalaire à chaque nœud. Chaque sortie est écrite une fois (par nœud) ⇒ indépendant du nombre de threads.

# Squared norm of a vector field, node by node.
norme2 = pyrucast.field.psca(vitesse, vitesse)  # one-component field, "psca"

integral(field, component, fespace=None) → float

Intègre un champ sur son support par la quadrature éléments finis, ∫_Ω f dΩ — le total d’une composante (p.ex. la résultante d’une densité de force distribuée) :

\[ \int_\Omega f \, d\Omega \;=\; \sum_{\text{cell}} \sum_g f(\text{cell}, g)\, |J|_g\, w_g . \]

• sur un NodeField : les valeurs nodales sont relevées aux points de Gauss par les fonctions de forme, ∫ Σ_i f_i N_i dΩ — fespace est requis ;
• sur un ElementField : les valeurs (déjà aux points de Gauss) sont intégrées directement — fespace est ignoré.

Comme xty, la somme flottante dépend du nombre de threads jusqu’au dernier ULP. En interne, la réduction parallèle sur les cellules passe par le driver kernel::reduce_cells.

# Resultant of a surface force density f_y on a plate (through N_i).
r_y = pyrucast.measure.integral(densite, "f_y", fespace=fes)
# Measure of the domain: ∫ 1 dΩ.
aire = pyrucast.measure.integral(champ_unite, "u", fespace=fes)

Somme et xtx

Pour une résultante de forces déjà nodales (sortie de internal_forces, réactions…), la résultante est une simple somme par nœud — exposée comme méthode, à côté de min / max :

Sur un champ aux nœuds, sum et xtx comptent une fois un nœud que plusieurs zones portent : il n’y a qu’une valeur (voir Champ). Elles regroupent la somme en parallèle : dépendantes du nombre de threads au dernier ULP (contrairement à min / max, exactes quel que soit l’ordre).

Appelées sans argument, min et max lisent le champ comme la liste plate de ses valeurs — toutes composantes confondues, et, au niveau agrégat, toutes zones confondues. C’est l’esprit de xtx, et la même mise en garde : sur un champ dont les composantes ne portent pas la même unité (sigma_xx à côté de sigma_xy), la réponse est « la plus petite valeur là-dedans », pas une grandeur physique — il faut alors nommer la composante. Pour la liste des extremums composante par composante, une compréhension suffit : {c: f.min(c) for c in f.components()}.

Par défaut xtx somme toutes les composantes. En passant components, on restreint la somme à celles-là (les autres sont ignorées) — utile pour mesurer la norme d’un résidu sur un sous-jeu de degrés de liberté. Une composante absente d’une zone y est simplement ignorée ; l’appel n’échoue que si aucune zone ne porte l’une des composantes demandées.

# Resultant of a nodal force field, component by component.
rx = forces.sum("f_x")
ry = forces.sum("f_y")
# Squared norm of the residual, for a convergence test.
r2 = pyrucast.measure.xtx(residu)
# The same norm, restricted to the translation components alone.
r2_uy = pyrucast.measure.xtx(residu, components=["f_y"])
# Extrema of a named component…
fy_max = forces.max("f_y")
# …or, without an argument, of the whole field, components pooled.
partout = forces.min()

À venir

Le module est conçu pour accueillir d’autres dérivations sur le même patron (champ, espace EF) → champ : projection L² vers les nœuds (project_to_nodes), mesures non linéaires de déformation (Green-Lagrange). Elles arriveront avec les premiers besoins.