Introduction

pyrucast est une librairie d’éléments finis dont le cœur est écrit en Rust et qui expose une API Python. Elle s’inspire des principes de cast3m : un noyau d’objets typés, accompagné de fonctions opérant sur ces objets.

📖 Référence API Rust (rustdoc). Cette documentation couvre les principes, l’architecture et les exemples ; la référence par item (signatures, types, modules) est générée par rustdoc et publiée à côté de ce book : https://pyrucast.github.io/pyrucast/rust/.

Philosophie

• Simplicité avant tout. Le code doit rester maintenable et éditable par un humain non expert ; on évite la sophistication gratuite.
• Dépendances minimales. Tout ajout de dépendance externe (Rust ou Python) requiert un accord explicite.
• Vérification continue. Chaque objet est livré avec des tests unitaires Rust, des doctests, des tests Python et un chapitre de cette documentation.

Modèle d’objets (arbre de dépendances)

Chaque structure ne dépend que des structures qui la précèdent. La plupart viennent par paire zone / agrégat (SubMesh/Mesh…) — c’est le motif Agrégat, avec sa composition par union |.

Coords ── Node                          (NodeId u32 stable ; Node = accesseur RAII)
   │
   ├── NodeField  (agrège des SubNodeField)        valeurs par nœud × composante
   │
   └── Mesh  (agrège des SubMesh, + ElementType)   géométrie
          └── FiniteElementSpace  (agrège des SubFiniteElementSpace)
                 │                  (+ Interpolation, QuadratureRule)
                 ├── ElementField  (agrège des SubElementField)   valeurs aux Gauss
                 └── Model  (agrège des SubModel : physiques + contraintes)
                        └── ops::matrix ──► Matrix   (matrice creuse, DOFs nommés)

Model + ElementField (matériau) ──► stiffness ──► Matrix
Matrix + NodeField (second membre) ──► solve ──► NodeField (solution)

Deux traits transverses factorisent le comportement commun : Aggregate (accès len/[i]/union |) et Field/SubField (composantes nommées, min/max, arithmétique), partagés entre NodeField et ElementField.

Résumé des rôles :

Premiers pas

Exemple minimal en Rust :

#[test]
fn un_maillage_minimal() -> 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(), &[1.0, 0.0])?;

    let mut sm = SubMesh::new(coords.clone(), ElementType::SEG2);
    sm.add_cell(&[a.id(), b.id()])?;

    // The aggregate does not carry the `Coords`: the submeshes are what hold it.
    // `Mesh::from_submesh(sm)` is the shortcut for the single-submesh case.
    let mut mesh = Mesh::empty();
    mesh.add_sub(Handle::new(sm))?;
    println!("{}", mesh); // Mesh: 1 submesh(es), 1 cell(s) total
    Ok(())
}

Exemple équivalent en Python :

import pyrucast

c = pyrucast.Coords(dim=2)
a = c.add_node([0.0, 0.0])
b = c.add_node([1.0, 0.0])

mesh = pyrucast.Mesh(c, "SEG2")
mesh.unit().add_cell([a, b])
print(mesh)  # Mesh: 1 submesh(es), 1 cell(s) total

Correspondance Rust ↔ Python

Cette page liste, par module, les structures (exposées en classes Python) et les fonctions libres (exposées en fonctions de module Python). Elle matérialise la règle de Conventions :

• une structure containers::…::Foo est exposée sous le même nom pyrucast.Foo (le wrapper PyO3 interne PyFoo est masqué) ;
• une fonction libre ops::<module>::f est exposée dans le sous-module de même nom : pyrucast.<module>.f. Le module Rust porte le nom du conteneur produit (mesh, node_field, element_field, matrix, model, coords), ou de l’activité quand il ne produit aucun conteneur (measure, geom, export) ; solver est l’exception nommée. Le miroir est sans exception : aucune fonction libre ne vit au top-level ;
• une surcharge d’opérateur Rust devient un dunder Python (Add → __add__, Index → __getitem__, …) ;
• un constructeur nommé Rust devient un classmethod / constructeur Python.

Les tableaux ci-dessous donnent la forme canonique — la fonction libre. Beaucoup de ces opérations sont aussi exposées comme méthode de leur sujet, selon la règle des trois conditions (Conventions) : premier argument sujet, retour conteneur, sens pour toute instance du type. La règle étant mécanique, la liste n’est pas recopiée ici — elle est vérifiée par un test (tests/python/test_method_exposure.py), qui lit le stub et échoue si une opération éligible perd sa méthode. Ce test porte aussi la liste des exclusions, chacune avec sa raison.

Deux points à connaître, illustrés plus bas : le nom peut changer entre les deux formes (matrix.stiffness(model, mats) / model.stiffness_matrix(mats), element_field.sub_material_field(sub, …) / sub_model.material_field(…)), et les cinématiques (deformation, beam_deformation, shell_deformation, thermal_strain) n’ont pas de méthode : elles exigent des composantes nommées, elles n’auraient pas de sens sur un champ quelconque. Même raison pour internal_forces, qui lit la contrainte de Voigt par nom.

filter_components et rename_component n’ont que la forme méthode (f.filter_components(["u_x"]), f.rename_component("U", "DX")) : un seul conteneur, de petits arguments, une vue dérivée — R1 en fait du vocabulaire du champ, pas un opérateur.

La complétude du miroir est elle aussi vérifiée par un test (tests/python/test_mirror_completeness.py), dans les deux sens : aucun opérateur Rust sans binding Python, et aucune fonction Python sans opérateur Rust. Les dérogations y vivent avec leur raison.

Source de vérité : le #[pymodule] de src/lib.rs (enregistrement des classes et fonctions) et le stub python/pyrucast/_pyrucast/__init__.pyi (signatures typées). Cette page en est un instantané, à régénérer à la main si l’API bouge.

Structures ↔ classes

Le nom de la classe Python est identique au nom de la structure Rust.

Quelques types Rust ne sont pas exposés en classes Python : ce sont des détails d’implémentation (SubModelKind, l’énum des physiques sous SubModel ; DofOrdering, l’ordonnancement des DOFs d’une SubMatrix ; SubValue, OutOfRange, ValueKind et Interpolated, internes à l’Evolution — une valeur tabulée se passe directement en scalaire ou en champ, la politique hors-plage en chaîne "error"/"clamp"/"extrapolate").

Les sous-objets Sub* (SubMesh, SubFiniteElementSpace, SubElementField, SubMatrix, SubModel) ne se construisent pas directement côté Python : ce sont des vues obtenues par indexation de leur parent (parent[i]). On construit toujours au niveau parent — Mesh(coords, type), FiniteElementSpace(mesh), ElementField(fes, comps), Matrix.block(...) — ou par l’opérateur qui rend ce parent (model.heat_conduction(fes)), et on compose plusieurs zones avec | (union — Rust : union). Voir la règle « Agrégats : un ou plusieurs » de CONVENTIONS.md.

Exception : SubEvolution. Seul sous-objet à la fois vue (via evolution[i]) et constructible directement — SubEvolution([(t0, v0), (t1, v1), …]) — car une courbe tabulée n’a pas de « parent » géométrique qui la définirait. On compose ensuite les courbes par zone avec | (cf. Évolution).

Fonctions ↔ fonctions

Toutes les fonctions ops prennent leurs conteneurs par référence et renvoient Result<T> (converti en exception Python RuntimeError). Les signatures ci-dessous omettent le & et le Result pour la lisibilité.

ops::mesh — construction et transformation de maillages

ops::node_field — opérateurs produisant un champ aux nœuds

flux et internal_forces (BSIG, ∫ Bᵀ σ) sont des assemblages, mais leur résultat est un vecteur nodal et non un opérateur : on se range par la sortie.

ops::coords — écriture dans le magasin de coordonnées

ops::element_field — opérateurs produisant un champ aux points de Gauss

ops::measure — réductions à un nombre

ops::field — opérateurs génériques

Leur produit est un conteneur — toujours — mais pas un conteneur déterminé : il dépend de l’argument. La règle « un module par conteneur produit » ne désigne donc pas un module, et ils se rangent par domaine.

La composition de zones passe par l’union (| Python / union Rust) ; l’arithmétique champ + scalaire et champ + champ par les opérateurs +, -, *, /, ** (cf. Opérateurs). Ni l’une ni l’autre n’est une fonction ops — merge(a, b) est juste un alias nommé de a | b.

ops::model — déclaration des physiques

Chaque opérateur rend un Model couvrant tout le support reçu (une zone par sous-espace) ; on compose les physiques hétérogènes avec |. Aucun n’a de forme méthode : le premier argument est le support que le modèle recouvre, pas un sujet qu’on transforme.

Les deux plis du catalogue. Rust nomme la symétrie et la loi par une enum (MaterialSymmetry, ElasticLaw, PlasticLaw, DamageLaw) ; Python n’expose pas ces enums, et replie donc la symétrie en mot-clé symmetry= et déplie les lois en une fonction chacune. Le catalogue est le même des deux côtés — les dérogations correspondantes sont enregistrées, avec leur raison, dans tests/python/test_mirror_completeness.py.

ops::matrix — assemblage des matrices

ops::solver — résolution

ops::export — export vers des formats externes

spill — le swap utilisateur

Les clés du dictionnaire sont threshold, count, largest, live et peak, en octets ; threshold vaut None quand rien ne déborde. Voir Calculs plus gros que la RAM.

archive — sauvegarde et relecture d’un graphe

Les deux seuls verbes qui restent au niveau racine : ils ne produisent aucun conteneur déterminé, mais un dictionnaire de ce qu’on leur a donné.

À l’écriture les types sont connus du compilateur, d’où la tranche de paires ; à la relecture non, d’où la table nommée dont on tire chaque objet avec son type attendu (objets.mesh("clef")?). Voir Sauvegarde et relecture.

ops::geom héberge locate_points (mapping iso-paramétrique inverse, sous le baignage) et project_points (projection au point le plus proche sur une surface, sous le contact) ; ces deux primitives sont internes (API Rust), pas encore exposées en Python — la seule dérogation de module entier du garde-fou de complétude (tests/python/test_mirror_completeness.py). Le nœud le plus proche, lui, n’est pas un opérateur : c’est la méthode mesh.nearest_node([x, y]), des deux côtés.

Opérateurs (dunders ↔ traits Rust)

Toutes les classes implémentent __repr__ (← Debug, vue structurelle) et __str__ (← Display, vue résumée façon cast3m) — voir Conventions. Les autres opérateurs, classe par classe :

Arithmétique sur les champs

L’arithmétique (f + s, f + g, …) existe au niveau zone et au niveau agrégat ; les opérations par composante et entre champs sont portées par les traits SubField / Field. Les dunders +, -, *, /, ** dispatchent selon l’opérande droite : un float déclenche l’arithmétique scalaire, un champ du même type l’arithmétique champ + champ (valeur à valeur).

f + s / f + g renvoie un nouveau champ ; += n’est pas surchargé. La composition de zones n’est pas sur + : c’est l’union | (union en Rust, cf. ci-dessous). L’opérateur + est entièrement réservé à l’arithmétique de champ — scalaire (f + 1.0) et champ + champ valeur à valeur (f + g via merge_components/merge_field) ; p. ex. deux champs constants valant 1 s’additionnent en un champ constant valant 2. Pour fusionner des zones avec vérification (et non additionner) : merge(a, b) ≡ a | b.

Algèbre de Matrix / SubMatrix

Matrix.__mul__ dispatche selon l’opérande droite, comme l’arithmétique de champ ci-dessus : un NodeField déclenche le produit matrice-vecteur (mul_field), un float la mise à l’échelle paresseuse du facteur (voir Matrice creuse). / n’existe que pour le facteur (Matrix n’a pas de division matrice-vecteur), et refuse un diviseur nul (ZeroDivisionError) ou non fini (ValueError). Comme pour les champs, tous ces opérateurs renvoient une nouvelle Matrix — jamais de mutation en place.

+ et - additionnent deux opérateurs ; ils ne sont pas la composition (|, ci-dessous), qui écarte un bloc dont elle tient déjà l’emplacement : k | k vaut k, k + k vaut 2k.

Indexation par clé

Protocole séquence — len(x), x[i], for _ in x

Union | (composition d’agrégats)

La composition d’agrégats est l’union : côté Python elle s’écrit | (comme set | set), côté Rust ce sont les méthodes nommées union / union_sub / union_sub_first / union_subs (renvoient Result<…>). Les sous-objets sont partagés (refcount), jamais copiés ; les contraintes de domaine (même Coords pour Mesh, etc.) restent vérifiées.

Sémantique d’union (uniforme pour tous les agrégats) :

1. Déduplication par handle : un sous-objet dont le Handle désigne un objet déjà présent (cf. Handle::same_object) n’est pas ajouté deux fois.
2. Finalisation (Aggregate::finalize) : par défaut un no-op ; les champs la surchargent pour fusionner les zones partageant un même support (voir plus bas).

Vaut pour les sept agrégats (Mesh, FiniteElementSpace, Model, Matrix, NodeField, ElementField, Evolution) plus Node.

Finalisation des champs (fusion par support)

Après l’union par handle, NodeField et ElementField fusionnent les sous-champs définis sur le même support (même Handle<SubMesh> pour NodeField, même Handle<SubFiniteElementSpace> pour ElementField) :

• le sous-champ fusionné porte l’union des composantes ;
• une composante définie par plusieurs sous-champs doit y avoir la même valeur partout (comparaison exacte), sinon | lève une erreur ;
• pour NodeField, 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.

Ces opérations sont aussi exposées en Rust : ops::node_field::consolidate (NodeField) et ops::element_field::consolidate (ElementField).

+ est réservé à l’arithmétique de champ

L’opérateur + (et -, *, /) reste l’arithmétique scalaire des sous-champs (subfield + 2.0 → ajoute la valeur à chaque composante) et a vocation à porter, à terme, l’addition réelle de champs (valeur au nœud = somme des deux). Il n’est jamais utilisé pour composer des agrégats — c’est | qui s’en charge, sans collision.

Installation et démarrage rapide

Cette page suffit pour utiliser pyrucast en Python : quelques commandes, un premier script, et de quoi vérifier que tout fonctionne. Pour développer sur la librairie (tests Rust, doctests, génération de la documentation, features Cargo), voir Compilation et tests.

Rust pur ? pyrucast s’utilise aussi comme bibliothèque Rust sans Python : pyo3 y est une dépendance optionnelle, donc un build par défaut ne tire ni pyo3 ni libpython. Voir Usage en Rust pur.

Prérequis

• Rust stable, installé via rustup.
• Python ≥ 3.11 (pour l’API Python ; inutile en Rust pur).
• Linux uniquement : les en-têtes Python — python3-dev (Debian/Ubuntu) ou python3-devel (Fedora/RHEL). pyo3 en a besoin pour l’édition de liens ; sur Windows l’installateur officiel les inclut déjà.

Compilation et installation

À la racine du dépôt cloné. pyo3 localise l’interpréteur Python via la variable VIRTUAL_ENV : activez toujours le venv avant cargo ou maturin, sinon la compilation échoue avec error: failed to run the Python interpreter at ....

Linux / macOS (bash)

python3 -m venv .venv
source .venv/bin/activate
pip install --upgrade pip maturin
maturin develop --release --features extension-module,viz-interactive

Windows (PowerShell)

python -m venv .venv
.\.venv\Scripts\Activate.ps1
pip install --upgrade pip maturin
maturin develop --release --features extension-module,viz-interactive

maturin develop compile le module Rust et l’installe dans le venv (mode editable). L’option --release compile en optimisé : recommandé pour tout usage réel (un build debug est typiquement 10× plus lent à l’exécution). Après toute modification du Rust, relancer simplement maturin develop --release — pas besoin de réinstaller ni de réactiver le venv tant qu’il reste actif. Les options –features extension-module,viz-interactive permettent d’activer la visualisation interactive.

Vérification immédiate

python -c "import pyrucast; c = pyrucast.Coords(2); n = c.add_node([0.0, 0.0]); print(c); print(n)"

Sortie attendue :

Coords: dim=2, configs=1 (active="default"), nodes=1 (0 collected), permutation: identity
<Node #0>

Premier script

Un exemple minimal — une Coords 2D, deux nœuds, un maillage POI1 :

import os

import pyrucast

c = pyrucast.Coords(dim=2)
a = c.add_node([0.0, 0.0])
b = c.add_node([1.0, 0.0])

mesh = pyrucast.Mesh(c, "SEG2")  # un sous-maillage
mesh.unit().add_cell([a, b])

# Without a screen (continuous integration, a remote session), `plot()` would fail:
# we fall back to a file. That is the condition winit tests itself.
ecran = os.environ.get("DISPLAY") or os.environ.get("WAYLAND_DISPLAY")
mesh.plot(save=None if ecran else "apercu.svg")

print(c)
print(mesh)  # Mesh: 1 submesh(es), 1 cell(s) total
mesh.dump()

À partir d’ici, le chapitre Introduction présente le modèle d’objets, et la section Objets détaille chaque brique (coordonnées, maillage, champs, modèle physique…). Les chapitres Conduction thermique et Mécanique déroulent des problèmes complets de bout en bout.

Aller plus loin

• Référence API Rust : cargo doc --no-deps --lib --open.
• Exemples Python complets et exécutables : dossier examples/ du dépôt (thermique, treillis, élasticité, poutres…).
• Développer sur pyrucast (tests, doc, features) : Compilation et tests.
• Tout construire d’un coup (prérequis + compilation + tests + les trois documentations + module avec visu interactive) : bash script/build.sh (Linux/macOS) ou .\script\build.ps1 (Windows) — voir Scripts « tout-en-un ».

Aspect informatique

Avant de détailler les objets un par un, ce chapitre donne la vue informatique de pyrucast : comment le code est organisé, comment les objets vivent en mémoire, et les quelques motifs transverses qu’on retrouve partout. Chaque section renvoie au chapitre qui en donne le détail.

Deux langages, un seul cœur

pyrucast est une librairie Rust exposée à Python.

• Le cœur de calcul (structures de données, maillage, champs, assemblage, solveur) est écrit en Rust : typage fort, pas de ramasse-miettes runtime, performances natives.
• Le binding Python est une fine couche pyo3, compilée par maturin en un module d’extension .so / .pyd. Chaque classe Python (pyrucast.Coords, pyrucast.Mesh, …) enveloppe un type Rust ; chaque fonction de module (pyrucast.matrix.stiffness, pyrucast.solver.solve, …) appelle une fonction Rust.

Le binding est un miroir 1:1 du Rust, sans logique propre : une méthode Rust reste une méthode, une fonction libre reste une fonction. La table de correspondance exhaustive est dans Correspondance Rust ↔ Python ; la règle qui décide méthode vs fonction libre est dans Conventions.

   Script Python  ──►  module pyrucast (.so, pyo3)  ──►  crate Rust pyrucast
        |                      |                                |
   pyrucast.Coords        wrapper PyCoords                 struct Coords
   pyrucast.matrix.stiffness     fonction de module              ops::matrix::stiffness

Les objets se désignent par des handles

Les objets ne se relient pas entre eux par des références Rust directes, mais par un Handle<T> : une référence comptée munie de son propre verrou (Arc<RwLock<T>>). Pas de registre à interroger, pas de session à passer — le handle est l’adresse de l’objet.

Trois propriétés en découlent, et elles structurent tout le reste de la librairie :

1. Libération automatique. Clone partage, Drop relâche. Quand le dernier handle d’un objet disparaît, l’objet est détruit — aucune fonction remove() à appeler.
2. Toujours valide. Détenir un handle maintient l’objet en vie : il n’y a pas de référence périmée, et read / write ne peuvent pas échouer.
3. L’identité, c’est le pointeur. same_object dit si deux références désignent le même objet — la base de l’union des agrégats.

L’accès passe par un guard (read / write) qui verrouille ce seul objet le temps de l’opération (RAII). Le détail — guards possédés, granularité du verrouillage, compteur par nœud — est dans Modèle mémoire.

Refcount à deux niveaux

La gestion de durée de vie opère à deux échelles indépendantes :

• au niveau objet : le Handle<T> décide si un objet entier (un Coords, un SubMesh…) est vivant ;
• au niveau interne : à l’intérieur d’un Coords, un second compteur par nœud décide si tel nœud est vivant. Le ramasse-miettes manuel Coords::gc() opère sur ce niveau-là.

C’est pourquoi un nœud reste protégé tant qu’un maillage ou un champ le référence, même si tous les Node utilisateurs ont disparu. Voir Coordonnées et Nœud.

Le motif agrégat / sous-objet

La plupart des conteneurs viennent par paires : un objet zone (Sub…) et son agrégat (une liste de zones partageant la même grammaire d’accès) :

Tous les agrégats exposent la même interface (len, [i], itération, unit()) et la même composition par union | (côté Rust : union). Les sous-objets ne se construisent pas directement : on construit au niveau parent, et on indexe (parent[i]) pour obtenir une vue sur une zone. Ce motif est factorisé dans le trait Aggregate — voir Agrégat.

Union |, pas +. Composer deux zones, c’est l’union (mesh_a | mesh_b), avec partage des sous-objets (refcount) et déduplication par handle. L’opérateur + est réservé à l’arithmétique des champs (cf. Champ).

Trois niveaux d’affichage

Chaque objet implémente trois vues, du plus court au plus complet :

• __str__ (Rust Display) — résumé une ligne, façon listing cast3m ;
• __repr__ (Rust Debug) — vue structurelle bornée, pour le développement ;
• dump() — contenu intégral (valeurs, topologie) imprimé sur la sortie standard, au-delà de ce que repr montre.

Détail dans Conventions.

Erreurs

Toute l’API publique renvoie Result<T, PyrucastError>. Côté Python, PyrucastError est converti automatiquement en RuntimeError. Il n’y a qu’un seul type d’erreur dans la librairie — voir Conventions.

Persistance portable

Un trait unique, Portable (serde + bincode), fixe le contrat d’octets : un format binaire identique Linux ↔ Windows. Au-dessus, pyrucast.save / pyrucast.load écrivent un graphe d’objets et le relisent en préservant le partage — deux champs sur un support restent deux champs sur un support. Voir Sauvegarde et relecture, Conventions et Modèle mémoire.

Pour le développeur

L’organisation des fichiers Rust (où vit chaque morceau) est décrite dans Arborescence.

Mailler une géométrie

Cette page est un guide de choix, pas une référence : elle compare les mailleurs entre eux et dit lequel prendre. La description de chaque opérateur, ses arguments et ses pièges vivent dans Opérateurs → Maillage.

La partie 3D viendra plus tard ; seul le 2D est traité ici.

Mailler en 2D

Quatre opérateurs remplissent l’intérieur d’un contour fermé. Ils prennent tous la même chose — un ou plusieurs contours orientés en SEG2 et une taille de maille visée — et rendent tous un maillage dont le bord est le contour, nœuds compris.

Comment lire les figures

Chaque figure montre le même contour maillé quatre fois. Il est construit une seule fois, puis dupliqué par translate : aucun mailleur ne bénéficie d’une discrétisation de bord différente des autres. La disposition ne change jamais — en haut triangulate_surface et pave_surface, en bas grid_surface et grid_surface2. Le contour est en bleu et ses nœuds en rouge.

Les figures et les chiffres de cette page sortent d’un seul script, examples/comparer_mailleurs_2d.py.

La qualité citée est la mean ratio du pire coin : 1 pour un coin droit à côtés égaux, 0 pour un coin plat, négatif pour une maille retournée. C’est la même mesure pour les triangles et pour les quadrangles, ce qui est la seule façon de comparer les quatre sur un pied d’égalité. On donne aussi le 5ᵉ centile, qui dit ce que valent les mailles médiocres et non la seule plus mauvaise.

En pratique cette mesure ne descend jamais à zéro dans les tableaux qui suivent, et ce n’est pas une chance : les quatre mailleurs refusent de rendre un maillage portant une maille retournée ou plate. Une telle maille a un jacobien négatif ou nul, qu’aucun code éléments finis n’intègre — ce n’est pas un maillage médiocre mais un maillage faux, et mieux vaut une erreur qui situe la maille fautive. Un mailleur qui échoue là où un autre passe est donc un résultat en soi, à lire dans les tableaux comme tel.

Une forme rectilinéaire posée sur la grille

C’est le cas le plus favorable aux mailleurs en grille, et la différence est franche.

Les deux mailleurs en grille rendent le maillage qu’on dessinerait à la main : toutes les mailles sont des rectangles et il n’y a pas un seul triangle. Le paveur, lui, rend une pelure d’oignon avec quatre coutures diagonales — sa faiblesse est là où deux de ses rangées se rencontrent, et cette ligne-là existe même sur un rectangle.

Une forme rectilinéaire qui ne tombe pas sur la grille

Dès que les cotes ne sont plus des multiples de la taille visée, grid_surface et grid_surface2 se séparent nettement.

grid_surface pose une ligne sur la coordonnée où repose chaque côté aligné, puis découpe entre deux lignes d’après le côté qui les enjambe : toutes ses lignes sont droites. grid_surface2 donne à chaque nœud du contour la ligne qui le traverse et laisse ses rangées plier pour aller chercher le contour — une même rangée peut alors rejoindre deux parois qui se font face à deux hauteurs différentes, ce qu’aucune droite ne sait faire.

Le L étiré à 1,02 est le seul de la série que grid_surface2 ne règle pas : les deux côtés y posent un nombre de nœuds différent sur la même portée, dix contre onze, et aucune disposition de lignes n’invente la rangée manquante. Le remède est dans le contour, pas dans le mailleur — voir la règle de discrétisation.

Un profil à sept angles rentrants

Le profil crénelé est décliné en deux versions dont seule la base change : coupée sous chaque barre, ou d’un seul tenant. C’est le meilleur révélateur de la sensibilité d’un mailleur à la discrétisation du contour.

grid_surface2 rend exactement le même nombre de mailles dans les deux cas. Les trois autres paient la base d’un seul tenant. C’est la propriété la plus utile de ce mailleur : il pardonne une discrétisation de contour que les autres font payer.

Une forme oblique ou courbe

Ici le classement s’inverse, et c’est le seul endroit où il le fait.

Une grille ne peut pas suivre une oblique : elle la découpe en escalier, et tout ce qui s’en approche est rendu au paveur frontal. Sur la maison, où l’oblique ne fait que le toit, ce qui reste de rectilinéaire suffit encore à grid_surface2 ; sur le cercle, qui n’a plus rien de droit, triangulate_surface est le meilleur des quatre.

Sur une forme sans direction dominante, grid_surface2 n’est pas un candidat. Il prend l’écartement de ses lignes sur le contour, et un cercle n’en dicte aucun : rien ne pose de ligne entre y = 0,24 et y = 1, le vide y vaut quinze fois l’écartement moyen, et le cœur s’effondre — 0,005, une maille au bord de la dégénérescence. grid_surface, lui, rend 0,288 avec un 5ᵉ centile de 0,796, de loin le meilleur des quatre. Sur une courbe, prenez grid_surface, ou triangulate_surface si la qualité du pire élément commande.

Le cercle mérite un mot de plus. pave_surface y tombe à 0,031 — ses rangées se rejoignent au centre en une étoile à quatre branches, et c’est cette ligne-là qui porte tout le défaut. triangulate_surface, lui, ne descend jamais sous 0,42, sur aucune des onze formes : c’est la signature d’un Delaunay raffiné, qui garantit un angle minimal et rien de plus. Il ne s’effondre jamais, mais il ne monte jamais non plus.

Ce qu’il faut retenir

Le coût en mailles est le classement le plus stable : sur les dix formes rectilinéaires ou peu obliques, grid_surface2 ≤ grid_surface < pave_surface < triangulate_surface, avec un facteur quatre entre les extrêmes à qualité au moins égale. Pour un calcul, c’est le facteur qui compte juste après la qualité. Le cercle est la seule exception, et il l’est dans les deux colonnes à la fois : grid_surface2 y rend plus de mailles pour une qualité effondrée, ce qui est la façon la plus nette de dire qu’il n’est pas fait pour ça.

En pratique :

• forme rectilinéaire — grid_surface2, et d’autant plus si ses côtés n’ont pas été coupés aux angles qui leur font face ;
• forme courbe — triangulate_surface si la qualité du pire élément commande, grid_surface s’il faut des quadrangles ; jamais grid_surface2, qui y perd son cœur ;
• forme franchement oblique — pave_surface, dont les rangées épousent la pente ; mais tant qu’il reste des parois droites autour de l’oblique, les mailleurs en grille tiennent — la maison en est l’exemple ;
• quadrangles obligatoires — pave_surface accepte all_quad=True et refuse par une erreur claire un contour dont la parité l’interdit ;
• triangles voulus — triangulate_surface, seul à en produire par construction.

Le front peut garder ses angles : relax

Les chiffres de pave_surface ci-dessus sont ceux du réglage par défaut. Entre deux rangées, le front est relaxé — c’est ce qui l’empêche de se plisser — et cette relaxation est un laplacien, qui arrondit les angles. Or un front ne perd des nœuds qu’à ses coins : une fois les coins arrondis, il garde tous ses nœuds pendant que son périmètre rétrécit, et le milieu du domaine sort plus fin que la taille demandée. Un carré 20 × 20 à la taille 1 sort à 600 mailles au lieu de 400, les plus intérieures à 0,43 de l’aire visée.

relax="along" garde le même déplacement mais projeté sur le front : l’écartement s’égalise encore, la forme n’est plus rabotée, et le carré sort comme les 400 carrés exacts qu’on dessinerait. relax="none" ne relaxe pas du tout. Aucun des trois ne gagne partout — sur une courbe, il n’y a pas d’angle à préserver et le front a tout à gagner à se redresser :

Le détail du mécanisme est sur la page Mailler des opérateurs. Les deux mailleurs en grille prennent le même réglage, qui n’y gouverne que la bande : la grille est posée droite quoi qu’il arrive.

Deux limites connues

Un vide que le contour ne borde pas n’est pas maillé en grille. grid_surface2 prend l’écartement de ses lignes sur le contour ; là où le contour ne dit rien — le triangle du toit de la maison, qu’aucune paroi verticale ne borde — il coupe le vide en mailles entières tant qu’il ne dépasse pas trois fois l’écartement moyen, et l’abandonne au front au-delà. Remplir un tel vide de rangées ne crée pas du cœur : ces rangées n’existent que pour être érodées par l’oblique qui les traverse, et elles morcellent le travail du front au lieu de lui donner une région propre. Mesuré sur la maison : l’abandon coûte un triangle là où le remplissage en coûtait cinq, et le 5ᵉ centile monte de 0,644 à 0,676.

C’est aussi ce qui coûte le cercle, dont les vides valent quinze fois l’écartement moyen. Aucun seuil ne sépare les deux cas — abandonner un vide de huit moyennes abandonne a fortiori un vide de quinze — et c’est un choix assumé : grid_surface2 sert les formes rectilinéaires, grid_surface sert les courbes.

Le cercle est un objectif instable. Sur une forme sans direction dominante, déplacer la grille d’un millième change la pire maille du simple au décuple. Les chiffres du cercle ci-dessus sont justes mais ne se prolongent pas : ne réglez aucun paramètre dessus.

Reproduire les figures

maturin develop --features extension-module,viz
python examples/comparer_mailleurs_2d.py

# Figures du livre :
PYRUCAST_IMG_DIR=book/src/img python examples/comparer_mailleurs_2d.py

Le script tient dans une fonction : elle prend un contour, le duplique trois fois par translate, maille chaque copie par une méthode et réunit le tout sur une figure.

def comparer(nom, contour, taille, fichier, noeuds=True):
    """Maille `contour` par les quatre méthodes et trace la figure.
    / Meshes `contour` with all four methods and draws the figure.

    FR — Le contour n'est construit qu'une fois : ses trois copies viennent de
    `translate`, qui rend un maillage neuf dans les mêmes `Coords` — c'est ce
    qui permet de les réunir sur une seule figure.
    """
    points = [n.position() for sub in contour for cell in sub for n in cell]
    xs, ys = [p[0] for p in points], [p[1] for p in points]
    dx, dy = (max(xs) - min(xs)) * 1.15, (max(ys) - min(ys)) * 1.25
    coins = [(0.0, 0.0), (dx, 0.0), (0.0, -dy), (dx, -dy)]

    tout, lignes = None, []
    for (titre, mailler), coin in zip(MAILLEURS, coins):
        copie = pc.mesh.translate(contour, list(coin))
        maillage = mailler(copie, taille)
        qs = sorted(
            qualite([n.position() for n in cell]) for sub in maillage for cell in sub
        )
        comptes = dict(zip(maillage.element_types(), maillage.cell_counts()))
        lignes.append(
            f"{titre:20s} {len(qs):6d} mailles, {comptes.get('TRI3', 0):5d} tri, "
            f"pire {qs[0]:.3f}, p5 {qs[len(qs) // 20]:.3f}"
        )
        for sub in copie:
            sub.face_color = (20, 90, 200)
        bloc = maillage | copie
        if noeuds:
            points_rouges = pc.mesh.to_poi1(copie)
            for sub in points_rouges:
                sub.face_color = (220, 30, 30)
            bloc = bloc | points_rouges
        tout = bloc if tout is None else tout | bloc

    print(f"\n{nom}   (taille visée {taille:g})")
    for ligne in lignes:
        print(f"   {ligne}")
    tout.plot(
        view=VUE,
        show_axes=False,
        save=os.path.join(OUT, fichier),
        title=f"{nom} — h={taille:g} — haut : triangulate, pave — bas : grid, grid2",
    )

Mailler en 3D

À venir.

Objets

Cette partie décrit, un par un, les objets du modèle pyrucast : pour chacun, son principe (à quoi il sert, comment il est structuré) et son interface (API Rust et Python). Les objets sont présentés dans leur ordre de dépendance — chacun ne s’appuie que sur les précédents.

