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)
   │
   └── Mesh  (agrège des SubMesh, + ElementType)   géométrie
          ├── NodeField  (agrège des SubNodeField, support POI1)   valeurs par nœud × composante
          └── 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)

Evolution  (agrège des SubEvolution)    scalaire, NodeField ou ElementField
                                        tabulé vs une variable

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 — | en Python, union en Rust) 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 mesh = mesh::line(&a, &b, 10, ElementType::SEG2)?;
    println!("{}", mesh); // Mesh: 1 submesh(es) [#…], 10 cell(s) total
    Ok(())
}

En Python, il suffit d’installer le paquet :

pip install pyrucast

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.line(a, b, 10)
print(mesh)  # Mesh: 1 submesh(es) [#…], 10 cell(s) total

Installation et démarrage rapide

Cette page suffit pour utiliser pyrucast en Python : l’installer, vérifier qu’il fonctionne, écrire un premier script. Deux voies : la version précompilée, sans rien compiler, ou la compilation depuis les sources. 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.

Version précompilée

pip install pyrucast

et pour passer à la dernière version :

pip install --upgrade pyrucast

PyPI distribue pyrucast en wheels — des paquets précompilés, prêts à installer — pour trois plateformes :

Seul prérequis : Python ≥ 3.11. Ni Rust ni compilateur. La wheel comprend la visualisation (plot), fenêtre interactive incluse.

Sur une autre plateforme (Linux ARM, musl, BSD…), pip ne trouve pas de wheel et compile lui-même les sources : il faut alors Rust (voir ci-dessous). Cette compilation automatique se fait sans la visualisation — voir Ce que porte chaque distribution. Pour l’avoir, compilez à la main, comme le décrit la partie Compilation ci-dessous.

Compilation

Pour suivre la version en développement, ou pour compiler avec la visualisation sur une plateforme sans wheel.

Prérequis

• Git, pour cloner le dépôt.
• Rust stable, installé via rustup.
• Python ≥ 3.11.
• Linux uniquement : la visualisation lie fontconfig et freetype, il faut leurs en-têtes — pkg-config libfontconfig1-dev libfreetype6-dev (Debian/Ubuntu) ou pkg-config fontconfig-devel freetype-devel (Fedora/RHEL). Windows et macOS passent par les polices du système.

Compiler et installer

On clone le dépôt, puis on compile à sa racine, dans un environnement virtuel : maturin develop refuse de s’en passer (« Couldn’t find a virtualenv or conda environment »).

Linux / macOS (bash)

git clone https://github.com/Pyrucast/pyrucast
cd pyrucast
python3 -m venv .venv
source .venv/bin/activate
pip install --upgrade pip maturin
maturin develop --release --features extension-module,viz-interactive

Windows (PowerShell)

git clone https://github.com/Pyrucast/pyrucast
cd pyrucast
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). --release compile en optimisé, ce qu’il faut pour tout usage réel : un build de débogage est bien plus lent à l’exécution. Les features font le reste : extension-module fait du binaire un module Python, et viz-interactive ajoute la visualisation, fenêtre interactive comprise (elle inclut viz, l’export PNG/SVG).

Après toute modification du Rust, relancer la même commande : sans --features …viz-interactive, le module recompilé perdrait la visualisation.

Vérifier l’installation

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>

pyrucast.__features__ dit ce que porte le module installé : viz y figure quand la visualisation est compilée.

Premier script

Une ligne de 10 éléments entre deux nœuds, par la fonction libre line du module pyrucast.mesh :

import pyrucast

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

ligne = pyrucast.mesh.line(a, b, 10)
print(ligne)  # Mesh: 1 submesh(es) [#…], 10 cell(s) total

Trois niveaux de listing

Tout objet pyrucast se liste à trois niveaux : print(x) donne le résumé d’une ligne, repr(x) la structure, x.dump() le contenu complet — les listings résumé et complet de Cast3M. Le détail, avec un exemple, est dans Aspect informatique.

Fonctions libres et méthodes

Les opérateurs de pyrucast sont des fonctions libres, rangées par ce qu’elles produisent : tout ce qui rend un maillage est dans pyrucast.mesh. Quand une fonction prend en premier l’objet qu’elle transforme et rend un nouvel objet, elle existe en règle générale aussi en méthode de cet objet — les deux écritures sont équivalentes :

plaque = pyrucast.mesh.extrude(ligne, [0.0, 1.0], 5)  # fonction libre
plaque = ligne.extrude([0.0, 1.0], 5)  # la même, en méthode

La méthode sert à composer. Pour la bordure d’une plaque tirée de la ligne, les fonctions libres s’imbriquent et se lisent de l’intérieur ; les méthodes s’enchaînent et se lisent dans l’ordre des opérations :

bord = pyrucast.mesh.border(pyrucast.mesh.extrude(ligne, [0.0, 1.0], 5))  # imbriqué
bord = ligne.extrude([0.0, 1.0], 5).border()  # enchaîné
print(bord)  # Mesh: 1 submesh(es) [#…], 30 cell(s) total

La règle complète — quelles fonctions ont leur méthode, et pourquoi — est dans Conventions.

Voir le maillage

Une fenêtre interactive s’ouvre, s’il y a un écran :

import os

# 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")
plaque.plot(save=None if ecran else "apercu.svg")

Aller plus loin

À 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.

• Tutoriel : un parcours guidé, du maillage au calcul thermique et mécanique — Premiers pas.
• Référence API Rust : publiée sur https://pyrucast.github.io/pyrucast/rust/, ou en local avec 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, vue de Python : comment les objets vivent, et les quelques motifs qu’on retrouve partout. Quand cela éclaire l’usage, un encart dit ce qu’il y a dessous, côté Rust. 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. Les classes Python (pyrucast.Coords, pyrucast.Mesh, …) enveloppent des types Rust, et les fonctions de module (pyrucast.matrix.stiffness, pyrucast.solver.solve, …) appellent des fonctions 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. Seules quelques parties sont écrites en Python, parce qu’elles n’ont rien à faire en Rust : les échanges avec gmsh et MED (mesh.from_gmsh, export.to_medcoupling…), qui appellent ces modules Python, et la couche pas à pas (pyrucast.stepping), qui orchestre les opérateurs. 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

Durée de vie des objets

En Python, rien à libérer à la main : un objet pyrucast vit tant qu’une variable ou un autre objet le référence — un maillage garde ses nœuds, un champ garde son support.

Les nœuds ont en plus leur propre compteur, dans les Coords : un nœud reste vivant tant qu’un maillage ou un champ l’utilise, même quand plus aucune variable ne le désigne. coords.gc() collecte ceux qui ne servent plus. Voir Coordonnées et Nœud.

Côté Rust, chaque objet vit derrière un Handle<T> — une référence comptée munie de son propre verrou — et la Coords porte un compteur par nœud. Le détail est dans Modèle mémoire.

Agrégats et zones

La plupart des conteneurs viennent par paires : une zone (Sub…) et son agrégat, une liste de zones qui s’utilise toujours de la même façon :

Un agrégat s’utilise comme une liste de zones : len(mesh), mesh[i], for zone in mesh, et mesh.unit() quand il n’en a qu’une. Les zones ne se construisent pas directement (sauf SubEvolution, une courbe tabulée) : on construit l’agrégat, et mesh[i] donne une vue sur une zone.

Deux agrégats se composent par union : mesh_a | mesh_b. Les zones y sont partagées, pas copiées, et une zone déjà présente n’est pas ajoutée deux fois. L’opérateur +, lui, additionne des valeurs — champs ou matrices — et ne compose jamais (cf. Champ).

Côté Rust, ce motif est factorisé dans le trait Aggregate, et l’union s’y écrit a.union(&b) — voir Agrégat.

Trois niveaux de listing

Tout objet pyrucast se liste à trois niveaux, du plus court au plus complet :

• print(x) ou str(x) — le résumé, une ligne : le listing résumé de Cast3M ;
• repr(x) — la structure : compteurs, dimensions, noms, identifiants, sans jamais le contenu ;
• x.dump() — le contenu : valeurs, coordonnées, connectivité, comme le listing complet de Cast3M. Il reste borné : au-delà de max_rows lignes (20 par défaut), la table est coupée et la coupure annoncée.

Sur des coordonnées qui portent les 11 nœuds d’une ligne de 10 éléments :

import pyrucast

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

print(c)  # résumé : Coords: dim=2, configs=1 (active="default"), nodes=11 …
print(repr(c))  # structure : Coords { dim: 2, frame: Cartesian, … }
c.dump()  # contenu : la table des 11 nœuds, coordonnées et références

Seule exception, Array — le tableau que rendent les fonctions d’échange — n’a pas de dump() : simple tampon sans homologue Rust, ce n’est pas un objet pyrucast à proprement parler.

Côté Rust, les trois niveaux sont les traits Display, Debug et Dump.

Erreurs

En Python, une erreur signalée par pyrucast lève une RuntimeError, dont le message dit ce qui ne va pas. Les cas que Python nomme déjà gardent leur exception : IndexError pour un indice hors bornes, ZeroDivisionError pour une matrice divisée par zéro, KeyboardInterrupt sur un Ctrl+C.

Côté Rust, toute l’API publique renvoie Result<T, PyrucastError> : un seul type d’erreur pour toute la librairie — voir Conventions.

Persistance portable

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. Le fichier est identique Linux ↔ Windows. Voir Sauvegarde et relecture.

Côté Rust, ce format repose sur un trait unique, Portable (serde + bincode) — voir Conventions et Modèle mémoire.

Pour le développeur Rust

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

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                            le magasin de nœuds, et l'un d'eux
   │
   └── Mesh                               agrège des SubMesh
          ├── NodeField                   champ aux nœuds, sur un support POI1
          └── FiniteElementSpace          agrège des SubFiniteElementSpace
                 ├── ElementField         champ aux points de Gauss
                 └── Model                agrège des SubModel : physiques et contraintes
                        └── Matrix        matrice creuse, sortie d'assemblage

Evolution                                 agrège des SubEvolution : valeur tabulée selon une variable

Une distinction structure toute cette liste et se lit dans l’arborescence : Coords est le magasin (src/coords.rs), Node, Cell et Element 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 (durée de vie des objets, motif zone/agrégat), voir Aspect informatique ; pour ce qu’il y a dessous côté Rust (Handle<T>, compteurs de références), Modèle mémoire. Pour les écrire dans un fichier et les relire tels qu’ils étaient, Sauvegarde et relecture.

Coordonnées (Coords)

Un Coords est le magasin des nœuds : il leur donne une identité et des coordonnées, en dimension fixée, dans une ou plusieurs configurations (référence, déformée…). C’est le premier objet de pyrucast — maillages, champs et espaces éléments finis s’y rattachent. On y crée les nœuds, et chacun se manipule par un Node, l’objet du chapitre suivant.

import pyrucast

c = pyrucast.Coords(dim=2)
a = c.add_node([0.0, 0.0])  # un pyrucast.Node
b = c.add_node([1.0, 0.0])
print(a.id, b.id)  # 0 1 — identifiants stables
print(c)  # Coords: dim=2, configs=1 (active="default"), nodes=2 …

ligne = pyrucast.mesh.line(a, b, 2)  # pose un nœud au milieu, id 2
print(c.refcount(a.id))  # 2 — la variable a, et la maille qui le contient
print(c.refcount(2))  # 2 — les deux mailles qui le contiennent
milieu = c.acquire(2)  # un Node pour un id existant
print(milieu.position())  # [0.5, 0.0]

Identité et durée de vie des nœuds

Chaque nœud créé reçoit un identifiant (node.id) stable et jamais réutilisé, valable pour toute la vie du Coords. C’est par lui que les maillages et les champs désignent leurs nœuds.

Un nœud ne se supprime pas : il n’existe aucune méthode pour le retirer, et ses coordonnées restent lisibles tant que le Coords existe. refcount(id) dit combien de références le désignent — chaque Node qui le tient, chaque maille qui le contient — et acquire(id) rend un nouveau Node pour un identifiant existant. Ces deux méthodes prennent un identifiant brut, pas un Node : passer un Node ajouterait la référence même qu’on cherche à compter.

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.

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é par add_node comme par 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 de plot la balaie pour dessiner le corps de révolution lui-même — voir Visualisation.

Plusieurs configurations

Une configuration est un jeu complet de coordonnées pour les mêmes nœuds : référence, déformée, prédite. add_config(nom) clone la configuration active sous un nouveau nom et rend son indice ; select(i) change la configuration active, et node.position() lit toujours dans celle-là.

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

Pour calculer dans une autre configuration — une matrice sur la déformée, une normale sur la paroi telle qu’elle est —, with coords.in_config(i): l’active le temps du bloc, puis rend la configuration d’avant, même si le bloc lève une exception. C’est la forme à préférer à une paire select(i) … select(0), qui laisse la mauvaise configuration active à la première exception. Les blocs s’imbriquent. La configuration active appartient au Coords, pas au bloc : un autre fil qui calcule sur le même Coords pendant le bloc la voit aussi.

import pyrucast

c = pyrucast.Coords(dim=2)
n = c.add_node([1.0, 0.0])
deformee = c.add_config("deformed")

# Le bloc voit `deformee` ; la configuration d'avant revient à la sortie,
# même après une exception.
with c.in_config(deformee):
    n.set_position([1.5, 0.0])
    x = n.position()

print(x)  # [1.5, 0.0]
print(n.position())  # [1.0, 0.0] — la référence, de nouveau active

Permutation pour le solveur

Une permutation optionnelle fixe l’ordre des nœuds pour le solveur sans toucher à leur identité : permutation[id] est la place du nœud id, et les matrices assemblées rangent leurs degrés de liberté dans cet ordre. Elle se pose à la main par set_permutation, avec une entrée par nœud créé (capacity) et aucune valeur en double ; clear_permutation revient à l’ordre des identifiants.

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])

# One entry per node created, no value twice.
print(c.capacity)  # 3
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

Inspecter un Coords

dim, is_axisymmetric, names() et active décrivent le Coords ; print(c) le résume en une ligne, repr(c) en donne la structure et c.dump() la table complète des nœuds, coordonnées et références — les trois niveaux de listing.

Pourquoi plusieurs configurations dans un seul Coords

L’autre conception possible serait un Coords par configuration. pyrucast garde un seul jeu d’identités par domaine, et plusieurs configurations dedans :

Des Coords distincts servent aux domaines vraiment indépendants.

Côté Rust

En Rust, un Coords se partage par un Handle<Coords>, une référence comptée munie de son verrou, et chaque nœud porte en plus son propre compteur dans le Coords : c’est lui que refcount lit, et que Node, les maillages et les champs incrémentent. Le détail est dans Modèle mémoire. Les mêmes opérations s’écrivent ainsi :

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();
    assert_eq!(axi.dim(), 2);
    assert!(axi.is_axisymmetric());
}

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

in_config y reçoit le bloc en fermeture, l’idiome Rust du with Python :

#[test]
fn calculer_dans_une_autre_configuration() -> pyrucast::Result<()> {
    use pyrucast::atoms::Node;
    use pyrucast::handle::Handle;
    use pyrucast::ops::coords as ops_coords;

    let coords = Handle::new(Coords::new(2)?);
    let n = Node::create_in(coords.clone(), &[1.0, 0.0])?;
    let deformee = coords.write().add_config("deformed");

    // Le bloc voit `deformee` ; la configuration d'avant revient à la sortie,
    // même si le bloc échoue ou panique.
    ops_coords::in_config(&coords, deformee, || n.set_position(&[1.5, 0.0]))?;
    let x = ops_coords::in_config(&coords, deformee, || n.position())?;

    assert_eq!(x, vec![1.5, 0.0]);
    assert_eq!(n.position()?, vec![1.0, 0.0]); // la référence, de nouveau active
    Ok(())
}

#[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: one entry per node created.
    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());
}

Nœud (Node)

Un Node désigne un nœud d’un Coords. Il ne stocke pas de coordonnées : il porte son Coords et l’identifiant du nœud, et lit tout dans le Coords — les coordonnées qu’il rend sont donc celles de la configuration active. C’est l’objet qu’on passe partout où une fonction attend un nœud : pyrucast.mesh.line(a, b, 10), pyrucast.mesh.arc(…), pyrucast.mesh.circle(…).

Créer et retrouver un nœud

Un Node naît toujours d’un Coords :

• coords.add_node([x, y, …]) crée le nœud et rend son Node ;
• coords.acquire(id) rend un Node pour un nœud qui existe déjà — un nœud posé par un mailleur, par exemple, dont on connaît l’identifiant.

Interface

• id — l’identifiant entier, stable dans son Coords ;
• position() — ses coordonnées dans la configuration active ;
• set_position([x, y, …]) — les réécrit, dans la configuration active ;
• coords() — le Coords auquel il appartient ;
• node | node — un Mesh POI1 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 d’une seule zone ;
• print(n), repr(n) et n.dump() — les trois niveaux de listing.

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), puis un point de plus.
poi = a | b
print(poi)  # Mesh: 1 submesh(es) [#…], 2 cell(s) total
poi = poi | c.add_node([2.0, 0.0])
print(poi.cell_count)  # 3

print(a)  # <Node #0>
a.dump()  # Node #0 @ (0.500, 0.500)

Un Node compte comme une référence

Chaque Node compte comme une référence à son nœud, celle que coords.refcount(id) additionne avec celles des maillages et des champs : add_node et acquire en ajoutent une, et le Node la rend quand Python le détruit. Une seconde variable sur le même objet (b = a) ne compte pas deux fois : c’est le même Node.

C’est ce qui le distingue d’un simple identifiant. Un entier ne dit pas à quel Coords il appartient et ne compte pour rien ; le Node porte son Coords, lit ses coordonnées dans la bonne configuration, et le nœud le compte parmi ses références tant qu’il existe. Porter son Coords lui permet aussi d’être refusé là où il n’a pas sa place : une maille, un champ, une matrice posés sur un autre Coords lèvent une erreur plutôt que de prendre le nœud de ce Coords qui porte le même numéro.

Côté Rust

En Rust, un Node est une paire (Handle<Coords>, NodeId) qui tient son compte par RAII : Clone ajoute une référence dans le Coords, Drop la rend. Il se crée par Node::create_in(coords, &[x, y, …]), et Node::acquire(coords, id) en rend un de plus pour un identifiant existant — le verbe que Python porte sur le Coords, Rust le porte sur le nœud. Il expose id(), position(), set_position(&[…]) et coords(), qui rend le Handle<Coords> porté.

#[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 id = n.id();
    let m = n.clone(); // une référence de plus
    assert_eq!(coords.read().refcount(id), 2);
    drop(n);
    drop(m); // chaque Drop rend la sienne
    assert_eq!(coords.read().refcount(id), 0);
    Ok(())
}

Le code interne peut manipuler des NodeId sans passer par Node, à condition de tenir lui-même le compte avec Coords::incref / Coords::decref : c’est ce que font les maillages et les champs. Le détail est dans Modèle mémoire.

Agrégat (Aggregate)

La plupart des conteneurs de pyrucast viennent par paires : une zone (SubMesh, SubFiniteElementSpace, SubNodeField, SubElementField, SubModel, SubMatrix, SubEvolution) et son agrégat (Mesh, FiniteElementSpace, NodeField, ElementField, Model, Matrix, Evolution), une liste de zones. Les sept agrégats s’utilisent exactement de la même façon et se composent de la même manière, par l’union |. Ce chapitre décrit ce contrat commun une fois pour toutes ; chaque chapitre d’objet n’en redonne que les spécificités.

Un agrégat ne copie jamais ses zones : il les partage. Deux agrégats peuvent porter la même zone, et la modifier par l’un se voit par l’autre.

Interface commune

agg[i] est une vue, agg[i:j] un agrégat. L’indice entier rend la zone elle-même ; le slicing rend un agrégat du même type, qui partage les zones retenues.

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

unit() plutôt que [0]. agg[0] prend sans rien dire la première de plusieurs zones ; agg.unit() exige qu’il n’y en ait qu’une et lève sinon. Là où une seule zone a un sens — ajouter une maille à un maillage qu’on vient de créer —, unit() dit ce qu’on suppose.

Les zones ne se construisent pas seules

Une zone s’obtient par agg[i] ; on ne la construit pas directement. On construit l’agrégat — Mesh(coords, type), FiniteElementSpace(mesh), ElementField(fes, composantes) — ou on le reçoit d’un opérateur (pyrucast.model.heat_conduction(fes)…), puis on indexe pour atteindre une zone. Seule exception, SubEvolution : une courbe tabulée n’a pas de parent qui la définirait, elle se construit directement (cf. Évolution).

Composition : l’union |

Composer deux agrégats, c’est l’union, a | b. Elle rend un agrégat neuf et laisse ses deux opérandes intacts, avec la même règle pour les sept agrégats :

1. Pas de doublon. Une zone déjà présente n’est pas ajoutée une seconde fois ; les zones gardent l’ordre de leur première apparition.
2. Partage, pas copie. Les zones retenues sont celles des opérandes, pas des copies.
3. Cohérence vérifiée. Ce qui doit être commun l’est, à l’ajout : par exemple, deux maillages sur des Coords différents se refusent (« mismatched Coords »).
4. Fusion des champs. Pour NodeField et ElementField, les zones posées sur le même support sont fusionnées en fin d’union (voir Champ, Champ aux nœuds).

L’union accepte agrégats et zones, dans les deux sens :

En place : add_sub / add_subs. Pour modifier un agrégat plutôt que d’en créer un, agg.add_sub(zone) ajoute une zone et agg.add_subs(other) toutes celles d’un autre agrégat, dans l’ordre. C’est de la concaténation : la cohérence est vérifiée, mais rien n’est dédoublonné ni fusionné. Pour composer, | reste la forme à préférer.

Les nœuds s’unissent aussi (cf. Nœud) :

| compose, + calcule. L’union s’écrit toujours |. L’opérateur + (et -, *, /) est l’arithmétique — des champs (champ + 2.0, champ + champ) et des matrices (k + m) — et ne compose jamais : un Mesh + Mesh lève un TypeError.

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

Côté Rust

Un agrégat est une liste de Handle<Sub> : il partage ses zones par le compte de références de leur Handle, ce qui explique qu’il ne les copie jamais.

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

Tout le contrat est porté par le trait Aggregate (src/aggregate.rs), qui le dérive de deux accesseurs, items() et items_mut() : len, is_empty, get(i), subset(indices) (sur lequel s’appuie le slicing Python), iter, unit, push, add_sub / add_subs, et l’union sous des noms explicites — a.union(&b), a.union_sub(&h) (zone en queue), a.union_sub_first(&h) (zone en tête), T::union_subs(&a, &b). « Pas de doublon » y veut dire : pas deux Handle sur le même objet (Handle::same_object). Deux points de personnalisation seulement : check_push, qui vérifie la cohérence à chaque ajout, et finalize, qui tourne en fin d’union et que les champs utilisent pour fusionner leurs zones.

C’est ce qui garantit qu’un sous-maillage, un sous-espace EF ou un sous-modèle s’indexent, se parcourent et se composent strictement de la même façon : aucun accès n’est réécrit par type. Un agrégat se câble par des macros : impl_aggregate!, qui pose aussi Index<usize> et IntoIterator, le plus souvent impl_aggregate_dump! pour le listing, et côté Python impl_aggregate_pymethods!.

Maillage (Mesh/SubMesh)

Un maillage est un agrégat de zones :

• une zone (SubMesh) regroupe des mailles d’un seul type d’élément : que des TRI3, ou que des QUA4… Chaque maille y est une liste de nœuds d’un Coords ;
• un maillage (Mesh) est une liste de zones, toutes sur le même Coords. C’est ce qui lui permet de mêler triangles et quadrangles : une zone par type.

Le type POI1 est le cas dégénéré : une maille, un nœud. Une zone POI1 est donc exactement une liste de nœuds, et c’est elle qui sert de support aux champs aux nœuds.

Cette page décrit l’objet maillage. Pour construire un maillage, on passe en général par les opérateurs de maillage (line, extrude, pave_surface…, cf. Mailler une géométrie) ; la construction à la main, maille par maille, sert surtout aux petits cas et aux tests.

Créer et lire un maillage

pyrucast.Mesh(coords, "TRI3") crée un maillage d’une zone vide de ce type (pyrucast.Mesh(coords), sans type, un maillage sans zone). mesh.unit() atteint cette zone, et zone.add_cell([n0, n1, n2]) y pose une maille et rend son indice. Les nœuds se donnent dans l’ordre local de l’élément, celui que décrit sa fiche, section « Élément de référence » ; un nombre de nœuds qui ne correspond pas au type est refusé, comme un nœud d’un autre Coords. Plusieurs types se composent par l’union |.

import pyrucast

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

# Mesh(coords, type) : un maillage d'une zone vide ; unit() l'atteint,
# add_cell y pose une maille (nœuds dans l'ordre local de l'élément) et rend
# son indice.
quad = pyrucast.Mesh(c, "QUA4")
quad.unit().add_cell([a, b, n3, n4])  # 0
tri = pyrucast.Mesh(c, "TRI3")
tri.unit().add_cell([b, n5, n3])

# Deux types d'éléments, deux zones, composées par l'union.
mesh = quad | tri
print(mesh.element_types())  # ['QUA4', 'TRI3']
print(mesh.cell_counts())  # [1, 1]
print(mesh.cell_count)  # 2

# Lire : zone i, maille j, nœud k.
cell = mesh[1][0]  # la même vue que mesh.cell(1, 0)
print(len(cell), cell.element_type)  # 3 TRI3
print(cell[1].position())  # [2.0, 0.5]
print(mesh.node(0, 0, 2).id)  # 2
print(mesh.nearest_node([1.9, 0.4]).id)  # 4

Une maille (Cell) est une vue : elle tient sa zone et son indice, et lit ses nœuds à la demande. len(cell) donne son nombre de nœuds, cell[k] le k-ième Node, cell.nodes() la liste ; cell.index et cell.element_type la situent. cell.dump() liste ses nœuds avec leurs coordonnées.

Les opérateurs de maillage s’appellent aussi en méthode : mesh.skin() est pyrucast.mesh.skin(mesh). Maillage, zone et maille s’affichent sur les trois niveaux de listing ; le dump() d’une zone donne sa table de connectivité.

Types d’éléments

Le type d’une zone se désigne par son nom.

Les huit derniers types sont quadratiques. Ils reprennent les sommets de leur parent linéaire, dans le même ordre, puis ajoutent les nœuds de milieu d’arête, dans l’ordre d’arêtes de la convention VTK ; QUA9 et HEX27 ajoutent encore leurs nœuds de face et leur nœud central. QUA8, PENTA15 et HEX20 sont sérendipité (nœuds d’arête seulement), les autres sont des Lagrange complets. Un espace éléments finis les accepte avec l’interpolation LAGRANGE2, qu’il prend de lui-même : pyrucast.FiniteElementSpace(mesh) choisit pour chaque zone le Lagrange de son degré (cf. Espace éléments finis).

Un maillage compte ses nœuds

Chaque fois qu’un nœud apparaît dans une maille, la zone prend une référence sur lui dans le Coords : un nœud partagé par quatre mailles en reçoit quatre, que coords.refcount(id) additionne avec celles des Node et des autres objets (cf. Coordonnées). La zone les rend toutes à sa destruction, quand plus rien ne la tient — ni variable Python, ni maillage, ni espace éléments finis, ni champ. Un add_cell refusé ne prend aucune référence.

Scellement

Un espace éléments finis, un champ ou une matrice indexent les mailles ou les nœuds de la zone sur laquelle ils sont posés. Ajouter une maille ensuite les laisserait sur une numérotation qui n’est plus la bonne. La zone est donc scellée dès qu’un tel objet la prend : sa connectivité est figée pour de bon, add_cell lève une RuntimeError (« submesh is sealed »), et zone.is_sealed vaut True. Ce qui scelle :

• la construction d’un espace éléments finis scelle ses zones ;
• une matrice scelle les zones POI1 qui portent ses lignes et ses colonnes ;
• un champ aux nœuds scelle sa zone support, une zone POI1.

Contenir une zone ne la scelle pas : un maillage, une union |, une tranche la partagent sans la figer, et la zone peut grossir tant qu’aucun de ces objets ne s’y est posé. De même, pyrucast.mesh.merge_nodes(…, in_place=True), qui réécrit la connectivité en place, refuse une zone scellée.

