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