Coords ── Node
   │
   ├── Mesh (agrège des SubMesh)
   │      └── FiniteElementSpace (agrège des SubFiniteElementSpace)
   │             ├── ElementField (champ aux points de Gauss)
   │             └── Model (agrège des SubModel : physiques)
   │                    └── Matrix (matrice creuse, sortie d'assemblage)
   └── NodeField (champ aux nœuds)

Evolution (agrège des SubEvolution : valeur tabulée vs variable, interpolée)

Une distinction structure toute cette liste et se lit dans l’arborescence : Coords est le magasin (src/coords.rs), Node et Cell sont des atomes insécables (src/atoms/), et tout le reste — les sept agrégats et leurs vues Sub* — sont des conteneurs divisibles (src/containers/), les seuls qui puissent être le sujet d’un opérateur. Voir Conventions.

Deux abstractions transverses traversent la liste et méritent d’être lues tôt :

• l’Agrégat — la grammaire commune (len, [i], union |) de tous les conteneurs « parent / zones » (Mesh, FiniteElementSpace, Model, Matrix, NodeField, ElementField, Evolution) ;
• le Champ (Field / SubField) — le contrat partagé entre NodeField et ElementField : composantes nommées, min/max, arithmétique scalaire et par composante.

Pour le contexte informatique (handles, refcount, motif zone/agrégat), voir Aspect informatique. Pour les écrire dans un fichier et les relire tels qu’ils étaient, Sauvegarde et relecture.

Coordonnées (Coords)

Un Coords héberge une ou plusieurs configurations (jeux de coordonnées) pour le même ensemble de nœuds, en dimension fixée. C’est le premier objet du modèle pyrucast — tous les autres (Mesh, NodeField, FE space…) viennent s’y greffer. L’accesseur utilisateur d’un nœud, le Node, fait l’objet du chapitre suivant.

Repère de révolution

Un Coords déclare aussi comment lire ses coordonnées : cartésien par défaut, ou axisymétrique — le plan méridien \( (r, z) \) d’un solide de révolution, avec \( x = r \ge 0 \) (rayon) et \( y = z \) (axe). C’est l’équivalent du OPTI MODE AXIS de Cast3M, et sa place est bien la géométrie : le repère change la mesure d’intégration elle-même,

\[ d\Omega = 2\pi r \, |J| \, d\xi, \]

donc rigidité, masse, conductivité, flux réparti, volumes et forces internes d’un coup — sur le corps comme sur ses bords. Le facteur \( 2\pi \) est celui de l’anneau complet : les masses, volumes et résultantes nodales sont ceux de la pièce de révolution entière.

use pyrucast::coords::Coords;

#[test]
fn le_repere_se_choisit_a_la_construction() {
    // Cartesian (the default): free dim.
    let plan = Coords::new(2).unwrap();
    assert!(!plan.is_axisymmetric());

    // Revolution: the dimension is necessarily 2, hence no argument.
    let axi = Coords::axisymmetric().unwrap();
    assert_eq!(axi.dim(), 2);
    assert!(axi.is_axisymmetric());
}

import pyrucast

c = pyrucast.Coords.axisymmetric()
assert c.dim == 2 and c.is_axisymmetric
c.add_node([1.0, 0.0])  # r = 1, z = 0
try:
    c.add_node([-1.0, 0.0])  # x is a radius: it must be ≥ 0
except RuntimeError as erreur:
    print(erreur)

Un rayon négatif est refusé à l’ajout (et au set_position) plutôt que de ressortir en |J| négatif au fond d’une intégrale. Tout espace éléments finis bâti sur ces Coords hérite du repère, si bien qu’un corps et son bord ne peuvent pas diverger.

Côté mécanique, l’axisymétrie ajoute la déformation orthoradiale \( \varepsilon_{\theta\theta} = u_r/r \), qui relève du modèle et non de la géométrie : voir Élasticité linéaire. La thermique n’a rien à changer.

Côté visualisation, un tracé axisymétrique montre par défaut la section méridienne ; l’option revolve la balaie pour dessiner le corps de révolution lui-même — voir Visualisation.

Identité d’un nœud

Chaque nœud créé reçoit un identifiant interne stable (NodeId), unique pour toute la vie du Coords : aucun id n’est jamais réutilisé, même après ramassage par le GC. C’est ce qui permet aux maillages et champs de référencer un nœud par son id sans s’inquiéter de la stabilité.

Politique de suppression : pas de suppression directe

Il n’existe aucune méthode remove_node. Un nœud référencé est protégé. Seul le ramasse-miettes Coords::gc() retire les nœuds dont le refcount interne est tombé à 0.

use pyrucast::handle::Handle;

#[test]
fn un_noeud_survit_tant_qu_on_le_tient() {
    let coords = Handle::new(Coords::new(2).unwrap());
    // add_node initializes refcount = 1; without a decrement, the node is protected.
    let id = coords.write().add_node(&[0.0, 0.0]).unwrap();
    assert_eq!(coords.write().gc(), 0);

    // Après décrément, gc ramasse.
    coords.write().decref(id).unwrap();
    assert_eq!(coords.write().gc(), 1);
}

Modèle de refcount à deux niveaux

        Handle<Coords>             ◀── refcount de l'Arc
                │                       (le Coords est-il vivant ?)
                ▼
        ┌──────────────────┐
        │      Coords      │
        └──────────────────┘
                │
                │ refcount par nœud
                ▼                   ◀── refcount sur chaque NodeId interne
        NodeId(0)  NodeId(1)  …        (le nœud est-il vivant ?)

Les deux niveaux sont indépendants :

• tant qu’un Handle<Coords> existe, le Coords reste en mémoire ;
• tant qu’au moins un Node (ou un objet aval comme Mesh / Field via incref/decref) référence un NodeId, ce nœud est protégé du GC.

Le détail du niveau objet est dans Modèle mémoire ; le niveau nœud, dans Nœud.

Plusieurs configurations

Utile pour basculer entre référence / déformée / prédite. La configuration active est désignée par index ; lire les coordonnées d’un nœud (node.position() en Python, Coords::position côté Rust) renvoie celles de la configuration active. add_config(name) clone la configuration active sous un nouveau nom.

Rust :

#[test]
fn une_seconde_configuration_clone_la_courante() {
    let coords = Handle::new(Coords::new(2).unwrap());
    let id = coords.write().add_node(&[0.0, 0.0]).unwrap();

    let c2 = coords.write().add_config("deformed");
    coords.write().select(c2).unwrap();
    // the following `set_position` now change the "deformed" configuration.
    coords.write().set_position(id, &[0.1, 0.05]).unwrap();

    coords.write().select(0).unwrap();
    assert_eq!(coords.read().position(id).unwrap(), vec![0.0, 0.0]);
    coords.write().select(c2).unwrap();
    assert_eq!(coords.read().position(id).unwrap(), vec![0.1, 0.05]);
}

Python :

import pyrucast

c = pyrucast.Coords(dim=2)
n = c.add_node([0.0, 0.0])

# Create a second configuration (a clone of the active one).
c2 = c.add_config("deformed")
print(c.names())  # ['default', 'deformed']

# Switch to the deformed configuration and change the coordinates.
c.select(c2)
n.set_position([0.1, 0.05])

# The coordinates read depend on the active configuration.
c.select(0)
print(n.position())  # [0.0, 0.0]  — configuration de référence
c.select(c2)
print(n.position())  # [0.1, 0.05] — configuration déformée
print(c.active)  # 1

Pourquoi plusieurs configurations dans un Coords plutôt que plusieurs Coords ?

cast3m suit historiquement la convention inverse : un objet de coordonnées = une seule configuration, et on multiplie les objets. Les deux modèles ont des compromis distincts.

Rien n’énumère les Coords vivants — chacun n’est atteignable que par les handles qui le désignent (cf. Modèle mémoire) : propager un add_node « à tous les Coords équivalents » n’aurait de toute façon aucun point d’entrée.

Choix pyrucast : un seul jeu d’identités par domaine géométrique, plusieurs configurations pour les variantes (référence / déformée / prédite). Des Coords distincts restent prévus pour les vrais domaines indépendants (sous-domaines, maillages adaptés de plus haute densité).

Permutation solveur

Une permutation optionnelle (Vec<u32>, longueur = capacity) sépare l’ordre solveur de l’identité : permutation[node_id] donne l’ordre solveur associé. Elle est posée par l’appelant aujourd’hui ; une renumérotation réduisant la bande/profil (Cuthill–McKee) la calculera. L’identité (NodeId) n’est jamais modifiée, dans les deux cas.

Rust :

#[test]
fn une_permutation_renumerote_pour_le_solveur() {
    let coords = Handle::new(Coords::new(2).unwrap());
    // Trois nœuds créés ; ids = 0, 1, 2.
    coords.write().add_node(&[0.0, 0.0]).unwrap();
    coords.write().add_node(&[1.0, 0.0]).unwrap();
    coords.write().add_node(&[0.5, 1.0]).unwrap();

    // Permutation set by hand (the automatic computation is still to be written).
    coords.write().set_permutation(vec![2, 0, 1]).unwrap();
    // permutation[0] = 2: the node with id 0 is at solver position 2.
    println!("{:?}", coords.read().permutation());

    // Back to the identity.
    coords.write().clear_permutation();
    assert!(coords.read().permutation().is_none());
}

Python :

import pyrucast

c = pyrucast.Coords(dim=2)
c.add_node([0.0, 0.0])
c.add_node([1.0, 0.0])
c.add_node([0.5, 1.0])

# Set a permutation by hand.
c.set_permutation([2, 0, 1])
print(c.permutation())  # [2, 0, 1]

# Back to the identity (None = identity).
c.clear_permutation()
print(c.permutation())  # None

API Python

import pyrucast

c = pyrucast.Coords(dim=2)
n = c.add_node([0.0, 0.0])  # n is a pyrucast.Node; refcount = 1
m = c.add_node([1.0, 0.0])

print(c)  # Coords: dim=2, configs=1 (active="default"), nodes=2 ...
n.set_position([0.5, 0.5])

# GC touches nothing as long as at least one Python Node exists.
assert c.gc() == 0

# del + collect force le Drop côté Rust et libère le refcount.
import gc as pygc

del n
pygc.collect()
assert c.gc() == 1

Méthodes d’inspection utiles : node_count() (nœuds vivants), capacity() (slots alloués, vivants + non encore collectés), is_alive(id) et refcount(id) — qui prennent un id brut (pas un Node) : un Node portant un refcount ne pourrait jamais être observé mort. acquire(id) rend un Node supplémentaire pour un id existant (refcount += 1). dump() imprime le contenu intégral (coordonnées, configurations) sur la sortie standard.

Nœud (Node)

Le Node est l’accesseur utilisateur d’un nœud d’un Coords. Il ne stocke pas de coordonnées : il porte un Handle<Coords> et un NodeId, et délègue tout au Coords (les coordonnées lues dépendent donc de la configuration active). C’est l’objet qu’on passe partout où une API attend un nœud.

Principe : un identifiant qui se protège

Un Node est conceptuellement une paire (Coords, NodeId), mais avec une propriété de plus : il maintient le refcount interne automatiquement, par RAII.

• Clone incrémente le refcount du nœud dans le Coords ;
• Drop le décrémente.

Tant qu’un Node (ou un objet aval — un SubMesh, un champ — qui a fait son propre incref) référence un NodeId, ce nœud est protégé du ramasse-miettes Coords::gc(). C’est le niveau interne du refcount à deux niveaux décrit dans Coordonnées.

#[test]
fn un_noeud_est_un_compteur_de_references() -> Result<()> {
    let coords = Handle::new(Coords::new(2)?);
    let n = Node::create_in(coords.clone(), &[1.0, 2.0])?;
    let m = n.clone(); // refcount = 2
    drop(n); // refcount = 1
    drop(m); // refcount = 0
    coords.write().gc(); // collecte
    Ok(())
}

Le code interne peut toujours manipuler directement les NodeId sans passer par Node, mais perd alors la protection automatique : il doit appeler manuellement Coords::incref / Coords::decref (c’est ce que font les maillages et les champs).

Création

Un Node ne se construit jamais « dans le vide » : il naît d’un Coords.

• Rust : Node::create_in(coords_handle, &[x, y, …]) crée le nœud et rend le Node (refcount = 1). Pour obtenir un Node supplémentaire sur un id déjà existant, Node::acquire(coords_handle, id) — côté Rust le verbe est porté par l’atome, côté Python par le magasin.
• Python : coords.add_node([x, y, …]) renvoie directement un pyrucast.Node. coords.acquire(id) rend un accesseur de plus.

Interface

Côté Python, le Node expose :

• la propriété id — l’identifiant entier stable dans son Coords ;
• position() — ses coordonnées dans la configuration active ;
• set_position([x, y, …]) — réécrit ses coordonnées (dans la configuration active) ;
• coords() — le Coords auquel il appartient ; filet de secours quand la poignée a été lâchée côté Python, comme Mesh.coords() ;
• l’union node | node → un Mesh POI1 unitaire sur les deux nœuds (la même union | que les agrégats, cf. Agrégat) ; et mesh | node ajoute un point à un Mesh POI1 unitaire ;
• les vues repr / str et dump().

import pyrucast

c = pyrucast.Coords(dim=2)
a = c.add_node([0.0, 0.0])
b = c.add_node([1.0, 0.0])

print(a.id)  # 0
print(a.position())  # [0.0, 0.0]
a.set_position([0.5, 0.5])

# Union de nœuds → maillage POI1 (deux points).
poi = a | b
print(poi)  # Mesh: 1 submesh(es), 2 cell(s) total

Côté Rust, Node expose id(), position(), set_position(&[…]), coords() (le Handle<Coords> porté), plus Clone/Drop qui gèrent le refcount comme décrit ci-dessus.

Pourquoi un accesseur protecteur plutôt qu’un simple id ?

Manipuler des NodeId bruts est possible mais dangereux : rien n’empêche un nœud d’être ramassé pendant qu’on tient encore son id. Le Node rend la protection automatique et locale — exactement comme un Rc/Arc le ferait pour un objet du heap, mais ici au niveau du nœud interne d’un Coords. Une alternative (modèle « cast3m pur » où un nœud n’existe qu’au sein d’un maillage) est discutée dans Modèle mémoire ; pyrucast conserve pour l’instant le Node protecteur.

Agrégat

La plupart des conteneurs de pyrucast viennent par paires : un objet zone (SubMesh, SubFiniteElementSpace, SubNodeField, SubElementField, SubModel, SubMatrix, SubEvolution) et son agrégat (Mesh, FiniteElementSpace, NodeField, ElementField, Model, Matrix, Evolution). Tous les agrégats partagent exactement la même grammaire d’accès et la même composition par union — factorisées dans le trait Rust Aggregate (src/aggregate.rs). Ce chapitre décrit ce contrat commun une fois pour toutes ; chaque chapitre d’objet n’en redonne que les spécificités.

Principe : une liste de handles

Un agrégat est, au fond, un Vec<Handle<Sub>> : une liste de handles vers des sous-objets adressés par handle. Il ne copie jamais ses zones — il en partage les handles (refcount). Le trait Aggregate dérive toute la mécanique d’accès de deux accesseurs seulement (items() / items_mut()), si bien qu’ajouter un nouvel agrégat ne duplique aucun code d’indexation.

   Mesh                FiniteElementSpace        Model
   ├── Handle<SubMesh> ├── Handle<SubFES>        ├── Handle<SubModel>
   ├── Handle<SubMesh> ├── Handle<SubFES>        └── Handle<SubModel>
   └── …               └── …

Interface uniforme

Tout agrégat expose, côté Python :

Côté Rust, le trait Aggregate fournit les mêmes : len, is_empty, get(i), subset(indices), iter, unit, push/add_sub/add_subs, plus Index<usize> et IntoIterator via la macro impl_aggregate_std_traits!.

agg[i] vue, agg[i:j] agrégat. L’indexation entière renvoie une vue sur une zone (un Sub…) ; le slicing renvoie un agrégat du même type que agg, dont les zones sont des handles partagés (pas de copie). Le slicing s’appuie sur subset côté Rust et préserve les invariants de l’agrégat (check_push, finalize).

Passer une seule zone à un opérateur. Les opérateurs (ops/) prennent des agrégats, pas des vues : mesh.invert(me[1]) lève un TypeError (« SubMesh object is not an instance of Mesh »). Pour n’en traiter qu’une, utiliser le slicing, qui rend un agrégat unitaire : mesh.invert(me[1:2]).

unit() vs [0]. agg[0] prend silencieusement la première de plusieurs zones ; agg.unit() exige qu’il n’y en ait qu’une et lève sinon. Aux frontières qui n’ont de sens que mono-zone (par exemple ajouter une cellule à un maillage qu’on vient de créer avec un seul sous-maillage), préférer unit() — plus honnête (cf. CONVENTIONS.md).

Les sous-objets ne se construisent pas seuls

Un Sub… obtenu par agg[i] est une vue : on ne le construit pas directement côté Python. On construit toujours au niveau parent (Mesh(coords, type), FiniteElementSpace(mesh), ElementField(fes, comps)) ou par un opérateur qui rend un parent (model.heat_conduction(fes)…), puis on indexe pour atteindre une zone.

Composition : l’union |

Composer deux agrégats, c’est l’union : a | b côté Python, a.union(&b) côté Rust. La sémantique est uniforme pour les sept agrégats :

1. Déduplication par handle. Une zone dont le Handle désigne un objet déjà présent (Handle::same_object) n’est pas ajoutée deux fois. L’ordre est celui de première apparition.
2. Partage, pas copie. Les zones retenues sont partagées (refcount), jamais dupliquées en mémoire.
3. Contraintes de domaine. Les invariants (même Coords pour un Mesh, etc.) restent vérifiés au moment de l’ajout (check_push).
4. Finalisation. Un crochet finalize() tourne en fin d’union : no-op pour la plupart des agrégats, mais les champs le surchargent pour fusionner les zones partageant un même support (voir Champ, Champ aux nœuds).

Les quatre formes de l’union :

L’union est donc acceptée dans les deux sens : sub | agg passe par le __ror__ de l’agrégat (une zone seule ne sait unir qu’une autre zone), et ne diffère de agg | sub que par l’ordre des zones. La déduplication et la finalisation sont identiques — une zone déjà présente est simplement déplacée en tête.

En place : add_sub / add_subs. L’union rend un agrégat neuf et laisse ses deux opérandes intacts. Pour modifier l’agrégat courant, les deux primitives bas niveau sont agg.add_sub(sub) (une zone) et agg.add_subs(other) (toutes les zones d’un autre agrégat, dans l’ordre). Elles vérifient check_push mais ne dédupliquent pas et n’appellent pas finalize : c’est de la concaténation. L’API d’usage pour composer reste | (cf. CONVENTIONS.md).

Plus, pour les nœuds (cf. Nœud) :

| compose, + calcule. L’union est toujours |. L’opérateur + (et -, *, /) est réservé à l’arithmétique des champs (field + 2.0, addition de champs…) — il n’est jamais utilisé pour composer des agrégats. Les deux ne se télescopent donc jamais. C’est un changement par rapport aux toutes premières versions, où la composition passait par + ; aujourd’hui + est entièrement libéré pour le calcul.

Exemple

import pyrucast

c = pyrucast.Coords(dim=2)
ns = [c.add_node(p) for p in [(0, 0), (1, 0), (1, 1), (0, 1)]]

tri = pyrucast.Mesh(c, "TRI3")
tri.unit().add_cell([ns[0], ns[1], ns[2]])

qua = pyrucast.Mesh(c, "QUA4")
qua.unit().add_cell([ns[0], ns[1], ns[2], ns[3]])

# Union of two meshes (zones shared by handle).
mesh = tri | qua
print(len(mesh))  # 2 sous-maillages
print(mesh)  # Mesh: 2 submesh(es), 2 cell(s) total

# In-place addition: `add_sub` for one zone, `add_subs` for all those of
# another aggregate (concatenation, without deduplication).
tri.add_subs(qua)
print(len(tri))  # 2 sous-maillages

Pourquoi un trait commun ?

Factoriser l’accès et l’union dans Aggregate garantit qu’un sous-maillage, un sous-espace EF ou un sous-modèle s’indexent, s’itèrent et se composent strictement de la même façon. Aucun __getitem__ n’est réécrit par type, aucun comportement d’union ne diverge d’un agrégat à l’autre, et un nouvel agrégat se câble en deux macros (impl_aggregate! + impl_aggregate_pymethods!). Le seul point de personnalisation est finalize() (et check_push()), que les champs exploitent pour leur fusion par support.

Maillage (Mesh / SubMesh)

Le maillage de pyrucast se compose à deux niveaux, selon le motif agrégat / zone :

• SubMesh : regroupe toutes les cellules d’un même ElementType. Stocke la connectivité à plat (un Vec<NodeId> de longueur cell_count × nodes_per_cell).
• Mesh : agrège plusieurs SubMesh liés à la même Coords.

Le cas POI1 est volontairement dégénéré : un sous-maillage POI1 est exactement une liste de nœuds (un nœud par cellule), ce qui sert de support naturel aux NodeField.

Ce chapitre décrit l’objet maillage (structure, types d’éléments, refcount). La construction de maillages par des générateurs (line, triangulate_surface, extrude…) relève des opérateurs de maillage.

Types d’éléments

L’enum ElementType liste les types supportés ; chaque variante porte sa propre méta-information (nombre de nœuds, dimension topologique, nom court cast3m).

Les huit derniers types sont quadratiques (Lagrange-2) : ils reprennent la numérotation des sommets de leur parent linéaire puis ajoutent les nœuds de milieu d’arête, dans l’ordre d’arêtes de la convention VTK (voir le rustdoc d’ElementType). 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 = quadrangle/hexaèdre bi-/tri-quadratiques, avec nœuds de face et central). Ils se posent avec l’interpolation LAGRANGE2 (cf. Espace éléments finis).

Ajouter un nouveau type d’élément est purement additif : un fichier src/atoms/element_kind/<nom>.rs et une variante — voir Ajouter un élément fini.

Cellule (Cell)

mesh.cell(submesh_idx, cell_idx) (ou mesh[submesh_idx][cell_idx]) renvoie une vue Cell sur une cellule : len(cell) donne son nombre de nœuds et cell[k] le k-ième Node. C’est l’accès lecture à la connectivité ; l’ajout passe par submesh.add_cell([...]).

Refcount sur les nœuds — interaction avec le GC

Chaque appel à SubMesh::add_cell incrémente le refcount interne de chaque nœud dans la Coords (cf. Coordonnées). Le Drop du SubMesh les décrémente. Tant qu’un SubMesh référence un nœud, le ramasse-miettes le protège — même si tous les Node utilisateurs ont disparu.

   Coords             ◀── refcount par NodeId
        │                          (le nœud est-il vivant ?)
        ├── Node(s) utilisateur(s) ── chacun +1
        └── SubMesh(s)             ── chacun +1 par cellule incidente

En cas d’échec partiel d’add_cell (par exemple un nœud déjà ramassé), les incréments déjà effectués pour la cellule courante sont annulés (rollback transactionnel à l’échelle d’une cellule).

Construire en bloc. Un add_cell prend le verrou d’écriture de la Coords et lâche les caches dérivés du sous-maillage : c’est le bon grain pour une maille posée à la main, mais sur un million de mailles c’est ce prix-là qu’on paie, pas la connectivité. Un opérateur qui produit un gros maillage passe donc par SubMesh::from_connectivity(coords, type, connectivité) (côté Rust), qui valide et incrémente tout le tableau en une seule prise de verrou — une unité par occurrence, comme toujours — et par Coords::add_nodes pour créer ses nœuds d’un coup. Tous les opérateurs de maillage l’empruntent : les balayages (sweep, extrude, revolve, sweep_solid), transfinite, les copies rigides (translate, rotate, symétries, copy), merge_nodes, convert, to_quadratic, skin, border, chain, orient, select, consolidate, barycenter, to_poi1, la lecture gmsh et les mailleurs frontaux. Un nœud dont la position n’est pas encore arrêtée — celle qu’un front déplace tant qu’il avance, celle qu’un relâchement va lisser — n’est même plus créé avant de l’être : il vit comme un rang dans un tableau, et ne devient un NodeId qu’à la fin.

Scellement (connectivité figée après consommation)

Par convention, un maillage n’est plus modifié une fois construit : un espace éléments finis, un champ ou une matrice indexent ses cellules, et lui ajouter une maille par la suite les laisserait dans un état incohérent.

Cette convention est désormais imposée. Dès qu’un objet autre qu’un maillage capture un SubMesh et indexe ses cellules (construction d’un SubFiniteElementSpace, d’un support de SubMatrix…), ce sous-maillage est scellé : sa connectivité est figée pour toujours. add_cell / add_cell_taking renvoient alors l’erreur MeshSealed. Un Mesh qui se contente de contenir le sous-maillage ne le scelle pas — il peut continuer à grossir tant qu’aucun consommateur ne s’y attache. On teste l’état avec is_sealed.

Un champ aux nœuds fait exception, parce qu’il n’indexe pas les cellules du maillage : il se pose sur le compagnon POI1 de la zone, et c’est ce compagnon-là qui est scellé. Le maillage donné en argument, lui, reste modifiable ; le modifier lâche son compagnon, si bien qu’un champ construit ensuite se retrouve sur un autre support — les champs déjà construits ne sont pas invalidés pour autant, ils restent définis sur le nuage de nœuds d’avant (voir Champ aux nœuds).

Pour repartir d’un maillage scellé et le modifier à nouveau, on en prend une copie profonde avec duplicate() : un SubMesh (ou Mesh) neuf, non scellé, avec la même connectivité (les nœuds sont partagés — même Coords —, seuls leurs refcounts augmentent). L’opérateur mesh.copy(m, new_nodes=…) est la même chose au niveau Mesh, avec le choix supplémentaire de recréer des nœuds neufs aux mêmes endroits plutôt que de partager ceux de l’original.

mesh = pyrucast.Mesh(c, "TRI3")
mesh.unit().add_cell([a, b, n3])

pyrucast.FiniteElementSpace(mesh)  # scelle mesh[0]
assert mesh[0].is_sealed
# mesh[0].add_cell([...])           # → RuntimeError (MeshSealed)

copie = mesh.duplicate()  # neuf, modifiable
copie.unit().add_cell([b, n3, n4])  # OK

API Rust

#[test]
fn un_maillage_se_compose_de_zones() -> 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(), &[1.0, 0.0])?;
    let c = Node::create_in(coords.clone(), &[0.5, 1.0])?;

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

    let sm_handle = Handle::new(sm);
    let mut mesh = Mesh::empty(); // l'agrégat ne porte pas la `Coords`
    mesh.add_sub(sm_handle)?;
    assert_eq!(mesh.cell_count(), 1);
    Ok(())
}

API Python

import pyrucast

c = pyrucast.Coords(dim=2)
a = c.add_node([0.0, 0.0])
b = c.add_node([1.0, 0.0])
n3 = c.add_node([0.5, 1.0])

# Mesh(coords, type) creates a mesh with a single submesh; unit() gives the
# view of it, add_cell adds a cell.
mesh = pyrucast.Mesh(c, "TRI3")
mesh.unit().add_cell([a, b, n3])
print(mesh)  # Mesh: 1 submesh(es), 1 cell(s) total
print(mesh.element_types())  # ['TRI3']
print(mesh.cell_counts())  # [1]

# Composer plusieurs zones : l'union | (jamais +).
quad = pyrucast.Mesh(c, "QUA4")
# … add_cell … ;  combined = mesh | quad

Durée de vie et refcount

Les SubMesh et Mesh portent un effet de bord dans leur Drop : ils décrémentent le refcount des nœuds qu’ils référencent dans la Coords. Cet effet a lieu exactement une fois, quand le dernier Handle sur le maillage disparaît. Le détail est dans le chapitre Modèle mémoire.

Visualisation

mesh.plot(...) trace le maillage (chaque sous-maillage avec sa propre couleur, ou coloré par un champ). Voir Visualisation.

Espace éléments finis (FiniteElementSpace)

Ce chapitre couvre FiniteElementSpace et son objet associé SubFiniteElementSpace : la couche éléments finis posée par-dessus la couche géométrique (Mesh / SubMesh). Le Mesh reste purement géométrique ; le FiniteElementSpace lui ajoute la formulation (fonctions de forme, points de Gauss, dérivées). Cette séparation permet de réutiliser un même maillage avec plusieurs formulations, et — plus tard — d’admettre un déplacement du maillage tout en réévaluant correctement les grandeurs physiques.

Architecture

Hiérarchie miroir de celle du maillage :

Mesh                                  FiniteElementSpace
├── SubMesh (ElementType)             ├── SubFiniteElementSpace (Interpolation, QuadratureRule)
├── SubMesh                           ├── SubFiniteElementSpace
└── ...                               └── ...

• SubFiniteElementSpace détient un Handle<SubMesh>, une Interpolation et une QuadratureRule. Il porte les tables de référence précalculées une fois pour toute (points de Gauss, fonctions de forme et dérivées de référence évaluées à ces points).
• FiniteElementSpace détient un Handle<Mesh> figé et un Vec<Handle<SubFiniteElementSpace>> en correspondance un-pour-un avec les sous-maillages. La topologie (connectivité, types d’éléments) est figée à la construction.

POI1 n’est pas un élément fini au sens classique : un sous-maillage POI1 est rejeté à la construction d’un SubFiniteElementSpace.

Le catalogue détaillé — fonctions de forme, dérivées et quadrature de chaque type — est en section Éléments finis supportés. Ce chapitre-ci décrit la machinerie commune (mapping isoparamétrique, Jacobien, gradients physiques) partagée par tous.

Conventions de l’élément de référence

Chaque ElementType fixe son repère de référence \( \xi \) et la numérotation locale de ses nœuds. Ces conventions sont aussi documentées dans le rustdoc de ElementType et reproduites ici pour référence centrale.

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

QUA8, HEX20, PENTA15 sont sérendipité (arêtes seulement) ; SEG3, TRI6, TET10, QUA9, HEX27 sont des Lagrange complets (QUA9/HEX27 = tenseurs Q2 complets, avec nœuds de face et central).

Ces conventions sont cohérentes avec celles déjà imposées ailleurs dans le code : orientation CCW des triangles produits par triangulate_surface, ordre HEX8/PENTA6 utilisé par extrude et sweep_solid, ordre des nœuds de milieu d’arête aligné sur VTK (export verbatim) et réaligné à la lecture gmsh.

Théorie : élément isoparamétrique

Transformation géométrique de référence

Soit un élément avec \( n \) nœuds, de coordonnées physiques \( \mathbf{x}_1, \dots, \mathbf{x}_n \in \mathbb{R}^{d_s} \) (avec \( d_s \) la dimension géométrique de la Coords). La transformation géométrique \( \chi : \hat{K} \to K \) qui envoie l’élément de référence \( \hat{K} \) sur l’élément physique \( K \) est interpolée par les mêmes 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 est interpolée exactement comme un champ scalaire. Avec une interpolation Lagrange-1, \( \chi \) est affine sur les simplexes (SEG2, TRI3, TET4) et tri-linéaire sur les tenseurs (QUA4, HEX8).

Fonctions de forme Lagrange-1

Une fonction de forme Lagrange est définie par la propriété de Kronecker : \( N_i(\xi_j) = \delta_{ij} \) aux nœuds de référence. Les formules explicites sur les cinq types supportés sont :

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 (tri-liné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) \]

Une propriété immédiate, vérifiée par les tests unitaires : \( \sum_i N_i(\xi) = 1 \) en tout \( \xi \) (partition de l’unité).

Fonctions de forme quadratiques (Lagrange-2)

L’interpolation Lagrange2 couvre les six types quadratiques. Pour les Lagrange complets, avec les coordonnées barycentriques \( L_i \) :

• 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 \).

Les sérendipité n’ont pas de nœud de face/intérieur :

• QUA8 : 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.
• QUA9 : produit tensoriel complet \( N_i(\xi,\eta) = \ell_{a}(\xi)\,\ell_{b}(\eta) \) des trois fonctions de Lagrange 1D \( \ell_{-}(t)=\tfrac12 t(t-1),\; \ell_0(t)=1-t^2,\; \ell_{+}(t)=\tfrac12 t(t+1) \).
• HEX27 : produit tensoriel 3D \( N_i(\xi,\eta,\zeta) = \ell_{a}(\xi)\,\ell_{b}(\eta)\,\ell_{c}(\zeta) \) des mêmes fonctions 1D.
• HEX20 : 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 permutations).
• PENTA15 : produit du TRI6 par un facteur quadratique en \( \zeta \), avec correction sérendipité aux sommets.

Toutes vérifient Kronecker (\( N_i(\xi_j)=\delta_{ij} \)), partition de l’unité, et leurs dérivées analytiques sont recoupées par différences finies dans les tests.

Fonctions de forme cubiques d’Hermite (Hermite-3)

Les deux familles précédentes sont C⁰ : le champ est continu d’un élément au suivant, sa dérivée ne l’est pas. Cela suffit à toute équation du second ordre, et à rien d’autre. 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 uniquement :