Le champ aux nœuds mérite un mot. Posé sur une zone POI1, il prend cette zone pour support et la scelle. Posé sur une autre zone, il n’en indexe pas les mailles : il se pose sur son nuage de nœuds, une zone POI1 tirée d’elle (voir ci-dessous), et c’est ce nuage-là qui est scellé. La zone de départ reste modifiable ; si on la modifie, les champs déjà posés restent valides sur le nuage d’avant, et les suivants se posent sur un nouveau nuage.

Pour repartir d’une zone scellée, on en prend une copie : mesh.duplicate() rend un maillage neuf, non scellé, de même connectivité, sur les mêmes nœuds (seules leurs références augmentent). L’opérateur pyrucast.mesh.copy(mesh, new_nodes=…) fait de même, avec en plus le choix de créer des nœuds neufs aux mêmes positions.

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

pyrucast.FiniteElementSpace(mesh)  # indexe les mailles : scelle mesh[0]
print(mesh[0].is_sealed)  # True
# mesh[0].add_cell([b, n3, n4])  # RuntimeError: submesh is sealed: …

copie = mesh.duplicate()  # neuve, non scellée, sur les mêmes nœuds
copie.unit().add_cell([b, n3, n4])
print(copie.cell_count, mesh.cell_count)  # 2 1

Ce qu’une zone garde en cache

Une zone garde deux résultats tirés de sa connectivité. Chacun se calcule à la première demande, puis resert tel quel.

Son nuage de nœuds. pyrucast.mesh.to_poi1(mesh) rend, pour chaque zone, la zone POI1 de ses nœuds distincts, dans leur ordre de première apparition. Ce nuage est calculé une fois par zone, puis le même est rendu à chaque appel. Tout ce qui pose un champ aux nœuds sur un maillage passe par lui : NodeField(mesh, …), les supports des matrices d’un modèle — donc le résultat de solve(K, f) et le produit K * x —, pyrucast.node_field.restrict(champ, mesh), les sorties de divergence et de sum_to_nodes.

C’est ce qui permet à ces champs de se combiner. Deux champs aux nœuds ne s’additionnent ou ne se soustraient nœud à nœud que s’ils sont posés sur la même zone support ; sur deux supports distincts, même portant les mêmes nœuds, ils restent côte à côte, en deux zones (cf. Champ aux nœuds). Grâce au cache, une différence solve(K, f) - restrict(g, mesh) tombe sur le même support et se calcule directement. Le cache épargne aussi le dédoublonnage des nœuds, que chaque construction de modèle referait sinon.

L’index de ses nœuds. Une table « identifiant de nœud → rang dans la zone ». C’est par elle qu’un champ aux nœuds trouve la ligne d’un nœud (champ.value(nœud, "T"), set_value) et qu’une matrice trouve la ligne et la colonne d’un degré de liberté. Construite à la première recherche, elle sert ensuite à tous les champs et à toutes les matrices posés sur cette zone ; en pratique, c’est sur les zones POI1 support qu’elle se construit.

Les deux caches sont oubliés dès que la connectivité change (add_cell, merge_nodes(…, in_place=True)), et recalculés à la demande suivante : ils répondent toujours pour l’état courant de la zone. Un champ déjà posé garde son nuage, et restrict le porte sur le nouveau :

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

# Le nuage de nœuds d'une zone se calcule une fois : chaque appel rend
# ensuite la même zone POI1.
nuage = pyrucast.mesh.to_poi1(tri)
print(len(nuage | pyrucast.mesh.to_poi1(tri)))  # 1 — une seule zone

# Deux champs posés sur tri partagent donc leur support : ils se combinent
# nœud à nœud.
f = pyrucast.NodeField(tri, ["T"])
g = pyrucast.NodeField(tri, ["T"])
f[0].set_value(a, "T", 1.0)
g[0].set_value(a, "T", 0.25)
print((f - g).value(a, "T"))  # 0.75

# Ajouter une maille oublie le cache : le nuage suivant est une autre zone.
tri.unit().add_cell([b, n4, n3])
print(len(nuage | pyrucast.mesh.to_poi1(tri)))  # 2
print(f.node_count)  # 3 — f reste sur le nuage d'avant
r = pyrucast.node_field.restrict(f, tri)  # le porte sur le nouveau
print(r.node_count)  # 4

Les caches ne sont pas écrits par pyrucast.save : après une relecture, le nuage se recalcule, et un champ relu ne partage donc plus son support avec les champs qu’on pose ensuite sur le maillage relu. restrict(champ, mesh) l’y ramène.

Visualisation

mesh.plot(…) trace le maillage, chaque zone dans sa couleur (face_color), ou coloré par un champ. Voir Visualisation.

Côté Rust

Une zone stocke sa connectivité à plat : un Vec<NodeId> où la maille i occupe les entrées [i·n, (i+1)·n), n étant le nombre de nœuds du type. Elle porte aussi son Handle<Coords>, sa couleur, son drapeau de scellement et ses deux caches. Mesh est l’agrégat de Handle<SubMesh> ; il ne porte pas de Coords lui-même et le retrouve par sa première zone, et son check_push refuse une zone d’un autre Coords.

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

• Références. SubMesh::add_cell incrémente le compteur de chaque nœud et annule les incréments déjà faits si l’un échoue ; add_cell_taking adopte des références que l’appelant possède déjà (un nœud tout juste créé) ; le Drop de la zone les rend toutes, une fois, quand le dernier Handle disparaît (cf. Modèle mémoire).
• Construction en bloc. add_cell prend le verrou d’écriture du Coords et vide les caches à chaque maille. SubMesh::from_connectivity(coords, type, connectivité) reçoit tout le tableau, le valide et incrémente ses nœuds en une seule prise de verrou, et Coords::add_nodes crée des nœuds d’un coup : c’est par là que les opérateurs de maillage produisent leurs zones.
• Scellement. Les consommateurs passent par la fonction libre mesh::seal(&handle), qui ne prend le verrou d’écriture que si la zone n’est pas déjà scellée : sceller une zone qu’un lecteur tient déjà ne bloque pas. add_cell, add_cell_taking et remap_nodes rendent alors PyrucastError::MeshSealed.
• Caches. Deux OnceLock, lus par node_index() et to_poi1(), vidés par chaque mutation de la connectivité et marqués #[serde(skip)].
• Types d’éléments. ElementType porte le nom, le nombre de nœuds et la dimension topologique ; ce que l’élément sait de lui-même (repère, fonctions de forme, quadrature) vit dans src/atoms/element_kind/, un fichier par type — cf. Ajouter un élément fini.

Espace éléments finis (FiniteElementSpace/SubFiniteElementSpace)

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

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

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

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

Construire un espace

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

import pyrucast

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

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

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

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

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

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

La construction est refusée pour :

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

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

Lire un espace

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

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

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

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

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

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

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

Les tableaux sont rendus à plat, ligne par ligne :

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

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

Déplacer le maillage

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

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

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

Ce qui est mis en cache

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

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

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

Sauvegarde

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

Théorie

Élément de référence

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

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

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

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

Transformation isoparamétrique

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

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

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

Fonctions de forme de Lagrange

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

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

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

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

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

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

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

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

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

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

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

Dérivées de référence

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

Fonctions de forme d’Hermite (HERMITE3)

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

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

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

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

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

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

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

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

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

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

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

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

Pas de base du tout : MODEL_EMBEDDED

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

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

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

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

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

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

Quadrature de Gauss

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

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

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

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

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

Jacobien et grandeurs physiques

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

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

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

Dimensions admises :

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

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

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

Côté Rust

SubFiniteElementSpace détient :

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

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

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

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

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

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

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

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

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

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

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

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

Le reste de l’API sert les physiques :

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

Champ (Field/SubField)

pyrucast a deux familles de champs : le champ aux nœuds (NodeField), qui donne une valeur à chaque nœud, et le champ aux points de Gauss (ElementField), qui en donne une à chaque point d’intégration de chaque maille. Les deux se lisent et se calculent de la même façon, et cette page décrit ce qu’elles ont en commun ; les deux pages suivantes donnent ce qui est propre à chacune.

En Python, il n’y a pas de classe Field : tout ce qui suit s’écrit à l’identique sur NodeField, sur ElementField et sur leurs zones (SubNodeField, SubElementField). Côté Rust, ce socle commun est écrit une seule fois, dans deux traits (voir Côté Rust).

Zones, supports et composantes

Un champ est un agrégat de zones. Chaque zone :

• est posée sur un support : un nuage de nœuds (une zone de maillage POI1) pour un champ aux nœuds, un sous-espace d’espace éléments finis pour un champ aux points de Gauss. zone.support() le rend, et champ.support() rend celui de chaque zone, dans l’ordre ;
• porte des composantes nommées ("T", "u_x", "E", "sigma_xx"…) : au moins une, chaque nom une seule fois ;
• range une valeur par composante et par ligne de son support : un nœud, ou un couple (maille, point de Gauss).

Les règles de cohérence

Tout champ respecte deux règles, que chaque opération garde :

1. Une composante par support. Sur un même support, une composante est portée par une seule zone. Plusieurs zones peuvent partager un support si leurs composantes diffèrent : les matériaux de deux physiques sur un même maillage, par exemple.
2. Une valeur par nœud. Un nœud que plusieurs zones portent y a la même valeur, pour chaque composante qu’elles ont en commun. Deux zones de maillage qui se touchent partagent leurs nœuds de frontière : le champ y a une valeur, rangée dans chacune. La règle ne concerne que le champ aux nœuds, car deux sous-espaces n’ont aucun point de Gauss en commun.

La première se vérifie à l’entrée. NodeField(maillage, …) refuse un maillage qui liste deux fois la même zone, et ElementField(espace, …) un espace qui liste deux fois le même sous-espace : l’union maillage | maillage garde chaque zone une fois. add_sub refuse une zone qui porterait une composante déjà portée sur son support.

La seconde se vérifie à l’union, et en sortie d’une opération entre deux champs quand c’est nécessaire (plus bas). Seule l’écriture dans une zone peut encore la rompre : champ[i].set_value(…) ne touche pas les autres zones qui portent le même nœud. Sur un champ aux nœuds, champ.check() la vérifie à la demande.

Lire un champ

Les zones d’un même champ n’ont pas forcément les mêmes composantes. L’exemple coupe un carré en deux zones de maillage, gauche et droite, qui partagent la diagonale : la température T est connue sur les deux, le flux q sur celle de droite seulement.

# T sur les deux zones, q sur celle de droite seulement.
t = pyrucast.NodeField(gauche, ["T"])
for noeud, valeur in zip(n[:3], (0.0, 10.0, 20.0)):
    t[0].set_value(noeud, "T", valeur)
tq = pyrucast.NodeField(droite, ["T", "q"])
for noeud, valeur in zip((n[1], n[3], n[2]), (10.0, 30.0, 20.0)):
    tq[0].set_value(noeud, "T", valeur)
tq[0].add_to_component("q", 5.0)
champ = t | tq

print(champ.components())  # ['T', 'q'] : les composantes de toutes les zones
print(champ.min("T"), champ.max("T"))  # 0.0 30.0 : cherchés dans les deux zones
print(champ.min())  # 0.0 : toutes les valeurs, toutes composantes confondues
print(champ.sum("q"))  # 15.0 : les trois nœuds de droite, à 5
print(len(champ["q"]), champ["q"].components())  # 1 ['q'] : la zone sans q disparaît

Une composante qu’aucune zone ne porte lève une erreur.

champ["c"] laisse de côté une zone qui ne porte aucune des composantes demandées. Une zone qui n’en porte pas d’autres est partagée avec le champ d’origine, sans copie : écrire dans l’une modifie l’autre. Une zone dont il faut retirer des composantes est recopiée. rename_component partage de même les zones qui ne portent pas a.

sum(c) ajoute une valeur par nœud : un nœud que deux zones partagent n’a qu’une valeur (deuxième règle), et compte une fois. Ce n’est pas une intégrale, qui pèse chaque valeur par la taille des mailles : c’est pyrucast.measure.integral. Deux autres réductions s’appliquent aux deux familles, et comptent de même un nœud une fois : pyrucast.measure.xtx(champ), la somme des carrés des valeurs, et pyrucast.measure.xty(a, b), la somme des produits des valeurs que a et b portent sur les mêmes supports et sous les mêmes noms (voir Opérateurs sur les champs).

Calculer avec des champs

Par un scalaire

champ + s, champ - s, champ * s, champ / s et champ ** s appliquent l’opération à chaque valeur, de chaque composante, de chaque zone, et rendent un nouveau champ ; le scalaire se place d’un côté ou de l’autre. -champ change le signe de chaque valeur. Il en va de même sur une zone.

plus_chaud = champ + 5.0  # un nouveau champ : chaque valeur + 5
ecart = 100.0 - champ  # le scalaire se place d'un côté ou de l'autre
carre = champ**2.0  # puissance, valeur par valeur
print(plus_chaud.max("T"), ecart.min("T"), carre.max("q"))  # 35.0 70.0 25.0
print((-champ).min("T"), (2.0**champ).max("q"))  # -30.0 32.0

garde = champ
champ += 1.0  # un nouveau champ, rangé sous le même nom
print(champ.max("T"), garde.max("T"))  # 31.0 30.0 : `garde` a l'ancien

champ += 1.0 construit lui aussi un nouveau champ, et le range sous le même nom : une autre variable qui désignait l’ancien le garde tel quel. Pour modifier un champ lui-même, on passe par les opérations par composante.

Aucune valeur n’est protégée : une division par zéro donne inf, la puissance fractionnaire d’un nombre négatif donne nan, comme avec numpy. La forme pow(champ, s, m), avec un modulo, est refusée.

Sur une composante, en place

add_to_component(c, s), sub_to_component, mul_to_component et div_to_component modifient le champ lui-même : seule la composante c change, sur toutes les zones qui la portent. Elles lèvent une erreur si aucune zone ne porte c, et div_to_component refuse une division par zéro.

champ.add_to_component("T", -1.0)  # en place : T seulement, sur ses deux zones
champ.mul_to_component("q", 2.0)
print(champ.max("T"), champ.max("q"))  # 30.0 12.0

Sur une zone, set_uniform(c, v) donne la valeur v à la composante c en chaque nœud ou point de Gauss.

Entre deux champs

a + b, a - b, a * b, a / b et a ** b combinent les valeurs de deux champs de la même famille. Les valeurs s’associent par support et par composante :

• une composante que les deux champs portent sur un même support se combine valeur à valeur, quelle que soit la zone de ce support qui la range de chaque côté ;
• une composante qu’un seul des deux porte sur un support passe telle quelle, quelle que soit l’opération : a - b rend b, et non -b, sur une composante que seul b porte. C’est ce qui permet de ne toucher qu’à une composante : u + du, où du ne porte que u_x, laisse u_y tel quel, et u * f, où f ne porte que u_x, ne multiplie que u_x ;
• le résultat garde les zones de a, dans l’ordre. Une composante que seul b porte rejoint la zone de a sur ce support quand a n’en a qu’une ; sinon elle garde une zone à elle, comme une zone de b dont le support n’existe pas dans a.

Les deux champs n’ont donc besoin ni des mêmes composantes, ni des mêmes zones :

# a porte T et q, b porte T et p, sur la même zone.
a = pyrucast.NodeField(gauche, ["T", "q"])
a[0].add_to_component("T", 10.0)
a[0].add_to_component("q", 1.0)
b = pyrucast.NodeField(gauche, ["T", "p"])
b[0].add_to_component("T", 3.0)
b[0].add_to_component("p", 4.0)

d = a - b  # T se combine ; q et p, d'un seul côté, passent tels quels
print(d.components())  # ['T', 'q', 'p']
print(d.value(n[0], "T"), d.value(n[0], "q"), d.value(n[0], "p"))  # 7.0 1.0 4.0

champ + zone combine zone, selon les mêmes règles, avec les zones du champ posées sur son support, et recopie les autres ; il lève une erreur si aucune ne l’est. Entre deux zones, les supports doivent être les mêmes, sinon c’est une erreur.

Le résultat respecte les règles de cohérence. Une composante peut toutefois ne pas avoir été obtenue de la même façon d’une zone à l’autre : calculée ici, passée telle quelle là. Deux zones de ce genre qui partagent un nœud n’y donneraient pas la même valeur. Dans ce seul cas, le résultat est vérifié, et l’opération lève une erreur si deux zones ne s’accordent pas en un nœud.

Même support veut dire même objet

Deux zones s’associent quand leurs supports sont le même objet, jamais en comparant leurs nœuds. Deux régions qui se touchent, construites chacune sur sa zone de maillage, ont des supports distincts, même là où elles partagent des nœuds : g + dr laisserait chaque zone telle quelle, et un nœud de la frontière vaudrait 1 dans l’une et 2 dans l’autre. L’opération lève donc une erreur, qui renvoie à restrict_like.

Pour additionner deux régions, on les ramène d’abord sur les mêmes supports. restrict(champ, maillage) pose un champ aux nœuds sur les zones de maillage, et deux champs restreints au même maillage ont les mêmes supports ; restrict_like(b, a) pose b sur les supports de a (voir Opérateurs sur les champs). Un nœud que le champ restreint ne portait pas y vaut 0.

# Deux régions qui se touchent par la diagonale, chacune sur sa zone.
g = pyrucast.NodeField(gauche, ["T"])
g[0].add_to_component("T", 1.0)
dr = pyrucast.NodeField(droite, ["T"])
dr[0].add_to_component("T", 2.0)

try:
    g + dr  # deux supports distincts : n1 vaudrait 1 dans une zone, 2 dans l'autre
except RuntimeError as erreur:
    message = str(erreur)
    print(message)  # incoherent NodeField: node 1, component T: 1 ≠ 2: …

carre = gauche | droite
somme = g.restrict(carre) + dr.restrict(carre)  # les mêmes supports des deux côtés
print(somme.value(n[0], "T"), somme.value(n[1], "T"), somme.value(n[3], "T"))
# 1.0 3.0 2.0 : la diagonale reçoit les deux

Pour un champ aux points de Gauss, le support est le sous-espace : deux espaces construits sur le même maillage sont deux supports (voir Champ aux points de Gauss).

Fonctions mathématiques et masques

pyrucast.field réunit les fonctions qui s’appliquent valeur par valeur et rendent un nouveau champ de la même famille : abs, sqrt, exp, log (népérien), log10, cos, sin, tan, cosh, sinh et tanh. Chacune existe aussi en méthode : champ.sqrt() donne pyrucast.field.sqrt(champ). Elles ne protègent rien non plus : log d’un nombre négatif donne nan.

champ >= s, champ > s, champ <= s et champ < s rendent un masque : un champ de mêmes zones et de mêmes composantes, qui vaut 1 où la condition tient et 0 ailleurs (c’est mask, voir Opérateurs sur les champs). == et != ne comparent pas les valeurs : champ == 0.0 vaut False.

import pyrucast as pc

amorti = pc.field.exp(champ * -0.1)  # e^(-0.1 v), valeur par valeur
print(round(amorti.min("T"), 4))  # 0.0498 : e^-3, au nœud le plus chaud
print(pc.field.cos(champ).max("T"))  # 1.0 : cos 0, au nœud à 0
print(champ.sqrt().max("T"))  # 5.477… : chaque fonction est aussi une méthode

chaud = champ >= 15.0  # 1 où la valeur atteint 15, 0 ailleurs
print(chaud.value(n[0], "T"), chaud.value(n[3], "T"))  # 0.0 1.0

L’union |

L’union compose des zones sans toucher aux valeurs ; c’est l’arithmétique qui combine des valeurs. Elle respecte les règles de cohérence, à sa manière pour chaque famille :

• champ aux nœuds : les zones posées sur un même support sont fondues en une. Un nœud doit avoir la même valeur dans toutes les zones qui le portent, pour chaque composante qu’elles ont en commun, sinon l’union lève une erreur ;
• champ aux points de Gauss : deux zones posées sur un même sous-espace restent côte à côte si leurs composantes sont différentes ; une même composante portée deux fois sur un sous-espace est refusée.

consolidate() fond, à la demande, les zones qui partagent un support (voir Champ aux nœuds et Champ aux points de Gauss).

Ce qui est mis en cache

Un champ ne garde aucun cache : une zone ne contient que le nom de ses composantes, ses valeurs et une référence à son support. Ce qui ne se calcule qu’une fois vit sur les supports, et c’est ce qui rend les opérations ci-dessus possibles ou rapides :

• le nuage de nœuds d’une zone de maillage. NodeField(maillage, …), restrict et les opérateurs qui produisent un champ aux nœuds posent leurs zones sur le nuage des nœuds de chaque zone de maillage. Ce nuage est calculé une fois, puis rendu tel quel à chaque demande, tant que la zone de maillage ne change pas. Deux champs construits sur le même maillage reçoivent donc le même support et se combinent nœud à nœud : c’est ce qui fait marcher a - b et restrict(…, carre) dans les exemples ci-dessus ;
• la table qui donne la ligne d’un nœud dans ce nuage, construite à la première lecture d’un nœud : value(nœud, c), set_value et les opérateurs trouvent la ligne d’un nœud sans parcourir la liste ;
• les tables de l’espace éléments finis (points et poids de Gauss, fonctions de forme et leurs dérivées) et la coloration de ses mailles, pour un champ aux points de Gauss (voir Champ aux points de Gauss).

Côté Rust

Le socle commun est écrit une seule fois, dans deux traits de src/containers/field.rs. Python n’en voit que le résultat : les méthodes et les opérateurs décrits plus haut, présents sur les quatre classes.

• SubField, la zone. Un type n’a que cinq choses à fournir : support(), components(), values(), values_mut() et same_support_with(…), une zone neuve, à zéro, sur le même support. Tout le reste en découle : component_index, min, max, sum, xtx, dot, pscal, map_all, add_to_component…, merge_components, select_components, rename_component. SubNodeField et SubElementField l’implémentent. SubNodeField fournit en plus row_points(), le nœud de chaque ligne, lu sur place dans le support : c’est ce qui permet de compter un nœud une fois (la première zone qui le porte le compte, les suivantes le sautent) et de vérifier la deuxième règle. Une zone aux points de Gauss garde la réponse par défaut, None.
• Field, l’agrégat. Il est implémenté d’office pour tout Aggregate dont la zone est un SubField : NodeField et ElementField le reçoivent sans une ligne de code : components, min, max, sum, map_subs, combine_scalar, map_all, merge_field, merge_subfield, dot_field, pscal_field, filter_components, rename_component. Les réductions comptent une fois un nœud que plusieurs zones portent ; merge_field, merge_subfield, dot_field et pscal_field associent les valeurs par support et par composante.
• Règles de cohérence. add_sub refuse une zone qui porterait une composante déjà portée sur son support : c’est le crochet check_add, qu’une union ne passe pas, puisque son finalize vérifie ou fond le résultat entier. Les opérateurs qui empilent des contributions partielles avant de les additionner nœud à nœud (divergence, forces internes et externes, sum_to_nodes) ne vérifient pas la deuxième règle en chemin : ce ne sont pas encore des champs.

Le contrat ne dit pas ce qu’est une ligne — un nœud pour SubNodeField, un couple (maille, point de Gauss) pour SubElementField. Il fixe seulement la disposition des valeurs dans un tableau plat, où la composante varie le plus vite :

   zone aux composantes [c0, c1]
   values = [ ligne0.c0, ligne0.c1, ligne1.c0, ligne1.c1, … ]

#[test]
fn le_contrat_des_champs_s_ecrit_une_fois() -> Result<()> {
    let coords = Handle::new(Coords::new(1)?);
    let a = Node::create_in(coords.clone(), &[0.0])?;
    let b = Node::create_in(coords, &[1.0])?;
    let support = mesh::poi1_from_nodes(&[a.clone(), b.clone()])?;

    let t = NodeField::new(&support, vec!["T".into()])?;
    t.add_to_component("T", 2.0)?; // Field : en place, sur chaque zone qui porte T
    let tq = NodeField::new(&support, vec!["T".into(), "q".into()])?;
    tq.add_to_component("T", 3.0)?;

    // Les opérateurs rendent un `Result` : combiner deux champs peut échouer.
    let s = (&t + &tq)?; // même support : T se combine, q passe tel quel
    assert_eq!(s.value(a.id(), "T")?, 5.0);
    assert_eq!((2.0 * &t)?.value(b.id(), "T")?, 4.0);

    // Rust n'a pas d'opérateur de puissance : on passe par les primitives.
    let carre = t.combine_scalar(f64::powf, 2.0)?;
    let puissance = t.merge_field(&tq, f64::powf)?; // 2³
    assert_eq!(carre.value(a.id(), "T")?, 4.0);
    assert_eq!(puissance.value(a.id(), "T")?, 8.0);

    // Une zone se lit par le trait `SubField`.
    let zone = t.get(0)?;
    assert_eq!(zone.read().component_index("T"), Some(0));
    assert_eq!(SubField::max(&*zone.read(), Some("T"))?, 2.0);
    Ok(())
}

• Opérateurs. Trois macros écrivent + - * / : impl_subfield_scalar_ops! (une zone et un scalaire, d’un côté ou de l’autre), impl_subfield_field_ops! (deux zones, par merge_components) et impl_field_ops! (deux agrégats par merge_field, un agrégat et un scalaire par combine_scalar). Combiner peut échouer, sauf une zone avec un scalaire : ces opérateurs rendent un Result. Une zone prise par valeur est modifiée sur place, sans copie ; prise par référence, elle est d’abord copiée.
• Puissance. Rust n’a pas d’opérateur ** : combine_scalar(f64::powf, s) et merge_field(&autre, f64::powf) font ce que fait ** en Python.
• Support. same_support compare les Handle par identité (Handle::same_object).
• Mêmes composantes exigées. check_same_components refuse deux zones dont les supports ou les composantes diffèrent. L’interpolation d’une Evolution l’appelle avant merge_components : deux valeurs tabulées d’un même champ doivent porter les mêmes composantes.
• Lecture par rang. Un noyau de calcul ne cherche pas une composante par son nom à chaque point : resolve_components traduit, une fois par zone, une liste de noms en rangs dans la zone, et resolve_optional_components marque ABSENT_COMPONENT les noms que la zone ne porte pas.
• Fonctions mathématiques. Les traits MapValues et Pscal donnent un nom commun aux quatre types (zone ou agrégat, nœuds ou points de Gauss) : ops::field::cos, ops::field::psca… sont écrits une seule fois.
• Vue de lecture. Field::view prend un verrou de lecture par zone et le garde le temps de la vue. Les opérateurs qui lisent beaucoup (gradient, deformation, divergence, restrict, integral, le tracé…) lisent ainsi les valeurs sur place, sans copie ; d’autres lectures peuvent avoir lieu en même temps, une écriture attend la fin de la vue.
• Parallélisme. min, max, sum, xtx, map_all… parcourent les valeurs en parallèle. min et max donnent le même résultat quel que soit le nombre de threads ; sum, xtx et dot peuvent différer au dernier chiffre.

Champ aux nœuds (NodeField/SubNodeField)

Un champ aux nœuds porte des valeurs par nœud : une température, un déplacement, des forces nodales, les positions des nœuds. Là où un champ aux points de Gauss donne une valeur à chaque point d’intégration de chaque maille, celui-ci en donne une à chaque nœud, que les mailles voisines partagent. C’est un agrégat de zones :

• une zone (SubNodeField) range les valeurs des nœuds d’un nuage de nœuds, une zone de maillage POI1 ;
• un champ (NodeField) est une liste de zones, qui n’ont pas forcément les mêmes composantes : un champ multiphysique porte T sur tout le domaine et ux, uy sur la partie solide seulement.

Les règles communes aux deux familles de champs (composantes, règles de cohérence, statistiques, arithmétique) sont décrites dans Champ ; cette page en donne ce qui est propre aux nœuds.

Créer et lire un champ

pyrucast.NodeField(maillage, ["T"]) crée une zone par zone du maillage, avec ces composantes, toutes les valeurs à zéro ; sur une seule zone de maillage (SubMesh), le champ a une zone. NodeField.with_components_per_submesh(maillage, [["T", "ux"], ["T"]]) donne une liste de composantes par zone. Une liste vide, un nom en double, un nombre de listes qui ne correspond pas aux zones ou un maillage qui liste deux fois la même zone sont refusés.

Le champ se lit nœud par nœud : champ.value(nœud, "T") lit la première zone qui porte ce nœud et cette composante, et lève une erreur si aucune ne les porte. champ.values(nœuds, "T") lit plusieurs nœuds d’un coup et rend leurs valeurs dans le même ordre ; nœuds est une liste de Node, ou une zone ou un maillage dont la connectivité est lue dans l’ordre (pour une zone POI1, ses points). Un nœud d’un autre Coords est refusé.

