Modèle physique (Model)

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

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

Architecture

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

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

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

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

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

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

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

Identification des DOFs : primal ≠ dual

Chaque physique déclare :

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

Elles sont presque toujours différentes :

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

Chargements complètement séparés du Model

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

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

Cette séparation a deux mérites :

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

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

Les physiques disponibles

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

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

Pour ajouter une physique, voir Ajouter une physique.

Ce que chaque physique calcule

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

Nature physique et filtrage

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

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

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

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

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

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

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

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

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

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

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

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

Règle invariante : un Model = une Matrice

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

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

API Rust

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

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

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

API Python

import pyrucast

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

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

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

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

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

Assemblage et résolution

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

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

Limitations actuelles

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