\[ \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 à la fois : à 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’annulent dans les deux. 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.

Deux conséquences structurelles

Quatre fonctions pour deux nœuds. Le nombre de fonctions de forme cesse d’être le nombre de nœuds ; c’est shape_count() qui le donne, et tous les tableaux de référence sont dimensionnés dessus. nodes_per_cell() reste ce qu’il était, et continue de dimensionner la géométrie.

L’élément devient sous-paramétrique. La géométrie d’un SEG2 est un segment droit, interpolée en Lagrange-1, quel que soit le champ qu’il porte. L’espace tabule donc deux bases :

Elles coïncident pour toute interpolation de Lagrange, et l’accesseur de champ se rabat alors sur la géométrique — d’où un coût mémoire nul dans le cas courant, et aucun consommateur existant à modifier.

Dérivées secondes

Hermite3 est la seule famille à tabuler \( \partial^2 N_i/\partial\xi^2 \), et c’est cohérent : la courbure est une grandeur primaire pour un élément C¹, alors qu’une base de Lagrange n’en a aucun usage. Sur SEG2 elles sont linéaires en \( \xi \) — donc la courbure varie dans l’élément, là où un SEG2 de Lagrange en donnerait une identiquement nulle.

Le passage au physique est une simple règle de chaîne, le terme \( \partial J/\partial\xi \) disparaissant 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₂ et H₄ portent une pente \( \partial w/\partial\xi \), tandis 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 est délibérément absent des fonctions de forme : c’est exactement le passage référence → physique que toute autre base traverse, donc il vit là où tous les autres jacobiens sont appliqués. Écrite autrement — avec un L dans la fonction de forme — la base cesserait d’être une grandeur de l’élément de référence, et ne pourrait plus être tabulée une fois par type d’élément.

Ce qui le vérifie

La raideur d’Euler-Bernoulli est le seul oracle capable de falsifier cette base :

\[ 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 classique, à la précision machine. Une erreur dans les fonctions, dans leurs dérivées secondes ou dans le facteur jacobien des pentes atterrit dans cette comparaison — ce qui vaut mieux que seize assertions sur la base seule. L’intégrande étant quadratique en \( \xi \), deux points de Gauss l’intègrent exactement.

Pas de base du tout : MODEL_EMBEDDED

Un élément structurel dont la matrice élémentaire est une forme fermée — une barre, un portique — n’évalue jamais une fonction de forme. Son intégrale a été faite une fois, au crayon, et seul le résultat figure dans le code.

Déclarer une interpolation de Lagrange sur un tel espace énonce quelque chose que la physique n’utilise pas, et parfois quelque chose de faux : la forme fermée de l’élément de portique exact vient de fonctions cubiques et quadratiques, pas linéaires.

MODEL_EMBEDDED le dit. C’est la troisième combinaison des deux bases :

La géométrie reste entièrement définie — coordonnées, jacobien, mesure — donc toute la plomberie d’assemblage est inchangée : champ matériau, sortie de comportement, layout, coloriage, dispersion parallèle. Seule l’interpolation de l’inconnue n’appartient pas à l’espace.

Le refus a des dents

Un accesseur de champ sur un tel espace erronne, en nommant la situation, plutôt que de rendre la base géométrique. La distinction compte parce que la base géométrique est disponible et paraîtrait plausible : c’est exactement le repli silencieux que cette variante existe pour empêcher.

Ce qui suppose que les consommateurs posent la question au bon endroit. Les opérateurs qui interpolent une inconnue sont donc passés à l’accesseur de champ :

Sur une poutre, le premier est le plus parlant : la charge répartie cohérente vaut qL/2 et qL²/12 en moment, ce qui suppose de connaître la base. Un repli linéaire donnerait qL/2 sans moment — plausible, et faux. Refuser est la bonne réponse.

La visualisation, elle, reste sur la base géométrique, et délibérément : ce qu’on colorie est ce qu’on dessine — un segment droit, une facette plane — donc la couleur doit varier le long du tracé. C’est une image, pas une valeur calculée.

Et ce que la formulation doit alors fournir elle-même

La reconstruction des efforts passe à la charge de la formulation, puisque l’espace ne peut plus interpoler. beam_deformation évalue donc les déformations depuis les fonctions de forme de l’élément, ce qui l’oblige à recevoir le matériau : Φ en dépend. C’est le prix, et il est juste — on ne peut pas reconstituer la courbure d’une poutre sans connaître sa raideur de cisaillement.

Dérivées de référence

Les dérivées de référence \( \partial N_i / \partial \xi_k \) suivent par dérivation directe. Pour Lagrange-1 sur les simplexes (TRI3, TET4) et sur SEG2, ces dérivées sont constantes sur l’élément. Sur QUA4 et HEX8 — et sur tous les types quadratiques — elles sont polynomiales en les coordonnées de référence.

En dérivant la partition de l’unité, on obtient également \( \sum_i \partial N_i / \partial \xi_k = 0 \) pour chaque direction de référence \( k \). Cette identité est aussi testée.

Stockage flat

Le buffer plat des dérivées de référence retourné par Interpolation::dshape_dxi(et, &xi) est de longueur \( n_\text{nodes} \times d_r \) (avec \( d_r = \dim \hat{K} \)) et row-major : \[ \mathtt{dN}[i \times d_r + k] = \frac{\partial N_i}{\partial \xi_k}(\xi) \]

Théorie : quadrature de Gauss

Pour intégrer une fonction \( f \) sur l’élément physique, on remonte à 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| \) le déterminant (au sens généralisé, défini ci-dessous) du Jacobien. Le couple \( (\xi_g, w_g) \) est la règle de quadrature.

pyrucast utilise une règle « par défaut » par type d’élément, calibrée pour intégrer exactement la matrice de masse Lagrange-1 sur un élément à géométrie droite. Les règles sont :

Les types quadratiques utilisent une règle exacte pour leur matrice de masse (degré 4) sur géométrie droite ; l’exactitude des règles custom TRI6 et TET10 est vérifiée par des tests d’intégration de monômes.

La somme des poids vaut le volume 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 (vérifié par les tests, pour tous les types à la fois).

Théorie : Jacobien et grandeurs physiques

Jacobien

Le Jacobien de la transformation \( \chi \) est la matrice de 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 \in {0, \dots, d_s - 1} \) parcourt les directions physiques et \( k \in {0, \dots, d_r - 1} \) les directions de référence. Le buffer plat retourné par SubFiniteElementSpace::jacobian(cell, g) suit la convention row-major \( \mathtt{J}[a \times d_r + k] = J{a,k} \).

Cas standard : \( d_s = d_r \)

Quand le maillage et son espace ambiant ont la même dimension (par exemple TRI3 dans une Coords 2D), \( J \) est carrée. Le déterminant ordinaire \( \det(J) \) mesure la dilatation locale du volume ; son valeur absolue intervient dans l’intégration. La fonction SubFiniteElementSpace::det_jacobian retourne \( |\det(J)| \).

La dérivation des fonctions de forme par rapport aux coordonnées physiques utilise 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 manifold : \( d_s > d_r \)

Un sous-maillage peut être plongé dans un espace de dimension supérieure : SEG2 dans une Coords 2D ou 3D (contour, courbe), TRI3 dans une Coords 3D (surface plongée). C’est exactement ce que produit triangulate_surface quand on lui donne un contour 3D plan.

Dans ce cas, \( J \) est rectangulaire (taille \( d_s \times d_r \)). Le déterminant standard n’a plus de sens, mais on peut définir la métrique tirée en arrière : \[ G(\xi) = J(\xi)^T\, J(\xi) \quad \text{de taille } d_r \times d_r \] \( G \) est symétrique définie positive (sous condition de non-dégénérescence). L’élément de mesure devient \[ d\mu = \sqrt{\det G}\, d\xi \] qui s’utilise comme \( |J| \) dans la quadrature. La fonction det_jacobian retourne ce \( \sqrt{\det G} \) — toujours positif par construction.

Le gradient tangent d’un champ sur la surface (la projection du vrai gradient sur l’espace tangent) est donné par la pseudo-inverse : \[ \nabla_s N_i = J\, G^{-1}\, \nabla_\xi N_i \] La fonction SubFiniteElementSpace::dn_dx retourne ces composantes, dans le repère ambiant à \( d_s \) dimensions. Pour \( d_s = d_r \), ces formules se réduisent à celles du cas standard (\( J G^{-1} = J^{-T} \)).

Combinaisons valides

Le couple \( (d_r, d_s) \) possible pour notre v0 :

\( d_s < d_r \) n’a pas de sens (impossible de définir un élément 2D dans un espace 1D) ; le constructeur de SubFiniteElementSpace rejette ce cas.

Stratégie de stockage : invariant vs variant à la déformation

Le maillage support d’un FiniteElementSpace est topologiquement figé, mais ses coordonnées peuvent évoluer (déplacement de maillage, mise à jour incrémentale). pyrucast scinde donc les grandeurs en deux catégories selon leur invariance vis-à-vis de cette déformation :

Ce choix donne deux propriétés importantes :

1. Empreinte mémoire indépendante du nombre de cellules. Un SubFiniteElementSpace ne stocke que de l’ordre de \( n_g \times n_\text{nodes} \times d_r \) flottants — quelques centaines au plus par sous-espace. À comparer aux GB qu’un précalcul des Jacobiens demanderait sur un maillage 3D fin.
2. Robustesse au déplacement. Réécrire les coordonnées dans la Coords (par exemple via Coords::set_position) suffit à mettre à jour automatiquement toutes les évaluations de \( J \), \( |J| \) et \( \partial N_i / \partial x_a \) — pas d’invalidation à signaler, pas de cache à reconstruire.

Le coût est CPU plutôt que mémoire : chaque appel à jacobian(cell, g) recalcule la somme \( J = \sum_i \mathbf{x}i\, \nabla\xi N_i \). En pratique, l’assemblage matrice-élémentaire procède cellule par cellule : on calcule \( J \), \( |J| \), \( \nabla_x N_i \) une fois par couple (cellule, Gauss), puis on les réutilise pour tous les termes intégrés. Le surcoût reste donc proportionnel à \( n_\text{cells} \times n_g \) — soit le minimum incompressible — et non à \( n_\text{cells} \times n_g \times n_\text{termes} \).

Si une mesure montrait un jour que ce recalcul devient un goulot d’étranglement, un cache invalidé sur incrément d’un compteur de version du Coords pourrait être ajouté sans changer l’API publique — déclenché par la mesure, pas par anticipation.

Validation à la construction

SubFiniteElementSpace::new rejette à la création :

• un sous-maillage de type POI1 (pas de repère de référence) ;
• un couple (ElementType, Interpolation) non supporté (par exemple TRI3 + Lagrange2 tant que Lagrange2 n’existe pas) ;
• un couple (ElementType, QuadratureRule) non supporté ;
• \( d_s < d_r \) (incompatible avec la définition du Jacobien).

FiniteElementSpace::with (et donc lagrange1, new) vérifie en plus :

• que le maillage contient au moins un sous-maillage ;
• que la longueur de la liste (interpolation, quadrature) correspond au nombre de sous-maillages ;
• que chaque SubFiniteElementSpace se construit sans erreur.

Le déterminant du Jacobien n’est pas vérifié à la construction : un élément dégénéré ou inversé ne sera détecté qu’à la première évaluation de det_jacobian. Cette validation paresseuse permet de construire un FE space sans toucher aux coordonnées et reste correcte vis-à-vis du déplacement de maillage ultérieur (un élément valide peut devenir dégénéré après un déplacement et inversement).

API Rust

Constructeur principal — Lagrange-1 partout, quadrature de Gauss par défaut :

#[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 :
    // affine mapping, det(J) = 4 = 2 × the triangle's physical area (1/2 × 2 × 2).
    for g in 0..s.gauss_count() {
        let dj = s.det_jacobian(0, g)?;
        assert!((dj - 4.0).abs() < 1e-12);
    }
    Ok(())
}

Constructeur explicite — utile pour mélanger les interpolations / quadratures par sous-maillage :

#[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(())
}

Évaluation des grandeurs sur une cellule :

#[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)?;
            // … use these buffers in the element matrix assembly …
            let _ = (n, dn, jac, det_j, dn_dx);
        }
    }
    Ok(())
}

Déplacement de maillage : exemple

Après modification des coordonnées dans la Coords, les évaluations à la volée reflètent automatiquement le nouvel état :

#[test]
fn deplacer_un_noeud_change_les_evaluations_a_venir() -> 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(), &[1.0, 0.0])?;
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    mesh.add_cell(&[a.id(), b.id()])?;
    let sub = FiniteElementSpace::lagrange1(&mesh)?.get(0)?;

    // Initial SEG2: nodes at x=0 and x=1 → |J| = 0.5 (length 1 over [-1,+1]).
    let dj_before = sub.read().det_jacobian(0, 0)?;
    assert!((dj_before - 0.5).abs() < 1e-12);

    // Stretch: the second node moves to x=4 → |J| = 2.0 (length 4 over [-1,+1]).
    coords.write().set_position(b.id(), &[4.0, 0.0])?;
    let dj_after = sub.read().det_jacobian(0, 0)?;
    assert!((dj_after - 2.0).abs() < 1e-12);
    Ok(())
}

API Python

L’objet est exposé sous le nom pyrucast.FiniteElementSpace, avec les sous-espaces accessibles via pyrucast.SubFiniteElementSpace. Les interpolations et règles de quadrature sont passées en chaînes de caractères ("LAGRANGE1", "GAUSS").

import pyrucast

c = pyrucast.Coords(dim=2)
n0 = c.add_node([0.0, 0.0])
n1 = c.add_node([2.0, 0.0])
n2 = c.add_node([0.0, 2.0])

mesh = pyrucast.Mesh(c, "TRI3")
mesh.unit().add_cell([n0, n1, n2])

# Default constructor: Lagrange1 + Gauss everywhere.
fes = pyrucast.FiniteElementSpace(mesh)
assert len(fes) == 1  # 1 sous-espace = 1 sous-maillage
sub = fes[0]  # typed view of subspace 0
assert sub.element_type == "TRI3"
assert sub.interpolation == "LAGRANGE1"
assert sub.quadrature == "GAUSS"
assert sub.gauss_count() == 3
assert sub.space_dim == 2
assert sub.ref_dim == 2

# Evaluations at a given Gauss point.
for g in range(sub.gauss_count()):
    print(sub.gauss_xi(g), sub.gauss_weight(g))
    print(sub.n_at_g(g))  # N_i(ξ_g), flat
    print(sub.dn_at_g(g))  # ∂N_i/∂ξ_j(ξ_g), flat

# Physical quantities (on the fly) on cell 0.
print(sub.jacobian(0, 0))  # J, flat row-major
print(sub.det_jacobian(0, 0))  # |J|, scalaire
print(sub.dn_dx(0, 0))  # ∂N_i/∂x_a, flat row-major

Variantes de construction :

# Same Lagrange1 + same Gauss for every submesh, explicitly.
fes = pyrucast.FiniteElementSpace(mesh, interpolation="LAGRANGE1", quadrature="GAUSS")

# "Class method" form equivalent to the default constructor.
fes = pyrucast.FiniteElementSpace.lagrange1(mesh)

# (Interpolation, quadrature) explicit per submesh.
fes = pyrucast.FiniteElementSpace.with_choices(mesh, [("LAGRANGE1", "GAUSS")])

Déplacement du maillage : le Jacobien reflète automatiquement les coordonnées courantes du Coords — pas de cache à invalider.

print(sub.det_jacobian(0, 0))  # |J| initial

# Moving a node → every evaluation to come sees the
# nouvelles coordonnées.
n1.set_position([4.0, 0.0])
print(sub.det_jacobian(0, 0))  # |J| recalculé

Limitations actuelles

• Quatre interpolations : Lagrange1 (types linéaires), Lagrange2 (types quadratiques SEG3, TRI6, QUA8, QUA9, TET10, PENTA15, HEX20, HEX27), Hermite3 (C¹, sur SEG2 seul) et ModelEmbedded (aucune base de champ : la formulation possède la sienne). Pour les familles de Lagrange le degré doit correspondre au type d’élément : un maillage quadratique se pose avec interpolation="LAGRANGE2" (le constructeur par défaut LAGRANGE1 refuse un type quadratique, et inversement). HERMITE3, lui, se pose sur un SEG2 sans en être le degré — c’est le cas sous-paramétrique. Les ordres supérieurs (Lagrange-3…) restent à venir.
• Deux quadratures : QuadratureRule::Gauss (la règle standard par défaut par ElementType) et QuadratureRule::Reduced (intégration réduite : un point au centroïde, exact pour les constantes — utilisée par exemple pour le terme de cisaillement de la poutre de Timoshenko, anti-verrouillage). Les variantes d’ordre supérieur viendront comme nouvelles variantes. Voir le tableau croisé quadrature × élément pour la compatibilité couple par couple.
• POI1 n’est pas un élément fini (pas de repère de référence). Un sous-maillage POI1 dans le maillage support fait échouer la construction du FiniteElementSpace.
• Pas encore de cache invalidable des grandeurs physiques (\( J \), \( |J| \), \( \nabla_x N_i \)) : tout est recalculé à la volée. Une optimisation à base d’invalidation par compteur de version pourra être ajoutée si la mesure le justifie, sans changement d’API.

Champ (Field / SubField)

Les deux familles de champs de pyrucast — le champ aux nœuds et le champ aux points de Gauss — partagent un contrat commun, capturé par deux traits Rust (src/containers/field.rs) :

• SubField — une zone : un bloc homogène de valeurs, des composantes nommées + un buffer plat ;
• Field — le niveau agrégat, blanket-implémenté pour tout Aggregate dont la zone est un SubField.

Ce chapitre décrit ce contrat (composantes, statistiques, arithmétique) ; les deux chapitres suivants donnent les spécificités de chaque famille.

Principe du trait

SubField : un bloc homogène

Une zone porte :

• une liste ordonnée de composantes nommées ("UX", "UY", "T", "k", "sigma_xx", …) — au moins une, noms uniques ;
• un buffer plat de f64 dans lequel l’indice de composante varie le plus vite (stride = nombre de composantes).

Le contrat est purement structurel : peu importe que les « lignes » du buffer soient des nœuds (SubNodeField) ou des couples (cellule, point de Gauss) (SubElementField) — dès qu’une zone fournit ses composantes et son buffer, le trait en dérive le reste : recherche d’une composante par nom, min/max par composante, opérations scalaires.

   SubField (composantes [c0, c1], stride = 2)
   values = [ ligne0.c0, ligne0.c1, ligne1.c0, ligne1.c1, … ]
                                     └─ composante varie le plus vite ─┘

Chaque zone connaît aussi son support (support() : un Handle<SubMesh> pour un champ aux nœuds, un Handle<SubFiniteElementSpace> pour un champ aux Gauss). Deux zones sont « sur le même support » (same_support) si leurs handles désignent le même slot — c’est la précondition pour les combiner.

Field : le repli sur les zones

Au niveau agrégat, Field replie les opérations sur les zones :

• components() — l’union des composantes des zones (ordre de première apparition) ; une composante peut n’exister que sur certaines zones ;
• min(c) / max(c) — repliés sur les zones qui définissent c (erreur si aucune) ; sans argument, sur tout le champ, toutes composantes et toutes zones confondues ;
• view() — une vue zéro-copie (un guard de lecture par zone), utilisée par les opérateurs qui font beaucoup de lectures (gradient, solveur, viz).

Opération arithmétique

L’arithmétique des champs se décline en trois familles, du plus simple au plus contraint.

1. Scalaire (broadcast)

field + s, field - s, field * s, field / s, field ** s renvoient un nouveau champ où l’opération est appliquée à toutes les valeurs de toutes les composantes de toutes les zones. Disponible au niveau zone (SubField) et au niveau agrégat (NodeField / ElementField, dunders Python __add__, …, __pow__).

scaled = mat * 1.1  # nouveau champ, toutes composantes × 1.1
shifted = u - 5.0  # nouveau champ
energy = u**2.0  # element-wise power (fractional exponent is fine)

+= n’est pas surchargé : f + s ne mute pas f. (Côté Rust, la version consommante est zéro-copie, la version par référence clone d’abord.)

Puissance — Python seulement. Rust n’a pas d’opérateur de puissance ; +/-/*// passent par Add/Sub/Mul/Div, mais ** est exposé par le seul dunder __pow__. Côté Rust, faire une puissance via les primitives génériques : combine_scalar(|a, b| a.powf(b), s) (scalaire) ou merge_components(other, |a, b| a.powf(b)) (binaire). La forme ternaire pow(x, y, z) (modulo) est refusée — elle n’a pas de sens sur des flottants.

2. Par composante (en place)

Pour ne toucher qu’une composante, sur toutes les zones qui la portent : add_to_component(c, s), sub_to_component, mul_to_component, div_to_component — en place, erreur seulement si aucune zone ne définit c (la division par zéro est refusée). Au niveau zone, set_uniform(c, v) force une composante à une valeur constante.

mat.mul_to_component("E", 0.95)  # scales "E" only

3. Binaire entre champs

Combiner deux champs valeur à valeur se décline selon le niveau.

Les opérateurs + - * / (et **), aussi bien de zone à zone (merge_components) que d’agrégat à agrégat (merge_field), sont par (support, composante) en union avec passthrough — les deux opérandes n’ont pas besoin du même jeu de composantes ni de la même décomposition :

• la sortie couvre l’union des supports des deux opérandes ;
• sur un support porté des deux côtés, les zones sont fusionnées composante par composante (merge_components) : une composante définie des deux côtés devient op(a, b) ; une composante d’un seul côté passe telle quelle (passthrough brut, pour tous les opérateurs — donc a - b sur une composante propre à b vaut b, pas -b) ;
• un support porté d’un seul côté voit sa (ses) zone(s) passer inchangées.

Cela suppose l’invariant de champ (au plus une zone par (support, composante), garanti par l’union |) ; la sortie l’hérite par construction.

« Même support » = le même objet, pas la même géométrie. L’appariement se fait par identité d’objet (Handle::same_object, via SubField::same_support), jamais en comparant les nœuds. Deux supports distincts qui portent les mêmes nœuds — ou qui n’en partagent que quelques-uns, comme deux régions de bord adjacentes assemblées chacune de son côté — comptent donc chacun comme « porté d’un seul côté » : les deux zones passent inchangées et se retrouvent côte à côte dans la sortie. Rien n’est sommé au nœud commun, et à la lecture agrégat (gather, value) c’est la première zone définissant (nœud, composante) qui l’emporte — l’autre contribution est perdue sans erreur. Pour additionner deux régions qui se touchent, il faut donc les ramener sur un support partagé au préalable : restrict(a, m) + restrict(b, m) (même maillage m ⇒ même support POI1 canonique, mis en cache), ou restrict_like(b, a) pour retomber sur celui de a (cf. Opérateurs sur les champs).

Tous les opérateurs passent par merge_components (zone) / merge_field (agrégat) / merge_subfield (maj ciblée d’une zone). Là où un jeu de composantes divergent est une erreur plutôt qu’un passthrough (p. ex. l’interpolation Evolution entre deux valeurs tabulées du même champ), on garde merge_components mais précédé du garde-fou SubField::check_same_components (même support et même jeu de composantes, sinon erreur).

La même mécanique vaut pour field ** field (puissance élément par élément, exposant pris dans le second champ). La division — et la puissance à exposant fractionnaire sur base négative — ne se protègent pas des cas limites à ce niveau (sémantique numpy : inf / nan).

4. Fonctions unaires (cos, exp, …)

Des fonctions mathématiques de base s’appliquent élément par élément, renvoyant un nouveau champ de même type (zone ou agrégat). Style numpy, exposées au top-level côté Python :

import pyrucast as pc

champ2 = pc.field.cos(champ1)  # cosine of every value
e = pc.field.exp(pc.field.abs(u) * -1.0)  # they compose freely
norme = pc.field.sqrt(sx**2.0 + sy**2.0)

Jeu disponible : abs, sqrt, exp, log (népérien), log10, cos, sin, tan, sinh, cosh, tanh. Sémantique non gardée comme le reste (log d’un négatif → nan). Côté Rust ce sont des fonctions nommées (ops::field::cos(&f), …, génériques via le trait MapValues) ; il n’y a pas de syntaxe cos(x) pour un opérateur en Rust, donc seul Python en profite à l’écriture. Tout repose sur la primitive map_all (cf. plus haut), sans logique nouvelle.

Interface (résumé)

Les deux familles concrètes ajoutent leurs accès indexés et leurs constructeurs propres :

• NodeField / SubNodeField — valeurs par nœud, lecture agrégat field.value(node, "c"), écriture par zone ;
• ElementField / SubElementField — valeurs par (cellule, point de Gauss).

Et l’union | ? L’union compose des zones (structure), l’arithmétique combine des valeurs. Pour un ElementField, l’union ne fusionne plus les zones : elle valide simplement qu’aucune composante n’est portée deux fois sur le même support (deux zones de même support à composantes disjointes restent côte à côte). Pour un NodeField, l’union finalise encore en fusionnant les zones de même support (et lève si elles divergent sur une valeur partagée). La fusion explicite reste offerte par node_field.consolidate / element_field.consolidate (cf. Opérateurs sur les champs).

Champ aux nœuds (NodeField / SubNodeField)

Un champ aux nœuds porte une ou plusieurs valeurs par nœud. Il suit la même grammaire d’agrégat que tous les conteneurs de pyrucast (cf. Conventions) :

• SubNodeField — les valeurs d’une zone : un bloc multi-composantes sur les nœuds d’un sous-maillage POI1 (cf. Maillage) ;
• NodeField — l’agrégat : une liste de SubNodeField, un par zone, avec éventuellement des composantes différentes d’une zone à l’autre.

C’est le miroir exact de ElementField / SubElementField côté valeurs aux nœuds. L’intérêt de l’agrégat : un champ multiphysique (par exemple T sur tout le domaine, UX/UY sur la zone solide seulement) se représente sans inventer de 0.0 pour les couples (nœud, composante) qu’aucune zone ne définit — rien n’est densifié.

Support : un sous-maillage POI1 par zone

Un SubMesh POI1 est, par construction, exactement une liste de nœuds (un nœud par cellule). Chaque SubNodeField s’appuie sur un support POI1 :

• construit depuis un SubMesh POI1, le champ partage le handle du support tel quel (aucun refcount par nœud supplémentaire : le SubMesh est l’unique propriétaire des increfs, garder son handle suffit à garder les nœuds en vie) ;
• construit depuis un SubMesh d’un autre type d’élément, le champ se pose sur le compagnon POI1 canonique de la zone — ses nœuds distincts dans l’ordre de première apparition —, matérialisé une fois et mémoïsé par le sous-maillage. Deux champs bâtis sur la même zone tombent donc sur le même support, et s’apparient (cf. same_support).

Le champ ne stocke aucun identifiant de nœud : un POI1 est la liste de nœuds, et c’est le support qui la détient, une fois pour tous les champs qui s’y posent. Recopier la liste dans chaque champ reviendrait à la stocker autant de fois qu’il y a de champs (et de pas de temps) sur ce support ; le champ ne porte que ses valeurs, et lit les nœuds dans le support quand il en a besoin.

Scellement : le compagnon, pas le maillage

C’est le support qui est scellé — le compagnon POI1 —, jamais le maillage donné en argument. Le champ indexe ses lignes par position dans ce support, donc ce support ne doit plus bouger ; mais la zone d’origine, elle, n’a aucune raison de geler : le champ n’y touche plus.

Modifier cette zone (add_cell, remap_nodes) lui fait lâcher son compagnon : le nuage de nœuds a changé, le cache n’y répond plus. Conséquence, un champ construit après la modification se pose sur un nouveau support et ne s’apparie plus avec les précédents. Ce n’est pas une invalidation : les champs d’avant gardent leur support en vie et restent parfaitement lisibles, simplement définis sur le maillage d’avant. Pour les ramener sur le nouveau, un restrict suffit.

   NodeField (agrégat)
   ├── SubNodeField zone 0 ── support POI1 ── values[i × ncomp + c]
   ├── SubNodeField zone 1 ── support POI1 ── values[...]
   └── …

Composantes nommées

Chaque zone porte ses noms de composantes ("UX", "UY", "T", …), rangés en row-major : la composante c du nœud i est à l’indice i × ncomp + c. Au moins une composante par zone, noms uniques, valeurs initialisées à 0.0. Au niveau agrégat, components() renvoie l’union des composantes des zones (ordre de première apparition).

Les caractéristiques communes à tous les champs (composantes, min, max, sum, arithmétique scalaire et par composante) sont portées par les traits Rust SubField (niveau zone) et Field (niveau agrégat) — partagés avec ElementField.

Nœuds d’interface : duplication, lecture, cohérence

Un nœud partagé par plusieurs zones (nœud d’interface) est stocké une fois par zone. Trois règles régissent cette duplication :

• lecture agrégat (field.value(nœud, comp)) : la première zone définissant le couple gagne — aucune vérification au fil de l’eau. La forme par lot field.values(nœuds, comp) lit une liste de valeurs dans le même ordre : nœuds est une liste de nœuds, un SubMesh POI1, ou un Mesh POI1 (ses points pris dans l’ordre de la connectivité) ; même règle « première zone » et même erreur qu’un nœud non défini ;
• écriture : il n’y a pas d’écriture au niveau agrégat ; toute mutation passe par les zones (field[i]), exactement comme ElementField ;
• cohérence à la demande : field.check() vérifie que toutes les zones stockant un même couple (nœud, composante) portent la même valeur (comparaison exacte) ; node_field.consolidate(field) fait cette vérification puis fusionne par support — les zones définies sur le même SubMesh (identité de handle) deviennent une seule zone portant l’union de leurs composantes (valeurs des composantes communes vérifiées), les supports distincts restent séparés.

Composition : union |

Comme pour tous les agrégats, a | b unit les zones (handles partagés, pas de copie ; déduplication par handle) — ce n’est pas une addition de valeurs. L’union finalise en fusionnant les zones de même support (voir node_field.consolidate ci-dessus) et lève si deux zones divergent sur une valeur partagée. L’arithmétique scalaire (f + 2.0, f * 0.5, …) vit au niveau zone (SubNodeField) sur +/*/… Le nommé merge(a, b) ≡ a | b.

API Rust

#[test]
fn un_champ_aux_noeuds_s_ecrit_par_zone_et_se_lit_par_agregat() -> 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(), &[1.0, 0.0])?;

    // Support : SubMesh POI1 contenant a et b.
    let sm = {
        let mut sm = SubMesh::new(coords.clone(), ElementType::POI1);
        sm.add_cell(&[a.id()])?;
        sm.add_cell(&[b.id()])?;
        Handle::new(sm)
    };

    // Champ de déplacement 2D mono-zone : composantes UX, UY.
    let u = NodeField::from_submesh(&sm, vec!["UX".into(), "UY".into()])?;

    // Écriture : via la zone. Lecture : via l'agrégat (ou la zone).
    u.get(0)?.write().set_value(a.id(), "UX", 1.5)?;
    assert_eq!(u.value(a.id(), "UX")?, 1.5);
    assert_eq!(u.value(b.id(), "UX")?, 0.0); // valeur par défaut

    // From a multi-zone mesh: one SubNodeField per submesh.
    let mesh = Mesh::from_submesh(SubMesh::new(coords, ElementType::POI1));
    let field = NodeField::new(&mesh, vec!["T".into()])?;
    assert_eq!(field.len(), mesh.len());
    field.check()?; // zones cohérentes aux interfaces
    Ok(())
}

API Python

import pyrucast

c = pyrucast.Coords(dim=2)
a = c.add_node([0.0, 0.0])
b = c.add_node([1.0, 0.0])

mesh = pyrucast.Mesh(c, "POI1")
mesh.unit().add_cell([a])
mesh.unit().add_cell([b])

# One SubNodeField per submesh of the support (Mesh or SubMesh).
u = pyrucast.NodeField(mesh, ["UX", "UY"])
print(u)  # NodeField: 1 subfield(s)
print(u.unit())  # SubNodeField: 2 node(s), 2 component(s) [UX, UY]

# Écriture via la zone, lecture via l'agrégat.
u[0][a, "UX"] = 1.5
print(u.value(a, "UX"))  # 1.5

# Batch read: a list of nodes (or a POI1 Mesh/SubMesh) → an ordered list.
print(u.values([a, b], "UX"))  # [1.5, 0.0]
print(u.values(mesh, "UX"))  # [1.5, 0.0]  — points of the POI1 mesh
print(u.min("UX"), u.max("UX"))  # 0.0 1.5
print(u.sum("UX"))  # 1.5  — Σ over the nodes (resultant of a force field)

# Components per zone (multiphysics):
f = pyrucast.NodeField.with_components_per_submesh(two_zone_mesh, [["T"], ["UX", "UY"]])
print(f.components())  # ['T', 'UX', 'UY']
f.check()  # interface coherence (raises otherwise)
g = pyrucast.node_field.consolidate(f)  # merge as tightly as possible

Refcount et durée de vie

Le champ ne fait aucune comptabilité par nœud : il garde un clone du Handle<SubMesh> de son support, et c’est le SubMesh qui possède les increfs par nœud dans la Coords (cf. Coords). Tant qu’une zone du champ est vivante, son support l’est aussi, donc ses nœuds aussi — même si tous les Node utilisateurs ont disparu.

La libération est automatique : quand le dernier Handle sur une zone disparaît, la zone est détruite, son clone du handle de support avec elle, et les nœuds redeviennent collectables si plus rien ne les retient (cf. Modèle mémoire).

Opérateurs consommant un champ

Les opérateurs détaillés sont décrits dans Opérateurs sur les champs (dérivations), Assemblage (flux, second membre) et Solveur. Tous consomment l’agrégat et résolvent les nœuds à travers les zones (règle premier-trouvé) :

Champ aux points de Gauss (ElementField / SubElementField)

Un champ aux points de Gauss porte une ou plusieurs valeurs par (cellule, point de Gauss) sur un espace éléments finis. Il suit la même grammaire d’agrégat que tous les conteneurs de pyrucast (cf. Agrégat) et le contrat de champ :

• SubElementField — les valeurs d’une zone : un bloc multi-composantes sur les cellules × points de Gauss d’un SubFiniteElementSpace ;
• ElementField — l’agrégat : une liste de SubElementField, un par sous-espace, avec éventuellement des composantes différentes d’une zone à l’autre.

C’est le miroir exact de NodeField côté points d’intégration. C’est l’objet sur lequel s’écrivent naturellement :