On écrit dans une zone : champ[i] ou champ.unit() la rendent, et ce sont des vues (écrire dans la zone, c’est écrire dans le champ). Dans une zone, un nœud se désigne par lui-même (zone[nœud, "T"], set_value) ou par sa ligne, son rang dans le support de la zone (get(ligne, rang), set, node_values). Les opérations par composante, comme add_to_component, s’appliquent au champ entier, en place, dans toutes les zones qui portent la composante.

import pyrucast

# Un carré coupé en deux triangles, deux zones de maillage qui partagent la
# diagonale a–d.
c = pyrucast.Coords(2)
a, b, d, e = (c.add_node(p) for p in ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]))
bas = pyrucast.Mesh(c, "TRI3")
bas.unit().add_cell([a, b, d])
haut = pyrucast.Mesh(c, "TRI3")
haut.unit().add_cell([a, d, e])
carre = bas | haut

# Une zone par zone de maillage, toutes les valeurs à zéro.
T = pyrucast.NodeField(carre, ["T"])
print(len(T), T.node_count)  # 2 4 — a et d sont rangés dans les deux zones

# En place, dans toutes les zones qui portent T ; puis dans une zone seulement.
T.add_to_component("T", 20.0)
T[1].set_value(e, "T", 25.0)  # e n'appartient qu'à la zone du haut

# Lire par le champ : la première zone qui porte le nœud répond.
print(T.value(e, "T"))  # 25.0
print(T.values([a, b, d, e], "T"))  # [20.0, 20.0, 20.0, 25.0]
print(T.sum("T"))  # 85.0 — a et d ne comptent qu'une fois

# Une zone est une vue : écrire dans zone, c'est écrire dans T.
zone = T[0]
print(zone)  # SubNodeField: 3 node(s), 1 component(s) [T]
print(zone[d, "T"], zone.get(2, 0))  # 20.0 20.0 — d est la ligne 2 de la zone

# Écrire un nœud partagé dans une seule zone rompt la deuxième règle.
zone[d, "T"] = 30.0
try:
    T.check()
except RuntimeError as err:
    print(err)  # incoherent NodeField: node 2, component T: 30 ≠ 20
T[1][d, "T"] = 30.0
T.check()  # les deux zones s'accordent de nouveau

# Une liste de composantes par zone.
mixte = pyrucast.NodeField.with_components_per_submesh(carre, [["T", "ux"], ["T"]])
print(mixte.components())  # ['T', 'ux']

Le support : un nuage de nœuds

Chaque zone est posée sur une zone de maillage POI1 : une liste de nœuds, un par maille. La ligne i de la zone est le i-ème nœud de ce support ; le champ ne range que ses valeurs et lit ses nœuds dans le support.

• Sur une zone POI1, le champ prend cette zone pour support, telle quelle. Il la refuse si elle liste deux fois un même nœud, puisqu’il range une ligne par nœud ; pyrucast.mesh.to_poi1 en donne les nœuds distincts.
• Sur une zone d’un autre type (TRI3, HEX8…), il prend son nuage de nœuds : la zone POI1 de ses nœuds distincts, dans leur ordre de première apparition, celle que rend pyrucast.mesh.to_poi1. La zone de maillage le calcule une fois et le garde en cache : deux champs construits sur le même maillage reçoivent le même support.

Le support est scellé, puisque le champ range ses lignes dans l’ordre de ses nœuds. La zone TRI3 d’où vient le nuage, elle, reste modifiable. Si on lui ajoute une maille, elle oublie son nuage : les champs déjà posés restent valides sur l’ancien, et un champ construit ensuite se pose sur un nouveau nuage. Les deux ne se combinent plus nœud à nœud ; pyrucast.node_field.restrict(champ, maillage) porte l’ancien sur le nouveau (voir Maillage).

champ.support() rend le support de chaque zone, dans l’ordre, et zone.support() celui d’une zone. C’est le même objet, pas une copie : un champ construit dessus se combine nœud à nœud avec le premier. Le champ garde son support en vie, mais pas la zone de maillage dont il a été tiré : il n’a plus besoin d’elle.

# Posé sur une zone TRI3, le champ prend son nuage de nœuds : une zone POI1,
# scellée. La zone TRI3 elle-même reste modifiable.
s = T.support()
print(s[0].element_type, s[0].cell_count)  # POI1 3
print(s[0].is_sealed, carre[0].is_sealed)  # True False

# Ce nuage est celui que rend to_poi1, le même objet : l'union ne garde
# que deux zones, pas quatre.
print(len(s | pyrucast.mesh.to_poi1(carre)))  # 2

# Une zone POI1 est elle-même le support.
points = pyrucast.mesh.poi1_from_nodes([b, d])
F = pyrucast.NodeField(points, ["F"])
print(len(F.support() | points), points[0].is_sealed)  # 1 True

Nœuds partagés entre zones

Deux zones de maillage qui se touchent partagent leurs nœuds de frontière, et un champ posé sur elles range ces nœuds dans chaque zone, comme a et d dans le premier exemple. La deuxième règle de cohérence veut qu’un tel nœud ait la même valeur dans toutes les zones qui le portent, pour chaque composante qu’elles ont en commun. Le champ s’appuie sur elle :

• value et values lisent la première zone qui porte le nœud ;
• sum, node_count, pyrucast.measure.xtx et xty comptent un nœud partagé une fois ;
• add_to_component et les autres opérations par composante modifient toutes les zones à la fois, et gardent donc la règle.

L’union et les opérations entre deux champs la vérifient quand c’est nécessaire. Écrire dans une zone, en revanche, ne touche pas les autres : T[0][d, "T"] = 30.0 laisse l’ancienne valeur de d dans la zone du haut. champ.check() compare alors les valeurs exactement et lève une erreur au premier désaccord.

Combiner des champs

Les opérations par un scalaire et entre deux champs suivent les règles de Champ : les valeurs s’associent par support et par composante, et deux zones s’associent quand leurs supports sont le même objet.

L’union a | b réunit les zones sans toucher aux valeurs :

• les zones posées sur un même support sont fondues en une, qui porte l’union de leurs composantes ; une composante portée des deux côtés doit y avoir les mêmes valeurs ;
• les zones posées sur des supports distincts restent côte à côte, et doivent s’accorder sur les nœuds qu’elles partagent.

Un désaccord lève une erreur. pyrucast.node_field.merge(a, b) est cette même union. add_sub ajoute une zone telle quelle, sans rien fondre, et refuse une zone qui porterait une composante déjà portée sur son support ; consolidate() fond ensuite les zones qui partagent un support, après la même vérification.

# Sur un même support, l'union fond les zones en une.
q = pyrucast.NodeField(carre, ["q"])
tout = T | q
print([z.components() for z in tout])  # [['T', 'q'], ['T', 'q']]

# Deux supports distincts restent côte à côte, mais doivent s'accorder sur
# les nœuds qu'ils partagent.
coin = pyrucast.NodeField(pyrucast.mesh.poi1_from_nodes([d]), ["T"])
try:
    T | coin
except RuntimeError as err:
    print(err)  # incoherent NodeField: node 2, component T: 30 ≠ 0
coin.unit()[d, "T"] = 30.0
print(len(T | coin))  # 3

# add_sub ajoute une zone telle quelle ; consolidate fond celles qui
# partagent un support.
w = pyrucast.NodeField(carre, ["T"])
w.add_sub(q[0])
print(len(w), len(w.consolidate()))  # 3 2

D’où viennent ces champs, et où ils vont

Les opérateurs qui produisent un champ aux nœuds sont rangés dans pyrucast.node_field (détail dans Opérateurs sur les champs et Assemblage) :

La résolution en produit un aussi : pyrucast.solver.solve(K, f) pose sa solution sur les supports des blocs de la matrice, ceux que restrict(…, maillage) donne sur le même maillage (voir Solveur).

Ces champs sont ensuite consommés :

• portés aux points de Gauss : interp_to_gauss, gradient, deformation, beam_deformation, shell_deformation, strain (dans pyrucast.element_field) ;
• comme second membre de solve, où un degré de liberté qu’aucune zone ne porte vaut zéro ;
• pour déplacer les nœuds : pyrucast.coords.set, displace, displace_from, qui traitent une fois un nœud partagé ;
• par pyrucast.model.update_contact, les réductions de pyrucast.measure (integral, xtx, xty), l’export (pyrucast.export.to_arrays(…, node_fields=[…]), export_vtk) et le tracé.

Ce qui est mis en cache

Un champ aux nœuds ne garde aucun cache : une zone ne contient que ses valeurs, les noms de ses composantes et son support. Il s’appuie sur deux caches des zones de maillage (détaillés dans Maillage) :

• le nuage de nœuds, gardé par la zone de maillage d’origine, calculé à la première demande puis rendu tel quel. C’est lui qui donne le même support à NodeField(maillage, …), à restrict(…, maillage), à la solution de solve et aux sorties de divergence et de sum_to_nodes sur ce maillage : ils se combinent nœud à nœud. Il épargne aussi le dédoublonnage des nœuds de la zone à chaque nouveau champ ;
• l’index des nœuds, gardé par le support : la table « nœud → ligne », construite à la première lecture par nœud (value, set_value, values, les opérateurs qui lisent le champ nœud par nœud). Elle sert ensuite à tous les champs posés sur ce support : trouver la ligne d’un nœud ne parcourt plus la liste.

Le support étant scellé, son index ne change plus. Le nuage, lui, est oublié dès que la zone de maillage d’origine change.

Visualisation

champ.plot(component="T") trace un nuage de points colorés : le support POI1 n’a pas de mailles, aucune surface ne s’en déduit. Pour colorer les mailles, on passe le champ au maillage : maillage.plot(field=champ, component="T"). Voir Visualisation.

Sauvegarde

pyrucast.save écrit le champ avec son support et son Coords. Des champs sauvegardés dans un même fichier partagent leur support à la relecture et se combinent nœud à nœud. Le nuage mis en cache n’est pas écrit : un champ construit sur le maillage relu se pose sur un nouveau nuage, et restrict(champ, maillage) y ramène un champ relu. Les nœuds relus appartiennent au Coords relu : value refuse un Node du script qui a écrit le fichier, et on lit avec le nœud relu, o["maillage"].coords().acquire(id).

Côté Rust

SubNodeField porte le Handle<SubMesh> de son support, ses noms de composantes et un Vec<f64> à plat, rangé ligne → composante :

values[ligne * component_count + c]

Les composantes d’un nœud sont donc contiguës : node_values(ligne) en rend la tranche sans copie.

Le Handle<SubMesh> du support est le compagnon POI1 lui-même, sans objet intermédiaire. SubMesh::to_poi1 garde ce compagnon dans un OnceLock et en rend un clone du Handle ; le champ garde ce clone. Un Handle étant un Arc<RwLock<SubMesh>>, le champ et la zone d’origine tiennent deux clones du même Arc, et le champ ne référence pas la zone d’origine.

NodeField est l’agrégat de Handle<SubNodeField>. Son check_push refuse une zone d’un autre Coords, son check_add une composante déjà portée sur le support de la zone (première règle), et son finalize, appelé à la fin de chaque union, est ops::node_field::consolidate.

#[test]
fn un_champ_aux_noeuds_s_ecrit_par_zone_et_se_lit_par_agregat() -> Result<()> {
    let coords = Handle::new(Coords::new(2)?);
    let n: Vec<Node> = [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]
        .iter()
        .map(|p| Node::create_in(coords.clone(), p))
        .collect::<Result<_>>()?;
    let mut tri = Mesh::from_submesh(SubMesh::new(coords, ElementType::TRI3));
    tri.add_cell(&[n[0].id(), n[1].id(), n[2].id()])?;

    // Une zone par zone de maillage, posée sur son nuage de nœuds : le
    // `Handle` même que `to_poi1` garde en cache, sans objet intermédiaire.
    let u = NodeField::new(&tri, vec!["ux".into(), "uy".into()])?;
    let zone = u.get(0)?;
    let nuage = tri.get(0)?.read().to_poi1()?;
    assert!(zone.read().support().same_object(&nuage));

    // Écriture par la zone, lecture par le champ.
    zone.write().set_value(n[1].id(), "ux", 1.5)?;
    assert_eq!(u.value(n[1].id(), "ux")?, 1.5);
    assert_eq!(u.value(n[0].id(), "ux")?, 0.0); // valeur initiale
    u.check()?;
    Ok(())
}

• Contrat commun. Les noms de composantes, min/max/sum, les opérations par composante et l’arithmétique viennent des traits SubField et Field (src/containers/field.rs), partagés avec le champ aux points de Gauss.
• Ligne d’un nœud. SubNodeField::index_of passe par SubMesh::node_index(), un OnceLock<HashMap<NodeId, usize>> du support ; index_of_with prend un verrou de lecture déjà tenu, pour une boucle qui cherche beaucoup de nœuds.
• Lecture par les opérateurs. NodeFieldView tient un verrou de lecture sur chaque zone et sur son support pendant tout un opérateur (gradient, restrict, project, le tracé…), et lit les valeurs en place.
• Un nœud une fois. first_rows parcourt les zones dans l’ordre et donne chaque nœud une fois, à la première zone qui le porte : c’est ce qu’utilisent node_ids, node_count et le déplacement des nœuds.

Champ aux points de Gauss (ElementField/SubElementField)

Un champ aux points de Gauss porte des valeurs par maille et par point de Gauss : là où un champ aux nœuds donne une valeur à chaque nœud, celui-ci en donne une à chaque point d’intégration de chaque maille. Deux mailles voisines ne partagent donc rien, et le champ peut sauter d’une maille à l’autre. C’est un agrégat de zones :

• une zone (SubElementField) porte les valeurs d’un sous-espace d’un espace éléments finis, c’est-à-dire d’une zone de maillage munie de sa quadrature ;
• un champ (ElementField) est une liste de zones, qui n’ont pas forcément les mêmes composantes.

On y range ce qui vit aux points où se calculent les intégrales :

• les propriétés matériau (module d’Young, coefficient de Poisson, conductivité, masse volumique…) ;
• les variables internes d’une loi de comportement (déformation plastique, endommagement…), gardées d’un pas de calcul au suivant ;
• les grandeurs dérivées d’une solution (déformations, contraintes, gradients, flux…).

On en construit rarement à la main : la plupart sortent des opérateurs (material_field, gradient, integrate_behavior…). Les règles communes aux deux familles de champs (composantes, statistiques, arithmétique) sont décrites dans Champ ; cette page en donne ce qui est propre aux points de Gauss.

Créer et lire un champ

pyrucast.ElementField(fes, ["E", "nu"]) crée une zone par sous-espace de fes, avec ces composantes, toutes les valeurs à zéro. ElementField.with_components_per_subspace(fes, [["E", "nu"], ["k"]]) donne une liste de composantes par sous-espace. Une liste vide, un nom en double ou un nombre de listes qui ne correspond pas aux sous-espaces sont refusés.

On lit et on écrit dans une zone : champ[i] ou champ.unit() la rendent, et ce sont des vues (écrire dans la zone, c’est écrire dans le champ). Un point s’y désigne par trois indices : la maille, dans l’ordre des mailles de la zone de maillage ; le point de Gauss dans la maille, dans l’ordre de la quadrature ; la composante, par son nom ou par son rang. Un indice hors bornes ou une composante inconnue lèvent une erreur.

import pyrucast

# Deux triangles et leur espace éléments finis : 3 points de Gauss par TRI3.
c = pyrucast.Coords(2)
a = c.add_node([0.0, 0.0])
b = c.add_node([1.0, 0.0])
d = c.add_node([1.0, 1.0])
e = c.add_node([0.0, 1.0])
mesh = pyrucast.Mesh(c, "TRI3")
mesh.unit().add_cell([a, b, d])
mesh.unit().add_cell([a, d, e])
fes = pyrucast.FiniteElementSpace(mesh)

# Une zone par sous-espace de fes, toutes les valeurs à zéro.
mat = pyrucast.ElementField(fes, ["E", "nu"])
z = mat.unit()  # la zone, une vue : écrire dans z, c'est écrire dans mat
print(z)  # SubElementField: 2 cell(s) × 3 gauss × 2 component(s) [E, nu]

# Écrire : toute la zone, une maille entière, un point.
z.set_uniform("E", 210e9)
z.set_uniform("nu", 0.3)
z.set_cell_uniform(1, "E", 70e9)  # la maille 1, à tous ses points de Gauss
z[1, 2, "nu"] = 0.33  # la maille 1, son point de Gauss 2

# Lire : par nom, par indice de composante, ou un point entier.
print(z.value(0, 0, "E"))  # 210000000000.0
print(z[1, 2, "nu"])  # 0.33
print(z.get(1, 0, 0))  # 70000000000.0
print(z.point_values(1, 2))  # [70000000000.0, 0.33]
print(mat.max("E"), mat.min("E"))  # 210000000000.0 70000000000.0

# Une liste de composantes par sous-espace (fes n'en a qu'un).
cond = pyrucast.ElementField.with_components_per_subspace(fes, [["k"]])
print(cond.components())  # ['k']

Les positions des points de Gauss ne sont pas dans le champ : leurs coordonnées dans l’élément de référence et leurs poids sont sur l’espace (fes[i].gauss_xi(g), fes[i].gauss_weight(g)), et leur position réelle se calcule à partir des nœuds.

Composantes

Chaque zone a ses propres composantes, au moins une, chacune nommée une seule fois. Au niveau du champ, components() rend leur union, dans l’ordre de première apparition : une composante peut n’exister que sur certaines zones. C’est le cas d’un champ matériau multiphysique, qui porte E et nu sur la zone de l’élasticité et k sur celle de la thermique.

champ["E"] et champ.filter_components(["E", "nu"]) rendent un nouveau champ réduit à ces composantes. Une zone qui n’en porte aucune disparaît ; une zone qui n’en porte pas d’autres est partagée avec le champ d’origine, sans copie, si bien qu’écrire dans l’une modifie l’autre ; une zone dont il faut retirer des composantes est recopiée.

Le support : un espace éléments finis

Chaque zone est posée sur un sous-espace de l’espace éléments finis, qui fixe ses mailles (celles de sa zone de maillage) et son nombre de points de Gauss par maille (celui de sa quadrature). Ces deux dimensions sont lues à la création de la zone et ne changent plus : l’espace a scellé sa zone de maillage, aucune maille ne peut y être ajoutée. Les nœuds, eux, peuvent bouger : les valeurs du champ restent ce qu’elles sont, et les points de Gauss suivent les mailles.

Deux champs se combinent point à point seulement s’ils sont posés sur le même espace, au sens du même objet. Deux appels à FiniteElementSpace(mesh) sur un même maillage créent deux espaces distincts, et leurs champs restent côte à côte au lieu de s’additionner (voir l’exemple ci-dessous). On construit donc un espace une fois et on le passe à tous les opérateurs ; material_field(modele, …) pose ses zones sur les sous-espaces des sous-modèles, ceux-là mêmes.

Un champ rend aussi le sien : champ.support() donne le sous-espace de chaque zone, dans l’ordre (champ.support()[i] est champ[i].support()), et zone.support() celui d’une zone. C’est le même objet, pas une copie : un champ construit dessus se combine point à point avec le premier. Deux zones posées sur un même sous-espace le citent deux fois ; s | s le garde une fois.

Un champ reste sur sa quadrature : aucun opérateur ne transporte des valeurs d’une quadrature à une autre. Pour obtenir un champ sur un autre espace, on l’interpole depuis les nœuds (interp_to_gauss), ou on l’intègre maille par maille vers un espace à un point par maille (cell_integral).

Un champ garde son espace en vie, et l’espace garde son maillage et son Coords : qu’un script ne tienne plus de variable sur l’espace ou le maillage ne gêne pas le champ. Il ne compte pas les nœuds pour autant ; c’est la zone de maillage qui les compte (voir Maillage).

Combiner des champs

Les opérations suivent les règles de Champ :

• une opération par un scalaire (champ * 1.1, champ + 2.0…) s’applique à toutes les valeurs et rend un nouveau champ ;
• a + b, a - b, a * b, a / b, a ** b entre deux champs combinent une composante que les deux portent sur un même sous-espace, même quand chaque côté la range dans une autre de ses zones sur ce sous-espace ; une composante ou une zone présente d’un seul côté passe telle quelle ;
• champ + zone ne touche que les zones du champ posées sur le même sous-espace ;
• entre deux zones, les sous-espaces doivent être les mêmes, sinon c’est une erreur ;
• abs, sqrt, exp, cos… s’appliquent valeur par valeur, et champ >= s rend un masque de 0 et de 1 (mask).

L’union a | b compose des zones sans toucher aux valeurs. Elle garde côte à côte deux zones posées sur le même sous-espace si leurs composantes sont disjointes, et refuse qu’une même composante y soit portée deux fois. consolidate() fond en une seule zone celles qui partagent un sous-espace. Contrairement au champ aux nœuds, l’union ne les fond pas d’elle-même.

# Deux champs posés sur le même espace se combinent point à point.
k1 = pyrucast.ElementField(fes, ["k"])
k2 = pyrucast.ElementField(fes, ["k"])
k1.unit().set_uniform("k", 2.0)
k2.unit().set_uniform("k", 0.5)
somme = k1 + k2
print(len(somme), somme.unit().value(0, 0, "k"))  # 1 2.5

# Un second espace, même sur le même maillage, est un autre support : les
# zones restent côte à côte, rien n'est additionné.
autre = pyrucast.ElementField(pyrucast.FiniteElementSpace(mesh), ["k"])
print(len(k1 + autre))  # 2

# L'espace d'un champ, rendu tel quel : un champ posé dessus s'y combine.
print(len(pyrucast.ElementField(k1.support(), ["k"]) + k1))  # 1

# L'union garde côte à côte deux zones aux composantes disjointes ;
# consolidate les fond en une.
tout = k1 | mat
print(len(tout), tout.components())  # 2 ['k', 'E', 'nu']
print(len(tout.consolidate()))  # 1

sum(c) additionne les valeurs de tous les points de Gauss, sans poids : ce n’est pas une intégrale. L’intégrale d’une composante sur chaque maille est cell_integral (exemple de la section suivante).

D’où viennent ces champs, et où ils vont

Les opérateurs qui produisent un champ aux points de Gauss sont rangés dans pyrucast.element_field (détail dans Opérateurs sur les champs, Construction et Comportement) :

material_field crée une zone par sous-modèle qui consomme du matériau, sur l’espace de ce sous-modèle, et ignore les autres (Dirichlet…). Il refuse une liste à laquelle manque une composante qu’une physique attend. Deux sous-modèles posés sur un même sous-espace qui lisent une même composante — une conduction et une élasticité qui lisent toutes deux rho — partagent une seule zone, puisqu’un sous-espace porte chaque composante dans une seule zone ; avec material_field_per_sub_model, ils doivent alors lui donner la même valeur.

Ces champs sont ensuite consommés par l’assemblage (stiffness, mass, tangent, internal_forces, external_forces, cf. Assemblage), ramenés aux nœuds (divergence, sum_to_nodes), tracés (plot) ou exportés (pyrucast.export.to_arrays(…, element_fields=[…])).

# Le matériau d'un modèle : une zone par sous-modèle qui en consomme.
elasticite = pyrucast.model.elasticity(fes, "plane_stress")
modele = elasticite | pyrucast.model.heat_conduction(fes)
materiau = pyrucast.element_field.material_field(
    modele, [("E", 210e9), ("nu", 0.3), ("k", 45.0)]
)
print([zone.components() for zone in materiau])  # [['E', 'nu'], ['k']]

# Une température T = x aux nœuds, portée aux points de Gauss.
T = pyrucast.NodeField(mesh, ["T"])
for n in (a, b, d, e):
    T[0].set_value(n, "T", n.position()[0])
T_g = pyrucast.element_field.interp_to_gauss(T, fes)
print([T_g.unit().value(0, g, "T") for g in range(3)])  # [0.5, 1.0, 0.5]
grad = pyrucast.element_field.gradient(T, fes)
print(grad.components())  # ['grad_T_x', 'grad_T_y']
print(grad.unit().point_values(0, 0))  # [1.0, 0.0]

# Une somme sur les points de Gauss n'est pas une intégrale : l'intégrale
# par maille se pose sur un espace à un point par maille.
un = pyrucast.ElementField(fes, ["un"]) + 1.0
centres = pyrucast.FiniteElementSpace(mesh, quadrature="REDUCED")
aires = pyrucast.element_field.cell_integral(un, centres)
print(un.sum("un"), aires.sum("un"))  # 6.0 1.0

Ce qui est mis en cache

Un champ aux points de Gauss ne garde aucun cache : il ne contient que ses valeurs, les noms de ses composantes et ses deux dimensions. Ce qui ne se calcule qu’une fois vit sur l’espace éléments finis qui le porte :

• les tables de référence (coordonnées et poids des points de Gauss, fonctions de forme et leurs dérivées en ces points), calculées à la création de l’espace ;
• la coloration des mailles : des groupes de mailles qui ne partagent aucun nœud. Elle sert aux opérateurs qui accumulent des contributions de mailles sur les nœuds (l’assemblage d’une matrice, divergence, sum_to_nodes, internal_forces, external_forces) : les mailles d’un même groupe s’écrivent en parallèle sans conflit. Elle est calculée au premier de ces calculs, puis resert à tous les suivants sur ce sous-espace.

C’est une raison de plus de réutiliser un même espace : un second espace sur le même maillage refait ses tables et sa coloration.

Visualisation

champ.plot(component="E") trace le champ sur son propre support : chaque zone retrouve sa zone de maillage par son sous-espace, sans copie. Le rendu reste par maille : dans chaque maille, des valeurs aux nœuds sont ajustées aux valeurs de Gauss par moindres carrés, sans moyenne avec les mailles voisines, si bien que les sauts d’une maille à l’autre (contraintes, flux) restent visibles. smooth règle la subdivision des mailles (0 : une couleur par maille). Voir Visualisation.

Sauvegarde

pyrucast.save écrit le champ avec son espace, son maillage et son Coords. Des champs sauvegardés dans un même fichier partagent toujours leur espace à la relecture et se combinent point à point ; un champ relu ne partage pas son espace avec ceux du script qui l’a écrit, ni avec ceux d’un autre fichier. Le scellement de la zone de maillage est relu avec elle ; la coloration, elle, n’est pas écrite et se recalcule au premier usage.

Côté Rust

SubElementField porte un Handle<SubFiniteElementSpace>, ses noms de composantes, ses deux dimensions et un Vec<f64> à plat, rangé maille → point de Gauss → composante :

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

Les composantes d’un point sont donc contiguës : point_values(maille, g) en rend la tranche sans copie, et c’est ce que lit un noyau de calcul. ElementField est l’agrégat de Handle<SubElementField> ; son finalize refuse une composante portée deux fois sur un même sous-espace.

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

• Contrat commun. Les noms de composantes, min/max/sum, les opérations par composante et l’arithmétique viennent des traits SubField et Field (src/containers/field.rs), partagés avec le champ aux nœuds.
• Support. SubField::same_support compare les Handle par identité (same_object) : c’est ce qui fait de deux espaces sur un même maillage deux supports distincts.
• Durée de vie. Le Handle de l’espace est un Arc : la zone le garde en vie, l’espace garde son SubMesh, qui garde son Coords et compte ses nœuds (cf. Modèle mémoire).
• Coloration. Un OnceLock sur SubFiniteElementSpace, rempli par coloring(…) et marqué #[serde(skip)].

Modèle physique (Model/SubModel)

Un modèle décrit le problème à résoudre : quelles physiques, sur quelles zones du maillage, et quelles contraintes les tiennent. Il ne porte aucune valeur de matériau : la conductivité, le module d’Young ou la densité d’un chargement arrivent dans un champ matériau au moment de l’assemblage, et un même modèle sert avec d’autres valeurs. C’est un agrégat :

• un sous-modèle (SubModel) est une physique sur une zone, c’est-à-dire sur un sous-espace d’un espace éléments finis, ou bien une contrainte ;
• un modèle (Model) est une liste de sous-modèles, de physiques éventuellement différentes : une conduction et une élasticité sur le même maillage, avec leurs appuis.

Le modèle est la troisième couche, posée sur les deux autres :

Déclarer un modèle