• les propriétés matériau (module d’Young, Poisson, conductivité, masse volumique…) évaluées là où les intégrales sont calculées ;
• les variables internes (déformation plastique, endommagement…) gardées de cellule en cellule et de point en point ;
• les grandeurs dérivées d’une solution (contraintes, déformations, flux…) pour le post-traitement, produites par les opérateurs de champ (gradient, deformation…) et le comportement (integrate_behavior).

Support : un sous-espace éléments finis par zone

Chaque SubElementField est attaché à un seul SubFiniteElementSpace (cf. Espace éléments finis), qui détermine :

• la liste des cellules concernées (via son SubMesh) ;
• le nombre de points de Gauss par cellule (via sa QuadratureRule).

Les trois dimensions d’une zone sont figées à la construction : cell_count (du SubMesh), gauss_count (de la quadrature), component_count (choisi par l’utilisateur). Le buffer interne est dimensionné une fois pour toutes et n’est jamais réalloué — la topologie du maillage sous-jacent doit rester figée pour la durée de vie du champ (les coordonnées, elles, peuvent évoluer ; cf. FiniteElementSpace).

Les coordonnées et poids des points de Gauss ne sont pas stockés dans le champ : ils restent sur le SubFiniteElementSpace comme données de référence.

   ElementField (agrégat)
   ├── SubElementField zone 0 ── support SubFiniteElementSpace ── values[…]
   ├── SubElementField zone 1 ── support SubFiniteElementSpace ── values[…]
   └── …

Composantes nommées

Chaque zone porte ses noms de composantes ("E", "nu", "sigma_xx", "plastic_strain"…) : au moins une, noms uniques, valeurs initialisées à 0.0. Au niveau agrégat, components() renvoie l’union des composantes des zones (ordre de première apparition) ; une composante peut n’exister que sur certaines zones.

Disposition mémoire

Les valeurs d’une zone sont rangées à plat, ligne-major, dans l’ordre cellule → gauss → composante :

values[cell_idx * gauss_count * component_count
       + g * component_count
       + c]

Cet ordre rend deux accès courants cache-friendly :

• lire toutes les composantes à un point de Gauss d’une cellule (par exemple (E, nu, rho) pendant l’assemblage) — component_count flottants contigus, exposés par point_values(cell, g) ;
• balayer tous les points de Gauss d’une cellule pour une composante donnée — gauss_count flottants régulièrement espacés.

Construction : au niveau agrégat

Comme tout agrégat, un ElementField se construit au niveau parent, à partir d’un FiniteElementSpace : une zone par sous-espace.

• ElementField(fes, components) — la même liste de composantes pour chaque sous-espace ;
• ElementField.with_components_per_subspace(fes, [...]) — une liste de composantes par sous-espace (multiphysique / multi-matériau).

Pour fabriquer un champ matériau prêt pour l’assemblage, l’opérateur material_field(model, [...]) est plus direct : il crée les zones nécessaires aux sous-modèles qui consomment du matériau et les remplit en un appel.

Refcount et cycle de vie

Chaque zone détient un Handle<SubFiniteElementSpace> (cloné, donc compté). Tant qu’une zone est vivante, son sous-espace ne peut pas être collecté ; à son Drop, le refcount du sous-espace décroît et la cascade descend jusqu’au SubMesh puis à la Coords. Un ElementField n’incrémente pas le refcount des nœuds : il n’a pas de support nodal direct — les nœuds restent protégés par le SubMesh du sous-espace.

API Rust

#[test]
fn un_champ_aux_points_de_gauss_porte_le_materiau() -> 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(), &[1.0, 0.0])?;
    let c = Node::create_in(coords.clone(), &[0.0, 1.0])?;
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords, ElementType::TRI3));
    mesh.add_cell(&[a.id(), b.id(), c.id()])?;
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // 2-D linear elasticity: two material properties, one zone (one subspace).
    let mat = ElementField::new(&fes, vec!["E".into(), "nu".into()])?;
    {
        let mut z = mat.get(0)?.write(); // la zone (SubElementField) — guard
        z.set_uniform("E", 210e9)?; // module d'Young constant
        z.set_uniform("nu", 0.3)?; // Poisson constant
        assert_eq!(z.value(0, 0, "E")?, 210e9);
    }

    // Components per subspace (multi-material):
    let mat2 = ElementField::with(
        &fes,
        &[vec!["E".into(), "nu".into()]], // une liste par sous-espace
    )?;

    // Statistics and arithmetic at the aggregate level.
    assert_eq!(Field::max(&mat, Some("E"))?, 210e9);
    let scaled = &mat * 1.1; // nouveau champ (référence : préserve `mat`)
    mat.mul_to_component("E", 0.95)?; // en place, seulement "E"
    let _ = (mat2, scaled);
    Ok(())
}

API Python

import pyrucast

# Maillage + FE space — préparation.
c = pyrucast.Coords(dim=2)
a = c.add_node([0.0, 0.0])
b = c.add_node([1.0, 0.0])
c2 = c.add_node([0.0, 1.0])
mesh = pyrucast.Mesh(c, "TRI3")
mesh.unit().add_cell([a, b, c2])
fes = pyrucast.FiniteElementSpace(mesh)

# Material field: one zone per subspace of `fes`.
mat = pyrucast.ElementField(fes, ["E", "nu"])
print(mat)  # ElementField: 1 subfield(s)
print(mat.unit())  # SubElementField: 1 cell(s) × 3 gauss × 2 component(s) [E, nu]

# Writing through the zone; reading through the zone (or the aggregate stats).
z = mat.unit()  # the only zone (an error if there were several)
z.set_uniform("E", 210e9)
z.set_uniform("nu", 0.3)
assert z.value(0, 0, "E") == 210e9

# Dictionary-like access on the zone — `sub[cell, gauss, "name"]`.
z[0, 2, "nu"] = 0.28
assert z[0, 2, "nu"] == 0.28

# Stats and arithmetic at the aggregate level.
print(mat.min("E"), mat.max("E"))  # 210000000000.0 210000000000.0
print(mat.sum("E"))  # Σ over the Gauss points
mat.mul_to_component("E", 0.95)  # en place, seulement "E"
scaled = mat * 1.1  # nouveau champ

# Components per subspace (multiphysics / multi-material).
ef = pyrucast.ElementField.with_components_per_subspace(fes, [["E", "nu"]])
print(ef.components())  # ['E', 'nu']

Le plus souvent, on ne construit pas le champ matériau à la main : on appelle l’opérateur material_field qui apparie les zones aux sous-modèles consommateurs de matériau et ignore les autres (Dirichlet, …).

Visualisation

element_field.plot(...) colore le champ sur son propre support : chaque zone retrouve son sous-maillage via son sous-espace EF (partagé, pas copié). Le rendu raisonne par élément — les valeurs nodales viennent d’un moindre carré local à l’élément sur les valeurs de Gauss, sans moyenne inter-éléments, de sorte que les discontinuités (flux, contraintes) restent visibles. Voir Visualisation.

Sérialisation

SubElementField (et donc ElementField) implémente Portable via serde comme tous les objets pyrucast : le buffer de valeurs et la liste de noms voyagent dans le format binaire portable Linux ↔ Windows. Le lien vers l’espace EF, lui, est un Handle : il devient un identifiant local au fichier, et l’espace est écrit avec le champ — voir Sauvegarde et relecture.

Limitations actuelles

• Pas de mécanisme de rééchantillonnage entre quadratures (« projeter ce champ Gauss-2-points sur un autre Gauss-3-points ») : le sous-espace est figé à la création de chaque zone.
• L’arithmétique binaire entre champs est stricte (même support, mêmes composantes ; cf. Champ) — il n’y a pas (encore) de combinaison tolérante avec rééchantillonnage ou complétion par zéro.

Modèle physique (Model)

Un Model est l’objet orchestrateur qui décrit un problème physique et produit les matrices (raideur, masse) à la demande de l’utilisateur. Il vient se poser sur la couche éléments finis (FiniteElementSpace), elle-même posée sur le maillage géométrique.

Géométrie         Mesh / SubMesh                       (purement géométrique)
Formulation EF    FiniteElementSpace / SubFiniteElementSpace      (interpolation, quadrature)
Physique          Model / SubModel / SubModelKind      (loi, matériaux, assemblage)

Architecture

Model
├── sub_models: Vec<SubModel>
├── primal_vars(): Vec<String>      # union — colonnes des matrices
├── dual_vars():   Vec<String>      # union — lignes des matrices
├── filter(Physics) -> Model        # sous-modèles d'une nature donnée
└── fespace() -> FiniteElementSpace # 1 sous-espace par sous-modèle de domaine
                                     # (contraintes exclues, sans dédup)

ops::matrix (opérateurs, pas des méthodes de Model)
├── stiffness(model, materials) -> Matrix   # K  (assemblé sur demande)
└── mass(model)                 -> Matrix   # M  (assemblé sur demande)

SubModel  (énum de stockage + dispatch — AUCUNE logique)
├── HeatConduction(HeatConduction)    # chaque variante enveloppe une struct…
├── Dirichlet(Dirichlet)             # …qui porte ses données + impl SubModelKind
├── Mpc(Mpc)                         # contrainte multi-points (relations linéaires)
├── fespace() -> Option<SubFiniteElementSpace>  # sous-espace intégré (None si contrainte)
└── as_kind(&self) -> &dyn SubModelKind   # l'unique match du module modèle

SubModelKind  (trait de base — le dénominateur commun, co-localisé par physique)
├── primal_vars / dual_vars
├── physics       # ensemble de natures : &[Physics] (Mechanical|Thermal|Constraint|Other)
├── as_domain / as_constraint  # seams de capacité (None par défaut) — cf. ci-dessous
├── matrix_element          # pont vers les noyaux élémentaires, qui vivent sur Domain
├── stiffness_layout        # Some ⇒ bloc CALCULÉ (scatter parallèle) ; None ⇒ littéral
├── contributions           # défaut : dérivé du layout ; contraintes rendent leurs C/Cᵀ littéraux
├── build_stiffness_blocks  # défaut : stiffness_layout + Domain::element_matrix
├── build_mass_blocks       # (défaut : vide)
└── label / display / render

Capacités (sous-traits, miroir des natures ; une struct n'a que la sienne) :
├── Domain      # matériau + comportement (heat, elasticity, poutres, …)
└── Constraint  # multiplicateurs de Lagrange (Dirichlet, MPC, embedded, contact)

L’énum SubModel ne sert qu’au stockage et à la sérialisation (bincode) ; il délègue chaque appel à l’impl SubModelKind de la variante via as_kind(). Tout le code générique (l’agrégat Model, l’assembleur, Dump) passe par ce seul point — ajouter une physique ne touche donc aucun de ces sites. Voir le chapitre Ajouter une physique.

Le Model est purement orchestrateur : il énumère les DOFs, dimensionne la Matrix, boucle sur les sub-models et accumule. Aucune logique physique ne vit chez lui.

Identification des DOFs : primal ≠ dual

Chaque physique déclare :

• ses variables primales (les inconnues — composantes du vecteur solution, colonnes de la matrice) ;
• ses variables duales (les conjuguées énergétiques — composantes du vecteur chargement, lignes de la matrice).

Elles sont presque toujours différentes :

Les DOFs de la Matrix sont identifiés par le couple (NodeID, nom_de_champ) (voir Matrix) : deux SubModels qui utilisent le même nom ("T") sur des nœuds différents ne se collisionnent pas, et la jonction se fait automatiquement quand ils partagent un même (NodeID, nom).

Chargements complètement séparés du Model

Le Model ne porte aucune logique de second membre. L’utilisateur :

1. lit model.dual_vars() pour connaître les noms de composantes du vecteur force ;
2. construit un NodeField avec ces composantes (forces de Neumann, sources de chaleur, valeurs imposées de Dirichlet aux nœuds-multiplicateurs, …) ;
3. compose plusieurs sources avec | (union des zones, dédupliquée et fusionnée par support) — le nommé merge en est l’alias ;
4. passe Matrix + NodeField au solveur.

Cette séparation a deux mérites :

• les chargements sont des données utilisateur, faciles à composer ;
• le Model reste une description compacte et indépendante du chargement (le même modèle peut être résolu avec plusieurs chargements en cascade).

Pour la part contrainte du second membre (les valeurs imposées aux nœuds-multiplicateurs), le helper model.constraint_rhs([(nœud, g), …]) construit ce NodeField tout seul : on désigne chaque relation par un nœud contraint (Dirichlet) ou un nœud-terme (MPC) et sa valeur g, et le helper retrouve le nœud-multiplicateur et la composante à renseigner (imposed_<v>, mpc_rhs). Voir Contraintes.

Les physiques disponibles

Chaque physique est une struct sous src/models/ implémentant le trait SubModelKind, enveloppée par une variante de l’énum SubModel. Leur détail (équations, matériau, comportement, exemples) est dans la partie Détails des physiques ; on n’en rappelle ici que les opérateurs qui les déclarent et leurs variables, vus du Model :

Toutes balaient tous les sous-espaces du fes (une zone par sous-espace), sauf dirichlet, mpc, embedded et contact qui sont des contraintes portées par des maillages fournis par l’utilisateur. Le matériau est toujours fourni à l’assemblage, pas au modèle (cf. ci-dessous).

Pour ajouter une physique, voir Ajouter une physique.

Ce que chaque physique calcule

Le tableau ci-dessus dit quelles variables porte chaque physique. Ce qu’elle sait produire — les genres de matrice qu’elle déclare, la voie par laquelle elle obtient sa tangente, son intégration de comportement et sa particularité de calcul — est rassemblé physique par physique dans Détails des physiques.

Nature physique et filtrage

Chaque physique déclare un ensemble de natures — sa classification grossière, orthogonale à l’axe de capacité Domain/Constraint. Elle répond à « quel champ de physique » là où les capacités répondent à « domaine ou contrainte » :

Côté Python, les mêmes natures sont des chaînes : "mechanical", "thermal", "constraint", "other", "diffusion", "radiation".

Diffusion est une nature à part entière bien que la loi de Fick partage l’opérateur de la conduction : les variables diffèrent (c/j contre T/q), et un problème couplé doit pouvoir sélectionner l’une sans l’autre. Partager un opérateur n’est pas partager une physique.

La nature d’une physique de base est entièrement déterminée par la variante : c’est une constante par physique, exposée par SubModelKind::physics() — un slice &'static [Physics] (comme label()), pas un champ stocké. Le type est un ensemble pour deux raisons :

• une physique couplée (par ex. un futur élément thermo-mécanique) porte plusieurs natures — [Mechanical, Thermal] ;
• un bloc de matrice monté à la main, hors assemblage, n’en porte aucune — l’ensemble vide, le cas « rien ». Physics::Other est la nature « autre » explicite, pour un bloc qu’on veut classer plutôt que laisser sans étiquette.

L’ensemble voyage avec chaque bloc assemblé jusqu’à la SubMatrix (posé par l’assembleur sur les deux chemins, calculé et littéral, donc le couple C/Cᵀ d’un Dirichlet est étiqueté aussi).

Deux sélecteurs symétriques en découlent — ils gardent les entités dont l’ensemble contient la nature (une physique couplée apparaît donc sous chacune) — tous deux à partage par compteur de références (pas de copie profonde) :

• model.filter(Physics::Mechanical) → un Model ne gardant que les sous-modèles au moins mécaniques ;
• k.filter(Physics::Mechanical) → une Matrix ne gardant que les blocs au moins mécaniques (non assemblée — relancer Matrix::assemble avant de résoudre).

k.physics() renvoie l’ensemble des natures présentes dans la matrice (dédupliqué) : une matrice agrégeant plusieurs physiques y expose plusieurs tags (par ex. [Thermal, Constraint]). Un bloc à l’ensemble vide n’est jamais sélectionné par une nature concrète — l’étiqueter Physics::Other le rend atteignable par filter(Physics::Other).

#[test]
fn filtrer_un_modele_et_sa_matrice_par_nature() -> Result<()> {
    let coords = Handle::new(Coords::new(1)?);
    let a = Node::create_in(coords.clone(), &[0.0])?;
    let b = Node::create_in(coords.clone(), &[1.0])?;
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    mesh.add_cell(&[a.id(), b.id()])?;
    let fes = FiniteElementSpace::lagrange1(&mesh)?;
    let model = model::heat_conduction(&fes)?;
    let materials = element_field::material_field(&model, &[("k", 1.0)])?;
    let k = matrix::stiffness(&model, &materials)?;

    let meca = model.filter(Physics::Mechanical); // sous-modèles au moins mécaniques
    let k_meca = k.filter(Physics::Mechanical); // blocs au moins mécaniques (non assemblés)
    let natures = k.physics(); // ex. [Thermal, Constraint]

    assert!(meca.is_empty() && k_meca.is_empty()); // ce modèle est thermique
    assert!(natures.contains(&Physics::Thermal));
    Ok(())
}

Règle invariante : un Model = une Matrice

matrix::stiffness(model, materials) et matrix::mass(model, materials) produisent chacune une seule Matrix couvrant l’ensemble des DOFs du Model (primaux ⊕ multiplicateurs). Les conditions limites n’ont pas de statut spécial — ce sont des sub-models comme les autres qui contribuent leurs entrées dans la même matrice globale.

Cette uniformité simplifie tout : le solveur reçoit une seule Matrix + un seul NodeField ; pas besoin de jongler avec un système saddle-point composé.

API Rust

#[test]
fn un_modele_se_declare_et_s_assemble() -> Result<()> {
    // 1-D: a [0, 1] mesh with a single SEG2.
    let coords = Handle::new(Coords::new(1)?);
    let a = Node::create_in(coords.clone(), &[0.0])?;
    let b = Node::create_in(coords.clone(), &[1.0])?;
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    mesh.add_cell(&[a.id(), b.id()])?;
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // Model: conduction (the material is supplied at assembly, not here) +
    // Dirichlet on the left. Constructors at the parent level (they sweep
    // `fes`'s subspaces), composed with `union` — a `SubModel` is never built by
    // hand (see CONVENTIONS.md).
    let hc = model::heat_conduction(&fes)?;
    // Mesh of the imposed nodes + support of the multipliers (barycenter
    // co-locates fresh nodes). The model creates no node itself.
    let imposed = mesh::poi1_from_nodes(std::slice::from_ref(&a))?;
    let multiplier = mesh::barycenter(&imposed)?;
    let dir = model::dirichlet(&hc, "T", &imposed, &multiplier, RelationSense::Equality)?;
    let model = hc.union(&dir)?;

    // Material k = 1, applied to the sub-models that need it (Dirichlet is
    // skipped automatically), then assembly.
    let materials = element_field::material_field(&model, &[("k", 1.0)])?;
    let k = matrix::stiffness(&model, &materials)?;
    assert_eq!(k.n_rows()?, 3); // 2 nœuds physiques + 1 multiplicateur
    Ok(())
}

API Python

import pyrucast

c = pyrucast.Coords(dim=1)
a = c.add_node([0.0])
b = c.add_node([1.0])
mesh = pyrucast.Mesh(c, "SEG2")
mesh.unit().add_cell([a, b])
fes = pyrucast.FiniteElementSpace(mesh)

# Model: conduction (material supplied at assembly) + Dirichlet on the left.
# Constructors at the parent level, composed with `|` — no SubModel by hand.
# The multipliers' mesh is built from the imposed nodes.
imposed = pyrucast.mesh.poi1_from_nodes([a])
multiplier = pyrucast.mesh.barycenter(imposed)
cible = pyrucast.model.heat_conduction(fes)

model = cible | pyrucast.model.dirichlet(cible, "T", imposed, multiplier)

# Material k = 1 (the Dirichlet sub-models are skipped automatically).
materials = pyrucast.element_field.material_field(model, [("k", 1.0)])

K = pyrucast.matrix.stiffness(model, materials)
print("primal_vars =", model.primal_vars())  # ['T', 'lambda_T']
print("dual_vars =", model.dual_vars())  # ['q', 'imposed_T']
print(K)  # Matrix: 3 row(s) × 3 col(s), …

Assemblage et résolution

L’assemblage (stiffness / mass) et la résolution (solve) sont des opérateurs : ils consomment le Model (et le matériau, le chargement) et sont décrits dans la partie Détail des opérateurs — Assemblage et Solveur. Le solveur est une LU creuse directe (faer), dont la factorisation est mise en cache sur la Matrix — factoriser une fois, résoudre souvent.

Des exemples complets et à solution analytique (assemblage + contraintes + lecture des inconnues et des réactions) sont déroulés dans Dirichlet (Poisson 1-D), Conduction thermique et Mécanique.

Limitations actuelles

• Physiques disponibles : HeatConduction et BoundaryTransfer (échange de surface / film) (thermique) ; Truss, Elasticity, Plasticity, Mazars, Timoshenko, Frame, Frame3d (mécanique) ; et les contraintes Dirichlet, Mpc, Embedded, Contact (contraintes). Toute nouvelle physique est une struct implémentant SubModelKind (une variante de l’énum SubModel
  • un bras de as_kind, rien d’autre — cf. Ajouter une physique). Le coût d’ajout est O(1) fichier, indépendant du nombre de physiques existantes.
• Toutes les physiques n’ont pas tous les genres de matrice : chacune déclare les MatrixKind qu’elle sait produire (raideur, masse, raideur géométrique, tangente cohérente). Assembler un genre qu’une physique n’a pas ne casse rien — elle ne contribue simplement pas.
• Pas de check de cohérence pré-assemblage : la consistance (matériau définit bien "k" pour HeatConduction, compatibilité des FE spaces entre sub-models, etc.) est vérifiée au moment de matrix.stiffness / matrix.mass, pas à l’ajout du sub-model. Si on découvre des cas où ça pose problème, un check eager est facile à ajouter.
• Deux back-ends de solveur, tous deux directs : LU creuse (défaut) et Cholesky, au choix de l’appelant via method=, toutes deux en faer et avec cache de factorisation. Cholesky exige une matrice qui se déclare symétrique — un triangle ne se lisant que d’un côté, elle ne se tromperait pas bruyamment sur une matrice qui ne l’est pas — et définie positive : un point-selle à multiplicateurs, symétrique mais indéfini, est refusé au pivot fautif, qu’on l’élimine d’abord ou qu’on le résolve en LU. Pas de méthode itérative ; SolveMethod reste le point d’extension prévu pour cela.

Matrice creuse (Matrix)

Matrix est le conteneur de sortie d’un assemblage : c’est ce que produisent les opérateurs matrix::stiffness(model, materials) / matrix::mass(model, materials) à partir d’un Model. Elle représente une matrice creuse dont les lignes et les colonnes sont identifiées par des DOFs nommés.

Identification des DOFs : (NodeId, nom_de_champ)

Chaque ligne et chaque colonne d’une Matrix est identifiée par un couple (NodeId, ChampId) :

• NodeId — l’identifiant stable d’un nœud dans la Coords.
• ChampId — un indice compact dans une petite table de noms portée par la matrice (typiquement 5–10 entrées). Les noms sont des chaînes comme "T", "q", "ux", "lambda_w".

Le type concret est DofId { node_id, field_idx }. Cette représentation est compacte (un u32 par champ partagé sur tous les DOFs qui le portent) et conserve l’information sémantique : à chaque entrée numérique de la matrice est attaché « quel inconnu, à quel nœud ».

Les jeux de DOFs de lignes et de colonnes sont indépendants :

• ils peuvent avoir des tailles différentes (matrice rectangulaire — par exemple le bloc Lagrange d’une condition de Dirichlet) ;
• ils peuvent porter des noms de champs différents (les lignes étiquetées par des duales q, les colonnes par des primales T).

Blocs bi-mode : littéral ou calculé

Une Matrix est un agrégat de blocs SubMatrix, et un bloc est de l’un de deux modes :

• littéral — il porte ses valeurs, stockées en COO (coordinate triplet list). Chaque add_entry(...) ajoute un triplet (ligne, colonne, valeur) ; plusieurs entrées au même couple s’accumulent (sommées à l’assemblage), l’ordre d’insertion étant sans effet. C’est le mode historique — celui des contraintes (blocs C / Cᵀ de Dirichlet) et de tout bloc monté à la main.
• calculé — il ne porte aucune valeur, seulement une recette { sous-modèle, sous-espace EF, matériau }. Ses entrées sont produites à l’assemblage par le noyau élémentaire du sous-modèle, dispersées directement dans la matrice globale. C’est le mode des physiques volumiques (raideur), qui évite de matérialiser un COO intermédiaire.

Un bloc calculé garde son lien vers sa physique via la recette ; la Matrix, elle, reste un simple sac de blocs et ne référence pas le Model.

Le bloc ne recopie pas sa liste de nœuds

Un bloc est posé sur deux supports POI1 (lignes et colonnes, souvent le même objet) et n’en garde aucune copie : il lit leur connectivité en place à chaque accès, conformément à la règle Zéro-copie. C’est sûr parce que les deux supports sont scellés à la construction du bloc — leur connectivité ne peut plus changer, donc la numérotation ne peut pas dériver. Le NodeId → position passe par la table que le support porte déjà (SubMesh::node_index), partagée avec tous ses autres consommateurs au lieu d’être refaite par bloc.

Corollaire à connaître si l’on monte un bloc à la main : les nœuds d’un support doivent être distincts, ce que produit to_poi1. Un nœud répété n’est pas rejeté, mais il adresse la mauvaise ligne — la table du support donne un rang dédoublonné, qui s’écarte de la position dès la première répétition.

Étiquette de nature physique (physics)

Chaque bloc porte en plus un ensemble de natures Vec<Physics> (Mechanical, Thermal, Constraint, Other) — l’assembleur le pose sur tout bloc qu’il émet, sur les deux chemins (calculé et littéral), donc le couple C/Cᵀ d’un Dirichlet est étiqueté lui aussi. C’est ce qui rend l’étiquette utilisable là où la recette manque (blocs littéraux). Le tag est un ensemble : vide pour un bloc monté à la main hors assemblage (le cas « rien »), et à plusieurs éléments pour une physique couplée.

Il alimente Matrix::filter(Physics) — le miroir de Model::filter — qui renvoie une Matrix ne gardant que les blocs dont l’ensemble contient la nature donnée (handles partagés, pas de copie). Le résultat n’est pas assemblé : relancer Matrix::assemble avant de résoudre. Un bloc à l’ensemble vide n’est sélectionné par aucune nature concrète ; l’étiqueter Physics::Other le rend atteignable. Matrix::physics() renvoie l’ensemble des natures présentes dans la matrice (dédupliqué — « plusieurs tags » au niveau de l’agrégat).

#[test]
fn filtrer_une_matrice_par_nature() -> Result<()> {
    let (model, materials, _, _) = barre()?;
    let k = matrix::stiffness(&model, &materials)?;

    let k_meca = k.filter(Physics::Mechanical); // blocs au moins mécaniques
    let natures = k.physics(); // ex. [Thermal, Constraint]

    assert!(k_meca.is_empty()); // ce modèle est thermique
    assert!(natures.contains(&Physics::Thermal));
    Ok(())
}

Assemblage : motif + scatter

Passer d’un agrégat de blocs à une matrice utilisable se fait en deux temps :

1. Motif creux (sparsité CSR) — l’union dédoublée des DDL des blocs (une liste globale, via table de hachage) et de leurs entrées. Il ne dépend que de la topologie (bloc calculé via la connectivité, bloc littéral via sa COO), pas des matériaux ; stiffness le mémoïse donc sur le Model et le réutilise d’un assemblage à l’autre.
2. Valeurs — dispersées (scatter) dans le CSR : un bloc calculé lance son noyau élémentaire (en parallèle, par coloration des cellules — voir Parallélisme) ; un bloc littéral recopie sa COO. Chaque bloc remappe sa numérotation locale (nœud, variable) vers l’index global via une table de traduction — O(nnz), sans recherche par entrée (NodeId est déjà l’index nœud global dense, et add_entry retrouve la position d’un nœud en O(1)).

L’ordre des DOFs dans row_dofs() / col_dofs() est l’ordre de première rencontre des blocs — sauf si le Coords porte une permutation (ordre solveur), auquel cas la liste globale suit cet ordre (tri stable). Reproductible dans les deux cas.

finalize vs ops::matrix

• Matrix::finalize() n’assemble que des blocs littéraux (somme des COO → CSR). Il refuse un bloc calculé : le noyau vit dans models, hors de containers, et l’y appeler créerait un cycle matrix ↔ kernel. Il renvoie alors vers ops::matrix.
• ops::matrix::stiffness(model, materials) construit les blocs (calculés pour les physiques volumiques, littéraux pour Dirichlet) et assemble, motif mémoïsé sur le Model.
• Matrix::assemble(&mut self) réassemble une matrice depuis ses blocs seuls, sans Model : c’est le chemin de composition — combiner une sous-matrice neuve (de provenance quelconque) à une matrice existante puis réassembler. La Matrix ne dépendant que de ses blocs, cette composabilité de base est ainsi préservée y compris en présence de blocs calculés.

Pour qui veut une matrice creuse d’une autre bibliothèque, des conversions à la demande existent — elles fabriquent un objet neuf et ne retiennent rien :

• Matrix::to_csr → nalgebra_sparse::CsrMatrix<f64>
• Matrix::to_csc → nalgebra_sparse::CscMatrix<f64>
• Matrix::to_coo → nalgebra_sparse::CooMatrix<f64>
• Matrix::to_dmatrix → nalgebra::DMatrix<f64>

Mais le solveur n’en emprunte aucune. La forme assemblée n’a pas besoin d’être convertie, seulement d’être regardée sous le bon angle : une CSR est, octet pour octet, la CSC de la transposée, et une CSC triée sans doublon est exactement ce que le LU creux de faer demande. ops::solver::lu lit donc les tableaux par Matrix::csr_arrays — empruntés, la sparsité restant celle du motif — transpose une fois par tri par comptage, et tend à faer une vue. Aucune copie de la matrice ne coexiste avec la factorisation, qui est le moment où la mémoire est la plus tendue.

Le produit matrice-vecteur, lui, tire parti de l’orientation lignes : lectures contiguës, un accumulateur par ligne, parallélisable sur les lignes sans atomique. C’est la raison pour laquelle la forme assemblée reste une CSR.

Facteur scalaire et somme de matrices

Le facteur

Chaque SubMatrix porte un facteur f64, 1.0 par défaut, ajusté par bloc * s, bloc / s et -bloc. Le facteur ne touche que ce champ — jamais les valeurs stockées (coo) — ce qui le rend utilisable aussi bien sur un bloc littéral que sur un bloc calculé (dont les valeurs n’existent qu’à l’assemblage, produites par le noyau élémentaire). Il est pris en compte partout où une valeur du bloc est lue ou émise : les accesseurs directs (get, dense, to_dmatrix, to_coo, to_csr, to_csc, mul_dense) et les deux passes d’assemblage global (Matrix::finalize et ops::matrix::scatter, calculé comme littéral). Seules les formes locales brutes (local_triplets, local_coo_arrays) restent non mises à l’échelle — ce sont des vues internes destinées au remappage global, chaque consommateur y applique le facteur lui-même.

&Matrix * s, &Matrix / s et -&Matrix mettent à l’échelle une matrice entière : chaque bloc est cloné dans un nouvel objet avec son facteur ajusté. C’est nécessaire car add_sub/union/filter/subset partagent les Handle<SubMatrix> (même objet, compté) plutôt que de les copier ; muter le facteur en place rescalerait silencieusement toute autre Matrix référençant le même bloc. La forme possédante (matrix * s sur une valeur, pas une référence) fait l’économie de cette copie pour tout bloc que personne d’autre ne tient — ce que Handle::is_sole_owner établit.

Le scalaire se lit des deux côtés (2.0 * k comme k * 2.0). Ces opérateurs sont infaillibles, à une exception près : une division refuse un diviseur nul ou non fini, qui rendrait non finie chaque valeur du résultat. Côté Rust elle interrompt l’exécution ; côté Python elle lève ZeroDivisionError ou ValueError.

La CSR assemblée suit, mise à l’échelle. Mettre tous les blocs à la même échelle met chaque entrée à cette échelle et laisse la sparsité intacte : seul le tableau des valeurs est parcouru, les tableaux d’indices et la table de noms sont des Arc partagés. Une matrice assemblée reste donc assemblée après * s, et n’a pas à repasser par un assemble() qui relancerait tous les noyaux élémentaires pour appliquer un scalaire.

(Σ v) · s n’est pas, au bit près, le Σ (v · s) que calculerait un réassemblage : l’addition flottante n’est pas associative. Les deux valent la même quantité à l’arrondi près, et chaque chemin reste reproductible.

#[test]
fn diviser_une_matrice_ne_reecrit_aucune_valeur() -> Result<()> {
    let (model, materials, _, _) = barre()?;
    let m = matrix::stiffness(&model, &materials)?;
    let a = m.row_mesh()?.node(0, 0, 0)?.id();
    let dt = 0.1;

    // Facteur 1/dt sur chaque bloc : aucune valeur stockée n'est réécrite. Et
    // `m` étant déjà assemblée, sa CSR suit, mise à l'échelle — donc pas de
    // réassemblage, donc aucun noyau élémentaire relancé pour un scalaire.
    let m_dt = &m / dt;

    assert_eq!(m.get(a, "q", a, "T"), m_dt.get(a, "q", a, "T") * dt); // m inchangée
    assert!(m_dt.to_csr().is_ok()); // utilisable telle quelle
    Ok(())
}

La somme

a + b rend une Matrix portant les blocs des deux opérandes, partagés et délibérément non dédoublonnés. Rien n’est calculé : c’est l’assembleur qui somme ce qui retombe sur le même (row, col) global (build_global_triplets, scatter_serial/scatter_parallel). Une somme coûte donc quelques incréments de compteur, ne touche aucune valeur, et laisse un bloc calculé calculé. Comme pour filter, le résultat n’est pas assemblé : assemble() avant de résoudre.

a - b nie les blocs de droite, ce qui les recopie (le facteur vit dans le bloc) ; a + b ne copie rien. Les deux opérateurs acceptent indifféremment une Matrix ou une SubMatrix de chaque côté, et -a nie une matrice entière.

#[test]
fn additionner_deux_matrices_puis_resoudre() -> Result<()> {
    let (model, materials, _, rhs) = barre()?;
    let k = matrix::stiffness(&model, &materials)?;
    let m = matrix::stiffness(&model, &materials)?;
    let dt = 0.1;

    // `+` porte les blocs des deux opérandes, partagés : rien n'est copié,
    // aucune valeur n'est touchée. C'est l'assembleur qui somme ce qui retombe
    // au même (ligne, colonne) global.
    let mut sys = &m / dt + &k;
    sys.assemble()?; // la somme n'est pas assemblée : requis avant de résoudre
    let u = solver::lu::solve(&sys, &rhs)?;

    assert!(u.node_count()? > 0);
    Ok(())
}

Aucun traitement particulier n’est nécessaire quand K et M n’ont pas le même ensemble de DOFs (cas courant : un Dirichlet/MPC n’entre que dans la matrice de raideur, jamais dans la masse) — la somme prend simplement l’union des DOFs des deux côtés, et les blocs de M ne contribuent rien aux DOFs qu’ils ne portent pas.

| compose, + additionne

C’est la seule chose qui les sépare, et elle ne se voit que sur des blocs partagés : l’union écarte un bloc dont elle tient déjà l’emplacement, la somme le compte à chaque fois qu’on le lui donne. Donc k | k vaut k, tandis que k + k vaut 2k.

Prendre | pour composer un opérateur à partir de morceaux distincts (une raideur et son bloc de Dirichlet), + pour additionner deux opérateurs.

#[test]
fn union_compose_somme_additionne() -> Result<()> {
    let (model, materials, _, _) = barre()?;
    let k = matrix::stiffness(&model, &materials)?;
    let a = k.row_mesh()?.node(0, 0, 0)?.id();
    let kaa = k.get(a, "q", a, "T");

    // `+` compte une contribution chaque fois qu'on la lui donne.
    let mut deux_fois = &k + &k;
    deux_fois.assemble()?;
    assert_eq!(deux_fois.get(a, "q", a, "T"), 2.0 * kaa);

    // `|` écarte un bloc dont il tient déjà l'emplacement : c'est l'outil pour
    // composer un opérateur à partir de morceaux distincts, pas pour additionner.
    let mut une_fois = k.union(&k)?;
    une_fois.assemble()?;
    assert_eq!(une_fois.get(a, "q", a, "T"), kaa);

    // `-` nie les blocs de droite (ce qui les recopie) ; `+` ne copie rien.
    let mut nulle = &k - &k;
    nulle.assemble()?;
    assert_eq!(nulle.get(a, "q", a, "T"), 0.0);
    Ok(())
}

Symétrie

Le dernier argument des constructeurs de SubMatrix déclare quelle part de la symétrie de la matrice ce bloc porte :

La propriété visée est celle du tableau assemblé — A[i][j] == A[j][i] sur la CSR — et rien d’autre. Un bloc la déclare et on le croit : un modèle sait ce qu’il écrit, et rien ici ne le vérifie. Déclarer juste est donc tout le travail du producteur, et l’agrégat additionne les déclarations sans les corriger. Un bloc vide, par exemple, est symétrique : il déclare Full, et la règle n’a pas d’exception à prévoir pour lui.

La numérotation suit la déclaration

Un tableau n’est symétrique que si le rang i désigne des DDL conjugués des deux côtés. Or les deux ordres globaux se construisaient par deux parcours indépendants des blocs, et rien ne les faisait tomber d’accord : une contrainte qui introduit deux nœuds neufs d’un coup — embedded, dont le nœud immergé n’appartient à aucune physique — les faisait découvrir en ordre inverse de chaque côté.

Quand la matrice se déclare symétrique, les deux ordres sont donc construits en un seul parcours conjugué. Ce parcours n’apprend jamais que q est le dual de T : il lit seulement quelles positions se correspondent, ce que la déclaration dit déjà. Un bloc Full est carré sur un support unique, donc sa ligne k fait face à sa propre colonne k ; une paire Half croise, la ligne de l’un faisant face à la colonne de l’autre.

Deux incohérences y sont refusées, jamais rattrapées : un même DDL dual déclaré conjugué à deux DDL primaux différents, et une paire dont les deux membres n’ont pas le même nombre de DDL.

Une matrice non symétrique garde les deux parcours indépendants : une matrice rectangulaire n’a pas de conjugué à apparier.

Pourquoi une moitié

Une contrainte de Dirichlet introduit deux blocs rectangulaires, C et Cᵀ (voir Contraintes). Aucun des deux ne peut être symétrique — un bloc rectangulaire ne l’est jamais — mais ensemble ils le sont. C’est une propriété du couple, que Half rend exprimable : le producteur qui écrit le même coefficient des deux côtés est celui qui les apparie.

L’identité de la paire est une empreinte de son contenu : les nœuds des deux supports, les quatre noms de variables, le coefficient. Déterministe, donc l’archive reste reproductible ; et si deux paires réellement distinctes venaient à partager une empreinte, elles seraient rejetées, jamais acceptées à tort — on oublie une symétrie, on n’en invente pas.

Ce que l’agrégat en conclut

Matrix::symmetric() est vrai si chaque bloc est Full, ou Half avec ses deux membres présents en nombres égaux. Compter les membres plutôt que les blocs permet à une contrainte déclarée deux fois (quatre blocs, deux de chaque) de tenir, tandis qu’une paire coupée par un subset tombe.

Le stockage n’est pas dédupliqué : une matrice symétrique porte quand même ses deux triangles. Mais la déclaration, elle, est consultée : elle décide si la CSR assemblée peut être tendue telle quelle à la factorisation comme sa propre CSC. Une matrice symétrique l’est — CSR(A) est au bit près CSC(Aᵀ), et Aᵀ = A — si bien que le retournement, et la seconde copie complète de la matrice qui va avec, disparaissent. C’est elle aussi qui autorise Cholesky (voir Résolution).

Cas d’usage typique : matrice de raideur du laplacien

#[test]
fn les_entrees_vivent_dans_un_bloc() -> Result<()> {
    // The entries live in a **block**, never in the aggregate: a block knows its
    // POI1 supports (rows and columns) and its variable names.
    let coords = Handle::new(Coords::new(1)?);
    let a = Node::create_in(coords.clone(), &[0.0])?;
    let b = Node::create_in(coords.clone(), &[1.0])?;
    let support = {
        let mut sm = SubMesh::new(coords.clone(), ElementType::POI1);
        sm.add_cell(&[a.id()])?;
        sm.add_cell(&[b.id()])?;
        Handle::new(sm)
    };

    let mut block = SubMatrix::new(
        support.clone(),  // support des lignes
        support.clone(),  // support des colonnes (carré ici)
        vec!["q".into()], // variables duales   → lignes
        vec!["T".into()], // variables primales → colonnes
        DofOrdering::NodesThenVars,
        Symmetry::Full, // symétrique à lui seul
    );

    // A simple 2-node model (a segment):
    //   K = [[ 2, -1], [-1,  2]]
    block.add_entry(a.id(), "q", a.id(), "T", 2.0)?;
    block.add_entry(a.id(), "q", b.id(), "T", -1.0)?;
    block.add_entry(b.id(), "q", a.id(), "T", -1.0)?;
    block.add_entry(b.id(), "q", b.id(), "T", 2.0)?;

    let mut k = Matrix::empty();
    k.add_sub(Handle::new(block))?;
    k.finalize()?; // requis avant tout usage solveur

    assert_eq!(k.n_rows()?, 2);
    assert_eq!(k.n_cols()?, 2);
    assert!(k.symmetric());
    Ok(())
}

Matrice rectangulaire : bloc Lagrange

Une contrainte de Dirichlet introduit, par sa nature, un bloc rectangulaire : lignes indexées par les nœuds-multiplicateurs (un par contrainte), colonnes par les nœuds primaires contraints.

#[test]
fn un_bloc_de_lagrange_est_rectangulaire() -> Result<()> {
    let coords = Handle::new(Coords::new(1)?);
    let a = Node::create_in(coords.clone(), &[0.0])?;
    let b = Node::create_in(coords.clone(), &[1.0])?;
    let support = {
        let mut sm = SubMesh::new(coords.clone(), ElementType::POI1);
        sm.add_cell(&[a.id()])?;
        sm.add_cell(&[b.id()])?;
        Handle::new(sm)
    };

    // 2 constraints: the multipliers m0/m1 tie the primary nodes a/b.
    // The block is rectangular as soon as the two supports differ — here they
    // have the same size, but they are two distinct node clouds.
    let m0 = Node::create_in(coords.clone(), &[0.0])?;
    let m1 = Node::create_in(coords.clone(), &[1.0])?;
    let mult_support = {
        let mut sm = SubMesh::new(coords.clone(), ElementType::POI1);
        sm.add_cell(&[m0.id()])?;
        sm.add_cell(&[m1.id()])?;
        Handle::new(sm)
    };
    let mut block = SubMatrix::new(
        mult_support,
        support.clone(),
        vec!["T".into()],
        vec!["T".into()],
        DofOrdering::NodesThenVars,
        // Un bloc rectangulaire ne peut pas être symétrique seul ; construit à
        // la main, il n'a pas de moitié qui lui réponde.
        Symmetry::None,
    );
    block.add_entry(m0.id(), "T", a.id(), "T", 1.0)?;
    block.add_entry(m1.id(), "T", b.id(), "T", 1.0)?;

    let mut c = Matrix::empty();
    c.add_sub(Handle::new(block))?;
    c.finalize()?;
    assert_eq!(c.n_rows()?, 2);
    assert_eq!(c.n_cols()?, 2);
    // "T" is interned once only in the name table even though it appears on the
    // row side AND the column side (the collision is settled by the distinct
    // `NodeId`: the multipliers are nodes in their own right).
    assert_eq!(c.field_names().len(), 1);
    Ok(())
}

API Rust — accès en lecture

#[test]
fn lire_une_matrice_assemblee() -> Result<()> {
    let (model, materials, _, _) = barre()?;
    let k = matrix::stiffness(&model, &materials)?;
    let a: NodeId = k.row_mesh()?.node(0, 0, 0)?.id();
    let x = NodeField::from_submesh(&k.col_mesh()?.get(0)?, vec!["T".into()])?;

    // Toutes ces lectures traversent l'état assemblé : elles rendent un
    // `Result` and fail until `finalize()` (or `assemble()`) has
    // été appelé.

    // Value at a coordinate (the sum of every COO entry at that point).
    let v: f64 = k.get(a, "q", a, "T");

    // Dense row-major view (a flat Vec, handy for Python).
    let d: Vec<f64> = k.dense()?;
    assert_eq!(d.len(), k.n_rows()? * k.n_cols()?);

    // Typed nalgebra dense view (column-major DMatrix), ready for LU/Cholesky.
    let m: nalgebra::DMatrix<f64> = k.to_dmatrix()?;

    // nalgebra-sparse sparse views, ready for the sparse solvers. `to_csr`
    // materializes; `csr_arrays` borrows the three arrays without copying.
    let csr: nalgebra_sparse::CsrMatrix<f64> = k.to_csr()?;
    let csc: nalgebra_sparse::CscMatrix<f64> = k.to_csc()?;
    let (offsets, cols, vals): (&[usize], &[usize], &[f64]) = k.csr_arrays()?;
    assert_eq!(offsets.len(), k.n_rows()? + 1);
    assert_eq!(cols.len(), vals.len());

    // Iteration over the raw triplets (insertion order preserved). An entry is a
    // 5-tuple `(row node, dual var, column node, primal var, value)` — the
    // variable names are already resolved there.
    for (row_node, row_var, col_node, col_var, value) in k.iter_entries() {
        let _ = (row_node, row_var, col_node, col_var, value);
    }

    // Matrix · field product: `x` is read at the *column* DOFs (**primal** vars),
    // the result is a `NodeField` on the *row* DOFs (**dual** vars) — `K · u = f`.
    // The `*` operator is its sugar.
    let y: NodeField = k.mul_field(&x)?;
    let y_sucre: NodeField = (&k * &x)?; // le même produit, en opérateur

    let _ = (v, m, csr, csc, y, y_sucre);
    Ok(())
}

API Python

import pyrucast

# The entries live in a **block**, never in the aggregate: a block knows its
# POI1 supports (rows and columns) and its variables.
k = pyrucast.Matrix.block(support, support, ["q"], ["T"], symmetry="full")
bloc = k[0]
bloc.add_entry(a, "q", a, "T", 2.0)
bloc.add_entry(a, "q", b, "T", -1.0)
bloc.add_entry(b, "q", a, "T", -1.0)
bloc.add_entry(b, "q", b, "T", 2.0)
k.finalize()  # required before any solver use

assert k.n_rows() == 2
assert k.n_cols() == 2
assert k.symmetric is True

# What it weighs, estimated: the four entries of the block, plus the assembled
# CSR. The factorization, larger than both, cannot be counted (faer keeps the
# size of its factors private).
assert bloc.memory_bytes() == 4 * 24
assert k.memory_bytes() > bloc.memory_bytes()

Ce que pèse une matrice

memory_bytes(), sur un bloc comme sur la matrice, estime les octets de tas occupés. Pour un bloc, ce sont ses entrées stockées : un indice de ligne, un indice de colonne et une valeur, soit 24 octets chacune — un bloc calculé n’en stocke aucune et répond 0. Pour la matrice, c’est la CSR assemblée (un indice de colonne et une valeur par terme non nul, un décalage par ligne, une clé de DDL par ligne et par colonne) plus ce que gardent ses blocs.

L’estimation apparaît dans l’affichage : Matrix: 2 sous-matrice(s), 3 row(s) × 3 col(s), symmetric, ~1.2 kB.

Ce qu’elle ne compte pas, et qui est pourtant le plus gros : la factorisation. Ses facteurs pèsent vingt à soixante-cinq fois la matrice, mais faer garde leur taille privée. Pour la mémoire réellement consommée par un solve, voir Calculs plus gros que la RAM.

Sérialisation

Matrix implémente Portable via serde (comme tous les objets pyrucast). Les triplets COO, la table de noms et les DOFs voyagent dans le format binaire portable Linux ↔ Windows. La CSR assemblée et la factorisation, elles, ne sont pas écrites : elles se reconstruisent (voir Sauvegarde et relecture).

Limitations actuelles

• Cache de motif non invalidé par les mutations profondes : le motif creux mémoïsé sur le Model est invalidé à l’ajout d’un sous-modèle (add_sub), mais pas si le maillage / l’espace EF sous-jacent change en place (remaillage) — reconstruire le modèle dans ce cas. Le chemin de composition m.assemble(), lui, reconstruit toujours le motif depuis les blocs.
• Pas de produit matrice-matrice : à venir avec les premiers besoins concrets (préconditionneurs, formulations couplées).
• La somme n’assemble pas de manière opportuniste : a + b rend une matrice non assemblée même quand les deux opérandes le sont. Fusionner leurs CSR — ce qui éviterait de relancer les noyaux élémentaires dans une boucle en temps à pas variable — est possible sans changer la sémantique (l’ordre des DDL d’une concaténation est exactement celui de a suivi des DDL que seule b apporte), mais demande une addition creuse complète : retable des variables, remappage et retri des colonnes de b, fusion ligne à ligne. À faire quand un intégrateur en temps le justifiera.
• La symétrie déclarée n’est pas vérifiée numériquement à l’assemblage, et ne doit pas l’être : c’est une déclaration du modèle, pas une mesure. Des tests unitaires confrontent la déclaration à la CSR réellement assemblée ; le calcul, lui, fait confiance.
• Une symétrie découpée en tranches de lignes n’est pas exprimable : deux blocs rectangulaires qui sont chacun une tranche de lignes d’une matrice symétrique ne peuvent rien déclarer — ni Full, qui suppose un bloc carré, ni Half, les deux n’étant pas transposés l’un de l’autre. Le tableau assemblé est symétrique, le drapeau répond false, et le calcul prend le chemin général. C’est le mode de défaillance voulu : on oublie une symétrie, on n’en invente pas.

Évolution (Evolution / SubEvolution)

Une évolution associe une série de valeurs à une variable (souvent le temps, mais pas nécessairement) et interpole linéairement entre les échantillons tabulés. C’est l’analogue de l’EVOLUTIO de Cast3M, généralisé : l’abscisse n’est pas forcément le temps et la valeur peut être un champ entier, pas seulement un réel.

Elle suit la même grammaire d’agrégat que tous les conteneurs de pyrucast (cf. Agrégat) :

• SubEvolution — une courbe tabulée : une liste d’abscisses (triées, strictement croissantes) et la liste des valeurs en regard. Une valeur est un scalaire, un SubNodeField (champ aux nœuds) ou un SubElementField (champ aux points de Gauss) — toutes du même type, et pour les champs sur le même support. Son interpolation en x rend une valeur.
• Evolution — l’agrégat : une liste de SubEvolution, une par zone, exactement comme un NodeField agrège des SubNodeField. Son interpolation en x interpole chaque courbe, puis regroupe les sous-champs résultants en un NodeField / ElementField. Pour des scalaires, elle rend une liste de flottants (il n’existe pas d’agrégat de flottant).

   Evolution (agrégat)
   ├── SubEvolution zone 0 ── abscisses [t₀, t₁, …] × valeurs [v₀, v₁, …]
   ├── SubEvolution zone 1 ── …
   └── …

Interpolation linéaire

Entre deux échantillons encadrants x_lo ≤ x ≤ x_hi, le résultat est le mélange v_lo·(1−t) + v_hi·t avec t = (x − x_lo) / (x_hi − x_lo). Pour les champs, le mélange réutilise l’arithmétique de champs (map_all + merge_components précédé de check_same_components, cf. Champ) : aucune logique numérique n’est dupliquée, et la compatibilité des supports/composantes des deux champs encadrants est vérifiée à ce moment-là. Une abscisse tombant exactement sur un échantillon rend la valeur telle quelle.

Types d’abscisse et d’ordonnée

Une évolution peut porter le type physique de ses axes :

• abscissa_type — le type de l’abscisse (p. ex. "T", "time"). Valable pour toutes les évolutions. Il sert à étiqueter les tracés (axe X d’une courbe, slider d’un champ) et, lorsqu’on interpole un champ, à choisir la composante du champ à lire (voir ci-dessous).
• ordinate_type — le type de la valeur, pour les évolutions scalaires uniquement (p. ex. "young"). Il étiquette l’axe Y et nomme la composante produite quand on interpole un champ. Le donner sur une évolution de champs est une erreur (un champ a déjà ses propres composantes).

se = pc.SubEvolution(
    [(0.0, 0.0), (100.0, 210e9)], abscissa_type="T", ordinate_type="young"
)

Interpoler un champ (courbe de transfert)

Une évolution scalaire à une seule courbe s’utilise comme une fonction de transfert y = f(x) : au lieu d’un scalaire, on lui passe un champ et elle rend un autre champ de même support, où chaque nœud / point de Gauss est l’interpolation de la valeur d’entrée sur la courbe.

• La composante lue dans le champ d’entrée est celle nommée comme l’abscissa_type — la correspondance de type est vérifiée : si le champ n’a pas de composante de ce nom, c’est une erreur.
• Le champ de sortie a une seule composante, nommée d’après l’ordinate_type (à défaut "value").
• La politique hors-plage s’applique valeur par valeur, comme pour un scalaire.

# Loi matériau E(T) : module d'Young fonction de la température.
loi = pc.Evolution(
    [(0.0, 0.0), (100.0, 210e9)], abscissa_type="T", ordinate_type="young"
)
young = loi.interpolate(temperature)  # temperature : NodeField de composante "T"
# young : NodeField de composante "young"

Côté Evolution (agrégat), l’appel exige une seule courbe scalaire (sans quoi le choix de la courbe serait ambigu) ; une SubEvolution s’interpole directement.

Politique hors plage

Chaque évolution porte une politique appliquée quand l’abscisse demandée sort de l’intervalle tabulé [x_min, x_max] :

La politique stockée peut être surchargée à l’appel : evol.interpolate(x, out_of_range="clamp").

Construction

Deux voies, le constructeur haut niveau n’étant que du sucre au-dessus du primitif bas niveau (motif model.heat_conduction(fes) / SubModel + |) :

• temps-major (haut niveau) — Evolution([(t0, champ0), (t1, champ1), …]) avec un NodeField / ElementField / flottant complet par pas ; les champs entiers sont transposés en une courbe par zone (zones appariées entre pas par leur support, identique d’un pas à l’autre) ;
• zone-major (bas niveau) — construire chaque SubEvolution depuis sa liste (abscisse, sous-champ), puis agréger avec |.

L’union | et le slicing réinitialisent la politique hors-plage de l’agrégat à "error".

API Rust

#[test]
fn une_evolution_interpole_scalaires_et_champs() -> Result<()> {
    // Courbe scalaire X→Y : 0→10, 1→20.
    let se = SubEvolution::new(
        vec![(0.0, SubValue::Scalar(10.0)), (1.0, SubValue::Scalar(20.0))],
        OutOfRange::Error,
    )?;
    match se.interpolate(0.5, None)? {
        SubValue::Scalar(v) => assert_eq!(v, 15.0),
        _ => unreachable!(),
    }
    // Hors plage : Error (défaut) lève ; surcharge Clamp → extrémité.
    assert!(se.interpolate(2.0, None).is_err());

    // Agrégat scalaire → liste de flottants.
    let e = Evolution::from_scalars(vec![(0.0, 10.0), (1.0, 20.0)], OutOfRange::Error)?;
    match e.interpolate(0.5, None)? {
        Interpolated::Scalars(v) => assert_eq!(v, vec![15.0]),
        _ => unreachable!(),
    }
    Ok(())
}

API Python

import pyrucast as pc

# Scalar curve (one SubEvolution).
se = pc.SubEvolution([(0.0, 10.0), (1.0, 20.0)])
print(se.interpolate(0.5))  # 15.0
print(se.interpolate(2.0, out_of_range="clamp"))  # 20.0 (otherwise: an error)

# Agrégat scalaire → liste de flottants.
e = pc.Evolution([(0.0, 10.0), (1.0, 20.0)])
print(e.interpolate(0.5))  # [15.0]

# Low level: composing per-zone curves with `|`.
agg = pc.SubEvolution([(0.0, 1.0), (1.0, 2.0)]) | pc.SubEvolution(
    [(0.0, 3.0), (1.0, 4.0)]
)
print(agg.interpolate(0.5))  # [1.5, 3.5]

# High level, time-major: one whole NodeField per step → interpolated NodeField.
ev = pc.Evolution([(0.0, champ_t0), (2.0, champ_t1)])
champ = ev.interpolate(1.0)  # NodeField halfway

# Courbe de transfert : passer un champ → champ (loi matériau E(T)).
loi = pc.Evolution(
    [(0.0, 0.0), (100.0, 210e9)], abscissa_type="T", ordinate_type="young"
)
young = loi.interpolate(temperature)  # composante "T" lue → composante "young"

Tracé

evolution.plot(...) visualise l’évolution : courbe X-Y pour des scalaires, champ + slider de valeur tabulée pour des champs. Voir Visualisation › Tracé d’une évolution.

À défaut de x_label / y_label explicites, les étiquettes reprennent l’abscissa_type (axe X d’une courbe, slider d’un champ) et l’ordinate_type (axe Y d’une courbe).

e = pc.Evolution([(0.0, 10.0), (1.0, 20.0), (2.0, 5.0)])
e.plot(save="courbe.svg", x_label="temps", y_label="T")  # courbe scalaire
ev.plot(save="frame.png", frame=1)  # tabulated field (one step)

Place dans le modèle

SubValue est un enum de stockage inline (scalaire / SubNodeField / SubElementField), comme SubModel l’est pour les physiques. SubEvolution s’adresse par un Handle<SubEvolution> et sérialise ses valeurs en ligne via le trait Portable ; les courbes sont donc portables comme tout autre objet. L’homogénéité du type de valeur est garantie à la construction (au sein d’une courbe) et par check_push (entre zones d’un même agrégat).

Détails des physiques

Chaque physique est une variante de SubModel : elle déclare ses variables (primales / duales), son matériau, et sait assembler sa rigidité K (et, le cas échéant, intégrer son comportement COMP). Le Model orchestre, mais ne porte aucune logique physique — voir Modèle physique pour la mécanique générique (DOFs, assemblage, chargements, solveur).

Chaque page suit le même plan standard, dans cet ordre :

1. Introduction — nature du modèle, éléments et degrés de liberté.
2. Équations continues résolues — forme forte (et faible) du problème.
3. Forme discrétisée — opérateurs discrets (B, N) et expressions des matrices élémentaires (K, masse/capacité, tangente…).
4. Variables et matériau — noms primal/dual, composantes matériau, comportement (COMP) et, le cas échéant, état interne.
5. Mise en donnée (Rust, testé) — exemple Rust exécuté via {{#include}}.
6. Exemple Python — l’équivalent haut niveau.
7. Compléments (facultatif) — extensions propres au modèle (masse & rigidité géométrique, thermomécanique, convection, régime transitoire…).

Certaines physiques omettent une rubrique sans objet (p. ex. pas d’exemple Rust dédié pour un modèle piloté par le comportement) ; l’ordre reste identique.

• Conduction thermique — -∇·(k∇T) = 0, l’exemple canonique, la conduction orientée (tenseur K) et la convection de surface (Robin / film, q·n = h(T − T_ext)), plus le rayonnement à l’infini (q·n = σε(T⁴ − T_∞⁴), non linéaire, à tangente cohérente).
• Diffusion (loi de Fick) — ∇·(D∇c) = 0, concentration c_<espèce> et flux de matière j_<espèce>, l’espèce étant nommée à la construction ; même opérateur que la conduction, nature distincte. Avec le transfert d’interface j·n = h(c₁ − c₂), qui laisse le champ sauter entre deux corps.
• Échanges (frontière et interface) — la loi h(a − b), partagée par l’échange avec une ambiante et le transfert entre deux maillages : une seule loi, dont la seule différence est que l’autre côté soit une donnée ou une inconnue. Générique en composantes, donc aussi bien un film thermique qu’une fondation élastique.
• Mécanique — barre, élasticité linéaire, poutres et portiques :
  • Barre / treillis
  • Élasticité linéaire
  • Élasticité orthotrope et anisotrope
  • Plasticité parfaite (von Mises)
  • Lois d’écoulement plastique — écrouissage isotrope, Drucker-Prager, Ottosen
  • Fluage et viscoplasticité — Norton, Lemaitre, Blackburn, Chaboche et sa variante endommageable
  • Endommagement de Mazars
  • Lois d’endommagement — Damage TC, SiC/SiC orthotrope, plasticité poreuse de Gurson
  • Poutre d’Euler-Bernoulli — 1-D, plan, spatial ; exacte aux nœuds
  • Poutre de Timoshenko
  • Portique 2D
  • Cadre 3D
  • Coques — Reissner-Mindlin ou Kirchhoff discret, six DDL par nœud
• Contraintes — conditions limites imposées par multiplicateurs de Lagrange :
  • Dirichlet
  • Multi-points (MPC)
  • Baignage (embedded)
  • Contact (nœud-surface)

Ce regroupement est la nature physique (Physics) que chaque variante déclare : Thermal (conduction), Radiation (rayonnement, en plus de Thermal), Diffusion (Fick), Mechanical (barre → cadre 3D) et Constraint (les contraintes de Lagrange). On sélectionne les sous-modèles d’une nature avec model.filter(Physics::Mechanical) (et les blocs d’une matrice avec k.filter(...)) — voir Nature physique et filtrage.

Ce que chacune sait produire

La liste ci-dessus dit où lire chaque physique. Les tableaux qui suivent disent ce qu’elle sait produire : les genres de matrice qu’elle déclare (MatrixKind), la voie par laquelle elle obtient sa tangente, la nature de son intégration de comportement, et la particularité de calcul qui la distingue.

Ils se lisent en colonnes. Une case remplie dit que la physique déclare ce terme et que l’assembleur l’appelle ; un tiret dit qu’elle n’y contribue rien — et ce n’est pas un manque. matrix.mass(...) sur un modèle contenant un échange de bord voit simplement qu’il n’a pas de masse.

Les tags donnés sous chaque nom sont les chaînes exactes que prennent les opérateurs (from_name) : c’est ce qu’on écrit, pas une paraphrase.

Thermique — primale T, duale q, filter("thermal")

Diffusion — primale c_<espèce>, duale j_<espèce>, filter("diffusion")

Mécanique — milieux continus, primales u, duales f, filter("mechanical")

Mécanique — structurel, efforts de section, filter("mechanical")

Contraintes — multiplicateurs de Lagrange, filter("constraint")

« Par perturbation » veut dire différences centrées

Pour les huit lois qui n’ont pas de tangente en forme fermée, D_alg est obtenu en perturbant la déformation et en relançant le retour, composante par composante :

\[ D_{ij} \simeq \frac{\sigma_i(\varepsilon + h\,e_j) - \sigma_i(\varepsilon - h\,e_j)}{2h}, \qquad h = 10^{-6}\,\lVert \varepsilon \rVert_\infty . \]

Six composantes de Voigt, deux évaluations chacune : douze appels au retour par point de Gauss. Le pas doit rester bien au-dessus du bruit du retour (celui d’Ottosen ou de Gurson converge à une tolérance, pas exactement) et bien en dessous de l’échelle de courbure de la surface ; 1e-6·‖ε‖ tient confortablement entre les deux. Les colonnes de cisaillement sont divisées par deux en sortie, ce qui transforme ∂σ/∂ε_ij en ∂σ/∂γ_ij — la convention ingénieur du reste du dépôt.

Le partage n’est pas une question de difficulté mais de vérifiabilité : seule la forme fermée de von Mises a été confrontée à une différence finie et validée. La dérivation analytique de Drucker-Prager, écrite d’abord, était fausse de 24 % — plausible, et fausse ; seul l’oracle numérique l’a dit. Une tangente obtenue par perturbation ne peut pas être mal dérivée, coûte douze évaluations d’une mise à jour bon marché, et laisse Newton converger. C’est un bon échange.

Le rayonnement, lui, a bien une tangente analytique (4σεT³ ∫NᵢNⱼ) : sa non-linéarité est une puissance scalaire d’une seule variable, pas une carte de projection.

Trois choses que ces tableaux ne disent pas

Les forces internes ne suivent pas toujours Bᵀσ. Le défaut est le noyau de la mécanique des milieux continus, f_i = ∫ ∂N_i/∂x · σ. Une physique dont la duale n’est pas un vecteur déplacement le redéfinit : la thermique et la diffusion appliquent Bᵀ à un flux scalaire, tandis que la convection, le rayonnement et le transfert d’interface pondèrent par N et non par Bᵀ — leur intégrande est une densité surfacique, pas une grandeur conjuguée d’un gradient.

Les structurels le redéfinissent aussi, pour une autre raison : leur B existe, mais ce n’est pas le gradient symétrique. Barre, poutres et coques intègrent chacun le transposé du leur — truss::internal_force_element, beam::internal_force_into, shell::b_into — et c’est le même B que leur rigidité, ce que mesure tests/internal_forces.rs : la loi structurelle étant linéaire, ∫ Bᵀσ vaut K·u exactement, à l’arrondi près.

Deux d’entre eux ont fallu payer un prix pour cela. Une poutre de Timoshenko porte Φ dans son état, son interpolation en dépendant. Une coque a gagné un effort conjugué, M_drill : le vrillage travaillait dans K sans figurer nulle part ailleurs.

La tangente est symétrisée. D_alg est réduit au triangle supérieur puis relu en miroir, un format qui ne peut structurellement pas porter une matrice non symétrique. Or un écoulement non associé — Drucker-Prager, dont la dilatance diffère du frottement — en a une. Elle est donc symétrisée, ce qui coûte à Newton son taux quadratique sur cette loi, et rien d’autre. Voir Lois d’écoulement plastique.

État absent ≠ état nul. Au-delà de ε_p et p, une loi déclare l’état qu’elle veut par internal_names() : la déformation primaire de Blackburn, la contrainte de rappel de Chaboche (un tenseur complet), la porosité de Gurson, les deux histoires de Damage TC. Le premier pas passe un vecteur vide — et non un vecteur de zéros — pour qu’une loi démarrant d’une constante matériau puisse faire la différence. Sans cela, un métal poreux démarrerait dense et ne s’endommagerait jamais.

Pour ajouter une physique, voir Ajouter une physique.

Conduction thermique

Cette page décrit la physique de conduction thermique (HeatConduction) et la convection de surface associée. Elle suit le plan standard des physiques puis la déroule sur un exemple complet — une ligne chauffée par une source à une extrémité et maintenue à température fixe à l’autre — comparé à la solution analytique.

Pour la mécanique générique du Model (orchestration, DOFs, assemblage), voir Modèle physique. Ici on se concentre sur le cas thermique.

Équations continues résolues

En régime stationnaire, la forme forte est

\[ -\nabla\cdot\big(k\,\nabla T\big) = 0, \]

et en régime transitoire, l’équation de la chaleur porte un terme de stockage :

\[ \rho\,c_p\,\frac{\partial T}{\partial t} - \nabla\cdot(k\,\nabla T) = Q. \]

La forme variationnelle de Galerkine (multiplication par une température virtuelle, intégration par parties) fait apparaître la rigidité (conduction) et, en transitoire, la capacité (stockage) — leurs formes discrètes ci-dessous.

Forme discrétisée

La conductivité donne, cellule par cellule, la matrice de rigidité :

\[ K_{ij} = \int_K k(x)\,\nabla N_i\cdot\nabla N_j\,dx \quad\approx\quad \sum_g k(\xi_g)\,(\nabla N_i\cdot\nabla N_j)\big|_g\,|J|_g\,w_g \]

(implémentée dans src/models/heat_conduction.rs). En notant \( B = [\nabla N_1, \dots, \nabla N_n] \) la matrice des gradients de forme (taille \( d\times n \)), on a aussi \( K = \int_\Omega k\, B^\top B\, d\Omega \). Le bloc local est écrit aux positions row = (NodeId_i, "q"), col = (NodeId_j, "T"). Pour un SEG2 de longueur \(L\) et \(k\) uniforme on retrouve la matrice analytique \((k/L)\,[[1,-1],[-1,1]]\).

En transitoire, le terme de stockage discrétise en une matrice de capacité (l’analogue thermique de la matrice de masse, Cast3M CAPA) :

\[ C_{ij} = \int_\Omega \rho\,c_p\,N_i\,N_j\,d\Omega \;\approx\; \sum_g \rho\,c_p\,N_i(\xi_g)\,N_j(\xi_g)\,|J|_g\,w_g, \]

assemblée par assemble.mass (matériau rho, cp) et concentrable en diagonale par lump. Le système semi-discret est \( C\,\dot T + K\,T = F \) ; l’intégration en temps (θ-schéma, Euler implicite (C/\Delta t + K)) se pilote dans la couche Python.

Variables et matériau

La conductivité peut être orientée — voir Conduction orthotrope et anisotrope plus bas ; "k" est alors remplacée par les constantes de la symétrie choisie.

Mise en donnée (Rust, testé)

Le pipeline est toujours le même :

1. Coords — l’espace des nœuds (dimension géométrique).
2. Mesh — les éléments (ici des SEG2 alignés sur \([0,1]\)).
3. FiniteElementSpace — l’interpolation (lagrange1).
4. Matériau — un ElementField portant la composante "k", fabriqué commodément par element_field::material_field(&model, &[("k", …)]) (les sous-modèles sans matériau, comme Dirichlet, sont ignorés).
5. Model — model::heat_conduction(&fes), composé par | (union) avec les conditions limites.
6. Conditions limites :
   • Dirichlet (T imposée) : un sous-modèle model::dirichlet qui impose la valeur via multiplicateurs de Lagrange. L’utilisateur fournit le maillage des nœuds imposés et le maillage support des multiplicateurs — typiquement fabriqué depuis le premier avec le mesher générique barycenter (nœuds neufs colocalisés). La valeur imposée \(u_d\) s’écrit dans le chargement au slot imposed_T du nœud-multiplicateur (cf. Modèle physique).
   • Neumann / source : une charge ponctuelle est une valeur du chargement sur la composante duale "q" au nœud concerné ; un flux réparti sur un bord (ou un volume) se transforme en charges nodales cohérentes par l’opérateur flux (analogue de FLUX/PRES de Cast3M).
7. Assemblage + résolution — matrix::stiffness puis le solveur solver::lu::solve (LU creuse directe, voir Modèle physique).

Exemple : ligne chauffée

Problème. Sur \([0,1]\), une source de chaleur (flux de Neumann \(Q\)) est appliquée en \(x=0\), et la température est imposée à \(T=20\) en \(x=1\).

Solution analytique. Sans génération volumique, \(T’’=0\) : le profil est linéaire. En notant \(Q\) le flux injecté et \(k\) la conductivité,

\[ u(x) = 20 + \frac{Q}{k}\,(1 - x). \]

De plus, le multiplicateur de Lagrange au nœud imposé (la réaction qui maintient \(T=20\)) vaut exactement \(Q\) : tout le flux injecté en \(x=0\) ressort en \(x=1\) — un bilan d’énergie discret.

Code. L’exemple ci-dessous est le test d’intégration tests/thermal_line.rs : il est compilé et exécuté à chaque cargo test, donc garanti à jour avec l’API.

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::node_field::{NodeField, SubNodeField};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

#[test]
fn thermal_line_recovers_analytical_solution() -> Result<()> {
    // ── Problem data ───────────────────────────────────────────────────────
    const K: f64 = 1.0; // conductivité
    const Q: f64 = 10.0; // source de chaleur (flux de Neumann) en x = 0
    const T_IMPOSED: f64 = 20.0; // température imposée en x = 1
    const N_ELEMS: usize = 4;
    let h = 1.0 / N_ELEMS as f64;

    // ── Mesh: a line of SEG2 on [0, 1] ─────────────────────────────────────
    let coords = Handle::new(Coords::new(1)?);
    let nodes: Vec<Node> = (0..=N_ELEMS)
        .map(|i| Node::create_in(coords.clone(), &[i as f64 * h]))
        .collect::<Result<_>>()?;
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for i in 0..N_ELEMS {
        mesh.add_cell(&[nodes[i].id(), nodes[i + 1].id()])?;
    }
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Modèle : conduction + Dirichlet T = 20 en x = 1 ────────────────────
    // The multipliers' support is built from the imposed node by the
    // `barycenter` mesher (a fresh co-located node). The model creates nothing.
    let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(std::slice::from_ref(
        nodes.last().unwrap(),
    ))?);
    let multiplier = mesh::barycenter(&imposed)?;
    let mult = multiplier.node(0, 0, 0)?.id();

    let conduction = model::heat_conduction(&fes)?;
    let dirichlet = model::dirichlet(&conduction, "T", &imposed, &multiplier, Default::default())?;
    let model = conduction.union(&dirichlet)?;

    // ── Material: uniform k (Dirichlet is skipped automatically) ───────────
    let materials = pyrucast::ops::element_field::material_field(&model, &[("k", K)])?;

    // ── Chargement : source Q en x = 0 (composante duale "q"), valeur imposée
    //    T = 20 at the multiplier node ("imposed_T" slot) ──────────────────
    let node0 = nodes[0].id();
    let mut load_sm = SubMesh::new(coords.clone(), ElementType::POI1);
    load_sm.add_cell(&[node0])?;
    load_sm.add_cell(&[mult])?;
    let load_sm = Handle::new(load_sm);
    let mut rhs = SubNodeField::from_poi1(&load_sm, vec!["imposed_T".into(), "q".into()])?;
    rhs.set_value(node0, "q", Q)?;
    rhs.set_value(mult, "imposed_T", T_IMPOSED)?;
    let rhs = NodeField::from_sub(rhs);

    // ── Assemblage + résolution ────────────────────────────────────────────
    let stiffness = pyrucast::ops::matrix::stiffness(&model, &materials)?;
    let solution = solve(&stiffness, &rhs)?;

    // ── Compared with the analytical solution u(x) = 20 + (Q/k)(1 − x) ─────
    let tol = 1e-10;
    for (i, node) in nodes.iter().enumerate() {
        let x = i as f64 * h;
        let expected = T_IMPOSED + (Q / K) * (1.0 - x);
        let got = solution.value(node.id(), "T")?;
        assert!(
            (got - expected).abs() < tol,
            "T(x={x}) : obtenu {got}, attendu {expected}"
        );
    }
    // La réaction (multiplicateur de Lagrange) équilibre le flux injecté : λ = Q.
    let reaction = solution.value(mult, "lambda_T")?;
    assert!(
        (reaction - Q).abs() < tol,
        "réaction λ : obtenue {reaction}, attendue {Q}"
    );

    Ok(())
}

Exemple Python

La version Python équivalente et documentée est dans le dépôt : examples/thermal_line_1d.py (lancer avec python examples/thermal_line_1d.py après maturin develop). Les compléments 2-D ci-dessous ont eux aussi leur variante Python (thermal_square_2d.py, thermal_convection_2d.py).

Compléments

Exemple : un carré

La généralisation 2-D du cas précédent : un carré \([0,1]^2\) (grille structurée de QUA4), chauffé par une source répartie sur le bord gauche (\(x=0\)) et maintenu à \(T=20\) sur le bord droit (\(x=1\)). Les bords haut et bas ne portent aucune condition : c’est la condition naturelle (flux nul, bord isolé).

Comme les bords latéraux sont isolés, le champ ne dépend pas de \(y\) : le carré redonne le profil de la ligne,

\[ u(x) = 20 + \frac{Q}{k}\,(1 - x), \]

et la réaction totale (somme des multiplicateurs sur le bord imposé) vaut le flux injecté \(Q\).

Mise en donnée d’un flux réparti. Une source répartie se transforme en charges nodales cohérentes \(f_i = \int_\Gamma \varphi\,N_i\,d\Gamma\) par l’opérateur flux — l’analogue de FLUX/PRES de Cast3M. On lui donne le bord (ici un maillage SEG2, intégré comme une ligne : la mesure vient du Jacobien manifold) et la densité de flux (une constante, ou un champ par éléments) ; il renvoie un NodeField sur la composante duale "q", prêt à composer (|) avec le reste du chargement. Sous le capot, pour un flux uniforme sur des éléments linéaires, un nœud intérieur du bord reçoit \(Q\,h\) et un coin \(Q\,h/2\) (somme \(Q\)) — mais on n’a plus à le calculer à la main.

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::node_field::{NodeField, SubNodeField};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

#[test]
fn thermal_square_recovers_analytical_solution() -> Result<()> {
    // ── Problem data ───────────────────────────────────────────────────────
    const K: f64 = 1.0; // conductivité
    const Q: f64 = 10.0; // flux de chaleur TOTAL injecté sur le bord gauche
    const T_IMPOSED: f64 = 20.0; // température imposée sur le bord droit
    const N: usize = 4; // N×N éléments QUA4
    let h = 1.0 / N as f64;

    // ── Mesh: a structured (N+1)×(N+1) grid of QUA4 on [0,1]² ──────────────
    let coords = Handle::new(Coords::new(2)?);
    let idx = |i: usize, j: usize| j * (N + 1) + i; // nœud colonne i, ligne j
    let mut grid: Vec<Node> = Vec::with_capacity((N + 1) * (N + 1));
    for j in 0..=N {
        for i in 0..=N {
            grid.push(Node::create_in(
                coords.clone(),
                &[i as f64 * h, j as f64 * h],
            )?);
        }
    }
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::QUA4));
    for j in 0..N {
        for i in 0..N {
            mesh.add_cell(&[
                grid[idx(i, j)].id(),
                grid[idx(i + 1, j)].id(),
                grid[idx(i + 1, j + 1)].id(),
                grid[idx(i, j + 1)].id(),
            ])?;
        }
    }
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Dirichlet T = 20 on the right edge (x = 1) ──────────────────────────
    let right_nodes: Vec<Node> = (0..=N).map(|j| grid[idx(N, j)].clone()).collect();
    let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(&right_nodes)?);
    let multiplier = mesh::barycenter(&imposed)?;
    let mults: Vec<Node> = (0..=N)
        .map(|j| multiplier.node(0, j, 0))
        .collect::<Result<_>>()?;

    let conduction = model::heat_conduction(&fes)?;
    let dirichlet = model::dirichlet(&conduction, "T", &imposed, &multiplier, Default::default())?;
    let model = conduction.union(&dirichlet)?;
    // ── Chargement ─────────────────────────────────────────────────────────
    // Source: uniform flux (density Q) on the left edge, turned into
    // consistent nodal loads by the `flux` operator (Cast3m FLUX) — no more
    // Q·h / Q·h/2 distribution by hand. The edge is a SEG2 mesh built on
    // the grid's nodes; it integrates as a line.
    let mut left_edge = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for j in 0..N {
        left_edge.add_cell(&[grid[idx(0, j)].id(), grid[idx(0, j + 1)].id()])?;
    }
    let left_fes = FiniteElementSpace::lagrange1(&left_edge)?;
    let model = model.union(&model::flux(&left_fes, &model, "q".into())?)?;

    let materials =
        pyrucast::ops::element_field::material_field(&model, &[("k", K), ("phi_q", Q)])?;
    let source = pyrucast::ops::node_field::external_forces(&model, &materials)?;

    // Imposed value T = 20 at the multiplier nodes' "imposed_T" slot.
    let mut imposed_sm = SubMesh::new(coords.clone(), ElementType::POI1);
    for m in &mults {
        imposed_sm.add_cell(&[m.id()])?;
    }
    let imposed_sm = Handle::new(imposed_sm);
    let mut imposed_load = SubNodeField::from_poi1(&imposed_sm, vec!["imposed_T".into()])?;
    for m in &mults {
        imposed_load.set_value(m.id(), "imposed_T", T_IMPOSED)?;
    }

    // Loading = the edge's flux + the imposed values (union of the zones).
    let rhs = source.union(&NodeField::from_sub(imposed_load))?;

    // ── Assemblage + résolution ────────────────────────────────────────────
    let stiffness = pyrucast::ops::matrix::stiffness(&model, &materials)?;
    let solution = solve(&stiffness, &rhs)?;

    // ── Compared with the analytical u(x) = 20 + (Q/k)(1 − x), ∀ y ─────────
    let tol = 1e-9;
    for j in 0..=N {
        for i in 0..=N {
            let x = i as f64 * h;
            let expected = T_IMPOSED + (Q / K) * (1.0 - x);
            let got = solution.value(grid[idx(i, j)].id(), "T")?;
            assert!(
                (got - expected).abs() < tol,
                "T(x={x}, y={}) : obtenu {got}, attendu {expected}",
                j as f64 * h
            );
        }
    }
    // The total reaction on the imposed edge balances the injected flux: Σλ = Q.
    let total_reaction: f64 = mults
        .iter()
        .map(|m| solution.value(m.id(), "lambda_T"))
        .sum::<Result<f64>>()?;
    assert!(
        (total_reaction - Q).abs() < tol,
        "réaction totale : obtenue {total_reaction}, attendue {Q}"
    );

    Ok(())
}