Chaque physique se déclare par une fonction de pyrucast.model, qui reçoit l’espace éléments finis et rend un modèle couvrant tout l’espace : un sous-modèle par sous-espace. Une contrainte reçoit, au lieu d’un espace, les maillages qu’elle relie, et rend un seul sous-modèle, quel que soit leur nombre de zones.

On compose les modèles avec |, qui rend un nouveau modèle : les sous-modèles des deux côtés, dans l’ordre, sans répéter un sous-modèle déjà présent. Model() est un modèle vide ; add_sub et add_subs ajoutent des sous-modèles en place. Tous les sous-modèles d’un modèle vivent dans un même Coords : en composer deux de Coords différents est refusé.

Certains termes agissent sur les inconnues d’une autre physique : les contraintes (dirichlet, mpc, embedded, contact), les échanges (boundary_transfer, interface_transfer), le rayonnement (radiation) et le chargement réparti (flux). Ils reçoivent cette physique, leur cible : elle donne la ligne où écrire leur effet (q pour une contrainte sur T, f_x pour une contrainte sur u_x), et garantit qu’une physique du problème assemble bien cette ligne.

import pyrucast

# Un carré de 2 × 2 quadrangles.
c = pyrucast.Coords(dim=2)
a = c.add_node([0.0, 0.0])
b = c.add_node([1.0, 0.0])
carre = pyrucast.mesh.extrude(pyrucast.mesh.line(a, b, 2), [0.0, 1.0], 2)
fes = pyrucast.FiniteElementSpace(carre)

# Chaque déclaration rend un modèle : un sous-modèle par sous-espace de `fes`.
thermique = pyrucast.model.heat_conduction(fes)
solide = pyrucast.model.elasticity(fes, "plane_stress")
# Une contrainte reçoit le modèle dont elle fige une inconnue : sa cible.
appui = pyrucast.model.dirichlet(pyrucast.mesh.poi1_from_nodes([a]), solide, "u_x")

modele = thermique | solide | appui
print([sub.label() for sub in modele])  # ['HeatConduction', 'Elasticity', 'Dirichlet']
print(modele.primal_vars())  # ['T', 'u_x', 'u_y', 'lambda']
print(modele.dual_vars())  # ['q', 'f_x', 'f_y', 'imposed']
print(modele.dual_of("u_y"))  # f_y

# `modele[i]` est une vue sur un sous-modèle.
sub = modele[1]
print(sub.physics(), sub.material_components())  # ['mechanical'] ['E', 'nu']
print(modele[2].fespace(), modele[2].has_behavior)  # None False
print(len(modele.fespace()))  # 2 : l'appui n'intègre rien

# L'appui a créé son nœud-multiplicateur, à la position du nœud imposé.
m = modele.multiplier_mesh().node(0, 0, 0)
print(m.position(), m == a)  # [0.0, 0.0] False

modele[i] rend une vue sur un sous-modèle. Une vue se lit (variables, natures, sous-espace, matériau attendu) mais ne se modifie pas : on change un modèle en le recomposant.

Les physiques disponibles

Le tableau donne, pour chaque déclaration, les variables qu’elle apporte au modèle. Les équations, le matériau et les exemples sont dans la page de chaque physique ; ce que chacune sait assembler est résumé dans Détails des physiques.

symétrie ("isotropic"), mesure ("linear"), sens ("=") et tol ont une valeur par défaut. Une poutre lit sa configuration (1-D, plan, espace) sur le maillage. Les composantes matériau dépendent de la physique et de ses options : sub.material_components() les donne.

interface_transfer rend un sous-modèle par paire de sous-espaces des deux côtés, pris dans l’ordre. update_contact(modele, déplacement) ne déclare rien : il rend le modèle avec ses contacts appariés de nouveau sur la configuration déplacée (voir Contact).

Inconnues et multiplicateurs

Chaque physique déclare :

• ses inconnues, ou variables primales (primal_vars()) : les composantes de la solution, et les colonnes de la matrice ;
• ses duales (dual_vars()) : les composantes du second membre, et les lignes de la matrice. T a pour duale q, u_x a pour duale f_x.

Une inconnue de la matrice est un couple (nœud, nom). Deux physiques qui emploient le même nom au même nœud partagent l’inconnue : c’est ainsi qu’un échange de bord sur T s’ajoute à la conduction, et que deux zones d’élasticité se raccordent sur leurs nœuds communs. Deux noms différents ne se rencontrent jamais : une conduction et une élasticité sur le même maillage forment, dans la même matrice, deux problèmes que rien ne lie.

Une contrainte impose ses relations par multiplicateurs de Lagrange (voir Contraintes). À sa construction, elle crée dans le même Coords un nœud par relation, le nœud-multiplicateur, à la position d’un nœud de la relation : le nœud imposé d’un Dirichlet, le nœud esclave d’un contact. Ce nœud porte l’inconnue lambda, le multiplicateur, dont l’opposé est la réaction, et la duale imposed, la ligne où s’écrit la valeur imposée. Toutes les contraintes emploient ces deux noms sans se mélanger, puisque chacune a ses propres nœuds. multiplier_mesh() rend ces nœuds, en zones POI1.

Assembler et résoudre

matrix.stiffness(modele, materiaux) rend une seule matrice pour tout le modèle : physiques et contraintes, inconnues et multiplicateurs. Les conditions aux limites sont des sous-modèles comme les autres. Le second membre réunit trois sources :

• les données que porte le modèle, rendues par node_field.external_forces(modele, materiaux) : un chargement réparti (flux, dont la densité phi_<duale> est lue dans le matériau), la partie ambiante d’un échange de bord, le jeu d’un contact. Un modèle qui n’en a pas rend un champ vide ;
• les valeurs imposées par une contrainte, rendues par contrainte.constraint_rhs([(nœud, g), …]) : on désigne chaque relation par un de ses nœuds, et la valeur g est écrite à son nœud-multiplicateur, dans la composante imposed. Le modèle doit porter une seule contrainte ; constraint_rhs_by_index désigne la relation par son rang (voir Contraintes) ;
• ce qu’on écrit soi-même : un champ aux nœuds sur les composantes de dual_vars(), comme une force ou un flux ponctuels.

Deux champs posés sur le même maillage s’additionnent avec +, composante par composante. | réunit des champs qui ne se recouvrent pas, comme les valeurs imposées, portées par les nœuds-multiplicateurs.

# Une barre [0, 1] de 4 SEG2 : conduction, source répartie, T imposée à gauche.
c = pyrucast.Coords(dim=1)
o = c.add_node([0.0])
l = c.add_node([1.0])
barre = pyrucast.mesh.line(o, l, 4)
fes = pyrucast.FiniteElementSpace(barre)
conduction = pyrucast.model.heat_conduction(fes)
source = pyrucast.model.flux(fes, conduction, "q")  # densité lue dans phi_q
gauche = pyrucast.model.dirichlet(pyrucast.mesh.poi1_from_nodes([o]), conduction, "T")
modele = conduction | source | gauche

materiaux = pyrucast.element_field.material_field(modele, [("k", 1.0), ("phi_q", 1.0)])
K = pyrucast.matrix.stiffness(modele, materiaux)

# Les données que porte le modèle : ici la source, ∫ phi_q·N dx.
donnees = pyrucast.node_field.external_forces(modele, materiaux)
# Un flux de 1 entrant en x = 1, écrit à la main sur les nœuds de la barre.
ponctuel = pyrucast.NodeField(barre, ["q"])
ponctuel.unit().set_value(l, "q", 1.0)
# La valeur imposée par la contrainte : T = 2 en x = 0.
imposees = gauche.constraint_rhs([(o, 2.0)])

T = pyrucast.solver.solve(K, (donnees + ponctuel) | imposees)
print(T.value(l, "T"))  # 3.5 = 2 + 1 + 1/2
multiplicateur = gauche.multiplier_mesh().node(0, 0, 0)
print(T.value(multiplicateur, "lambda"))  # 2.0 : la réaction de l'appui

Les opérateurs sont détaillés dans Assemblage et Solveur. Des exemples complets à solution analytique sont déroulés dans Dirichlet (Poisson 1-D), Conduction thermique et Mécanique.

Une physique qui n’a pas un genre de matrice n’y contribue pas : matrix.mass sur une conduction avec échange de bord ne prend que la capacité de la conduction.

Nature physique et filtrage

Chaque sous-modèle déclare un ensemble de natures (physics()), qui dit à quelle physique il appartient :

Un flux et les échanges prennent la nature de leur cible : sur une conduction, ils sont "thermal", sur une diffusion, "diffusion". Une contrainte porte "constraint" et la nature de la physique dont elle fige les inconnues : un Dirichlet sur u_x est ["constraint", "mechanical"]. Une relation qui lie les inconnues de deux physiques différentes ne porte que "constraint" : aucun problème d’une seule physique ne saurait la résoudre.

filter(nature) garde les sous-modèles qui portent cette nature, et exclude(nature) ceux qui ne la portent pas. Les deux se composent, et les sous-modèles sont partagés, pas copiés. Comme une contrainte porte la nature de ce qu’elle fige, filter("mechanical") rend la mécanique avec ses appuis, un problème qui se résout tel quel, et filter("mechanical").exclude("constraint") la mécanique seule :

# L'appui fige un déplacement : il est une contrainte, et mécanique aussi.
print(appui.unit().physics())  # ['constraint', 'mechanical']
mecanique = modele.filter("mechanical")  # l'élasticité et son appui
print([sub.label() for sub in mecanique])  # ['Elasticity', 'Dirichlet']
seule = mecanique.exclude("constraint")  # l'élasticité seule
print(seule.primal_vars())  # ['u_x', 'u_y']

Les natures disent ce qu’un sous-modèle fige, pas ce qu’il est : un contact et un Dirichlet sur les déplacements portent tous deux ["constraint", "mechanical"]. label() nomme la physique ("Dirichlet", "Contact"…) ; c’est par lui que le pas à pas reconnaît le contact.

Chaque bloc d’une matrice assemblée garde les natures du sous-modèle qui l’a produit : une matrice se sélectionne de la même façon (voir Matrice).

Ce qui est vérifié, et quand

À la déclaration :

• l’espace convient à la physique : une poutre d’Euler-Bernoulli demande un espace "HERMITE3", une poutre de Timoshenko un espace "MODEL_EMBEDDED" ;
• la variable qu’une contrainte fige est une inconnue de sa cible ;
• la composition refuse un sous-modèle d’un autre Coords.

À l’assemblage : le champ matériau porte, sur le sous-espace de chaque sous-modèle, les composantes qu’il demande. Le matériau n’est connu qu’à ce moment-là.

Ce qui est mis en cache

Un modèle garde, pour chaque genre de matrice (raideur, masse, raideur géométrique, tangente), le motif de sa matrice : quelles lignes et quelles colonnes portent un terme, et à quelle place de la matrice creuse chaque terme élémentaire vient s’ajouter. Il est calculé au premier assemblage de ce genre, puis repris tel quel : les assemblages suivants ne calculent que les valeurs.

Ce motif ne dépend que des inconnues que relient les sous-modèles, pas du matériau, de la position des nœuds ou de l’état mécanique. Un calcul non linéaire, qui réassemble sa tangente à chaque itération, le réutilise donc à chaque fois. Le calculer représente plus de la moitié du premier assemblage : sur un million de quadrangles en élasticité plane (4 cœurs), le premier assemblage prend 1,9 s, les suivants 0,7 s.

Le motif est oublié quand on ajoute un sous-modèle (add_sub, add_subs). Un nouveau modèle, rendu par |, filter, exclude, une tranche modele[i:j] ou update_contact, part sans motif : dans une boucle, on garde le modèle dans une variable plutôt que de le recomposer à chaque tour.

Sauvegarde

Un modèle se sauvegarde avec tout ce qu’il désigne : ses espaces, leurs maillages, le Coords, les nœuds-multiplicateurs de ses contraintes. Sauvegardés dans un même fichier, un modèle et son maillage restent liés à la lecture : le modèle relu est posé sur le maillage relu. Le motif des matrices n’est pas écrit ; il se recalcule au premier assemblage.

Côté Rust

Model est l’agrégat de Handle<SubModel>. Il garde aussi le motif de chaque genre de matrice, dans matrix_patterns, un tableau de OnceLock<Arc<AssemblyPattern>> à une case par MatrixKind. Son check_push refuse un sous-modèle d’un autre Coords, et son post_push vide les motifs. ops::matrix::assemble_kind demande ses blocs à chaque sous-modèle, puis prend le motif du modèle, construit au premier appel par ops::scatter::build_pattern, et y répartit les valeurs en parallèle.

SubModel est une enum à une variante par physique, chacune enveloppant la struct de src/models/ qui porte ses données. L’enum sert à ranger et à sauvegarder : bincode écrit le rang de la variante, donc une nouvelle variante se place à la fin. Le code commun (l’agrégat, l’assemblage, l’affichage) passe par as_kind(), qui rend la struct de la variante en &dyn SubModelKind.

SubModelKind (src/models/mod.rs) est ce que toute physique déclare : ses variables, ses natures (physics(), un &'static [Physics]), son nom (label()), ses blocs pour chaque genre de matrice, ses termes de forces internes et externes, son affichage. Trois traits la complètent selon la physique, rendus par as_domain, as_constraint et as_behavior (None par défaut) : Domain pour le matériau et les matrices élémentaires, Constraint pour les nœuds-multiplicateurs et les relations, Behavior pour la loi de comportement. Le chapitre Ajouter une physique détaille ces traits. Un bloc de matrice monté à la main reçoit ses natures par SubMatrix::set_physics, que Python n’expose pas.

#[test]
fn un_modele_se_declare_et_s_assemble() -> Result<()> {
    // Une barre [0, 1] d'un seul 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 maillage = Mesh::from_submesh(SubMesh::new(coords.clone(), ElementType::SEG2));
    maillage.add_cell(&[a.id(), b.id()])?;
    let fes = FiniteElementSpace::lagrange1(&maillage)?;

    // La conduction, puis l'appui qui fige T en `a` ; `union` est le `|` de
    // Python.
    let conduction = model::heat_conduction(&fes, MaterialSymmetry::Isotropic)?;
    let impose = mesh::poi1_from_nodes(std::slice::from_ref(&a))?;
    let appui = model::dirichlet(&impose, &conduction, "T", RelationSense::Equality)?;
    let modele = conduction.union(&appui)?;

    // Le matériau arrive à l'assemblage ; l'appui n'en demande pas.
    let materiaux = element_field::material_field(&modele, &[("k", 1.0)])?;
    let k = matrix::stiffness(&modele, &materiaux)?;
    assert_eq!(k.n_rows()?, 3); // 2 nœuds + 1 nœud-multiplicateur

    // L'appui fige une température : il est thermique aussi.
    let thermique = modele.filter(Physics::Thermal);
    let seule = thermique.exclude(Physics::Constraint);
    assert_eq!((thermique.len(), seule.len()), (2, 1));
    Ok(())
}

Matrice creuse (Matrix/SubMatrix)

Une matrice est ce que rend un assemblage : pyrucast.matrix.stiffness(modele, materiaux) rend une seule matrice pour tout le modèle, physiques et contraintes comprises, et mass, geometric et tangent font de même (voir Assemblage). Elle est creuse : seuls ses termes non nuls sont gardés. C’est un agrégat :

• un bloc (SubMatrix) est la part d’un sous-modèle : la conduction d’une zone, ou l’une des deux moitiés d’une contrainte ;
• une matrice (Matrix) est la liste de ses blocs, plus ses valeurs assemblées : la matrice que lit le solveur.

Lire une matrice

Une ligne est un couple (nœud, duale), une colonne un couple (nœud, inconnue) : la ligne (0, "q") est le flux au nœud 0, la colonne (0, "T") la température au même nœud. row_dofs() et col_dofs() les donnent dans l’ordre de la matrice, chaque nœud par son identifiant.

import pyrucast

# Une barre [0, 1] de deux SEG2 : conduction, T imposée à gauche.
c = pyrucast.Coords(dim=1)
o = c.add_node([0.0])
l = c.add_node([1.0])
barre = pyrucast.mesh.line(o, l, 2)  # crée le nœud du milieu, d'id 2
fes = pyrucast.FiniteElementSpace(barre)
conduction = pyrucast.model.heat_conduction(fes)
appui = pyrucast.model.dirichlet(pyrucast.mesh.poi1_from_nodes([o]), conduction, "T")
modele = conduction | appui
materiaux = pyrucast.element_field.material_field(
    modele, [("k", 1.0), ("rho", 1.0), ("cp", 1.0)]
)

K = pyrucast.matrix.stiffness(modele, materiaux)
print(
    K
)  # Matrix: 3 sub-matrix(es) [#…, #…, #…], 4 row(s) × 4 col(s), symmetric, ~296 B
print(K.row_dofs())  # [(0, 'q'), (2, 'q'), (1, 'q'), (3, 'imposed')]
print(K.col_dofs())  # [(0, 'T'), (2, 'T'), (1, 'T'), (3, 'lambda')]
print(K.get(o, "q", o, "T"))  # 2.0
K.dump()
#               (0,T)   (2,T)   (1,T)  (3,lambda)
#       (0,q)   2.000  -2.000   0.000       1.000
#       (2,q)  -2.000   4.000  -2.000       0.000
#       (1,q)   0.000  -2.000   2.000       0.000
# (3,imposed)   1.000   0.000   0.000       0.000

Les lignes et les colonnes sont rangées dans l’ordre où les blocs les présentent, ou selon la permutation du Coords s’il en porte une. Dans une matrice symétrique, la ligne i et la colonne i se répondent : même nœud, la duale en ligne et son inconnue en colonne (q et T, imposed et lambda). Une matrice peut aussi être rectangulaire : un bloc de contrainte a une ligne par relation et une colonne par nœud que les relations relient.

get rend 0 pour une case que la matrice ne contient pas. entries() rend toutes les valeurs, ligne par ligne, chacune en tuple (nœud, duale, nœud, inconnue, valeur). dump() affiche la matrice en tableau, une lecture à réserver aux petites matrices.

Les blocs

K[i] est une vue sur un bloc. Un bloc a deux supports, les nœuds de ses lignes et ceux de ses colonnes, en zones POI1, et ses noms de lignes et de colonnes. La matrice ci-dessus a trois blocs : la conduction, puis les deux moitiés de l’appui, C (ligne imposed, colonne T) et Cᵀ (ligne q, colonne lambda).

for bloc in K:
    print(bloc.row_dofs(), bloc.col_dofs(), bloc.physics(), bloc.entry_count)
# [(0, 'q'), (2, 'q'), (1, 'q')] [(0, 'T'), (2, 'T'), (1, 'T')] ['thermal'] 0
# [(3, 'imposed')] [(0, 'T')] ['constraint', 'thermal'] 1
# [(0, 'q')] [(3, 'lambda')] ['constraint', 'thermal'] 1

Un bloc est de l’une de deux sortes :

• calculé : il ne garde aucune valeur, seulement de quoi les calculer, son sous-modèle, son sous-espace et son matériau. Ses valeurs sont calculées à l’assemblage et écrites directement dans la matrice. Les physiques produisent des blocs calculés : la conduction n’a aucune valeur écrite (entry_count) ;
• écrit : il garde ses valeurs, une par une. Ce sont les blocs des contraintes, et ceux qu’on écrit soi-même.

Les valeurs se lisent donc sur la matrice, pas sur ses blocs : sur un bloc calculé, get et entries lèvent une erreur. entry_count, sur un bloc comme sur la matrice, compte les valeurs écrites dans les blocs.

Pour écrire un bloc soi-même, Matrix.block rend une matrice d’un seul bloc, vide, qu’on remplit par add_entry. Deux valeurs écrites dans la même case s’additionnent. assemble() donne ensuite ses valeurs à la matrice.

Un bloc est scellé dès qu’une matrice qui le contient reçoit ses valeurs : par assemble(), par un assembleur ou par une somme. add_entry est alors refusé, puisque les valeurs de la matrice ont été calculées à partir de ses entrées. On écrit donc un bloc avant d’assembler, comme on remplit une zone de maillage avant qu’un espace ou un champ ne la scelle.

noeuds = barre.to_poi1()  # les trois nœuds de la barre, chacun une fois
B = pyrucast.Matrix.block(noeuds, noeuds, ["q"], ["T"], symmetry="full")
b = B.unit()
b.add_entry(o, "q", o, "T", 1.0)
b.add_entry(l, "q", l, "T", 1.0)
b.add_entry(o, "q", o, "T", 0.5)  # s'ajoute à la précédente
B.assemble()
print(B.get(o, "q", o, "T"))  # 1.5
print(b.memory_bytes)  # 72 : trois valeurs écrites, 24 octets chacune

Les nœuds d’un support doivent être distincts, comme ceux que rend to_poi1 : un support qui répète un nœud est refusé, comme pour un champ aux nœuds. Quand un bloc porte plusieurs noms, ordering dit comment ses lignes se suivent : nœud par nœud ("nodes_then_vars", par défaut) ou nom par nom ("vars_then_nodes"). symmetry vaut "full" ou "none" (par défaut), voir Symétrie.

Assemblée ou non

Une matrice a des valeurs une fois assemblée. Les assembleurs rendent une matrice assemblée, sur la configuration active au moment de l’appel. Ce qui réunit des blocs rend une matrice sans valeurs : |, filter, exclude, une tranche K[i:j]. add_sub et add_subs, qui changent la matrice en place, effacent ses valeurs.

assemble() calcule les valeurs à partir des blocs seuls, sans modèle, sur la configuration active : un bloc calculé refait son calcul, un bloc écrit apporte ses valeurs.

U = K.filter("thermal").exclude("constraint")  # la conduction sans son appui
print(U)  # Matrix: 1 sub-matrix(es) [#…], 3 row(s) × 3 col(s), symmetric, ~0 B
try:
    U.get(o, "q", o, "T")
except RuntimeError as erreur:
    print(
        erreur
    )  # Matrix::get: this matrix carries a computed block that was never assembled…
U.assemble()  # calcule les valeurs, sur la configuration active
print(U.get(o, "q", o, "T"))  # 2.0

Sur une matrice sans valeurs, K * x et solve lèvent une erreur. get et entries lisent les valeurs écrites dans les blocs, et lèvent une erreur dès qu’un bloc est calculé.

Une matrice est assemblée une fois : ses valeurs ne bougent plus, et un second assemble() ne fait rien. Les valeurs d’une autre configuration s’obtiennent en refaisant la matrice.

Calculer avec des matrices

• 2.0 * K, K * 2.0, K / dt et -K rendent une matrice neuve, assemblée si K l’était, aux valeurs multipliées. Ses blocs sont des copies qui portent le facteur (factor) : K n’est pas touchée. Diviser par zéro lève ZeroDivisionError, par un nombre infini ou nan, ValueError.
• A + B et A - B additionnent les valeurs telles qu’elles ont été assemblées, sur la réunion des lignes et des colonnes des deux côtés : le résultat est assemblé et se résout tel quel. Rien n’est recalculé, si bien que deux matrices assemblées sur deux configurations différentes s’additionnent telles quelles. Un opérande dont un bloc calculé n’a jamais été assemblé (une sortie de | ou de filter) lève une erreur ; un opérande fait de blocs écrits seulement est assemblé au passage. Chaque côté peut être une matrice ou un bloc.
• Les deux côtés n’ont pas besoin des mêmes lignes. Dans M / dt + K, l’appui n’entre que dans K, et M ne contribue rien à ses lignes. A + B range d’abord les lignes de A, puis celles que seule B apporte.
• K | K est K : les mêmes blocs une seule fois, sans valeurs. K + K vaut 2K. On compose une matrice de morceaux distincts avec |, on additionne deux matrices avec +.
• K * x, où x est un champ aux nœuds, est le produit : x est lu aux colonnes (une valeur absente vaut 0), et le résultat est un champ posé sur les lignes.
• lump() rend la matrice diagonale dont chaque terme est la somme de sa ligne : la masse concentrée. Elle garde les lignes et les colonnes de la matrice d’origine, avec un bloc diagonal par support de lignes, qui porte ses noms et ses natures : elle se résout et se multiplie comme elle.

Deux matrices ne se multiplient pas entre elles.

M = pyrucast.matrix.mass(modele, materiaux)  # la capacité : l'appui n'y entre pas
dt = 0.1
A = M / dt + K  # assemblée, prête à résoudre
print(A.row_dofs() == K.row_dofs())  # True

print((K + K).get(o, "q", o, "T"))  # 4.0 : chaque opérande compte
une_fois = K | K  # les mêmes blocs, une seule fois, sans valeurs
une_fois.assemble()
print(une_fois.get(o, "q", o, "T"))  # 2.0

# Le produit par un champ lit T aux colonnes et rend q aux lignes.
x = pyrucast.NodeField(barre, ["T"])
x.unit().set_uniform("T", 1.0)
print((K * x).value(l, "q"))  # 0.0 : une température uniforme ne crée aucun flux

Natures physiques

Chaque bloc garde les natures du sous-modèle qui l’a produit (voir Modèle), et K.physics() réunit celles de tous les blocs. filter(nature) garde les blocs qui portent cette nature, exclude(nature) ceux qui ne la portent pas ; les blocs sont partagés, et le résultat est sans valeurs. Un bloc écrit à la main n’a pas de nature : filter ne le garde jamais, exclude le garde toujours.

Symétrie

K.symmetric dit si la matrice assemblée est symétrique. Les valeurs ne sont pas comparées : chaque bloc déclare sa part, et la matrice réunit les déclarations.

• Un bloc de raideur élastique ou de masse est symétrique à lui seul (is_symmetric).
• Les blocs C et Cᵀ d’une contrainte ne le sont pas seuls, mais forment ensemble une symétrie : la matrice n’est symétrique que s’ils sont présents tous les deux.
• Un bloc écrit à la main déclare symmetry="full" ou "none".

print([bloc.is_symmetric for bloc in K])  # [True, False, False]
print(K.symmetric)  # True : C et Cᵀ sont présents tous les deux
print(K[0:2].symmetric)  # False : Cᵀ a été laissé de côté

Le solveur lit la matrice colonne par colonne. Une matrice symétrique s’y prête telle quelle, une autre est d’abord recopiée dans ce sens. Seule une matrice symétrique accepte method="cholesky". Une déclaration fausse n’est pas détectée : déclarer "full" un bloc qui ne l’est pas fait résoudre le système de la matrice transposée.

Résoudre

K.solve(second_membre) est pyrucast.solver.solve(K, second_membre). La matrice porte tout ce qu’il faut : les lignes des relations, leurs coefficients et leur sens. Les options sont décrites dans Solveur.

impose = appui.constraint_rhs([(o, 1.0)])  # T = 1 en x = 0
T = K.solve(impose)  # factorise K, et garde la factorisation dans K
print(T.value(l, "T"))  # 1.0
T = K.solve(2.0 * impose)  # reprend la factorisation : la descente-remontée seule
print(T.value(l, "T"))  # 2.0

Ce qui est mis en cache

Le premier solve sur une matrice la factorise, c’est-à-dire la décompose en un produit de matrices triangulaires, l’étape coûteuse, et garde le résultat dans la matrice. Les solve suivants sur la même matrice le reprennent et ne font que la descente-remontée, bien moins chère : plusieurs chargements, un transitoire à pas de temps constant, des itérations qui ne changent pas la matrice.

La matrice n’a qu’une case pour cela. Son contenu dépend de la voie du solveur (voir Solveur) :

• une factorisation de la matrice, LU ou Cholesky, quand on résout avec les multiplicateurs sans inégalité, ou par élimination une matrice sans relation ;
• une condensation, pour l’élimination (relations="eliminate") : la table T qui exprime chaque inconnue éliminée par les autres, la matrice sans ses contraintes, et la factorisation du système réduit K̂ = TᵀKT ; ou, si des inégalités restent, ce système réduit et l’état de l’ensemble actif qui le résout ;
• l’état de l’ensemble actif, quand la matrice porte des inégalités : les inégalités, le dernier statut convergé (d’où repart la résolution suivante) et ses factorisations LU — celle de la base de Schur et ses colonnes, ou celle du dernier statut.

Le contenu n’est repris que s’il répond à la demande : la même voie, et, pour une factorisation, la même method. Sinon il est recalculé et prend la place.

Ce qui est gardé est oublié quand la matrice change (add_sub, add_subs). Une matrice neuve part sans : 2.0 * K, M / dt + K ou un nouvel assemblage factorisent de nouveau. Dans une boucle, on construit la matrice une fois et on la garde dans une variable. solve(…, cache=False) ne garde rien.

Une matrice rendue par un assembleur partage sa structure, c’est-à-dire quelles cases sont non nulles et l’ordre des lignes et des colonnes, avec le motif que le modèle garde pour ce genre de matrice (voir Modèle) : deux raideurs du même modèle n’ont qu’une structure en mémoire, et chacune ne possède que ses valeurs.

Ce que pèse une matrice

memory_bytes estime les octets occupés. Un bloc compte 24 octets par valeur écrite (deux indices et la valeur), et 0 s’il est calculé. Une matrice compte ses valeurs assemblées et leur structure (un indice de colonne et une valeur par terme non nul, un début par ligne, une clé par ligne et par colonne), plus ce que gardent ses blocs. Une structure partagée entre matrices est comptée dans chacune. L’estimation apparaît dans l’affichage (~296 B plus haut).

La factorisation n’est pas comptée, alors qu’elle est en général bien plus grosse que la matrice : faer, la bibliothèque qui la calcule, ne donne pas la taille de ses facteurs. Pour la mémoire d’un solve, voir Calculs plus gros que la RAM.

Sauvegarde

Une matrice se sauvegarde avec ses blocs, et avec ce que désignent ses blocs calculés : sous-modèle, sous-espace, matériau. Ni ses valeurs assemblées ni sa factorisation ne sont écrites : à la relecture, assemble() recalcule les valeurs depuis les blocs, sur la configuration active. Une somme de matrices assemblées sur deux configurations est donc recalculée sur une seule.

Côté Rust

Matrix est l’agrégat de Handle<SubMatrix>, avec deux membres de plus :

• assembled, les valeurs assemblées et la numérotation des lignes et des colonnes, None tant que la matrice n’est pas assemblée. Les valeurs sont rangées ligne par ligne (forme CSR) : trois tableaux, le début de chaque ligne, la colonne de chaque terme et sa valeur. Pour une matrice rendue par un assembleur, les deux premiers et la numérotation sont des Arc partagés avec le motif du modèle (AssemblyPattern) ; seules les valeurs lui appartiennent.
• factorization, ce que garde le solveur, dans un parking_lot::Mutex<Option<Arc<dyn Any + Send + Sync>>>. dyn Any accepte n’importe quel type, si bien que containers ne dépend pas du solveur : cached_factorization::<T>() le rend s’il est de type T, et store_factorization le range. Le Mutex permet de le remplir depuis solve(&Matrix), qui ne reçoit qu’une référence partagée.

check_push refuse un bloc posé sur un autre Coords, et post_push vide assembled et factorization.

Une ligne ou une colonne est une DofKey : un u64 qui porte l’identifiant du nœud dans ses 32 bits hauts, et dans ses 32 bits bas le rang de son nom dans la table de noms de la matrice (dof_vars()). row_dofs() rend la forme lisible, NamedDof = (NodeId, String).

SubMatrix porte :

• row_support et col_support, deux Handle<SubMesh> POI1 scellés à la construction du bloc. Le bloc ne recopie pas leur liste de nœuds : il la lit en place et trouve le rang d’un nœud par la table node_index du support (voir Maillage). D’où la règle des nœuds distincts : cette table donne le rang parmi les nœuds sans répétition, et les constructeurs (new, computed, from_coo) refusent un support qui répète un nœud ou qui n’est pas un POI1 ;
• dual_vars, primal_vars et ordering (DofOrdering) : la ligne d’un couple (nœud, nom) est rang_du_nœud * n_noms + rang_du_nom en NodesThenVars, rang_du_nom * n_nœuds + rang_du_nœud en VarsThenNodes ;
• coo, les valeurs d’un bloc écrit, en liste de triplets (ligne, colonne, valeur) (nalgebra_sparse::CooMatrix), vide pour un bloc calculé ;
• recipe, Some(ComputedRecipe) pour un bloc calculé : le sous-modèle, les sous-espaces, le matériau, le MatrixKind et les KernelInputs, ce que le calcul lit en plus du matériau (la contrainte pour la raideur géométrique, la déformation et l’état précédent pour la tangente) ;
• symmetry, physics, factor, et row_sense, le sens (égalité ou inégalité) des lignes imposed, que lit le solveur ;
• sealed, posé quand une matrice qui contient le bloc reçoit ses valeurs (assemblage ou somme), et lu par add_entry. Il se pose comme celui d’une zone de maillage : la fonction seal lit d’abord, et ne prend le verrou d’écriture que si le bloc n’est pas encore scellé. Il n’est pas sauvegardé : une matrice relue n’a pas de valeurs, et son assemblage scelle de nouveau ses blocs.

Sur un bloc calculé, toute lecture de valeurs (get, iter_entries, to_dmatrix, to_coo, to_csr, to_csc) rend une erreur.

Assemblage. ops::matrix::stiffness, comme mass, geometric et tangent, demande ses blocs à chaque sous-modèle, prend le motif du modèle (construit au premier appel par ops::scatter::build_pattern), calcule les valeurs en parallèle par coloriage des mailles (voir Parallélisme) et les remet à la matrice par set_assembled. Matrix::assemble est défini dans ops::matrix et non avec la matrice : calculer un bloc calculé demande le calcul élémentaire, qui vit dans models, dont containers ne dépend pas. Il construit le motif depuis les blocs à chaque appel.

Facteur. Il s’applique partout où une valeur du bloc est lue ou émise (get, iter_entries, les conversions, l’assemblage), sauf dans local_triplets et local_coo_arrays, les valeurs brutes que chaque lecteur multiplie lui-même. &k * s copie chaque bloc ; k * s, qui prend k, multiplie en place les blocs que personne d’autre ne tient (Handle::is_sole_owner). Les valeurs assemblées sont multipliées aussi, la structure restant partagée : pas de réassemblage. (Σ v) · s n’est pas, au bit près, le Σ (v · s) d’un réassemblage, l’addition des nombres à virgule flottante n’étant pas associative. En Rust, une division par zéro ou par un nombre non fini arrête le programme (panic).

Somme. &a + &b rend un Result<Matrix>. Deux matrices assemblées sur le même motif additionnent leurs valeurs case à case ; sinon, les lignes et les colonnes sont renumérotées sur leur réunion.

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

    // `/` ne peut pas échouer et rend une `Matrix` ; `+` rend un `Result`,
    // puisqu'un opérande peut ne pas avoir de valeurs.
    let systeme = (&m / dt + &k)?;
    let t = solver::solve(&systeme, &rhs)?;

    assert!(t.node_count() > 0);
    Ok(())
}