Version Python : examples/thermal_square_2d.py (lancer avec python examples/thermal_square_2d.py après maturin develop).

Conduction orthotrope et anisotrope

Un matériau feuilleté, fibré ou laminé ne conduit pas la chaleur de la même façon dans toutes les directions. La conductivité devient alors un tenseur K, et la rigidité

\[ K_{ij} = \int_\Omega \nabla N_i^{\mathsf T}\, \mathbf{K}\, \nabla N_j \, d\Omega \]

dont le cas isotrope K = k·I redonne le produit scalaire habituel.

C’est le même axe de symétrie matériau qu’en mécanique (chapitre Élasticité orthotrope), avec un tenseur d’ordre 2 au lieu de 4 :

Le repère est donné par des vecteurs — V1X, V1Y en 2-D, V1X…V1Z, V2X…V2Z en 3-D — comme MATE 'DIRECTION' V1 V2 de Cast3M. Ils sont orthonormalisés en interne.

model = pyrucast.model.heat_conduction(fes, symmetry="orthotropic")
materials = pyrucast.element_field.material_field(
    model,
    [("k_1", 12.0), ("k_2", 3.0), ("k_3", 12.0), ("V1X", cos_a), ("V1Y", sin_a)],
)

La conductivité isotrope reste lue au point de Gauss, donc variable à l’intérieur d’une maille ; les constantes orientées sont lues par maille, comme les modules mécaniques.

L’exemple Rust est un test de patch, qui est ce qu’appelle une conductivité orientée : un champ de température linéaire est harmonique pour n’importe quel tenseur constant, donc l’imposer au bord doit le reproduire à l’intérieur quelle que soit K. Le test ne s’arrête pas là — il relit le flux produit et le compare à K·∇T calculé à la main, ce qui est le seul moyen de prendre la rotation en défaut :

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::element_field::ElementField;
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::model::Model;
use pyrucast::containers::node_field::{NodeField, SubNodeField};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::symmetry::MaterialSymmetry;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

/// `N×N` QUA4 grid on the unit square.
const N: usize = 3;

#[test]
fn orthotropic_conduction_passes_the_linear_patch_test() -> Result<()> {
    const K1: f64 = 12.0; // along the first material axis
    const K2: f64 = 3.0; // transverse
    let theta = 30.0_f64.to_radians();
    let (c, s) = (theta.cos(), theta.sin());

    let (grid, fes, _coords) = unit_square()?;
    let (model, multipliers) = patch_model(&grid, &fes, MaterialSymmetry::Orthotropic)?;
    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[
            ("k_1", K1),
            ("k_2", K2),
            ("k_3", K1),
            ("V1X", c),
            ("V1Y", s),
        ],
    )?;

    let solution = solve_patch(&model, &materials, &grid, &multipliers)?;

    // The linear field must be reproduced exactly, tensor or no tensor.
    let h = 1.0 / N as f64;
    let tol = 1e-9;
    for j in 0..=N {
        for i in 0..=N {
            let x = i as f64 * h;
            let got = solution.value(grid[j * (N + 1) + i].id(), "T")?;
            assert!((got - x).abs() < tol, "T({x}) = {got}");
        }
    }

    // …and the flux must be the first column of the **rotated** tensor.
    let expect_xx = K1 * c * c + K2 * s * s;
    let expect_yx = (K1 - K2) * c * s;
    let (fx, fy) = uniform_flux(&model, &solution, &fes, &materials, &grid)?;
    assert!(
        (fx - expect_xx).abs() < 1e-9,
        "flux_x = {fx}, expected {expect_xx}"
    );
    assert!(
        (fy - expect_yx).abs() < 1e-9,
        "flux_y = {fy}, expected {expect_yx}"
    );
    Ok(())
}

Avec ∇T = (1, 0), le flux est la première colonne de K : K_xx = k₁cos²θ + k₂sin²θ et K_yx = (k₁ − k₂)·cosθ·sinθ. Le terme extra-diagonal n’est non nul que si le matériau est à la fois anisotrope et désaligné — précisément le cas qu’une rotation fausse manquerait.

Rayonnement à l’infini (Stefan-Boltzmann)

Une surface qui échange avec un environnement lointain à \(T_\infty\) rayonne

\[ q\cdot n = \sigma\,\varepsilon\,\big(T^4 - T_\infty^4\big) \]

où \(\sigma\) est la constante de Stefan-Boltzmann et \(\varepsilon\) l’émissivité. Primale "T", duale "q" — les mêmes degrés de liberté que la conduction, donc un bord rayonnant se couple directement dans sa rigidité, comme la convection. Et comme elle, il n’a besoin d’aucune normale : la direction est déjà consommée en écrivant q·n, il ne reste sous l’intégrale qu’un scalaire et la mesure de surface.

Ce qui change par rapport à la convection : c’est non linéaire

La loi de Newton est linéaire en T, si bien que la convection ne contribue qu’une matrice de film constante. T⁴ ne l’est pas, d’où trois termes :

Linéariser la rigidité autour de \(T_\infty\) plutôt qu’autour de l’état courant est ce qui la laisse être une matrice constante : c’est l’opérateur dont on part pour une boucle de Newton, et à lui seul une itération de Picard tout à fait utilisable. La tangente porte la vraie non-linéarité : elle évalue 4σεT³ à la température courante, au point de Gauss, quand on la demande — comme le D_alg plastique, et pour la même raison : personne d’autre ne la lirait.

Deux natures

Le rayonnement déclare [Thermal, Radiation]. Un bord rayonnant fait partie du problème thermique — filter("thermal") doit le rendre — tandis que filter("radiation") isole le terme non linéaire à part, pour l’assembler ou l’inspecter seul. C’est le premier usage du caractère ensembliste de physics().

Unités

sigma vaut par défaut la constante SI, et T est alors une température absolue (Kelvin) : une puissance quatrième n’a aucune invariance permettant de translater une origine. Dans un autre système d’unités, fournir sigma comme composante matériau.

conduction = pyrucast.model.heat_conduction(volume)
model = conduction | pyrucast.model.radiation(bord, conduction)
materials = pyrucast.element_field.material_field(
    model, [("k", 20.0), ("emis", 0.8), ("T_inf", 300.0)]
)

Ce que ça vaut comme vérification

Deux choses se contrôlent sans acrobatie analytique : le flux rayonné doit valoir exactement σε(T⁴ − T_∞⁴) fois l’aire, et la tangente doit être la dérivée du résidu. Une loi en T⁴ est précisément là où une tangente incohérente se cache — Newton ramperait au lieu de converger quadratiquement — d’où sa comparaison à une différence finie :

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::element_field::ElementField;
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::model::Model;
use pyrucast::containers::node_field::{NodeField, SubNodeField};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::radiation::STEFAN_BOLTZMANN;
use pyrucast::models::Physics;
use pyrucast::ops::element_field;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

const EMIS: f64 = 0.8; // emissivity
const T_INF: f64 = 300.0; // far-field temperature (K)
const T_WALL: f64 = 500.0; // temperature imposed on the far side (K)

#[test]
fn the_radiated_flux_matches_stefan_boltzmann() -> Result<()> {
    let (fixture, materials) = radiating_square()?;

    // The boundary is at the uniform wall temperature: interpolate it to the
    // Gauss points, integrate the law, and scatter back to the nodes.
    let temperature = uniform_temperature(&fixture, T_WALL)?;
    let at_gauss = element_field::interp_to_gauss(&temperature, &fixture.boundary_fes)?;
    let state =
        element_field::behavior::integrate(&fixture.radiation, &at_gauss, None, &materials, None)?;
    let reaction = pyrucast::ops::node_field::internal_forces(
        &fixture.radiation,
        &state,
        &temperature,
        &materials,
    )?;

    // The radiating edge has unit length, so the total flux is the density.
    let expected = STEFAN_BOLTZMANN * EMIS * (T_WALL.powi(4) - T_INF.powi(4));
    let total: f64 = fixture
        .edge
        .iter()
        .map(|n| reaction.value(n.id(), "q").unwrap_or(0.0))
        .sum();
    assert!(
        (total - expected).abs() < 1e-9 * expected.abs(),
        "radiated flux {total}, expected {expected}"
    );
    Ok(())
}

Convection de surface (Robin / film)

Le modèle BoundaryTransfer (src/models/boundary_transfer.rs) ajoute un échange convectif avec un fluide à température ambiante \(T_\text{ext}\) sur un bord : la loi de Newton du refroidissement

\[ q\cdot n = h\,\big(T - T_\text{ext}\big) \]

où \(h\) est le coefficient d’échange (film). Injectée dans le terme de bord de la forme faible de la conduction, elle se scinde en deux ingrédients :

\[ \underbrace{K_{ij} = h \int_\Gamma N_i\,N_j\,d\Gamma}{\text{matrice de film (raideur)}} \qquad \underbrace{f_i = h\,T\text{ext} \int_\Gamma N_i\,d\Gamma}_{\text{charge (second membre)}} \]

On le construit contre la conduction qu’il refroidit, en lui passant les couples de variables à échanger — ceux de la conduction, ce qui fait que le terme se couple directement dans sa raideur. La conduction, elle, lui donne sa nature thermique, et refuse un couple qu’elle n’assemble pas :

conduction = pyrucast.model.heat_conduction(bord_fes)
film = pyrucast.model.boundary_transfer(bord_fes, conduction, [("T", "q")])

Ce modèle n’a rien de thermique : la même loi décrit un transfert de masse en surface ou une fondation élastique, selon les composantes qu’on lui donne, et il partage son noyau avec le transfert d’interface. Voir Échanges pour la loi commune, la structure en quatre blocs et le choix entre un échange et une contrainte.

Mise en donnée. Le modèle fournit la matrice de film ; la part externe \(h\,T_\text{ext}\) est un chargement, bâti avec le même opérateur flux que la source (densité \(h\,T_\text{ext}\)). Le terme de film rend la matrice définie : un problème purement Neumann + convection est bien posé sans Dirichlet.

Exemple. Une dalle \([0,1]^2\) chauffée par un flux \(Q\) sur le bord gauche et refroidie par convection sur le bord droit (haut/bas isolés). Tout le flux ressort par convection, d’où le profil linéaire

\[ T(x) = T_\text{ext} + \frac{Q}{h} + \frac{Q}{k}\,(1 - x). \]

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

#[test]
fn thermal_convection_recovers_analytical_solution() -> Result<()> {
    // ── Problem data ───────────────────────────────────────────────────────
    const K: f64 = 2.0; // conductivité
    const Q: f64 = 10.0; // densité de flux injectée sur le bord gauche
    const H: f64 = 5.0; // coefficient d'échange (film) sur le bord droit
    const T_EXT: f64 = 20.0; // température ambiante du fluide
    const N: usize = 4; // N×N éléments QUA4
    let step = 1.0 / N as f64;

    // ── Mesh: a structured (N+1)×(N+1) grid of QUA4 on [0,1]² ──────────────
    let coords = Handle::new(Coords::new(2)?);
    let idx = |i: usize, j: usize| j * (N + 1) + i; // nœud colonne i, ligne j
    let mut grid: Vec<Node> = Vec::with_capacity((N + 1) * (N + 1));
    for j in 0..=N {
        for i in 0..=N {
            grid.push(Node::create_in(
                coords.clone(),
                &[i as f64 * step, j as f64 * step],
            )?);
        }
    }
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::QUA4));
    for j in 0..N {
        for i in 0..N {
            mesh.add_cell(&[
                grid[idx(i, j)].id(),
                grid[idx(i + 1, j)].id(),
                grid[idx(i + 1, j + 1)].id(),
                grid[idx(i, j + 1)].id(),
            ])?;
        }
    }
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Modèle : conduction (volume) + convection (bord droit x = 1) ───────
    // Le bord droit est un maillage SEG2 bâti sur les nœuds de la grille ;
    // it integrates as a line (film matrix h ∫ N_i N_j dΓ).
    let mut right_edge = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for j in 0..N {
        right_edge.add_cell(&[grid[idx(N, j)].id(), grid[idx(N, j + 1)].id()])?;
    }
    let right_fes = FiniteElementSpace::lagrange1(&right_edge)?;

    let conduction = model::heat_conduction(&fes)?;
    let convection =
        model::boundary_transfer(&right_fes, &conduction, vec![("T".into(), "q".into())])?;
    let model = conduction.union(&convection)?;

    // Material: k for the conduction, h and the ambient for the convection (each
    // sub-model takes the component it requires from the supplied list).
    // ── Chargement ─────────────────────────────────────────────────────────
    // Source: uniform flux (density Q) on the left edge, as nodal loads
    // cohérentes via `flux`.
    let mut left_edge = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for j in 0..N {
        left_edge.add_cell(&[grid[idx(0, j)].id(), grid[idx(0, j + 1)].id()])?;
    }
    let left_fes = FiniteElementSpace::lagrange1(&left_edge)?;
    let model = model.union(&model::flux(&left_fes, &model, "q".into())?)?;

    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[("k", K), ("h_T", H), ("a_ext_T", T_EXT), ("phi_q", Q)],
    )?;

    // Both given terms — the left edge's source and the convection's external
    // part h·T_ext — belong to the model, which returns them together. Nothing
    // left to union by hand, hence nothing left to forget.
    let rhs = pyrucast::ops::node_field::external_forces(&model, &materials)?;

    // ── Assembly + solve (K made definite by the film term) ────────────────
    let stiffness = pyrucast::ops::matrix::stiffness(&model, &materials)?;
    let solution = solve(&stiffness, &rhs)?;

    // ── Compared with the analytical T(x) = T_ext + Q/h + (Q/k)(1 − x), ∀ y ─
    let tol = 1e-9;
    for j in 0..=N {
        for i in 0..=N {
            let x = i as f64 * step;
            let expected = T_EXT + Q / H + (Q / K) * (1.0 - x);
            let got = solution.value(grid[idx(i, j)].id(), "T")?;
            assert!(
                (got - expected).abs() < tol,
                "T(x={x}, y={}) : obtenu {got}, attendu {expected}",
                j as f64 * step
            );
        }
    }

    // Energy balance: all the injected flux leaves by convection, so the right
    // edge's temperature is exactly T_ext + Q/h.
    let t_right = solution.value(grid[idx(N, 0)].id(), "T")?;
    assert!(
        (t_right - (T_EXT + Q / H)).abs() < tol,
        "T(x=1) : obtenu {t_right}, attendu {}",
        T_EXT + Q / H
    );

    Ok(())
}

Version Python : examples/thermal_convection_2d.py (lancer avec python examples/thermal_convection_2d.py après maturin develop).

Diffusion (loi de Fick)

Introduction

La diffusion d’une espèce dans un milieu — humidité dans un béton, hydrogène dans un acier, chlorures dans un enrobage — obéit à la loi de Fick, dont l’opérateur est celui de la conduction thermique.

C’est pourtant une physique distincte, et pyrucast la traite comme telle. La variable primale est la concentration c, la duale le flux de matière j, et sa nature déclarée est Physics::Diffusion. Partager un opérateur n’est pas partager une physique : dans un problème couplé thermo-diffusif, on doit pouvoir écrire model.filter("diffusion") et n’obtenir que la partie diffusive, sans traîner la thermique avec.

Le modèle vit sur n’importe quel espace EF volumique (2-D ou 3-D, linéaire ou quadratique), avec un degré de liberté scalaire par nœud.

Équations continues résolues

La première loi de Fick relie le flux au gradient de concentration :

\[ \mathbf j = -\,\mathsf D\,\nabla c \]

et la conservation de l’espèce, en régime transitoire avec un coefficient de stockage \( \varphi \) (la porosité, pour une espèce diffusant dans un solide poreux) :

\[ \varphi\,\frac{\partial c}{\partial t} + \nabla\!\cdot\mathbf j = 0 \qquad\Longleftrightarrow\qquad \varphi\,\frac{\partial c}{\partial t} - \nabla\!\cdot(\mathsf D\,\nabla c) = 0 . \]

En stationnaire c’est l’équation de Laplace pondérée par \( \mathsf D \). C’est la même équation que la conduction thermique, à un changement de noms près (\( c \leftrightarrow T \), \( \mathsf D \leftrightarrow k \), \( \varphi \leftrightarrow \rho c_p \)) — d’où un modèle qui partage la totalité du noyau de conduction thermique, et n’en diffère que par ses variables et sa nature physique.

La forme faible, après intégration par parties, s’écrit : trouver c tel que pour tout δc admissible,

\[ \int_\Omega \varphi\,\delta c\,\dot c\; d\Omega

• \int_\Omega \nabla \delta c \cdot \mathsf D\,\nabla c\; d\Omega = -\int_{\partial\Omega} \delta c\;\mathbf j\!\cdot\!\mathbf n\; d\Gamma . \]

Forme discrétisée

\[ K_{ij} = \int_\Omega \nabla N_i^\top\,\mathsf D\;\nabla N_j\; d\Omega \quad \text{(rigidité de diffusion — Cast3M \texttt{COND})}, \] \[ C_{ij} = \int_\Omega \varphi\,N_i\,N_j\; d\Omega \quad \text{(stockage — Cast3M \texttt{CAPA})}. \]

D est un tenseur, dont le cas isotrope D = D·I redonne le produit scalaire habituel ∇N_i · ∇N_j. Les trois symétries matériau décrites au chapitre Élasticité orthotrope s’appliquent identiquement, avec un tenseur d’ordre 2 au lieu de 4 :

Le repère est donné par les vecteurs V1X, V1Y (2-D) ou V1X…V1Z, V2X…V2Z (3-D), exactement comme en mécanique.

Variables et matériau

Le comportement (COMP) rend le flux sous forme faible D·∇c, en composantes j_x, j_y(, j_z). Comme en thermique, c’est l’opposé du flux physique de Fick : ce choix garantit ∫ Bᵀ·j = K·c, donc l’accord entre le comportement et la rigidité dans le cas linéaire. Les composantes sont nommées d’après la variable duale (j_*) et non flux_*, afin qu’un modèle portant à la fois conduction et diffusion garde deux champs de flux non ambigus.

L’entrée du comportement est le gradient grad_c_x, …, tel que le produit l’opérateur gradient sur un champ dont la composante est c.

Mise en donnée (Rust, testé)

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::node_field::{NodeField, SubNodeField};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::fick::{dual_var, primal_var};
use pyrucast::models::Physics;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

/// The diffusing species — every name of this physics carries it.
const SPECIES: &str = "H2";

#[test]
fn fick_line_recovers_the_linear_profile() -> Result<()> {
    const D: f64 = 2.0; // diffusivity
    const J: f64 = 10.0; // injected species flux at x = 0
    const C_IMPOSED: f64 = 1.0; // concentration imposed at x = 1
    const N_ELEMS: usize = 4;
    let h = 1.0 / N_ELEMS as f64;

    // ── Mesh: a line of SEG2 on [0, 1] ─────────────────────────────────────
    let coords = Handle::new(Coords::new(1)?);
    let nodes: Vec<Node> = (0..=N_ELEMS)
        .map(|i| Node::create_in(coords.clone(), &[i as f64 * h]))
        .collect::<Result<_>>()?;
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for i in 0..N_ELEMS {
        mesh.add_cell(&[nodes[i].id(), nodes[i + 1].id()])?;
    }
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Modèle : diffusion + Dirichlet c = 1 en x = 1 ──────────────────────
    let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(std::slice::from_ref(
        nodes.last().unwrap(),
    ))?);
    let multiplier = mesh::barycenter(&imposed)?;
    let mult = multiplier.node(0, 0, 0)?.id();

    let diffusion = model::fick(&fes, SPECIES)?;
    let dirichlet = model::dirichlet(
        &diffusion,
        &primal_var(SPECIES),
        &imposed,
        &multiplier,
        Default::default(),
    )?;
    let model = diffusion.union(&dirichlet)?;

    // ── Matériau : diffusivité uniforme ────────────────────────────────────
    let materials =
        pyrucast::ops::element_field::material_field(&model, &[(&format!("D_{SPECIES}"), D)])?;

    // ── Loading: flux J at x = 0, imposed concentration at the multiplier
    let node0 = nodes[0].id();
    let mut load_sm = SubMesh::new(coords.clone(), ElementType::POI1);
    load_sm.add_cell(&[node0])?;
    load_sm.add_cell(&[mult])?;
    let load_sm = Handle::new(load_sm);
    let mut rhs = SubNodeField::from_poi1(
        &load_sm,
        vec![
            format!("imposed_{}", primal_var(SPECIES)),
            dual_var(SPECIES),
        ],
    )?;
    rhs.set_value(node0, &dual_var(SPECIES), J)?;
    rhs.set_value(mult, &format!("imposed_{}", primal_var(SPECIES)), C_IMPOSED)?;
    let rhs = NodeField::from_sub(rhs);

    // ── Assemblage + résolution ────────────────────────────────────────────
    let stiffness = pyrucast::ops::matrix::stiffness(&model, &materials)?;
    let solution = solve(&stiffness, &rhs)?;

    // ── Compared with the analytical profile c(x) = 1 + (J/D)(1 − x) ───────
    let tol = 1e-10;
    for (i, node) in nodes.iter().enumerate() {
        let x = i as f64 * h;
        let expected = C_IMPOSED + (J / D) * (1.0 - x);
        let got = solution.value(node.id(), &primal_var(SPECIES))?;
        assert!(
            (got - expected).abs() < tol,
            "c(x={x}) : {got} ≠ {expected}"
        );
    }
    // Mass balance: the reaction at the imposed edge balances the injected flux.
    let reaction = solution.value(mult, &format!("lambda_{}", primal_var(SPECIES)))?;
    assert!((reaction - J).abs() < tol, "réaction λ : {reaction} ≠ {J}");
    Ok(())
}

Exemple Python

"""1-D Fick diffusion — a bar fed with a species, compared with the analytical solution.

Problème
--------
On the segment [0, 1]:

  * en x = 0 : un **flux d'espèce** imposé ``J`` (Neumann) ;
  * at x = 1: an **imposed concentration** ``c = 1`` (Dirichlet).

In the steady regime without a volume source the profile is linear ::

    c(x) = 1 + (J / D) * (1 - x)

and the Lagrange multiplier at the imposed node is exactly ``J``: everything
entering at x = 0 leaves at x = 1 (mass balance).

The operator is the thermal conduction one; what changes is the **physics**.
The primal is the concentration ``c``, the dual the flux ``j``, and the kind
declared is ``"diffusion"`` — so that a coupled thermo-diffusive model splits
with ``model.filter(...)``, which the second part
de l'exemple montre.

This is the Python equivalent of the Rust integration test ``tests/fick.rs``.

Lancement
---------
Once the extension is built in the venv ::

    maturin develop --features extension-module
    python examples/diffusion_1d.py
"""

import pyrucast

# ── Problem data ────────────────────────────────────────────────────────────
SPECIES = "H2"  # the diffusing species — every name carries it
D = 2.0  # diffusivité
J = 10.0  # flux d'espèce injecté en x = 0
C_IMPOSED = 1.0  # concentration imposée en x = 1
N_ELEMS = 4
K = 5.0  # thermal conductivity, for the coupled part


def ligne(n_elems):
    """A line of ``n_elems`` SEG2 on [0, 1], with its nodes."""
    c = pyrucast.Coords(1)
    h = 1.0 / n_elems
    nodes = [c.add_node([i * h]) for i in range(n_elems + 1)]
    mesh = pyrucast.Mesh(c, "SEG2")
    for i in range(n_elems):
        mesh.unit().add_cell([nodes[i], nodes[i + 1]])
    return c, nodes, pyrucast.FiniteElementSpace(mesh), h


def profil_stationnaire() -> None:
    c, nodes, fes, h = ligne(N_ELEMS)

    # ── Modèle : diffusion + concentration imposée en x = 1 ──────────────────
    imposed = pyrucast.Mesh(c, "POI1")
    imposed.unit().add_cell([nodes[-1]])
    multiplier = pyrucast.mesh.barycenter(imposed)
    mult = multiplier.node(0, 0, 0)

    cible = pyrucast.model.fick(fes, SPECIES)

    model = cible | pyrucast.model.dirichlet(cible, f"c_{SPECIES}", imposed, multiplier)
    materials = pyrucast.element_field.material_field(model, [(f"D_{SPECIES}", D)])

    # ── Loading: flux J at x = 0, imposed value at the multiplier ────────────
    load = pyrucast.Mesh(c, "POI1")
    load.unit().add_cell([nodes[0]])
    load.unit().add_cell([mult])
    rhs = pyrucast.NodeField(load, [f"imposed_c_{SPECIES}", f"j_{SPECIES}"])
    rhs[0].set_value(nodes[0], f"j_{SPECIES}", J)
    rhs[0].set_value(mult, f"imposed_c_{SPECIES}", C_IMPOSED)

    # ── Assemblage + résolution ─────────────────────────────────────────────
    stiffness = pyrucast.matrix.stiffness(model, materials)
    solution = pyrucast.solver.solve(stiffness, rhs)

    print("Diffusion de Fick 1-D")
    print(f"  D = {D}, flux injecté J = {J}, c(1) = {C_IMPOSED}")
    print()
    print("     x      c calculé    c analytique")
    print("  " + "-" * 36)
    for i, node in enumerate(nodes):
        x = i * h
        attendu = C_IMPOSED + (J / D) * (1.0 - x)
        obtenu = solution.value(node, f"c_{SPECIES}")
        print(f"  {x:5.3f}   {obtenu:10.6f}   {attendu:12.6f}")
        assert abs(obtenu - attendu) < 1e-10

    reaction = solution.value(mult, f"lambda_c_{SPECIES}")
    print()
    print(f"  Bilan de matière : réaction = {reaction:.6f}, flux injecté = {J}")
    assert abs(reaction - J) < 1e-10


def couplage_avec_la_thermique() -> None:
    """Diffusion and conduction on the same mesh: two distinct physics."""
    _c, _nodes, fes, _h = ligne(3)
    model = pyrucast.model.fick(fes, SPECIES) | pyrucast.model.heat_conduction(fes)

    # A single material field carries both sets: the assembler resolves each zone
    # through the components its physics requires (`D` here, `k` there).
    materials = pyrucast.element_field.material_field(
        model, [(f"D_{SPECIES}", D), ("k", K)]
    )
    pyrucast.matrix.stiffness(model, materials)

    print()
    print("Modèle couplé diffusion + thermique")
    print(f"  sous-modèles          : {len(model)}")
    print(f"  filter('diffusion')   : {len(model.filter('diffusion'))}")
    print(f"  filter('thermal')     : {len(model.filter('thermal'))}")
    print(f"  filter('mechanical')  : {len(model.filter('mechanical'))}")
    assert len(model.filter("diffusion")) == 1
    assert len(model.filter("thermal")) == 1
    assert len(model.filter("mechanical")) == 0


def main() -> None:
    profil_stationnaire()
    couplage_avec_la_thermique()


if __name__ == "__main__":
    main()

Compléments

Coexister avec la thermique

Les deux physiques peuvent vivre sur le même maillage sans se gêner :

model = pyrucast.model.fick(fes, "H2") | pyrucast.model.heat_conduction(fes)
materials = pyrucast.element_field.material_field(model, [("D_H2", 2.0), ("k", 5.0)])
k = pyrucast.matrix.stiffness(model, materials)

len(model.filter("diffusion"))  # 1
len(model.filter("thermal"))  # 1

Un seul champ matériau porte les deux jeux de coefficients. L’assembleur résout la zone de chaque physique par les composantes qu’elle exige (D ici, k là) — il n’y a rien à consolider à la main. Et parce que les deux natures sont distinctes, filter les sépare de nouveau après coup, aussi bien sur le modèle que sur la matrice assemblée.

Les degrés de liberté restent séparés (c d’un côté, T de l’autre) : le système est bloc-diagonal. Un vrai couplage — une diffusivité fonction de la température, ou une thermodiffusion — se pilote depuis Python, en réassemblant la partie diffusive à chaque pas avec un champ matériau recalculé.

Transfert à travers une interface

Deux corps qui se touchent ne partagent pas forcément leurs nœuds. Un contact imparfait, un revêtement, un joint, une membrane laissent le champ sauter à la traversée, tandis qu’un flux la franchit proportionnellement à ce saut :

\[ j\cdot n = h\,\big(c_1 - c_2\big) \]

h_c_H2 est le coefficient de transfert (son inverse est la résistance de contact) : un par grandeur transférée, nommé d’après elle.

Ce modèle n’a rien de diffusif non plus : on lui passe [("T", "q")] et une conduction pour cible pour une résistance de contact, les couples de déplacement et une élasticité pour un joint collé de raideur finie — la nature vient de la cible. La loi commune, sa structure en quatre blocs dont deux hors-diagonale, l’exigence de conformité des deux côtés et le critère qui départage un échange d’une contrainte MPC sont dans Échanges.

corps = pyrucast.model.fick(gauche, "H2") | pyrucast.model.fick(droite, "H2")
model = corps | pyrucast.model.interface_transfer(
    face_gauche, face_droite, corps, [("c_H2", "j_H2")]
)
materials = pyrucast.element_field.material_field(
    model, [("D_H2", 2.0), ("h_c_H2", 5.0)]
)

Ce que ça vaut comme vérification

Deux carrés côte à côte, un flux q injecté d’un côté, la concentration imposée de l’autre : le profil est linéaire par morceaux avec une chute q/D dans chaque carré et un saut q/h à l’interface. C’est ce saut qui distingue une interface d’un nœud partagé, et il est porté entièrement par les blocs hors-diagonale. Quand h → ∞, le saut s’efface et l’on retrouve le corps continu.

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::element_field::ElementField;
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::model::Model;
use pyrucast::containers::node_field::{NodeField, SubNodeField};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::interface_transfer::DEFAULT_TOL;
use pyrucast::models::Physics;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

/// The diffusing species — every name of the Fick physics carries it.
const SPECIES: &str = "H2";

const D: f64 = 2.0; // diffusivity, both squares
const Q: f64 = 10.0; // flux density injected at x = 0
const C_RIGHT: f64 = 1.0; // concentration imposed at x = 2

#[test]
fn an_interface_law_makes_the_field_jump() -> Result<()> {
    const H: f64 = 5.0; // transfer coefficient
    let (geom, solution) = solve_two_squares(H)?;

    let c = |n: &Node| solution.value(n.id(), &format!("c_{SPECIES}"));
    let far_left = c(&geom.left[0])?; // (0, 0)
    let left_face = c(&geom.left[1])?; // (1, 0), left side
    let right_face = c(&geom.right[0])?; // (1, 0), right side
    let far_right = c(&geom.right[1])?; // (2, 0)

    let tol = 1e-10;
    assert!((far_right - C_RIGHT).abs() < tol, "c(2) = {far_right}");
    // Slope q/D over the unit width of each square.
    assert!(
        (far_left - left_face - Q / D).abs() < tol,
        "left square: {far_left} → {left_face}"
    );
    assert!(
        (right_face - far_right - Q / D).abs() < tol,
        "right square: {right_face} → {far_right}"
    );
    // …and the jump across the interface is q/h — the exchange law itself.
    let jump = left_face - right_face;
    assert!(
        (jump - Q / H).abs() < tol,
        "jump = {jump}, expected {}",
        Q / H
    );
    Ok(())
}

Régime transitoire

La matrice de stockage s’assemble avec matrix.mass(...), qui exige alors la composante poro. L’intégration en temps est orchestrée en Python, comme pour la thermique transitoire — le noyau Rust fournit K et C, pas la boucle.

Bilan de matière

Comme en thermique, le multiplicateur de Lagrange d’une concentration imposée est le flux d’espèce qui traverse la frontière. C’est la vérification la plus directe d’un calcul de diffusion, et c’est ce que contrôle le test d’intégration ci-dessus : la réaction au bord imposé égale exactement le flux injecté à l’autre bout.

Échanges (frontière et interface)

Introduction

Deux physiques, une seule loi. Un échange fait traverser un flux proportionnellement à un écart :

\[ q\cdot n = h\,\big(a - b\big) \]

et toute la différence entre les deux tient à de quel côté du signe égal vit le milieu d’en face :

Le bloc hors-diagonale est le second membre rendu implicite. C’est pourquoi les deux partagent leur noyau (src/models/transfer.rs) au lieu d’en porter chacun une copie : le terme de frontière est le terme d’interface avec les deux côtés sur la même maille.

Ce qui est transféré appartient à l’appelant

Aucun des deux ne sait ce qu’il transporte. On lui donne une liste de couples (primale, duale) — la forme que prennent déjà embedded et contact, les deux autres lois qui lient des maillages — et tout le reste en découle :

Un coefficient par grandeur transférée, nommé d’après elle. C’est ce qui permet une raideur par direction, et ce qu’un h unique ne pouvait pas exprimer.

Nommer les DDL de la physique de volume est aussi ce qui fait que le terme s’y couple directement : un boundary_transfer sur ("T", "q") entre dans la raideur d’un heat_conduction sans adaptateur, parce que la matrice est indexée par le couple (nœud, nom de champ).

La cible donne la nature

Des noms de variables libres ne disent pas quelle nature physique porte l’échange ("thermal", "diffusion", "mechanical") — c’est pourtant elle que model.filter(...) sélectionne. Le modèle dont l’échange nomme les DDL, lui, la connaît. Les deux lois se construisent donc contre ce modèle, passé en argument target après le support, exactement comme une charge répartie ou une contrainte :

conduction = model.heat_conduction(volume)
film       = model.boundary_transfer(peau, conduction, [("T", "q")])      # thermique
corps      = model.fick(gauche, "H2") | model.fick(droite, "H2")
joint      = model.interface_transfer(face_g, face_d, corps, [("c_H2", "j_H2")])  # diffusion

Chaque couple doit être assemblé par la cible : un de ses sous-modèles compte la primale parmi ses inconnues et l’apparie à la duale. La vérification se fait une fois, à la construction, et refuse trois erreurs qui passaient jusqu’ici :

Pour deux natures, on construit deux échanges. Le rayonnement reçoit lui aussi sa cible — il écrit dans la ligne q d’une conduction —, mais ses natures restent les siennes, [Thermal, Radiation] : la cible n’y sert qu’à prouver que la ligne est assemblée.

Forme discrétisée

Échange avec une ambiante

La forme faible du terme de bord se scinde en deux ingrédients, dont le sous-modèle ne porte que le premier :

\[ \underbrace{K_{ij} = h \int_\Gamma N_i\,N_j\,d\Gamma}{\text{matrice de film (raideur)}} \qquad \underbrace{f_i = h\,a\text{ext} \int_\Gamma N_i\,d\Gamma}_{\text{charge (second membre)}} \]

La part ambiante est un chargement, et le sous-modèle le porte : a_ext est une composante matériau exigée, à côté de son coefficient, et le terme se récupère par external_forces. Elle a longtemps été bâtie à la main avec l’opérateur flux, ce qui laissait l’oublier — et un ambiant oublié ne se voit pas : il se lit comme un ambiant nul, donc comme une paroi qui échange avec le zéro absolu. Le terme de film rend par ailleurs la matrice définie : un problème purement Neumann + échange est bien posé sans Dirichlet.

Aucune normale à choisir. Elle est déjà consommée en passant de \( q\cdot n \) à \( h(a - a_\text{ext}) \) ; ce qui reste sous l’intégrale est un scalaire fois la mesure \( d\Gamma = \sqrt{\det(J^\top J)} \), une magnitude indépendante de l’orientation du maillage de bord — contrairement à une pression ou à un flux signé.

Échange entre deux maillages

Ni l’une ni l’autre n’a de loi de comportement. Le h·a d’un bord et le h·(a₁−a₂) d’une interface sont le coefficient de leur propre opérateur appliqué en un point, pas un comportement : le même h a déjà bâti ∫h NᵀN. Elles ont longtemps dû s’en inventer une, parce que Domain en exigeait une ; depuis que la loi est une capacité à part, elles n’en déclarent plus.

Côté résidu, l’intégrale d’une différence est la différence des intégrales : l’interface rend ∫h(a₁−a₂)N sur A et son opposé sur B, chacune dispersée sur son propre espace. Ce qui est couplé, c’est la matrice — lignes sur A, colonnes sur B — pas le vecteur, qui ne produit qu’un nombre par nœud. La conformité exigée des deux côtés sert à savoir où lire a₂ : maille pour maille, les deux interpolations s’alignent indice pour indice et le saut est une soustraction. Sans elle il faudrait localiser les points de Gauss de l’un dans l’autre — ce qui n’est plus la même méthode, mais un mortar.

Quand l’autre côté est une inconnue, la forme faible devient

\[ \int_\Gamma h\,(a_1 - a_2)\,(\delta a_1 - \delta a_2)\; d\Gamma, \]

qui se développe en une structure \( 2\times2 \) sur les DDL des deux côtés :

\[ \begin{bmatrix} +K & -K \\ -K & +K \end{bmatrix}, \qquad K_{ij} = h \int_\Gamma N_i\,N_j\; d\Gamma . \]

Les deux blocs diagonaux sont exactement deux termes de frontière, un par côté. Les deux autres ont leurs lignes sur un maillage et leurs colonnes sur l’autre : c’est le genre de contribution Coupling, dont ce modèle est le seul utilisateur (voir Ajouter une physique).

Les quatre sortent d’une seule fonction, transfer::exchange_matrix : les diagonaux avec la même maille des deux côtés et le signe +, les croisés avec la maille en vis-à-vis et le signe −. Le signe est donc porté par le noyau et non par un facteur qu’il faudrait faire circuler dans l’assembleur. Chaque bloc de couplage pris seul est non symétrique ; leur réunion l’est — exactement comme la paire C / Cᵀ de Dirichlet.

Les grandeurs transférées ne se couplent pas entre elles : la chaleur qui traverse un joint n’y pousse pas l’hydrogène. Seule la diagonale en indice de variable est écrite, donc le coût reste linéaire en nombre de composantes.

Conformité

Les deux côtés doivent être conformes : même type d’élément, même nombre de mailles, maille i face à maille i, et nœud local k face au nœud local k. C’est vérifié géométriquement à la construction — les nœuds appariés doivent être colocalisés — et signalé plutôt qu’approché. Une interface non conforme est un problème de maillage ; la rattraper par une projection silencieuse fabriquerait des flux faux sans le dire.

Les deux maillages sont donc géométriquement confondus mais numérotés séparément : c’est cette duplication des nœuds qui laisse le champ sauter, là où un nœud partagé l’interdirait.

Mise en donnée (Rust, testé)

/// A boundary exchange on the displacement is an **elastic foundation**: the
/// face rests on a distributed spring, and the analytical `u = q/(E/L + h)`
/// comes out for stiffnesses spanning four decades.
#[test]
fn a_boundary_exchange_on_displacements_is_an_elastic_foundation() -> Result<()> {
    for h in [0.0, 1e3, 1e4, 1e5, 1e6] {
        let u = free_face_displacement(4, h)?;
        let exact = Q / (E / L + h);
        assert!(
            (u - exact).abs() < 1e-9 * exact.abs().max(1e-12),
            "h = {h}: face displacement {u}, exact {exact}"
        );
    }
    Ok(())
}

Exemple Python

# Thermal film: enters the stiffness of the conduction it cools, and
# en prend la nature.
conduction = pyrucast.model.heat_conduction(fes)
model = conduction | pyrucast.model.boundary_transfer(peau, conduction, [("T", "q")])
materials = pyrucast.element_field.material_field(
    model, [("k", 5.0), ("h_T", 12.0), ("a_ext_T", 20.0)]
)

# Contact resistance between two meshes.
corps = pyrucast.model.heat_conduction(gauche) | pyrucast.model.heat_conduction(droite)
joint = pyrucast.model.interface_transfer(face_gauche, face_droite, corps, [("T", "q")])

# Elastic foundation: the same law, on displacements.
plaque = pyrucast.model.elasticity(fes, "plane_stress")
appui = pyrucast.model.boundary_transfer(
    semelle, plaque, [("u_x", "f_x"), ("u_y", "f_y")]
)

# No kind in argument: each took its target's own.
joint[0].physics()  # ["thermal"]
appui[0].physics()  # ["mechanical"]

Compléments

Ce que la généralisation apporte

Rien dans la loi n’est thermique ni diffusif, donc deux modèles sortent sans une ligne de mécanique :

• une fondation élastique de Winkler — un échange de frontière sur les déplacements. Une barre poussée sur sa face libre, cette face reposant sur un ressort réparti, c’est la barre et le ressort en parallèle sous la traction appliquée : \( u = q / (E/L + h) \), ce que le test vérifie à 1e-9 sur quatre décades de raideur ;
• un joint collé de raideur finie — le même sur une interface.

Échange ou contrainte ?

Un interface_transfer est la régularisation par pénalité de ce qu’un MPC impose exactement :

mpc([(T, maillage₁, +1), (T, maillage₂, −1)]) = 0   ⟺  T₁ = T₂, exactement
interface_transfer(maillage₁, maillage₂), h → ∞      ⟶  T₁ = T₂, à 1/h près

Le critère de choix tient en une phrase : si h vient d’une mesure, c’est de la physique ; s’il a été choisi « assez grand », il fallait une contrainte. Un joint imparfait, une résistance thermique de contact, une couche adhésive ont un h mesurable ; « lier deux surfaces » n’en a pas.

Ce que ça vaut comme vérification

Le film thermique est confronté à une solution analytique — une dalle chauffée d’un côté, refroidie de l’autre, dont tout le flux ressort par convection (voir Conduction thermique).

L’interface l’est au saut : deux carrés côte à côte, un flux q injecté d’un côté, la concentration imposée de l’autre, et un saut q/h à la traversée (voir Diffusion). C’est ce saut qui distingue une interface d’un nœud partagé, et il est porté entièrement par les blocs hors-diagonale ; quand h → ∞ il s’efface et l’on retrouve le corps continu.

S’y ajoute la structure : la raideur assemblée reste symétrique alors qu’aucun bloc de couplage ne l’est, ce qui est la vérification que les quatre blocs atterrissent où ils doivent.

Ce qui n’est pas couvert

Pas d’échange entre maillages non coïncidents : la liaison exacte a sa version non conforme avec embedded et ses poids d’interpolation, l’échange fini n’en a pas. Le mécanisme existe, il n’a simplement jamais été branché là.

Mécanique

Les physiques mécaniques, chacune décrite sur sa page : équations résolues, exemple de mise en donnée Rust testé, et exemple Python.

Comme pour la thermique, ce sont des variantes de SubModel : chacune déclare ses variables, son matériau, et assemble sa rigidité K (et, le cas échéant, son comportement COMP).

Convention de nommage : primal = déplacements u_x, u_y, u_z (les inconnues) ; dual = forces nodales f_x, f_y, f_z (second membre / réactions).

• Barre / treillis — élément à effort axial (1-D/2-D/3-D).
• Élasticité linéaire — continuum 2-D (CP/DP) et 3-D.
• Élasticité orthotrope et anisotrope — la symétrie matériau, repère donné par vecteurs.
• Plasticité parfaite (von Mises) — J2 sans écrouissage, retour radial (état interne εᵖ, p).
• Lois d’écoulement plastique — la loi comme attribut : écrouissage isotrope, Drucker-Prager, Ottosen.
• Fluage et viscoplasticité — les lois dépendantes du temps, qui exigent dt.
• Endommagement de Mazars — endommagement isotrope du béton, deux variables (état interne κ).
• Lois d’endommagement — la loi comme attribut : Damage TC, SiC/SiC orthotrope, Gurson.
• Poutre d’Euler-Bernoulli — sans cisaillement transverse, interpolation d’Hermite (1-D / plan / spatial).
• Poutre de Timoshenko — flexion + cisaillement, intégration réduite (anti-verrouillage).
• Portique 2D — poutre orientée (axial + flexion + cisaillement), transformation local→global.
• Coques — surface à six DDL par nœud, Reissner-Mindlin à cisaillement sous-intégré ou Kirchhoff discret (DKT/DKQ) sans cisaillement.
• Cadre 3D — space frame 6 DOF/nœud (axial + torsion + flexion 2 plans), orientation automatique.

Barre / treillis

Élément SEG2 à 2 nœuds transmettant uniquement l’effort axial (treillis). Fonctionne à l’identique en 1-D, 2-D et 3-D.

Équations continues résolues

La barre suit la loi axiale 1-D et son équilibre le long de l’axe s :

\[ N = E\,A\,\varepsilon, \qquad \varepsilon = \frac{du}{ds}, \qquad \frac{dN}{ds} + f = 0, \]

où N est l’effort normal, ε la déformation axiale (dérivée du déplacement le long de la barre) et f la charge axiale répartie. L’orientation est déduite des coordonnées des nœuds via le cosinus directeur c = (x_B − x_A)/L.

Forme discrétisée

Avec l’interpolation linéaire SEG2, la déformation est constante par élément. Projetée sur les directions physiques par c, la rigidité élémentaire globale (en d dimensions) s’écrit

\[ K_e = \frac{E\,A}{L} \begin{bmatrix} c\,c^\top & -c\,c^\top \\ -c\,c^\top & c\,c^\top \end{bmatrix}, \]

écrite aux positions (NodeId_i, f_a) × (NodeId_j, u_b). En 1-D, c = 1 et l’on retrouve (EA/L)[[1,−1],[−1,1]].

Variables et matériau

• primal : u_x, u_y(, u_z) — dual : f_x, f_y(, f_z).
• matériau : E (module d’Young), A (section) ; rho facultatif (masse).
• comportement (COMP) : effort axial N = E·A·(cᵀ ε c), à partir de la déformation ε (op deformation).

⚠️ Une barre n’a aucune raideur transversale : pour un système bien posé il faut bloquer les DOFs transverses (treillis triangulé, appuis), sinon la matrice est singulière.

Mise en donnée (Rust, testé)

Barre horizontale encastrée à gauche, appuyée transversalement à droite, force axiale F au bout ⇒ u_x = F·L/(E·A). Le code est le test d’intégration tests/truss.rs, exécuté à chaque cargo test :

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::model::Model;
use pyrucast::containers::node_field::{NodeField, SubNodeField};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

#[test]
fn truss_bar_recovers_axial_elongation() -> Result<()> {
    const E: f64 = 210.0e9; // Young's modulus (Pa)
    const A: f64 = 1.0e-4; // section area (m²)
    const L: f64 = 2.0; // length (m)
    const F: f64 = 1000.0; // axial force at the right end (N)

    // ── Mesh: one horizontal SEG2 bar ──────────────────────────────────────
    let coords = Handle::new(Coords::new(2)?);
    let n0 = Node::create_in(coords.clone(), &[0.0, 0.0])?;
    let n1 = Node::create_in(coords.clone(), &[L, 0.0])?;
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    mesh.add_cell(&[n0.id(), n1.id()])?;
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Modèle : barre + appuis (Dirichlet homogènes) ──────────────────────
    // Homogeneous (u = 0) BCs: the imposed value defaults to 0, so we only need
    // to introduce the constraint. A bar has no transverse stiffness, hence
    // `u_y` is clamped at both nodes to make the system well-posed.
    let clamp = |target: &Model, node: &Node, var: &str| -> Result<Model> {
        let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(std::slice::from_ref(node))?);
        let multiplier = mesh::barycenter(&imposed)?;
        model::dirichlet(target, var, &imposed, &multiplier, Default::default())
    };
    let mut model = model::truss(&fes)?;
    model = model.union(&clamp(&model, &n0, "u_x")?)?;
    model = model.union(&clamp(&model, &n0, "u_y")?)?;
    model = model.union(&clamp(&model, &n1, "u_y")?)?;

    // ── Matériau E, A (Dirichlet ignoré automatiquement) ───────────────────
    let materials = pyrucast::ops::element_field::material_field(&model, &[("E", E), ("A", A)])?;

    // ── Loading: axial force F at the right node ───────────────────────────
    let mut load_sm = SubMesh::new(coords.clone(), ElementType::POI1);
    load_sm.add_cell(&[n1.id()])?;
    let load_sm = Handle::new(load_sm);
    let mut rhs = SubNodeField::from_poi1(&load_sm, vec!["f_x".into()])?;
    rhs.set_value(n1.id(), "f_x", F)?;
    let rhs = NodeField::from_sub(rhs);

    // ── Assemblage + résolution ────────────────────────────────────────────
    let stiffness = pyrucast::ops::matrix::stiffness(&model, &materials)?;
    let solution = solve(&stiffness, &rhs)?;

    // ── Comparaison à l'analytique : u_x = F·L / (E·A) ─────────────────────
    let expected = F * L / (E * A);
    let ux = solution.value(n1.id(), "u_x")?;
    assert!(
        (ux - expected).abs() < 1e-10 * expected,
        "u_x = {ux}, attendu {expected}"
    );
    // The left node is clamped, the right end does not move transversally.
    assert!(solution.value(n0.id(), "u_x")?.abs() < 1e-18);
    assert!(solution.value(n1.id(), "u_y")?.abs() < 1e-18);

    Ok(())
}

Exemple Python

"""Bar / truss — a bar in tension, compared with the analytical solution.

Physique
--------
A 2-node `SEG2` element transmitting the axial force only. Law: `N = E·A·ε`
with `ε = du/ds` (axial strain along the bar). Global stiffness
`K_e = (E·A/L)·[[c⊗c, -c⊗c], [-c⊗c, c⊗c]]`, where `c` is the direction cosine
(derived from the nodes' coordinates) — works in 1-D/2-D/3-D.

Problème
--------
A horizontal bar of length `L`, clamped on the left (`u_x = u_y = 0`),
transversally supported on the right (`u_y = 0`), axial force `F` on the right.
A bar having no transverse stiffness, `u_y` is blocked at both nodes.
Solution analytique : `u_x = F·L / (E·A)`.

Lancement ::

    maturin develop --features extension-module
    python examples/truss.py
"""

import pyrucast

E, A, L, F = 210.0e9, 1.0e-4, 2.0, 1000.0


def _clamp(target, node, var):
    """Homogeneous Dirichlet (u = 0) on `var` at node `node`."""
    imposed = pyrucast.mesh.poi1_from_nodes([node])
    multiplier = pyrucast.mesh.barycenter(imposed)
    return pyrucast.model.dirichlet(target, var, imposed, multiplier)


def main() -> None:
    c = pyrucast.Coords(2)
    n0 = c.add_node([0.0, 0.0])
    n1 = c.add_node([L, 0.0])
    mesh = pyrucast.mesh.line(n0, n1, 1)  # un seul SEG2 (mailleur `line`)
    fes = pyrucast.FiniteElementSpace(mesh)

    model = pyrucast.model.truss(fes)
    model = model | _clamp(model, n0, "u_x")
    model = model | _clamp(model, n0, "u_y")
    model = model | _clamp(model, n1, "u_y")  # no transverse stiffness

    materials = pyrucast.element_field.material_field(model, [("E", E), ("A", A)])

    load = pyrucast.mesh.poi1_from_nodes([n1])
    rhs = pyrucast.NodeField(load, ["f_x"])
    rhs[0].set_value(n1, "f_x", F)

    solution = pyrucast.solver.solve(pyrucast.matrix.stiffness(model, materials), rhs)

    ux = solution.value(n1, "u_x")
    expected = F * L / (E * A)
    print(f"u_x (bout) = {ux:.6e}   (analytique F·L/E·A = {expected:.6e})")
    assert abs(ux - expected) < 1e-10 * expected
    print("OK : élongation axiale conforme à F·L/(E·A).")


if __name__ == "__main__":
    main()

Compléments

Masse & rigidité géométrique

En plus de la rigidité, la barre assemble :

• masse consistante M = (ρAL/6)[[2,1],[1,2]] sur chaque composante de translation (rho composante matériau optionnelle) — pyrucast.matrix.mass ;
• rigidité géométrique K_g = (N/L)·(I − c⊗c) transverse, sous l’effort axial N (sortie n du comportement) — pyrucast.matrix.geometric.

Élasticité linéaire

Continuum en petites déformations : 2-D (TRI3 / QUA4) ou 3-D (TET4 / HEX8), les cas 2-D couvrant aussi l’axisymétrie (solide de révolution maillé dans son plan méridien).

Équations continues résolues

Sur le domaine \( \Omega \), avec \( b \) les efforts volumiques :

\[ \underbrace{\nabla\cdot\sigma + b = 0}{\text{équilibre}}, \qquad \underbrace{\sigma = \mathbb{D} : \varepsilon}{\text{loi de Hooke}}, \qquad \underbrace{\varepsilon = \tfrac12(\nabla u + \nabla u^\top)}_{\text{cinématique}}. \]

La forme faible (multiplication par un déplacement virtuel \( v \), intégration par parties) s’écrit : trouver \( u \) tel que pour tout \( v \),

\[ \int_\Omega \varepsilon(v) : \mathbb{D} : \varepsilon(u)\,d\Omega = \int_\Omega v\cdot b\,d\Omega + \int_{\Gamma_N} v\cdot t\,d\Gamma, \]

où \( t \) est la traction imposée sur le bord de Neumann \( \Gamma_N \).

Forme discrétisée

En convention de Voigt (déformation ingénieur \( \gamma = 2\varepsilon \)), le champ discret \( u_h = \sum_i N_i u_i \) donne \( \varepsilon = B\,u_e \), avec la matrice déformation-déplacement \( B \) bâtie des dérivées physiques \( \partial N_i/\partial x_a \) (voir dn_dx). En 2-D (\( \varepsilon = [\varepsilon_{xx}, \varepsilon_{yy}, \gamma_{xy}]^\top \)), le bloc du nœud \( i \) est

\[ B_i = \begin{bmatrix} \partial_x N_i & 0 \\ 0 & \partial_y N_i \\ \partial_y N_i & \partial_x N_i \end{bmatrix}, \]

et en 3-D (\( \varepsilon = [\varepsilon_{xx}, \varepsilon_{yy}, \varepsilon_{zz}, \gamma_{yz}, \gamma_{xz}, \gamma_{xy}]^\top \)),