Symétrie. Symmetry vaut Full, Half(PairId) ou None. Les deux moitiés d’une paire reçoivent ensemble leur identité par mint_pair, calculée sur leur contenu (nœuds des deux supports, noms, coefficient) : elle ne dépend pas de l’exécution, et deux sauvegardes des mêmes objets donnent les mêmes octets. Deux paires distinctes qui tomberaient sur la même identité seraient rejetées, jamais acceptées à tort. Matrix::symmetric() est vrai si chaque bloc est Full, ou Half avec ses deux membres présents en nombres égaux : une contrainte déclarée deux fois reste symétrique.

Dans une matrice symétrique, lignes et colonnes sont numérotées en un seul parcours, pour que la ligne i et la colonne i se répondent : la ligne k d’un bloc Full fait face à sa colonne k, et une paire Half se croise. Deux incohérences y sont refusées : un même degré de liberté dual associé à deux primaux différents, et une paire dont les membres n’ont pas la même taille. Deux blocs rectangulaires qui sont chacun une tranche de lignes d’une matrice symétrique ne peuvent rien déclarer : la matrice répond false, et le solveur prend le chemin général.

Solveur. csr_arrays() prête les trois tableaux sans copie. La factorisation de faer lit la matrice par colonnes (forme CSC) ; or les tableaux CSR d’une matrice sont exactement les tableaux CSC de sa transposée. Une matrice symétrique est donc passée telle quelle, une autre est transposée une fois. to_csr, to_csc, to_coo et to_dmatrix fabriquent des objets nalgebra neufs.

#[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()])?;

    // Une valeur.
    let v: f64 = k.get(a, "q", a, "T")?;

    // Les trois tableaux de la matrice creuse, prêtés sans copie…
    let (debuts, colonnes, valeurs): (&[usize], &[usize], &[f64]) = k.csr_arrays()?;
    assert_eq!(debuts.len(), k.n_rows()? + 1);
    assert_eq!(colonnes.len(), valeurs.len());
    // …ou des objets neufs, de nalgebra et nalgebra-sparse.
    let plein: nalgebra::DMatrix<f64> = k.to_dmatrix()?;
    let csr: nalgebra_sparse::CsrMatrix<f64> = k.to_csr()?;
    let csc: nalgebra_sparse::CscMatrix<f64> = k.to_csc()?;

    // Le produit par un champ : `x` est lu aux colonnes, le résultat est posé
    // sur les lignes.
    let y: NodeField = (&k * &x)?;

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

Un bloc écrit se construit par SubMatrix::new, puis add_entry ; SubMatrix::set_physics, que Python n’expose pas, lui donne des natures, par exemple Physics::Other pour que filter puisse le retrouver.

#[test]
fn ecrire_un_bloc() -> 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)
    };

    let mut bloc = SubMatrix::new(
        support.clone(),  // support des lignes
        support.clone(),  // support des colonnes
        vec!["q".into()], // variables duales : les lignes
        vec!["T".into()], // variables primales : les colonnes
        DofOrdering::NodesThenVars,
        Symmetry::Full, // symétrique à lui seul
    )
    .unwrap();
    // K = [[2, -1], [-1, 2]]
    bloc.add_entry(a.id(), "q", a.id(), "T", 2.0)?;
    bloc.add_entry(a.id(), "q", b.id(), "T", -1.0)?;
    bloc.add_entry(b.id(), "q", a.id(), "T", -1.0)?;
    bloc.add_entry(b.id(), "q", b.id(), "T", 2.0)?;

    let mut k = Matrix::empty();
    k.add_sub(Handle::new(bloc))?;
    k.assemble()?; // assemble les blocs écrits

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

Évolution (Evolution/SubEvolution)

Une évolution associe des valeurs à des abscisses et interpole linéairement entre elles. L’abscisse est souvent le temps, mais pas forcément : une température, une déformation, ou plusieurs variables à la fois. La valeur est un nombre, un champ aux nœuds ou un champ aux points de Gauss.

C’est un agrégat :

• une SubEvolution est une courbe : la liste de ses abscisses et la liste des valeurs en regard. L’interpoler en x rend une valeur ;
• une Evolution réunit des SubEvolution, une par zone, comme un NodeField réunit ses zones. L’interpoler interpole chaque courbe et regroupe les résultats : un NodeField ou un ElementField pour des champs, une liste de flottants pour des nombres.

Une courbe

SubEvolution(échantillons) prend une liste de couples (abscisse, valeur). Les échantillons sont triés à la construction ; deux à la même abscisse sont refusés. Entre deux échantillons x₀ ≤ x ≤ x₁, la valeur est v₀ + (v₁ − v₀)·(x − x₀)/(x₁ − x₀) ; sur un échantillon, c’est sa valeur.

import pyrucast

# Une température qui monte de 20 à 120 °C en 10 s, puis reste à 120 °C.
montee = pyrucast.SubEvolution([(0.0, 20.0), (10.0, 120.0), (60.0, 120.0)])
print(montee.interpolate(2.5))  # 45.0
print(montee.interpolate(10.0))  # 120.0 : un échantillon, rendu tel quel
print(montee.abscissas())  # [0.0, 10.0, 60.0]
print(len(montee))  # 3

Hors plage

Une abscisse en dehors de la plage tabulée suit la politique hors plage de l’évolution :

La politique se donne à la construction (out_of_range="clamp"), se lit par out_of_range, et s’impose pour un seul appel : interpolate(x, out_of_range="clamp").

try:
    montee.interpolate(70.0)  # politique "error", par défaut
except RuntimeError as erreur:
    print(erreur)  # evolution: abscissa 70 is outside the tabulated range [0, 60]
print(montee.interpolate(70.0, out_of_range="clamp"))  # 120.0
print(montee.interpolate(-1.0, out_of_range="extrapolate"))  # 10.0

bornee = pyrucast.SubEvolution([(0.0, 20.0), (10.0, 120.0)], out_of_range="clamp")
print(bornee.out_of_range, bornee.interpolate(70.0))  # clamp 120.0

Des champs comme valeurs

Evolution([(t₀, champ₀), (t₁, champ₁), …]) prend un champ entier par abscisse, un NodeField ou un ElementField. Elle en fait une courbe par zone : les zones sont celles du premier champ, et chaque champ suivant doit porter une zone sur chacun de ces supports ; une zone qu’il aurait en plus est ignorée. Toutes les valeurs d’une courbe ont le même support et les mêmes composantes.

L’évolution copie les valeurs des champs ; le support, lui, est partagé. Modifier un champ ensuite ne change pas l’évolution, et le champ interpolé est posé sur le même support que les champs d’origine : il se combine avec eux zone à zone.

shared_abscissas() rend les abscisses, communes à toutes les zones, et frames() les champs tabulés, un par abscisse.

c = pyrucast.Coords(1)
noeuds = [c.add_node([float(i)]) for i in range(3)]
support = pyrucast.mesh.poi1_from_nodes(noeuds)
froid = pyrucast.NodeField(support, ["T"])
froid.unit().set_uniform("T", 20.0)
chaud = pyrucast.NodeField(support, ["T"])
chaud.unit().set_uniform("T", 120.0)

chauffe = pyrucast.Evolution([(0.0, froid), (10.0, chaud)])
t = chauffe.interpolate(2.5)  # un NodeField, sur le même support
print(t.value(noeuds[0], "T"))  # 45.0
print(chauffe.shared_abscissas())  # [0.0, 10.0]
print(len(chauffe.frames()))  # 2 : les champs tabulés, un par abscisse

Une abscisse à plusieurs variables

Une abscisse peut être un tuple de deux à quatre nombres : la courbe devient une table. Chaque échantillon est alors ((T, φ), valeur), et interpolate((T, φ)) rend la valeur en ce point.

# Limite d'élasticité σ_y (MPa) selon la température T et la fraction φ :
# deux valeurs de φ mesurées à 20 °C, trois à 500 °C.
sigma_y = pyrucast.SubEvolution(
    [
        ((20.0, 0.0), 300.0),
        ((20.0, 1.0), 450.0),
        ((500.0, 0.0), 200.0),
        ((500.0, 0.5), 260.0),
        ((500.0, 1.0), 320.0),
    ],
    abscissa_type=("T", "phi"),
    ordinate_type="sigma_y",
)
print(sigma_y.interpolate((300.0, 0.4)))  # 294.666…

La table se lit variable après variable. Les échantillons sont rangés par leur première variable : ici, un groupe à 20 °C (φ = 0 et 1), un autre à 500 °C (φ = 0, 0,5 et 1). Chaque groupe est une courbe le long de la variable suivante. En (300, 0,4) :

   φ
   1  ● 450                               ● 320
      │                                   │
 0,5  │                                   ● 260
 0,4  ◆ 360 ─────────────── ★ ─────────── ◆ 248
      │                   294,7           │
   0  ● 300                               ● 200
      └────────────────────┬──────────────┴───── T
      20                  300            500

  ●  échantillons de la table
  ◆  étape 1 : chaque groupe est lu le long de φ, en φ = 0,4
  ★  étape 2 : les deux résultats sont pesés le long de T, en T = 300

1. On cherche les deux groupes qui encadrent T = 300 : 20 °C et 500 °C.
2. Chacun est interpolé le long de φ en 0,4 : entre 300 et 450 à 20 °C, ce qui donne 360 ; entre 200 et 260 à 500 °C (φ de 0 à 0,5), ce qui donne 248.
3. On interpole entre ces deux résultats le long de T : t = (300 − 20)/(500 − 20) ≈ 0,583, et 360 + 0,583 × (248 − 360) ≈ 294,7.

Pour deux variables, avec T_a ≤ T ≤ T_b les groupes encadrants et f_a, f_b leurs courbes en φ :

\[ f(T, \varphi) = (1-t)\, f_a(\varphi) + t\, f_b(\varphi), \qquad t = \frac{T - T_a}{T_b - T_a}. \]

Avec trois ou quatre variables, chaque groupe se lit de la même façon le long de la variable suivante. Ce que cela donne :

• les groupes n’ont pas à partager leurs abscisses : la table ci-dessus mesure deux valeurs de φ à 20 °C et trois à 500 °C ;
• sur une grille complète, le résultat est l’interpolation multilinéaire (bilinéaire pour deux variables). Une fonction linéaire en chaque variable séparément, comme a + bT + cφ + dTφ, est retrouvée exactement, grille complète ou non ;
• un groupe où la variable suivante ne prend qu’une valeur est constant le long d’elle : une propriété mesurée pour une seule valeur de φ à une température donnée ne dépend pas de φ à cette température ;
• l’ordre des variables compte : la première range les groupes, et c’est le long d’elle que se fait le dernier mélange. Sur une table dont les groupes ne partagent pas leurs abscisses, échanger l’ordre des variables change le résultat ;
• hors de la plage d’une variable, la politique hors plage s’applique à cette variable, dans le groupe où elle est lue.

abscissas() rend les abscisses triées : des nombres pour une courbe, des tuples pour une table, rangés par première variable puis par la suivante.

Types d’abscisse et d’ordonnée

Une courbe peut porter le nom de ce que mesurent ses axes :

• abscissa_type, le nom de l’abscisse ("T", "t"), un tuple de noms pour une table (("T", "phi")), un par variable ;
• ordinate_type, le nom de la valeur ("E"), pour une courbe de nombres seulement : un champ a déjà ses composantes, et le donner à une évolution de champs est une erreur.

Ils se donnent à la construction et se lisent par abscissa_type() et ordinate_type() ; ceux d’une Evolution sont ceux de sa première courbe. Ils étiquettent les tracés et relient une courbe à un champ (section suivante).

young = pyrucast.SubEvolution(
    [(20.0, 210e3), (500.0, 170e3)], abscissa_type="T", ordinate_type="E"
)
print(young.abscissa_type(), young.ordinate_type())  # T E
print(sigma_y.abscissa_type())  # ('T', 'phi')

Appliquer une courbe à un champ

Une courbe de nombres s’applique à un champ comme cos : interpolate reçoit un champ au lieu d’une abscisse, et rend un champ sur le même support, dont chaque nœud ou point de Gauss est la courbe évaluée en ses propres valeurs.

• La courbe lit les composantes nommées comme son abscissa_type, une par variable. Sans abscissa_type, ou si le champ n’a pas ces composantes, c’est une erreur.
• Elle écrit une composante, nommée d’après son ordinate_type, ou "value" à défaut.
• La politique hors plage s’applique point par point.

etat = pyrucast.NodeField(support, ["T", "phi"])
for noeud, (T, phi) in zip(noeuds, [(20.0, 0.0), (260.0, 0.4), (500.0, 1.0)]):
    etat.unit().set_value(noeud, "T", T)
    etat.unit().set_value(noeud, "phi", phi)

E = young.interpolate(etat[0])  # lit la composante "T", écrit "E"
print(E.components(), E.value(noeuds[1], "E"))  # ['E'] 190000.0

Une Evolution applique toutes ses courbes : k courbes donnent un champ à k composantes, chacune lisant ses propres variables. Deux courbes qui écriraient la même composante sont refusées. C’est ainsi qu’on calcule des propriétés matériau qui dépendent de l’état, par exemple E(T) et σ_y(T, φ) aux points de Gauss en un seul appel, avant de les réunir au besoin aux propriétés constantes ; il n’y a pas d’opérateur matériau à part.

lois = young | sigma_y  # deux courbes : E(T) et σ_y(T, φ)
materiau = lois.interpolate(etat)  # chaque nœud, en ses propres T et φ
print(materiau.components())  # ['E', 'sigma_y']
print(materiau.value(noeuds[1], "sigma_y"))  # 304.0

Composer des courbes

| réunit des SubEvolution en une Evolution, une courbe par zone. Toutes portent le même genre de valeur : des nombres, ou des champs aux nœuds, ou des champs aux points de Gauss. evolution[i] rend la courbe i.

Une Evolution interpolée applique sa politique hors plage à toutes ses courbes. Celle que rendent | et une tranche evolution[i:j] est "error" ; chaque courbe garde la sienne, qui sert quand on l’interpole seule.

zone_a = pyrucast.SubEvolution([(0.0, 1.0), (1.0, 2.0)], out_of_range="clamp")
zone_b = pyrucast.SubEvolution([(0.0, 3.0), (1.0, 4.0)])
deux = zone_a | zone_b
print(deux.interpolate(0.5))  # [1.5, 3.5] : un flottant par courbe
print(deux.out_of_range)  # error : l'agrégat applique la sienne à ses courbes
print(deux.interpolate(2.0, out_of_range="clamp"))  # [2.0, 4.0]

Le long d’une ligne

pyrucast.evolution.along(field, line, components=None) (aussi field.along(line)) rend les valeurs d’un champ aux nœuds le long d’une ligne : une courbe par composante, la valeur en fonction de la distance parcourue. line est un maillage POI1, lu dans l’ordre de ses mailles : cet ordre est le trajet. L’abscisse d’un point est la longueur du trajet jusqu’à lui — la somme des distances en ligne droite d’un point au suivant —, 0 au premier point ; son ordonnée, la valeur de la composante en ce nœud. components=None prend toutes les composantes du champ, dans son ordre ; chacune donne une courbe, dont l’ordinate_type est son nom.

Un champ d’un autre maillage, ou un champ aux points de Gauss, se porte d’abord sur les points de la ligne par pyrucast.node_field.project. Deux points consécutifs confondus sont refusés : l’abscisse se répéterait.

# Une température relevée le long d'un trajet de trois points : la courbe porte
# la valeur en fonction de la distance parcourue — 0, 3, puis 3 + 4 = 7.
plan = pyrucast.Coords(2)
trajet = [plan.add_node(p) for p in ([0.0, 0.0], [3.0, 0.0], [3.0, 4.0])]
ligne = pyrucast.mesh.poi1_from_nodes(trajet)
releve = pyrucast.NodeField(ligne, ["T"])
for noeud, t in zip(trajet, [20.0, 50.0, 90.0]):
    releve[0].set_value(noeud, "T", t)
le_long = pyrucast.evolution.along(releve, ligne)  # aussi releve.along(ligne)
print(le_long[0].abscissas())  # [0.0, 3.0, 7.0]

Tracer

plot() dessine une évolution de nombres en courbe, une ligne par zone, et une évolution de champs en champ sur le maillage, avec un curseur pour choisir l’abscisse dans la fenêtre interactive (voir Visualisation). Sans x_label ni y_label, les axes et le curseur prennent l’abscissa_type et l’ordinate_type. Une table ne se dessine pas en courbe.

courbe = pyrucast.Evolution(
    [(0.0, 20.0), (10.0, 120.0)], abscissa_type="t", ordinate_type="T"
)
courbe.plot(save="courbe.svg")  # axes étiquetés "t" et "T"
chauffe.plot(save="chauffe.png", frame=1)  # le champ du pas tabulé n° 1

Sauvegarde

Une évolution se sauvegarde comme tout objet, avec sa politique hors plage et ses types. Les champs qu’elle tabule gardent leur support, partagé avec les autres objets du même fichier.

Côté Rust

SubEvolution porte :

• abscissas, un Vec<f64> plat : dims nombres par échantillon, les échantillons triés dans l’ordre de leurs tuples (première variable, puis deuxième…) ;
• dims, le nombre de variables, de 1 à MAX_ABSCISSAS (4) ;
• values, un Vec<SubValue>, une valeur par échantillon. SubValue vaut Scalar(f64), Node(SubNodeField) ou Element(SubElementField) : les champs sont rangés dans la courbe, et leur support est un Handle partagé ;
• out_of_range (OutOfRange), abscissa_type (Vec<String>) et ordinate_type (Option<String>).

SubEvolution::new prend des échantillons dont l’abscisse implémente Abscissa : un f64, un [f64; N] ou un Vec<f64>.

Evolution est l’agrégat de Handle<SubEvolution>, plus sa politique hors plage. check_push refuse une courbe dont les valeurs ne sont pas du même genre que celles des autres. Evolution::interpolate_at rend un Interpolated : Scalars(Vec<f64>), Node(NodeField) ou Element(ElementField). from_scalars, from_node_fields et from_element_fields sont les constructeurs que Python appelle.

Interpolation. Une interrogation pèse au plus deux échantillons par variable, donc au plus 16. gather les trouve par recherche dichotomique, variable après variable, et range leurs poids dans un tableau de taille fixe, sur la pile : les poids ne demandent aucune allocation. Pour des champs, la somme pondérée reprend l’arithmétique des champs (map_all, puis merge_components après check_same_components).

Appliquer à un champ. SubEvolution::interpolate_field et Evolution::interpolate_node_field / interpolate_element_field résolvent une fois, avant de parcourir le champ, quelles composantes chaque courbe lit et quelle politique elle applique. Les lignes du champ sont ensuite traitées en parallèle, sans allocation ; chaque ligne du résultat est écrite une fois, depuis sa propre ligne d’entrée, si bien que le résultat ne dépend pas du nombre de fils.

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

Détails des physiques

Une physique est un sous-modèle (SubModel) que construit un opérateur de pyrucast.model. L’opérateur reçoit son support, un espace éléments finis ou les maillages que relie une contrainte, et rend un modèle qui le couvre en entier : un sous-modèle par zone. Les physiques se composent par | : une conduction et la contrainte qui fige sa température sur un bord forment un seul modèle.

Cette page dit ce qu’ont en commun toutes les physiques : leur tag, leur type, leurs variables, et ce que chacune sait produire. Chaque page qui suit décrit une physique selon le même plan (voir le plan des pages).

Les tags

Chaque sous-modèle porte un ou plusieurs tags, qui disent à quelle physique il appartient. sub.physics() les donne, et un modèle se restreint à un tag (voir Modèle). Les pages sont rangées par tag :

Les échanges et le flux imposé n’ont pas de tag propre : ils prennent celui de la physique qu’ils visent, thermal sur une conduction, mechanical sur une élasticité. Une contrainte porte constraint et le tag de la physique dont elle fige les inconnues.

Les types de modèle

Une physique est de l’un de trois types, selon ce qu’elle sait faire :

Ce que les accesseurs d’un sous-modèle répondent :

Seul un domaine avec comportement se prête aux opérateurs de comportement de pyrucast.element_field : strain, integrate_behavior, et la tangente pyrucast.matrix.tangent (voir Comportement).

Les variables

Chaque physique nomme ses variables :

• les primales (primal_vars()) sont les inconnues, les colonnes de la matrice : T, u_x… ;
• les duales (dual_vars()) sont les lignes, où se lit le second membre : q, f_x… La duale d’une primale est ce qui travaille avec elle : un flux de chaleur avec une température, une force avec un déplacement ;
• le matériau est un champ aux points de Gauss qui porte les composantes requises (material_components()) et, au besoin, des composantes facultatives ;
• l’entrée du comportement est ce que strain(modele, solution, materiaux) calcule depuis la solution : un gradient, une déformation, des courbures ;
• la sortie du comportement est ce que rend integrate_behavior(modele, entrée, materiaux, prev, dt) : le flux, la contrainte ou les efforts, puis l’état que la loi relira au pas suivant (prev) ;
• les forces internes (pyrucast.node_field.internal_forces) sont ce que la sortie exerce sur les duales : ∫ Bᵀσ, où B est l’opérateur qui produit l’entrée depuis la solution. Chaque page dit comment sa physique les calcule.

Une contrainte n’a qu’une primale, lambda (le multiplicateur), et une duale, imposed (la ligne de la relation, où s’écrit la valeur imposée).

Ce que chaque physique produit

Les tableaux disent, pour chaque physique, les matrices qu’elle apporte et comment s’obtient sa tangente. Un tiret dit qu’elle n’y contribue rien : une masse assemblée sur un modèle qui contient un échange ignore simplement l’échange.

La tangente est ce que rend pyrucast.matrix.tangent. « = raideur » dit que la physique répond à une demande de tangente par sa raideur ; « analytique » qu’elle a une formule ; « perturbation » qu’elle la calcule par différences finies sur sa loi ; « aucune » qu’elle refuse la demande.

thermal

diffusion

Le tag de la cible

mechanical

constraint

Le plan des pages

Chaque physique est décrite dans cet ordre :