\[ B_i = \begin{bmatrix} \partial_x N_i & 0 & 0 \\ 0 & \partial_y N_i & 0 \\ 0 & 0 & \partial_z N_i \\ 0 & \partial_z N_i & \partial_y N_i \\ \partial_z N_i & 0 & \partial_x N_i \\ \partial_y N_i & \partial_x N_i & 0 \end{bmatrix}. \]

La rigidité élémentaire est alors, intégrée par quadrature de Gauss,

\[ K_e = \int_{\Omega_e} B^\top D\, B\, d\Omega \;\approx\; \sum_g B(\xi_g)^\top D\, B(\xi_g)\,|J(\xi_g)|\,w_g, \]

écrite aux positions (NodeId_i, f_a) × (NodeId_j, u_b) (ordre des DOFs nœud-majeur). Le second membre nodal cohérent d’une traction de bord est \( f_i = \int_{\Gamma_N} N_i\,t\,d\Gamma \) (opérateur flux).

Matrice constitutive D

Le modèle fixe \( D \) (isotrope, module d’Young \( E \), coefficient de Poisson \( \nu \)) :

• plane_stress (contraintes planes, \( \sigma_{zz}=0 \)), avec \( c = \dfrac{E}{1-\nu^2} \) :

\[ D = c\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \tfrac{1-\nu}{2} \end{bmatrix}; \]

• plane_strain (déformations planes, \( \varepsilon_{zz}=0 \), \( \sigma_{zz}\neq 0 \)), avec \( c = \dfrac{E}{(1+\nu)(1-2\nu)} \) :

\[ D = c\begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \tfrac{1-2\nu}{2} \end{bmatrix}; \]

• full_3d (3-D), même \( c \), avec le module de cisaillement \( G = c\,\tfrac{1-2\nu}{2} \) :

\[ D = \begin{bmatrix} c(1-\nu) & c\nu & c\nu & & & \\ c\nu & c(1-\nu) & c\nu & & & \\ c\nu & c\nu & c(1-\nu) & & & \\ & & & G & & \\ & & & & G & \\ & & & & & G \end{bmatrix} \quad (\text{ordre } [xx, yy, zz, yz, xz, xy]). \]

• axisymmetric (solide de révolution), même \( c \) : les trois directions normales \( r, z, \theta \) étant orthogonales, le bloc normal est l’isotrope 3×3 et \( rz \) est le seul cisaillement,

\[ D = c\begin{bmatrix} 1-\nu & \nu & \nu & \\ \nu & 1-\nu & \nu & \\ \nu & \nu & 1-\nu & \\ & & & \tfrac{1-2\nu}{2} \end{bmatrix} \quad (\text{ordre } [rr, zz, \theta\theta, rz]). \]

Axisymétrie

Un solide de révolution se maille dans son plan méridien \( (r, z) \) sur des Coords déclarées axisymétriques (\( x = r \ge 0 \), \( y = z \)) — voir Coordonnées. Deux choses changent, et elles ont deux origines distinctes :

1. la mesure d’intégration, portée par la géométrie : \( d\Omega = 2\pi r\,|J|\,d\xi \). Elle vaut pour toutes les intégrales — rigidité, masse, conductivité, flux réparti, volumes, forces internes, y compris sur les sous-maillages de bord SEG2, dont \( \int 2\pi r\,N \) donne directement l’effort sur l’anneau. Rien à écrire : c’est CellGeom::det_j_w qui l’applique, en un seul point ;
2. la déformation orthoradiale, portée par le modèle : \( \varepsilon_{\theta\theta} = u_r / r \), que le gradient méridien ne peut pas exprimer. Elle ajoute une quatrième composante de Voigt et une ligne à \( B \) :

\[ B_i = \begin{bmatrix} \partial_r N_i & 0 \\ 0 & \partial_z N_i \\ N_i / r & 0 \\ \partial_z N_i & \partial_r N_i \end{bmatrix} \quad (\varepsilon = [\varepsilon_{rr}, \varepsilon_{zz}, \varepsilon_{\theta\theta}, \gamma_{rz}]^\top). \]

Les points de Gauss étant intérieurs à la maille, \( r > 0 \) même pour un élément qui touche l’axe : le terme \( N_i/r \) reste fini, sans traitement particulier de l’axe.

Nommage (convention Cast3M) : les composantes s’appellent sigma_xx, sigma_yy, sigma_zz, sigma_xy et eps_xx, eps_yy, eps_zz, eps_xy, où zz désigne l’orthoradial \( \theta\theta \) — le plan méridien n’occupant que xx, yy et xy, il n’y a pas de collision.

Le modèle et le repère doivent s’accorder dans les deux sens : une géométrie de révolution refuse plane_stress / plane_strain, et axisymmetric refuse une géométrie cartésienne. Sans cela on mélangerait silencieusement une loi plane avec la mesure \( 2\pi r \).

La thermique n’a rien de spécifique à faire : le flux \( q = -k\nabla T \) est déjà purement méridien, et le facteur \( 2\pi r \) suffit à produire le profil logarithmique d’un cylindre creux. La plasticité et Mazars supportent l’axisymétrie, leur état interne étant déjà stocké en 3-D complet. En revanche barre et portique la refusent : un segment du plan méridien engendre une coque de révolution, que leurs noyaux ne modélisent pas.

Un maillage de bord (SEG2 en 2-D) est par ailleurs refusé comme domaine par les trois physiques de milieu continu : B y serait bâti sur le gradient tangent et \( B^\top D B \) serait déficient en rang dans la direction normale. Un bord porte des charges (flux, convection), il n’est pas un massif.

Validation : tests/axisymmetric.rs (Lamé, patch test de dilatation uniforme, \( \int B^\top\sigma = K u \), volume et masse de révolution, conduction logarithmique) et tests/python/test_axisymmetric.py.

Convergence sur Lamé

La solution de Lamé \( u_r = c_1 r + c_2/r \) comporte un terme rationnel : aucune base de Lagrange, de quelque degré que ce soit, ne la reproduit exactement. Les éléments quadratiques gagnent un ordre, pas l’exactitude. La solution ne dépendant que de \( r \), le problème discret est une EDO 1-D et les valeurs nodales sont superconvergentes en \( O(h^{2p}) \) :

(erreur relative maximale sur \( u_r \)). Les contraintes, une dérivée plus bas, passent de \( O(h) \) à \( O(h^2) \).

Le cas exact existe néanmoins : lorsque \( c_2 = 0 \) — dilatation uniforme \( u_r = c\,r \) — l’état de déformation est constant et même Q1 le reproduit à la précision machine (c’est le patch test de la suite de validation).

Matrice de masse

Pour la dynamique, la masse consistante (composante matériau rho) est

\[ M_e = \int_{\Omega_e} \rho\,N^\top N\, d\Omega \;\approx\; \sum_g \rho\,N(\xi_g)^\top N(\xi_g)\,|J(\xi_g)|\,w_g, \]

où \( N \) place \( N_i \) sur chaque composante de translation — assemblée par assemble.mass, et concentrable en diagonale par lump.

Variables et matériau

• primal : u_x, u_y(, u_z) — dual : f_x, f_y(, f_z).
• matériau : E (Young), nu (Poisson) ; facultatif alpha (dilatation thermique, cf. thermomécanique), rho (masse) — accepté par le champ matériau mais jamais exigé pour un assemblage purement élastique.
• comportement (COMP) : σ = D ε (convention tenseur → ingénieur γ = 2ε), à partir de la déformation ε (op deformation).
• modèles : plane_stress, plane_strain, axisymmetric (2-D) et full_3d (3-D).

Mise en donnée (Rust, testé)

Carré unité en contraintes planes : appuis u_x = 0 (gauche) et u_y = 0 (bas), traction S sur le bord droit appliquée en charges nodales cohérentes par l’opérateur flux (composante f_x). Solution exacte u_x = (S/E)·x, u_y = −(ν S/E)·y. Code = test tests/elasticity.rs (le fichier contient aussi un test 3-D sur un cube HEX8) :

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::model::Model;
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::tensor::Kinematics;
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

#[test]
fn elasticity_unit_square_uniaxial_tension() -> Result<()> {
    const E: f64 = 210.0; // Young's modulus
    const NU: f64 = 0.3; // Poisson's ratio
    const S: f64 = 2.0; // traction on the right edge
    const N: usize = 2; // N×N QUA4 grid
    let h = 1.0 / N as f64;

    // ── QUA4 mesh on [0,1]² ────────────────────────────────────────────────
    let coords = Handle::new(Coords::new(2)?);
    let idx = |i: usize, j: usize| j * (N + 1) + i;
    let mut grid: Vec<Node> = Vec::new();
    for j in 0..=N {
        for i in 0..=N {
            grid.push(Node::create_in(
                coords.clone(),
                &[i as f64 * h, j as f64 * h],
            )?);
        }
    }
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::QUA4));
    for j in 0..N {
        for i in 0..N {
            mesh.add_cell(&[
                grid[idx(i, j)].id(),
                grid[idx(i + 1, j)].id(),
                grid[idx(i + 1, j + 1)].id(),
                grid[idx(i, j + 1)].id(),
            ])?;
        }
    }
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Modèle : élasticité plane stress + appuis (rollers) ────────────────
    // The support receives the model it constrains: the dual is read there, and
    // contraints s'y vérifient.
    let roller = |target: &Model, nodes: &[Node], var: &str| -> Result<Model> {
        let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(nodes)?);
        let multiplier = mesh::barycenter(&imposed)?;
        model::dirichlet(target, var, &imposed, &multiplier, Default::default())
    };
    let left: Vec<Node> = (0..=N).map(|j| grid[idx(0, j)].clone()).collect();
    let bottom: Vec<Node> = (0..=N).map(|i| grid[idx(i, 0)].clone()).collect();
    let mut model = model::elasticity(&fes, Kinematics::PlaneStress)?;
    model = model.union(&roller(&model, &left, "u_x")?)?;
    model = model.union(&roller(&model, &bottom, "u_y")?)?;

    // ── Loading: traction S on the right edge (consistent nodal loads, on the
    //    f_x component). The load is a sub-model: it joins the model, its density
    //    the material. ──────────────────────────────────────────────────────
    let mut right_edge = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    for j in 0..N {
        right_edge.add_cell(&[grid[idx(N, j)].id(), grid[idx(N, j + 1)].id()])?;
    }
    let right_fes = FiniteElementSpace::lagrange1(&right_edge)?;
    model = model.union(&model::flux(&right_fes, &model, "f_x".into())?)?;

    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[("E", E), ("nu", NU), ("phi_f_x", S)],
    )?;
    let rhs = pyrucast::ops::node_field::external_forces(&model, &materials)?;

    // ── Assemblage + résolution ────────────────────────────────────────────
    let stiffness = pyrucast::ops::matrix::stiffness(&model, &materials)?;
    let solution = solve(&stiffness, &rhs)?;

    // ── Compared with the analytical u_x = (S/E)·x, u_y = −(ν S/E)·y ───────
    let tol = 1e-10;
    for j in 0..=N {
        for i in 0..=N {
            let (x, y) = (i as f64 * h, j as f64 * h);
            let ux = solution.value(grid[idx(i, j)].id(), "u_x")?;
            let uy = solution.value(grid[idx(i, j)].id(), "u_y")?;
            assert!((ux - S / E * x).abs() < tol, "u_x({x},{y}) = {ux}");
            assert!((uy + NU * S / E * y).abs() < tol, "u_y({x},{y}) = {uy}");
        }
    }
    Ok(())
}

Exemple Python

"""Élasticité linéaire — traction d'un carré (contraintes planes).

Physique
--------
Continuum en petites déformations : équilibre `∇·σ = 0`, loi de Hooke
`σ = D : ε`, kinematics `ε = ½(∇u + ∇uᵀ)`. The stiffness is
`K = ∫ Bᵀ D B dΩ` (B : matrice déformation-déplacement en Voigt, D : matrice
constitutive isotrope, ici en contraintes planes).

Problème
--------
Carré unité, appuis `u_x = 0` (bord gauche) et `u_y = 0` (bord bas), traction
`S` on the right edge applied as consistent nodal loads by the `flux` operator
(on the `f_x` component). Exact (uniaxial) solution:
`u_x = (S/E)·x`, `u_y = -(ν·S/E)·y`.

Lancement ::

    maturin develop --features extension-module
    python examples/elasticity.py
"""

import pyrucast

E, NU, S, N = 210.0, 0.3, 2.0, 2


def _clamp(target, nodes, var):
    imposed = pyrucast.mesh.poi1_from_nodes(nodes)
    multiplier = pyrucast.mesh.barycenter(imposed)
    return pyrucast.model.dirichlet(target, var, imposed, multiplier)


def main() -> None:
    h = 1.0 / N
    c = pyrucast.Coords(2)

    def idx(i, j):
        return j * (N + 1) + i

    # An N×N grid of QUA4 by sweeping two SEG2 lines (`sweep`).
    bottom = pyrucast.mesh.line(c.add_node([0.0, 0.0]), c.add_node([1.0, 0.0]), N)
    top = pyrucast.mesh.line(c.add_node([0.0, 1.0]), c.add_node([1.0, 1.0]), N)
    mesh = pyrucast.mesh.sweep(bottom, top, N)

    # Nodes laid out by idx(i, j) (i along x, j along y) by reading the
    # connectivité QUA4 : maille (cy, cx) = cy*N + cx, nœuds locaux 0..3.
    grid = [None] * ((N + 1) * (N + 1))
    for cy in range(N):
        for cx in range(N):
            cell = cy * N + cx
            grid[idx(cx, cy)] = mesh.node(0, cell, 0)
            grid[idx(cx + 1, cy)] = mesh.node(0, cell, 1)
            grid[idx(cx + 1, cy + 1)] = mesh.node(0, cell, 2)
            grid[idx(cx, cy + 1)] = mesh.node(0, cell, 3)
    fes = pyrucast.FiniteElementSpace(mesh)

    left = [grid[idx(0, j)] for j in range(N + 1)]
    bottom = [grid[idx(i, 0)] for i in range(N + 1)]
    model = pyrucast.model.elasticity(fes, "plane_stress")
    model = model | _clamp(model, left, "u_x")
    model = model | _clamp(model, bottom, "u_y")

    # Traction S on the right edge → consistent nodal loads (the flux op).
    right = pyrucast.Mesh(c, "SEG2")
    for j in range(N):
        right.unit().add_cell([grid[idx(N, j)], grid[idx(N, j + 1)]])
    right_fes = pyrucast.FiniteElementSpace(right)
    model = model | pyrucast.model.flux(right_fes, model, "f_x")
    materials = pyrucast.element_field.material_field(
        model, [("E", E), ("nu", NU), ("phi_f_x", S)]
    )
    rhs = pyrucast.node_field.external_forces(model, materials)

    solution = pyrucast.solver.solve(pyrucast.matrix.stiffness(model, materials), rhs)

    print(f"{'x':>5} {'y':>5} {'u_x':>12} {'u_y':>12}")
    tol = 1e-10
    for j in range(N + 1):
        for i in range(N + 1):
            x, y = i * h, j * h
            ux = solution.value(grid[idx(i, j)], "u_x")
            uy = solution.value(grid[idx(i, j)], "u_y")
            print(f"{x:5.2f} {y:5.2f} {ux:12.6e} {uy:12.6e}")
            assert abs(ux - S / E * x) < tol
            assert abs(uy + NU * S / E * y) < tol
    print("\nOK: uniaxial field matching u_x=(S/E)x, u_y=-(νS/E)y.")


if __name__ == "__main__":
    main()

Compléments

Thermomécanique non couplée

Première brique de thermomécanique : une température imposée ΔT engendre une déformation thermique de libre dilatation ε_th = α·(T − T_ref), d’où des contraintes mécaniques — sans rétroaction de la mécanique sur le thermique. En petites déformations, la rigidité K reste l’élastique ; le terme thermique n’agit que sur le second membre et sur la contrainte réelle :

\[ \sigma = D : (\varepsilon(u) - \varepsilon_{th}), \qquad f_{th} = \int_\Omega B^\top D\, \varepsilon_{th}\, d\Omega. \]

Aucune physique nouvelle : on compose les briques existantes. alpha est fourni au champ matériau (composante facultative) ; la température, portée aux points de Gauss par interp_to_gauss, alimente thermal_strain (EPTH) ; la charge thermique sort de integrate_behavior + internal_forces (BSIG) ; enfin la contrainte réelle se relit sur deformation(u) − ε_th.

Deux régimes sur une barre chauffée valident les fermetures analytiques : bord en x encastré aux deux bouts ⇒ σ_xx = −E·α·ΔT ; appuis simples ⇒ dilatation libre u = α·ΔT·(x, y) sans contrainte. Code = test tests/thermoelastic_bar.rs :

use pyrucast::aggregate::Aggregate;
use pyrucast::atoms::{ElementType, Node};
use pyrucast::containers::field::Field;
use pyrucast::containers::finite_element_space::FiniteElementSpace;
use pyrucast::containers::mesh::{Mesh, SubMesh};
use pyrucast::containers::model::Model;
use pyrucast::containers::node_field::{NodeField, SubNodeField};
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::tensor::Kinematics;
use pyrucast::ops::element_field::{deformation, interp_to_gauss, thermal_strain};
use pyrucast::ops::mesh;
use pyrucast::ops::model;
use pyrucast::ops::node_field::internal_forces;
use pyrucast::ops::solver::lu::solve;
use pyrucast::Result;

#[test]
fn thermoelastic_constrained_bar_stress() -> Result<()> {
    const E: f64 = 210_000.0;
    const NU: f64 = 0.3;
    const ALPHA: f64 = 1e-5;
    const T_REF: f64 = 20.0;
    const DT: f64 = 100.0;
    const NX: usize = 4;
    const NY: usize = 2;
    const L: f64 = 4.0;
    const H: f64 = 1.0;
    let (hx, hy) = (L / NX as f64, H / NY as f64);

    // ── QUA4 mesh on [0,L]×[0,H] ───────────────────────────────────────────
    let coords = Handle::new(Coords::new(2)?);
    let idx = |i: usize, j: usize| j * (NX + 1) + i;
    let mut grid: Vec<Node> = Vec::new();
    for j in 0..=NY {
        for i in 0..=NX {
            grid.push(Node::create_in(
                coords.clone(),
                &[i as f64 * hx, j as f64 * hy],
            )?);
        }
    }
    let mut mesh = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::QUA4));
    for j in 0..NY {
        for i in 0..NX {
            mesh.add_cell(&[
                grid[idx(i, j)].id(),
                grid[idx(i + 1, j)].id(),
                grid[idx(i + 1, j + 1)].id(),
                grid[idx(i, j + 1)].id(),
            ])?;
        }
    }
    let fes = FiniteElementSpace::lagrange1(&mesh)?;

    // ── Modèle : élasticité + deux bords en x encastrés + appui u_y en bas ──
    let clamp = |target: &Model, nodes: &[Node], var: &str| -> Result<Model> {
        let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(nodes)?);
        let multiplier = mesh::barycenter(&imposed)?;
        model::dirichlet(target, var, &imposed, &multiplier, Default::default())
    };
    let left: Vec<Node> = (0..=NY).map(|j| grid[idx(0, j)].clone()).collect();
    let right: Vec<Node> = (0..=NY).map(|j| grid[idx(NX, j)].clone()).collect();
    let bottom: Vec<Node> = (0..=NX).map(|i| grid[idx(i, 0)].clone()).collect();
    let mut model = model::elasticity(&fes, Kinematics::PlaneStress)?;
    model = model.union(&clamp(&model, &left, "u_x")?)?;
    model = model.union(&clamp(&model, &right, "u_x")?)?;
    model = model.union(&clamp(&model, &bottom, "u_y")?)?;

    // `alpha` supplied through the material field — an optional elastic component.
    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[("E", E), ("nu", NU), ("alpha", ALPHA)],
    )?;

    // ── Imposed temperature T = T_ref + ΔT everywhere, carried at the Gauss points
    let support = Handle::new(SubMesh::poi1_from_nodes(&grid)?);
    let mut t_nodal = SubNodeField::from_poi1(&support, vec!["T".into()])?;
    for n in &grid {
        t_nodal.set_value(n.id(), "T", T_REF + DT)?;
    }
    let t_elem = interp_to_gauss(&NodeField::from_sub(t_nodal), &fes)?;

    // ── Charge thermique f_th = ∫ Bᵀ D ε_th (BSIG de σ_th = D:ε_th) ─────────
    let eps_th = thermal_strain(&t_elem, &materials, &fes, T_REF)?;
    let sig_th =
        pyrucast::ops::element_field::behavior::integrate(&model, &eps_th, None, &materials, None)?;
    let f_th = internal_forces(&model, &sig_th, &NodeField::empty(), &materials)?;

    // ── Assemblage + résolution ────────────────────────────────────────────
    let solution = solve(
        &pyrucast::ops::matrix::stiffness(&model, &materials)?,
        &f_th,
    )?;

    // ── Déplacement propre (u_x, u_y) puis σ = D:(ε(u) − ε_th) ─────────────
    let disp_support = Handle::new(SubMesh::poi1_from_nodes(&grid)?);
    let mut disp = SubNodeField::from_poi1(&disp_support, vec!["u_x".into(), "u_y".into()])?;
    for n in &grid {
        disp.set_value(n.id(), "u_x", solution.value(n.id(), "u_x")?)?;
        disp.set_value(n.id(), "u_y", solution.value(n.id(), "u_y")?)?;
    }
    let eps = deformation(&NodeField::from_sub(disp), &fes)?;
    let eps_mech = eps.merge_field(&eps_th, |a, b| a - b)?;
    let sigma = pyrucast::ops::element_field::behavior::integrate(
        &model, &eps_mech, None, &materials, None,
    )?;

    // ── Vérification : σ_xx = −E·α·ΔT, σ_yy = 0 ────────────────────────────
    let expected = -E * ALPHA * DT;
    let tol = 1e-6 * expected.abs();
    let sub = sigma.get(0)?.read();
    for cell in 0..sub.cell_count() {
        for g in 0..sub.gauss_count() {
            assert!((sub.value(cell, g, "sigma_xx")? - expected).abs() < tol);
            assert!(sub.value(cell, g, "sigma_yy")?.abs() < tol);
        }
    }
    Ok(())
}

"""Thermomécanique non couplée — barre chauffée (contraintes planes).

Physique
--------
An imposed temperature ΔT generates a free thermal strain
`ε_th = α·(T − T_ref)`. En petites déformations la rigidité reste élastique ;
the thermal term acts only on the right-hand side (equivalent thermal load
`f_th = ∫ Bᵀ D ε_th`) and on the real stress
`σ = D:(ε(u) − ε_th)`. Uncoupled: the mechanics does not feed back on the thermics.

Bricks composed by hand (no "all-in-one" operator):
`interp_to_gauss` (T nodale → points de Gauss), `thermal_strain` (EPTH),
`integrate_behavior` (σ = D:ε), `internal_forces` (BSIG), `solve`, puis
`deformation` and a field subtraction for ε_mech = ε(u) − ε_th.

`alpha` travels through the material field (`material_field`) as an
**optional** component of the elasticity, beside `E`/`nu`.

Two regimes on the same bar
------------------------------
- **bloquée** (deux bords en x encastrés) : `σ_xx = −E·α·ΔT`, `σ_yy = 0` ;
- **libre** (appuis simples) : dilatation `u = α·ΔT·(x, y)`, `σ ≈ 0`.

Lancement ::

    maturin develop --features extension-module
    python examples/thermoelastique_barre.py
"""

import pyrucast

E, NU, ALPHA = 210_000.0, 0.3, 1e-5
T_REF, DT = 20.0, 100.0
NX, NY, L, H = 4, 2, 4.0, 1.0


def _clamp(target, nodes, var):
    imposed = pyrucast.mesh.poi1_from_nodes(nodes)
    multiplier = pyrucast.mesh.barycenter(imposed)
    return pyrucast.model.dirichlet(target, var, imposed, multiplier)


def _bar():
    """An NX×NY grid of QUA4 on [0,L]×[0,H]. Returns (coords, grid, fes, idx)."""
    c = pyrucast.Coords(2)
    hx, hy = L / NX, H / NY

    def idx(i, j):
        return j * (NX + 1) + i

    grid = [c.add_node([i * hx, j * hy]) for j in range(NY + 1) for i in range(NX + 1)]
    mesh = pyrucast.Mesh(c, "QUA4")
    for j in range(NY):
        for i in range(NX):
            mesh.unit().add_cell(
                [
                    grid[idx(i, j)],
                    grid[idx(i + 1, j)],
                    grid[idx(i + 1, j + 1)],
                    grid[idx(i, j + 1)],
                ]
            )
    return c, grid, pyrucast.FiniteElementSpace(mesh), idx


def _uniform_temperature(c, grid, fes, value):
    """A temperature field 'T' = value everywhere, carried at the Gauss points."""
    t_mesh = pyrucast.Mesh(c, "POI1")
    for node in grid:
        t_mesh.unit().add_cell([node])
    t_nodal = pyrucast.NodeField(t_mesh, ["T"])
    for node in grid:
        t_nodal[0].set_value(node, "T", value)
    return pyrucast.element_field.interp_to_gauss(t_nodal, fes)


def _displacement(solution, c, grid):
    """Extracts a clean (u_x, u_y) field (without the Lagrange multipliers)."""
    u_mesh = pyrucast.Mesh(c, "POI1")
    for node in grid:
        u_mesh.unit().add_cell([node])
    u = pyrucast.NodeField(u_mesh, ["u_x", "u_y"])
    for node in grid:
        u[0].set_value(node, "u_x", solution.value(node, "u_x"))
        u[0].set_value(node, "u_y", solution.value(node, "u_y"))
    return u


def _solve_thermal(model, materials, fes, c, grid):
    """ε_th → charge thermique → u → σ = D:(ε(u) − ε_th)."""
    eps_th = pyrucast.element_field.thermal_strain(
        _uniform_temperature(c, grid, fes, T_REF + DT), materials, fes, T_REF
    )
    sig_th = pyrucast.element_field.integrate_behavior(model, eps_th, materials)
    # A **load**, not a residual: the divergence of the prescribed tensor, hence
    # the geometric operator. They are then renamed into dual rows — that is
    # where, and only where, these numbers become forces.
    f_th = (
        pyrucast.node_field.divergence(sig_th, "sigma")
        .rename_component("div_sigma_x", "f_x")
        .rename_component("div_sigma_y", "f_y")
    )
    solution = pyrucast.solver.solve(pyrucast.matrix.stiffness(model, materials), f_th)
    u = _displacement(solution, c, grid)
    sigma = pyrucast.element_field.integrate_behavior(
        model, pyrucast.element_field.deformation(u, fes) - eps_th, materials
    )
    return u, sigma


def main() -> None:
    # ── Régime bloqué : σ_xx = −E·α·ΔT ──────────────────────────────────────
    c, grid, fes, idx = _bar()
    left = [grid[idx(0, j)] for j in range(NY + 1)]
    right = [grid[idx(NX, j)] for j in range(NY + 1)]
    bottom = [grid[idx(i, 0)] for i in range(NX + 1)]

    model = pyrucast.model.elasticity(fes, "plane_stress")
    model = (
        model
        | _clamp(model, left, "u_x")
        | _clamp(model, right, "u_x")
        | _clamp(model, bottom, "u_y")
    )
    materials = pyrucast.element_field.material_field(
        model, [("E", E), ("nu", NU), ("alpha", ALPHA)]
    )

    _u, sigma = _solve_thermal(model, materials, fes, c, grid)
    expected = -E * ALPHA * DT
    sub = sigma[0]
    sxx = sub.value(0, 0, "sigma_xx")
    print(f"Bloquée : σ_xx = {sxx:12.4f}  (attendu {expected:.4f} = −E·α·ΔT)")
    assert abs(sxx - expected) < 1e-6 * abs(expected)

    # ── Régime libre : dilatation u = α·ΔT·(x, y), σ ≈ 0 ────────────────────
    c, grid, fes, idx = _bar()
    left = [grid[idx(0, j)] for j in range(NY + 1)]
    bottom = [grid[idx(i, 0)] for i in range(NX + 1)]

    model = pyrucast.model.elasticity(fes, "plane_stress")
    model = model | _clamp(model, left, "u_x") | _clamp(model, bottom, "u_y")
    materials = pyrucast.element_field.material_field(
        model, [("E", E), ("nu", NU), ("alpha", ALPHA)]
    )

    u, sigma = _solve_thermal(model, materials, fes, c, grid)
    tip = grid[idx(NX, NY)]
    ux, uy = u.value(tip, "u_x"), u.value(tip, "u_y")
    print(
        f"Libre   : u(coin) = ({ux:.6e}, {uy:.6e})  (attendu ({ALPHA * DT * L:.6e}, {ALPHA * DT * H:.6e}))"
    )
    assert abs(ux - ALPHA * DT * L) < 1e-9 and abs(uy - ALPHA * DT * H) < 1e-9
    assert abs(sigma[0].value(0, 0, "sigma_xx")) < 1e-6

    print("\nOK: blocked bar → σ_xx = −E·α·ΔT; free bar → expansion without stress.")


if __name__ == "__main__":
    main()

Élasticité orthotrope et anisotrope

Introduction

L’élasticité linéaire suppose le matériau isotrope : deux constantes, aucune direction privilégiée. Beaucoup de matériaux n’obéissent pas à cette hypothèse — un composite tissé, un bois, une tôle laminée, un monocristal ont des raideurs différentes selon la direction.

pyrucast traite cela comme un axe à part entière, la symétrie matériau, porté par le même modèle elasticity et orthogonal à l’hypothèse cinématique (contraintes planes, déformations planes, axisymétrie, massif) :

C’est la convention de Cast3M, où ISOTROPE / ORTHOTROPE / ANISOTROPE qualifie le matériau d’une formulation plutôt que de nommer un modèle différent. Les degrés de liberté ne changent donc pas : déplacement u_x, u_y(, u_z) en primal, force nodale f_x, … en dual, exactement comme en isotrope. Seule la matrice de Hooke D change.

Équations continues résolues

Les mêmes qu’en élasticité linéaire — équilibre ∇·σ + f = 0 en petites déformations, avec ε = ½(∇u + ∇uᵀ). C’est la loi de comportement qui se généralise :

\[ \sigma_{ij} = C_{ijkl}\,\varepsilon_{kl} \]

où \( C \) est le tenseur d’élasticité d’ordre 4. Ses symétries — mineures \( C_{ijkl} = C_{jikl} = C_{ijlk} \), qui viennent de celles de \( \sigma \) et \( \varepsilon \), et majeure \( C_{ijkl} = C_{klij} \), qui vient de l’existence d’un potentiel élastique \( W = \tfrac12\,\varepsilon : C : \varepsilon \) — le réduisent d’un tenseur à 81 composantes à une matrice \( 6\times6 \) symétrique en notation de Voigt, soit 21 constantes dans le cas général.

L’orthotropie est le cas où le matériau possède trois plans de symétrie orthogonaux. Dans ses axes propres, la souplesse S = C⁻¹ se découple : les termes normaux ne sont couplés qu’entre eux, et chaque cisaillement est isolé.

\[ S = \begin{bmatrix} 1/E_1 & -\nu_{21}/E_2 & -\nu_{31}/E_3 & & & \\ -\nu_{12}/E_1 & 1/E_2 & -\nu_{32}/E_3 & & & \\ -\nu_{13}/E_1 & -\nu_{23}/E_2 & 1/E_3 & & & \\ & & & 1/G_{23} & & \\ & & & & 1/G_{13} & \\ & & & & & 1/G_{12} \end{bmatrix} \]

Les relations de réciprocité \( \nu_{ji}/E_j = \nu_{ij}/E_i \) rendent la matrice symétrique, d’où neuf constantes seulement : trois modules d’Young, trois coefficients de Poisson, trois modules de cisaillement. Le bloc normal et les trois cisaillements sont découplés, ce qui est la définition même de l’orthotropie : une traction selon un axe propre ne produit aucun cisaillement.

Attention — un jeu de constantes n’est pas physique par construction : S doit rester définie positive, ce qui impose ν_ij² < E_i/E_j. pyrucast le vérifie en inversant S et erronne si elle est singulière, plutôt que d’assembler une raideur non définie positive en silence.

Le repère d’orthotropie

Les constantes sont données dans les axes matériau, qui ne coïncident pas avec les axes globaux. Il faut donc dire où ils pointent.

pyrucast suit Cast3M et les décrit par des vecteurs, pas par des angles d’Euler (MATE 'DIRECTION' V1 V2). Ils voyagent dans le champ matériau comme n’importe quel autre coefficient :

Ils sont orthonormalisés en interne (Gram-Schmidt) : V2 n’a besoin d’être que grossièrement perpendiculaire à V1, c’est le plan qu’ils engendrent qui compte. Des vecteurs plutôt que des angles, parce qu’il n’y a aucune convention à retenir, aucun cas de blocage de cardan — et surtout parce que le repère varie alors naturellement d’une maille à l’autre (un composite bobiné, une pièce courbe), en passant par le canal matériau existant.

Un V1 nul, ou un V2 parallèle à V1, est un repère dégénéré : il est signalé, jamais complété arbitrairement.

Forme discrétisée

La chaîne est celle de l’élasticité — K_e = Σ_g Bᵀ D B |J| w, avec le même opérateur B. Ce qui change tient en trois étapes, faites une fois par maille :

1. construire D dans les axes matériau, où l’orthotropie est diagonale ;
2. le tourner vers les axes globaux ;
3. le réduire au modèle cinématique (bloc [xx, yy, xy] en déformations planes, sa condensation statique sur ε_zz en contraintes planes, le bloc [rr, zz, θθ, rz] en axisymétrie, le 6×6 complet en massif).