1. En bref : ce qu’elle modélise, son opérateur Python, son tag, son type, et les éléments qu’elle accepte ;
2. Théorie : les équations, de la forme locale à la matrice élémentaire ;
3. Variables : primales, duales, matériau, entrée et sortie du comportement, état ;
4. Ce qu’elle produit : matrices, tangente, forces internes et externes ;
5. Exemple en Python ;
6. Côté Rust.

Une page qui regroupe plusieurs lois d’une même physique (les lois plastiques, de fluage, d’endommagement) donne l’en-bref et les variables communes une fois, puis la théorie et le matériau loi par loi.

Côté Rust

Une physique est un struct de src/models/ qui implémente le trait SubModelKind : ses variables (primal_vars, dual_vars), ses tags (physics), son nom (label), et ses capacités, qu’elle expose par trois méthodes :

• as_domain() rend Some(&dyn Domain) pour un domaine : le sous-espace où lire le matériau (material_fespace), les composantes requises et facultatives, les matrices élémentaires ;
• as_behavior() rend Some(&dyn Behavior) pour un domaine qui porte une loi : son entrée (deformation_reads, produite par strain_operator), sa sortie (behavior_output_components), l’état qu’elle relit (state_reads), sa tangente au point (tangent_point, tangent_source), et requires_dt pour une loi qui dépend du pas de temps ;
• as_constraint() rend Some(&dyn Constraint) pour une contrainte : ses nœuds multiplicateurs et ses relations.

Le code générique (assembleur, comportement, sauvegarde) ne passe que par ce trait : SubModel::as_kind est le seul match sur les physiques. Les matrices qu’une physique apporte se déclarent par matrix_layout(kind), une par genre (MatrixKind : raideur, masse, géométrique, tangente) ; une contrainte ou un transfert d’interface rend ses blocs par contributions. Pour en écrire une, voir Ajouter une physique.

Conduction thermique

En bref

La conduction de la chaleur dans un solide : la température T est l’inconnue, et la chaleur circule des zones chaudes vers les zones froides.

symmetry choisit la forme de la conductivité : "isotropic" (un scalaire k), "orthotropic" ou "anisotropic" (un tenseur orienté par des axes matériau).

Théorie

En régime transitoire, l’équation de la chaleur s’écrit

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

avec \(\mathbf K\) le tenseur de conductivité, \(\rho\,c_p\) la capacité volumique et \(Q\) une source. En régime permanent, le premier terme disparaît. Multipliée par une température virtuelle \(v\) et intégrée par parties, elle devient

\[ \int_\Omega \nabla v\cdot\mathbf K\,\nabla T\,d\Omega + \int_\Omega v\,\rho\,c_p\,\dot T\,d\Omega = \int_\Omega v\,Q\,d\Omega + \int_\Gamma v\,\bar q\,d\Gamma , \]

où \(\bar q\) est le flux qui entre par le bord. Un bord sans condition est isolé : rien n’y entre.

Avec \(T \approx \sum_j N_j\,T_j\), chaque maille donne une matrice de conduction et une matrice de capacité :

\[ K_{ij} = \int_e \nabla N_i^{\mathsf T}\,\mathbf K\,\nabla N_j\,d\Omega , \qquad C_{ij} = \int_e \rho\,c_p\,N_i\,N_j\,d\Omega , \]

intégrées aux points de Gauss, où le matériau est lu : la conductivité et la capacité peuvent varier à l’intérieur d’une maille. Le système semi-discret est \(C\,\dot T + K\,T = F\) ; un pas d’Euler implicite résout \((C/\Delta t + K)\,T_{n+1} = F + (C/\Delta t)\,T_n\).

Conductivité orientée. Un matériau fibré ou feuilleté conduit mieux dans certaines directions. Le tenseur \(\mathbf K\) est alors construit dans les axes matériau, diagonal (orthotropic, conductivités principales k_1, k_2, k_3) ou complet (anisotropic, k_11…k_33), puis tourné dans les axes globaux. Les axes sont des vecteurs du champ matériau : V1X, V1Y en 2-D, le second axe étant la normale dans le plan ; V1X…V1Z et V2X…V2Z en 3-D, le troisième étant V1 × V2. Ils sont normés en interne et peuvent varier d’une maille à l’autre. En 2-D, k_3 (la direction normale au plan) est demandé aussi : le tenseur est construit en 3-D, puis réduit au plan.

Variables

La sortie est 𝐊·∇T, le flux de la forme faible. Le flux de Fourier, celui qui va du chaud vers le froid, est son opposé : −𝐊·∇T.

Ce qu’elle produit

Exemple

Une barre de conductivité 2, chauffée par un flux de 10 à gauche, à 20 à droite. Le profil est linéaire, \(T(x) = 20 + Q\,(1 - x)/k\) :

import pyrucast

# Une barre [0, 1] de quatre SEG2 : un flux Q = 10 entre à gauche, T = 20 à droite.
c = pyrucast.Coords(1)
gauche = c.add_node([0.0])
droite = c.add_node([1.0])
barre = pyrucast.mesh.line(gauche, droite, 4)
fes = pyrucast.FiniteElementSpace(barre)

conduction = pyrucast.model.heat_conduction(fes)
appui = pyrucast.model.dirichlet(
    pyrucast.mesh.poi1_from_nodes([droite]), conduction, "T"
)
modele = conduction | appui
materiaux = pyrucast.element_field.material_field(
    modele, [("k", 2.0), ("rho", 1.0), ("cp", 3.0)]
)

source = pyrucast.NodeField(barre, ["q"])  # sur les nœuds de la barre
source.unit().set_value(gauche, "q", 10.0)
second_membre = source | appui.constraint_rhs([(droite, 20.0)])

K = pyrucast.matrix.stiffness(modele, materiaux)
T = K.solve(second_membre)
print(T.value(gauche, "T"))  # 25.0 = 20 + Q·L/k
print(T.value(droite, "T"))  # 20.0

Le gradient et le flux aux points de Gauss, et les forces internes, qui rendent la source :

gradient = pyrucast.element_field.strain(modele, T, materiaux)
print(gradient.components())  # ['grad_T_x']
print(gradient[0].value(0, 0, "grad_T_x"))  # -5.0

repos = pyrucast.element_field.initial_state(
    modele, materiaux
)  # l'état de départ de la loi
flux = pyrucast.element_field.integrate_behavior(modele, gradient, materiaux, repos)
print(flux[0].value(0, 0, "flux_x"))  # -10.0 = k·∂T/∂x

f_int = pyrucast.node_field.internal_forces(modele, flux, T, materiaux)
print(f_int.value(gauche, "q"))  # 10.0 : la source que la conduction équilibre

Un pas de temps, depuis une barre à 20 partout :

C = pyrucast.matrix.mass(modele, materiaux)  # la capacité ρ·cp ; l'appui n'y entre pas
T0 = pyrucast.NodeField(barre, ["T"])
T0.unit().set_uniform("T", 20.0)

dt = 0.01
A = C / dt + K  # Euler implicite : (C/Δt + K)·T₁ = F + (C/Δt)·T₀
T1 = A.solve(second_membre + (C / dt) * T0)
print(round(T1.value(gauche, "T"), 3))  # 20.306 : la chaleur commence à entrer

Une conductivité orthotrope, d’axe principal (0,6 ; 0,8) :

c2 = pyrucast.Coords(2)
bas = pyrucast.mesh.line(c2.add_node([0.0, 0.0]), c2.add_node([1.0, 0.0]), 2)
haut = pyrucast.mesh.line(c2.add_node([0.0, 1.0]), c2.add_node([1.0, 1.0]), 2)
carre = pyrucast.mesh.sweep(bas, haut, 2)
plaque = pyrucast.FiniteElementSpace(carre)

orientee = pyrucast.model.heat_conduction(plaque, symmetry="orthotropic")
print(orientee[0].material_components())  # ['k_1', 'k_2', 'k_3', 'V1X', 'V1Y']
materiaux = pyrucast.element_field.material_field(
    orientee,
    [("k_1", 12.0), ("k_2", 3.0), ("k_3", 3.0), ("V1X", 0.6), ("V1Y", 0.8)],
)
K = pyrucast.matrix.stiffness(orientee, materiaux)

Côté Rust

HeatConduction (src/models/heat_conduction.rs) garde son sous-espace et sa MaterialSymmetry. Les matrices élémentaires sont des fonctions pures, element_stiffness et element_capacity, qui lisent le matériau ligne par ligne aux points de Gauss par les indices résolus une fois par zone (ElementLayout). Le tenseur orienté est construit par src/models/symmetry.rs, le même code que pour la diffusion et l’élasticité.

Le même problème que l’exemple, en Rust :

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::symmetry::MaterialSymmetry;
use pyrucast::ops::model;
use pyrucast::ops::solver::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 constraint mints its multiplier beside the imposed node (a fresh
    // co-located node); its value is written there.
    let imposed = Mesh::from_submesh(SubMesh::poi1_from_nodes(std::slice::from_ref(
        nodes.last().unwrap(),
    ))?);

    let conduction = model::heat_conduction(&fes, MaterialSymmetry::Isotropic)?;
    let dirichlet = model::dirichlet(&imposed, &conduction, "T", Default::default())?;
    let mult = dirichlet.multiplier_mesh()?.node(0, 0, 0)?.id();
    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" 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".into(), "q".into()])?;
    rhs.set_value(node0, "q", Q)?;
    rhs.set_value(mult, "imposed", 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")?;
    assert!(
        (reaction - Q).abs() < tol,
        "réaction λ : obtenue {reaction}, attendue {Q}"
    );

    Ok(())
}

Le test de la conduction orientée impose un champ linéaire au bord, qui doit être reproduit à l’intérieur pour tout tenseur constant, puis relit le flux et le compare à 𝐊·∇T calculé à la main :

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::model;
use pyrucast::ops::solver::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(())
}

Rayonnement

En bref

Le rayonnement d’une surface vers un environnement lointain à la température \(T_\infty\) (loi de Stefan-Boltzmann). Il se pose sur le bord d’une conduction et partage ses variables : il entre directement dans sa matrice.

fespace est l’espace du bord, target la conduction qu’il refroidit. Une cible qui n’assemble pas T avec sa duale q est refusée : le rayonnement écrirait dans une ligne que personne ne résout.

Les deux tags permettent de garder le rayonnement dans le problème thermique (tag thermal) ou de le prendre seul (tag radiation) pour l’assembler ou l’inspecter à part.

Théorie

Le flux qui quitte la surface vaut

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

où \(\sigma\) est la constante de Stefan-Boltzmann et \(\varepsilon\) l’émissivité. Aucune normale n’entre dans le calcul : la direction est déjà dans \(q\cdot n\), il ne reste qu’un scalaire et la mesure de la surface.

La loi n’est pas linéaire en \(T\), d’où trois termes :

Linéarisée autour de \(T_\infty\), la raideur est une matrice constante : c’est l’opérateur de départ d’une boucle de Newton, et à elle seule un opérateur de point fixe utilisable. La tangente porte la non-linéarité : elle est évaluée à la température de chaque point de Gauss.

sigma vaut par défaut la constante SI, 5,670374419·10⁻⁸ W·m⁻²·K⁻⁴, et T est alors une température absolue, en kelvins : une puissance quatrième ne permet pas de déplacer l’origine. Dans un autre système d’unités, on donne sigma dans le matériau.

Variables

Ce qu’il produit

Dans un calcul

Le flux se calcule aux points de Gauss de la face. La température y est interpolée à partir de ses valeurs aux nœuds, \( T_g = \sum_j N_j(x_g)\,T_j \), élevée à la puissance quatrième en ce point, et le flux \( \sigma\varepsilon(T_g^4 - T_\infty^4) \) est intégré avec le poids \( w_g \) de chaque point (la mesure de la surface comprise) :

\[ R_i = \sum_g w_g\, N_i(x_g)\, \sigma\varepsilon\,(T_g^4 - T_\infty^4), \qquad (K_r)_{ij} = \sum_g w_g\, 4\sigma\varepsilon\,T_g^3\, N_i(x_g)\, N_j(x_g) \]

\( R \) est le résidu, que rendent les forces internes, et \( K_r \) sa dérivée exacte par rapport aux températures nodales : la tangente.

Le calcul pas à pas annule le résidu par corrections successives. À chaque pas de temps, pyrucast.stepping part de la température du pas précédent et calcule le résidu de tout le problème — conduction, capacité, chargements, rayonnement. Il résout la correction avec un opérateur qui contient la tangente \( K_r \), prise une fois au début du pas, puis recommence avec le résidu de la nouvelle température, jusqu’à ce que la correction tombe sous precision (\( 10^{-4} \) par défaut) relativement à la température. Le résidu suit donc la température à chaque itération ; l’opérateur, lui, reste celui du début du pas. Comme il contient la dérivée de \( T^4 \), une ou deux corrections suffisent quand le pas de temps est raisonnable (algorithme détaillé dans Calcul pas à pas, Algo 2).

La raideur, film linéarisé autour de \( T_\infty \), est une matrice constante : employée seule à chaque itération, elle donne une itération de point fixe, qui converge plus lentement que la tangente.

Une autre méthode consiste à écrire le rayonnement comme une convection dont le coefficient dépend de la température. Puisque \( T^4 - T_\infty^4 = (T^2 + T_\infty^2)(T + T_\infty)(T - T_\infty) \), le flux vaut \( h_r\,(T - T_\infty) \), avec

\[ h_r = \sigma\varepsilon\,(T^2 + T_\infty^2)(T + T_\infty). \]

On calcule \( h_r \) aux nœuds, avec la température nodale de l’itération en cours, on l’interpole aux points de Gauss, et chaque itération résout un problème de convection : \( \int_\Gamma h_r N_i N_j \) dans la matrice, \( \int_\Gamma h_r\,T_\infty N_i \) au second membre. On recommence avec la nouvelle température jusqu’à ce qu’elle ne bouge plus.

Deux choses changent. \( h_r \) est une sécante, pas la tangente \( 4\sigma\varepsilon T^3 \) : chaque itération corrige moins, et il en faut davantage pour la même précision. Et le problème discret n’est pas tout à fait le même : \( h_r \) interpolé depuis les nœuds n’est pas \( h_r \) calculé avec la température des points de Gauss. Les deux méthodes convergent vers la même solution quand le maillage s’affine ; sur un maillage donné, leurs solutions convergées diffèrent de l’ordre de l’erreur de discrétisation.

Exemple

Une plaque à 1000 K qui rayonne par son bord supérieur vers 300 K. On lit le flux au bord, puis l’écart entre tangente et raideur au coin du bord : un seul segment, de longueur L = 0,5, touche ce coin, où \(\int N_0^2 = L/3\) :

conduction = pyrucast.model.heat_conduction(plaque)
rayonnement = pyrucast.model.radiation(pyrucast.FiniteElementSpace(haut), conduction)
print(rayonnement[0].physics())  # ['thermal', 'radiation']
modele = conduction | rayonnement
materiaux = pyrucast.element_field.material_field(
    modele, [("k", 20.0), ("emis", 0.8), ("T_inf", 300.0)]
)

T = pyrucast.NodeField(carre, ["T"])
T.unit().set_uniform("T", 1000.0)  # en kelvins
entree = pyrucast.element_field.strain(modele, T, materiaux)  # ∇T dedans, T au bord
repos = pyrucast.element_field.initial_state(modele, materiaux)
sortie = pyrucast.element_field.integrate_behavior(modele, entree, materiaux, repos)
print(sortie[1].value(0, 0, "flux"))  # 44995.55… = σε(T⁴ − T∞⁴)

K = pyrucast.matrix.stiffness(modele, materiaux)  # film linéarisé en T∞
Kt = pyrucast.matrix.tangent(
    modele, materiaux, entree, repos
)  # film à la température T
coin = haut.node(0, 0, 0)
print(
    Kt.get(coin, "q", coin, "T") - K.get(coin, "q", coin, "T")
)  # 29.425… = 4σε(T³ − T∞³)·L/3

Côté Rust

Radiation (src/models/radiation.rs) garde son sous-espace de bord. Son noyau de raideur évalue 4σεT∞³ au point de Gauss, son comportement (integrate_point) le flux, et tangent_point la dérivée 4σεT³. STEFAN_BOLTZMANN est la constante par défaut.

Le test compare le flux rayonné à σε(T⁴ − T∞⁴) fois l’aire, et la tangente à une différence finie du résidu :

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::symmetry::MaterialSymmetry;
use pyrucast::models::Physics;
use pyrucast::ops::element_field;
use pyrucast::ops::model;
use pyrucast::ops::solver::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,
        &element_field::behavior::initial_state(&fixture.radiation, &materials)?,
        &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(())
}

Diffusion (loi de Fick)

En bref

La diffusion d’une espèce dans un milieu : de l’humidité dans un béton, de l’hydrogène dans un acier, des chlorures dans un enrobage. La concentration est l’inconnue, et l’espèce va des zones riches vers les zones pauvres.

species nomme l’espèce : toutes ses grandeurs portent ce nom (c_H2, j_H2, D_H2), si bien que deux espèces vivent sur le même maillage sans se confondre. Ce qui appartient au milieu ne porte pas le nom de l’espèce : le coefficient de stockage poro et les axes matériau. symmetry choisit la forme de la diffusivité, comme pour la conduction.

L’équation est celle de la conduction, mais la physique est distincte : son tag est diffusion, et non thermal, pour qu’un problème qui couple les deux puisse prendre l’une sans l’autre.

Théorie

La première loi de Fick relie le flux de l’espèce au gradient de concentration, et la conservation de l’espèce donne, avec un coefficient de stockage \(\varphi\) (la porosité pour une espèce qui diffuse dans un solide poreux) :

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

C’est l’équation de la chaleur avec \(c\) pour \(T\), \(\mathsf D\) pour \(\mathbf K\) et \(\varphi\) pour \(\rho\,c_p\). Sa forme faible et sa discrétisation sont les mêmes :

\[ K_{ij} = \int_e \nabla N_i^{\mathsf T}\,\mathsf D\,\nabla N_j\,d\Omega , \qquad C_{ij} = \int_e \varphi\,N_i\,N_j\,d\Omega . \]

La diffusivité \(\mathsf D\) est un scalaire (isotropic) ou un tenseur orienté par les axes matériau (orthotropic, anisotropic), construit comme la conductivité.

Variables

Pour l’espèce H2 :

La sortie est 𝖣·∇c, le flux de la forme faible ; le flux de Fick est son opposé. Elle est nommée d’après la duale de l’espèce, si bien qu’un modèle qui porte une conduction et plusieurs espèces garde des flux distincts.

Ce qu’elle produit

Exemple

Une barre de diffusivité 0,5 et de longueur 2, où un flux d’hydrogène de 1 entre à gauche, la concentration étant nulle à droite. Le multiplicateur de l’appui est le flux qui sort, égal au flux qui entre :

import pyrucast

# Une barre [0, 2] : un flux d'hydrogène de 1 entre à gauche, c = 0 à droite.
c = pyrucast.Coords(1)
gauche = c.add_node([0.0])
droite = c.add_node([2.0])
barre = pyrucast.mesh.line(gauche, droite, 4)
fes = pyrucast.FiniteElementSpace(barre)

fick = pyrucast.model.fick(fes, "H2")
print(fick[0].primal_vars(), fick[0].dual_vars())  # ['c_H2'] ['j_H2']
appui = pyrucast.model.dirichlet(pyrucast.mesh.poi1_from_nodes([droite]), fick, "c_H2")
modele = fick | appui
materiaux = pyrucast.element_field.material_field(modele, [("D_H2", 0.5)])

source = pyrucast.NodeField(barre, ["j_H2"])
source.unit().set_value(gauche, "j_H2", 1.0)
c_h2 = pyrucast.matrix.stiffness(modele, materiaux).solve(
    source | appui.constraint_rhs([(droite, 0.0)])
)
print(round(c_h2.value(gauche, "c_H2"), 9))  # 4.0 = j·L/D
print(
    round(c_h2.value(appui.multiplier_mesh().node(0, 0, 0), "lambda"), 9)
)  # 1.0 : le flux qui sort

gradient = pyrucast.element_field.strain(modele, c_h2, materiaux)
repos = pyrucast.element_field.initial_state(
    modele, materiaux
)  # l'état de départ de la loi
flux = pyrucast.element_field.integrate_behavior(modele, gradient, materiaux, repos)
print(flux.components(), flux[0].value(0, 0, "j_H2_x"))  # ['j_H2_x'] -1.0 = D·∂c/∂x

Deux espèces et une conduction sur le même maillage forment un système à trois blocs indépendants :

melange = (
    pyrucast.model.fick(fes, "H2")
    | pyrucast.model.fick(fes, "O2")
    | pyrucast.model.heat_conduction(fes)
)
print(melange.primal_vars())  # ['c_H2', 'c_O2', 'T']
print(len(melange.filter("diffusion")), len(melange.filter("thermal")))  # 2 1
materiaux = pyrucast.element_field.material_field(
    melange, [("D_H2", 0.5), ("D_O2", 0.1), ("k", 2.0)]
)
K = pyrucast.matrix.stiffness(melange, materiaux)  # trois blocs, sans couplage

Une diffusivité qui dépend de la température se calcule par une évolution appliquée au champ de température, avant chaque assemblage. Un transfert d’espèce à travers une interface entre deux corps est un transfert d’interface.

Côté Rust

Fick (src/models/fick.rs) garde son sous-espace, sa MaterialSymmetry et le nom de l’espèce. Ses noyaux, element_stiffness et element_storage, ont la forme de ceux de la conduction et partagent avec eux la construction du tenseur orienté (symmetry::transport_tensor_from). Les noms se construisent depuis l’espèce : fick::primal_var("H2") rend c_H2, fick::dual_var rend j_H2.

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::symmetry::MaterialSymmetry;
use pyrucast::models::{Physics, IMPOSED, MULTIPLIER};
use pyrucast::ops::model;
use pyrucast::ops::solver::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 diffusion = model::fick(&fes, SPECIES, MaterialSymmetry::Isotropic)?;
    let dirichlet = model::dirichlet(
        &imposed,
        &diffusion,
        &primal_var(SPECIES),
        Default::default(),
    )?;
    let mult = dirichlet.multiplier_mesh()?.node(0, 0, 0)?.id();
    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![IMPOSED.to_string(), dual_var(SPECIES)])?;
    rhs.set_value(node0, &dual_var(SPECIES), J)?;
    rhs.set_value(mult, IMPOSED, 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, MULTIPLIER)?;
    assert!((reaction - J).abs() < tol, "réaction λ : {reaction} ≠ {J}");
    Ok(())
}

Échanges et flux imposé

Trois physiques sans comportement, qui s’appliquent à une autre physique, leur cible, et en prennent les variables et le tag :

• l’échange de bord, entre un bord et un milieu extérieur connu ;
• le transfert d’interface, entre deux maillages qui se font face ;
• le flux imposé, une densité donnée sur un bord ou un volume.

Les deux échanges suivent la même loi : un flux proportionnel à un écart,

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

Pour l’échange de bord, \(b\) est une donnée, la valeur ambiante ; pour le transfert d’interface, c’est une inconnue, le champ de l’autre maillage. Les deux partagent leur noyau, transfer::exchange_matrix.

Ce qui est échangé est donné par l’appelant. On passe la liste des variables échangées (["T"], ["c_H2"], ["u_x", "u_y"]). Chacune doit être une inconnue de la cible, qui lui donne sa duale, la ligne où le flux est écrit. Le coefficient \(h\) est un par variable, nommé d’après elle : h_T, h_c_H2, h_u_x. L’échange prend le tag de la cible : thermal sur une conduction, diffusion sur une diffusion, mechanical sur une élasticité. Un nom inconnu de la cible, ou des variables de deux tags dans un même échange, sont refusés à la construction.

Échange de bord

En bref

Un échange entre un bord et un milieu extérieur : la convection avec un fluide (loi de Newton), le transfert de masse en surface, ou, sur les déplacements, une fondation élastique.

Théorie

Sur le bord \(\Gamma\), le flux sortant vaut \(h\,(a - a_\text{ext})\). Dans la forme faible, il donne une matrice de film et un second membre :

\[ K_{ij} = \int_\Gamma h\,N_i\,N_j\,d\Gamma , \qquad f_i = \int_\Gamma h\,a_\text{ext}\,N_i\,d\Gamma . \]

Aucune normale n’entre dans le calcul : elle est déjà dans \(q\cdot n\), et il ne reste sous l’intégrale qu’un scalaire et la mesure du bord, qui ne dépend pas de l’orientation des mailles. Le film rend la matrice définie : un problème avec un flux imposé et un échange, sans aucune température imposée, a une solution unique.

Variables

Pour la variable T d’une conduction :

Ce qu’il produit

Exemple

Un carré de conductivité 2, avec une source de 10 sur le bord gauche (un flux imposé) et un film h = 5 vers 20 sur le bord droit. Tout ce qui entre ressort par le film, d’où \(T = T_\text{ext} + Q/h\) à droite et \(T_\text{ext} + Q/h + Q/k\) à gauche :

import pyrucast

# Un carré [0, 1]² : ses bords gauche et droit sont les deux lignes balayées.
c = pyrucast.Coords(2)
gauche = pyrucast.mesh.line(c.add_node([0.0, 0.0]), c.add_node([0.0, 1.0]), 2)
droite = pyrucast.mesh.line(c.add_node([1.0, 0.0]), c.add_node([1.0, 1.0]), 2)
carre = pyrucast.mesh.sweep(gauche, droite, 2)

conduction = pyrucast.model.heat_conduction(pyrucast.FiniteElementSpace(carre))
source = pyrucast.model.flux(pyrucast.FiniteElementSpace(gauche), conduction, "q")
film = pyrucast.model.boundary_transfer(
    pyrucast.FiniteElementSpace(droite), conduction, ["T"]
)
print(source[0].physics(), film[0].physics())  # ['thermal'] ['thermal']
modele = conduction | source | film
materiaux = pyrucast.element_field.material_field(
    modele, [("k", 2.0), ("phi_q", 10.0), ("h_T", 5.0), ("a_ext_T", 20.0)]
)

second_membre = pyrucast.node_field.external_forces(
    modele, materiaux
)  # φ à gauche, h·T_ext à droite
T = pyrucast.matrix.stiffness(modele, materiaux).solve(second_membre)
print(round(T.value(droite.node(0, 0, 0), "T"), 9))  # 22.0 = T_ext + Q/h
print(round(T.value(gauche.node(0, 0, 0), "T"), 9))  # 27.0 = T_ext + Q/h + Q/k

Le même échange sur les déplacements est une fondation élastique : un ressort réparti par direction.

plaque = pyrucast.model.elasticity(pyrucast.FiniteElementSpace(carre), "plane_stress")
semelle = pyrucast.model.boundary_transfer(
    pyrucast.FiniteElementSpace(gauche), plaque, ["u_x", "u_y"]
)
print(semelle[0].physics())  # ['mechanical']
print(semelle[0].material_components())  # ['h_u_x', 'h_u_y', 'a_ext_u_x', 'a_ext_u_y']

Transfert d’interface

En bref

Un échange entre deux maillages qui se font face sans partager leurs nœuds : une résistance de contact, un revêtement, un joint, une membrane. Le champ saute à la traversée, et le flux qui passe est proportionnel au saut.

side_a et side_b sont les espaces des deux faces. Elles doivent être conformes : même type de maille, même nombre de mailles, la maille i face à la maille i, et leur nœud local k à moins de tol l’un de l’autre. C’est vérifié à la construction, et une interface non conforme est refusée. Un sous-modèle est créé par paire de zones des deux côtés, prises dans l’ordre.

Théorie

La forme faible du terme d’interface est

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

qui donne quatre blocs sur les inconnues des deux côtés :

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

Les deux blocs diagonaux sont deux films, un par côté. Les deux autres ont leurs lignes sur un maillage et leurs colonnes sur l’autre ; chacun seul n’est pas symétrique, leur paire l’est. Les variables échangées ne se couplent pas entre elles : la chaleur qui traverse un joint ne pousse pas l’hydrogène.

C’est la version par pénalité de ce qu’un multi-points impose exactement : quand \(h\) grandit, le saut tend vers zéro. Si \(h\) est une donnée physique (une résistance de contact mesurée), c’est un transfert ; s’il n’est là que pour lier deux surfaces, c’est une contrainte.

Variables

Pour la variable T :

Ce qu’il produit

Exemple

Deux carrés côte à côte, un joint h = 4 entre eux, un flux de 1 entrant à gauche et T = 0 à droite. Le saut à la traversée vaut \(Q/h\) :

# Deux carrés côte à côte ; la ligne x = 1 existe deux fois, avec des nœuds distincts.
c = pyrucast.Coords(2)
x0 = pyrucast.mesh.line(c.add_node([0.0, 0.0]), c.add_node([0.0, 1.0]), 2)
x1_g = pyrucast.mesh.line(c.add_node([1.0, 0.0]), c.add_node([1.0, 1.0]), 2)
x1_d = pyrucast.mesh.line(c.add_node([1.0, 0.0]), c.add_node([1.0, 1.0]), 2)
x2 = pyrucast.mesh.line(c.add_node([2.0, 0.0]), c.add_node([2.0, 1.0]), 2)
corps = pyrucast.model.heat_conduction(
    pyrucast.FiniteElementSpace(pyrucast.mesh.sweep(x0, x1_g, 2))
) | pyrucast.model.heat_conduction(
    pyrucast.FiniteElementSpace(pyrucast.mesh.sweep(x1_d, x2, 2))
)
joint = pyrucast.model.interface_transfer(
    pyrucast.FiniteElementSpace(x1_g), pyrucast.FiniteElementSpace(x1_d), corps, ["T"]
)
source = pyrucast.model.flux(pyrucast.FiniteElementSpace(x0), corps, "q")
fixe = x2.to_poi1()  # les trois nœuds de x = 2
appui = pyrucast.model.dirichlet(fixe, corps, "T")
modele = corps | joint | source | appui
materiaux = pyrucast.element_field.material_field(
    modele, [("k", 1.0), ("h_T", 4.0), ("phi_q", 1.0)]
)

second_membre = pyrucast.node_field.external_forces(
    modele, materiaux
) | appui.constraint_rhs([(fixe.node(0, i, 0), 0.0) for i in range(3)])
T = pyrucast.matrix.stiffness(modele, materiaux).solve(second_membre)
saut = T.value(x1_g.node(0, 0, 0), "T") - T.value(x1_d.node(0, 0, 0), "T")
print(round(saut, 9))  # 0.25 = Q/h : le champ saute à la traversée

Flux imposé

En bref

Une densité donnée, intégrée en forces nodales : une source de chaleur sur un bord ou dans un volume, une traction, une pression. Elle ne dépend pas de la solution, et n’entre dans aucune matrice.

dual est la ligne de la cible où la densité est écrite : "q" pour une source de chaleur, "f_x" pour une traction en x. "n" désigne une pression : la densité agit selon la normale de la maille et s’écrit dans toutes les lignes de déplacement (f_x, f_y, f_z). La normale est celle de la configuration active : sur une configuration déformée, la pression suit la paroi.

Théorie

\[ f_i = \int \varphi\,N_i\,d\Gamma , \qquad f_{i,a} = \int \varphi\,n_a\,N_i\,d\Gamma \;\text{pour une pression.} \]

Variables

La densité est une composante du champ matériau : uniforme, ou variable d’un point de Gauss à l’autre.

Ce qu’il produit

Exemple

La source de l’échange de bord ci-dessus est un flux imposé : model.flux(gauche, conduction, "q") et la composante phi_q du matériau.

Côté Rust

BoundaryTransfer (src/models/boundary_transfer.rs) et InterfaceTransfer (src/models/interface_transfer.rs) s’appuient sur src/models/transfer.rs, qui porte leur noyau commun, exchange_matrix, la vérification des variables contre la cible et le choix du tag (owner_physics). Ils implémentent Domain sans Behavior. Ce qu’ils lisent pour leurs forces internes est la primale interpolée aux points de Gauss (residual_input), et non un état.

Le transfert d’interface rend ses quatre blocs par contributions : deux blocs ordinaires et deux blocs de couplage (Contribution::Coupling), dont les lignes et les colonnes vivent sur deux sous-espaces différents.

Flux (src/models/flux.rs) ne contribue à aucune matrice ; il déclare seulement son terme de forces externes (external_force_contribution).

Le test de la fondation élastique compare le déplacement d’une barre poussée sur un ressort réparti à \(u = q / (E/L + h)\), sur quatre décades de raideur :

/// 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(())
}

Le test du transfert d’interface vérifie le profil de part et d’autre et le saut \(q/h\) :

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::symmetry::MaterialSymmetry;
use pyrucast::models::{Physics, IMPOSED};
use pyrucast::ops::model;
use pyrucast::ops::solver::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(())
}

Mécanique

Les physiques de tag mechanical. Leurs inconnues sont les déplacements u_x, u_y, u_z (et les rotations r_x, r_y, r_z des poutres et des coques), leurs duales les forces nodales f_x, f_y, f_z (et les moments m_x, m_y, m_z).

Milieux continus, une même cinématique (B, la déformation tirée du déplacement) et une loi qui change :

• Élasticité linéaire — elasticity, et sa symétrie matériau orthotrope ou anisotrope ;
• Plasticité de von Mises — plasticity_perfect, plasticity_isotropic ;
• Drucker-Prager et Ottosen — des surfaces qui voient la pression ;
• Fluage et viscoplasticité — Norton, Blackburn, Lemaitre, Chaboche, Lemaitre-Chaboche ;
• Endommagement de Mazars ;
• Autres lois d’endommagement — Damage TC, SiC/SiC, Gurson ;
• Grandes transformations — les mesures de déformation green_lagrange et hencky.

Éléments de structure, dont les efforts sont intégrés sur la section :

• Barre — truss, l’effort axial seul ;
• Poutre d’Euler-Bernoulli — bernoulli ;
• Poutre de Timoshenko — timoshenko, le cisaillement en plus ;
• Coques — shell, six inconnues par nœud.

Ce que partagent les milieux continus

Une loi de milieu continu se déclare par pyrucast.model.<loi>(fespace, kinematics, …, strain_measure="linear").

kinematics est l’hypothèse cinématique :

L’axisymétrie demande des Coords de révolution (voir Coordonnées), et réciproquement : un modèle plan est refusé sur une géométrie de révolution. La déformation suit la convention de l’ingénieur : le glissement \(\gamma_{xy} = 2\varepsilon_{xy}\) est rangé dans eps_xy.

Les lois plastiques et de fluage intègrent leur état en 3-D, quelle que soit la cinématique, et rendent en plus de la contrainte la déformation totale et, pour les modèles plans, sigma_zz. La sortie d’un pas de toute loi non linéaire est complète : elle sert telle quelle de prev au pas suivant.

Le matériau accepte deux composantes facultatives : rho (la masse) et alpha (la dilatation thermique, lue par pyrucast.element_field.thermal_strain).

Barre

En bref

Une barre ne transmet que l’effort axial : un treillis, un câble tendu, un tirant. Elle se pose sur des segments, en 1-D, 2-D ou 3-D.

Une barre n’a aucune raideur transversale : chaque nœud doit être tenu dans les directions que la barre ne porte pas (par d’autres barres, un treillis triangulé, ou des appuis), sinon la matrice est singulière. Tout autre élément qu’un SEG2 est refusé, comme une géométrie de révolution.

Théorie

Le long de l’axe \(s\) de la barre, de section \(A\) :

\[ N = E\,A\,\varepsilon, \qquad \varepsilon = \frac{du_s}{ds}, \qquad \frac{dN}{ds} + f = 0 . \]

La direction de la barre est le vecteur unitaire \(c = (x_B - x_A)/L\), tiré des coordonnées des nœuds. Avec l’interpolation linéaire, la déformation est constante dans l’élément, et la raideur s’écrit directement dans les axes globaux, sans matrice de repère :

\[ K_e = \frac{E\,A}{L} \begin{bmatrix} c\,c^{\mathsf T} & -c\,c^{\mathsf T} \\ -c\,c^{\mathsf T} & c\,c^{\mathsf T} \end{bmatrix}. \]

En 1-D, \(c = 1\) et l’on retrouve \((EA/L)\,[[1, -1], [-1, 1]]\). La masse est \((\rho A L/6)\,[[2, 1], [1, 2]]\) sur chaque direction, et la raideur géométrique sous l’effort \(N\) vaut \((N/L)\,[[I, -I], [-I, I]]\) : elle raidit toutes les directions, l’axiale comprise.

Avec strain_measure="green_lagrange", la barre suit les grands déplacements : voir Grandes transformations.

Variables

Selon la dimension des Coords :

n est l’effort axial \(E\,A\,c^{\mathsf T}\varepsilon\,c\), lu le long de la barre dans la déformation du point.

Ce qu’elle produit

Exemple

Une barre de longueur 2, de module 200 000 et de section 0,01, tirée par une force de 1000 : \(u = F\,L/(E\,A) = 1\).

import pyrucast

# Une barre de longueur 2, fixée à gauche, guidée en y à droite, tirée par F = 1000.
c = pyrucast.Coords(2)
a = c.add_node([0.0, 0.0])
b = c.add_node([2.0, 0.0])
ligne = pyrucast.mesh.line(a, b, 1)
barre = pyrucast.model.truss(pyrucast.FiniteElementSpace(ligne))
appuis = (
    pyrucast.model.dirichlet(pyrucast.mesh.poi1_from_nodes([a]), barre, "u_x")
    | pyrucast.model.dirichlet(pyrucast.mesh.poi1_from_nodes([a]), barre, "u_y")
    | pyrucast.model.dirichlet(pyrucast.mesh.poi1_from_nodes([b]), barre, "u_y")
)
modele = barre | appuis
materiaux = pyrucast.element_field.material_field(modele, [("E", 200e3), ("A", 0.01)])

force = pyrucast.NodeField(ligne, ["f_x"])
force.unit().set_value(b, "f_x", 1000.0)
u = pyrucast.matrix.stiffness(modele, materiaux).solve(force)  # valeurs imposées nulles
print(u.value(b, "u_x"))  # 1.0 = F·L/(E·A)

deformation = pyrucast.element_field.strain(modele, u, materiaux)
repos = pyrucast.element_field.initial_state(
    modele, materiaux
)  # l'état de départ de la loi
effort = pyrucast.element_field.integrate_behavior(
    modele, deformation, materiaux, repos
)
print(effort.components(), effort[0].value(0, 0, "n"))  # ['n'] 1000.0

Côté Rust

Truss (src/models/truss.rs) calcule ses matrices par des fonctions pures, element_stiffness, element_mass et element_geometric, et ses forces internes par internal_force_element, qui applique le même \(B\) que la raideur : pour une loi linéaire, \(\int B^{\mathsf T}\sigma\) vaut exactement \(K\,u\).

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::element_field::StrainMeasure;
use pyrucast::ops::model;
use pyrucast::ops::solver::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))?);
        model::dirichlet(&imposed, target, var, Default::default())
    };
    let mut model = model::truss(&fes, StrainMeasure::Linear)?;
    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(())
}

Élasticité linéaire

En bref

Un solide qui se déforme réversiblement, proportionnellement à la charge (loi de Hooke), en petites déformations.

kinematics est l’hypothèse cinématique (voir Mécanique), symmetry la symétrie du matériau (voir Orthotropie), strain_measure la mesure de déformation (voir Grandes transformations). Cette page décrit le cas isotrope en petites déformations.

Un maillage de bord (des segments en 2-D) est refusé comme domaine : \(B\) y serait bâti sur le seul gradient tangent, et la raideur n’aurait aucune rigidité dans la direction normale.

Théorie

Sur le domaine \(\Omega\), avec \(b\) les forces de volume :

\[ \nabla\cdot\sigma + b = 0, \qquad \sigma = \mathbb D : \varepsilon, \qquad \varepsilon = \tfrac12\big(\nabla u + \nabla u^{\mathsf T}\big). \]

La forme faible : pour tout déplacement virtuel \(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 \(\Gamma_N\).

Avec \(u \approx \sum_i N_i\,u_i\), la déformation s’écrit \(\varepsilon = B\,u_e\). En 2-D, \(\varepsilon = [\varepsilon_{xx}, \varepsilon_{yy}, \gamma_{xy}]\) et 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}; \]

en 3-D, \(\varepsilon = [\varepsilon_{xx}, \varepsilon_{yy}, \varepsilon_{zz}, \gamma_{yz}, \gamma_{xz}, \gamma_{xy}]\) avec un bloc de six lignes de même forme. La raideur et la masse de l’élément sont

\[ K_e = \int_{\Omega_e} B^{\mathsf T} D\,B\,d\Omega , \qquad M_e = \int_{\Omega_e} \rho\,N^{\mathsf T} N\,d\Omega , \]

intégrées aux points de Gauss, où le matériau est lu.

La matrice de Hooke \(D\), avec \(E\) le module d’Young et \(\nu\) le coefficient de Poisson :

• plane_stress (\(\sigma_{zz} = 0\)), avec \(c = E/(1-\nu^2)\) : \(D = c\,[[1, \nu, 0], [\nu, 1, 0], [0, 0, (1-\nu)/2]]\) ;
• plane_strain (\(\varepsilon_{zz} = 0\)), avec \(c = E/((1+\nu)(1-2\nu))\) : \(D = c\,[[1-\nu, \nu, 0], [\nu, 1-\nu, 0], [0, 0, (1-2\nu)/2]]\) ;
• full_3d, même \(c\) : le bloc normal \(3\times3\) de \(c(1-\nu)\) sur la diagonale et \(c\nu\) ailleurs, et \(G = c\,(1-2\nu)/2\) sur les trois cisaillements ;
• axisymmetric, même \(c\) : le bloc normal isotrope sur \([rr, zz, \theta\theta]\) et \(G\) sur \(rz\).

Axisymétrie. Un solide de révolution se maille dans son plan méridien \((r, z)\). Deux choses changent. La mesure d’intégration devient \(2\pi r\,|J|\), pour toutes les intégrales (raideur, masse, forces) ; c’est la géométrie qui l’apporte, pas la physique. Et la déformation orthoradiale \(\varepsilon_{\theta\theta} = u_r / r\) ajoute une ligne à \(B\), \([N_i/r,\; 0]\) ; elle est rangée dans eps_zz. Les points de Gauss étant intérieurs aux mailles, \(r > 0\) même pour une maille qui touche l’axe.

Dilatation thermique. Une température impose une déformation libre \(\varepsilon_{th} = \alpha\,(T - T_\text{ref})\), sans effet de la mécanique sur la thermique. La raideur ne change pas ; le terme thermique agit sur la contrainte et le second membre :

\[ \sigma = D : (\varepsilon(u) - \varepsilon_{th}), \qquad f_{th} = \int_\Omega B^{\mathsf T} D\,\varepsilon_{th}\,d\Omega . \]

alpha est une composante facultative du matériau ; pyrucast.element_field.thermal_strain(temperature, materiaux, fespace, t_ref) calcule \(\varepsilon_{th}\) aux points de Gauss depuis une température portée aux points de Gauss par interp_to_gauss ; integrate_behavior puis internal_forces en tirent la charge thermique.

Variables

Isotrope, en petites déformations, selon la cinématique :

En axisymétrie, eps_zz et sigma_zz sont les composantes orthoradiales.

Ce qu’elle produit

La raideur géométrique est décrite avec les opérateurs d’assemblage.

Exemple

Une plaque carrée en contraintes planes, tenue en u_x sur son bord gauche et en u_y à son coin inférieur gauche, tirée par une traction \(S = 100\) sur son bord droit. La solution est exacte : \(u_x = S\,x/E\), \(u_y = -\nu\,S\,y/E\) :

gauche, droite, carre = _carre()  # un carré [0, 1]² de 2×2 QUA4
plaque = pyrucast.model.elasticity(pyrucast.FiniteElementSpace(carre), "plane_stress")
coin = gauche.node(0, 0, 0)
modele = (
    plaque
    | pyrucast.model.dirichlet(gauche.to_poi1(), plaque, "u_x")
    | pyrucast.model.dirichlet(pyrucast.mesh.poi1_from_nodes([coin]), plaque, "u_y")
    | pyrucast.model.flux(pyrucast.FiniteElementSpace(droite), plaque, "f_x")
)
materiaux = pyrucast.element_field.material_field(
    modele,
    [("E", 200e3), ("nu", 0.3), ("phi_f_x", 100.0)],  # traction S = 100 à droite
)
u = pyrucast.matrix.stiffness(modele, materiaux).solve(
    pyrucast.node_field.external_forces(modele, materiaux)
)
haut_droit = droite.node(0, 1, 1)
print(round(u.value(haut_droit, "u_x"), 12))  # 0.0005 = S/E
print(round(u.value(haut_droit, "u_y"), 12))  # -0.00015 = −ν·S/E

eps = pyrucast.element_field.strain(modele, u, materiaux)
repos = pyrucast.element_field.initial_state(modele, materiaux)
sigma = pyrucast.element_field.integrate_behavior(modele, eps, materiaux, repos)
print(sigma.components())  # ['sigma_xx', 'sigma_yy', 'sigma_xy']
print(round(sigma[0].value(0, 0, "sigma_xx"), 9))  # 100.0

Côté Rust

Elasticity (src/models/elasticity.rs) délègue sa cinématique à un Continuum (src/models/continuum/), qui porte \(B\) pour chaque hypothèse et le noyau de forces internes ; la matrice \(D\) vient de src/models/symmetry.rs, l’isotropie gardant ses formes fermées.

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::symmetry::MaterialSymmetry;
use pyrucast::models::tensor::Kinematics;
use pyrucast::ops::element_field::StrainMeasure;
use pyrucast::ops::model;
use pyrucast::ops::solver::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)?);
        model::dirichlet(&imposed, target, var, 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,
        MaterialSymmetry::Isotropic,
        StrainMeasure::Linear,
    )?;
    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")?)?;

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

Sur un solide de révolution, tests/axisymmetric.rs vérifie la solution de Lamé d’un cylindre épais : avec QUA8, l’erreur sur les déplacements nodaux décroît à l’ordre 4 environ et celle des contraintes à l’ordre 2, sans jamais s’annuler, \(1/r\) n’étant pas un polynôme.

La barre chauffée vérifie la dilatation thermique : encastrée aux deux bouts, \(\sigma_{xx} = -E\,\alpha\,\Delta T\) ; libre, une dilatation sans contrainte :

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::symmetry::MaterialSymmetry;
use pyrucast::models::tensor::Kinematics;
use pyrucast::ops::element_field::{deformation, interp_to_gauss, thermal_strain, StrainMeasure};
use pyrucast::ops::model;
use pyrucast::ops::node_field::internal_forces;
use pyrucast::ops::solver::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)?);
        model::dirichlet(&imposed, target, var, 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,
        MaterialSymmetry::Isotropic,
        StrainMeasure::Linear,
    )?;
    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, StrainMeasure::Linear)?;
    let sig_th = pyrucast::ops::element_field::behavior::integrate(
        &model,
        &eps_th,
        &pyrucast::ops::element_field::behavior::initial_state(&model, &materials)?,
        &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, StrainMeasure::Linear)?;
    let eps_mech = eps.merge_field(&eps_th, |a, b| a - b)?;
    let rest = pyrucast::ops::element_field::behavior::initial_state(&model, &materials)?;
    let sigma = pyrucast::ops::element_field::behavior::integrate(
        &model, &eps_mech, &rest, &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(())
}

Élasticité orthotrope et anisotrope

En bref

Un composite tissé, un bois, une tôle laminée ou un monocristal n’ont pas la même raideur dans toutes les directions. C’est l’option symmetry de l’élasticité, indépendante de la cinématique :

Les variables ne changent pas : seule la matrice de Hooke \(D\) change. Les lois non linéaires (plasticité, fluage, endommagement de Mazars et Damage TC) gardent une élasticité isotrope ; la loi SiC/SiC a ses propres axes.

Théorie

La loi de comportement est \(\sigma_{ij} = C_{ijkl}\,\varepsilon_{kl}\). Les symétries de \(\sigma\) et \(\varepsilon\) et l’existence d’un potentiel élastique réduisent \(C\) à une matrice \(6\times6\) symétrique en notation de Voigt : 21 constantes dans le cas général.

Orthotropie. Dans ses axes propres, la souplesse \(S = C^{-1}\) découple les termes normaux des cisaillements :

\[ 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}, \qquad \frac{\nu_{ji}}{E_j} = \frac{\nu_{ij}}{E_i} . \]

Les relations de réciprocité ne laissent que neuf constantes. \(S\) doit être définie positive : un jeu de constantes qui la rend singulière est refusé, plutôt que d’assembler une raideur fausse.

Les axes matériau sont donnés par des vecteurs du champ matériau : V1X, V1Y en 2-D (le second axe est la normale dans le plan, le troisième la direction hors du plan, l’orthoradiale en axisymétrie) ; V1X…V1Z et V2X…V2Z en 3-D (le troisième est \(V_1 \times V_2\)). Ils sont orthonormalisés en interne : \(V_2\) n’a besoin d’être qu’à peu près perpendiculaire à \(V_1\). Un \(V_1\) nul, ou un \(V_2\) parallèle à \(V_1\), est refusé. Comme toute composante matériau, les axes peuvent varier d’une maille et d’un point de Gauss à l’autre : un composite bobiné, une pièce courbe.

En chaque point de Gauss, \(D\) est construit dans les axes matériau, tourné dans les axes globaux par le tenseur d’ordre 4,

\[ C'_{pqrs} = R_{pi}\,R_{qj}\,R_{rk}\,R_{sl}\,C_{ijkl}, \qquad R = \big[\,V_1\ \ V_2\ \ V_1\times V_2\,\big], \]

puis réduit à la cinématique : le bloc [xx, yy, xy] en déformations planes, sa condensation sur \(\varepsilon_{zz}\) en contraintes planes, le bloc [rr, zz, θθ, rz] en axisymétrie. En cisaillement de l’ingénieur, le passage Voigt ↔ tenseur ne porte aucun facteur, et la rotation d’ordre 4 n’a aucune convention à retenir.

Variables

Les neuf constantes orthotropes sont demandées même en 2-D : le tenseur est construit en 3-D avant d’être réduit. Les trois jeux de composantes sont disjoints, si bien qu’une zone isotrope et une zone orthotrope partagent un maillage sans rien consolider. L’entrée et la sortie du comportement sont celles de l’élasticité isotrope.

Exemple

Les constantes d’un matériau quatre fois plus raide selon son axe 1 :

gauche, droite, carre = _carre()
plaque = pyrucast.model.elasticity(
    pyrucast.FiniteElementSpace(carre), "plane_stress", symmetry="orthotropic"
)
print(plaque[0].material_components())
# ['E_1', 'E_2', 'E_3', 'nu_12', 'nu_13', 'nu_23', 'G_12', 'G_13', 'G_23', 'V1X', 'V1Y']
angle = 0.0  # l'axe raide, selon x
constantes = [
    ("E_1", 40e3),
    ("E_2", 10e3),
    ("E_3", 10e3),
    ("nu_12", 0.25),
    ("nu_13", 0.25),
    ("nu_23", 0.3),
    ("G_12", 5e3),
    ("G_13", 5e3),
    ("G_23", 4e3),
    ("V1X", math.cos(angle)),
    ("V1Y", math.sin(angle)),
]

Sur la plaque en traction de la page Élasticité, le fichier de tests vérifie les deux bornes : \(u_x = S/E_1\) quand l’axe raide suit la traction, \(S/E_2\) quand il lui est perpendiculaire.

Côté Rust

src/models/symmetry.rs construit les tenseurs orientés pour l’élasticité, la conduction et la diffusion ; MaterialSymmetry est l’énumération des trois symétries. Le test ci-dessous étire un carré orthotrope selon son premier axe matériau ; le même fichier nourrit une loi orthotrope de constantes isotropes, puis une loi anisotrope du tenseur isotrope, qui doivent se comporter comme l’isotrope quel que soit le repère, et refuse un repère dégénéré.

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;
use pyrucast::coords::Coords;
use pyrucast::handle::Handle;
use pyrucast::models::symmetry::MaterialSymmetry;
use pyrucast::models::tensor::Kinematics;
use pyrucast::ops::element_field::StrainMeasure;
use pyrucast::ops::model;
use pyrucast::ops::solver::solve;
use pyrucast::Result;

/// Traction on the right edge.
const S: f64 = 2.0;
/// `N×N` QUA4 grid on the unit square.
const N: usize = 2;

#[test]
fn orthotropic_square_stretches_along_its_first_material_axis() -> Result<()> {
    const E1: f64 = 200.0; // stiff direction — aligned with global x
    const E2: f64 = 50.0; // compliant transverse direction
    const NU12: f64 = 0.25;

    let (grid, fes, _coords) = unit_square()?;
    let model = clamped_model(&grid, &fes, MaterialSymmetry::Orthotropic)?;

    // The material axes travel through the material field like any other
    // coefficient: `V1` is the first orthotropy direction, here the global x.
    let materials = pyrucast::ops::element_field::material_field(
        &model,
        &[
            ("E_1", E1),
            ("E_2", E2),
            ("E_3", E2),
            ("nu_12", NU12),
            ("nu_13", NU12),
            ("nu_23", 0.25),
            ("G_12", 30.0),
            ("G_13", 30.0),
            ("G_23", 30.0),
            ("V1X", 1.0),
            ("V1Y", 0.0),
            ("phi_f_x", S),
        ],
    )?;

    let solution = solve_traction(&model, &materials)?;

    // σ_xx = S with σ_yy = σ_zz = 0 ⇒ the compliance row gives ε directly.
    let tol = 1e-10;
    let h = 1.0 / N as f64;
    for j in 0..=N {
        for i in 0..=N {
            let (x, y) = (i as f64 * h, j as f64 * h);
            let id = grid[j * (N + 1) + i].id();
            let ux = solution.value(id, "u_x")?;
            let uy = solution.value(id, "u_y")?;
            assert!((ux - S / E1 * x).abs() < tol, "u_x({x},{y}) = {ux}");
            assert!((uy + NU12 * S / E1 * y).abs() < tol, "u_y({x},{y}) = {uy}");
        }
    }
    Ok(())
}

Plasticité de von Mises

En bref

Un métal qui se déforme de façon irréversible au-delà d’une limite d’élasticité. Deux lois, au critère de von Mises :

Théorie

L’élastoplasticité en petites déformations partage la déformation en une part élastique et une part plastique :

\[ \varepsilon = \varepsilon^e + \varepsilon^p, \qquad \sigma = D : (\varepsilon - \varepsilon^p), \qquad f(\sigma, p) \le 0, \qquad \dot\varepsilon^p = \dot\lambda\,\frac{\partial f}{\partial\sigma}, \]

avec les conditions de Kuhn-Tucker \(\dot\lambda \ge 0\), \(f \le 0\), \(\dot\lambda\,f = 0\). Le critère de von Mises porte sur la contrainte équivalente \(q = \sqrt{\tfrac32\,s : s}\), où \(s\) est le déviateur de \(\sigma\) :

\[ f(\sigma, p) = q - \big(\sigma_y + H\,p\big), \qquad \dot p = \dot\lambda , \]

\(p\) étant la déformation plastique cumulée et \(H\) le module d’écrouissage, nul pour la loi parfaite.

Le retour radial. Un pas part de l’état \(A\) (\(\varepsilon^p_A\), \(p_A\)) et de la déformation totale \(\varepsilon_B\) :

1. prédiction élastique, \(\sigma^{tr} = D : (\varepsilon_B - \varepsilon^p_A)\), et sa contrainte équivalente \(q^{tr}\) ;
2. si \(q^{tr} \le \sigma_y + H\,p_A\), le pas est élastique ;
3. sinon, la normale étant colinéaire au déviateur, la consistance est affine et se résout sans itérer :

\[ \Delta p = \frac{q^{tr} - (\sigma_y + H\,p_A)}{3\mu + H}, \qquad s_B = s^{tr}\Big(1 - \frac{3\mu\,\Delta p}{q^{tr}}\Big), \qquad \varepsilon^p_B = \varepsilon^p_A + \frac{3\,\Delta p}{2\,q^{tr}}\,s^{tr}, \]

\(\mu\) étant le module de cisaillement ; la pression ne change pas.

Le calcul est mené en 3-D quelle que soit la cinématique. En déformations planes, les composantes hors du plan de la déformation sont nulles ; en contraintes planes, la condition \(\sigma_{zz} = 0\) est résolue par une itération locale autour du retour radial ; en axisymétrie, la déformation orthoradiale est mesurée.

La tangente cohérente est la dérivée exacte du retour,

\[ D_\text{alg} = K\,I\otimes I + 2\mu\,\theta\,\mathbb I_\text{dev} - 2\mu\,\bar\theta\;\hat n\otimes\hat n, \qquad \theta = \frac{\sigma_y + H\,p_B}{q^{tr}}, \qquad \bar\theta = \frac{3\mu}{3\mu + H} - (1 - \theta), \]

avec \(\hat n = s^{tr}/\lVert s^{tr}\rVert\) et \(K\) le module de compressibilité. Le facteur \(\theta\) rend compte du pas fini : c’est ce qui donne à Newton sa convergence quadratique.

La raideur reste la raideur élastique \(\int B^{\mathsf T} D\,B\) : un opérateur constant, celui d’un Newton modifié. La tangente \(\int B^{\mathsf T} D_\text{alg} B\) est celle d’un Newton complet.

Variables

Selon la cinématique :

En axisymétrie, eps_zz et sigma_zz sont les composantes orthoradiales.

La sortie du comportement est la contrainte, suivie de l’état :

La déformation plastique et la déformation totale ont toujours six composantes, quelle que soit la cinématique. Au premier pas, prev est l’état de repos, element_field.initial_state(modele, materiaux), où tout est nul. Ensuite, la sortie d’un pas se passe telle quelle comme prev au pas suivant.

Ce qu’elle produit

Exemple

Une maille étirée à \(\varepsilon_{xx} = 2\cdot10^{-3}\), au-delà de \(\sigma_y/E\). La contrainte équivalente reste à \(\sigma_y\) :

fes, u = _un_qua4(2e-3)  # une maille, ε_xx = 2·10⁻³ : au-delà de σ_y/E
modele = pyrucast.model.plasticity_perfect(fes, "plane_stress")
materiaux = pyrucast.element_field.material_field(
    modele, [("E", 210e3), ("nu", 0.3), ("sigma_y", 250.0)]
)

eps = pyrucast.element_field.strain(modele, u, materiaux)
repos = pyrucast.element_field.initial_state(
    modele, materiaux
)  # premier pas : l'état de repos
etat = pyrucast.element_field.integrate_behavior(modele, eps, materiaux, repos)
s = {nom: etat[0].value(0, 0, nom) for nom in ("sigma_xx", "sigma_yy", "sigma_xy")}
q = math.sqrt(
    s["sigma_xx"] ** 2
    + s["sigma_yy"] ** 2
    - s["sigma_xx"] * s["sigma_yy"]
    + 3 * s["sigma_xy"] ** 2
)
print(round(q, 3))  # 250.0 : la contrainte équivalente est plafonnée à σ_y
print(etat[0].value(0, 0, "p") > 0)  # True : de la déformation plastique

# Le pas suivant repart de cet état.
etat = pyrucast.element_field.integrate_behavior(modele, eps, materiaux, prev=etat)
Kt = pyrucast.matrix.tangent(
    modele, materiaux, eps, prev=etat
)  # le module algorithmique

Avec écrouissage, elle vaut \(\sigma_y + H\,p\) :

modele = pyrucast.model.plasticity_isotropic(fes, "plane_stress")
materiaux = pyrucast.element_field.material_field(
    modele, [("E", 210e3), ("nu", 0.3), ("sigma_y", 250.0), ("H", 1000.0)]
)
eps = pyrucast.element_field.strain(modele, u, materiaux)
etat = pyrucast.element_field.integrate_behavior(
    modele, eps, materiaux, pyrucast.element_field.initial_state(modele, materiaux)
)
s = {nom: etat[0].value(0, 0, nom) for nom in ("sigma_xx", "sigma_yy", "sigma_xy")}
q = math.sqrt(
    s["sigma_xx"] ** 2
    + s["sigma_yy"] ** 2
    - s["sigma_xx"] * s["sigma_yy"]
    + 3 * s["sigma_xy"] ** 2
)
print(round(q - (250.0 + 1000.0 * etat[0].value(0, 0, "p")), 3))  # 0.0 : q = σ_y + H·p

La boucle de Newton complète (pas de charge, résidu \(F_\text{ext} - F_\text{int}\), résolution, report de l’état) s’écrit en Python avec ces opérateurs ; examples/plasticite_poutre_console.py en est un exemple, une console cisaillée jusqu’à plastifier à l’encastrement.

Côté Rust

Plasticity (src/models/plasticity.rs) porte la loi d’écoulement en attribut : une même physique sert les dix lois de cette page, de la page Drucker-Prager et Ottosen, de la page Fluage et la loi de Gurson. Les deux lois de von Mises vivent dans src/models/plasticity/von_mises.rs ; H = 0 redonne la loi parfaite par le même chemin de code. examples/plasticite_poutre_console.rs est la même console en Rust, et sert de banc de parallélisme.

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::plasticity::law::PlasticLaw;
use pyrucast::models::tensor::Kinematics;
use pyrucast::ops::element_field::{
    behavior::{initial_state, integrate},
    deformation, material_field, StrainMeasure,
};
use pyrucast::ops::matrix::tangent;
use pyrucast::ops::model;
use pyrucast::ops::node_field::internal_forces;
use pyrucast::Result;

const AXES: [&str; 3] = ["x", "y", "z"];

/// von Mises with linear isotropic hardening.
const ISOTROPIC: &[(&str, f64)] = &[
    ("E", 70_000.0),
    ("nu", 0.3),
    ("sigma_y", 200.0),
    ("H", 5_000.0),
];
/// A frictional, mildly dilatant material — `ψ < α`, so the flow is
/// non-associated.
const DRUCKER: &[(&str, f64)] = &[
    ("E", 20_000.0),
    ("nu", 0.2),
    ("friction", 0.3),
    ("k", 30.0),
    ("psi", 0.1),
];
/// Ottosen's classic concrete set, for a tensile/compressive strength ratio of
/// about 0.1.
const OTTOSEN: &[(&str, f64)] = &[
    ("E", 30_000.0),
    ("nu", 0.2),
    ("a", 1.2759),
    ("b", 3.1962),
    ("k_1", 11.7365),
    ("k_2", 0.9801),
    ("sigma_c", 30.0),
];

#[test]
fn isotropic_hardening_satisfies_its_consistency_condition() -> Result<()> {
    let cube = Cube::new(PlasticLaw::Isotropic, ISOTROPIC)?;
    // Well past yield, so the step is definitely plastic.
    let s = cube.stress(&uniaxial(0.02))?;
    let (q, p) = (von_mises(&s.sigma), s.p);
    assert!(p > 0.0, "the step must be plastic (p = {p})");
    let expected = 200.0 + 5_000.0 * p;
    assert!(
        (q - expected).abs() < 1e-6 * expected,
        "q = {q}, expected σ_y + H·p = {expected}"
    );
    Ok(())
}

Drucker-Prager et Ottosen

En bref

Deux lois plastiques pour les matériaux plus résistants en compression qu’en traction : sols, roches, poudres, béton. Leur surface voit la pression, que von Mises ignore.

Elles suivent le cadre de la plasticité de von Mises : partition de la déformation, critère \(f \le 0\), écoulement et conditions de Kuhn-Tucker. Elles en diffèrent par la surface \(f\), le potentiel \(g\) de l’écoulement (\(\dot\varepsilon^p = \dot\lambda\,\partial g/\partial\sigma\)), et la tangente.

Théorie

Les surfaces s’écrivent sur les invariants de la contrainte :

\[ I_1 = \operatorname{tr}\sigma, \qquad J_2 = \tfrac12\,s : s, \qquad J_3 = \det s, \qquad q = \sqrt{3J_2}, \qquad \cos 3\theta = \frac{3\sqrt3}{2}\,\frac{J_3}{J_2^{3/2}} , \]

\(\theta \in [0, \pi/3]\) étant l’angle de Lode, la direction de la contrainte dans le plan déviatorique. Von Mises ne voit que \(q\) ; Drucker-Prager ajoute \(I_1\), donc la pression ; Ottosen ajoute \(\theta\).

Drucker-Prager

Le cône

\[ f(\sigma) = \alpha\,I_1 + \beta\,q - k \]

a pour pente \(\alpha\) (friction), pour cohésion \(k\) et pour poids déviatorique \(\beta\) (beta, 1 par défaut). L’écoulement n’est pas associé : le potentiel a sa propre pente, la dilatance \(\psi\) (psi), et son propre poids déviatorique \(\delta\) (delta, 1 par défaut), \(g(\sigma) = \psi\,I_1 + \delta\,q\). \(\psi = \alpha\) redonne un écoulement associé, \(\psi = 0\) un écoulement sans variation de volume.

Sur le flanc du cône, \(q\) et \(I_1\) sont affines en le multiplicateur, et le retour se résout sans itérer. Un retour qui donnerait \(q_B < 0\) a dépassé l’axe hydrostatique : la contrainte se pose alors au sommet du cône, \(s_B = 0\), \(\sigma_{m,B} = k/(3\alpha)\), où la tangente est nulle : un corps entièrement à son sommet ne porte plus rien.

Écrouissage. La cohésion croît avec la déformation plastique cumulée, \(dk = H\,dp\) (H, 0 par défaut ; négatif, il adoucit), jusqu’à une surface ultime \(\alpha_u I_1 + \beta_u q = k_u\) (friction_ult, beta_ult, k_ult, par défaut égales aux valeurs initiales). La surface courante interpole les deux :

\[ \lambda = \operatorname{clamp}\!\Big(\frac{H\,p}{k_u - k}, 0, 1\Big), \qquad (\alpha, \beta, k)(p) = (1-\lambda)\,(\alpha, \beta, k) + \lambda\,(\alpha_u, \beta_u, k_u), \]

si bien que \(k(p) = k + H\,p\) tant que la surface ultime n’est pas atteinte. La consistance n’est alors plus linéaire : elle est résolue par un Newton encadré, dont la solution sans écrouissage est le premier itéré. Sans écrouissage, le retour reste fermé.

Ottosen

Une surface à quatre paramètres, faite pour le béton :

\[ f(\sigma) = a\,\frac{J_2}{\sigma_c^2} + \lambda(\theta)\,\frac{\sqrt{J_2}}{\sigma_c} + b\,\frac{I_1}{\sigma_c} - 1, \]

\[ \lambda(\theta) = \begin{cases} k_1 \cos\!\big[\tfrac13 \arccos(k_2 \cos 3\theta)\big] & \text{si } \cos 3\theta \ge 0,\\ k_1 \cos\!\big[\tfrac\pi3 - \tfrac13 \arccos(-k_2 \cos 3\theta)\big] & \text{si } \cos 3\theta < 0 . \end{cases} \]

\(\sigma_c\) est la résistance en compression simple. Le terme en \(J_2\) courbe les méridiens, le terme en \(I_1\) les incline, et \(\lambda(\theta)\) fait de la section déviatorique un triangle arrondi qui s’ouvre vers la compression.

Aucun retour fermé n’existe. Le retour passe par l’algorithme du plan sécant, qui n’a besoin que de \(f\) : partant du prédicteur, on linéarise le critère à l’itéré courant, on en tire un multiplicateur, on corrige, jusqu’à \(|f|\) sous une tolérance. La normale \(\partial f/\partial\sigma\) est obtenue par différences centrées sur \(f\), avec un pas de \(10^{-6}\,\sigma_c\).

La tangente par perturbation

Les deux lois obtiennent leur tangente par différences centrées sur le retour lui-même, une colonne par composante de déformation :

\[ 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 , \]

soit douze retours par point de Gauss. Les colonnes de cisaillement sont divisées par deux, pour passer de \(\partial\sigma/\partial\varepsilon_{ij}\) à \(\partial\sigma/\partial\gamma_{ij}\).

La tangente est symétrisée : seul son triangle supérieur est gardé. Celle d’un écoulement non associé, comme Drucker-Prager, ne l’est pas ; la symétriser coûte à Newton sa convergence quadratique sur cette loi. Celle d’Ottosen est doublement numérique (la normale, puis la tangente) : moins précise, elle suffit à faire converger Newton.

Variables

Selon la cinématique :

En axisymétrie, eps_zz et sigma_zz sont les composantes orthoradiales.

La sortie du comportement est la contrainte, suivie de l’état :

Ce qu’elles produisent

Exemple

La tangente de Drucker-Prager, assemblée sur une maille plastifiée :

modele = pyrucast.model.drucker_prager(fes, "plane_strain")
materiaux = pyrucast.element_field.material_field(
    modele, [("E", 20e3), ("nu", 0.2), ("friction", 0.3), ("k", 10.0), ("psi", 0.1)]
)
eps = pyrucast.element_field.strain(modele, u, materiaux)
repos = pyrucast.element_field.initial_state(modele, materiaux)
etat = pyrucast.element_field.integrate_behavior(modele, eps, materiaux, repos)
Kt = pyrucast.matrix.tangent(
    modele, materiaux, eps, repos
)  # par perturbation, symétrisée
print(Kt.symmetric)  # True

Côté Rust

Les deux lois vivent dans src/models/plasticity/drucker_prager.rs et src/models/plasticity/ottosen.rs, et déclarent leur tangente par TangentSource::Numerical. tests/plastic_laws.rs vérifie que le retour atterrit sur la surface, l’effondrement au sommet du cône, l’écart entre traction et compression d’Ottosen, et confronte chaque tangente à une différence finie des forces internes.

Fluage et viscoplasticité

En bref

Des lois où l’écoulement prend du temps : la contrainte peut sortir du domaine élastique, et c’est l’écart qui fixe la vitesse à laquelle elle y revient. Le résultat d’un pas dépend de sa durée.

integrate_behavior et matrix.tangent exigent le pas de temps dt pour ces lois, et lèvent une erreur sans lui : intégrer une loi visqueuse comme si elle était instantanée donnerait un nombre plausible et faux.

Théorie

Le cadre est celui de la plasticité, sans les conditions de Kuhn-Tucker : la vitesse d’écoulement n’est plus l’inconnue d’une condition de consistance, elle est donnée par une loi,

\[ \dot\varepsilon^{vp} = \dot p\,\frac32\,\frac{s}{q}, \qquad \dot p = g(q, p, \dots). \]

L’écoulement est déviatorique, sans variation de volume, et radial.

Un solveur pour toutes. Le pas est intégré par Euler implicite : la vitesse est évaluée en fin de pas. La direction étant celle du prédicteur, le déviateur ne fait que se contracter, \(q(\Delta p) = q^{tr} - 3\mu\,\Delta p\), et le pas se ramène à une équation scalaire,

\[ R(\Delta p) = \Delta p - \Delta t\,g\big(q(\Delta p), p_A + \Delta p, \dots\big) = 0 , \]

dont la racine est unique dans \([0, q^{tr}/3\mu]\). Les lois de vitesse étant raides, elle est cherchée par un Newton encadré, doublé d’une dichotomie.

Norton. \(\dot p = (q/K)^n\), sans seuil : toute contrainte flue. \(K\) est la contrainte de référence, \(n\) la sensibilité.

Lemaitre. \(\dot p = (q/K)^N\,p^{-M}\), \(M > 0\) : la déformation accumulée ralentit l’écoulement ; à contrainte constante, \(p(t) \propto t^{1/(1+M)}\). \(p\) a un plancher minuscule, pour que \(p^{-M}\) reste fini au premier pas.

Blackburn. La déformation primaire est une variable à part, p_prim :

\[ \dot p_\text{prim} = r\,\big(\varepsilon_\infty(q) - p_\text{prim}\big), \qquad \varepsilon_\infty(q) = A\,\sinh(\alpha q), \qquad \dot p = \dot p_\text{prim} + B\,\sinh(\beta q) . \]

Le primaire tend vers son asymptote à la vitesse \(r\) puis s’éteint ; le secondaire persiste. En Euler implicite, le primaire s’inverse en un pas, \(p^B_\text{prim} = (p^A_\text{prim} + \Delta t\,r\,A\,\sinh(\alpha q))/(1 + \Delta t\,r)\).

Chaboche. Une surface translatée par une contrainte de rappel \(X\) et agrandie par un écrouissage isotrope \(R\) :

\[ f = J(\sigma - X) - R - k, \qquad \dot p = \Big\langle \frac fK \Big\rangle^n, \qquad \dot X = \tfrac23 C\,\dot\varepsilon^{vp} - \gamma\,X\,\dot p, \qquad \dot R = b\,(Q - R)\,\dot p , \]

où \(\langle\cdot\rangle\) est la partie positive : sous la surface, le comportement est élastique. \(X\) produit l’effet Bauschinger et sature à \(C/\gamma\) ; \(R\) sature à \(Q\). La direction est gelée au prédicteur, et les deux variables sont implicites en \(\Delta p\) : \(X_B = (X_A + \tfrac23 C\,\Delta p\,\hat n)/(1 + \gamma\,\Delta p)\), \(R_B = (R_A + b\,Q\,\Delta p)/(1 + b\,\Delta p)\).

Lemaitre-Chaboche remplace la contrainte par la contrainte effective \(\tilde\sigma = \sigma/(1 - D)\), et fait croître l’endommagement avec l’écoulement :

\[ \dot D = \Big(\frac YS\Big)^s \dot p, \qquad Y = \frac{\tilde\sigma_\text{eq}^2}{2E} . \]

Un matériau endommagé flue plus vite, ce qui l’endommage davantage : c’est le fluage tertiaire. \(D\) ne décroît jamais et s’arrête à D_c.

La tangente est obtenue par perturbation et symétrisée, comme celle de Drucker-Prager.

Variables

Selon la cinématique :

En axisymétrie, eps_zz et sigma_zz sont les composantes orthoradiales.

La sortie du comportement est la contrainte, suivie d’un état commun aux cinq lois, puis de l’état propre à chacune :

Ce qu’elles produisent

Exemple

Une déformation maintenue sous Norton : sans dt, l’intégration est refusée ; avec, la contrainte se relaxe d’un pas à l’autre.

modele = pyrucast.model.creep_norton(fes, "plane_strain")
materiaux = pyrucast.element_field.material_field(
    modele, [("E", 150e3), ("nu", 0.3), ("K", 400.0), ("n", 5.0)]
)
eps = pyrucast.element_field.strain(modele, u, materiaux)  # déformation maintenue
repos = pyrucast.element_field.initial_state(modele, materiaux)
try:
    pyrucast.element_field.integrate_behavior(modele, eps, materiaux, repos)
except RuntimeError as erreur:
    print(erreur)  # Plasticity: this law is rate-dependent and needs a time increment…

etat = pyrucast.element_field.integrate_behavior(modele, eps, materiaux, repos, dt=1.0)
sigma_1 = etat[0].value(0, 0, "sigma_xx")
etat = pyrucast.element_field.integrate_behavior(modele, eps, materiaux, etat, dt=1.0)
print(etat[0].value(0, 0, "sigma_xx") < sigma_1)  # True : la contrainte se relaxe

Les paramètres sont des composantes du matériau. Une dépendance à la température se calcule par une évolution appliquée au champ de température.

Côté Rust

Les cinq lois vivent dans src/models/plasticity/viscous.rs, et déclarent requires_dt. Le test vérifie que Norton suit sa loi de vitesse sur un pas court, qu’une déformation maintenue se relaxe d’autant plus que le pas est long, que le primaire de Blackburn sature, que la contrainte de rappel de Chaboche plafonne sous \(C/\gamma\), que l’endommagement ne dépasse pas \(D_c\), et que toutes refusent d’intégrer sans dt.

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::plasticity::law::PlasticLaw;
use pyrucast::models::tensor::Kinematics;
use pyrucast::ops::element_field::{
    behavior::{initial_state, integrate},
    deformation, material_field, StrainMeasure,
};
use pyrucast::ops::model;
use pyrucast::Result;

const AXES: [&str; 3] = ["x", "y", "z"];

/// Norton: a steady creep with a strongly non-linear stress dependence.
const NORTON: &[(&str, f64)] = &[("E", 150_000.0), ("nu", 0.3), ("K", 400.0), ("n", 5.0)];

#[test]
fn norton_creep_follows_its_rate_law() -> Result<()> {
    let cube = Cube::new(PlasticLaw::CreepNorton, NORTON)?;
    // A step short enough that the stress barely relaxes, so the closed-form
    // rate at the trial stress is the right comparison.
    let dt = 1e-6;
    let state = cube.step(&uniaxial(2e-3), &cube.rest()?, dt)?;
    let q = von_mises(&state.sigma);
    let expected = dt * (q / 400.0_f64).powf(5.0);
    assert!(
        (state.p - expected).abs() < 1e-3 * expected,
        "Δp = {}, expected dt·(q/K)^n = {expected}",
        state.p
    );
    Ok(())
}

Endommagement de Mazars

En bref

Le béton qui se fissure : il perd de sa raideur, sans déformation permanente. Un seul scalaire d’endommagement \(D \in [0, 1)\) affaiblit la contrainte.

Théorie

La loi est sécante : la contrainte est la contrainte élastique affaiblie,

\[ \sigma = (1 - D)\,D_\text{el} : \varepsilon . \]

L’endommagement est piloté par la déformation équivalente, bâtie sur les parts positives des déformations principales,

\[ \tilde\varepsilon = \sqrt{\textstyle\sum_I \langle\varepsilon_I\rangle_+^2}, \qquad \kappa = \max\big(\varepsilon_{d0},\ \max_t \tilde\varepsilon\big), \]

où \(\varepsilon_{d0}\) est le seuil. \(\kappa\) ne décroît jamais : le matériau ne guérit pas. L’endommagement mêle une branche de traction et une branche de compression,

\[ D_t = 1 - \frac{\varepsilon_{d0}(1 - A_t)}{\kappa} - \frac{A_t}{e^{B_t(\kappa - \varepsilon_{d0})}}, \qquad D = \alpha_t\,D_t + \alpha_c\,D_c , \]

\(D_c\) s’écrivant de même avec \(A_c\) et \(B_c\). Les poids \(\alpha_t\), \(\alpha_c\) viennent du partage de la contrainte effective \(D_\text{el} : \varepsilon\) en parts positive et négative, sur ses valeurs principales.

Le calcul est mené sur le tenseur 3-D : en déformations planes, \(\varepsilon_{zz} = 0\) ; en contraintes planes, \(\varepsilon_{zz} = -\nu\,(\varepsilon_{xx} + \varepsilon_{yy})/(1 - \nu)\) ; en axisymétrie, l’orthoradiale est mesurée.

La raideur est la raideur élastique non endommagée : l’opérateur d’itération d’un Newton modifié. La loi n’a pas de tangente.

Variables

Selon la cinématique :

En axisymétrie, eps_zz et sigma_zz sont les composantes orthoradiales.

La sortie du comportement est la contrainte, suivie de l’état :

La contrainte ne dépend que de la déformation totale et de \(\kappa\) : il n’y a pas d’incrément. Au premier pas, prev est l’état de repos, element_field.initial_state(modele, materiaux), où kappa vaut 0 : la loi prenant le maximum avec le seuil, cela revient à \(\kappa = \varepsilon_{d0}\).

Ce qu’elle produit

Exemple

Une maille étirée à deux fois le seuil :

fes, u = _un_qua4(2e-4)  # ε_xx = 2·10⁻⁴, au-delà du seuil ε_d0 = 10⁻⁴
modele = pyrucast.model.mazars(fes, "plane_stress")
materiaux = pyrucast.element_field.material_field(
    modele,
    [
        ("E", 30e3),
        ("nu", 0.2),
        ("eps_d0", 1e-4),
        ("A_t", 0.8),
        ("B_t", 20e3),
        ("A_c", 1.4),
        ("B_c", 1900.0),
    ],
)
eps = pyrucast.element_field.strain(modele, u, materiaux)
etat = pyrucast.element_field.integrate_behavior(
    modele, eps, materiaux, pyrucast.element_field.initial_state(modele, materiaux)
)
print(
    round(etat[0].value(0, 0, "kappa"), 6)
)  # 0.0002 : la plus grande déformation équivalente atteinte
print(0.0 < etat[0].value(0, 0, "damage") < 1.0)  # True

Côté Rust

Damage (src/models/damage.rs) porte sa loi en attribut (DamageLaw), comme Plasticity : Mazars, Damage TC et SiC/SiC. Mazars vit dans src/models/damage/mazars.rs ; la décomposition spectrale passe par nalgebra. tests/damage_laws.rs vérifie que Mazars ne rend qu’un endommagement scalaire.

Damage TC, SiC/SiC et Gurson

En bref

Trois lois où un matériau perd sa raideur en se dégradant :

Damage TC et SiC/SiC sont des lois d’endommagement, comme Mazars ; Gurson est une loi plastique.

Théorie

Le cadre de l’endommagement. La contrainte effective \(\tilde\sigma = D_\text{el} : \varepsilon\) est celle du matériau sain sous la même déformation ; la contrainte est sa version dégradée. Un moteur scalaire \(\tau\), construit sur \(\varepsilon\) ou \(\tilde\sigma\), alimente une variable d’histoire \(r = \max(r_0, \max_t \tau)\), qui ne décroît jamais, et l’endommagement est une fonction croissante de \(r\). Il n’y a pas de déformation permanente, et la mise à jour se fait en un pas, sans itérer.

Damage TC

Les parts positive et négative de la contrainte effective, découpées sur ses valeurs principales, sont dégradées chacune par sa variable :

\[ \sigma = (1 - d^+)\,\tilde\sigma^+ + (1 - d^-)\,\tilde\sigma^- . \]

Une fissure ouverte en traction qui se referme reprend donc toute la raideur en compression : c’est l’effet unilatéral, que le scalaire de Mazars ne peut pas rendre. Chaque variable a son moteur et son histoire :

\[ \tau^+ = \sqrt{\tilde\sigma^+ : \varepsilon}, \qquad \tau^- = \sqrt{\sqrt3\,\big|K\,\tilde\sigma^-_\text{oct} + \tilde\tau^-_\text{oct}\big|}, \qquad r_0^+ = \frac{f_t}{\sqrt E}, \quad r_0^- = \frac{f_c}{\sqrt E}, \]

l’énergie élastique de la part positive en traction, une mesure octaédrique (contrainte normale \(\tfrac13\operatorname{tr}\tilde\sigma^-\) et cisaillement) sensible au confinement en compression, avec \(K = 0{,}171\). L’adoucissement est exponentiel en traction, durcissant puis adoucissant en compression :

\[ d^+ = 1 - \frac{r_0^+}{r^+}\,e^{A_t(1 - r^+/r_0^+)}, \qquad d^- = 1 - \frac{r_0^-}{r^-}(1 - A_c) - A_c\,e^{2(1 - r^-/r_0^-)} . \]

SiC/SiC

Dans un composite tissé, la matrice fissure dans des plans normaux aux torons, et les fibres continuent de porter la charge : la raideur chute par direction. Chaque direction matériau \(i\) a son endommagement, piloté par la déformation dans les axes matériau \(\varepsilon^\text{mat} = R^{\mathsf T}\varepsilon R\) :

\[ \kappa_i = \max_t \langle\varepsilon^\text{mat}_{ii}\rangle_+, \qquad d_i = d_{\max,i}\Big(1 - e^{-(\kappa_i - \varepsilon_{0,i})/\varepsilon_{c,i}}\Big) \ \text{si } \kappa_i > \varepsilon_{0,i} . \]

Une direction comprimée n’est pas dégradée. La dégradation s’applique terme à terme dans les axes matériau, \(\sigma^\text{mat}_{ij} = \sqrt{1-d_i}\sqrt{1-d_j}\;\tilde\sigma^\text{mat}_{ij}\), ce qui garde l’opérateur symétrique. Chaque \(d_i\) sature à \(d_{\max,i}\) au lieu d’atteindre 1 : les fibres restent. Les axes sont ceux de l’orthotropie, les vecteurs V1 (et V2 en 3-D) du matériau ; l’élasticité, elle, est isotrope.

Gurson

La porosité \(f\) rétrécit la surface de charge :

\[ \Phi = \Big(\frac{q}{\sigma_y}\Big)^2 + 2\,q_1 f^* \cosh\!\Big(\frac{3\,q_2\,\sigma_m}{2\,\sigma_y}\Big) - \big(1 + q_3\,f^{*2}\big), \]

qui redonne von Mises à \(f = 0\). Le cosh rend la surface sensible à la pression : les cavités croissent en traction triaxiale. Au-delà d’une porosité critique \(f_c\), une porosité effective \(f^*\), bilinéaire en \(f\), accélère la coalescence et atteint \(1/q_1\) à la porosité de rupture \(f_f\). La porosité croît par conservation de la masse, \(\dot f = (1 - f)\,\operatorname{tr}\dot\varepsilon^p\), à partir de sa valeur initiale \(f_0\) ; la germination n’est pas modélisée. L’écoulement n’est donc pas isochore. Le retour se fait par plan sécant à normale numérique, comme celui d’Ottosen, et la tangente par perturbation.